Thermodynamic theory of dislocation/grain boundary interaction

Abstract

The thermodynamic theory of dislocation/grain boundary interaction, including dislocation pileup against, absorption by, and transfer through the grain boundary, is developed for non-uniform plastic deformations in polycrystals. The case study is carried out on the boundary conditions affecting work hardening of a bicrystal subjected to plane constrained shear for three types of grain boundaries: (i) impermeable hard grain boundary, (ii) grain boundary that allows dislocation transfer without absorption, and (iii) grain boundary that absorbs dislocations and allows them to pass later.

1 Introduction

The interaction between dislocations and grain boundaries plays an essential role in the mechanical response of polycrystalline materials to loading. Several outcomes of this interaction are possible. At low dislocation density, the grain boundary blocks the movement of dislocations and forces them to pileup against it, with the by-product of linearly increased kinematic work hardening. As the dislocation density exceeds some critical value, a certain amount of dislocations is absorbed by the grain boundary, causing an increase in the surface dislocation density and the misorientation between adjacent grains. Dislocations can also be emitted at a later stage of straining in the neighboring grain or transferred through the grain boundary. In the presence of grain boundaries, non-uniform plastic deformation occurs, leading to non-redundant (geometrically necessary) dislocations, so the need for continuum models that include the gradient of plastic deformation becomes evident. Berdichevsky and Sedov [1] were the first to introduce the curl of plastic deformation (Nye–Bilby–Kröner’s tensor) into the continuum dislocation theory (CDT) and associate it with the density of non-redundant dislocations. They also proposed the thermodynamics of dislocations, using this density of non-redundant dislocations as a state variable in the free energy and giving the variational formulation of CDT (see the further developments in [2, 3]). The main drawback of this thermodynamic approach is the absence of configurational entropy and its dual quantity, the effective temperature of dislocations, introduced a few decades later by Langer et al. [4]. Together with the density of redundant dislocations, this configurational entropy enables ones to build up the meaningful thermodynamics of dislocations satisfying two universal laws for plastic flows discovered in [5, 6]. Recently, one of us extended CDT by introducing the density of redundant dislocations as well as the configurational entropy as additional state variables [7]. He has shown that the new theory, called thermodynamic dislocation theory (TDT), can predict both kinematic and isotropic work hardening. This TDT has been applied to model grain boundaries as hard obstacles that block dislocations, leading to dislocation pileup, kinematic work hardening, and Bauschinger effect [8, 9]. Piao and Le [10] further applied the theory to the problem of interaction between screw dislocations and permeable grain boundary based on the experiments performed by Kondo et al. [11]. There, an interfacial dissipation potential is proposed and additional boundary conditions for the grain boundary are derived to model screw dislocation/grain boundary interaction.

Alternative phenomenological models that include the interaction between dislocations and grain boundaries have been proposed in the context of strain gradient plasticity by Fleck et al. [12]; Aifantis [13]; Gao et al. [14]; Gurtin [15], and a large number of investigations have been carried out within this approach. Fredriksson and Gudmundson [16] proposed an energetic model according to which the plastic work is stored as interfacial energy that depends on the jump of plastic slip at the grain boundary, while Fleck and Willis [17] argued that material hardening due to the grain boundary is mainly carried by energy dissipation. Note that both models allow for a discontinuity in the plastic slip at the grain boundary. Gurtin [18] introduced a Burgers’ tensor at the grain boundary, analogous to the Nye–Bilby–Kröner’s tensor for bulk crystals, to determine both the interfacial energy and dissipation. One merit of this tensorial measure is that it allows the theory to satisfy the slip transfer criteria. Van Beers et al. [19] proposed a similar theory that takes into account the criteria of slip transfer, but differs from the formulation of Gurtin [18]. Inter-grain interaction modules from two theories are compared by Gottschalk et al. [20]; Bayerschen et al. [21]. Erdle and Böhlke [22] developed the continuum model of dislocation/grain boundary interaction based on the micromorphic strain gradient plasticity. We also mention the references [23,24,25,26,27,28,29] which use the same strain gradient plasticity. Similar to the CDT developed in [1,2,3], this standard strain gradient plasticity ignores completely configurational entropy, so none of the work done within this approach is consistent with thermodynamics of dislocations proposed in [4,5,6].

The aim of this work is to model the interaction of edge dislocations with the grain boundary in polycrystals within the framework of TDT. In contrast to the interaction between screw dislocations and grain boundaries studied in our previous paper [10], the processes of dislocation pileup against, absorption by, and transfer through the grain boundary may occur in this situation. Using the variational formulation together with the proposed surface energy and dissipation potential, we will derive the interface boundary conditions for these processes. The case study is carried out on the boundary conditions affecting work hardening of a bicrystal subjected to plane constrained shear for three types of grain boundaries: (i) impermeable hard grain boundary, (ii) grain boundary that allows dislocation transfer without absorption, and (iii) grain boundary that absorbs dislocations and allows them to pass later. Note that, according to various experimental observations reported, e.g., in [30, 31], the grain boundary could be migrated during this dislocation absorption, but this possible migration is ignored in this continuum model. We show especially the effect of dislocation absorption on work hardening.

The paper is organized as follows. After this brief introduction, we lay down the kinematics and thermodynamic framework of TDT in Sect. 2. In Sect. 3, we analyze the effects of dissipation potential and interfacial energy on the work hardening of materials and derive the conditions at the grain boundary for various processes during plastic deformation. In Sect. 4, the application of the developed theory to the problem of a bicrystal subjected to plane constrained shear is illustrated and its numerical treatment is presented. Detailed numerical simulations are then performed and the results are discussed. We conclude the paper with a brief summary in Sect. 5.

2 Formulation for bulk crystals

2.1 Kinematics

In the small strain theory, we ignore the differences between the Lagrangian and Eulerian coordinates and assume that the total displacement gradient (distortion) ${\mathbf {u}}\nabla$ can be decomposed additively into elastic and plastic counterparts

$$\begin{aligned} {\mathbf {u}} \nabla =\varvec{\beta }^\text {e}+\varvec{\beta }, \end{aligned}$$

(1)

where ${\mathbf {u}}({\mathbf {x}})$ is a displacement field, ${\mathbf {x}}=(x_1,x_2,x_3)$ a position vector of a material point of the body, and $\nabla$ the nabla operator. The incompatible plastic distortion of the crystal having n active slip systems takes the form

$$\begin{aligned} \varvec{\beta }({\mathbf {x}})= \sum _{\alpha =1}^{n} \beta ^{\alpha }({\mathbf {x}}) {\mathbf {s}}^{\alpha } \otimes {\mathbf {m}}^{\alpha }, \end{aligned}$$

(2)

with ${\mathbf {s}}^{\alpha }$ and ${\mathbf {m}}^{\alpha }$ being the pair of unit vectors denoting the slip direction and the normal to the slip planes of the corresponding $\alpha$-th slip system. The plastic slip of the $\alpha$-th slip system, $\beta ^{\alpha }({\mathbf {x}})$, is assumed to be continuously differentiable function except at the interfaces or grain boundaries. Since ${\mathbf {s}}^{\alpha }$ and ${\mathbf {m}}^{\alpha }$ are mutually orthogonal, the trace of the plastic distortion always vanishes, tr$\varvec{\beta }=0$, which indicates volume-preserving plastic deformation due to the conservative motion of dislocations. The total strain $\varvec{\varepsilon }$ and the plastic strain $\varvec{\varepsilon }^\text {p}$ are the symmetric part of the displacement gradient and the plastic distortion, respectively, so that

$$\begin{aligned} \varvec{\varepsilon }=\frac{1}{2}( \nabla {\mathbf {u}}+{\mathbf {u}}\nabla ), \quad \varvec{\varepsilon }^\text {p}=\frac{1}{2}(\varvec{\beta }+\varvec{\beta }^T). \end{aligned}$$

(3)

Accordingly, the elastic strain equals

$$\begin{aligned} \varvec{\varepsilon }^\text {e}=\varvec{\varepsilon }-\varvec{\varepsilon }^\text {p}. \end{aligned}$$

(4)

The net Burger’s vector ${\mathbf {B}}$ of non-redundant dislocations whose lines cross the area ${\mathcal {S}}$ bounded by the closed curve $\partial {\mathcal {S}}$ is given by

$$\begin{aligned} {\mathbf {B}}= \oint _{\partial {\mathcal {S}}} \varvec{\beta } \cdot \mathrm {d} {\mathbf {x}}= \int _{\mathcal {S}} \varvec{\alpha }\cdot {\mathbf {n}} \mathrm {d}a, \end{aligned}$$

(5)

where $\varvec{\alpha }$ is the Nye–Bilby–Kröner’s dislocation tensor [32,33,34],

$$\begin{aligned} \varvec{\alpha }= -\varvec{\beta }\times \nabla , \quad (\alpha _{ij}=\epsilon _{jkl}\beta _{il,k}). \end{aligned}$$

(6)

In particular, if there is only one active slip system such that ${\mathbf {s}}$ and ${\mathbf {m}}$ lie in the $(x_1,x_2)$-plane, while the dislocation lines are parallel to the $x_3$-axis, then $\varvec{\alpha }$ under the plane strain has only two nonzero components $\alpha _{i3}=s_i(\partial _s \beta )$, where $\partial _s$ denotes the derivative in the ${\mathbf {s}}$-direction. Consequently, the scalar density of the non-redundant dislocations is given by

$$\begin{aligned} \rho ^\text {g}=\frac{|\mathrm {d}{\mathbf {B}}|}{b\mathrm {d}a}=\frac{1}{b}|\partial _s \beta |, \end{aligned}$$

(7)

with b being the magnitude of Burger’s vector.

2.2 Thermodynamic framework

Suppose the material body in form of a slab of constant depth L occupies the 3-D domain ${\mathcal {V}}$, and its boundary $\partial {\mathcal {V}}$ consists of two non-intersecting surfaces, $\partial _k$ and $\partial _s$, on which the displacement field $\mathbf {u(x)}$ and the traction-free condition are specified, respectively. Assuming that no body force acts on the crystals, the energy functional reads

$$\begin{aligned} I=\int _{{\mathcal {V}}} \psi (\varvec{\varepsilon }^\text {e}, \rho ^\text {g}, \rho ^\text {r}, \chi ) \mathrm {d}^3x, \end{aligned}$$

(8)

where $\mathrm {d}^3x=\mathrm {d}x_1\mathrm {d}x_2\mathrm {d}x_3$ denotes the volume element and $\psi$ is the density of the Helmholtz free energy. The dislocation network is assumed to be two-dimensional, such that the dislocation lines are parallel to the $x_3$-axis and the Burgers’ vector of edge dislocations lies in the $(x_1,x_2)$-plane. We let $\rho ^\text {r}$ denote the density of redundant dislocations whose resultant Burgers’ vector vanishes, and $\rho ^\text {g}$ the density of non-redundant dislocations (see subsection 2.1 ). The sum of the two types of dislocation densities is the total dislocation density $\rho =\rho ^\text {r}+\rho ^\text {g}$. In the spirit of the dislocation mediated plasticity [4], which decomposes the thermodynamic system of dislocated crystal into two subsystems and traces it back to statistical aspects of the defects, $\chi$ is the effective temperature characterizing the slow rearrangement of configuration of atoms during the motion of dislocations, while the ordinary temperature T characterizes the fast vibration of atoms of the kinetic vibrational subsystem. In some cases, such as thermal softening [35,36,37] or adiabatic shear banding [38], the ordinary temperature plays an important role. However, in the present study we set T to be constant and drop it from the list of arguments of the free energy density $\psi$. The state variables $\varvec{\varepsilon }^\text {e}$ and $\rho ^\text {g}$ are dependent variables expressed by ${\mathbf {u}}$ and $\beta$, while $\rho ^\text {r}$ and $\chi$ are independent state variables. We decompose $\psi$ into the part due to the elastic strain, $\psi _\text {e}$, and the remaining parts

$$\begin{aligned} \psi (\varvec{\varepsilon }^\text {e}, \rho ^\text {g}, \rho ^\text {r}, \chi )=\psi _\text {e}(\varvec{\varepsilon }^\text {e})+\psi _\text {r}(\rho ^\text {r})+\psi _\text {m}(\rho ^\text {g})+\psi _\text {c}(\chi , \rho ), \end{aligned}$$

(9)

where $\psi _\text {r}$ is the self-energy density of redundant dislocations, $\psi _\text {m}$ is the energy density of the non-redundant dislocations, while $\psi _\text {c}$ is the configurational heat associated with the configurational entropy [7, 35]. The additive decomposition of the energy of dislocations into two parts $\psi _\text {r}$ and $\psi _\text {m}$ is possible because the interaction between dislocation dipole and dislocation has a much shorter range compared to the interaction between dislocations of the same sign. These contributions are given below

$$\begin{aligned} \begin{aligned}&\psi _\text {e}(\varvec{\varepsilon }^\text {e})=\frac{1}{2} \lambda (\varepsilon ^\text {e}_{kk})^2+\mu \varepsilon ^\text {e}_{ij}\varepsilon ^\text {e}_{ij}, \quad \psi _\text {r}(\rho ^\text {r})=\gamma _\text {D} \rho ^\text {r},\\&\psi _\text {m}(\rho ^\text {g})=\gamma _\text {D} \rho ^\text {g}_{ss} \ln \left( \frac{1}{1-\frac{\rho ^\text {g}}{\rho ^\text {g}_{ss}}}\right) ,\\&\psi _\text {c}(\chi , \rho )=-{\bar{\chi }}(-\rho \ln (a^2 \rho )+\rho ), \end{aligned} \end{aligned}$$

(10)

where $\lambda$ and $\mu$ are Lam$\acute{\text {e}}$ constants, and $\gamma _\text {D}$ the dislocation energy per unit length. Since $-\rho \ln (a^2 \rho )+\rho$ is the entropy per unit area of the two-dimensional dislocation network (see the detailed derivation in [39]), the “two-dimensional” temperature ${\bar{\chi }}=\chi /L$ is introduced to get the configurational heat having the unit of energy density [7]. Note that a is a length scale of the order of atomic spacing. The defect energy $\psi _\text {m}$ describing the interactions of non-redundant dislocations captures the kinematic hardening. The logarithmic term, originated from [40], ensures that the defect energy increases linearly for small dislocation density and tends to infinity as $\rho ^\text {g}$ approaches the saturated density of non-redundant dislocations $\rho ^\text {g}_{ss}$. We assume that $\rho ^\text {g}_{ss}=k_0 \rho _{ss}$ is a fraction of $\rho _{ss}$ by coefficient $k_0$, where $\rho _{ss}=(1/a^2)e^{-\gamma _\text {D}/{\bar{\chi }}}$ is the steady-state dislocation density determined by minimizing the free energy of the configurational subsystem [4]. Note that $k_0$ is characterized by the interaction between dislocations of the same sign, which is similar to the interaction between charges [41]. We will see in the subsequent section the role of $\rho ^\text {g}_{ss}$ in the contribution of interfacial energy to the work hardening.

The applied loading can lead to nucleation, propagation and movement of dislocations, which result in plastic deformations in the crystal. In these processes, dislocations always experience resistance, which causes energy dissipation. The dissipation potential in the bulk has the form

$$\begin{aligned} D_\text {b}({\dot{\beta }}, {\dot{\rho }}, {\dot{\chi }})=\tau _\text {i} {\dot{\beta }} +\frac{1}{2}d_{\rho }{\dot{\rho }}^2+\frac{1}{2}d_{\chi }{\dot{\chi }}^2. \end{aligned}$$

(11)

The first term is the plastic power, where $\tau _\text {i}$ is the internal stress, while the last two terms represent the dissipation due to the multiplication of dislocations and the increase of the configurational temperature [7]. Langer et al. [4] developed the kinetics of dislocation depinning that dominantly controls the plastic slip rate. Averaging this kinetic equation over the representative volume element, we have established in [42] that

$$\begin{aligned} \dot{{\bar{\beta }}}=\frac{q(T,\tau _\text {i},\rho ^\text {r})}{t_0}=\frac{b}{t_0}\sqrt{\rho }\left[ f_\text {P}(T,\tau _\text {i},\rho ^\text {r})-f_\text {P}(T,-\tau _\text {i},\rho ^\text {r}) \right] \end{aligned}$$

(12)

where $t_0$ is a microscopic time scale of the order of the inverse Debye frequency, and

$$\begin{aligned} f_\text {P}(T, \tau _\text {i}, \rho ^\text {r}) = \exp \left[ -\frac{T_\text {P}}{T} \exp \left( -\frac{\tau _\text {i}}{\tau _\text {T}(\rho ^\text {r})}\right) \right] . \end{aligned}$$

(13)

Note that $T_\text {P}$ is the activation temperature of the potential well of a pinning site and $\tau _\text {T}(\rho ^\text {r})=\mu _\text {T} b \sqrt{\rho ^\text {r}}$ is the Taylor stress with $\mu _\text {T}$ being proportional to the shear modulus $\mu$. It was shown that this kinetic equation for $\dot{{\bar{\beta }}}$ can correctly describe the rate-reducing behavior of isotropic work hardening both when the crystal is loaded on one direction and in the load reversal process [10]. Provided there is only one slip system inclined at an angle $\phi$ to the $x_1$-axis, then, using similar arguments as in [8], we lay down one governing equation for ${\dot{\tau }}_\text {i}$

$$\begin{aligned} {\dot{\tau }}_\text {i}=\mu \frac{q_0}{t_0}\left( \cos 2\phi -\frac{q(T,\tau _\text {i},\rho ^\text {r})}{q_0}\right) \end{aligned}$$

(14)

where $q_0/t_0$ denotes the total shear strain rate. This equation plays the role of the kinematic constraint for the internal stress $\tau _\text {i}$. The remaining equations in the bulk can be derived from the following variational principle: The true displacement field $\check{{\mathbf {u}}}\mathbf {(x)}$, the true plastic slips ${\check{\beta }}({\mathbf {x}})$, the true density of redundant dislocations ${\check{\rho }}^\text {r}({\mathbf {x}})$, and the true configurational temperature ${\check{\chi }}({\mathbf {x}})$ obey the variational equation

$$\begin{aligned} \delta I +\int _{{\mathcal {V}}} \left( \frac{\partial D_\text {b}}{\partial {\dot{\beta }}} \delta \beta +\frac{\partial D_\text {b}}{\partial {\dot{\rho }}} \delta \rho +\frac{\partial D_\text {b}}{\partial {\dot{\chi }}} \delta \chi \right) \mathrm {d}^3x=0 \end{aligned}$$

(15)

for all variations of admissible fields ${\mathbf {u}}({\mathbf {x}})$, $\beta ({\mathbf {x}})$, $\rho ^\text {r}({\mathbf {x}})$ and $\chi ({\mathbf {x}})$. Focusing on the contribution of the grain boundary to work hardening, we assume that the whole external boundary $\partial {\mathcal {V}}$ admit dislocation to reach it. In this case, the plastic slip can be varied arbitrarily at $\partial {\mathcal {V}}$. The variational equation (15) allows one to derive the remaining governing equations. From the variation of I with respect to the displacement field $\mathbf {u(x)}$, we obtain the quasi-static equilibrium equation and the boundary condition [43]

$$\begin{aligned} \begin{aligned}&\sigma _{ij,j}=0 \quad \hbox { in}\ {\mathcal {V}},\\&\sigma _{ij}n_j=0 \quad \hbox { on}\ \partial _s, \end{aligned} \end{aligned}$$

(16)

with $\varvec{\sigma }=\partial \psi _\text {e}/\partial \varvec{\varepsilon }^\text {e}$ being the Cauchy stress tensor and ${\mathbf {n}}$ the outward unit normal vector to the external boundary. In Eq. (15), the vanishing variation with respect to the plastic slip $\beta ({\mathbf {x}})$ yields the balance of microforces acting on dislocations [7] and the natural boundary condition at the entire external boundary as

$$\begin{aligned} \begin{aligned}&\tau -\tau _\text {b}-\tau _\text {i}=0 \quad \text { in}\ {\mathcal {V}},\\&\frac{\partial \psi _\text {m}}{\partial \rho ^\text {g}}=\gamma _\text {D} \quad \text {on}\ \partial {\mathcal {V}}, \end{aligned} \end{aligned}$$

(17)

where $\tau =s_i \sigma _{ij}m_j$ is the resolved shear stress acting on the slip system, while the back stress is

$$\begin{aligned} \tau _\text {b}=-\frac{1}{b}\left( \frac{\partial \psi _\text {m}}{\partial \rho ^\text {g}} \text {sign}\beta _{,s}\right) _{,s}=-\frac{1}{b^2}\frac{\partial ^2\psi _\text {m}}{\partial (\rho ^\text {g})^2}\beta _{,ss}. \end{aligned}$$

(18)

Note that $(.)_{,s}=s_i \partial _i(.)$ is the shorthand notation for $\partial _{s}(.)$. Following the suggestions by Le [7] for functions of $d_{\rho }$ and $d_{\chi }$, the remaining equations obtained from Eq. (15) for ${\dot{\chi }}$ and ${\dot{\rho }}$ can be cast into the form

$$\begin{aligned} \begin{aligned}&{\dot{\chi }}={\mathcal {K}}_{\chi } e_\text {D} \tau _\text {i} \frac{q(T,\tau _\text {i},\rho ^\text {r})}{t_0}\left( 1-\frac{\chi }{\chi _0}\right) , \\&{\dot{\rho }}={\mathcal {K}}_{\rho } \frac{\tau _\text {i}}{a^2 \, \nu (T, \rho ^\text {r}, q_0)^2} \frac{q(T,\tau _\text {i},\rho ^\text {r})}{t_0} \left( 1-\frac{\rho }{\rho _{ss}(\chi )}\right) . \end{aligned} \end{aligned}$$

(19)

Here, $\chi _0$ is a constant denoting the steady-state configurational temperature, while the steady-state dislocation density at a given effective temperature is $\rho _{ss}(\chi )=(1/a^2)e^{-\gamma _\text {D}/{\bar{\chi }}}$. ${\mathcal {K}}_{\chi }$ is a factor inversely proportional to the effective specific heat and ${\mathcal {K}}_{\rho }$ is a energy conversion coefficient [44]. The function $\nu (T, \rho ^\text {r}, q_0)$ is defined as

$$\begin{aligned} \nu (T,\rho ^\text {r}, q_0)= \ln \left( \frac{T_\text {P}}{T}\right) -\ln \left[ \ln \left( \frac{b\sqrt{\rho ^\text {r}}}{q_0}\right) \right] . \end{aligned}$$

(20)

3 Interface boundary conditions

As a prototype, we study a simple problem of plane deformation of a bicrystal whose single grain boundary, which is a plane perpendicular to the y-axis, lies at $y=y_\text {g}$ (see Fig. 1). The extension to the case of polycrystals with multiple grain boundaries can be done without complications. Since the tensorial index notation is no longer required, the coordinates $(x_1,x_2,x_3)$ are changed to (x, y, z) for simplicity. For the sake of definiteness, we assume that only edge dislocations move to and interact with the grain boundary. If the plastic slip $\beta$ depends only on y due to the sample geometry and specific boundary conditions, then the non-redundant dislocation density is proportional to $\beta _{,y}$, $\rho ^\text {g}=\varkappa | \beta _{,y}|/ b$, where $\varkappa$ depends on the slip system, as will be discussed in Sect. 4.

3.1 Interfacial energy and dissipation

Dislocations are long-lived and well-defined defects and have various interactions with the grain boundary, e.g., dislocation pileup, absorption, and transmission. Within the continuum approach, it should be possible to mimic these interactions by the interfacial energy and dissipation that reflect the underlying physics at the nano-/submicron scale. The low-angle grain boundary itself can be regarded as an array of dislocations [45] whose interfacial energy depends on the surface dislocation density [46, 47]. It also acts as an additional barrier to the movement of dislocations and can later absorb or release them. Therefore, in addition to bulk dissipation, interfacial dissipation must also be taken into account. In order to choose the appropriate argument for the interfacial energy, the state variables of the surface must first be recognized. Although the plastic distortion and the corresponding plastic slip are not the state variables in the bulk, the jump in the plastic slip representing the surface dislocation density is the state variable at the grain boundary (cf. with the Burgers tensor introduced in [18]). We propose the densities of interfacial energy and dissipation potential as follows

$$\begin{aligned} \psi _\text{ s }(\beta \bigr |_{y_\text{ g }\pm 0})=\psi _\text{ s }(|[\![\beta ]\!]_{\varGamma }|), \end{aligned}$$

(21)

and

$$\begin{aligned} D_\text{ s }({\dot{\beta }} \bigr |_{y_\text{ g }\pm 0})=\zeta _\text{ y } \frac{|\langle \!\langle {\dot{\beta }}\rangle \!\rangle _{\varGamma }|}{b}+\frac{1}{2} \zeta _\text{ a }\frac{([\![{\dot{\beta }}]\!]_{\varGamma })^2}{b^2}. \end{aligned}$$

(22)

Here, the vertical line followed by $y_\text {g}\pm 0$ indicates the limits of the preceding expression as y approaches the grain boundary position $y_\text {g}$ from above and below, while $\langle \!\langle . \rangle \!\rangle _{\varGamma }$ and $[\![. ]\!]_{\varGamma }$ denote the mean value and the jump of a variable at the interface $\varGamma$, respectively. Thus, according to these definitions

$$\begin{aligned} \langle \!\langle {\dot{\beta }}\rangle \!\rangle _{\varGamma }=\frac{1}{2}({\dot{\beta }}\bigr |_{y_\text{ g }+0}+{\dot{\beta }}\bigr |_{y_\text{ g }-0}), \quad [\![{\dot{\beta }}]\!]_{\varGamma }=({\dot{\beta }}\bigr |_{y_\text{ g }+0}-{\dot{\beta }}\bigr |_{y_\text{ g }-0}).\end{aligned}$$

(23)

Eq. (21) requires that only the surface dislocation density $|[\![\beta ]\!]_{\varGamma }|/b$, regarded as the state variable, can be the argument of $\psi _\text {s}$. The first term on the right-hand side of Eq. (22) is the dissipation by the rate of the mean plastic slip at the interface, where $\zeta _\text {y}$ plays a role of the yield surface tension. The second term in (22), quadratic in $[\![{\dot{\beta }}]\!]_{\varGamma }$, describes rate-dependent dissipation due to the acoustic emission of waves generated during the absorption of dislocations by the grain boundary.

The boundary conditions involving interfacial energy and dissipation can be derived from the following variational equation

$$\begin{aligned}&\delta \left( I+\int _{\varGamma }\psi _\text{ s }(|[\![\beta ]\!]_{\varGamma }|)\mathrm {d}a\right) +\int _{{\mathcal {V}}} \left( \frac{\partial D_\text{ b }}{\partial {\dot{\beta }}} \delta \beta +\frac{\partial D_\text{ b }}{\partial {\dot{\rho }}} \delta \rho +\frac{\partial D_\text{ b }}{\partial {\dot{\chi }}} \delta \chi \right) \mathrm {d}^3x \nonumber \\&\quad + \int _{\varGamma }\left( \frac{\partial D_\text{ s }}{\partial [\![{\dot{\beta }}]\!]_{\varGamma }}\delta [\![\beta ]\!]_{\varGamma }+\frac{\partial D_\text{ s }}{\partial \langle \!\langle {\dot{\beta }}\rangle \!\rangle _{\varGamma }} \delta \langle \!\langle \beta \rangle \!\rangle _{\varGamma }\right) \mathrm {d}a=0,\end{aligned}$$

(24)

which extends Eq. (15) to the case involving the grain boundary, where $\mathrm {d}a$ denotes the area element, while I and $D_\text {b}$ remain exactly the same as in (8) and (11), respectively. Edge dislocations may either penetrate the grain boundary or be absorbed by it through dissociation of the dislocations. The latter process requires the continuity of displacement field but causes a jump in plastic slip $[\![\beta ]\!]_{\varGamma }$ at the interface. This is revealed by the volume integral containing the variation of the non-redundant dislocation density $\delta \rho ^\text {g}$ during the integration by parts, permitting the discontinuity to enter into the interface as

$$\begin{aligned}&\int _{{\mathcal {V}}} \left( \frac{\partial \psi _\text{ m }}{\partial \rho ^\text{ g }}+\frac{\partial \psi _{\chi }}{\partial \rho ^\text{ g }}+\frac{\partial D_\text{ b }}{\partial {\dot{\rho }}^\text{ g }}\right) \delta \rho ^\text{ g } \mathrm {d}^3x =-\int _{{\mathcal {V}}} \frac{\varkappa }{b}\left( \frac{\partial \psi _\text{ m }}{\partial \rho ^\text{ g }}-\gamma _\text{ D }\right) _{,y}\text{ sign }\beta _{,y} \delta \beta \mathrm {d}^3x \nonumber \\&\quad +\int _{\partial _s} \frac{\varkappa }{b}\left( \frac{\partial \psi _\text{ m }}{\partial \rho ^\text{ g }}-\gamma _\text{ D }\right) \text{ sign }\beta _{,y} \delta \beta \mathrm {d}a -\int _{\varGamma } [\![\frac{\varkappa }{b}\left( \frac{\partial \psi _\text{ m }}{\partial \rho ^\text{ g }}-\gamma _\text{ D }\right) \text{ sign }\beta _{,y} \delta \beta ]\!]_{\varGamma } \mathrm {d}a. \end{aligned}$$

(25)

The governing equation with respect to $\rho ^\text {r}$ from Eq. (15), given as [7, 48]

$$\begin{aligned} \gamma _\text {D}+\chi \ln (a^2 \rho )/L+d_{\rho } {\dot{\rho }}=0, \end{aligned}$$

(26)

is utilized for $\gamma _\text {D}$ in the derivation of Eq. (25). Expanding the interface integral term on the right-hand side of Eq. (25) by means of the jump identity [49]

$$\begin{aligned} \begin{aligned}{}[\![p \, q]\!]_{\varGamma }=[\![p]\!]_{\varGamma }\langle \!\langle q\rangle \!\rangle _{\varGamma }+ [\![q]\!]_{\varGamma }\langle \!\langle p\rangle \!\rangle _{\varGamma } \end{aligned}\end{aligned}$$

(27)

and substituting it into Eq. (24), we transform the latter to

$$\begin{aligned}&\int _{\varGamma } \left\{ -\langle \!\langle \frac{\varkappa }{b}\left( \frac{\partial \psi _\text{ m }}{\partial \rho ^\text{ g }}-\gamma _\text{ D }\right) \text{ sign }\beta _{,y}\rangle \!\rangle _{\varGamma } \delta [\![\beta ]\!]_{\varGamma } -[\![\frac{\varkappa }{b}\left( \frac{\partial \psi _\text{ m }}{\partial \rho ^\text{ g }}-\gamma _\text{ D }\right) \text{ sign }\beta _{,y}]\!]_{\varGamma } \delta \langle \!\langle \beta \rangle \!\rangle _{\varGamma }\right. \nonumber \\&\quad \left. + \frac{\partial \psi _\text{ s }}{\partial [\![\beta ]\!]_{\varGamma }} \delta [\![\beta ]\!]_{\varGamma } +\zeta _\text{ y } \frac{1}{b}\text{ sign }\langle \!\langle {\dot{\beta }}\rangle \!\rangle _{\varGamma }\delta \langle \!\langle \beta \rangle \!\rangle _{\varGamma } +\zeta _\text{ a } \frac{1}{b^2}[\![{\dot{\beta }}]\!]_{\varGamma } \delta [\![\beta ]\!]_{\varGamma }\right\} \mathrm {d}a=0. \end{aligned}$$

(28)

To derive consequences from (28), we need to analyze the variations of plastic slip on both sides of the grain boundary, which may be subject to constraints depending on the accumulated dislocation densities. We therefore consider two different processes.

Pileup process: When dislocations move toward the grain boundary under the Peach–Koehler force, they first pile up against it. Since the dislocations cannot initially penetrate the grain boundary acting as an obstacle, the plastic slips on both sides of the grain boundary are subject to homogeneous Dirichlet boundary condition which fulfills Eq. (28) due to $\delta \beta \bigr |_{y_\text {g}-0}=\delta \beta \bigr |_{y_\text {g}+0}=0$. So we have during the pileup process

$$\begin{aligned} {\dot{\beta }}(y_\text {g}\pm 0, t)=0. \end{aligned}$$

(29)

The time derivative in Eq. (29) ensures the continuity of plastic slip at the interface with respect to time in the reversal loading [10].

Absorption and transmission processes: As more dislocations of the same sign pile up against the interface, the back stress near the pileup sites is increased. When the accumulated non-redundant dislocation density and with it the back stress in the vicinity of the grain boundary exceeds a threshold, the latter can no longer block dislocations. When the grain boundary starts to absorb dislocations, the jump in plastic slip as well as the misorientation between two adjacent grains will increase. The plastic slip at the grain boundary as well as its jump $[\![\beta ]\!]_{\varGamma }$ may change during the process of dislocation absorption, and consequently, $\delta \langle \!\langle \beta \rangle \!\rangle _{\varGamma }$ and $\delta [\![\beta ]\!]_{\varGamma }$ can be chosen independently and arbitrarily. Eq. (28) then implies

$$\begin{aligned} \varkappa \left( \frac{\partial \psi _\text {m}}{\partial \rho ^\text {g}}-\gamma _\text {D}\right) \Bigr |_{y_\text {g}+0}+\varkappa \left( \frac{\partial \psi _\text {m}}{\partial \rho ^\text {g}}-\gamma _\text {D}\right) \Bigr |_{y_\text {g}-0}=\zeta _\text {y}, \end{aligned}$$

(30)

and

$$\begin{aligned} {}[\![{\dot{\beta }}]\!]_{\varGamma }=\frac{b^2}{\zeta _\text{ a }}\left( [\![\frac{\varkappa }{2b}\frac{\partial \psi _\text{ m }}{\partial \rho ^\text{ g }}]\!]_{\varGamma } -\frac{\partial \psi _\text{ s }}{\partial [\![\beta ]\!]_{\varGamma }} \right) , \end{aligned}$$

(31)

provided $\text {sign}\beta _{,y}|_{y_\text {g}\pm 0}=\pm 1$, $\text{ sign }\langle \!\langle {\dot{\beta }}\rangle \!\rangle _{\varGamma }=1$ and $\text{ sign }[\![{\dot{\beta }}]\!]_{\varGamma }=1$.

For simplicity, let us assume the symmetry such that $\varkappa _l=\varkappa _u=\varkappa$. Let $\rho _\text {y}$ be the root of the equation $\partial \psi _m/\partial \rho ^\text {g}=\zeta _\text {y}/2\varkappa +\gamma _\text {D}$ which exists and is unique due to the monotonicity of function $\partial \psi _m/\partial \rho ^\text {g}$. It turns out that this dislocation density plays the role of the first threshold. If the dislocation density on the one side of the grain boundary exceeds this threshold value, Eq. (30) can only be satisfied if

$$\begin{aligned} {\left\{ \begin{array}{ll} \rho ^\text {g}|_{y_\text {g}+0}=\varkappa | \beta _{,y}|/ b|_{y_\text {g}+0}=\rho _\text {y}+\rho _\text {j}, \\ \rho ^\text {g}|_{y_\text {g}-0}=\varkappa | \beta _{,y}|/ b|_{y_\text {g}-0}=\rho _\text {y}-\rho _\text {j}. \end{array}\right. } \end{aligned}$$

(32)

Indeed, taking into account that $\rho _\text {j}$ is much less than $\rho _\text {y}$, we can expand function $\partial \psi _\text {m}/\partial \rho ^\text {g}$ into the Taylor series near $\rho _\text {y}$ and neglect terms of order higher than $\rho _\text {j}$. Then it is easy to see that (32) satisfies Eq. (30).

But there is another threshold that does not permit dislocations to reach the grain boundary even if the dislocation density on the one side exceeds $\rho _\text {y}$. The reason lies in the second term on the right-hand side of (31). Since $\psi _\text {s}$ is a function of $|[\![\beta ]\!]_{\varGamma }|$, its derivative is equal to

$$\begin{aligned} \frac{\partial \psi _\text{ s }}{\partial [\![\beta ]\!]_{\varGamma }}=\frac{\partial \psi _\text{ s }}{\partial |[\![\beta ]\!]_{\varGamma }|}\text{ sign }[\![\beta ]\!]_{\varGamma }.\end{aligned}$$

(33)

Because the largest absolute value of this derivative is achieved at $[\![\beta ]\!]_{\varGamma }=\vartheta _0$, with $\vartheta _0$ being the initial misalignment, as long as the first term in the bracket on the right-hand side is less than this value, $[\![\beta ]\!]_{\varGamma }$ cannot increase. To simplify (31), we take into account Eq. (32) and the smallness of $\rho _\text {j}$. Expanding $\partial \psi _\text {m}/\partial \rho ^\text {g}$ into the Taylor series near $\rho _\text {y}$, we transform (31) to

$$\begin{aligned} {}[\![{\dot{\beta }}]\!]_{\varGamma }=\frac{b^2}{\zeta _\text{ a }} \left( \frac{\varkappa }{b}\frac{\partial ^2\psi _\text{ m }}{\partial (\rho ^\text{ g})^2}(\rho _\text{ y}) \rho _\text{ j } -\frac{\partial \psi _\text{ s }}{\partial |[\![\beta ]\!]_{\varGamma }|}\text{ sign }[\![\beta ]\!]_{\varGamma } \right) . \end{aligned}$$

(34)

For the density of interfacial energy, we will consider the modified Read–Shockley’s interfacial energy [46],

$$\begin{aligned} \psi _\text {s}=\frac{\gamma _\text {D}}{4\pi (1-\nu )b}\vartheta \ln \frac{e\vartheta _\text {m}}{\vartheta +\delta }, \end{aligned}$$

(35)

where $\vartheta =|[\![\beta ]\!]_{\varGamma }|$. A small number $\delta$ is added to the denominator of the logarithm to make the derivative of $\psi _\text {s}$ finite at $\vartheta =0$. Finally, we need an equation for $\rho _\text {y}$ which is similar to Eq. (14) for $\tau _\text {i}$. As such, we propose

$$\begin{aligned} \rho _\text{ y }=\rho _{\text{ y }0}+\kappa _\text{ a } |[\![\beta ]\!]_{\varGamma }|. \end{aligned}$$

(36)

4 Application

4.1 Plane constrained shear

Suppose that a bicrystal slab is subjected to plane constrained shear by a given displacement $u_y(t)={\dot{\gamma }}t y$ specified at the upper and lower surfaces (see Fig. 1) and that the system is driven at a constant shear rate $\dot{\gamma }=q_0/t_0$. Let the green colored plane $\varGamma$ be the plane of the grain boundary that is parallel to the (x, y)-plane and is located at $y=y_\text {g}$. This grain boundary divides the bicrystal occupying ${\mathcal {V}}$ into two perfectly bonded single crystals occupying ${\mathcal {V}}^-$ and ${\mathcal {V}}^+$, i.e., ${\mathcal {V}}= {\mathcal {V}}^- \cup {\mathcal {V}}^+\cup \varGamma$. Let the cross section of this bicrystal perpendicular to the z-axis be a rectangle of width c and height h, having the same shape and size over the entire length L. We assume $h\ll c \ll L$ to neglect the end effect near $x=0$ and $x=c$. Edge dislocations can occur during the plastic deformation, and we assume that only one slip system is activated in each single crystal, colored blue and red, respectively. The slip directions are perpendicular to the z-axis and inclined at angles $\phi _u$ and $\phi _l$ to the plane of the grain boundary. Given the geometry of the specimen as well as boundary conditions, we assume that the displacements $(u_x,u_y)$ as well as the plastic slip $\beta$ depends only on y: $u_x=u_x(y)$, $u_y=u_y(y)$, $\beta =\beta (y)$. The displacement $u_z$ vanishes.

Fig. 1

Plane constrained shear of bicrystal

When the plastic deformation occurs, the edge dislocations move on the active slip system with the slip direction ${\mathbf {s}}=(\cos \phi , \sin \phi , 0)^T$ and the normal to the slip plane ${\mathbf {m}}=(-\sin \phi , \cos \phi , 0)^T$. The plastic distortion tensor is given by

$$\begin{aligned} \varvec{\beta }=\beta (y,t) \begin{pmatrix} -\sin \phi \cos \phi &{} \cos ^2\phi &{}0\\ -\sin ^2\phi &{} \sin \phi \cos \phi &{}0\\ 0&{}0&{}0 \end{pmatrix}. \end{aligned}$$

(37)

The plastic and elastic strain tensor becomes

$$\begin{aligned} \varvec{\varepsilon }^\text {p}&=\frac{1}{2}\beta (y,t) \begin{pmatrix} -\sin 2\phi &{} \cos 2\phi &{}0\\ \cos 2\phi &{} \sin 2\phi &{}0\\ 0&{}0&{}0 \end{pmatrix}, \end{aligned}$$

(38)

$$\begin{aligned} \varvec{\varepsilon }^\text {e}&=\frac{1}{2} \begin{pmatrix} \beta \sin 2\phi &{} u_{x,y}-\beta \cos 2\phi &{}0\\ u_{x,y}-\beta \cos 2\phi &{} 2u_{y,y}-\beta \sin 2\phi &{}0\\ 0&{}0&{}0 \end{pmatrix}, \end{aligned}$$

(39)

where the displacement field ${\mathbf {u}}=(u_x, u_y,0)$ is yet unknown. The nonzero components of Nye–Bilby–Kröner’s tensor are $\alpha _{xz}=\beta _{,y}\sin \phi \cos \phi$ and $\alpha _{yz}=\beta _{,y}\sin ^2\phi$. Therefore, the density of the non-redundant dislocation per unit area perpendicular to the z-axis equals

$$\begin{aligned} \rho ^\text {g}=\frac{1}{b}|\varvec{\alpha } \cdot \varvec{e_z }|=\varkappa |\beta _{,y}|/b, \quad \text {where} \quad \varkappa =|\sin \phi |. \end{aligned}$$

(40)

The total energy functional of a crystal becomes

$$\begin{aligned} I=cL\int _0^h [\psi _\text{ e }(\mathbf {\varepsilon }^\text{ e})+\gamma _\text{ D }\rho ^\text{ r }+\psi _\text{ m }(\rho ^\text{ g})+\psi _c] \mathrm {d}{y} + cL \psi _\text{ s }|_{y=y_\text{ g }},\end{aligned}$$

(41)

where $\psi _\text {e}$ is the energy density due to the elastic strain

$$\begin{aligned} \psi _\text {e}=\frac{1}{2}\lambda u^2_{y,y}+\frac{1}{2}\mu (u_{x,y}-\beta \cos 2\phi )^2+\frac{1}{4}\mu \beta ^2\sin ^22\phi +\mu \left( u_{y,y}-\frac{1}{2}\beta \sin 2\phi \right) ^2. \end{aligned}$$

(42)

Note that the angle $\phi (y)$ is piecewise constant:

$$\begin{aligned} \phi (y)= {\left\{ \begin{array}{ll} \phi _\text {l} \quad \text {for } 0<y<y_\text {g}, \\ \phi _\text {u} \quad \text {for } y_\text {g}<y<h. \end{array}\right. } \end{aligned}$$

(43)

Taking the variation of functional I with respect to $u_x$ and $u_y$ at fixed $\beta (y)$, we obtain the equilibrium equations which, after integration, lead to

$$\begin{aligned} \begin{aligned} u_{x,y}&=\gamma +(\beta -\langle \beta \rangle ) \cos 2\phi ,\\ u_{y,y}&=\kappa (\beta -\langle \beta \rangle ) \sin 2\phi , \end{aligned} \end{aligned}$$

(44)

with $\kappa =\frac{\mu }{\lambda +2\mu }$ and $\langle \beta \rangle =\frac{1}{h}\int _0^h \beta \mathrm {d}y$. Inserting Eq. (44) into Eq. (41), we reduce the energy functional to

$$\begin{aligned}&I=cL\int _0^h\Bigg [\frac{1}{2}\mu \kappa \langle \beta \rangle ^2 \sin ^2 2\phi +\frac{1}{2}\mu (\langle \beta \rangle \cos 2\phi -\gamma )^2+\frac{1}{2}\mu (1-\kappa )\beta ^2\sin ^22\phi +\gamma _\text{ D }\rho ^\text{ r } \nonumber \\&\quad +\gamma _\text{ D } k_0 \rho _{ss} \ln \left( \frac{1}{1-\frac{\rho ^\text{ g }}{k_0 \rho _{ss}}} \right) -\chi (-\rho \ln (a^2 \rho )+\rho )/L\Bigg ]\mathrm {d}y+ cL \psi _\text{ s }\Bigl |_{y=y_\text{ g }}. \end{aligned}$$

(45)

Varying this functional with respect to $\beta$ and substituting into (24), we obtain the equilibrium equation

$$\begin{aligned} \tau -\tau _\text {b}-\tau _\text {i}=0, \end{aligned}$$

(46)

with the resolved and back shear stresses being

$$\begin{aligned} \tau =-\mu (\kappa \langle \beta \rangle \sin ^2 2\phi +(\langle \beta \rangle \cos 2\phi -\gamma )\cos 2\phi +(1-\kappa )\beta \sin ^2 2\phi ), \end{aligned}$$

(47)

and

$$\begin{aligned} \tau _\text {b}=-\frac{C_1}{(1-C_2 |\beta _{,y}|)^2}\beta _{,yy}, \quad C_1=\frac{\gamma _\text {D}}{k_0 \rho _{ss} b^2}\sin ^2 \phi , \quad C_2=\frac{1}{k_0 \rho _{ss} b}|\sin \phi |. \end{aligned}$$

(48)

At the top and bottom surfaces of the sample, the boundary conditions $\partial \psi _m/\partial \rho ^\text {g}=\gamma _\text {D}$ must be fulfilled. It is obvious that these imply the Neumann conditions $\beta _{,y}(0,t)=\beta _{,y}(h,t)=0.$ Besides, at $y=y_\text {g}$, the boundary conditions (31) (or (34), alternatively), (32), and (36) must be fulfilled.

4.2 Rescaled boundary value problem

To facilitate numerical integration, we introduce the following rescaled variables and quantities

$$\begin{aligned} \begin{aligned}&{\tilde{y}}=\frac{y}{b}, \quad {\tilde{\rho }}=a^2 \rho , \quad {\tilde{\rho }}^\text {r}=a^2 \rho ^\text {r}, \quad {\tilde{\rho }}^\text {g}=b^2\rho ^\text {g}=\varkappa |\beta _{,{\tilde{y}}}|, \\&{\tilde{\chi }}=\frac{{\bar{\chi }}}{\gamma _\text {D}}, \quad \theta =\frac{T}{T_\text {P}}, \quad {\tilde{\tau }}=\frac{\tau }{\mu }, \quad {\tilde{\tau }}_\text {i}=\frac{\tau _\text {i}}{\mu }, \quad {\tilde{\tau }}_\text {b}=\frac{\tau _\text {b}}{\mu }, \end{aligned} \end{aligned}$$

(49)

where the variable ${\tilde{y}}$ changes from zero to ${\tilde{h}}=h/b$. The plastic slip rate is rewritten in the form

$$\begin{aligned} \frac{q(T,\tau _\text {i},\rho ^\text {r})}{t_0}=\frac{b}{a}\frac{{\tilde{q}}(\theta , {\tilde{\tau }}_\text {i}, {\tilde{\rho }}^\text {r})}{t_0}. \end{aligned}$$

(50)

in which

$$\begin{aligned} {\tilde{q}}(\theta , {\tilde{\tau }}_\text {i},{\tilde{\rho }}^\text {r})=\sqrt{{\tilde{\rho }}^\text {r}}({\tilde{f}}_\text {P}(\theta , {\tilde{\tau }}_\text {i},{\tilde{\rho }}^\text {r})-{\tilde{f}}_\text {P}(\theta , -{\tilde{\tau }}_\text {i}, {\tilde{\rho }}^\text {r})). \end{aligned}$$

(51)

We set ${\tilde{\mu }}_\text {T}=(b/a)\mu _\text {T}=\mu r$ and assume that r is independent of temperature and strain rate. Then

$$\begin{aligned} {\tilde{f}}_\text {P}(\theta , {\tilde{\tau }}_\text {i},{\tilde{\rho }}^\text {r})=\exp \left( -\frac{1}{\theta } \exp \left( -\frac{{\tilde{\tau }}_\text {i}}{r\sqrt{{\tilde{\rho }}^\text {r}}}\right) \right) . \end{aligned}$$

(52)

We define ${\tilde{q}}_0=(a/b)q_0$ such that $q/q_0={\tilde{q}}/{\tilde{q}}_0$, and function ${\tilde{\nu }}$ becomes

$$\begin{aligned} {\tilde{\nu }}(\theta , {\tilde{\rho }}^\text {r}, {\tilde{q}}_0)=\ln \left( \frac{1}{\theta }\right) -\ln \left( \ln \left( \frac{\sqrt{{\tilde{\rho }}^\text {r}}}{{\tilde{q}}_0}\right) \right) . \end{aligned}$$

(53)

The dimensionless steady-state quantities are

$$\begin{aligned} {\tilde{\rho }}_{ss}({\tilde{\chi }})=\exp \left( -\frac{1}{{\tilde{\chi }}}\right) ,\quad {\tilde{\chi }}_0=\frac{\chi _0}{e_\text {D}}. \end{aligned}$$

(54)

Using ${\tilde{q}}$ instead of q as the dimensionless measure of plastic strain rate means that we are effectively rescaling $t_0$ by a factor b/a. We assume that $(a/b)t_0=10^{-12}\,$s. Since the shear rate ${\dot{\gamma }}=q_0/t_0$ is constant, we can replace the time derivative by the derivative with respect to $\gamma$ so that $t_0 \partial /\partial t \rightarrow q_0 \partial /\partial \gamma$. Using the introduced dimensionless variables, the governing equations in the bulk, Eqs. (14), (19) and (46), are transformed to the following system of partial differential equations

$$\begin{aligned} \begin{aligned}&{\tilde{\tau }}-{\tilde{\tau }}_\text {b}-{\tilde{\tau }}_\text {i}=0, \\&\frac{\partial {\tilde{\tau }}_\text {i}}{\partial \gamma }=\cos 2\phi -\frac{{\tilde{q}}(\theta , {\tilde{\tau }}_\text {i},{\tilde{\rho }}^\text {r})}{{\tilde{q}}_0}, \\&\frac{\partial {\tilde{\chi }}}{\partial \gamma }={\mathcal {K}}_{\chi } {\tilde{\tau }}_\text {i} \frac{{\tilde{q}}(\theta ,{\tilde{\tau }}_\text {i},{\tilde{\rho }}^\text {r})}{{\tilde{q}}_0}\left( 1-\frac{{\tilde{\chi }}}{{\tilde{\chi }}_0}\right) ,\\&\frac{\partial {\tilde{\rho }}}{\partial \gamma }={\mathcal {K}}_{\rho } {\tilde{\tau }}_\text {i}\frac{{\tilde{q}}(\theta , {\tilde{\tau }}_\text {i},{\tilde{\rho }}^\text {r})}{{\tilde{\nu }}(\theta , {\tilde{\rho }}^\text {r}, {\tilde{q}}_0)^2{\tilde{q}}_0} \left( 1-\frac{{\tilde{\rho }}}{{\tilde{\rho }}_{ss}({\tilde{\chi }})}\right) , \end{aligned} \end{aligned}$$

(55)

with the rescaled resolved shear stress and back stress being

$$\begin{aligned} \begin{aligned}&{\tilde{\tau }}=-(\kappa \langle \beta \rangle \sin ^2 2\phi +(\langle \beta \rangle \cos 2\phi -\gamma )\cos 2\phi +(1-\kappa )\beta \sin ^2 2\phi ), \\&{\tilde{\tau }}_b=-\frac{k_1}{(1-k_2 |\beta _{,{\tilde{y}}}|)^2} \beta _{,{\tilde{y}}{\tilde{y}}}, \end{aligned} \end{aligned}$$

(56)

where

$$\begin{aligned} k_1=\frac{\gamma _\text {D}}{\mu b^2} k \sin ^2\phi , \quad k_2=k|\sin \phi |, \quad k=\frac{a^2}{k_0 {\tilde{\rho }}_{ss}b^2}. \end{aligned}$$

(57)

During the absorption and transmission process, the density of non-redundant dislocation $\rho ^\text {g}$ near the grain boundary hardly reaches the saturated density $\rho ^\text {g}_{ss}$. Therefore, using the Taylor expansion, we can approximate the defect energy $\psi _\text {m}$ as follows

$$\begin{aligned} \psi _{\text {m}}=\gamma _\text {D}{\tilde{\psi }}_\text {m}, \quad {\tilde{\psi }}_\text {m}=\rho ^\text {g}_{ss}\left( \frac{\rho ^\text {g}}{\rho ^\text {g}_{ss}} +\frac{1}{2}\left( \frac{\rho ^\text {g}}{\rho ^\text {g}_{ss}}\right) ^2\right) . \end{aligned}$$

(58)

By this approximation, the equilibrium equation for microforces (55)$_1$, a stiff second-order quasi-linear differential equation, can be transformed into an elliptic differential equation

$$\begin{aligned} A\beta ''+B\beta +C=0, \end{aligned}$$

(59)

which is numerically more stable than the first one when the jump condition is implemented. The prime denotes the derivative of a function with respect to ${\tilde{y}}$, and

$$\begin{aligned} \begin{aligned} A=\frac{\gamma _\text {D}}{\mu b^2} \frac{a^2}{b^2} \frac{1}{k_0 {\tilde{\rho }}_{ss}} \sin ^2\phi , \quad B=(1-\kappa )\sin ^2 2\phi , \quad \\ C=-(\kappa \langle \beta \rangle \sin ^22\phi +(\langle \beta \rangle \cos 2\phi -\gamma )\cos 2\phi +{\tilde{\tau }}_{\text {i}}), \end{aligned} \end{aligned}$$

(60)

At the external boundaries, Eq. (55)$_1$ is subjected to the boundary condition

$$\begin{aligned} \beta _{,{\tilde{y}}}(0,{\tilde{f}})=0, \quad \beta _{,{\tilde{y}}}({\tilde{h}},{\tilde{f}})=0. \end{aligned}$$

(61)

From (31), we get the rescaled interface condition at the grain boundary as

$$\begin{aligned} \begin{aligned}&\frac{\partial [\![\beta ]\!]_{\varGamma }}{\partial \gamma }=\frac{1}{{\tilde{\zeta }}_\text{ a }} \left( \frac{\varkappa }{2} [\![\frac{\partial {\tilde{\psi }}_{\text{ m }}}{\partial {\tilde{\rho }}^{\text{ g }}} ]\!]_{\varGamma } - \frac{\text{ sign }[\![\beta ]\!]_{\varGamma }}{4\pi (1-\nu _0)}\log \frac{\vartheta _{\text{ m }}}{|[\![\beta ]\!]_{\varGamma }|+\delta }\right) , \\ {}&[\![\beta _{,{\tilde{y}}}]\!]_{\varGamma }=\frac{2}{\varkappa }{\tilde{\rho }}_{\text{ y }}. \end{aligned}\end{aligned}$$

(62)

where ${\tilde{\zeta }}_\text {a}=\zeta _\text {a}t_0/(q_0b\,\gamma _D)$ and ${\tilde{\rho }}_{\text {y}}=b^2\rho _{\text {y}}$. The dimensionless form of Eq. (36) reads

$$\begin{aligned} {\tilde{\rho }}_\text{ y }={\tilde{\rho }}_{\text{ y }0}+{\tilde{\kappa }}_\text{ a } |[\![\beta ]\!]_{\varGamma }|,\end{aligned}$$

(63)

where ${\tilde{\kappa }}_\text {a}=b^2\kappa _\text {a}$.

We also need to pose the initial conditions for $\beta$, ${\tilde{\tau }}_i$, ${\tilde{\rho }}$ and ${\tilde{\chi }}$. For the internal shear stress, we assume that the specimen is not plastically pre-deformed, i.e., ${\tilde{\tau }}_\text {i}({\tilde{y}},0)=0$. As for the plastic slip, two scenarios can be studied: (i) vanishing initial misalignment and (ii) nonzero initial misalignment. The numerical simulations will be performed for the first case where $\beta ({\tilde{y}},0)=0$. The rescaled dislocation density ${\tilde{\rho }}({\tilde{y}},0)$ and the rescaled disorder temperature ${\tilde{\chi }}({\tilde{y}},0)$ are assigned reasonable initial values based on various numerical simulations. It should be noted that, from the physical point of view, these initial values depend on the sample preparation.

4.3 Numerical implementation

We discretize the boundary value problem formulated above by a finite difference method called the line method. The latter replaces the spatial derivatives by finite differences, which allow us to transform partial differential equations into a system of ordinary differential equations with respect to $\gamma$. We illustrate this method for the case ${\tilde{y}}_\text {g}={\tilde{h}}/2$. For the case ${\tilde{y}}_\text {g}\ne {\tilde{h}}/2$, the discretization can be done in a similar way. Let the interval $0< {\tilde{y}}< {\tilde{h}}$ be divided into 2n equal subintervals. Then the step in ${\tilde{y}}$ along the grid becomes $\varDelta {\tilde{y}}={\tilde{h}}/2n$. Because of the jump in $\beta$ at the interface, we distinguish the left and right nodes located at the same point ${\tilde{y}}_\text {g}=n \varDelta {\tilde{y}}$. The coordinates of the nodes on the left interval $0< {\tilde{y}}< {\tilde{h}}/2$ are

$$\begin{aligned} {\tilde{y}}_i=i \varDelta {\tilde{y}},\quad i=0,\ldots ,n, \end{aligned}$$

(64)

while those on the right interval ${\tilde{h}}/2< {\tilde{y}}< {\tilde{h}}$ are

$$\begin{aligned} {\tilde{y}}_i=(i-1)\varDelta {\tilde{y}},\quad i=n+1,\ldots , 2n+1. \end{aligned}$$

(65)

Thus, altogether we have $2(n+1)$ nodes. The spatial derivatives of the plastic slip, $\beta _{,{\tilde{y}}}$ and $\beta _{,{\tilde{y}}{\tilde{y}}}$, for all internal nodes not lying at the external boundaries or the interface, are approximated as

$$\begin{aligned} \frac{\partial \beta ({\tilde{y}}_i,\gamma )}{\partial {\tilde{y}}}=\frac{\beta _{i+1}-\beta _{i-1}}{2\varDelta {\tilde{y}}}, \quad \frac{\partial ^2 \beta ({\tilde{y}}_i,\gamma )}{\partial {\tilde{y}}^2}=\frac{\beta _{i+1}-2\beta _i+\beta _{i-1}}{\varDelta {\tilde{y}}^2}, \end{aligned}$$

(66)

where $\beta _i=\beta ({\tilde{y}}_i,\gamma )$. We decompose the integral $\langle \beta \rangle$ in the resolved shear stress into two integrals that are calculated using the trapezoidal rule. Thus,

$$\begin{aligned} \langle \beta \rangle =\frac{1}{2n}\left( \frac{1}{2}\beta _0+\beta _1+\cdots +\frac{1}{2}\beta _n+\frac{1}{2}\beta _{n+1}+\beta _{n+2}+\cdots +\frac{1}{2}\beta _{2n+1}\right) . \end{aligned}$$

(67)

After these replacements, Eq. (55)$_1$ becomes an ordinary differential equation in all internal nodes.

The first derivatives at the external nodes are replaced by

$$\begin{aligned} \frac{\partial \beta ({\tilde{y}}_0,\gamma )}{\partial {\tilde{y}}}=\frac{\beta _{1}-\beta _{0}}{\varDelta {\tilde{y}}}, \quad \frac{\partial \beta ({\tilde{y}}_{2n+1},\gamma )}{\partial {\tilde{y}}}=\frac{\beta _{2n+1}-\beta _{2n}}{\varDelta {\tilde{y}}}. \end{aligned}$$

(68)

Thus, the boundary conditions (61) become

$$\begin{aligned} \beta _{0}=\beta _{1}, \quad \beta _{2n+1}=\beta _{2n}. \end{aligned}$$

(69)

The first left and right derivatives at the interface are replaced by

$$\begin{aligned} \frac{\partial \beta ({\tilde{y}}_n,\gamma )}{\partial {\tilde{y}}}=\frac{\beta _{n}-\beta _{n-1}}{\varDelta {\tilde{y}}}, \quad \frac{\partial \beta ({\tilde{y}}_{n+1},\gamma )}{\partial {\tilde{y}}}=\frac{\beta _{n+2}-\beta _{n+1}}{\varDelta {\tilde{y}}}. \end{aligned}$$

(70)

while the jump in $\beta$ is

$$\begin{aligned} {}[\![\beta ]\!]_{\varGamma }=\beta _{n+1}-\beta _n. \end{aligned}$$

(71)

These transform Eqs. (62) into two differential algebraic equations.

Since the three remaining equations of (55) do not contain derivatives, they are satisfied node by node. Besides, as the unknown functions are continuous at the interface, we have

$$\begin{aligned} {\tilde{\tau }}_{\text {i}n}={\tilde{\tau }}_{\text {i}n+1},\quad {\tilde{\rho }}_{n}={\tilde{\rho }}_{n+1},\quad {\tilde{\chi }}_{n}={\tilde{\chi }}_{n+1}. \end{aligned}$$

(72)

It is easy to see that the system of partial differential equations and boundary conditions becomes a system of $4(2n+1)$ ordinary differential algebraic equations (DAE) for $4(2n+1)$ unknowns, so the problem is well posed. For the numerical discretization we choose $n=200$, while for the integration of DAEs with respect to the time-like variable $\gamma$ a step size $\varDelta \gamma =10^{-4}$ is chosen. The DAE system is then integrated by the MATLAB integrator ode15s.

4.4 Numerical results

In the numerical simulations, we keep the temperature and shear rate constant at $T=298\,$K and ${\tilde{q}}_0=10^{-13}$. The parameters for the thermodynamic dislocation theory for copper are chosen as presented in [4]:

$$\begin{aligned} r=0.0323, \quad \theta =0.0073, \quad {\mathcal {K}}_{\chi }=350, \quad {\mathcal {K}}_{\rho }=96.1, \quad {\tilde{\chi }}_0=0.25. \end{aligned}$$

(73)

We choose the following parameters for the copper bicrystal

$$\begin{aligned} h=5.1 \, \mathrm{\mu } \text {m}, \, b=0.255\, \text {nm}, \, a=10b, \, \mathrm{\mu }=50 \, \text {GPa}, \, \nu _0=0.34, \, \phi =30^{\circ }. \end{aligned}$$

(74)

The active slip systems of two grains in the bicrystal are constrained to symmetry. The initial conditions for dislocation density and configuration temperature are ${\tilde{\rho }}(0)=6.25 \times 10^{-5}$, ${\tilde{\chi }}(0)=0.23$ and we assume $\gamma _\text {D}=\mu b^2$. The critical density is set as ${\tilde{\rho }}_{\text {y}0}= 1.3 \times 10^{-5}$ and the other two interfacial dissipation parameters are given as ${\tilde{\zeta }}_\text {a}=1.5$ and ${\tilde{\kappa }}_\text {a}=5 \times 10^{-4}$. The parameters for the defect energy and interfacial energy are chosen as follows: $k_0=0.2$, $\vartheta _{\text {m}}=0.192$, $\delta =1\times 10^{-7}$. Unless some parameters are changed for the parameter study, they are kept as default values.

Fig. 2

(a) Distribution of plastic slips at $\gamma =8*10^{-3}$ for $k_0=0.1$, $k_0=0.2$, and $k_0=0.5$. (b) Corresponding hardening behaviors in the pileup process

The saturation density of the non-redundant dislocations $\rho ^\text {g}_{ss}$ as well as the defect energy $\psi _{\text {m}}$ is characterized by the factor $k_0$, so that its value affects the distribution of plastic slip and work hardening. Figure 2a shows the distributions of plastic slip in the pileup process for three different values of $k_0$. When $k_0$ is small, the plastic slip decreases gradually near the grain boundary, while it decreases more steeply to zero for a higher value of $k_0$. Since the slope of plastic slip is proportional to the accumulated dislocation density, these distributions imply that a smaller $k_0$ leads to a larger dislocation accumulation zone and a lower accumulated density of non-redundant dislocations near the grain boundary. This is because a small $k_0$ makes $\rho ^{\text {g}}_{ss}$ small, so a lower density of non-redundant dislocations is required for the backstress to have an effect on the balance of microforces. The hardening behavior in the pileup process is shown in Fig. 2b. It can be seen that the hardening rate is larger for smaller $k_0$. The reason is that the wider dislocation accumulation zone and the slower gradient increase of plastic slip resulting from the small $k_0$ lead to a higher resolved shear stress and consequently a larger hardening rate. A smaller value of $k_0$ also leads to a longer length of the accumulation phase. This is because the accumulated $\rho ^{\text {g}}$ is relatively low for small $k_0$, and therefore, it takes longer to reach the required critical density $\rho _{\text {y}}$.

Fig. 3

(a) Plastic slips during the pileup process, (b) evolution of the misorientation at the grain boundary

At different shear strains $\gamma$, the distribution of plastic slip in the processes of dislocation pileup is shown in Fig. 3a. The Neumann boundary conditions $\beta _{,{\tilde{y}}}(0)= \beta _{,{\tilde{y}}}({\tilde{h}})=0$ cause the curve of plastic slip to remain flat on two sides, while a groove is formed at the position of the grain boundary due to the Dirichlet interface condition. A dislocation pileup zone is located in the center and two dislocation-free zones are located at the boundary layers. In this process, the homogeneous Dirichlet interface condition is applied, and the dislocations cannot penetrate into the grain boundary. Therefore, the groove sinks with increasing strain, but the tip remains attached to the bottom (Fig. 3a). Figure 3b shows the evolution of the misalignment at the grain boundary, which is obtained from Eq. (31). When the non-redundant dislocation density $\rho ^g$ reaches the critical value $\rho _{\text {y}0}$ at $\gamma =0.011$, the evolution equation is activated. The grain boundary starts absorbing dislocations and the misorientation increases only when the strain exceeds the second threshold value at $\gamma =0.0151$, up to which the accumulation process continues. The dislocation absorption process stops at $\gamma =0.0403$ when the slope of the misorientation changes.

Fig. 4

(a) Plastic slips for absorption process, (b) plastic slip omitting the jump

The distribution of plastic slip during the dislocation absorption process is shown in Fig. 4a. In this process, the grain boundary is no longer impenetrable but absorbs dislocations, increasing the misorientation and causing non-redundant dislocations to accumulate. These are two factors that cause the distribution of plastic slip to take a different form than in the pileup process, as shown in Fig. 4a. Here, it can be observed that the plastic slip on two sides of the grain boundary converges not to a single point, but to two points apart, and the distance between them increases, which means that the misorientation of the grain boundary increases. To see the effect of the additionally accumulated non-redundant dislocations, we omit the jump from the plastic slip distribution at $\gamma =0.02$ and $\gamma =0.03$ (solid lines) and compare it with that at the beginning of the absorption process (dashed lines) in Fig. 4b. It can be seen that the tip leaves the bottom, and as the shear stress increases, the discrepancy between the solid and dashed lines also increases, indicating that the groove does not stop sinking.

Fig. 5

(a) Evolution of the mean density of the non-redundant dislocations, (b) evolution of the non-redundant dislocation density at two sides of grain boundary

In the pileup process, the rate of accumulated dislocation density is the highest among the different processes because all dislocations are blocked by the grain boundary. When positive and negative dislocations enter the grain boundary and increase the misorientation in the absorption process, the ability of the grain boundary to block dislocations from transfer is enhanced, and no additional dislocations are accumulated after this process. This behavior is shown in Fig. 5a. In the pileup process, the positive and negative non-redundant dislocation densities are symmetrically distributed on both sides of the grain boundary, while in the absorption process, a jump in the non-redundant density occurs near the grain boundary, as shown in Fig. 5b. Under the interfacial condition Eq. (32), the same amount of dislocation density is shifted from one side to the other.

Fig. 6

Rescaled stress–strain curves: (a) with impenetrable grain boundary (blue), with dislocation absorption process included (red), without dislocation absorption process (black), (b) partition of work hardening

In Fig. 6a, the stress–strain curve of the present model (red curve) is shown, where P1 and P2 indicate the first and second thresholds and P3 the end of the absorption process. For comparison we plot also the stress–strain curves obtained by two other models. The blue curve represents the hardening behavior of bicrystals with impenetrable grain boundary, which blocks all non-redundant dislocations and therefore hardens the most. The black curve describes the hardening behavior of bicrystals without absorption process. When $\rho ^{\text {g}}$ reaches the critical density $\rho _{\text {y}0}$, the pileup process ends at P1 and the traversal process follows. These two curves represent the limits above and below which the hardening with absorption can vary due to different interfacial dissipation. Figure 6b shows the partition of work hardening into isotropic and kinematic hardening. Isotropic work hardening results from the internal stress required by the dislocations to overcome the resistance due to pinning. Kinematic work hardening is caused by back stresses describing the interaction between non-redundant dislocations. In the non-uniform plastic deformation of small size crystals, kinematic work hardening tends to contribute strongly, while in the present study the contribution is small because bicrystals have only one grain boundary where non-redundant dislocations accumulate. It is worth mentioning that for crystals with multiple grain boundaries, the kinematic work hardening rate is higher [10].

Fig. 7

(a) Influence of ${\tilde{\zeta }}_\text {a}$ to the evolution of jump in plastic slip, (b) influence of ${\tilde{\kappa }}_\text {a}$ to the hardening behavior

Among the interfacial dissipation parameters, $\rho _{\text {y}0}$ determines the first threshold and ${\tilde{\zeta }}_\text {a}$ the second, so they determine the onset of the dislocation absorption process. The size of the jump in the plastic slip and the length of the absorption process are controlled by ${\tilde{\zeta }}_\text {a}$, as shown in Fig. 7a. Note that the jump exhibits a linear dependence on the shear strain. The kinematic strain hardening in the absorption process (Fig. 7b) is influenced by the additional accumulated non-redundant dislocations originating from the increased misorientation, so that ${\tilde{\kappa }}_\text {a}$ is the key controlling parameter.

Fig. 8

(a) Distribution of plastic slips in pileup process for $\phi =30^{\circ }$, $\phi =35^{\circ }$ and $\phi =55^{\circ }$, (b) corresponding hardening behaviors

Finally, we show the distribution of plastic slip at $\gamma =0.01$ with $\kappa _a=2\times 10^{-4}$ for $\phi =30^{\circ }$, $\phi =35^{\circ }$ and $\phi =55^{\circ }$ in Fig. 8a. It can be seen that for $0<\phi <45^{\circ }$ and $135^{\circ }<\phi <180^{\circ }$ the plastic slip is positively distributed, while for $45^{\circ }<\phi <135^{\circ }$ it is the reverse. The dashed curve is symmetrical to the black curve about the horizontal coordinate axis. We see that the blue and dashed curves coincide in the dislocation-free zone. This is because the resolved shear stress $\tau$ and the residual stress $\tau _{\text {i}}$ for $\phi =35^{\circ }$ are identical in magnitude but opposite in sign to those for $\phi =55^{\circ }$. The distributions of plastic slip are different near the grain boundary because the accumulation of non-redundant dislocation density is inconsistent due to the different slip systems. Figure 8b shows the strain hardening behavior for three different $\phi$, where the asterisks indicate the beginning and end of the absorption processes. It can be seen that the slip system affects not only the hardening rate but also the length of the absorption process $l_{\gamma }$, as indicated in the bar graph.

5 Conclusion

In the framework of thermodynamic dislocation theory for non-uniform plastic deformation, we have developed the interface conditions for various interactions of dislocations with the grain boundary, e.g., dislocation accumulation, absorption, and transfer. Interfacial dissipation in terms of rates for the mean and jump of plastic slip has been proposed. An evolution equation for grain boundary misorientation has been developed. The proposed interface dissipation yields two thresholds for the dislocation absorption process. The dislocation absorption leads to asymmetry in the distribution of plastic slip and non-redundant dislocation density on two sides of the grain boundary. It turns out that due to the additional accumulated density of non-redundant dislocations, kinematic work hardening also develops during the process of dislocation absorption. The extent of misorientation and the length of the dislocation absorption process are determined not only by the interfacial dissipation, but also by the slip system of the grains.