Gradient Crystal Plasticity: A Grain Boundary Model for Slip Transmission

Interaction between dislocations and grain boundaries (GBs) in the forms of dislocation absorption, emission, and slip transmission at GBs significantly affects size-dependent plasticity in fine-grained polycrystals. Thus, it is vital to consider those GB mechanisms in continuum plasticity theories. In the present paper, a new GB model is proposed by considering slip transmission at GBs within the framework of gradient polycrystal plasticity. The GB model consists of the GB kinematic relations and governing equations for slip transmission, by which the influence of geometric factors including the misorientation between the incoming and outgoing slip systems and GB orientation, GB defects, and stress state at GBs are captured. The model is numerically implemented to study a benchmark problem of a bicrystal thin film under plane constrained shear. It is found that GB parameters, grain size, grain misorientation, and GB orientation significantly affect slip transmission and plastic behaviors in fine-grained polycrystals. Model prediction qualitatively agrees with experimental observations and results of discrete dislocation dynamics simulations.


Introduction
Fine-grained polycrystalline materials with grain sizes ranging from hundreds of nanometers to tens of microns can be widely used as structural or functional materials in small-scaled engineering such as microdevices and microelectromechanical systems. Both their precision manufacturing and processing and their safety and reliability in practical applications entail a deep understanding of the plastic behavior in fine-grained polycrystals. As found in various experiments [1][2][3][4][5][6], if the average grain size decreases to the micrometer range or even smaller, the plastic behavior in crystals such as macroscopic yield strength and strain hardening rate become strongly size-dependent. Microscopically, plastic deformation in crystals mainly originates from the collective behavior of a vast number of dislocations, implying that size-dependent plasticity is closely related to the distinctive dislocation activities in small-scaled crystals. Since the grain boundaries (GBs) to volume ratio in fine-grained polycrystals is extremely high, interaction between dislocations and GBs may be an important, and in some situations, dominant deformation mechanism therein. From numerous in situ experiments as reviewed in [7] and atomistic simulations as reviewed in [8], dislocations interact with GBs in several different ways. Due to the crystallographic incompatibility between adjacent grains, GBs act as obstacles, impeding dislocations nucleated within the grain interior, which results in dislocation pile-ups at GBs [9]. In response to stress concentrations resulting from dislocation pile-ups, GBs may act as dislocation sinks or sources, absorbing the leading dislocation in the pile-up or emitting new dislocations onto the slip planes in the adjacent grains [9][10][11]. Under certain conditions, dislocations are transferred through GBs in a direct or indirect manner with the generation of residual GB defects for the conservation of Burgers vectors, which is also known as slip transmission [9,[12][13][14]. Various factors such as the initial GB structure, the GB orientation, the grain misorientation, and the stress state near the GB affect dislocation-GB interactions [7,15]. As illustrated by experiments [16,17] and discrete dislocation dynamics (DDD) simulations [18][19][20], dislocation-GB interactions can significantly affect macroscopic plastic behaviors in fine-grained polycrystals. Moreover, due to the GB constraint, plastic deformation in fine-grained polycrystals is strongly heterogeneous, which results in the accumulation of geometrically necessary dislocations (GNDs) stored to accommodate lattice curvature [21]. There exists a close connection between GNDs and size-dependent plasticity in crystals [22]. Conventional plasticity theories fail to capture dislocation-GB interactions and the influence of GNDs and are incapable of predicting size-dependent plasticity in fine-grained polycrystals. Thus, it is vital to properly model the above microscopic dislocation-relevant mechanisms, especially dislocation-GB interaction mechanisms at the continuum level, so as to construct non-classical plasticity theories aimed at the accurate characterization of size-dependent plasticity in fine-grained polycrystals.
In the last three decades or so, many efforts have been devoted to develop various strain gradient plasticity theories with consideration of the influence of plastic strain gradient or GNDs, which are capable of predicting some experimentally observed size effects [23][24][25][26][27]. In higher-order strain gradient plasticity theories [28][29][30][31][32][33][34][35], in addition to the traditional force balance equations, microscopic governing equations involving higher-order stresses/microstresses or back stresses resulting from interactions between GNDs are introduced. Consequently, additional microscopic boundary conditions at GBs/surfaces are required for the completeness of the mathematical framework. This provides the possibility to incorporate the GB/surface-dislocation interactions into those continuum theories. Within the framework of high-order strain gradient plasticity theories, some GB/surface models have been proposed. Two frequently used idealized GB/surface conditions are the microhard or microclamped model (corresponding to impenetrable GBs/surfaces) and the microfree model (corresponding to freely penetrable GB/surfaces). Additionally, some attention has been paid to develop intermediate GB/surface models for realistic GBs/surfaces with finite resistance against dislocation gliding. From the viewpoint of thermodynamics, the GB/surface-dislocation interaction is of both energetic and dissipative nature since emission, absorption, and transmission of dislocations may lead to the storage of residual defects/surface steps at GBs/surfaces, and such processes are resistive. Thus, the existing intermediate GB/surface models can be classified into those of dissipative ones [36][37][38][39], energetic ones [40][41][42][43][44][45], and both energetic and dissipative ones [25,[46][47][48][49][50][51][52]. The majority of these existing surface/GB models are phenomenological. A few physically based GB/surface models [46][47][48]52] capture GB/surface effects to some extent.
Due the complication of interaction between dislocations and GBs, further modeling efforts are still needed. In fact, different forms of dislocation-GB interaction such as GB absorbing, emitting, and transmitting dislocations may affect the macroscopic plastic behaviors in fine-grained polycrystals in different ways [53]. Therefore, it is necessary to distinguish these different GB mechanisms in GB models.
Based on the above background, the aim of the present work is to model slip transmission at the GB by considering the underlying physical mechanism. Within the framework of finite deformation gradient polycrystal plasticity, a new GB model for slip transmission is proposed. The GB model consisting of the GB kinematic relations and microscopic force balance equations for slip transmission captures the influence of geometric factors including the misorientation between the incoming and outgoing slip systems and GB orientation, GB defects, and stress state at GBs. The model is applied to study the plastic behavior of a bicrystal thin film under plane constrained shear. Results predicted by the model qualitatively agrees with those from experimental observations and DDD simulations.

Bulk Kinematics: Finite Deformation Gradient Polycrystal Plasticity
Following [54], in finite deformation plasticity, the deformation gradient tensor F is assumed to be multiplicatively decomposed into an elastic part F e and a plastic part F p , Consider that plastic deformation induces no change to volume, one obtains detF p = 1 with det being the determinant operator. In the following, if necessary, quantities in the intermediate and those in the current configuration will be identified by the subscripts i and c, respectively. The velocity gradient tensor L is defined as where [∇a] ij = a i,j represents the gradient of a vector a, u is the displacement vector, the superposed dot denotes time derivative, and the superscript −1 denotes the inverse of a tensor. The elastic velocity gradient L e and the plastic velocity gradient L p are expressed as Studying polycrystals, in each grain, a set of slip systems is introduced. For each slip system α, a unit vector of slip direction s (α) and a unit vector of slip plane normal m (α) attached to the lattice space in the intermediate configuration are defined. s (α) and m (α) are assumed to be unaltered from the reference configuration to the intermediate configuration. Following [55], the plastic velocity gradient consisting of contributions from all active slip systems reads withγ (α) being the plastic slip rate of the slip system α.
The involved strain measures are the Green-Lagrange strain tensor E and the right Cauchy-Green stretch tensor C, and their elastic parts E e and C e , which are expressed as The superscript T denotes the transpose of a tensor, and I is the second order identity tensor. For the slip system α, the evolutions of the edge and screw GND densities (The GND density characterizes the net density of polarized dislocations. Once the slip rate gradient is given, the edge and screw GND densities can be evaluated by Equations (8) and (9), respectively). ρ where ∇ i (·) = F p−T · ∇ r (·), b (α) is the magnitude of Burgers vector, and p (α) = s (α) × m (α) . In addition, the initial conditions are ρ and ρ , with ρ ge(α) 0i and ρ gs(α) 0i being the initial GND densities.

Kinematic Relations at GBs: Slip Transmission
As illustrated in Figure 1, slip transmission occurs in a way that dislocations on the incoming slip systems in grain A are transferred through the GB Ω AB onto the corresponding outgoing slip systems in grain B, during which the transmission of each dislocation produces a residual GB defect for the conservation of Burgers vector. In the following, subscripts A and B are used to identify quantities in grains A and B respectively. The outward normal unit vectors at the GB of grains A and B are denoted as N GB A and N GB B in the intermediate configuration.
In continuum crystal plasticity, it is assumed that slip planes and dislocations are continuously distributed. Thus, we consider the slip transmission as a point-level process at GBs. Generally, a specified incoming slip system may correspond to several different outgoing slip systems, and vice versa. Thus, for an incoming slip system α, the total number of incoming dislocations per unit time n

(α)
A transmitted through a point x GB on the GB Ω AB is expressed as where n in(αβ) AB denotes the number of incoming dislocations per unit time during the slip transmission process between the incoming slip system α and the outgoing slip system β. Accordingly, for an outgoing slip system β the total number of outgoing dislocations per unit time n where n out(αβ) AB is the number of outgoing dislocations per unit time during the slip transmission process between the incoming slip system α and the outgoing slip system β. As illustrated in Figure 1, for a specified slip transmission process between the incoming slip system α and the outgoing slip system β, each incoming dislocation may result in one outgoing dislocation. Therefore, the number of incoming dislocations per unit time equals the number of the outgoing dislocations per unit time, i.e., Given the number of dislocations per unit time n (α) gliding through a point x and the average distance between slip planes L (α) , the plastic slip rateγ (α) for a slip system α at that point can be expressed aṡ From Equations (10)-(13), the plastic slip ratesγ B for the incoming slip system α and the outgoing slip system β at the GB point x GB are written aṡ γ (α) withγ (αβ) B contributed by the slip transmission process between the incoming slip system α and the outgoing slip system β. Equation (14) represents the GB kinematic relation in the present model. By taking advantage of the balance of Burgers vector at the GB [46] and the GB kinematic relation in Equation (14), the rate of the density of residual GB defects G

(αβ)
ABi contributed by the slip transmission between the incoming slip system α and the outgoing slip system β in the intermediate configuration can be expressed asĠ

Balance Equations and Boundary Conditions
In the present work, the work-conjugate gradient crystal plasticity framework is adopted, and, hence, the balance equations and boundary conditions are derived via the principle of virtual power. Without loss of generality, a bicrystal domain consisting of grains A and B separated by the GB Ω AB is considered. The total volume and the total surface area of the domain are V i = V Ai V Bi and S i = S Ai S Bi . Following [31], in the intermediate configuration, the second Piola-Kirchhof stress tensor S e power-conjugate to the rate of the Green-Lagrange strain tensorĖ e , the microforce π (α) i power-conjugate to the slip rateγ (α) and the microstress ξ (α) i power-conjugate to the slip rate gradient ∇ iγ (α) for each slip system α contribute to the internal power in the grain interiors. Following [52], the contribution to the internal power by the surface microforce η S(α) i power-conjugate to the slip rateγ at the surface is considered to account for the influence of dislocation absorption by surfaces. In addition, to incorporate the influence of slip transmission between the incoming slip system α and the outgoing slip system β into the present theory, the contribution by the GB microforce η AB at the GB is introduced. It is assumed that no GB sliding or opening occurs.
In other words, displacements are continuous at the GB. Therefore, there is no macroscopic power contribution at the GB. Thus, the total internal power P int in the intermediate configuration is where the second Piola-Kirchholf stress tensor S e is defined as with σ being the Cauchy stress tensor, and J = detF e . In the absence of body forces, the traditional traction T i power-conjugate to the velocityu at the surface and the microtraction Ξ i power-conjugate to the plastic slip rateγ (α) at the surface contribute to the external power P ext , which gives According to the principle of virtual power, the variation of the internal power with respect to the velocityu and the plastic slip rateγ (α) equals that of the external power, resulting in Based on Equations (2), (3), (5) and (7), the variation ofĖ e is expressed as Substituting Equation (22) into Equation (21) and using the GB kinematic relation (14) and the divergence theorem, one can rewrite Equation (21) as where M = C e · S e is the Mandel stress tensor, and N S i is the surface outward normal unit vector. Considering that Equation (23) should be satisfied for arbitrary δu, one obtains the following traditional balance of momentum and the corresponding standard traction condition in the intermediate configuration, Div i (F e · S e ) = 0 (24) and at the part of surface S T i with prescribed traction. Further, given the validity of Equation (23) for arbitrary δγ (α) and δγ (αβ) ABi , microscopic balance equations in the bulk and the associated microscopic boundary conditions at the surface/GB are expressed as with τ being the Schmid stress, and As indicated in [52], Equation (27) is regarded as the governing equation for dislocation absorption by surfaces. In the present GB model, Equation (28) acts as the governing equation for the slip transmission process between the incoming slip system α in grain A and the outgoing slip system β in grain B at the GB, in which the normal components of the microstresses ξ (α) Ai · N GB Ai and ξ (β) Bi · N GB Bi from both grains are treated as the driving force for slip transmission.
In some existing GB models [46,47], each slip system possesses an independent microscopic force balance equation at the GB, and, hence, the correlation between plastic slips at one side of the GB and those at the other side is not directly reflected. However, during slip transmission, the plastic slip for the incoming slip system and that for the outgoing slip system should be coupled, which is captured in the present work. In the present GB model, an incoming slip system α in grain A and the corresponding outgoing slip system β in grain B involved in a specified slip transmission process share a GB governing equation, and the GB kinematic relation (14) serves as an additional boundary condition. In addition, since microstresses from both grains are involved in the GB governing Equation (28), the intergranular interaction between dislocations from the two adjacent grains at the GB is naturally considered.

Bulk Constitutive Relations: Gradient Crystal Plasticity
The elastic behavior of the crystal is captured by a compressible neo-Hookean material model, and, hence the hyperelastic constitutive relation for the second Piola-Kirchhof stress S e is given by where J = detF e = √ detC e , λ = 2νµ/[1 − 2ν], µ is the shear modulus, and ν is Poisson's ratio. The microstress ξ (α) i generally consisting of an energetic part and a dissipative part [56][57][58] is assumed to be purely energetic for the sake of simplicity. The constitutive relation for ξ (α) i based on the elastic interaction between GNDs within the same slip system and those from different slip systems suggested in [59] is adopted here, i.e., with where l denotes the radius of the domain within which the interaction between GNDs is considered. A visco-plastic power-law relation is adopted for the dissipative microforce π where sgn() is the sign function, m b is the rate-sensitivity exponent in the bulk, andγ b 0 is the bulk reference slip rate. The slip resistance R b(α) i resulting from the interaction between statistically stored dislocations (SSDs) has the following linear hardening form, where γ is the initial slip resistance, and H b(αβ) is the local hardening modulus. Similar to some existing works [52,58], in the present model, the evolution of the SSD density which quantifies the total non-polarized dislocation density is not explicitly considered. Instead, the influence of SSDs is reflected by the dependence of the slip resistance R b(α) i on the accumulated plastic slip which acts as a measure of the SSD density.
Inserting Equation (32) into Equation (26), one can rewrite the microscopic balance equation aṡ which acts as the evolution equation for the plastic slip in the bulk.

Surface Constitutive Relation: Dislocation Absorption by Surfaces
The surface model with consideration of the influence of dislocation absorption by surfaces in [52] is summarized here. In order to consider energetic and dissipative surface effects respectively, the surface microforce η which is derived by considering the change of the surface energy due to the formation of surface steps after dislocation absorption by surfaces. In Equation (35), Γ 1 denotes the surface free energy per unit surface step area, and Γ (αβ) 2 is the energetic surface hardening modulus measuring the interaction strength between surface steps from slip systems α and β. The dissipative microforce η Sdis(α) i accounting for the dissipative resistance against dislocation absorption by surfaces has a visco-plastic power-law form, whereγ S 0 is the reference slip rate at the surface, and m S is the surface rate-sensitivity exponent. The surface slip resistance R S(α) i is expressed as where γ S(β) acc = γ (β) dt is the accumulated slip at the surface, R S(α) 0 is the initial surface slip resistance depending on the initial density of surface defects, and H S(αβ) is the dissipative surface hardening modulus. The surface parameters may depend on the initial surface state such as surface coatings, oxide layers or initial surface defects.
The higher-order traction Ξ (α) i is ignored in the following. Then, by taking advantage of Equation (36), the microscopic surface boundary condition Equation (27) is rewritten aṡ which is the governing equation for dislocation absorption by surfaces.

GB Constitutive Relation: Slip Transmission
The key novelty of the present work is the consideration of slip transmission. To capture energetic and dissipative GB effects due to slip transmission, the GB microforce η To consider the energetic GB effect, it is assumed that the power expended by η GBen(αβ) ABi at the GB equals the increase rate of GB energyψ GB i induced by the accumulation of GB defects resulting from slip transmission, namely A simple quadratic form with respect to GB Burgers tensor G AB is adopted for the GB energy density such thatψ GB i is expressed aṡ where is the energetic GB hardening modulus. In Equation (41), the interaction between GB defects from different slip transmission processes are considered. Combing Equation (40) with Equation (41), one obtains the following constitutive relation for the energetic GB microforce η The dissipative GB microforce η GBdis(αβ) ABi captures the resistance against slip transmission and is assumed to possess the following visco-plastic power-law relation whereγ GB 0 is the reference slip rate at the GB, and m GB is the GB rate-sensitivity exponent. The GB resistance R GB(αβ) AB is expressed as where G is the initial GB resistance, and H

GB(αβϕω) AB
is the dissipative GB hardening modulus. In Equation (44), the interaction between different slip transmission processes is considered.
By substituting Equation (43) into Equation (28), the microscopic GB force balance Equation (28) is rewritten asγ which is the governing equation for the slip transmission process between the incoming slip system α in grain A and the outgoing slip system β in grain B at the GB, whereγ AB is regarded as the measure of the rate of slip transmission. The present GB model captures the influence of important factors including the misorientation between the incoming and outgoing slip systems, the GB orientation, GB defects and stress states at the GB which may affect slip transmission.

Numerical Example: A Bicrystal Thin Film Under Plane Constrained Shear
To illustrate influences of GB effects due to slip transmission on plastic behaviors in fine-grained polycrystals, the theory is implemented to study a plane strain benchmark problem of plane constrained shear of a bicrystal thin film. As depicted in Figure 2, the thin film consists of two grains A and B separated by a tilt GB parallel to the surfaces. Grain I (I = A, B) has a thickness h I in x 2 -direction and is infinitely long in x 1 -direction. For each grain, without loss of generality, two slip systems with each of them defined by a slip direction unit vector s Following [38], the misorientation angle ∆θ = θ (α) B and, following [47], the rotation angle of the GB relative to the symmetry plane θ GB = π − θ (1) A − θ (2) B are introduced to measure the geometric mismatch between the grains. If θ GB = 0, the GB is symmetric. Crystallographic parameters including the grain misorientation, GB orientation, and grain size are varied to illustrate their influence on the plastic behavior of the thin film. Therefore, the considered bicrystal thin film is not specified to a certain material. The proposed model applies to general fine-grained polycrystalline metallic materials where the plasticity is mediated by dislocation activities. The thin film suffers an externally prescribed shear rateγ ext with its lower surface being constrained such that the macroscopic boundary conditions are put down as where ∆u 1 and ∆u 2 denote the displacement increments after time increment ∆t, and h = h A + h B . The continuity condition of the displacement at the GB reads In the present plane strain problem, only edge GNDs are involved. Consider that slip transmission occurs at the GB. The corresponding relations between incoming and outgoing slip systems can be determined by taking advantage of slip transmission criteria [60,61], which is not elaborately pursued here. It is assumed that the incoming slip systems 1 and 2 in grain A correspond to the outgoing slip systems 1 and 2 in grain B respectively. Therefore, the microscopic boundary conditions at the GB arė AB ξ (2) Bi · N GB Bi − η AB ξ (2) Bi · N GB Bi − η where Equations (48) and (50) are the GB kinematic relations, and Equations (49) and (51) are the governing equations for slip transmission. In addition, in order to mimic the infinite extension of the thin film in x 1 -direction, a representative part of the thin film with length L is considered, and periodic boundary conditions are imposed for displacements and GND densities, i.e., The elastic parameters and the magnitude of the Burgers vector representative of aluminum are taken, i.e., µ = 26.3 GPa, ν = 0. 33  To simplify the discussion, GB parameters for the two slip transmission process are assumed to be the same, i.e., λ In addition, the average distance between slip planes of the incoming slip system is assumed to equal that of the corresponding outgoing slip system such that D

Influence of Energetic and Dissipative GB Parameters
In this section, the influences of the energetic and dissipative GB parameters are addressed.
A symmetric tilt GB with θ (1) The values of grain thickness are taken as h A = h B = 1 µm. For different values of the energetic GB hardening modulus λ GB , stress-strain curves, the evolution of the rate of slip transmission, the distribution of plastic slip at γ ext = 0.01, and the distribution of GND density at γ ext = 0.01 are plotted in Figure 3, where R GB 0 = 5 N/m and H GB = 0. Since the results of slip systems 1 and 2 are similar, only those of the slip system 1 are shown. In addition, the results for impenetrable GBs at which plastic slips vanish are also given for comparison. The corresponding results for different values of the dissipative GB hardening modulus and the initial GB slip resistance are displayed in Figure 4, where λ GB = 0. In Figures 3a and 4a, on each stress-strain curve (curves with symbols), two yielding points are observed. The first one is the traditional bulk yielding point. The second one denotes the onset of slip transmission after which the rate of slip transmission immediately increases from zero to a steady-state value (see Figures 3b and 4b). From Figure 4a,b, the critical load for the onset of slip transmission is governed by the initial GB resistance R GB 0 in a way that a larger R GB 0 gives rise to a larger critical load. As slip transmission weakens the dislocation pile-up at the GB, the strain hardening rate decreases after the onset of slip transmission. With the increase of the energetic and dissipative GB hardening modulus λ GB and H GB , it is more and more difficult for dislocations to penetrate the GB. Consequently, the steady-state value of the rate of slip transmission decreases, and the strain hardening rate at the stage with slip transmission increases. As shown in Figures 3c,d and 4c,d, for a smaller GB hardening modulus and/or a smaller initial GB resistance, the curve of distribution of plastic slip is smoother at the GB, and the GND density at the GB is smaller, indicating that slip transmission occurs more easily. In addition, if the initial GB slip resistance and/or the GB hardening modulus approach infinity, the present GB model reduces to the impenetrable GB model. The tendency of the stress-strain curves and the dependence of the strain hardening rate on the GB hardening modulus are consistent with those predicted by existing models [37][38][39]47,49].

Influence of Grain Size
To investigate the influence of grain size, for different values of grain thickness, the stress-strain curve, the evolution of the total average GND densityρ ge = h 0 ρ ge(1) + ρ ge(2) dx 2 /h, the evolution of the rate of slip transmission, and the distribution of plastic slip at γ ext = 0.01 are plotted in Figure 5a B − 60 • = 50 • with ∆θ = 20 • and θ GB = 0 • . By comparing the results of thin films with different values of grain thickness, it is seen that the plastic behavior is significantly size-dependent. From Figure 5a,b, for a specified load γ ext , a smaller grain thickness results in a larger flow stress and a larger total average GND density. The size-dependence of the total average GND density in Figure 5b is qualitatively consistent with the corresponding result by the DDD simulation in [20]. As displayed in Figure 5a,c, for the thin film with a smaller grain thickness, the critical load for the onset of slip transmission is smaller, and the steady-state value of the rate of slip transmission is larger, implying that the occurrence of slip transmission is easier. From Figure 5d, with the decrease of grain thickness, the plastic slip decreases in the middle of grains but increases near the GB. This is due to the fact that in the middle, the plastic slip is governed by dislocation interaction which is stronger if the grain thickness is smaller, while, near the GB, the plastic slip is dominated by slip transmission which is easier to occur in the film with a smaller grain thickness.

Influence of Grain Misorientation and GB Orientation
Firstly, the influence of grain misorientation is examined. To this end, four different values of ∆θ, i.e., 2 • , 5 • , 10 • , and 20 • are considered. The corresponding orientation angles of the slip systems are θ (1) A = 120 • + ∆θ/2, and θ The GB parameters are taken as λ GB = 20,000 N/m, R GB 0 = 15 N/m, and H GB = 10,000 N/m. In Figure 6, for different ∆θ, the evolution of the rate of slip transmission is plotted. It is seen that the steady-state value of the rate of slip transmission increases with the decrease of grain misorientation angle ∆θ. This indicates that slip transmission is more likely to occur if the grain misorientation angle is smaller, which is in accordance with experimental observations [62].  Figure 7 plots the grain size-dependence of the stress-strain curves for those four grain misorientation angles. Comparing the different cases, one can conclude that the influence of grain misorientation on the size-dependence of strain hardening rate after the onset of slip transmission is vital. In the case with a smaller grain misorientation angle (∆θ = 2 • or ∆θ = 5 • ), as slip transmission easily occurs, the dislocation pile-up at the GB is significantly relieved. Consequently, the strain hardening rate obviously decreases compared to the previous stage without slip transmission. Conversely, in the case with a larger grain misorientation angle (∆θ = 10 • or ∆θ = 20 • ), as slip transmission is difficult to proceed due to the strong GB resistance, the change of the strain hardening rate is trivial. In addition, for a given grain size, the larger the grain misorientation, the larger the flow stress at a specified load. These results are qualitatively consistent with the corresponding DDD simulation results in [18].
Then, the influence of the GB orientation is illustrated. To this end, a symmetric GB (θ GB = 0 • ) and an asymmetric GB (θ GB = 10 • ) are considered. The corresponding orientation angles of the slip systems are θ (1) Figure 8a,b, the evolution of the rate of slip transmission and the distribution of plastic slip at γ ext = 0.01 for the two cases are plotted. The GB orientation significantly affects slip transmission. For the case with a symmetric GB (θ GB = 0 • ), the critical loads for the onset of the two slip transmission processes are almost the same, and the difference in the steady-state value of the rate of slip transmission is small. For the case with an asymmetric GB (θ GB = 10 • ), the distribution of plastic slip and the evolution of the rate of slip transmission for the slip transmission process 1 obviously differ from those for the slip transmission process 2. From Figure 8b, the slip transmission process 2 occurs first. After the onset of the slip transmission process 1, the rate of slip transmission for the slip transmission process 2 suffers a sudden increase.

Influence of Surface Constraint
In the above discussion, in order to illustrate the GB influence separately, the microhard (impenetrable) surface which is completely constrained is considered. However, if the surface is permitted to absorb dislocations, slip transmission processes at the GB are also affected. Thus, the influence of surface constraint is investigated here. A microhard surface and a microfree surface freely penetrable by dislocations representing two extreme cases are considered. It is pointed out that the actual surface being in-between should be modeled by the surface model in Section 2.
The orientation angles of the slip systems are θ The GB parameters are λ GB = 20,000 N/m, R GB 0 = 15 N/m, and H GB = 10,000 N/m. Stress-strain curves, the evolution of the total average GND density, the evolution of the rate of slip transmission, and the distribution of plastic slip at γ ext = 0.01 for the two cases are shown in Figure 9. From Figure 9a,c, for the case with a microhard surface, the critical load for the onset of slip transmission is smaller and the steady-state value of the rate of slip transmission is larger, indicating the slip transmission is easier to occur. It is attributed to the fact that in the case with a microhard surface, since dislocations are not absorbed by surfaces, dislocation pile-up at the GB is stronger, giving rise to a larger driving force for slip transmission. Due to the surface constraint, in the case with a microhard surface, the decrease of the strain hardening rate and the increase rate of average GND densities after slip transmission is trivial compared to the case with a microfree surface. From Figure 9d, as the governing mechanism for the development of plastic slip near the surface and that near the GB are dislocation absorption by surfaces and slip transmission respectively, in the case with a microhard surface, plastic slip near the surface is smaller but that near the GB is larger. It reflects the competition between dislocation absorption by surfaces and slip transmission.

Summary
In this work, a GB model with consideration of slip transmission is proposed within the framework of finite deformation gradient polycrystal plasticity. In the GB model, for each slip transmission process, there are a GB kinematic relation between the plastic slip of the incoming slip system and that of the outgoing slip system at the GB and a governing equation for slip transmission, which constitute the GB microscopic boundary conditions. Both the energetic and the dissipative GB effects are considered by introducing an energetic and a dissipative GB microforce for each slip transmission process. By properly constructing energetic and dissipative GB constitutive relations, the important factors including grain misorientation, the GB orientation, GB defects and stress state at the GB which may affect slip transmission are captured in the present GB model. The main advantage of the present model over the existing models is the consideration of underlying physical mechanisms of slip transmission.
The mathematical framework of the present rate-dependent model is a fully coupled, strongly nonlinear initial boundary value problem with non-standard boundary conditions, the numerical implementation of which is done via a dual-mixed nonlinear finite element method. The treatment of microscopic GB conditions is emphasized. A benchmark problem of a bicrystal thin film under plane constrained shear is studied. It is found that the GB parameter, grain size, grain misorientation, GB orientation, and surface constraint can significantly affect the slip transmission and the plastic behavior of thin films. Particularly, in thin films with a smaller grain size and/or a smaller grain misorientation, slip transmission occurs more easily. These results predicted by the present model qualitatively agree with some experimental observations and results of DDD simulations in the literature.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
To realize the numerical implementation of the model, a dual-mixed finite element formulation put forward by [63] is adopted. The balance of momentum (24) and the evolution Equations (8) and (9) of the GND densities are treated as governing equations such that the displacement u and the GND densities ρ ge(α) i and ρ gs(α) i act as global variables. The plastic slip γ (α) is evaluated locally at the integration points by the evolution equation of the plastic slip in the bulk and that at the GB/surface. The treatment for the balance of momentum is standard. Thus, attention is restricted to the discretization and implementation of the evolution equations of the GND densities. As the evolution equations of the edge GND density and the screw GND density are similar, only the former is explained for the sake of brevity. To obtain its weak forms, the governing equation of the edge GND density is multiplied by the weight function δρ ge(α) i and then integrated over the volume, followed by the application of divergence theorem, which yields where the GB kinematic relation (14) is considered. For the spatial discretization, the volume of the body is divided into finite elements, in each of which, following the standard Galerkin approach, the unknown fields of the displacement and GND densities and the weight functions are approximated by their nodal values multiplied by the shape functions. In addition, the adjacent elements on the two sides of the GB are decoupled by placing double nodes at a GB nodal point. To satisfy the continuity condition of the displacement at the GB, the displacement of those node pairs are forced to be equal. For the benchmark problem in the present work, 4-nodes bilinear quadrilateral elements with 2 × 2 full integration are chosen for the spatial discretization of both the displacement and the GND densities. To evaluate the GB/surface integration, two additional integration points are placed on each surface/GB edge of elements at the GB/surface. For the time integration within a time increment [t n , t n+1 ], the implicit backward Euler scheme is adopted such thatγ , ∆t = t n+1 − t n and the subscript n identifies the quantities at the nth time step.
For the incorporation of the GB/surface contributions, the evaluation of the increment of the plastic slip at the GB/surface integration points is required. As discussed in [52], at the surface, the increment of the plastic slip ∆γ (α) is determined by the governing Equation (38) for dislocation absorption by surfaces The plastic slips required to evaluate the terms on the right hand side are extrapolated from the bulk. Accordingly, ∆γ where the plastic slips required to calculate the quantities on the right hand side are extrapolated from the interior of grain A or grain B. Since quantities from both grains are involved in Equation (A4), the adjacent elements on the two sides of the GB are coupled, by which intergranular interactions between dislocations are considered. By combing the spatial and temporal discretization, the weak forms of the governing equations reduce to a highly nonlinear, strongly coupled system of equations at each time step, which is solved numerically by a Newton-Raphson solution strategy.