Application of rigorous interface boundary conditions in mesoscale plasticity simulations

The interactions between dislocations and interface/grain boundaries, including dislocation absorption, transmission, and reflection, have garnered significant attention from the research community for their impact on the mechanical properties of materials. However, the traditional approaches used to simulate grain boundaries lack physical fidelity and are often incompatible across different simulation methods. We review a new mesoscale interface boundary condition based on Burgers vector conservation and kinetic dislocation reaction processes. The main focus of the paper is to demonstrate how to unify this boundary condition with different plasticity simulation approaches such as the crystal plasticity finite element (CPFEM), continuum dislocation dynamics (CDD), and discrete dislocation dynamics (DDD) methods. In DDD and CDD, plasticity is simulated based on dislocation activity; in the former, dislocations are described as discrete lines while in the latter in terms of dislocation density. CPFEM simulates plasticity in terms of slip on each slip system, without explicit treatment of dislocations; it is suitable for larger scale simulations. To validate our interface boundary condition, we implemented simulations using both the CPFEM method and a two-dimensional CDD model. Our results show that our compact and physically realistic interface boundary condition can be easily integrated into multiscale simulation methods and yield novel results consistent with experimental observations.


Introduction
Interfaces play a central role in the plastic deformation of polycrystalline materials.Plastic deformation commonly arises from the motion of lattice dislocations, while grain boundaries (GBs) and other interfaces act as barriers to the motion of dislocations within grains.Dislocation pile-ups near interfaces can induce back stresses within the grains leading to significant strengthening.However, introducing more interfaces into the materials can lead to the classic strength-ductility trade-off; i.e., where the materials become stronger but more brittle (Lu et al., 2004;Zhou et al., 2019).To overcome this dilemma, materials with nanoscale coherent twin boundaries or heterostructured materials have been the focus of much recent research (Lu et al., 2009a,b;Zhu and Wu, 2023).The improved performance of these materials may be attributed to tuned dislocation-interface interactions.For example, coherent twin boundaries can both block dislocations from within the grain as well as serve as dislocation glide planes; stress concentrations near dislocation slip plane-twin boundary intersections can be partially relaxed by such glide.This demonstrates the importance of a quantitative understanding of the interactions between interfaces and plasticity in the grain interior as a means of optimizing strength and ductility through interfacial engineering.
In situ transmission electron microscopy (TEM) observations show that dislocations from the grain interior may be absorbed into GBs and transformed (partially) into a set of disconnections on the GB (Kondo et al., 2016;Malyar et al., 2017).These direct observations provide guidance on how dislocations from the grain interior interact with interfaces.Based on numerous experimental results, the phenomenological Lee-Robertson-Birnbaum (LRB) criteria were proposed (Lee et al., 1989(Lee et al., , 1990b,a) ,a) to predict the tendency for slip transmission across GBs.However, TEM observations primarily capture individual dislocation-GB interaction events.Nonetheless, it remains a challenge to quantitatively predict the degree of slip transfer, disconnection activity (giving rise to GB sliding), dislocation pile-up, and generation of dislocations on complementary slip systems (giving rise to local hardening) for any particular GB.
Atomistic simulation (molecular dynamics, MD) can capture dislocation/interface interaction with atomic resolution (Dewald andCurtin, 2006, 2007;Zhang et al., 2014), consisted with in situ TEM observations as well as help determine parameters required for larger length and time scale simulation methods (Aragon et al., 2022).However, severe limitations on the length and time scales accessible to MD make them ineffective for simulating deformation at typical experimentally-relevant scales.This implies the need to perform mesoscale simulations with methods such as discrete dislocation dynamics (DDD) (LeSar and Capolungo, 2020), continuum dislocation dynamics (CDD) (El-Azab and Po, 2020) and crystal plasticity finite element method (CPFEM) (Roters et al., 2010).
Classically, most mesoscale simulation methods used to describe plasticity in polycrystalline materials view GBs as blockers of dislocation motion.The corresponding GB-interface boundary condition (BC) for dislocation density evolution is a Neumann BC (i.e., no dislocations crossing the GBs); e.g., see (Schulz et al., 2014;Dunne et al., 2007;Jiang et al., 2019).However, as discussed above, the interactions between lattice dislocations and interfaces is not this simple.Here we examine the application of a more robust interface BC for dislocation/interface interactions in the context of different mesoscale simulation methodologies.
In continuum dislocation dynamics (CDD), the fundamental degrees of freedom are a set of fields defined throughout the simulation region; the specific choice of field variables differs amongst CDD methods.Achaya proposed the field dislocation mechanics (FDM) based on geometrically necessary dislocations (GNDs), as described through the Nye tensor α (Acharya, 2001(Acharya, , 2003)).The dislocation velocity field depends on both the GND and statistically stored dislocation (SSD) densities; this is addressed phenomenologically.The corresponding grain boundary BCs are formulated by balancing Burgers vector content in patches of a finitely deforming body (Acharya, 2007).Hochrainer et al. (Hochrainer et al., 2007) proposed a higher-dimensional model which accounts for both the dislocation density and line-curvature fields, later simplifying the approach through phenomenological approximations for the evolution of the average curvature (Hochrainer et al., 2014).More recently, another CDD was proposed in which dislocations are represented by vector fields on each slip system (Leung et al., 2015;Xia and El-Azab, 2015).The local dislocation line direction is incorporated into the dislocation density and the dislocation density evolution laws at each material point to close the governing equations.GBs are considered through traction and dislocation flux boundary conditions (El-Azab, 2000;Xia and El-Azab, 2015).In this approach, GBs are no longer impenetrable to dislocations and include Burgers vector fluxes through a source term in the governing dislocation dynamics equations.Although this type of interface BC accommodates such issues as Burgers vector conservation, it is not simply implemented in the corresponding methods and is not compatible with the other CDD methods.
Crystal plasticity models describe plasticity on a larger scale, providing a means for direct comparisons with macroscopic experimental observations and are commonly implemented within finite element methods, CPFEM (Roters et al., 2010).Alternative implementations, including the micromechanical self-consistent method (Lebensohn and Tomé, 1993) and the full-field Fourier-based method (Lebensohn et al., 2012).At their core, crystal plasticity models, invoke a constitutive law that (in its local formulation) determines the plastic strain rate at each material point based on the current local material state and the applied load.CPFEM typically models dislocation-GB interactions through BCs that differ depending on whether the constitutive laws incorporate strain gradients.The first type of interface BC introduces an additional slip resistance into the rate equation for plastic slip within the framework of conventional CPFEM (Ma et al., 2006a,b;Mayeur et al., 2015).Here, special interface elements with larger slip resistance than those within grains are introduced; the proposed additional resistance is based on semi-empirical, rather than physical (dislocation reaction) considerations.The second type of interface BC is developed within the framework of strain gradient theory.The free energy of slip transfer at the interface is constructed by introducing an additional interface potential to the potential energy functional of the total system (Gudmundson, 2004;Aifantis et al., 2006;Aifantis and Willis, 2005).By setting the first variation of the potential energy functional to zero, jump conditions for high order tractions across interfaces are obtained; this method can be extended by incorporating more elaborate interface free energy functionals incorporating plastic strain jumps and the average plastic strain across the interface (Fleck and Willis, 2009).Such BCs do not explicitly account for crystallographic misorientation across the GBs.Gurtin proposed a quadratic free energy form for the GB that is essentially an interface defect tensor (Gurtin, 2008).The interface defect tensor measures the GND density at the GBs, analogous to the Nye tensor within the grain interior.The flow rules at the interface are derived based on the second law of thermodynamics in the form of a free-energy imbalance.This interface BC couples grain interiors to GBs by establishing a microscopic force balance between them.Compared with the above methods, this BC can explicitly account for GB crystallography through slip-interaction moduli.This BC was numerically implemented within a strain gradient crystal plasticity framework ( Özdemir and Yalçinkaya, 2014;Van Beers et al., 2013).These treatments of GBs can be used for the continuum description of dislocation-GB interactions within a thermodynamically consistent crystal plasticity framework, although the interactions should be viewed as a kinetic rather than a thermodynamic process.Note that the second type of BCs is only applicable to strain gradient crystal plasticity and is not compatible with other methods.
The discrete dislocation dynamics (DDD) approach for crystal plasticity tracks the motion of individual dislocation lines within an ensemble of dislocations evolving (Bulatov and Cai, 2006).A two-dimensional (2D) DDD model was developed to simulate grain boundary (GB) sliding and the transmission of lattice dislocations across GBs by introducing a Frank-Read source into the grains adjacent to the GB (Quek et al., 2014(Quek et al., , 2016)).While providing valuable insights into the importance of GB sliding during plastic deformation, this model is limited to 2D and relatively small numbers of discrete dislocations.Other phenomenological laws based upon the LRB criteria have been proposed to account for individual dislocation-GB interaction in 3D (Cho et al., 2020;Zhang et al., 2021).In contrast to the treatment of dislocations as "points" in 2D DDD, dislocations are modeled as discrete line segments gliding along slip planes under a driving force in 3D DDD.The kinetic process for such dislocation-GB interactions can be simulated based on prescribed empirical rules.For example, the process of dislocation transmission through GBs was achieved through a Frank-Read source dislocation bowing out model in GB adjacent grains when the resolved shear stress on the incoming dislocation exceeded a GB transmission strength (Zhou and LeSar, 2012).Another dislocation-GB interaction model assumed that dislocations can be absorbed into GBs once the resolved shear stress reached a critical value (Zhang et al., 2021) and dislocations can be emitted from GBs when specific rules (dependent on the residual Burgers vector and the resolved shear stress) are satisfied.Although these dislocation-GB interaction models are intuitive and easily implemented within DDD schematics, the empirical criteria are not rigorous and not compatible with CDD methods or the strain gradient plasticity finite element method.
The dislocation-interface BCs discussed above tend to be limited to specific simulation methods.In addition, most of these BCs rely on empirical criteria and are not based on rigorous microscopic mechanisms that quantitatively describes reactions between lattice dislocations and interfaces.Our proposed mesoscale interface boundary condition (BC) is more general and is based upon the kinetics of dislocation reactions at interfaces (Yu et al., 2023).Its compact form allows for easy combination with multiscale simulation methods.The main challenge is how to relate the input (dislocation density) and output (dislocation flux) quantities to the field variables of the simulation.In this paper, we first review the mesoscale interface BC in Section 2 and demonstrate the basic idea of its application.We then discuss how to implement this interface BC with different CDD simulation methods in Section 3 and CPFEM in Section 4, respectively.For each simulation method, we provide numerical examples that demonstrate the validity of our interface BC.Finally, we discuss how our interface BC can be naturally extended to DDD in Section 5.

The Mesoscale Interface Boundary Condition
When dislocations from the grain interior approach a GB, they tend to be blocked, pile up, and or absorbed into the GB.In order to release the accumulated energy of the GBs, dislocations may be emitted from the GBs into the same (reflection) or adjacent (transmission) grain (Kacher and Robertson, 2012;Kacher et al., 2014;Kacher and Robertson, 2014).These phenomena are related through dislocation reactions at the interface as shown in the schematic Fig. 1.We explicitly consider an incoming dislocation slip system in one grain, outgoing slip systems in the original and opposite grains, and slip along the interface.We briefly recall the formulation of the mesoscale interface BC (Yu et al., 2023).
We first label the slip system k by superscripts "(k)", such that the slip plane unit normal and slip direction are denoted n (k) and s (k) , respectively (l (k) ≡ n (k) × s (k) is defined to form a local Cartesian (3) ( 4) ξ (1)   n I

Interface
Figure 1: Schematic of single reaction process between the interface and lattice dislocations.
coordinate system for each slip system).When the dislocations approach an interface, the dislocation line ξ (k) will align with intersection slip plane-interface intersection, ξ ≡ n (k) ×n I (see Fig. 1).In a continuum, the dislocation density ρ (k) at the interface is defined as the number of dislocation lines threading a unit area with unit normal ξ (k) I .The dislocation flux J (k) is the number of dislocation lines that cross a unit area with unit normal ξ (k) I per unit time.To conserve Burgers vector, we balance all four Burgers vector fluxes J (k) (k = 1, 2, 3, 4) through a single reaction process (all fluxes are considered positive away form the slip plane-interface intersection should be concealed out at the interface): (1) Applying linear response theory, we write the dislocation fluxes as linearly proportional to the driving force (i.e., the Peach-Koehler force exerted) acting on the dislocation density.Combining this with the conservation of Burgers vector condition Eq. ( 1) allows us to write the interface BC as where the subscript "I" represents quantities evaluated at a point on the interface, the Burgers vectors are 4) , the resolved shear stresses (RSS) are and 3) , ρ (4) T are the generalized dislocation fluxes and densities.The matrix c ≡ c (234) , c (314) , c (124) , c (132) T only depends on the crystallographic orientation of the slip systems meeting at the interface and c (ijk Finally, the scalar quantity κ is a reaction constant which is related to the microscopic characteristics of the interface and depends on the detailed atomic structure and bonding at the interface and in dislocation cores. We now consider the most general case where there are N (N > 4) slip systems (see Fig. 1).An arbitrary reaction (denoted "(n)") occurs among any four of the N slip systems.Similar to Eq. ( 2), the kinetics of Reaction "(n)" can be described by κ (n) is similar to the κ in Eq. ( 2) and represents the reaction constant for Reaction "(n)" (each reaction will, in general, have a unique reaction constant).There are relations similar to Eq. ( 3) for each reaction.The total Burgers vector flux, due to all combinations of four reactions where This general form of the interface BC is applicable to reactions involving any number of slip systems.
In our interface boundary condition (BC) Eq. ( 4), the input and output quantities are the dislocation density ρ (k) and dislocation flux J (k) on each slip system, respectively.The main objective is to establish a connection between these two quantities using field variables appropriate for the different continuum plasticity simulation methods described above.
The first step is to determine the dislocation density at the interface and treat it as an input to our interface BC, as illustrated in Fig. (1).Subsequently, the output dislocation flux J (k) on each slip system is used to update the plastic strain rate for the CPFEM or the dislocation flux for the CDD method at the interface.The proposed interface BC can be employed to determine the outgoing slip system based on maximum energy dissipation rate in the DDD method.Details of the application of the interface BC to various plasticity methods are discussed below.

Application to Continuum Dislocation Dynamics
Growing interest in manipulating material properties through microstructure control, grain boundary engineering and interest in nanoscale devices through miniaturization of devices, has increased the community's focus on physically motivated, dislocation-based continuum theories of plasticity (Acharya, 2001;Leung and Ngan, 2016;Hochrainer et al., 2007).In recent years, several continuum approaches have been introduced that can simulate the motion of curve dislocation lines.As described above, several continuum dislocation dynamics (CDD) methods have been developed, based on different choices of field variables, to describe the coarse-grained dislocation density field.We consider three such CDD schemes here and show how to combined these with our general interface boundary condition.
Field dislocation mechanics, a continuum dislocation plasticity approach proposed by Achaya and coworkers (Acharya, 2001(Acharya, , 2003)), focuses on the Nye tensor α (dislocation density tensor) on each slip system as the primary state variable.The Nye tensor is a local measure of the geometrically necessary dislocation (GND) density: where b (k) represents the Burgers vector of k th slip system and ρ k) is the density of dislocations on that slip system with unit tangent vector ξ (k) , i.e., ρ (k) .The Nye tensor for slip system k th evolves as where v (k) is the dislocation velocity along its local normal ξ (k) × n (k) and g (k) is a source term on that slip system.The dislocation density where this slip plane intersects the interface plane is I .We can substitute this dislocation density ρ (k) at the interface into the interface BC, Eq. ( 4) and ρ (k) v (k)  in Eq. ( 6) can be updated from the dislocation fluxes J (k) calculated using the interface BC.Note that the updated dislocation flux only applies at the interface while the dislocation flux within the grain may be calculated using any appropriate bulk kinematic law.
The continuum dislocation-density function dynamics (CDDFD) approach of Ngan and co-workers (Leung and Ngan, 2016; Ngan, 2017) considers the dislocation line orientation (relative to the Burgers vector).Let ϱ (k) (r, θ) be the density of dislocations with line direction e (k) θ = cos θs (k) +sin θl (k) at r on slip system k.ϱ (k) (r, θ)dθ represents the dislocation density with orientations in the range θ to θ + dθ at r.The dislocation character density ϱ (k) (r, θ) evolves as where v (k) represents the dislocation velocity in direction e represents the dislocation rotation rate.Within the grain, the dislocation velocity v (k) can be calculated using an appropriate driving-force velocity relation.Once the v (k) field is obtained, the rotational velocity is evaluated as v θ , where e (k) θ is the dislocation line direction of orientation/character θ.When the dislocations move to the interface r I , the dislocation line will rotate to be collinear with the slip plane/interface intersection ξ (k) I as shown in Fig. 3.The dislocation density at the interface on slip system k is ϱ has orientation/character θ I .Substituting this dislocation density into the interface BC Eq. ( 4), we obtain the dislocation flux, Hochrainer and coworkers develop a higher-dimensional continuum dislocation dynamics (hdCDD) that includes the dislocation curvature k(x, φ) as a second fundamental field variable in addition to dislocation density ρ(x, φ) (Hochrainer et al., 2007).To illustrate the main idea, consider a two-dimensional case where the dislocations have the same Burgers vector b (k) moving on k th slip plane.ρ (k) (x, φ) represents the dislocation density at x on plane k, φ is the angle between the dislocation line direction and Burgers vector b (k) (similar to θ (Fig. 3) in the dislocation-density function dynamics approach (Leung and Ngan, 2016;Ngan, 2017)).
The original evolution equations for the scalar fields, dislocation density ρ (k) (x, φ) and curvature k (k) (x, φ), have proven to be computationally costly for simulations.This led to the development of simplified approaches based on Fourier expansions of the local variables ρ (k) (x, φ) and curvature k (k) (x, φ) (Hochrainer et al., 2014;El-Azab and Po, 2020).The zeroth-and first-order coefficients of ρ (k) (x, φ) in the Fourier expansion corresponding to the total scalar dislocation density ρ

Slip plane
Figure 3: Schematic of dislocation-density function dynamics The classical Nye tensor is simply recovered as k) .The evolution of the scalar dislocation density and curvature fields can be expressed as where n (k) is the unit normal of slip plane , κ ⊥ ≡ κ (k) × n (k) and k(k) is the average curvature.A phenomenological evolution equation for average curvature k(k) is required to close the governing equations Eq. ( 9).Since it is not directly related to our interface BC's application we will not discuss it further.
Similar to the FDM of Achaya in Eq. ( 6), the dislocation density at the interface is αξ k) .The dislocation flux magnitude is obtained by substituting this into our interface BC Eq. ( 4).Again, the fluxes ρ ⊥ in Eq. ( 9) at the interface can be updated by J (k) .Since the dislocation lines align with the interface/slip plane intersection line direction ξ (k) I (at the interface), the dislocation curvature at the interface should be set to k(k)

Numerical simulations
As shown above, the interface BC Eq. ( 4) can be easily implemented within different CDD methods; we now present such an application based upon a relatively simple two-dimensional (2D) CDD model.
Consider the lamellar configuration illustrated in Fig. 4(a) (periodic along the x and y-directions).There are two phases, α and β, of width L α and L β , separated by interfaces (denoted by black dashed lines) in each period.The slip direction and plane normal are s (k) and n (k) for (k) = α or β.The dislocation density vector is defined as k) is the dislocation line direction).All dislocation lines are straight, perpendicular to the plane, in the two-dimensional model and its direction is always; i.e, where k) is the angle between the slip direction s (k) and e x .
Defining the dislocation flux as k) , the evolution of dislocation density simply becomes the dislocation flux divergence: ρ The magnitude of the dislocation velocity may be described by a power law: where dislocation densities of opposite signs are denoted by subscripts "±", τ * (k) is the slip resistance for k th slip system, and the constant n depends on the range of stress.The dislocation density satisfies the balance: ρ(k) + ρ(k),gen +/− .As proposed by Arsenlis and co-workers (Arsenlis et al., 2004), dislocation annihilation occurs when opposite signed dislocations come within a capture radius r c ; i.e., the annihilation rates are ρ(k),ann When a stress is applied, a dislocation pair (opposite signs) is emitted from a source at the rate (Kocks et al., 1975) ρ(k),gen , where η and m are constants.The net dislocation density is ρ − .The total shear stress τ (k)  in Eq. ( 11) has contributions from both external σ ext and internal (associated with all other dislocations) σ int stresses.In the simulations, we employ reduced variables: ρ ≡ ρL 2 , x ≡ x/L, t ≡ v * t/L, ṽ ≡ v/v * , τ ≡ τ /K, where L ≡ L α + L β is the width of the whole cell in x and K ≡ µ/[2π(1 − ν)].For simplicity, we omit the tilde in the reduced quantities below.The parameters for these simulations are listed in Table .1; since we assume all slip systems have the same properties, the superscript "(k)" representing different slip systems are omitted.
We first consider a simple case where there is only one slip system within each phase (θ α = 30 • and θ β = 15 • ) under an external shear stress σ xy with one operable dislocation source at the centre of phase α. Figure 4(a) shows the equilibrium dislocation density for the case in which the interface is impenetrable (i.e., a Neumann BC); the reaction constant κ (n) = 0 in Eq. (3).Here, the dislocations glide along the slip plane and pile up near the interface.When dislocation reactions can occur at the interface (κ (n) ̸ = 0), with the interface boundary condition Eq. ( 4) applied, the model approaches a state as shown in Fig. 4(b).For this case (constant external stress with a Robin BC), there is no equilibrium or steady state; i.e., the dislocation point source continues to operate such that dislocations continue to cross the GBs.In this case, dislocations flow across the interface from the α phase into β. the interface is much larger in the case of the impenetrable interface than when the Reaction BC is applied; the Reaction boundary condition makes the interface penetrable.When the interface is impenetrable, the total resolved shear stress (RSS) at the dislocation source in α tends to zero as a result of the back stress from the pile up cancelling the stress at the source -effectively shutting down the dislocation source.On the other hand, when the Reaction BC is applied at the interface, the penetrability of the interface makes the pile up weaker and unable to cancel the applied stress -hence the stress at the source asymptotes to a non-zero value.Figure 4(d)-(f) illustrate the evolution of the dislocation distribution for a case where there are two slip systems in α phase (θ α1/α2 = ±30 • ) and one in β (θ β = 15 • ).The simulation results demonstrate that the dislocation pile up from the α 1 slip system not only can transmit into β phase, but also can be reflected back into α 2 slip system.(Note that there is an effect of the periodic BC in the simulation cell.)These simulation results are in excellent agreement with experimental observations.To illustrate the generality of our interface BC, we consider a more complex example involving three grains (denoted α, β, γ).Each grain has two slip systems with the following orientations: In these simulation, the external stress is σ xy .The dislocation sources in β are less active since the resolved shear stress on the β slip system is small relative to the other grains.For the same reason, it is difficult for slip to propagate from α and γ into β; see Fig. 5(a).In this scenario, only a small number of dislocations react at GBs and most of the slip is confined to the individual grains.The primary mechanism responsible for the plastic deformation of this tricrystal is the motion of dislocations within the grains, as depicted in Fig. 5(b).

Application to the Crystal Plasticity Finite Element Method
Crystal plasticity methods are robust, widely-used computational tools for examining the relationship between mechanical properties and material structure which have proven especially effective in micromechanics applications, including strain hardening in single crystals and texture evolution in polycrystalline aggregates (Roters et al., 2010).Crystal plasticity models are based upon the idea that plastic flow occurs via slip or shear on activated slip systems characterized by a critical resolved shear stress and/or a hardening law.Within the framework of small deformation, the plastic distortion is the sum of slip γ (k) on individual slip systems: Conventional CPFEM has no description of dislocation physics.A Burgers tensor can be defined to link dislocation slip to the plastic strain gradient (Gurtin et al., 2010): where G T e represents the Burgers vector (per unit area) for infinitesimal closed circuits on any plane with unit normal e.For any slip system k, s (k) and l (k) ≡ n (k) × s (k) form an orthonormal basis for slip plane Π (k) .Since ∇γ (k) × n (k) is orthogonal to n (k) , it can be expanded in terms of s (k) and l (k) , Hence we can rewrite the Burgers tensor as where k) is the Burgers tensor for slip system k.
(G (k) ) T e is the contribution of slip system k to the Burgers vector (per unit area) with unit normal e.When a dislocation contacts the interface, it aligns with the interface-slip plane intersection ξ (k) I .The total dislocation density of slip system k near the interface can be expressed as where the Burgers vector and slip direction s (k) have the same directions.Substituting Eq. ( 14) into Eq.( 15), the left hand side of Eq. ( 15) can be written as where we decompose the unit line direction vector as ξ k) .Comparing Eq. ( 15) and Eq. ( 16), we find the dislocation density of k th slip system near the interface as We relate the dislocation flux to crystal plasticity variables by taking the time derivative of Eq. ( 17) where the dislocation flux is ).The magnitude and direction of the flux l l (k) (i.e., perpendicular to the dislocation line ).Then the interface BC Eq. ( 2) can be expressed as where the dislocation density on each slip system is obtained by Eq. ( 17) from which the plastic strain rate on each slip system is determined.The interface BC implies that the slip systems are more likely to accumulate plastic deformation when there is a large plastic strain gradient in the corresponding slip system near the interface.The plastic strain rate on slip system k is γ k) .In accordance with the conservation of Burgers vector condition Eq. ( 1), the total net plastic strain rate at each reaction point must be zero.This implies that our interface BC implicitly incorporates material compatibility at the interface during dislocation reactions.

Numerical simulations
We incorporate crystal plasticity as follows.The slip rate on slip system α, γα , is formulated as a function of the resolved shear stress τ α and the critical shear stress: γα = f (τ α , τ α c ) and the evolution of the critical shear stress is dependent on the total shear γ and shear rate γα : τ α c = g(γ, γ).We employ the rate-dependent kinetic law for FCC metal slip systems proposed by Hutchinson (Hutchinson, 1976): where γ0 and m are material parameters that determine the reference shear rate and rate sensitivity.The strain hardening is characterized by the evolution of the strength through the incremental relation: where h αβ are the slip hardening moduli and the sum is over n activated slip systems.Here h αα and h αβ are the self-and latent-hardening moduli.A simple form for the slip hardening moduli is (Peirce et al., 1982) where h 0 and τ 0 are the initial hardening modulus and yield stress, τ s is the stage I threshold stress and γ is the Taylor cumulative shear strain on all slip system.q αβ is a measure of latent hardening; it is commonly set to 1.0 for coplanar slip systems α and β, and between 0 and 1 for other slip systems.Other types of models for slip hardening by replacing Eq. ( 20) with other expressions to describe, for example, three stage hardening of crystalline materials (Bassani and Wu, 1991;Zarka, 1975).This crystal plasticity model was incorporated into an Abaqus user material subroutine (UMAT) and applied to simulate the deformation of a 100 mm long bar of square cross section (20 mm × 20 mm) of FCC copper with a set of activated {111}⟨110⟩ slip systems (see Table 2).The bar is a bicrystal (Grains α and β); see Fig. 6(a).The GB is modeled as a thin region with a single {100}⟨001⟩ slip system.The slip rate of the GB elements is calculated through the interface BC: γI = J I b I = ρ I v I b I .The copper bicrystal bar is subject to uniaxial tensile stress of 200 MPa.The displacements of all nodes at the left end plane of the bar are constrained as shown in Fig. 6(a).In this example, we assume that all dislocation reactions have the same rates κ (n) in Eq. (4).
Figure 6 shows the von Mises stress, dislocation density for one of the slip systems, and displacements.The left and right columns correspond to the Reaction (Robin) (κ (n) ̸ = 0) and Impenetrable interface (κ (n) = 0) BC cases.As depicted in Fig. 6(b), the von Mises stress contour reveals stress concentrations near the grain boundary (GB).Figure 6(c) demonstrates the accumulation of dislocations with opposite signs near the GB for one of the activated slip systems.Figure 6(d1) shows that the displacement field is discontinuous across the Reaction BC GB due to GB sliding.Comparing the Reaction and Impenetrable cases, we see that dislocation reactions near the interface relax the high stress at the interface (Impenetrable BC case) by allowing some dislocation transmission.This is responsible for the smaller von Mises stress concentrations near the Reaction BC GB than in the Impenetrable BC case.It is especially interesting to note that the Reaction BC results in GB sliding, even when the bar is axially loaded, as shown in Fig. 6(e).Because the Reaction BC allows for some slip transfer across the GB than when the GB is impenetrable, the Reaction BC bicrystal shows a lower yield strength and is less hardening than when it is impenetrable, as illustrated in Fig. 6(f).

Application to Discrete Dislocation Dynamics
When an individual dislocation approaches an interface, it may be blocked, transmitted, or reflected onto another slip system.In most DDD simulations, dislocations were able to move on new slip systems according to empirical, geometric, or thermodynamic rules -often drawn from experimental or molecular simulation observations (Lu et al., 2019;Aragon et al., 2022).Here, we explicitly account for dislocation reactions at the interface based on our Reactive BC interface model and all possible slip systems.In DDD simulations, each dislocation is discrete in the reaction process.Assuming that there is single dislocation line from slip system i approaching the interface as shown in Fig. 7(a) .The dislocation density of incoming slip system is ρ (i) = 1/A e , where A e is the unit area near the interaction.There are many potential outgoing slip system from which to choose, corresponding to the many possible reactions amongst the different slip systems.For an arbitrary reaction, denoted by "n", we can substitute dislocation density ρ (i) into the interface BC Eq. ( 2) and obtain the dislocation fluxes of reaction n, J n = (J The fluxes of all (C 3 N −1 ≡ m) possible reactions with slip system i involved can be calculated following the same procedure (1) , J (1) , J (1) , J (1) ) (J (2) , J (2) , J In an energetically favorable structure, the nucleated dislocation should form such that it dissipates the store elastic energy at the fastest rate (this is equivalent to the maximum rate of entropy production) (Cho et al., 2020).We then determine the outgoing slip system following two steps: 1. Calculate the reaction rates of all possible reactions.
2. Choose the reaction with the maximum power dissipation amongst all possible reactions.Since the interface Reaction BC is based upon linear response theory, it corresponds to evolution in the direction that maximizes the rate of entropy production.The magnitude of the local entropy production rate can be expressed as the product of the driving force and dislocation flux for each slip system (Yu et al., 2023).The local entropy production rate for reaction "n" is the sum of the entropy production rate over all four slip systems involving in this reaction ( ṡ . The dominant reaction is that "n" with the largest entropy production rate (amongst all reactions ṡ(n) ) and the outgoing system k with the largest entropy production rate ṡ(k) (n) (amongst the three "outgoing" slip systems).
Once the outgoing slip system is determined, the dislocation line MN on slip system i rotates by β to align with slip plane k as Fig. 7(b) shows.The residual Burgers vector left at the interface is b r = b (i) − b (k) .The velocity of the emitted dislocation on the outgoing slip plane k is v k) , where dislocation density on slip system k is ρ (k) = 1/A e , the constant unit area A e cancels.The new segment will bow out along slip plane k as shown in Fig. 7(c) .The new node will be created and the distant between line MN and new node P is d (k) = v (k) dt, where dt is the time step of the iteration.

Discussion
In the simulations discussed above, we made the simplifying assumption that the interface is flat, while interfaces are often curved.To address this, the distribution of disconnections (line defects with Burgers vector and step character) in the interface should be considered (Han et al., 2022).Specifically, both the shape of any interface and the relative displacements between the phases/grains can be rigorously described in terms of the distribution of disconnection steps and Burgers vectors, respectively.Since the step heights are not involved in interfacial Burgers vector reactions, the interface boundary condition (BC) can still be applied to a curved interface by disregarding the disconnection steps at each point.However, some modifications are necessary when applying the interface BC to a curved interface since its fundamental character varies from point-to-point.First, the geometry term (c ⊗ c) in Eq. ( 2) is related to {s I , n I }, and its tangential and normal directions will change along the curved interface plane.Second, the driving force along the interface, τ I , and the reaction constant, κ, may vary with position along the curved interface.
The variation of reaction constants is necessary for different dislocation reactions.In our theory, the reaction constant reflects the energy dissipation rate for the dislocation reactions near the interface; this is controlled by microscopic mechanism (the atomic structure of the interface).This is similar to the Gurtin theory for GBs (Gurtin, 2008) which is based on a reduced dissipation inequality of the form where K is a tensor internal force defined over the region and Ġ measures the rate of defect accumulation at the interface.In order to satisfy the above inequality, the K should take the form, K = F : Ġ, where F is a fourth-order positive definite tensor.Gurtin proposed a simple dissipative constitutive relation in which F is proportional to the identity tensor, i.e, F = F I, where m > 0 is a rate-sensitivity modulus and F > 0 is a hardening-softening parameter that reflects the energy dissipation rate, similar to κ (n) in our theory.This parameter can only be determined based upon a microscopic mechanism (and atomistic simulation).The climbing-image nudged elastic band (CINEB) method (Zhu et al., 2007) may be employed to calculate the energy barrier Q for specific dislocation reactions.Then, by employing the equation κ = κ 0 exp(−Q/k B T ), where Q and k B T are the energy barrier and thermal energy, we can, in principle, determine the reaction constant κ.However, considering the numerous reactions (and reaction constants) occurring among all slip systems, this is not a practical solution.
In polycrystalline materials, plastic deformation is most frequently observed as a macroscopic average behavior across a large (compared to the atomic scale) spatial, polycrystalline region.Therefore, for simplicity, we can assume that all reaction constants may be "lumped" together, resulting in a single constant, as done in this paper.However, determining the specific reaction constant for simulations poses a question.One possible approach is to treat it as a parameter and calibrate it by fitting to experimental stress-strain data.Another approach for determining the reaction constant is to measure the activation volume Ω ≡ −(∂Q/∂σ) T .At constant temperature T and strain rate ε, the activation energy Q may be expressed as (Zhu et al., 2008): where N and ν 0 are known parameters.For a tensile test, the activation volume is Ω = k B T ∂ ln ε/∂σ, where ε and σ represent the strain rate and stress.The aforementioned approaches represent two extreme cases.On the one hand, it is impractical to calculate each reaction constant individually due to the numerous reactions occurring in polycrystalline materials and the computational cost of determining each.On the other hand, assuming identical reaction constants for all GBs oversimplifies the analysis.In order to obtain more accurate reaction constants without excessive computational resources, other improvements may be implemented.One possible approach is to utilize the nudged elastic band (NEB) method to calculate the activation energy for various types of GBs.Initially, the calculated GBs can be categorized into three distinct types: high-angle, low-angle, and coherent-twin GBs, as their activation energies differ significantly.These can be grouped, thereby reducing the number of activation energies determinations required.Subsequently, these reaction constants may be assigned to different GBs based on their respective proportions in the polycrystalline material.

Conclusion
In this paper, we proposed and demonstrated how our rigorous, conservation of Burgers vector-based, interface BC can be combined with different simulation methods to examine plastic deformation in polycrystalline and/or multiphase materials: (i) In continuum dislocation dynamics, the dislocation flux at the interface may be updated by the interface boundary condition.The numerical examples show that the reflection and transmission of dislocations at the interface can be captured by our interface boundary condition.
(ii) In the crystal plasticity finite element method, the plastic strain rate of one slip system near the interface will be coupled with the gradient of plastic strain of all activated slip systems.A bicrystal model simulation results show that the pile-ups of dislocation and stress concentration near the interface are weakened by the permeability of the interface to dislocations.Interface sliding is a natural consequence of interfacial dislocation reactions and is captured by this approach.
(iii) For discrete dislocation dynamics, possible dislocation reactions at the interface are diverse.The outgoing slip system can be determined based upon a maximum rate of entropy production criterion, which can be expressed as the product of the dislocation flux and driving force on each slip system participating in the interfacial reactions.
The successful implementation of our mesoscale interface boundary condition within the continuum dislocation dynamics and crystal plasticity finite element method has been validated through simulations.This makes it possible to apply the interface boundary condition to a wide-range of numerical plasticity approaches following similar procedures.

Figure 2 :
Figure 2: Flow chart showing the procedure for the application of the mesoscale interface BC to different plasticity simulation methods.
(k) t (x) and the geometrically necessary density b

ξ
(k) = e x × e y .Hence, the three different CDD schemes become identical and the evolution of dislocation density can be expressed as (Maxwell-Faraday equation):

Figure 4 :
Figure 4: The dislocation density distribution under an applied shear stress σxy with (a) Impenetrable interfaces (i.e., a Neumann BC) and (b) a Reaction BC (i.e., Robin BC) for a single slip system within each phase (orientations as indicated) where the dashed lines represent the interfaces.It should be noted that for Impenetrable BC case, it reached equilibrium state while there was no equilibrium or steady state for a Reaction BC case.(c) Solid and dashed lines represent the stress at the dislocation source τ α and the dislocation density at the interface for the two types of BCs (Impenetrable and Reaction); (d)-(f) The dislocation distribution under uniaxial loading for the case of a Reaction BC for the case of two slip systems in α and one in β at times, t = 0.4, 0.6, 0.8.

Figure 5 :
Figure 5: A three grain lamellar microstructure with two slip systems and multiple sources within each grain, subject to an external shear stress σxy.(a) The dislocation density distribution.(b) The corresponding von Mises plastic strain ϵ p .

Figure 6 :
Figure 6: (a) Schematic of finite element model; (b) The contours of Von Mises stress, (c) dislocation density for one of slip systems and (d) displacement are plotted, the left and right-hand sides correspond to the Robin and Neumann BC respectively; (e) The inset shows the sliding between two grains near the interface under Robin condition; (f) Comparison of engineering stress-strain curve under Robin and Neumann BC condition.

Figure 7 :
Figure 7: Schematic of transmissions of dislocation in DDD.(a) Single dislocation line from i th slip plane is approaching interface; (b) the dislocation line rotate by β angle; (c) dislocation nucleation of k th slip system.

Table 1 :
Parameters employed in the 2D CDD simulations.

Table 2 :
CPFEM parameters for copper