A method of coupling discrete dislocation plasticity to the crystal plasticity finite element method

A method of concurrent coupling of planar discrete dislocation plasticity (DDP) and a crystal plasticity finite element (CPFE) method was devised for simulating plastic deformation in large polycrystals with discrete dislocation resolution in a single grain or cluster of grains for computational efficiency; computation time using the coupling method can be reduced by an order of magnitude compared to DDP. The method is based on an iterative scheme initiated by a sub-model calculation, which ensures displacement and traction compatibility at all nodes at the interface between the DDP and CPFE domains. The proposed coupling approach is demonstrated using two plane strain problems: (i) uniaxial tension of a bi-crystal film and (ii) indentation of a thin film on a substrate. The latter was also used to demonstrate that the rigid substrate assumption used in earlier DDP studies is inadequate for indentation depths that are large compared to the film thickness, i.e. the effect of the plastic substrate modelled using CPFE becomes important. The coupling method can be used to study a wider range of indentation depths than previously possible using DDP alone, without sacrificing the indentation size effect regime captured by DDP. The method is general and can be applied to any problem where finer resolution of dislocation mediated plasticity is required to study the mechanical response of polycrystalline materials, e.g. to capture size effects locally within a larger elastic/plastic boundary value problem.


Introduction
In order to model most well-known material size effects, discrete resolution of defects in the crystal lattice is required. For example, grain size (Hall-Petch) and specimen size (micro-pillar and micro-film) strengthening cannot be captured by conventional crystal plasticity theories [1], however these effects arise naturally from discrete dislocation dynamics methods that capture the interactions between dislocations and interfaces [2]. However, large-scale boundary value problems remain computationally intractable for dislocation dynamics with respect to spatial and/or temporal scales. Embellishment of continuum methods with physically-based equations capturing phenomena occurring at smaller scales, as in strain gradient plasticity theories [3] and geometrically-necessary dislocation based crystal plasticity formulations [4] is one way of approaching the problem, however the mechanisms causing size effects are to some extent prescribed when taking that approach. Multi-scaling aims to alleviate the computational burden associated with small scale modelling methods by using a higher resolution method only in an area of interest, where it is required.
Multi-scale modelling generally employs either unidirectional (sub-modelling) or bi-directional (concurrent modelling) information passing. The former benefits from simplicity of implementation and computational efficiency, whereas the latter is more accurate, but at a cost. There is now a wide range of methods that have been established over recent decades, for example coupling electronic structure to atomistic modelling [5], atomistic to discrete dislocation dynamics [6], atomistic to continuum [7], etc. A natural method for capturing localised size effect behaviour in large polycrystal boundary value problems is concurrent coupling of discrete dislocation plasticity (DDP) to the crystal plasticity finite element (CPFE) method, a continuum method that accounts for the crystal structure and homogenises the effect of crystallographic slip.
A planar method for concurrent coupling of DDP to CPFE in a single crystal was developed in [8]. The coupling boundary between the DDP and CPFE regions was treated as transparent to dislocations, and the compatibility of tractions and displacement was ensured in an average sense. Tension and mode I crack growth problems were analysed; in the latter, computation time was reduced by up to a factor of 14 with only a minor reduction in accuracy. 3D DDP was coupled to strain gradient crystal plasticity via the solution of separate models with bi-directional information passing in [9]. The so-called multi-scale dislocation dynamics plasticity (MDDP) approach was also used in [10] to investigate micropillar compression, and in [11] to study spherical nano-indentation. Micro-pillar compression was also investigated in [12] using a hybrid DDP/finite element method (FEM) coupling method.
In this paper we present a methodology for concurrent coupling of planar DDP to crystal plasticity; the method is useful for, for example, extending a well-established relationship between hardness and indentation depth for a thin film on a substrate to include the indentation size effect regime, for a wide range of coating and substrate materials. In the proposed method, compatibility of displacements and tractions is ensured at the coupling boundary through an iterative boundary condition passing scheme; this is achieved at each node in the coupling boundary rather than in an average sense as done in prior studies [8].
Small-scale indentation is being increasingly used to measure yield stress and Young's modulus for films and small material volumes. It has been established experimentally that continuum relationships for hardness applicable to large scale conical (or pyramidal) and spherical indentation [13,14] fail when indentation depths are sufficiently small (nano/micro); the observation is that sufficiently smaller indentation depths correspond to larger indentation pressures [15][16][17][18][19]. The effect of the substrate material on the response becomes important at larger indentation depths, which can be characterized by empirical relationships [20,21].
At present no empirical relationship exists that is effective over the full range of indentation depths encountered in typical experiments.
Among the early models of the indentation size effect, 'Nix and Gao' developed an analytical formula that predicted size-dependent hardness based on geometrically necessary dislocations and the strain gradient [22]. Some further work based on strain gradients was conducted in [23][24][25][26], in order to study the relationship between the indentation size effect and the material length scale parameters. In a different approach, Guo et al [27] developed the Taylorbased nonlocal theory of plasticity based on general non-local plasticity theory [28].
However, continuum methods require fitting a length scale parameter to experimental data, e.g. the formula developed by Nix and Gao works well but only captures the experimental behaviour [29] after suitably fitting a length scale parameter in the model to the data. On the other hand, the length scales in DDP, which captures dislocation mediated size effects naturally, are measurable and have clear origins in the microstructure. The indentation size effect has been studied in detail using DDP methods. The work reported in [30][31][32] primarily focused on size effects in nano-indentation by sinusoidal and wedge-shaped indenters on crystals characterised by a low density of sources. The work was devoted to the investigation of the effect of microstructure, including nucleation source density, and number and orientation of pre-defined slip planes, on the indentation pressure [33][34][35]. Furthermore, the effect of indenter shape on the indentation pressure was studied using numerical simulations [36,37], which have also been successfully corroborated by experimental data [17]. The damage due to dislocation activity in indentation processes was simulated in [38], which aimed to provide the link between the dislocation structure and the macroscopic fatigue phenomenon. Recently, Nicola and co-workers [39][40][41][42] conducted DDP simulations of the flattening of surface roughness obtained by nano-imprinting thin single crystal films. In all the aforementioned studies, the dislocation activity is studied in simplified configurations, often with the film is adhering to a rigid substrate, and these studies cannot be easily extended to scenarios in which the crystal films are adhering to or resting on polycrystalline substrates characterised by similar elastic/plastic response.
Another important area in which coupled DDP/CPFE simulations may shed light on important mechanisms is the study of coating size effects across different length scales. In particular, Korsunsky et al [20] developed a formula that provides the hardness variation of a coated system as a function of indentation depth; it accounts for the effect of the substrate on the hardness of a thin film during the indentation process by looking at the total deforming volume and the work of indentation [21]. This method captures the transition between coating-dominated and substrate-dominated response well for many film-substrate systems, such as Cr-Al and Ni-Cu; however, it does not account for the indentation size effect observed at small indentation depths. The DDP/CPFE coupled framework developed here can, therefore, be used to explore the behaviour of coated systems across different length scales, for example by adding information about the indentation size effect regime at the smallest scales explored in other studies (e.g. Balint et al [34]) to the formula proposed in [20]. This is the subject of a forthcoming investigation.

Crystal plasticity finite element (CPFE) formulation
The CPFE implementation adopted in this project is the rate-independent formulation described in [43], which is based on the prior work of Asaro [44] and Anand and Kothari [45]. Elastic anisotropy and different crystal structures have since been incorporated to study energy dissipation at the meso-scale and localized plasticity and residual stresses (e.g. [46][47][48]). The formulation is based on the standard two-term multiplicative decomposition of the deformation gradient into the elastic and plastic parts, where plastic slip occurs on individual slip systems when the resolved shear stress on the slip system exceeds the current critical resolved shear stress, τ c . The rate of change of the resolved shear stress on the αth slip system is taken to evolve as [43]: where h αβ is the hardening modulus matrix, γ β is the shearing rate γ on the th β active slip system, and P α is given by: where s α and n α are the slip direction and slip plane normal of the th α active slip system, respectively, C is the fourth order tensor of elastic constants and D is the second order deformation rate tensor calculated from the velocity gradient L: and σ is the Cauchy stress and β α is computed as: with s n n s 1 2 .
T T The matrix of hardening moduli is defined as: with q representing the latent hardening parameter, and the self-hardening function, h, is defined as: where H F , H R , and H exp are constitutive variables used to capture the initial and final hardening modulus and p is the accumulated plastic strain. An example of how these parameters can be selected to match the experimental stress-strain response is reported in [48]. For the purpose of this article, the material parameters (critical resolved shear stress and hardening) are obtained by matching the plastic behaviour to the DDP model in uniaxial tension (see section 2.3.1), and are reported in table 1.
Rearranging equation (1) into a set of general linear equations, the unknown γ α can be found using a singular-value decomposition to avoid potential issues related to matrix inversion [49][50][51]. This enables the model to update stresses and strains throughout the simulation and capture the evolution of the accumulated plastic deformations along the slip systems.
The CPFE implementation discussed above is versatile and allows simulating both 2D and 3D real polycrystalline structures obtained using e.g. EBSD or randomly distributed grains or statistically equivalent grain structures modelled using Voronoi tessellations [46,52].

Discrete dislocation plasticity formulation
The plane strain, isotropic, quasi-static DDP formulation of Van der Giessen and Needleman [53] is used here. A fcc crystal structure is assumed, with the plane of the simulation taken perpend icular to the 101 [ ] crystal direction to satisfy the plane strain constraint, 0 33 13 23 Slip planes are lines in the planar formulation, defined by the intersection of the slip planes with the plane of the simulation, and dislocations are confined to glide along those lines. The material is assumed to be initially dislocation-free; all dislocations originate from Frank-Read sources. The plane of the simulation is that which intersects the pure edge segments of the incipient loop originating from a Frank-Read source, and the screw parts of the loop are effectively pinned by the plane strain condition and the planar loading, as depicted in figure 1. In 2D, the trapped dislocation line of the Frank-Read source is perpendicular to the simulation plane, hence is a point, and the dislocation loop originating from the source is a dipole of straight line edge dislocations with Burgers vector b. The applied loading is such that the glide component of the Peach-Koehler force is large on the edge segments of the dislocation loop, and small on the screw segments, hence the latter are effectively pinned and the loop expands via the edge segments.
Boundary conditions are satisfied using the superposition principle first employed by Lubarda et al [54], and Van der Giessen and Needleman [53] (see figure 2 and equation (8)). Displacement, stress and strain are decomposed into the infinite medium dislocation fields (˜) and a correction (^) that ensures the boundary conditions are satisfied; the former is obtained via the superposition of known analytical fields and the latter is obtained via a finite element solution of a problem where singularities are absent and the effect of the dislocations is mediated by the corrected boundary conditions.
Point Frank-Read sources are randomly distributed on the slip planes (lines) to a certain density. Source strengths nuc τ are taken from a normal distribution with specified mean and standard deviation; the latter captures the effect of the statistical distribution of trapped dislocation line lengths in the material. All sources nucleate dipoles with Burgers vector b at the equilibrium spacing, L nuc , given by: where μ is the shear modulus and ν is Poisson's ratio, which corresponds to an exact balance between the resolved shear stress and the dislocation line tension. This occurs when the nucleation stress nuc τ is reached and maintained for the nucleation time t nuc ; the latter is the time taken for the loop to reach its unstable configuration (Benzerga [55]). This in general depends on the applied stress, but is relatively insensitive to stress at low stresses, hence is taken to be a constant 10 ns in these simulations.
The Ith dislocation glides with velocity V I ( ) according to a linear mobility law:  where B is the drag coefficient and f I ( ) is the Peach-Koehler force on the I-th dislocation, given by: where n I ( ) is the slip plane unit normal, and dislocation self stress is omitted. A cut-off glide velocity v 20 m s cutoff 1 = − was found to prevent exceedingly small time steps without affecting the results; this is the same value used elsewhere, e.g. [56,57]. The introduction of a cut-off velocity aims to improve the simulation efficiency without inducing side effects, as described in [58]. Annihilation of a pair of dislocations of opposite sign on the same slip plane occurs when the distance between them is less than the critical annihilation distance L e = 6b. Obstacles are modelled as points that pin dislocations that attempt to glide over them. An obstacle releases a pinned dislocation when the Peach-Koehler force on the obstacle exceeds b obs τ , where obs τ is the obstacle strength. In the example problems studied below, the DDP region is modelled as one single crystal. The material properties and dislocation structures used in the simulations are reported in table 2. The material is initially stress-and dislocation-free with slip systems set as 54.7, 0, 1, 2, 3 with respect to the horizontal axis (this configuration corresponds to the plane strain slip system orientations of fcc crystals). Dislocation sources are randomly distributed on slip planes that are spaced b 100 apart, with a density nuc ρ in the specimen. The obstacles are randomly distributed with a density . A time step of t 0.5 ns ∆ = is used to resolve the motion and interactions of dislocations, hence a high loading rate is applied, e.g. U L / 600 s 1 = − is prescribed on the right-hand boundary of the DDP region in the tension case in order to obtain a strain of 0.0006 ε = in 2000 time steps.

General coupling framework
The method is based on a concurrent coupling of a sub-model of DDP inside a CPFE region (see figure 3). The fine scale area of interest is the DDP region (Ω DD ), while the remainder is modelled using a CPFE description (Ω CPFE ). The aim of the proposed coupling strategy is to achieve compatibility of stresses and displacements at a local level (i.e. at each node) instead of exchanging average information along their mutual boundaries.

Calibration of the material parameters.
The first requirement before coupling is to match the material behaviour. It is necessary to match both the elastic and plastic response of the two regions by modulating the CPFE inputs to match the stress-strain curve obtained using DDP in uniaxial tension; the latter is taken to be the actual material response, since it is based upon fundamental quantities rather than phenomenological input. Although the CPFE model is rate-independent and the DDP is rate-dependent, the behaviour of the DDP for a given rate (that which is applied) is effectively calibrated into the CPFE material response. Figure 4 shows the comparison between the two stress-strain responses after the CPFE calibration was performed. The parameters used in the dislocation simulations are summarized in table 2. They are also used in the indentation problem. The CPFE properties obtained from the uniaxial tension calibration are also applied to the CPFE region in the indentation problem, assuming the average strain rate under indentation, for the rate applied, is comparable to that experienced under uniaxial. The best match between the two descriptions was achieved using the CPFE parameters tabulated in table 1.

Interface coupling.
The key feature of the proposed coupling strategy is the passing of boundary conditions between two sub-regions via the coupling interface. The two subregions, Ω DD and Ω CPFE in figure 3, communicate via the coupling boundaries, Γ DD and Γ CPFE ; in particular, displacements and stresses are passed at nodal points, which are shared by the two finite element descriptions used for the correction field ( ). in the DDP simulation and  the CPFE description, in a forward-iterative manner. At every loading step, which is defined by updating the external applied loads and/or displacements, T and u respectively, along the external boundaries of the two domains, Ω DD and Ω CPFE , one field variable (either stress or displacement) is passed between the two coupling boundaries Γ DD and Γ CPFE until the other field variable converges between two consecutive iterations to within an error. The coupling interface in this method, which is intended for polycrystals, always coincides with a grain boundary and/or interface between plastically dissimilar materials, which is assumed in this study to be impenetrable to dislocations. Indeed, this feature will introduce convergence issues when a significant number of dislocations pile up at the coupling interface. Hence, in order for the stresses and displacements passed between the models to converge, it is sometimes beneficial to use a buffer zone in the DD region near the coupling interface. This buffer zone stops dislocations just before they reach the boundary, which effectively samples the dislocation fields outside of their singular cores, hence improving convergence; this is similar to other averaging and interpolation schemes proposed in the literature [6]. Furthermore, a running-average algorithm has been adopted to smooth the traction field passed to the CPFE model in the indentation case where significant numbers of dislocations pile up at the coupling interface. This helps to achieve convergence in the indentation case. It should be noted that the buffer zone and smoothing of the traction field is not always required, e.g. in the uniaxial tension case. They are used when needed as a compromise between efficiency of the convergence and retention of sufficient discrete detail.
A coupling strategy at the boundary that separates sub-regions with different material length scales is required not only to smooth the transition between the different descriptions but also to impose rules for the transfer of dislocations across the coupling boundary. Curtin et al developed a technique to deal with the passage of dislocations between the two regions Macroscopic stress-strain responses predicted by DDP and CPFE for the Al crystal after calibration of the CPFE constitutive parameters to produce the best match between the two simulations. [8]. Although their technique could be employed here, we align the interface between the DD and CPFE regions with an interface or grain boundary, i.e. dislocations are pinned at the interface and do not pass through it.
Advantages of the proposed coupling strategy include: (i) the flexibility to use any material description via the use of a commercial code (Abaqus) with user defined material subroutines; (ii) node-wise coupling of tractions and displacements rather than in an average sense; (iii) bi-directional coupling, rather than a one-directional passing of boundary conditions from the primary model to the sub-model.

Implementation of the coupling strategy and demonstration
The coupling strategy has been implemented to study two problems, namely uniaxial tension and indentation. In the following subsections we discuss the details of the algorithm implementation for these two problems.

Uniaxial tension: problem definition.
A 2D uniaxial tension problem, schematically shown in figure 5, is used as a vehicle to illustrate the implementation of the coupling procedure. The specimen is sub-divided into two rectangular L H 6 m 2 m μ μ × = × regions, modelled by CPFE and DDP, respectively. The CPFE sub-model is made of a single fcc crystal, represented by the fully 3D fcc slip systems although subject to plane strain, whose axes are oriented at zero degrees with respect to the reference system, i.e. parallel to the x-axis. In the DDP sub-model, also representative of a fcc crystal structure, the slip systems are ( ) φ α = ±°= α 54.7 , 0, 1, 2, 3. The dislocation structure properties are summarized in table 1. A uniform displacement in the x-direction is applied along the right (CPFE) edge of the sample; the other end of the specimen (left edge of the DDP domain) is constrained to not move in the x-direction, and one node is also prevented from moving in the y-direction to prevent translation of the sample. The demonstration consists of 9 coupling steps (loading increments), achieving 0.054 m μ final displacement (0.45% strain) on the right boundary.

Coupling strategy for the uniaxial tension problem.
After the initial calibration of the CPFE model, which is carried out to match the elasto-plastic response of the discrete dislocation model, the DD model is replaced by an equivalent elasto-plastic continuum description to initialise the simulation (the start-up model; see figure 6). This allows displacing the right hand boundary by U x until the onset of plasticity is exceeded, while ensuring continuity at boundary A during loading of the specimen; it generates an initial stress profile from the CPFE region to start the coupling, and amounts to an initial guess. The DD model is also separately initialised and set to receive the traction distribution at the interface DD Γ as the boundary condition on its right face. Once the start-up model exceeds yield, the tractions are passed to the DD model. The response of the DD model to the applied tractions is then computed. The displacements at the boundary between the two domains (the right face of the DD region) are then passed to the CPFE model and applied at CPFE Γ . The resulting tractions at the interface between the two models are again passed to the DD model, and the cycle is repeated until convergence is reached for the level of imposed strain. The above procedure is then repeated at predefined strain targets until the maximum value of imposed displacement, U x,max , is reached.

Indentation: problem definition.
A 2D plane strain analysis of indentation of a polycrystalline material is conducted using the coupling strategy outlined above. The key region of the polycrystalline material, i.e. where dislocation activity is likely to be high, is located underneath the tip of the indenter and is modelled using DDP. The surrounding region is model led using the CPFE description, which is modelled here as a single fcc crystal whose axes are oriented at zero degrees with respect to the reference system, i.e. parallel to the x-axis. The multi-scale modelling of the indentation problem is schematically represented in figure 7.
The indenter is modelled as a wedge-shaped rigid body with indentation angle 5 ω =° with respect to the x-axis. This wedge angle is constrained to a smaller value than that used in most experimental studies (e.g. Berkovich or Vickers) because larger indenter angles would lead to large contact area jumps in the discrete dislocation calculations [34].
The dimensions of the entire sample are . The size of the surrounding CPFE region is varied accordingly. The finite element mesh in the DDP region is highly refined with smallest element size μ μ × 0.04 m 0.5 m in order to model the contact accurately and ensure the effect of the singular (˜) fields is fully captured on the boundary of the ( ) problem. The finite element mesh in the CPFE region is characterised by increasingly smaller elements towards the DDP region and a much coarser description away from the DDP region. Another feature of the CPFE mesh is that the node positions and element sizes are exactly the same as the DDP region at the coupling boundary.
Perfect sticking (no slip) is assumed at the boundary between the DDP surface and the indenter. Thus, the boundary conditions are: where S contact denotes the part of the upper surface in contact with the indenter and δ is the indentation depth. The contact between the rigid wedge and the film is calculated based on the deformed top surface; the contact area 'A' (length) is defined as the difference between the x-coordinates of the rightmost and leftmost nodes in the contact. This gives results that are generally different from those based on the nominal contact area, ; this is due to material sink-in or pile-up of the deformed surface. The effect of surface roughness is neglected [59]. The other boundary conditions are: Here, T n i i j j σ = is the surface traction on a surface with outward normal vector n j . The total reaction force on the film in response to the indenter is computed as: The indentation pressure p is defined as the ratio of the indentation force to the actual contact length. A relatively high value of indentation displacement rate ˙0 .4 m s 1 δ = − was used in order to reduce the computational cost without compromising the accuracy of the results. The calculations consist of 30 coupling steps (loading increments); all the simulations reported below have a maximum applied indentation depth of 300 nm .

Coupling strategy for the indentation problem.
The coupling strategy for the indentation problem is illustrated in figure 8. The whole process initiates with a DDP calculation  alone on the region underneath the tip of the indenter, which is displaced vertically at the prescribed displacement rate. All degrees of freedom are set to zero along the coupling boundary. This initialization produces a preliminary stress field distribution on the coupling boundary. The stress obtained from the DDP initial simulation step is passed to the CPFE region at the nodes distributed along the coupling boundary. CPFE calculations are then performed and the displacements at the coupling interface are re-evaluated and passed back to the DDP region as new boundary conditions for the coupling boundary. The coupling loop ends when the difference between the tractions at the coupling boundary is lower than a given threshold. The coupling process will advance to the next step by increasing the indenter penetration, and will terminate when the indenter reaches the prescribed maximum indentation depth. In this case it was advantageous to apply a buffer zone since a large number of dislocations piled up at the coupling boundary. The height of the buffer zone was set to h b 200 0.05 m buffer μ = = , and its width was identical to the width of the DDP sub-region l 30 m buffer μ = .

The uniaxial tension case
The uniaxial tension case was studied first to test the methodology of the proposed concurrent coupling strategy. A comparison was performed between the response obtained using two coupled solutions: (i) the first approach employs a simple non-iterative and unidirectional coupling, whereby the tractions are passed from the CPFE solver to the DDP region at different times during the monotonic loading of the specimens while the CPFE boundaries are not affected by the evolution of the DDP domain; (ii) the second solution is obtained making use of the proposed fully coupled iterative scheme. Stress-strain curves are plotted to illustrate the differences between the two schemes. Figure 9 depicts the macroscopic stress-strain curves for the tensile specimen shown in figure 5, as predicted using the different coupling routes. The response obtained using the DDP and CPFE solutions after the initial calibration is performed are also reported on the plot and serve as references. The dashed blue curve shows the macroscopic stress-strain behaviour calculated in the DDP domain using the unidirectional coupling scheme. Although the traction distributions are correctly matched at the DDP/CPFE boundary, compatibility of displacements is not ensured using the sub-model approach and there is no feedback from the DDP region to the CPFE; the effect is small in this simple tension calculation, which is shown only to demonstrate the algorithm, but can be significant in other problems (see section 3.2 on the indentation case). The dashed black line shows the solution obtained using the iterative coupling strategy proposed by the authors. The solution is deemed convergent when the average values of the traction components at the interface between the two domains do not vary by more than 2%. The fully coupled formulation imposes displacement and traction matching at the interface between the DDP and the CPFE domains, hence produces a more accurate representation of the stress-strain response of the entire specimen than the unidirectional submodel approach; again, the difference between the methods is considerably more pronounced in more complicated problems, which will be demonstrated later.
A comparison of the DDP response for the different solution methods at the coupling interface, in terms of the evolution of dislocation density and the macroscopic stress-strain behaviour is reported in figure 10. It can be observed that the implementation of the coupling has an appreciable effect on the dislocation response, even in this simple tension problem, which also varies for the two different coupling schemes. In particular, both a greater decrease in the density of dislocations and relaxation of the stress are observed in the in DDP region, compared to the DD-only calculation.
A contour plot of the σ 11 component of stress for the entire tensile specimen, is shown in the inset of figure 10 for a strain 0.25% ε = . This shows both the evolution of dislocation structures in the DDP region and the associated plastic activity and the transition between the DDP and the CPFE domains, with the continuity of stresses across the boundary portrayed.

The indentation case
The results in figure 11 for indentation pressure versus depth show the effect of a substrate modelled using the proposed coupling method. For the coupling calculations with film thickness h 10 m μ = modelled using DDP, on a crystal plasticity substrate, the deviations in the indentation pressure are small compared to the curves obtained using a DDP model alone for the film (on a rigid substrate). The rigid substrate assumption in the DDP-only model isn't necessary, but is a simplifying assumption that has been made in prior studies (e.g. [34]) so is used here for comparison purposes; including a DDP model of the substrate is possible but would lead to considerably higher computational cost for the dimensions considered here compared to the coupling method. For a film thickness h 2 m μ = , the difference between the two simulations is much more significant and the presence of the CPFE substrate leads to a lower indentation pressure (softer) with increasing indentation depth. The indentation pressure curve for a h = 2 μm film on an elastic substrate is also included for reference in figure 11, which lies in between the DD-only (with a rigid substrate) result and that of the coupled calculation with CPFE substrate. A least-squares error estimate is used to quantify the difference between the same films modelled using DDP only (with a rigid substrate) and the coupled DDP/CPFE approach. The results are shown in table 3 for three DDP region film thicknesses.
The results obtained using the coupling strategy proposed here begin to differ considerably from those obtained using DDP alone (on a rigid substrate) when the DDP film thickness decreases below 5 μm; the solution obtained for thicker films is further validation of the coupling methodology in capturing the known DDP limit. Hence it is capable of capturing the effect of a non-rigid substrate (in this case characterised by a CPFE description matching the elasto-plastic behaviour of the film) when performing indentation, without the associated computational cost of modelling the entire system using DDP. This effectively expands the potential applications of DDP in modelling indentation and other problems to potentially much larger systems where the details of dislocations activity are only relevant in a region of interest; therefore the coupling method can predict the hardness of thin films and multi-layers on substrates, with elastic/plastic material properties that differ from those of the film, more accurately than was possible in [34]. Figure 12 shows the normal stress component in the coupled simulation with a 0.2 μm indentation depth. The dislocation activity within the domain is clearly visible, as is the continuity in the stress across the DDP/CPFE boundary and the pinning of the dislocations at the boundary between the two domains. It is evident that the coupling method improves the accuracy of the calculations compared to the stress and dislocation structure obtained in the simulations where the substrate is assumed to be a rigid boundary or a purely elastic substrate. The normalized dislocation density evolution versus the discrete dislocation film thickness h is shown in table 4. Both total and mobile dislocations are unaffected by the coupling approach when the DDP domain size is sufficiently large. However, the coupling effect, which manifests itself as the influence of the CPFE substrate in these simulations, becomes more apparent on the dislocation density when the DDP domain is reduced in size.
The deformed surface profiles for the h 2 m μ = film simulated using the DDP-only (with a rigid substrate) and coupled methods are shown in figure 13. Results are shown when the applied indentation depth 0.3 m δ µ = . The pile-up of material at the contact surface obtained from the coupled calculation is less significant than that from the DDP-only method (with a rigid substrate), further indicating the effect of the CPFE substrate in acting to relieve stress and thus reduce the dislocation activity in the film.
The multi-scale coupling methodology proposed here is strongly motivated by the need to improve the simulation efficiency by avoiding the use of DDP in the whole domain of the problem. Results obtained using the coupling method show that computation time can be reduced by a factor of 5.6 to 9 for the indentation simulations depending on the film thickness, which dictates the height of the DDP region, compared to corresponding simulations obtained by modelling the entire 300 μm × 100 μm domain with DDP.

Conclusion
A concurrent coupling framework between CPFE and DDP models linked via a scripted interface has been proposed, and a uniaxial bicrystal tension problem was used to demonstrate and verify the coupling strategy. The coupling used a CPFE material model implemented in  a general-purpose commercial finite element package (Abaqus) to showcase its wide applicability; any material description could be used for the primary model via the standard features of Abaqus or user defined material subroutines. A plane strain analysis of the indentation of a single crystal film bonded to a polycrystalline substrate by a rigid wedge-shaped indenter was also done. Although the computations were performed within a small strain framework, contact was modelled by considering the deformed surface profile of the films. It was demonstrated that the proposed coupling method efficiently captures the effect of a sub-domain with discrete dislocation resolution within a crystal plasticity model for incorporating size effects in a boundary value problem with much greater computational efficiency than DDP alone. Furthermore, prior DDP analyses of indentation [34,[39][40][41][42] relied on a rigid substrate condition to simplify the computation. It was demonstrated that for an indentation depth of 0.3 m μ , the response of films of thickness 2 m μ or less exhibited a significant dependence on the substrate: a crystal plasticity representation of the substrate led to appreciably smaller indentation pressures and diminished stress and dislocation activity in the film compared to the rigid and elastic substrate representations typically used in DDP analyses of indentation.
Using the proposed method, it is possible to embed a DDP region in any domain that can be modelled using Abaqus. The inclusion of schemes for capturing and correctly modelling the dislocation flux at the DDP/CPFE interface (following e.g. [58]) will be the subject of future investigations. Further efforts will also focus on the use of the proposed technique to develop a physically-based relationship to explain the transition in hardness of thin films on substrates across the scales, i.e. dislocation dominated size effects at smaller indentation  depths transitioning to the layer-dominated and substrate-dominated hardness observed at larger indentation depths, as described in previous studies [20]; nano-sliding and fretting problems will also be modelled using this method in a forthcoming study.