Diffuse-interface polycrystal plasticity: expressing grain boundaries as geometrically necessary dislocations

Abstract

The standard way of modeling plasticity in polycrystals is by using the crystal plasticity model for single crystals in each grain, and imposing suitable traction and slip boundary conditions across grain boundaries. In this fashion, the system is modeled as a collection of boundary-value problems with matching boundary conditions. In this paper, we develop a diffuse-interface crystal plasticity model for polycrystalline materials that results in a single boundary-value problem with a single crystal as the reference configuration. Using a multiplicative decomposition of the deformation gradient into lattice and plastic parts, i.e. F(X,t)=F ^L(X,t)F ^P(X,t), an initial stress-free polycrystal is constructed by imposing F ^L to be a piecewise constant rotation field R ⁰(X), and F ^P=R ⁰(X)^T, thereby having F(X,0)=I, and zero elastic strain. This model serves as a precursor to higher order crystal plasticity models with grain boundary energy and evolution.

Introduction

When a polycrystalline material is deformed, its microstructure generally experiences a reorientation of the crystal lattices of each grain towards a preferential distribution of orientations known as crystallographic texture. The study of texture evolution is important because textured metals typically exhibit plastic anisotropy, which plays a significant role on mechanical properties. Predicting the evolution of deformation-induced texture and the accompanying plastic anisotropy is the subject of polycrystal plasticity models (Beaudoin et al. 1993; Sarma and Dawson 1996; Kok et al. 2002; Estrin 2002). These models are typically formulated assuming that the microstructure of the polycrystal is associated with a representation of microscopic crystals whose individual responses, on average, determine the macroscopic response of the polycrystal. At the level of each grain, plastic deformation occurs by the standard mechanism of dislocation slip, and so (i) constitutive equations that relate dislocation motion to crystal deformation must be defined, and (ii) an averaging scheme that relates the response of individual crystals to the macroscopic stress-strain response of the polycrystal must also be defined. For single crystals, a multiplicative kinematic decomposition of the deformation gradient into elastic and plastic parts is typically used. This decomposition adequately describes the distinctly different kinematical mechanisms that operate during the plastic deformation of a crystal. It was formally introduced in continuum plasticity (Nemat-Nasser 1979; Simo 1988; Reina and Conti 2014), and then applied to describe the kinematics of single crystals (Asaro 1983; Lubarda 2004; Roters et al. 2010a). A feature of this decomposition is that it introduces an intermediate configuration between the reference and current configurations which is obtained by unloading the crystal to a stress-free state. The elasto-viscoplastic constitutive equations are generally written relative to this relaxed configuration.

Many numerical procedures have been proposed to integrate the crystal constitutive equations (Kalidindi et al. 1992; Cuitino and Ortiz 1993; Kuchnicki et al. 2006), generally implicit and semi-implicit procedures which are developed differently by particular selection of the primary variables (stresses (Harewood and McHugh 2007), shear rates (Zikry 1994), plastic deformation gradient (Rice 1971), etc.). Polycrystal plasticity models appear in various levels of sophistication. Along the venerable Sachs and Taylor models –in which the aggregate deformation (Sachs 1928; Kocks 1970a) or stress (Taylor and Quinney 1932; Hutchinson 1964) is computed by averaging from the individual crystal values–, self-consistent models have been developed and applied that express the global deformation in terms of linearized viscoplastic moduli that must be adjusted self-consistently (Lebensohn and Tomé 1993; Lebensohn et al. 2007; Acta Materialia 2012; International Journal of Plasticity 2013; Knezevic et al. 2014a). Models that spatially resolve grain boundaries (GB) have started to gain traction recently thanks to a higher efficiency of numerical solvers and a wider availability of computational resources. Roters et al. have provided a comprehensive review of the different variants of such approaches Roters et al. (2010b), which enable the calculation of the fine spatial features of strain and stress fields, including grain shape changes and nonuniform deformation. Some of these advances have also been discussed by Knezevic et al. (2014b).

However, in the above models, grain boundary processes –which are known to be relevant at high stresses and temperatures– cannot be captured by construction. For example, fundamental grain boundary properties such as energies and mobilities are extraneous to spatially-resolved standard (poly)crystal plasticity models.

The aim of this paper is to present a framework that preserves the ability to model intra-grain plasticity, while at the same time enabling a straightforward generalization to include grain boundary processes. To this end, we develop a ‘diffuse’-interface crystal plasticity model for polycrystalline materials based on a representation of grain boundaries as a special subclass of geometrically necessary dislocations (in the sense defined by Cermelli and Gurtin (2001, 2002)). In this model, with a single crystal as the reference configuration, a stress-free polycrystal is constructed by imposing a piecewise constant rotation field and its transpose as the lattice and plastic distortions respectively. To make the resulting model numerically tractable, we regularize the piecewise constant rotation field, resulting in a diffuse interface model, that preserves the zero-stress character of the grain boundaries. Our main intent here is to introduce the model and its potential, and perform a verification exercise before launching into more ambitious undertakings where grain boundary phenomena can be properly modeled. In the following sections we lay out the essential theoretical developments of our model and provide a verification exercise of the numerical implementation.

Classical crystal plasticity for single crystals

For reference, in this section, we introduce the framework of crystal plasticity for single crystals as a starting point. A body is represented as an open subset $\mathcal {B}$ of the three-dimensional Euclidean space $\mathbb {R}^{3}$. Let $\mathcal {B}^{0} \subset \mathbb {R}^{3}$ represent the reference configuration of the body. The position of an arbitrary material point in the reference configuration is denoted by X. A time-dependent deformation map is given by a one-to-one function y(X,t), such that detF≠0, where

$$\begin{array}{*{20}l} \boldsymbol{F}(\boldsymbol{X},t) := \nabla \boldsymbol{y}(\boldsymbol{X},t) \end{array} $$

(1)

is the gradient of the deformation map. In the theory of crystal plasticity, there exists a decomposition of the deformation gradient given by

$$\begin{array}{*{20}l} \boldsymbol{F} = \boldsymbol{F}^{\mathrm{L}} \boldsymbol{F}^{\mathrm{P}}, \end{array} $$

(2)

where F ^L and F ^P are lattice¹ and plastic components of F respectively, and detF ^P=1. In this paper, F ^P represents the deformation gradient of an infinitesimal material element, attributed to dislocation slip through its volume. Since such a process renders the lattice invariant, it follows that F ^P leaves the lattice undeformed. F ^L represents the deformation of the material due to the deformation of its underlying lattice. Note that F ^L(X,t) and F ^P(X,t)need not be gradients of a deformation map. Instead, since F ^L and F ^P are invertible, they represent deformation of an infinitesimally small neighborhood of X at time t. In other words, F ^PdX represents the deformation of a differential material element d X. The collection of all deformed differential material elements is referred to as lattice configuration. In this sense, F ^P maps the reference configuration to the lattice configuration, and F ^L maps the lattice configuration to the deformed configuration.

As is customary, dislocations move on slip systems α=1,2,…,A, where each α defines a glide direction s ^α and a slip plane normal to m ^α. These two are vectors in the lattice configurations such that

$$\begin{array}{*{20}l} |\boldsymbol{s}^{\alpha}| = |\boldsymbol{m}^{\alpha}| = 1; \quad \boldsymbol{s}^{\alpha} \cdot \boldsymbol{m}^{\alpha} = 0; \quad \boldsymbol{s}^{\alpha}, \boldsymbol{m}^{\alpha} = \text{constant}. \end{array} $$

(3)

Evolution of F ^P is governed by slip rates v ^α(X,t) on individual slip systems via the flow rule

$$\begin{array}{*{20}l} \dot{\boldsymbol{F}}^{\mathrm{P}} = \boldsymbol{L}^{\mathrm{P}} \boldsymbol{F}^{\mathrm{P}}, \end{array} $$

(4)

where

$$\begin{array}{*{20}l} \boldsymbol{L}^{\mathrm{P}}(\boldsymbol{X},t) := \sum_{{\alpha=1}^{A}} v^{\alpha}(\boldsymbol{X},t) \boldsymbol{s}^{\alpha} \otimes \boldsymbol{m}^{\alpha}. \end{array} $$

(5)

If the free energy density, denoted by ψ, depends on the lattice Lagrangian strain

$$\begin{array}{*{20}l} \mathbb{E}^{\mathrm{L}}:=\left(\left(\boldsymbol{F}^{\mathrm{L}}\right)^{\mathrm{T}} \boldsymbol{F}^{\mathrm{L}} - \boldsymbol{I}\right)/2, \end{array} $$

(6)

then the evolution equations of crystal plasticity are given by the flow rule in (4), along with the following macroscopic and microscopic force balance equations:

• Macroscopic force balance

  $$\begin{array}{*{20}l} \text{Div}~\boldsymbol{P}(\boldsymbol{X},t) &= \boldsymbol{0}, \quad \boldsymbol{X} \in \mathcal{B}^{0}, t>0, \end{array} $$

  (7a)
  $$\begin{array}{*{20}l} \boldsymbol{u} &= \boldsymbol{u}^{0}~\text{on}~\partial \mathcal{B}^{0}, \end{array} $$

  (7b)

  where

  $$\begin{array}{*{20}l} \boldsymbol{P} := \boldsymbol{F}^{\mathrm{L}} \psi,_{{\mathbb{E}}^{\mathrm{L}}} \left(\boldsymbol{F}^{\mathrm{P}}\right)^{\mathrm{-T}}, \end{array} $$

  (8)

  is the first Piola–Kirchhoff stress tensor, and $\psi,_{{\mathbb {E}}^{\mathrm {L}}}$ denotes the derivative of ψ with respect to ${\mathbb {E}}^{\mathrm {L}}$.

• Microscopic force balance for each slip system α

  $$\begin{array}{*{20}l} b^{\alpha} v^{\alpha}(\boldsymbol{X},t) = \psi,_{{\mathbb{E}}^{\mathrm{L}}} \boldsymbol{m}^{\alpha} \cdot \boldsymbol{C}^{\mathrm{L}} \boldsymbol{s}^{\alpha}, \end{array} $$

  (9)

  where b ^α≥0 is the inverse of the mobility associated with the slip v ^α, and C ^L=(F ^L)^T F ^L is the right Cauchy-Green strain tensor.

The non-negativity of the inverse mobilities is a necessary condition for thermodynamic consistency. The expression on the right-hand-side of (9) is commonly referred to as the resolved shear stress. See the work by Gurtin for a thermodynamically consistent derivation of (7) and (9) (Gurtin 2000, 2008). In standard crystal plasticity, a stress-free single crystal at t=0 is modeled using the initial conditions

$$\begin{array}{*{20}l} \boldsymbol{F}^{\mathrm{P}}(\boldsymbol{X},0) = \boldsymbol{F}^{\mathrm{L}}(\boldsymbol{X},0) \equiv \boldsymbol{I}. \end{array} $$

(10)

Note that, the above initial conditions are also used for polycrystals, with the difference that L ^P is evolved in a piecewise way in each grain due to the different orientation of the slip systems, and the free energy density given by $\psi \left (\boldsymbol {R}^{\mathrm {T}} {\mathbb {E}}^{\mathrm {L}} \boldsymbol {R}\right)$, where R is a piecewise constant rotation field describing the initial orientation of grains.²

In the next section, we first present a diffuse-interface polycrystal plasticity model which operates at a length scale where all grain boundaries are resolved explicitly. In contrast with assumption (10), the proposed framework gives us access to grain boundary dislocation densities, thus enabling us to model grain boundary energies.

Polycrystal plasticity

Consider a sharp-interface polycrystal, i.e. one where the orientation of the lattice is constant in the interior of one grain and has a jump discontinuity along the grain boundary. In this context, crystal plasticity is studied by having the stress-free polycrystal as the reference configuration. Due to the variation in orientation of the grains, the elastic and plastic response of each grain is different. Therefore, the elastic moduli and the slip systems (s ^α and m ^α) are piecewise constant, with jump discontinuities along the grain boundaries. If the polycrystal is stress-free at t=0, then the initial conditions are identical to (10). Thus, within this framework, polycrystal plasticity is identical to single crystal plasticity with the caveat that the elastic moduli, s ^α and m ^α are piecewise constant. While this model is remarkably simple, it is not straightforward to generalize it to model grain boundary-mediated deformation, such as shear-induced grain boundary motion, grain shrinkage and rotation, grain boundary sliding, etc. These phenomena can become important during plastic deformation at high stresses and/or temperatures, such as during recovery, recrystallization, and grain growth. In the following section, we present an alternate framework that lays the foundation to model polycrystal plasticity with grain boundary evolution.

Diffuse-interface polycrystal plasticity

The success of single crystal plasticity in describing the materials deformation lies in precisely identifying the independent mechanisms involved, and attributing them appropriately to the evolution of F ^P. For example, the rate of change is F ^P due to dislocation slip is identified with the slip rate projected on each slip system by way of the Schmid tensor. Similarly, additional mechanisms such as dislocation climb are built into the evolution law for F ^P (Weertman 1955; Thomson and Balluffi 1962). In addition to dislocations, a grain boundary sweeping through a material also results in plastic distortion. For example, consider a circular grain with lattice orientation θ ₂ embedded in a larger grain with orientation θ ₁. The misorientation of |θ ₂−θ ₁| results in a grain boundary energy. In order minimize the internal energy, the circular grain shrinks. As the circular grain boundary sweeps through the material, the lattice in the swept region rotates from an initial configuration of θ ₁ to θ ₂, while the rest of the lattice remains unchanged. If F ^P is equal to identity during this process, then this results in an incompatible F. This conclusively suggests that F ^P≢I in the swept area. In other words, grain boundary motion always results in plastic distortion.

Therefore, in the spirit of modeling plasticity due to bulk dislocations, plasticity due to grain boundary motion may thus be modeled by identifying the mechanism for the accompanying plastic distortion, and include it in the evolution law for F ^P. Identifying the pertinent GB-mediated plastic mechanisms is highly non-trivial. For example, recent atomistic simulations have revealed that for certain misorientations the interior grain not only shrinks but also rotates with no dislocation activity in the bulk. This suggests that unlike dislocation slip, there is no unique fundamental evolution law for F ^P that can be attributed to the motion of a grain boundary with a given misorientation. Therefore, we take an alternate approach to modeling plasticity due to grain boundary motion.

The central idea behind this approach is to identify dislocations as the basic defect carriers, and build grain boundaries as continuum aggregates of dislocations. Therefore, any motion of grain boundary is viewed as a collective motion of dislocations that form the boundary. The most important advantage of this approach is plastic distortion due to grain boundary motion emerges from the original flow rule given in (4) without identifying any new mechanisms. This approach can model phenomena such as shear-induced grain boundary motion, grain boundary sliding and grain rotation (Admal and Marian 2017). We next build a framework of polycrystal plasticity based on the idea described above.

Let R ⁰(X)∈SO(3), a step function in the space of special orthogonal tensor fields, represent the lattice rotation field in the polycrystal, with piecewise-constant values in each grain and smooth transitions across grain boundaries. In contrast to (10), the initial state of the polycrystal is chosen to be:

$$\begin{array}{*{20}l} \boldsymbol{F}^{\mathrm{L}}(\boldsymbol{X},0) = \boldsymbol{R}^{0}(\boldsymbol{X}), && \boldsymbol{F}^{\mathrm{P}}(\boldsymbol{X},0) = \boldsymbol{R}^{0}(\boldsymbol{X})^{\mathrm{T}}, \end{array} $$

(11)

resulting in

$$\begin{array}{*{20}l} \boldsymbol{F}(\boldsymbol{X},0) \equiv \boldsymbol{I}. \end{array} $$

(12)

The decomposition given in (11) is the central idea of the current framework, and we now describe its physical significance. Figure 1 demonstrates the decomposition given in (11) for the construction of a grain boundary in a bicrystal. Recall that F ^P deforms the material leaving the lattice fixed as shown in Fig. 1. On the other hand, F ^L deforms the lattice resulting in a total deformation gradient F that is compatible. Comparing the reference and the final configurations, Fig. 1 seems contradictory since the material is shown to be deformed although F≡I. We now discuss the correct mathematical interpretation that resolves this contradiction.

Fig. 1

In two-dimensions, the above construction results in exactly two non-zero components (G ₃₁ and G ₃₂) of G. In particular, for a symmetric tilt boundary oriented as shown above, G ₃₂≡0 when θ is a step function

We begin by noting that F ^P(X,0)=R ⁰(X)^T qualifies to be a plastic distortion due to dislocation slip, since a rotation can always be expressed as a product of three shear deformation tensors (Tanaka et al. 1986; Paeth 1986; Toffoli and Quick 1997).³ Interpreting the three resulting shear deformations as lattice-invariant shears obtained due to dislocation slips, the rotation tensor F ^P(X,0) may be interpreted as a lattice-invariant deformation. Since an arbitrary rotation rotates the material, it may seem contradictory for it to leave the lattice invariant (except of course when the rotation belongs to the point group of the lattice). The correct mathematical interpretation of a “lattice-invariant” rotation is given using the notion of weak-convergence discussed in Appendix B. In short, weak convergence represents convergence of functions/distributions on the “average”. In Appendix B, we show that, for a sequence of lattice constants a ⁱ→0 (as i→∞), F ^P(X,0) has to be viewed as a weak-limit of a sequence of deformations (F ^P)ⁱ that leave the a ⁱ-lattice invariant. Therefore, interpreting F ^P(X,0)=R ^T(X) and F≡I for a discrete lattice in an average sense resolves the apparent contradiction described in the previous paragraph.

An important consequence of the decomposition given in (11) is that the resulting polycrystal is stress-free since the Lagrangian strain, defined in (6), is equal to zero. Therefore, Eq. (11) describes a polycrystalline state which is obtained from a reference single crystal by the right amount of slip in each grain such that grains undergo relative rotation but the polycrystal remains stress free.

An advantage of the above construction is that we have immediate access to the grain boundary dislocation density content in the form of the geometrically necessary dislocation density G tensor defined as

$$\begin{array}{*{20}l} \boldsymbol{G} = \boldsymbol{F}^{\mathrm{P}}~\text{Curl}~\boldsymbol{F}^{\mathrm{P}}, \end{array} $$

(13)

where Curl denotes the curl of a tensor field with respect to the material/reference coordinate.⁴ For a given normal n in the lattice configuration, the vector G ^T n measures the net Burgers vector of dislocation lines per unit area passing through a plane of normal n.

Using the decomposition of F discussed above, we can now, in principle, study a polycrystal under a single boundary-value problem. Numerically, the problem still does not enjoy the nice characteristics of its single crystal counterpart as F ^L and F ^P are discontinuous. In order to overcome this challenge, we introduce a smooth-interface version of the above sharp-interface model. This can be achieved by constructing a stress-free diffuse interface crystal plasticity at t=0 with F ^P a smoothened step function in the space of rotation fields. This alteration ensures that all the resulting fields are smooth.

Numerical implementation

In this section, we discuss a three-dimensional numerical implementation of tensile tests of polycrystals of varying textures using the diffuse-interface model introduced in “Diffuse-interface polycrystal plasticity” section. The main aim of this section is to demonstrate the robustness of the diffuse-interface model.

We implement a simpler version of a crystal plasticity model for body-centered cubic (bcc) Fe used by Barton et al. (2013) that incorporates the role of latent hardening into the mobility variable in (9). The microscopic force balance we use in this implementation is given by

$$\begin{array}{*{20}l} v^{\alpha}(\boldsymbol{X},t) = v^{\alpha}_{0} \left (\frac{\tau^{\alpha}}{g^{\alpha}}\right)^{1/m}, \end{array} $$

(14)

where $v^{\alpha }_{0}$ the references shear rate, g ^α(X,t) is the slip system strength that captures the operating hardening mechanism, $\tau ^{\alpha }(\boldsymbol {X},t)=\psi,_{{\mathbb {E}}^{\mathrm {L}}} \boldsymbol {m}^{\alpha } \cdot \boldsymbol {C}^{\mathrm {L}} \boldsymbol {s}^{\alpha }$ is the resolved shear stress, and m=0.05 is the strain-rate sensitivity exponent. The slip strength g ^α depends on the network dislocation density ρ _n (X,t) via Taylor hardening:

$$\begin{array}{*{20}l} g^{\alpha} = g_{0} + b\mu_{0}\sqrt{h_{n} \rho_{n}}, \end{array} $$

(15)

where the constant g ₀=90 MPa refers to the slip strength in a single crystal, μ ₀=86 GPa is the rigidity modulus of iron, and h _n =0.125. The network dislocation density ρ _n in (15) evolves according to the Kocks–Mecking type evolution model (Mecking and Kocks 1981):

$$\begin{array}{*{20}l} \dot \rho_{n} = v \left(k_{1} \sqrt{\rho_{0} \rho_{n}} - k_{2} \rho_{n}\right), \end{array} $$

(16)

with

$$\begin{array}{*{20}l} k_{2} = k_{20} \left(\frac{v_{k0}}{v} \right)^{\frac{1}{n}}. \end{array} $$

(17)

The variable $v = \sum _{\alpha } |v^{\alpha }|$ is the aggregate slip rate, ρ ₀=10¹² m ⁻² is the reference network bulk dislocation density, and the Kocks-Mecking parameters k ₁, k ₂₀, and v _k0 are equal to 450, 14 and 10¹⁰ s ⁻¹ respectively. Finally, the elastic free energy ψ is taken to be of the form:

$$\psi\left({\mathbb{E}}^{\mathrm{L}}\right)=\frac{1}{2}\mathbb{C}~\mathbb{E}^{\mathrm{L}}\cdot{\mathbb{E}}^{\mathrm{L}} $$

where $\mathbb {C}$ is the elasticity matrix, which for a cubic material is fully characterized by three independent elastic constants whose values for Fe are: C ₁₁=228 MPa, C ₁₂=132 MPa, and C ₄₄=116 MPa (Barton et al. 2013).

Finite element implementation

The finite element method is used to solve the resulting system of equations in (4), (7), and (16), with the displacement field u, the plastic distortion F ^P, and the bulk network dislocation density ρ _n as unknowns. In particular, the three displacement variables u ₁, u ₂ and u ₃ are interpolated using the Lagrange quadratic finite elements, while ρ _n is interpolated using the Lagrange linear finite elements. Recall that at t=0, F ^P is a field in SO(3). This implies that it satisfies the condition of orthogonality, i.e. (F ^P)^T F ^P(X,0)≡I. On the other hand, the components of F ^P interpolated using the Lagrange finite elements cannot satisfy the orthogonality constraint in the interior of the finite elements. Therefore, the components of F ^P cannot be interpolated using the Lagrange finite elements. Instead, using the polar decomposition, F ^P is expressed as R ^P U ^P, where R ^P(X,t)∈SO(3), and U ^P(X,t) is the resulting positive-definite plastic stretch tensor. Using the angle-axis representation for rotation tensors, R ^P is expressed in terms of a vector $\boldsymbol {q} \in {\mathbb {R}}^{3}$:

$$\begin{array}{*{20}l} \boldsymbol{R}(\boldsymbol{q}) = \boldsymbol{I} + \frac{\sin |\boldsymbol{q}|}{|\boldsymbol{q}|} \boldsymbol{W} + \frac{1}{2} \left [ \frac{\sin(|\boldsymbol{q}|/2)}{(|\boldsymbol{q}|/2)} \right ]^{2} \boldsymbol{W}^{2}, \end{array} $$

(18)

where W is the skew-symmetric matrix associated with q, and |q| and q/|q| represent the angle and axis of the rotation tensor. Lagrange linear finite element interpolation is then chosen for the variables U ^P and q, from which F ^P is locally computed as F ^P=R ^P(q)U ^P. This method guarantees that F ^P(X,0) can describe an exact plastic rotation field without numerical artifacts due to the interpolation method.

To simulate tensile tests of polycrystals with different textures, we impose the boundary conditions shown in Fig. 2. The initial conditions for u, U ^P and ρ _n are chosen to be

$$\begin{array}{*{20}l} \boldsymbol{u}(\partial \mathcal{B}^{0}) = \boldsymbol{0}, \quad \boldsymbol{U}^{\mathrm{P}}(\mathcal{B}^{0}) \equiv \boldsymbol{I},\text{ and}~\quad \rho_{n}(\mathcal{B}^{0}) \equiv 20 \rho_{0}, \end{array} $$

(19)

Fig. 2

A schematic of the geometry of the polycrystal domain and the imposed boundary conditions. The length of the cubic domain L=3 micrometers, and the imposed strain rate $\dot \epsilon _{11}=100 s^{-1}$

respectively. The remaining initial conditions for q, which defines the texture of the polycrystal, is discussed in “Construction of polycrystals with different textures” section.

The system of Eqs. (4), (7), and (16) are evolved in a segregated manner using the MUMPS direct solver, and BDF (Backward Differential Formula) time stepping algorithm implemented in COMSOL5.2. In particular, due to the highly nonlinear nature of (4) expressed in q and U, we enable the “automatic highly nonlinear (Newton)” option to obtain well-behaved solutions. On the other hand, we rely on the default “constant (Newton)” option for solving (7), and (16). All simulations were performed on a finite element mesh with 99883 elements, and 595620 degrees of freedom.

Construction of polycrystals with different textures

In this section we describe the generation of diffuse interface polycrystals of different textures. The grain orientations are outputted in the form of a smoothened rotation vector field q(X,0) which serves as an initial condition along with those given in (19).

A stress-free polycrystal with N grains is constructed by randomly choosing N points, $\boldsymbol {\mathcal {P}}^{1},\dots, \boldsymbol {\mathcal {P}}^{N}$, within the domain, and constructing a corresponding diffuse Voronoi tessellation. The grain orientations are prescribed by associating random rotation vectors q ¹,…,q ^N to each grain. The diffuse tessellation is constructed using a grid of size 100×100×100, and assigning each grid point to a grain based on the Voronoi construction, i.e. a grid point p _i is associated with a grain α if

$$\begin{array}{*{20}l} \text{dist}\left(\boldsymbol{p}^{i},\boldsymbol{\mathcal{P}}^{\alpha}\right) < \text{dist}\left(\boldsymbol{p}^{i}, \boldsymbol{\mathcal{P}}^{\beta}\right), \quad \forall \beta \ne \alpha, \end{array} $$

(20)

where $\text {dist}(\boldsymbol {p}^{i},\boldsymbol {\mathcal {P}}^{\beta })$ is the distance between p ⁱ and $\boldsymbol {\mathcal {P}}^{\beta }$. Finally, the rotation vector q ^α is associated to the grid point p ⁱ. The polycrystal is outputted in the form of the rotation vector field on the grid which is then interpolated as a smooth vector field q(X,0) using the nearest neighbor interpolation implemented in COMSOL5.2. Therefore, the texture of the resulting collection of grains depends on the distribution of the initial collection of N random points, and the grain boundary “thickness” is inversely proportional to the resolution of the grid. The pseudocode for the above algorithm is described in Algorithm 1.

We study textures with (i) a log normal distribution of grain sizes⁵, (ii) elongated grains, and (iii) flat grains. The size of a grain α is defined as

$$\begin{array}{*{20}l} \text{size}(\alpha) = \min_{\beta}\left\{\text{dist}\left(\boldsymbol{\mathcal P}^{\alpha},\boldsymbol{\mathcal{P}}^{\beta}\right): \beta \in \{1,\dots,N\}, \beta \ne \alpha\right\}. \end{array} $$

(21)

Grains with a log normal distribution of sizes are generated by sampling the initial N points from a log normal distribution based on Algorithm 1 described in Appendix A.

We use the standard Euclidean metric for dist in (20) in the generation of the texture with log normal distribution of grain sizes. On the other hand, elongated and flat grains with aspect ratio equal to 4 are generated by sampling the initial N points from a Dirac probability measure supported on 0.2L, and using the metric

$$\begin{array}{*{20}l} \text{dist}(\boldsymbol{x},\boldsymbol{y}) &= \left(\frac{x_{1}-y_{1}}{sx}\right)^{2}+ \left(\frac{x_{2}-y_{2}}{sy}\right)^{2}+ \left(\frac{x_{3}-y_{3}}{sz}\right)^{2}, \end{array} $$

(22)

with the scales sx=1, sy=1, sz=4 for elongated grains, and sx=0.25, sy=0.25, and sz=1 for flat grains. Polycrystals with the three textures studied in this paper are shown in Fig. 3, with the colors obtained by plotting the q ₃ component, indicating different grains.

Fig. 3

Polycrystals with different textures with a log normal grain size distribution with parameters σ=0.15 and μ=− log105, (b) and c flat and elongated grains both with an aspect ratio of 4

Results

In this section, we present our results of the simulated tensile tests on polycrystals of varying textures. Figure 4 shows a color plot of the grain boundary dislocation density for a polycrystal with log normal grain size distribution, calculated using (13). Note that the field G is not available in the classical polycrystal plasticity implementation. In the proposed model, the initialization (11) allows to construct a kinematically consistent grain boundary structure which evolves in time as a consequence of slip in each grain.

Fig. 4

A color plot of the norm of the geometrically necessary dislocation density tensor G:=F ^PCurlF ^P expressed in units of m ⁻¹

Figure 5 shows plots of the intermediate stress $\psi,_{\mathbb {E}}^{\mathrm {L}}$ versus the first axial component of the total strain E:=(F ^T F−I)/2 for different textures and loading orientations. In addition, Fig. 5 also shows the variation of the normalized dislocation density h with respect to the total strain. We have verified that the plots shown in Fig. 5 are insensitive to further mesh refinement. In addition, since the computed properties are aggregates, as expected, we ensured that the results are not sensitive to grain boundary thickness. We expect that local properties such as stress concentration will be sensitive to the choice of grain boundary thickness. We also compute the Taylor factor, which is known to be 2.9 for an equiaxed bcc random polycrystal Kocks (1970b).

Fig. 5

Plots of the intermediate stress $\psi _{,{\mathbb {E}}^{\mathrm {L}}_{11}}$ and the normalized dislocation density ρ _n /ρ ₀ versus the first axial component of total strain E:=(F ^T F−I)/2 for polycrystals of different textures and loading orientations

The Taylor factor M is defined as the ratio of the aggregate microscopic shear rate in a polycrystal to the macroscopic shear rate. It is defined using the following equivalence of the power supplied by external loads to the power dissipated due to slip:

$$\begin{array}{*{20}l} \boldsymbol{P} \cdot \dot{\boldsymbol{F}} = \sum_{\alpha} \tau^{\alpha} v^{\alpha}. \end{array} $$

(23)

Assuming there exists a constant critical resolved shear stress τ ^c>0 for every slip system at the which a crystal slips, (23) can be simplified to

$$\begin{array}{*{20}l} \boldsymbol{P} \cdot \dot{\boldsymbol{F}} = \tau^{\mathrm{c}} \sum_{\alpha} |v^{\alpha}|. \end{array} $$

(24)

The Taylor factor M is defined as

$$\begin{array}{*{20}l} M &:= \left \langle \frac{\boldsymbol{P} \cdot \dot{\boldsymbol{F}}}{\tau^{\mathrm{c}} |\dot{\boldsymbol{F}}|} \right \rangle, \end{array} $$

(25a)

$$\begin{array}{*{20}l} &= \left \langle \frac{\sum_{\alpha} |v^{\alpha}|}{|\dot{\boldsymbol{F}}|} \right \rangle, \end{array} $$

(25b)

where we have used (24) to arrive at the last equality, and 〈·〉 denotes spatial average. The traditional definition of the Taylor factor given in (25) cannot be used in a straightforward manner in our implementation since the slip does not occur precisely at a critical load. In fact, when implemented, (25a) and (25b) neither agree, nor converge with time. On the other hand, by factoring out $\langle \sum _{\alpha } \tau ^{\alpha }\rangle$ instead of τ ^c in (23), we show that the following two definitions for M given by

$$\begin{array}{*{20}l} M &:= \left \langle \frac{\boldsymbol{P} \cdot \dot{\boldsymbol{F}}}{\left(\sum_{\alpha} \tau^{\alpha}\right) |\dot{\boldsymbol{F}}|} \right \rangle, \end{array} $$

(26a)

$$\begin{array}{*{20}l} &= \left \langle \frac{\sum_{\alpha} \tau^{\alpha} |v^{\alpha}|}{\left(\sum_{\alpha} \tau^{\alpha}\right)|\dot{\boldsymbol{F}}|} \right \rangle \end{array} $$

(26b)

are not only consistent with each other, but also converge to a constant value as shown in Fig. 6. The converged values of the Taylor factors for different textures are listed in Table 1.

Fig. 6

Plots of the variation of the Taylor factor computed using (26) for polycrystals of different textures and loading orientations

Table 1 Taylor factors for polycrystals of different textures

Discussion and conclusions

Most metals and alloys in usable form display an internal microstructure characterized by a collection of grains with different lattice orientation separated by grain boundaries. Metals deformation, particularly at high temperatures and stresses, such as during hot working, involves not just intragranular plasticity but also plasticity controlled by grain boundary mechanisms. Standard formulations of crystal plasticity decouple both types of deformation, probably due to our good deal of understanding about low temperature processes, e.g. cold working, which tends to dominate our thinking of plasticity. Indeed, this decoupling has been the governing principle behind the development of new methodologies to study recrystallization in metals (Singer-Loginova and Singer 2008; Steinbach 2009; Takaki and Tomita 2010; Abrivard et al. 2012; Kamachali 2013).

However, during dynamic phenomena such as continuous dynamic recrystallization, bulk and grain-boundary plastic processes can occur simultaneously, and therefore the underlying plasticity model must be capable of capturing both types of deformations concurrently. This is the motivation behind the present work: to devise a computational model that combines bulk and grain boundary plasticity by design within the same framework. Our purpose at the moment is simply to demonstrate that our formulation is capable of rendering the same response as standard crystal plasticity models for conventional problems in polycrystal plasticity. Only after fulfilling this step can we truly apply our methodology to phenomena involving grain boundary processes. We have undertaken this verification exercise by solving the same problem, standard Taylor hardening in body-centered cubic Fe, using both methodologies, and comparing the results obtained. To explore the capabilities of our model further, we have considered several different textures and misorientation ranges and have calculated the associated Taylor factors. In all cases, our results agree with those obtained using standard polycrystal plasticity.

In summary, we have developed a diffuse-interface model for polycrystalline materials deformation that expresses grain boundaries as a special class of geometrically necessary dislocations, such that the stress-free nature of the polycrystalline structures obtained is naturally recovered. We have tested the robustness of the method by simulating tensile tests and calculating Taylor factors for polycrystals of varying textures. Our model provides a pathway from which grain boundary energies and mobilities can eventually be obtained directly from dislocation densities, which opens the door to integrated models of intragranular and grain boundary-governed plasticity such as recrystallization in hot working.

Endnotes

¹ In the literature, it is more common to refer to the lattice distortion as an elastic distortion using the notation F ^E. Since the decomposition given in (2) is purely geometric in nature (as opposed to energetic), we prefer the term “lattice distortion” denoted by F ^L, a terminology adopted by Clayton (2010).

² It is important to note that, within the framework of crystal plasticity, a constitutive response function of the form $\psi (\boldsymbol {R}^{\mathrm {T}} {\mathbb {E}}^{\mathrm {L}} \boldsymbol {R})$ with a non-constant R does not imply F ^L=R and F ^P≡I as this would result in an incompatible F.

³ For example, a rotation by angle θ about the z-axis can be decomposed multiplicatively into three shear deformations as

$$\begin{array}{*{20}l} \left[ \begin{array}{ll} \cos \theta & -\sin\theta \\ \sin \theta & \cos \theta \end{array}\right] = \left[ \begin{array}{lc} 1 & -\tan\left(\frac{\theta}{2}\right) \\ 0 & 1 \end{array}\right] \left[ \begin{array}{cc} 1 & 0\\ \sin\theta & 1 \end{array}\right] \left[ \begin{array}{cc} 1 & -\tan\left(\frac{\theta}{2}\right) \\ 0 & 1 \end{array}\right]. \end{array} $$

⁴ The curl of a tensor field T is defined as

$$\begin{array}{*{20}l} (\text{Curl}~T)\boldsymbol{n} := \text{Curl}~(\boldsymbol{T}^{\mathrm{T}} \boldsymbol{n}), \notag \end{array} $$

where n is an arbitrary constant vector, and the curl on the right-hand-side of the above equation is the curl of a vector field defined as (v) _i =ε _ijk v _j,k , for any vector field v.

⁵ The distribution of a random variable whose logarithm is distributed normally is called a log normal distribution. The cumulative distribution function of a log normal random variable with parameters σ and μ is given by

$$\begin{array}{*{20}l} \Phi \left (\frac{\ln x-\mu}{\sigma} \right), \end{array} $$

(27)

where Φ is the cumulative distribution function of the standard normal distribution.

Appendix A: Algorithm to generate polycrystals with different textures

In this section, we describe the algorithm used to generate the different polycrystal textures simulated in this paper. Algorithm 1 is able to generate a polycrystal with a given cumulative distribution function f for grain sizes. In addition, grains of desired aspect ratio can be generated using the scales sx, sy and sz as given in Algorithm 1.

The variable maxiter has to be set by trial and error until a satisfactory distribution of grain size is obtained relative to the distribution f. A very high or a low value skews the resulting distribution away from f. Intuitively, increasing maxiter increases the number of tries to pack more grains such that the distribution of grain size is consistent with the given distribution. But as the number of grains increases, the correlation between sizes of adjacent grains increases resulting in a distribution away from the desired distribution. For example, a value of maxiter =1000 is used to generate the texture shown in Fig. 3. From Fig. 7, which compares the texture’s grain size distribution resulting from Algorithm 1 to a randomly generated log normal distribution of numbers, we conclude that the two distributions are reasonably close.

Fig. 7

A comparison of the histogram plots of the probability density functions of a log-normal distribution, and a grain size distribution resulting from Algorithm ?? with the parameter maxiter =1000

Appendix B: Interpretation of F ^P=R ^T and F=I using the notion of weak convergence

In this section, we use the notion of weak-convergence to arrive at a physical interpretation of the decomposition given in (11), and depicted in Fig. 1 for a discrete lattice. Recall the apparent contradiction we arrive at by interpreting F ^P=R ^T in an absolute sense for a discrete lattice. On the one hand, F ^P=R ^T should be a lattice-invariant deformation, while on the other hand an arbitrary rotation need not preserve the lattice. We will now show that, for a discrete lattice, F ^P=R ^T and F=I should be viewed in an average sense using the notion of weak convergence Rudin (2006).

Definition 1

A sequence of distributions Λ ⁱ converges weakly to a distribution Λ if

$$\begin{array}{*{20}l} {\lim}_{i \to \infty} \Lambda_{i}(\phi) = \Lambda \phi \end{array} $$

(28)

for all ϕ in the space of smooth functions with compact support, denoted by $C_{\mathrm {c}}^{\infty }$.

Given a constant rotation R, we will now construct a sequence of deformations (F ^P)ⁱ that converge weakly to R. Each (F ^P)ⁱ leaves a lattice with lattice constant a ⁱ unchanged, and a ⁱ→0 as i→∞. In other words, (F ^P)ⁱ converges to R ^T on an “average” as the lattice constant tends to zero. Assuming a square lattice, $\left (\boldsymbol {F}^{\mathrm {P}}\right)^{i}(\boldsymbol {X}) := \nabla \tilde {\boldsymbol {x}}^{i} (\boldsymbol {X})$, where

$$\begin{array}{*{20}l} \tilde{\boldsymbol{x}}^{i}(\boldsymbol{X}) = \left \lfloor \frac{\boldsymbol{R}^{\mathrm{T}} \boldsymbol{X}}{a^{i}} \right \rfloor a^{i}, \end{array} $$

(29)

and ⌊·⌋ denotes the floor function. The deformation given by (29) ensures that the lattice remains unchanged. Note that (F ^P)ⁱ should be viewed as a distribution since $\tilde {\boldsymbol {u}}^{i}$ is a piecewise constant vector field. It can be easily shown that $\tilde {\boldsymbol {x}}^{i}(\boldsymbol {X})$ uniformly converges to R ^T X as the lattice constant a ⁱ→0. On the other hand, (F ^P)ⁱ does not converge, pointwise or uniformly, to R ^T. Instead, it converges weakly to R ^T. This can be easily demonstrated using the divergence theorem. For an arbitrary $\phi \in C_{\mathrm {c}}^{\infty }$, we have

$$\begin{array}{*{20}l} {\lim}_{i\to 0} \int_{\Omega} (\boldsymbol{F}^{\mathrm{P}})^{i} \phi \, d\boldsymbol{X} &= -{\lim}_{i\to 0} \int_{\Omega} \tilde{\boldsymbol{x}}^{i} \otimes \nabla \phi \, d\boldsymbol{X} \notag \\ &= -\int_{\Omega} \boldsymbol{R}^{\mathrm{T}} \boldsymbol{X} \otimes \nabla \phi \, dX \notag \\ &= \int_{\Omega} \boldsymbol{R}^{\mathrm{T}} \phi \, d\boldsymbol{X}, \end{array} $$

(30)

where we have used the divergence theorem along with ϕ=0 on ∂ Ω to arrive at the first and last equalities, and the uniform convergence of $\tilde {\boldsymbol {x}}^{i}$ to interchange the limit and the integral signs in the first equality. By the definition of weak convergence, (30) implies (F ^P)ⁱ→R ^T weakly. Assuming F ^L=R, it can be similarly shown that the sequence F ⁱ:=F ^L(F ^P)ⁱ converges weakly to the identity.