A computational approach towards modelling dislocation transmission across phase boundaries

ABSTRACT

To study the nanoscopic interaction between edge dislocations and a phase boundary within a two-phase microstructure the effect of the phase contrast on the internal stress field due to the dislocations needs to be taken into account. For this purpose a 2D semi-discrete model is proposed in this paper. It consists of two distinct phases, each with its specific material properties, separated by a fully coherent and non-damaging phase boundary. Each phase is modelled as a continuum enriched with a Peierls–Nabarro (PN) dislocation region, confining dislocation motion to a discrete plane, the glide plane. In this paper, a single glide plane perpendicular to and continuous across the phase boundary is considered. Along the glide plane bulk induced shear tractions are balanced by glide plane shear tractions based on the classical PN model. The model's ability to capture dislocation obstruction at phase boundaries, dislocation pile-ups and dislocation transmission is studied. Results show that the phase contrast in material properties (e.g. elastic stiffness, glide plane properties) alone creates a barrier to the motion of dislocations from a soft to a hard phase. The proposed model accounts for the interplay between dislocations, external boundaries and phase boundary and thus represents a suitable tool for studying edge dislocation–phase boundary interaction in two-phase microstructures.

KEYWORDS:

• Dislocations
• dislocation pile-ups
• Peierls–Nabarro model
• phase boundary
• plasticity

1. Introduction

Grain and phase boundaries have been shown to strongly influence the ductile behaviour of polycrystalline microstructures, as described e.g. by the Hall–Petch relation [1]. Internal boundaries constitute barriers against dislocation motion, resulting in dislocation obstruction and dislocation pile-ups. Eventually, under sufficiently high shear stress, a variety of events may occur: transmission of a dislocation into the neighbouring grain, dislocation reflection, dislocation nucleation in the neighbouring grain or dislocation migration into the boundary. The aim of this work is to gain a better understanding of the interaction between edge dislocations and phase boundaries related to dislocation obstruction, pile-up and transmission. This involves the non-local interaction between the individual dislocations and also between dislocations and phase boundary. While analytical solutions are limited to simple cases, most non-trivial studies require a numerical treatment. To address the latter, various approaches have been pursued in the literature. Employing Discrete Dislocation Dynamics (DDD) [2,3] crack growth within single crystals [4] and along bi-material interfaces [5,6] has been studied. This method is capable of modelling discrete dislocation activity in domain sizes $>10\mu \mathrm{m}$ and represents an efficient tool to predict corresponding macro-scale material behaviour. However, DDD requires case-specific constitutive relations for dislocation–interface interactions (e.g. dislocation transmission). Furthermore, DDD is based on the analytical solution of single dislocations – often according to the Volterra model – which is questionable in the vicinity of phase and grain boundaries where the change in dislocation core can play a dominant role.

Molecular dynamics (MD), on the contrary, accounts for the full atomistic complexity of dislocation behaviour. It is thus capable of fully capturing the local interaction between dislocations and the interface. Yet, the utilisation of a proper interatomic potential is essential. This can be a problem within the dislocation core, as pointed out by Das and Gavini [7]. Recently, MD was utilised to study dislocation pile-ups, by Wang [8] in Al with $\mathrm{\Sigma }11$ asymmetrical tilt grain boundaries and by Elzas and Thijsse [9] in iron-precipitate microstructures with coherent, semi-coherent and non-coherent interfaces. The complexity of the simulations, however, entails a high computational cost. Hence, many simulations are limited to a small number of dislocations. To model larger systems, the quasicontinuum method (QCM) [10] has been developed. It preserves atomistic accuracy where required and approximates with continuum modelling elsewhere. A particular version of QCM is the coupled atomistic discrete dislocation (CADD) method. It was proposed by Shilkrot et al. [11] and treats the continuum domain with a DDD framework. Both methods were recently adopted to study dislocation–grain boundary interactions in [12] (QCM) and [13–15] (CADD). Besides the complex implementation of the pad-region connecting the continuum and the atomistics, the proper choice of the interatomic potential remains a crucial matter.

A full continuum approach that retains the discrete nature of dislocations was introduced by Acharya in [16] based on the early work of Kröner [17,18], Mura [19], Fox [20] and Willis [21]. In this theory, dislocations are viewed as sources of elastic incompatibility and are modelled with a dissipative dislocation transport equation. Following De Wit [22], Fressengeas [23] included disclinations to account for rotational defects such as grain boundaries. Recently, Zhang et al. [24,25] incorporated the effect of lattice periodicity and stacking fault energy to capture dislocation core effects and studied arrays of edge dislocations on one or multiple slip planes in homogeneous and bi-metallic media. The phase boundary is considered as impenetrable and hence introduces an artificial influence on dislocations interacting with the phase boundary.

Another approach was established by Shen and Wang [26] with a phase-field methodology. The total energy is decomposed in the stacking fault energy of the glide plane, the elastic strain energy following the microelasticity theory of Khachaturyan and Shatalov [27] and a gradient energy term for the dislocation core. The Ginzburg–Landau relation describes the evolution of the dislocation structure. Mianroodi et al. [28,29] extended this framework to the atomistic determination of the gradient term via energy scaling and showed an improved accuracy of the dislocation core structure in comparison with Peierls–Nabarro modelling. In a recent study, Zeng et al. [30] studied transmission of a single dislocation across a bi-metallic interface. Considering a lattice mismatch between both phases, coherency stresses and residual dislocations at the phase boundary were incorporated. However, basing the coherency stress purely on the mismatch of material properties, no change in coherency stress due to dislocation interaction is included – as one would expect from the change of the boundary structure as a result of approaching dislocations. Also, no spreading out of a screw dislocation core into the interface is possible – as observed by Anderson and Li [31] and by Shehadeh et al. [32]. One of the major drawbacks of this method is computational cost. Hence, domain sizes in the order of 100 nm are the largest feasible, resulting in relatively strong image forces arising from the commonly applied periodic boundary conditions.

To investigate the interaction of a single screw dislocation with a slipping and coherent bi-metallic interface, Anderson and Li proposed in [31] a simple and computationally efficient approach based on the Peierls–Nabarro (PN) model [1,33,34]. Following the analytical approach of Paracheco and Mura [35], they described the screw dislocation by a series of infinitesimally small dislocations with the two-phase stress field by Head [36] and balanced it against the glide plane tractions of the Frenkel sinusoidal [37]. Shehadeh et al. [32] extended this model for studying screw dislocation transmission across a non-coherent Cu–Ni interface. The Generalized Stacking Fault (GSF) energy surface [38,39] is utilised for accurate glide plane and phase boundary tractions and additional coherency stresses describe the lattice mismatch between both phases. Similar to the phase-field approach by Zeng et al. [30], a change of the coherency stresses due to the local change in phase boundary structure is not accounted for. Additionally, due to its simplicity, an extension towards higher complexities such as phase boundary decohesion or multiple slip systems is not straightforward.

Other approaches based on the PN model concentrated only on studying dislocations in homogeneous media and not on dislocations in two-phase microstructures: Schoeck [40,41] studied the dissociation of dislocations and the energy of a dislocation dipole considering the GSF energy surface and using a variational approach with an ansatz for the dislocation profile. Rice [42] suggested a methodology to analyse the nucleation of perfect and partial dislocations from stressed crack tips based on the J-integral. Xu et al. [43,44] extended the variational boundary integral method of Xu and Ortiz [45] towards the PN dislocation. By modelling the dislocation as a continuous distribution of infinitesimal dislocations [46] they analysed dislocation nucleation from crack tips under mixed mode loading. Later, Xu and Argon [47] adapted this method towards nucleation of dislocation loops in perfect crystals via a centred circular perturbation and Xu and Zhang [48] towards dislocation nucleation from crystal surfaces. Xiang et al. [49] introduced with the generalised PN model an extension towards 3D, including not only a planar disregistry but also a normal opening of the glide plane. The elastic stress field of a curved dislocation was assessed with Mura's formula. They then studied the core structure of perfect and partial dislocation loops in Cu and Al under different values of applied stress.

In this paper, the PN model is utilised in 2D plane strain modelling of edge dislocation–phase boundary interaction in two-phase microstructures. Considered is a single glide plane which intersects a fully coherent phase boundary. The glide plane is modelled in accordance with the PN model as an interface that splits the domain into two linear elastic media. Along the glide plane, bulk induced shear tractions balance glide plane shear tractions. The latter originate from the disregistry, the displacement across the glide plane, by a periodic shear traction–disregistry law based on the Frenkel sinusoidal. Individual phases are modelled by adopting phase specific material properties for the elastic medium and the glide plane. The full model is discretised spatially using the Finite Element Method (FEM) and in time by a backward Euler scheme. A highly non-linear problem results, which is solved numerically with an implicit Newton–Raphson algorithm.

Different to previous applications of the PN model, the dislocation induced bulk stresses are not given by an analytical solution, but are rather the result of solving, numerically, the governing equations. The model accounts for the interplay between dislocations, external boundary conditions and the phase boundary, while automatically taking the material properties into account. With only a few parameters it is a rather straightforward model that captures the mechanics of the edge dislocation–phase boundary interaction in an adequate manner. It presents therefore a suitable tool for studying dislocation transmission across phase boundaries.

In this study, a two-phase microstructure is considered, consisting of a soft phase that is flanked by a harder phase. The glide plane is assumed to be perpendicular to the phase boundaries throughout the two phases. On it, at the centre of the soft phase, lies a dislocation source that nucleates dislocation dipoles under sufficiently high shear stress. At the outer boundary of the problem domain a shear deformation is applied which triggers dislocation nucleation and drives the dislocations towards the phase boundary. This configuration is employed to study dislocation obstruction at the phase boundaries, pile-up formation and dislocation transmission into the harder phase.

After the problem statement and model formulation in Section 2, the numerical treatment is detailed in Section 3. In Section 4 the model is verified by a comparison with analytical solutions and a convergence study is performed. The influence of phase contrast and dislocation pile-ups on dislocation transmission is then studied in Section 5. The paper finishes with conclusions in Section 6.

2. Problem statement and model formulation

2.1. Problem description

To study the interaction between edge dislocations and a phase boundary the domain Ω with its boundary $\mathrm{\partial }\mathrm{\Omega }$, as sketched in Figure 1, is considered. Each material point in Ω is mapped by the position vector $\overrightarrow{x}$ in the two-dimensional Euclidean point space $R^2$ with the basis vectors ${\overrightarrow{e}}_x$ and ${\overrightarrow{e}}_y$. Domain Ω consists of Phase A that is flanked by Phase B, each with distinct material properties. Between both phases lies the shared phase boundary ${\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$ which is assumed normal to ${\overrightarrow{e}}_x$ (1) $\mathrm{\Omega }={\mathrm{\Omega }}^{\mathrm{A}}\cup {\mathrm{\Omega }}^{\mathrm{B}}$(1) (2) $\mathrm{\partial }\mathrm{\Omega }=\left(\mathrm{\partial }{\mathrm{\Omega }}^{\mathrm{A}}\setminus {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}\right)\cup \left(\mathrm{\partial }{\mathrm{\Omega }}^{\mathrm{B}}\setminus {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}\right)$(2) where ${\mathrm{\Omega }}^{\mathrm{A}}$ and ${\mathrm{\Omega }}^{\mathrm{B}}$ are the regions occupied by Phases A and B, respectively. In each phase $i\in \{\mathrm{A},\mathrm{B}\}$ lies a single glide plane ${\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i$ normal to ${\overrightarrow{e}}_y$ with the slip direction ${\overrightarrow{e}}_x$, splitting ${\mathrm{\Omega }}^i$ into ${\mathrm{\Omega }}_{\pm }^i$ (3) ${\mathrm{\Omega }}^i={\mathrm{\Omega }}_{+}^i\cup {\mathrm{\Omega }}_{-}^i$(3) (4) $\mathrm{\partial }{\mathrm{\Omega }}^i=\left(\mathrm{\partial }{\mathrm{\Omega }}_{+}^i\setminus {\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i\right)\cup \left(\mathrm{\partial }{\mathrm{\Omega }}_{-}^i\setminus {\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i\right)$(4)

Figure 1. Idealised dislocation motion in a two-phase microstructure with a single glide plane and a centred edge dislocation source.

The glide plane is initially defect-free and carries straight edge dislocation lines in the direction ${\overrightarrow{e}}_x\times {\overrightarrow{e}}_y$ only. Dislocations are nucleated from a dislocation source centred in Phase A. It emits a dislocation pair with Burgers vector $\pm b{\overrightarrow{e}}_x$ (a dislocation dipole) when the local resolved shear stress exceeds the critical value ${\tau }_{\mathrm{n}\mathrm{u}\mathrm{c}}$, representing the activation stress of a Frank–Read source. The resolved shear stress results from a shear deformation applied on the remote boundary $\mathrm{\partial }\mathrm{\Omega }$ and the dislocation induced internal stress field. A freshly nucleated dipole starts moving under the local Peach–Koehler force towards the phase boundary, where it will interact with this interface.

For the present study, the glide plane ${\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}$ is taken perpendicular to and continuous across a fully coherent phase boundary ${\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$. Interfacial decohesion is inhibited. Only dislocation motion from a soft Phase A to a stiff Phase B is considered. Accordingly, two possible dislocation–phase boundary interactions may exist: (1) dislocation obstruction and (2) dislocation transmission into Phase B. Evidently, dislocation obstruction may lead to the formation of a dislocation pile-up in which dislocations exert repelling forces on each other.

2.2. Governing equations

Under the absence of body forces the equilibrium of domain Ω is governed by the momentum balance in ${\mathrm{\Omega }}_{\pm }^i$, along ${\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i$ and along ${\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$:(5) $\mathrm{d}\mathrm{i}\mathrm{v}{\sigma }^i=\overrightarrow{0}\mathrm{i}\mathrm{n}{\mathrm{\Omega }}_{\pm }^i$(5) (6) ${\overrightarrow{T}}_{\mathrm{g}\mathrm{p},-}^i={\overrightarrow{T}}_{\mathrm{g}\mathrm{p},+}^i\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i$(6) (7) ${\overrightarrow{T}}_{\mathrm{p}\mathrm{b}}^{\mathrm{A}}={\overrightarrow{T}}_{\mathrm{p}\mathrm{b}}^{\mathrm{B}}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$(7) where ${\sigma }^i$ is the stress tensor in ${\mathrm{\Omega }}_{\pm }^i$, ${\overrightarrow{T}}_{\mathrm{g}\mathrm{p},\pm }^i={\sigma }_{\pm }^i\cdot {\overrightarrow{e}}_y$ are the tractions acting on $\mathrm{\partial }{\mathrm{\Omega }}_{\pm }^i\cap {\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i$ and ${\overrightarrow{T}}_{\mathrm{p}\mathrm{b}}^i={\sigma }^i\cdot {\overrightarrow{e}}_x$ the tractions acting on $\mathrm{\partial }{\mathrm{\Omega }}^i\cap {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$. The perfectly bonded and non-damaging phase boundary ${\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$ is enforced by the displacement continuity(8) ${\overrightarrow{u}}^{\mathrm{A}}={\overrightarrow{u}}^{\mathrm{B}}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$(8)

Following the concept of the classical Peierls–Nabarro (PN) model [1,33,34], linear elastic material behaviour is adopted in each subdomain ${\mathrm{\Omega }}_{\pm }^i$ at time t:(9) ${\sigma }^i={}^4{\mathbf{C}}^i:\epsilon$(9) (10) $\epsilon =\frac{1}{2}\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{u}+{\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{u}\right)}^T\right)$(10) where ${}^4{\mathbf{C}}^i$ is the fourth-order elasticity tensor associated with Phase i. Along the glide plane ${\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i$ the normal displacement is continuous, i.e.(11) $\left[\left[\overrightarrow{u}\right]\right]\cdot {\overrightarrow{e}}_y=\left({\overrightarrow{u}}_{+}-{\overrightarrow{u}}_{-}\right)\cdot {\overrightarrow{e}}_y=0$(11) but a tangential displacement jump Δ is allowed for that originates from dislocation glide:(12) $\mathrm{\Delta }=\left[\left[\overrightarrow{u}\right]\right]\cdot {\overrightarrow{e}}_x=\left({\overrightarrow{u}}_{+}-{\overrightarrow{u}}_{-}\right)\cdot {\overrightarrow{e}}_x$(12) Associated with Δ is a glide plane shear traction $T^i\left(\mathrm{\Delta }\right)$ which needs to balance the tangential part of the glide plane traction vector ${\overrightarrow{T}}_{\mathrm{g}\mathrm{p},\pm }^i$, i.e. $T^i\left(\mathrm{\Delta }\right)={\overrightarrow{T}}_{\mathrm{g}\mathrm{p},\pm }^i\cdot {\overrightarrow{e}}_x$. Note that the normal part of the traction vectors ${\overrightarrow{T}}_{\mathrm{g}\mathrm{p},\pm }^i$ follows directly from Equation (6(6) ${\overrightarrow{T}}_{\mathrm{g}\mathrm{p},-}^i={\overrightarrow{T}}_{\mathrm{g}\mathrm{p},+}^i\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i$(6) ) combined with the requirement of normal displacement continuity (11(11) $\left[\left[\overrightarrow{u}\right]\right]\cdot {\overrightarrow{e}}_y=\left({\overrightarrow{u}}_{+}-{\overrightarrow{u}}_{-}\right)\cdot {\overrightarrow{e}}_y=0$(11) ).

In this paper a periodic shear traction–disregistry law is employed as follows(13) $T^i\left(\mathrm{\Delta }\right)={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\dot{\mathrm{\Delta }}$(13) with the traction amplitude ${\hat{T}}^i$ and the damping coefficient ${\eta }^i$. The first term in Equation (13(13) $T^i\left(\mathrm{\Delta }\right)={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\dot{\mathrm{\Delta }}$(13) ) is the classical Frenkel sinusoidal function [37] and represents the periodic restoring force acting between the initially aligned atoms above and below the glide plane. We propose to add the second term in Equation (13(13) $T^i\left(\mathrm{\Delta }\right)={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\dot{\mathrm{\Delta }}$(13) ), with the damping coefficient ${\eta }^i$, to the classical expression. Using this term, the dynamic evolution of the dislocations is regularised. We find this to be necessary for numerical stability of the simulations, which without this term frequently fail to converge. As shown in Appendix 1, η can be related to the drag coefficient commonly used in DDD. The choice of a Frenkel sinusoidal function follows the classical PN model. As shown by Rice [42], Δ not only includes the inelastic, but also the elastic disregistry between the two adjacent atomic layers. Although this implies that the elastic response of that layer is accounted for twice, qualitative studies remain possible and valuable. A physically more realistic shear traction-disregistry relation is based on the Generalized Stacking Fault energy [39], which can be obtained from atomistic calculations [38,50,51] or from density-functional theory [52,53].

In the classical PN model a single static dislocation in a homogeneous and infinite medium is considered. The analytical solution for the disregistry profile can then be obtained as [1](14) $\mathrm{\Delta }\left(x\right)=-\frac{b}{\pi }{\tan }^{-1}\left(\frac{x-x_1}{\zeta }\right)+\frac{b}{2}$(14) where $x_1$ denotes the dislocation position and ζ the dislocation core half-width (in an isotropic crystal [39])(15) $\zeta =\frac{\mu b}{4\pi \left(1-\nu \right)\hat{T}}$(15) with the shear modulus μ and Poisson's ratio ν. For a two-phase microstructure no analytical solution is available and a numerical solution strategy is thus required.

2.3. Boundary conditions

Considering the model in Figure 1, it is apparent that the solution of the problem has the symmetry property $\overrightarrow{u}(-\overrightarrow{x})=-\overrightarrow{u}\left(\overrightarrow{x}\right)$ with respect to the dislocation source (at $\overrightarrow{x}=\overrightarrow{0}$). It is hence possible to reduce the model domain to the half-model presented in Figure 2 with model height $2H$, total length L and length D of Phase A. The symmetry boundary condition reads(16) $\overrightarrow{u}(0,y)=-\overrightarrow{u}(0,-y)\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{p}\mathrm{s}}=\left\{\right(0,y\left)\right|-H<y<H\}$(16) The shear deformation is imposed on $\mathrm{\partial }\mathrm{\Omega }\setminus {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{s}}$ in form of displacements that are consistent with a purely linear elastic model response – neglecting the glide plane – which would result under an applied shear load τ. Correspondingly, τ is used to parametrise the Dirichlet boundary conditions(17) $u_y=\frac{\tau }{{\mu }^{\mathrm{A}}}x\mathrm{o}\mathrm{n}\mathrm{\partial }{\mathrm{\Omega }}^{\mathrm{A}}\setminus \left({\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}\cup {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{s}}\right)$(17) (18) $u_y=\frac{\tau }{{\mu }^{\mathrm{B}}}(x-D)+\frac{\tau }{{\mu }^{\mathrm{A}}}D\mathrm{o}\mathrm{n}\mathrm{\partial }{\mathrm{\Omega }}^{\mathrm{B}}\setminus {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$(18) (19) $u_x=0\mathrm{o}\mathrm{n}\mathrm{\partial }\mathrm{\Omega }\setminus {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{s}}$(19)

Figure 2. Geometrical two-phase model for edge dislocation–phase boundary interaction with point-symmetry at x = 0, y = 0.

3. Numerical implementation

3.1. Finite element discretisation

The weak form of the balance of momentum (5(5) $\mathrm{d}\mathrm{i}\mathrm{v}{\sigma }^i=\overrightarrow{0}\mathrm{i}\mathrm{n}{\mathrm{\Omega }}_{\pm }^i$(5) ) and (6(6) ${\overrightarrow{T}}_{\mathrm{g}\mathrm{p},-}^i={\overrightarrow{T}}_{\mathrm{g}\mathrm{p},+}^i\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}^i$(6) ) with the constitutive relations (9(9) ${\sigma }^i={}^4{\mathbf{C}}^i:\epsilon$(9) ), (10(10) $\epsilon =\frac{1}{2}\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{u}+{\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{u}\right)}^T\right)$(10) ) and (13(13) $T^i\left(\mathrm{\Delta }\right)={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\dot{\mathrm{\Delta }}$(13) ) and the phase boundary constraints (8(8) ${\overrightarrow{u}}^{\mathrm{A}}={\overrightarrow{u}}^{\mathrm{B}}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$(8) ) and (7(7) ${\overrightarrow{T}}_{\mathrm{p}\mathrm{b}}^{\mathrm{A}}={\overrightarrow{T}}_{\mathrm{p}\mathrm{b}}^{\mathrm{B}}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}$(7) ) can be derived from the Generalized Clapeyrons formula [54, Chapter 5.14] as(20) ${\int }_{\mathrm{\Omega }}\left(\overrightarrow{\mathrm{\nabla }}{\overrightarrow{u}}_{\ast }\right):{}^4\mathbf{C}:\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{u}\right)\mathrm{d}\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}}{\mathrm{\Delta }}_{\ast }T\mathrm{d}\mathrm{\Gamma }=0$(20) Here, ${\overrightarrow{u}}_{\ast }$ is the test displacement and ${\mathrm{\Delta }}_{\ast }=\left[\right[{\overrightarrow{u}}_{\ast }\left]\right]\cdot {\overrightarrow{e}}_x$ the test disregistry. In order to solve Equation (20(20) ${\int }_{\mathrm{\Omega }}\left(\overrightarrow{\mathrm{\nabla }}{\overrightarrow{u}}_{\ast }\right):{}^4\mathbf{C}:\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{u}\right)\mathrm{d}\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}}{\mathrm{\Delta }}_{\ast }T\mathrm{d}\mathrm{\Gamma }=0$(20) ) numerically for the displacement field $\overrightarrow{u}(\overrightarrow{x},t)$ a discretisation in time and space is required. For the temporal discretisation the backward Euler method is adopted. Let $\overrightarrow{u}(\overrightarrow{x},t)$ be a known solution. An approximation of the velocity $\dot{\overrightarrow{u}}$ is given by(21) $\dot{\overrightarrow{u}}\approx \frac{\overrightarrow{u}\left(\overrightarrow{x},t+\theta \right)-\overrightarrow{u}\left(\overrightarrow{x},t\right)}{\theta }=\frac{\overrightarrow{u}-{\overrightarrow{u}}^t}{\theta }$(21) where θ is the time step. Hence, the shear traction-disregistry law (13(13) $T^i\left(\mathrm{\Delta }\right)={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\dot{\mathrm{\Delta }}$(13) ) reads in discretised form(22) $T^i={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\frac{\mathrm{\Delta }-{\mathrm{\Delta }}^t}{\theta }$(22) This leads to a non-linearity of the weak form which is solved implicitly with a Newton–Raphson scheme at each time $t+\theta$:(23) $\begin{array}{rl}& {\int }_{\mathrm{\Omega }}\left(\overrightarrow{\mathrm{\nabla }}{\overrightarrow{u}}_{\ast }\right):{}^4\mathbf{C}:\left(\overrightarrow{\mathrm{\nabla }}\delta \overrightarrow{u}\right)\mathrm{d}\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}}{\mathrm{\Delta }}_{\ast }\delta T\mathrm{d}\mathrm{\Gamma }=\\ & -{\int }_{\mathrm{\Omega }}\left(\overrightarrow{\mathrm{\nabla }}{\overrightarrow{u}}_{\ast }\right):{}^4\mathbf{C}:\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{u}\right)\mathrm{d}\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}}{\mathrm{\Delta }}_{\ast }T\mathrm{d}\mathrm{\Gamma }\end{array}$(23) where $\overrightarrow{u}$ and T denote the computed state variables at the previous iteration k and $\delta \overrightarrow{u}$ is the iterative correction to be solved for. $\delta T$ is related to $\delta \overrightarrow{u}$ via Equation (12(12) $\mathrm{\Delta }=\left[\left[\overrightarrow{u}\right]\right]\cdot {\overrightarrow{e}}_x=\left({\overrightarrow{u}}_{+}-{\overrightarrow{u}}_{-}\right)\cdot {\overrightarrow{e}}_x$(12) ) and(24) $\delta T=M\delta \mathrm{\Delta }$(24) with, from Equation (22(22) $T^i={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\frac{\mathrm{\Delta }-{\mathrm{\Delta }}^t}{\theta }$(22) ),(25) $M={\left.\frac{\mathrm{d}T}{\mathrm{d}\mathrm{\Delta }}\right|}_{\mathrm{\Delta }={\mathrm{\Delta }}_k}={\hat{T}}^i\frac{2\pi }{b^i}\cos \left(\frac{2\pi {\mathrm{\Delta }}_k}{b^i}\right)+\frac{{\eta }^i}{\theta }$(25)

For the spatial discretisation use is made of the FEM with Lagrange shape functions. Adopting the Voigt notation, Equation (23(23) $\begin{array}{rl}& {\int }_{\mathrm{\Omega }}\left(\overrightarrow{\mathrm{\nabla }}{\overrightarrow{u}}_{\ast }\right):{}^4\mathbf{C}:\left(\overrightarrow{\mathrm{\nabla }}\delta \overrightarrow{u}\right)\mathrm{d}\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}}{\mathrm{\Delta }}_{\ast }\delta T\mathrm{d}\mathrm{\Gamma }=\\ & -{\int }_{\mathrm{\Omega }}\left(\overrightarrow{\mathrm{\nabla }}{\overrightarrow{u}}_{\ast }\right):{}^4\mathbf{C}:\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{u}\right)\mathrm{d}\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}}{\mathrm{\Delta }}_{\ast }T\mathrm{d}\mathrm{\Gamma }\end{array}$(23) ) can be rewritten as(26) $\left({\int }_{\mathrm{\Omega }}{\underset{\_}{B}}^T\underset{\_}{C}\underset{\_}{B}\mathrm{d}\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}}{\underset{\_}{P}}^T{\underset{\sim }{\overline{N}}}^T\underset{\_}{M}\underset{\sim }{\overline{N}}\underset{\_}{P}\mathrm{d}\mathrm{\Gamma }\right)\delta \underset{\sim }{u}=-{\underset{\sim }{f}}^{\mathrm{i}\mathrm{n}\mathrm{t}}$(26) with(27) ${\underset{\sim }{f}}^{\mathrm{i}\mathrm{n}\mathrm{t}}={\int }_{\mathrm{\Omega }}{\underset{\_}{B}}^T\underset{\_}{C}\underset{\_}{B}\mathrm{d}\mathrm{\Omega }\underset{\sim }{u}+{\int }_{{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}}{\underset{\_}{P}}^T{\underset{\sim }{\overline{N}}}^T\underset{\sim }{T}\mathrm{d}\mathrm{\Gamma }$(27) Here, $\underset{\sim }{N}$ and ${\displaystyle \underset{\sim }{\overline{N}}}$ are the shape functions for bulk and glide plane, respectively, $\underset{\_}{B}$ contains the spatial derivatives of the shape functions in ${\displaystyle \underset{\sim }{N}}$, $\underset{\_}{C}$ is the Voigt elasticity matrix, $M$ is the glide plane tangent stiffness, ${\displaystyle \underset{\sim }{u}}$ are the nodal displacements and ${\displaystyle \underset{\sim }{T}}$ are the nodal glide plane tractions. $\underset{\_}{P}$ is a matrix such that in accordance with Equation (12(12) $\mathrm{\Delta }=\left[\left[\overrightarrow{u}\right]\right]\cdot {\overrightarrow{e}}_x=\left({\overrightarrow{u}}_{+}-{\overrightarrow{u}}_{-}\right)\cdot {\overrightarrow{e}}_x$(12) )(28) $\underset{\sim }{\mathrm{\Delta }}=\underset{\_}{P}\underset{\sim }{u}$(28) holds, e.g. for the linear quadrilateral interface element(29) $\left[\begin{array}{c}{\mathrm{\Delta }}^{1-4}\\ {\mathrm{\Delta }}^{2-3}\end{array}\right]=\left[\begin{array}{cccc}-1& 0& 0& 1\\ 0& -1& 1& 0\end{array}\right]\left[\begin{array}{c}{u}_x^1\\ u_x^2\\ u_x^3\\ u_x^4\end{array}\right]$(29)

with nodes 1 and 2 of ${\mathrm{\Omega }}_{-}$ initially coincident with nodes 4 and 3 of ${\mathrm{\Omega }}_{+}$. Tie constraints for the description of the point-symmetry (16(16) $\overrightarrow{u}(0,y)=-\overrightarrow{u}(0,-y)\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{p}\mathrm{s}}=\left\{\right(0,y\left)\right|-H<y<H\}$(16) ) and the glide plane constraint ${\mathrm{\Delta }}_n=0$ are adopted in a standard FEM manner.

Each domain ${\mathrm{\Omega }}_{\pm }^i$ is discretised by bilinear quadrilateral elements with four Gauss integration points. For the glide plane linear interface elements of zero thickness with two Gauss integration points are used.

3.2. Dislocation nucleation

As stated in Section 2, a dislocation source is used to nucleate dislocation dipoles each time the resolved shear stress ${\tau }^{\mathrm{r}\mathrm{e}\mathrm{s}}$ at the position of the source exceeds the critical value ${\tau }_{\mathrm{n}\mathrm{u}\mathrm{c}}$. The numerical implementation entails: (1) checking the condition ${\tau }^{\mathrm{r}\mathrm{e}\mathrm{s}}\ge {\tau }_{\mathrm{n}\mathrm{u}\mathrm{c}}$; (2) introducing the positive part of the dipole by perturbing the solution such that the disregistry reflects its presence. The position x where the dislocation (dipole) is introduced is selected to be close to its stable equilibrium position. The latter is motivated by the computational expense that would otherwise be required to propagate a dislocation (dipole) from its nucleation position towards its final equilibrium position. Focusing on the process of dislocation transmission and not dislocation motion towards the pile-up, a nucleation close to the static equilibrium is desirable.

For the nucleation criterion, the (converged) mean resolved shear stress ${\tau }_{\mathrm{m}}^{\mathrm{r}\mathrm{e}\mathrm{s}}$ is computed in the vicinity of the dislocation source using Equation (22(22) $T^i={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\frac{\mathrm{\Delta }-{\mathrm{\Delta }}^t}{\theta }$(22) ):(30) ${\tau }_{\mathrm{m}}^{\mathrm{r}\mathrm{e}\mathrm{s}}=\frac{1}{L_0}{\int }_0^{L_0}{T}^{\mathrm{A}}\mathrm{d}x$(30) The disregistry rate $(\mathrm{\Delta }-{\mathrm{\Delta }}^t)/\theta$ plays only a negligible role since either a potential dislocation is very close to the source and the sinusoidal term is highly negative (no dislocation nucleation) or the dislocation is distant and the disregistry rate insignificant (potential dislocation nucleation). The averaging distance $L_0$ in (30(30) ${\tau }_{\mathrm{m}}^{\mathrm{r}\mathrm{e}\mathrm{s}}=\frac{1}{L_0}{\int }_0^{L_0}{T}^{\mathrm{A}}\mathrm{d}x$(30) ) is defined according to the critical loop size for a Volterra dislocation dipole [2]:(31) $L_0=\frac{1}{2}\frac{{\mu }^{\mathrm{A}}{b}^{\mathrm{A}}}{2\pi (1-\nu ){\tau }_{\mathrm{n}\mathrm{u}\mathrm{c}}}$(31) Once the nucleation criterion is fulfilled, a positive dislocation is introduced (and hence due to the point-symmetry also the dipole) by perturbing the solution and letting the system equilibrate. For numerical stability, this is done in a two-step procedure. In the first step, the disregistry profile of a single dislocation in an infinite medium, according to Equation (14(14) $\mathrm{\Delta }\left(x\right)=-\frac{b}{\pi }{\tan }^{-1}\left(\frac{x-x_1}{\zeta }\right)+\frac{b}{2}$(14) ), is added to the existing disregistry solution:(32) $\delta u_{x,\pm }=\pm \frac{b^{\mathrm{A}}}{2}\left(-\frac{1}{\pi }{\tan }^{-1}\left(\frac{x-L_{\mathrm{n}\mathrm{u}\mathrm{c}}}{{\zeta }^{\mathrm{A}}}\right)+\frac{1}{2}\right)\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}$(32) where $\delta u_{x,\pm }$ are the displacements belonging to $\mathrm{\partial }{\mathrm{\Omega }}_{\pm }\cap {\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}$, and $L_{\mathrm{n}\mathrm{u}\mathrm{c}}$ the nucleation position of the new dislocation (see below). The resulting disregistry profile is kept fixed while the surrounding elastic region is allowed to equilibrate. In the second step the disregistry is released and additional equilibrium iterations allow the model to find an equilibrium solution featuring the new dislocation. Once such a new equilibrium state of the complete system has been established, the simulation proceeds to the next time step by incrementing the applied displacements on the external boundary.

To find a suitable nucleation position $L_{\mathrm{n}\mathrm{u}\mathrm{c}}$, the position x is sought at which the current glide plane shear tractions $T\left(x\right)$ and the shear stresses ${\sigma }_{xy}^{\pm }\left(x\right)$ of the to be nucleated positive/negative dislocation along the glide plane (i.e. for y = 0), situated at $\pm x$, are in equilibrium:(33) $T\left(x\right)+{\sigma }_{xy}^{+}\left(x\right)+{\sigma }_{xy}^{-}\left(x\right)=0$(33) Neglecting the influence of the left boundary, we employ as an approximation for ${\sigma }_{xy}^{\pm }\left(x\right)$ the stress field for a Volterra dislocation in a two-phase elastic medium as given by Head [55, Section 10]. Mathematically, Equation (33(33) $T\left(x\right)+{\sigma }_{xy}^{+}\left(x\right)+{\sigma }_{xy}^{-}\left(x\right)=0$(33) ) may have n+1 stable equilibrium positions, where n is the number of dislocations (dislocation dipoles) introduced so far, and one unstable equilibrium position in $L_0\le x\le D$. However, for n > 0, only one stable equilibrium position lies between $L_0$ and the position of the most recently introduced dislocation, $x_n$. Hence, an upper bound for the dislocation nucleation position is introduced as(34) $L_{\mathrm{m}\mathrm{a}\mathrm{x}}=\left\{\begin{array}{ll}D& \mathrm{f}\mathrm{o}\mathrm{r}n=0\\ x_n& \mathrm{f}\mathrm{o}\mathrm{r}n>0\end{array}\right.$(34) The unstable equilibrium position characterises the point at which a small perturbation leads either to annihilation of the dipole or to dislocation motion towards the phase boundary into the stable equilibrium position. Thus, the nucleation position is chosen as the position $x\in [L_0,L_{\mathrm{m}\mathrm{a}\mathrm{x}}]$ where Equation (33(33) $T\left(x\right)+{\sigma }_{xy}^{+}\left(x\right)+{\sigma }_{xy}^{-}\left(x\right)=0$(33) ) holds and that is the closest to the phase boundary (x = D).

4. Verification and convergence study

4.1. Single dislocation in a homogeneous medium

The verification of the implemented Peierls–Nabarro continuum model is accompanied by a mesh convergence study in a single dislocation framework. The numerically obtained stress field of the dislocation in a homogeneous material is compared with the analytical solution presented in [1], i.e. a true verification step. Considered is a model of dimensions $L\times 2H=100b\times 100b$. These dimensions, here and in what follows, have been selected sufficiently large for the results to depend only marginally on them. Uniform meshes are employed, with element sizes of $h=b/8,b/4,b/2,b,2b$. For the linear elastic behaviour of the single Phase A considered, a shear modulus of ${\mu }^{\mathrm{A}}=80.77\mathrm{G}\mathrm{P}\mathrm{a}$ and a Poisson's ratio of ${\nu }^{\mathrm{A}}=0.3$ are adopted. The glide plane amplitude is chosen in alignment with the small strain limit of linear elasticity as ${\hat{T}}^{\mathrm{A}}={\mu }^{\mathrm{A}}{b}^{\mathrm{A}}/2\pi d^{\mathrm{A}}$, where $d^{\mathrm{A}}$ is the interplanar spacing [1]. For simplicity $b^{\mathrm{A}}$ and $d^{\mathrm{A}}$ are assumed to be equal such that ${\hat{T}}^{\mathrm{A}}={\mu }^{\mathrm{A}}/2\pi$. The Burgers vector is $b=b^{\mathrm{A}}=0.25\mathrm{n}\mathrm{m}$. Following Appendix 1, the viscosity is defined by the relationship ${\eta }^{\mathrm{A}}=B/b$ in terms of the Drag coefficient $B$. The critical nucleation stress is set to ${\tau }_{\mathrm{n}\mathrm{u}\mathrm{c}}=0.048{\hat{T}}^{\mathrm{A}}$.

With the application of a constant shear load of $\tau =0.054{\hat{T}}^{\mathrm{A}}$ dislocation nucleation is triggered, enforced here at $L_{\mathrm{n}\mathrm{u}\mathrm{c}}\approx 40.5b$. Subsequently, the dislocation moves in accordance with the BVP towards the rigid surface at x = L and settles in its equilibrium position $x_1$. Here, the dislocation position $x_1$ is defined as the coordinate where $\mathrm{\Delta }=b/2$. Figure 3(a) shows the final dislocation position $x_1$ as a function of the element size h. The convergence behaviour is presented in Figure 3(b) where the absolute error with respect to the finest discretisation considered, $ϵ=\left|x_1\right(h=b/8)-x_1(h\left)\right|$, is plotted against the mesh size h. A strong influence of the discretisation on the final dislocation position is evident for coarse discretisations. Mesh sizes of $h\ge b$ lead to an absolute error of $ϵ\approx b$ and above. For element sizes $h\le b/2$, however, the dislocation position quickly converges, resulting in an accuracy on the order of ${10}^{-6}b$ or better. A (too) coarse discretisation is unable to represent the gradient of the disregistry and leads hence to a large deviation of the dislocation position.

Figure 3. Influence of the discretisation on the dislocation equilibrium position $x_1$.

Studying the dislocation induced stress field, one recovers the well-known 2D stress field, as shown in Figure 4 for $h=b/8$. To visualise the discretisation influence on the stress field representative line plots are extracted. In Figure 5(a,b) the shear stress ${\sigma }_{xy}$ is plotted as a function of the position $x-x_1$ along the x-axis at y = 0 and $y=b/16$, respectively. In both figures, the analytical solution for the infinite domain is added as a reference. The numerical curves have been obtained by evaluating the stresses along the line $y=0^{+}$, i.e. along the lower boundaries of the elastic elements just above the glide plane. No intra- or extrapolation has been employed, other than that of the finite element discretisation of the displacement field. For ${\sigma }_{xy}(x,y=0)$, the glide plane shear traction T provides a rather accurate measure. The shear stress is captured adequately for mesh sizes of $h\le b/2$. For the discretisation h = b, small discrepancies emerge, while for h = 2b the results become too inaccurate. This is in alignment with the convergence study of the dislocation position and shows that the proper quantification of the shear tractions plays a major role for accuracy of the dislocation evolution. The larger deviations that are observed for ${\sigma }_{xy}(x,y=b/16)$ (Figure 5 (b)) illustrate the significant impact of the spatial discretisation on the shear stress ${\sigma }_{xy}$. Besides a vertical offset, also stress jumps across the elements arise. Both offset and stress jumps are reduced by refining the mesh size.

Figure 4. 2D stress field of a single dislocation computed with element size $h=b/8$.

Figure 5. Shear stress ${\sigma }_{xy}(x-x_1,y)$ at a) y = 0 and b) $y=b/16$ for different discretisation sizes, along with the analytical solution.

Similar stress jumps are observed for ${\sigma }_{xx}(x=x_1,y)$ and ${\sigma }_{yy}(x=x_1,y)$, with significant irregularities in the latter one, as visualised in Figure 6(a,b), respectively. For the sake of clarity only the stresses for $h=b/8,b/2,2b$  are presented. Note that ${\sigma }_{xx}$ exhibits smaller fluctuations due to its smoother stress gradient in comparison to ${\sigma }_{xy}$ and ${\sigma }_{yy}$. Nevertheless, the solution converges towards the analytical solution for increasingly fine meshes. Different to the glide plane traction T and the related dislocation position $x_1$ where an element size of $h\le b/2$ suffices, a rather fine mesh with $h\le b/8$ is required for a high accuracy in the stress response.

Figure 6. Normal stresses a) ${\sigma }_{xx}(x=x_1,y)$ and b) ${\sigma }_{yy}(x=x_1,y)$ for different discretisation sizes along with the analytical solution.

4.2. Two-phase microstructure and dislocation transmission

For the study of dislocation obstruction and dislocation transmission, a two-phase microstructure is considered with the grain size of Phase A 2D = L. The properties of Phase A are identical to those introduced in Section 4. For Phase B the following material properties are adopted: ${\mu }^{\mathrm{B}}=7{\mu }^{\mathrm{A}}/6$, ${\hat{T}}^{\mathrm{B}}=7{\hat{T}}^{\mathrm{A}}/6$, ${\nu }^{\mathrm{B}}={\nu }^{\mathrm{A}}$, $b^{\mathrm{B}}=b$ and ${\eta }^{\mathrm{B}}={\eta }^{\mathrm{A}}$. The shear load is applied through a stress rate $\dot{\tau }$, here given by $\dot{\tau }=2.59\cdot {10}^{-5}{\hat{T}}^{{\mathrm{A}}^2}/B$. This rather high loading rate is necessary to sufficiently activate the viscous term in (13(13) $T^i\left(\mathrm{\Delta }\right)={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\dot{\mathrm{\Delta }}$(13) ), in order to stabilise the simulations; for significantly lower loading rates, divergence of the Newton iterations was sometimes observed.

Figure 7(a) shows the dislocation evolution for $h=b/2$ in terms of the dislocation position $x_1$ (horizontal axis) as a function of the applied shear load τ (vertical axis). Phase B is visualised by the grey domain. After nucleation at $\tau ={\tau }_{\mathrm{n}\mathrm{u}\mathrm{c}}$ the dislocation moves under the influence of the applied shear load towards the phase boundary, where it is obstructed. The physics governing the dislocation obstruction stems from the introduced phase contrast itself. The stiffer Phase B induces an image stress field that repels the approaching dislocation. Naturally, the shear load needs to be increased further to overcome the repulsion and to push the dislocation closer to the interface. Eventually, the applied shear load reaches a critical value, the external transmission stress ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$, here ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}\approx 0.15{\hat{T}}^{\mathrm{A}}$. At this point, at which the dislocation is situated at the phase boundary ($x_1=D$), the local shear stress acting on the dislocation is in equilibrium with the transmission resistance. A small increase of the shear load drives the dislocation rapidly into Phase B. The slow and linear motion of the dislocation prior to transmission resembles quasi-static solutions. It is therefore expected that the transmission process is only slightly influenced by the magnitude of the applied stress rate. To study the local stress build-up during dislocation transmission, the maximum normal stress $max\left({\sigma }_{xx}\right)$ at the boundary x = D is evaluated. Following the analytical solution the maximum of ${\sigma }_{xx}$ occurs in ${\mathrm{\Omega }}_{-}$ at y = 0 [1]. This point is shared by two elements, one lying in Phase A and one in Phase B. Hence, two stress evaluations are available at the respective nodal position, extrapolated from the element in Phase A and Phase B, here denoted as ${\sigma }_{xx,\mathrm{A}}$ and ${\sigma }_{xx,\mathrm{B}}$. Figure 7(b) displays these stresses with the elemental averaged nodal stress ${\sigma }_{xx,\mathrm{A}\mathrm{B}}=({\sigma }_{xx,\mathrm{A}}+{\sigma }_{xx,\mathrm{B}})/2$ (vertical axis) as a function of the dislocation position (horizontal axis). Again, the grey area indicates Phase B. The deviation visible in ${\sigma }_{xx,\mathrm{A}\mathrm{B}}$ originates from the spatial discretisation with a coarse element size of $h=b/2$. Clearly, a high interface stress emerges only when a dislocation reaches the direct vicinity of the interface, with its maximum at $x_1=D$.

Figure 7. (a) Dislocation position $x_1$ with respect to the applied shear load τ; (b) nodal stresses extrapolated from the element in Phase A/B ${\sigma }_{xx,\mathrm{A}/\mathrm{B}}(D,0)$ and the element averaged nodal stress ${\sigma }_{xx,\mathrm{A}\mathrm{B}}$ as a function of the dislocation position.

The mesh convergence results for dislocation transmission are presented in Figure 8. Figure 8(a) plots the external transmission stress ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ (vertical axis) against the element size h (horizontal axis). Although h = b exhibits an inaccuracy in terms of the equilibrium position (see Figure 3) it provides with ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}\approx 0.15{\hat{T}}^{\mathrm{A}}$ a reliable identification of the external transmission stress. In Figure 8(b) the convergence of the maximum interface stresses ${\sigma }_{xx,\mathrm{A}/\mathrm{B}}(x=D,y=0)$ and their average ${\sigma }_{xx,\mathrm{A}\mathrm{B}}(x=D,y=0)$ are depicted as a function of h. For a good approximation, an element size of $h=b/8$ is required – equivalent to the stress field approximation in a homogeneous medium (see Figures 5 and 6). With ${\sigma }_{xx,\mathrm{A}/\mathrm{B}}=4.4{\hat{T}}^{\mathrm{A}}$ the normal stress clearly exceeds the conventional approximation of the theoretical strength E/10. This is due to the subatomic discretisation in the continuum model where the conventional approximation no longer applies.

Figure 8. Influence of the mesh size h on a) the external transmission stress ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ and b) the maximum nodal stresses extrapolated from the element in Phase A/B ${\sigma }_{xx,\mathrm{A}/\mathrm{B}}(D,0)$ and the element averaged nodal stress ${\sigma }_{xx,\mathrm{A}\mathrm{B}}$.

5. Influence of phase contrast and pile-ups on dislocation transmission

5.1. Transmission in a single dislocation case

As illustrated in Section 4.2, a phase contrast provides a natural source for dislocation obstruction. The aim of this section is to investigate the influence of different phase contrasts on the external transmission stress ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ in a single dislocation framework. For this purpose model dimensions of $L\times 2H=400b\times 500b$ and D = 200b are adopted. The element size is set to h = b which provides adequate results for ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ (see Section 4.2). Material properties are: ${\mu }^{\mathrm{A}}=80.77\mathrm{G}\mathrm{P}\mathrm{a}$, ${\nu }^{\mathrm{B}}={\nu }^{\mathrm{A}}=0.3$, ${\hat{T}}^{\mathrm{A}}={\mu }^{\mathrm{A}}/2\pi$, $b^{\mathrm{B}}=b^{\mathrm{A}}=b=0.25\mathrm{n}\mathrm{m}$, ${\eta }^{\mathrm{A}}={\eta }^{\mathrm{B}}=B/b$ $B={10}^{-4}\mathrm{P}\mathrm{a}\mathrm{s}$ and ${\tau }_{\mathrm{n}\mathrm{u}\mathrm{c}}=0.012{\hat{T}}^{\mathrm{A}}$. The shear modulus ${\mu }^{\mathrm{B}}$ and the glide plane traction amplitude ${\hat{T}}^{\mathrm{B}}$ are varied in the case studies examined. The shear load is applied at the rate $\dot{\tau }=0.32\cdot {10}^{-5}{\hat{T}}^{{\mathrm{A}}^2}/B$.

In terms of phase contrast, three selected cases are considered (k > 1):

• Case I: ${\mu }^{\mathrm{B}}=k{\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}={\hat{T}}^{\mathrm{A}}$

• Case II: ${\mu }^{\mathrm{B}}={\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}=k{\hat{T}}^{\mathrm{A}}$

• Case III: ${\mu }^{\mathrm{B}}=k{\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}=k{\hat{T}}^{\mathrm{A}}$

These cases are chosen such that two different sources of dislocation obstruction are triggered. One relates to the repelling long-range image stress field induced by the difference in elastic properties. The second source stems from an increased traction amplitude $\hat{T}$ in Phase B which causes a short-range image stress field. The results are shown in Figure 9 in terms of the external transmission stress ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ as a function of the phase contrast k.

Figure 9. Influence of phase contrast k on the external transmission stress for Case I: ${\mu }^{\mathrm{B}}=k{\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}={\hat{T}}^{\mathrm{A}}$; Case II: ${\mu }^{\mathrm{B}}={\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}=k{\hat{T}}^{\mathrm{A}}$; Case III: ${\mu }^{\mathrm{B}}=k{\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}=k{\hat{T}}^{\mathrm{A}}$. The dashed line represents the scaling factor for the image stress of an edge dislocation on itself ${\tau }_{xy}(x-D)/b$ [55].

For the pure elasticity contrast, in Case I, a trend is obtained that is comparable to the scaling factor for the image stress of an edge dislocation at x on itself as derived for a Volterra dislocation by Head [55, Section 10]. This trend is included in Figure 9 as a dashed curve. Despite its premise of a Volterra dislocation it presents a fairly good approximation for ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{I}}$. Both Cases II and III exhibit a linear trend between ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ and k with ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{I}\mathrm{I}\mathrm{I}}>{\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{I}\mathrm{I}}>{\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{I}}$. This shows that for a single dislocation the glide plane contrast and the related repulsive short-range image stress field has a higher impact on dislocation obstruction than the pure elasticity contrast with its long-range image stress field. Furthermore, since ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{I}\mathrm{I}\mathrm{I}}>{\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{I}\mathrm{I}}+{\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{I}}$, the influence of glide plane and elasticity contrast cannot be decoupled, i.e. the elasticity contrast increases the short-range stress contribution that is invoked by the higher traction amplitude.

5.2. Transmission in a dislocation pile-up case

So far the case of a single dislocation dipole was considered. This is next extended to multiple dislocations, allowing for the formation of dislocation pile-ups, which will impact dislocation transmission. Model settings are used as in the single dislocation case. Phase contrasts are adopted in alignment with Cases I–III, as defined in Section 5.1, and k = 2.

During the evolution process dislocations are nucleated at positions following from Section 3.2. Figure 10 displays the dislocation evolution and the eventually triggered dislocation transmission for the different phase contrasts. No upper limits for the number of dislocations are used. The position of each individual dislocation $x_j$ (horizontal axis), where $\mathrm{\Delta }=(j-1/2)b$, is plotted by a single line against the applied shear load τ (vertical axis). Upon increasing the applied shear load, dislocations are nucleated one by one and pile-up at the phase boundary. Here, the absence of a long-range image stress field in Case II becomes apparent as the leading dislocation moves directly into the phase boundary. By further increasing the shear load, the leading dislocation gets ultimately transmitted into Phase B. Transmission is recorded for Case I with four dislocations at ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}\approx 0.120{\hat{T}}^{\mathrm{A}}$, for Case II with four dislocations at ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}\approx 0.116{\hat{T}}^{\mathrm{A}}$ and for Case III with seven dislocations at ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}\approx 0.221{\hat{T}}^{\mathrm{A}}$. The nucleation of dislocations in the vicinity of existing dislocations triggers a momentarily strong driving force on the current pile-up – i.e. a small kink in the remaining dislocation evolution curves.

Figure 10. Dislocation position with respect to the applied shear load for each dislocation for Case I: ${\mu }^{\mathrm{B}}=2{\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}={\hat{T}}^{\mathrm{A}}$ (solid line); Case II: ${\mu }^{\mathrm{B}}={\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}=2{\hat{T}}^{\mathrm{A}}$ (dashed line); Case III: ${\mu }^{\mathrm{B}}=2{\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}=2{\hat{T}}^{\mathrm{A}}$ (dash-dot line). The numbers in the diagram refer to the dislocation number j.

In all cases, the phase contrast and the grain size $2D$ induce an upper limit on the number of dislocations in the pile-up. While the phase contrast defines the resistance against dislocation transmission and relates hence to the external transmission stress ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$, the grain size limits the number of stable dislocation dipoles that can be accommodated in Phase A at $\tau ={\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$. For Case II, where no elasticity related long-range image stress field is present, an analytical expression for this upper limit can be found. Following [1] the number N of positive dislocations in a dipole pile-up of length l within a homogeneous medium reads:(35) $N=\frac{(1-\nu )l\tau }{\mu b}$(35) with a resolved shear stress acting on the leading dislocation of(36) ${\tau }_1^{\mathrm{r}\mathrm{e}\mathrm{s}}=\frac{\pi }{4}N\tau$(36) Combining Equations (35(35) $N=\frac{(1-\nu )l\tau }{\mu b}$(35) ) and (36(36) ${\tau }_1^{\mathrm{r}\mathrm{e}\mathrm{s}}=\frac{\pi }{4}N\tau$(36) ) yields(37) $N=\sqrt{\frac{4}{\pi }\frac{(1-\nu )l{\tau }_1^{\mathrm{r}\mathrm{e}\mathrm{s}}}{\mu b}}$(37) Projecting this relation upon the presented model, the natural upper limit of N can be approximated when a dipole pile-up is fully formed at the point of transmission. Accordingly, ${\tau }_1^{\mathrm{r}\mathrm{e}\mathrm{s}}$ is taken equal to ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ of a single dislocation and l is the grain size of Phase A, $2D$. In the present case N = 4 results (after rounding down), which is in agreement with the simulation. If the size of Phase A is so large that transmission occurs before the dipole pile-up is fully formed, a mixture between single pile-up and dipole pile-up is captured. A straightforward analytical expression is then not possible. An analytical expression also fails for Cases I and III where the presence of dislocations gives rise to the repulsive long-range image stress field and hence an increased dislocation obstruction.

The influence of the number of dislocations on the external transmission stress ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ is studied by varying the number of dislocations $N\le 4$ for Cases I and II and $N\le 7$ for Case III, by de-activating the dislocation source once these respective numbers of dislocations have been introduced. The results are presented in Figure 11. Although a strong decay in external transmission stress is obvious, the expected reduction towards zero is not apparent. The reason is twofold. (i) The Dirichlet boundary conditions (17(17) $u_y=\frac{\tau }{{\mu }^{\mathrm{A}}}x\mathrm{o}\mathrm{n}\mathrm{\partial }{\mathrm{\Omega }}^{\mathrm{A}}\setminus \left({\mathrm{\Gamma }}_{\mathrm{p}\mathrm{b}}\cup {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{s}}\right)$(17) )–(19(19) $u_x=0\mathrm{o}\mathrm{n}\mathrm{\partial }\mathrm{\Omega }\setminus {\mathrm{\Gamma }}_{\mathrm{p}\mathrm{s}}$(19) ) are applied under the assumption of linear elasticity, neglecting the influence of slip along ${\mathrm{\Gamma }}_{\mathrm{g}\mathrm{p}}$. While this effect is marginal for a single dislocation, its impact grows with the number of dislocations. Thus, ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ is increasingly overestimated. (ii) At the point of transmission, the dislocation pile-ups are not in their static equilibrium, yet. The resolved shear stress on the leading dislocation is lower than for a fully developed dislocation pile-up. Hence, for dislocation transmission the applied shear load needs to be increased further than expected for a fully developed pile-up. This effect, although minor, cannot be avoided as without the viscous regularisation the simulations are numerically unstable.

Figure 11. External transmission stress with respect to the number of dislocations for Case I: ${\mu }^{\mathrm{B}}=2{\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}={\hat{T}}^{\mathrm{A}}$; Case II: ${\mu }^{\mathrm{B}}={\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}=2{\hat{T}}^{\mathrm{A}}$; Case III: ${\mu }^{\mathrm{B}}=2{\mu }^{\mathrm{A}}$ and ${\hat{T}}^{\mathrm{B}}=2{\hat{T}}^{\mathrm{A}}$.

Comparing Cases I and II, a growing influence of the elasticity induced long-range image stress field becomes obvious. As more dislocations are introduced, the repulsive long-range image stress field stacks up. Consequently, dislocation obstruction gets stronger than for the pure glide plane contrast where dislocation obstruction arises from the short-range image stress field caused by solely the transmitting dislocation. Hence, Case I exhibits a smaller decay in external transmission stress than Case II.

6. Conclusion

A numerical 2D model for edge dislocation–phase boundary interaction has been presented. It consists of two linear elastic phases with an embedded glide plane according to the Peierls–Nabarro model. To solve the governing equations, the model is discretised spatially by FEM and in time by a backward Euler scheme. This approach provides a natural interplay between dislocations, phase boundary and external boundary conditions with only a limited number of parameters. It is therefore a suitable method for studying the stress build-up at phase boundaries due to dislocation pile-ups.

This paper considered a soft Phase A that is flanked by a harder Phase B. At the centre of Phase A, a dislocation source emits edge dislocation dipoles under a sufficiently high local shear stress. A single glide plane, perpendicular to a fully coherent and non-damaging phase boundary, was analysed. A shear load was applied in the form of Dirichlet boundary conditions.

To verify the presented model, the stress field of a single dislocation within a homogeneous material was studied. It could be shown that for a sufficiently fine discretisation ($h\le b/8$), the numerical solution converges towards the analytical solution. Similarly, the adequate spatial discretisation for capturing dislocation motion and dislocation transmission in a two-phase microstructure was determined, with $h\le b/2$.

For a first qualitative study, the influence of phase contrast and the number of dislocations in a dipole pile-up on dislocation transmission was determined. Due to the presence of a phase contrast, a resistance against dislocation motion – from soft to hard phase – emerges. Different sources of resistance were identified as (i) a long-range image stress field due to the difference in elastic properties and (ii) a short-range image stress field arising from the increased traction amplitude of the PN model.

In order to differentiate between the obstruction due to the long- and short-range stress fields, a transmission study was performed exploring three different kinds of phase contrast: pure phase contrast in elastic properties, pure phase contrast in the glide plane traction amplitude and a combination of both. The following insights were gained. (i) Dislocation obstruction due to pure elasticity contrast scales in a comparable manner as the image shear stress of a Volterra edge dislocation on itself as derived by Head [55]. (ii) For a pure contrast in the PN traction amplitude, obstruction due to the short-range image stress field scales linearly and is stronger than obstruction due to a pure long-range image stress field. (iii) For combined phase contrasts the short-range image stress field is enhanced by the elasticity contrast.

A study of the influence of the number of dislocations in a dipole pile-up followed. An increasing number of dislocations strongly decreases – as expected – the applied shear load needed for dislocation transmission. A pure elasticity contrast exhibited a smaller decrease in external transmission stress than a pure glide plane contrast. This confirms an enhanced repulsive long-range stress field for an increasing number of dislocations. For all phase contrasts the external transmission stress did not converge asymptotically towards zero as expected from the theory. Two reasons could be identified as follows. (i) The Dirichlet boundary conditions neglected the influence of slip along the glide plane. Accordingly, the obtained external transmission stress was somewhat overestimated, an effect that increases with a the number of dislocations. A larger model domain would reduce this effect. (ii) At the point of transmission, the dislocation pile-ups were not in their static equilibrium yet. For fully evolved dislocation pile-ups, the external transmission stress is expected to drop – an effect that increases with the number of dislocations.

Despite the simplifications made, the model is a suitable tool to gain insight into edge dislocation–phase boundary interaction in two-phase microstructures. Its strength is in particular that it incorporates the repulsion of dislocations by the phase boundary and the resistance against dislocation transmission exclusively through the adopted mechanics. No additional constitutive relations for dislocation–interface interactions, e.g. a transmission criterion, are required – cf. DDD where this would be necessary. On the other hand, the model can be better controlled and is computationally significantly less expensive than MD. We believe that the modelling could be extended without much difficulty to multiple glide planes, inclined glide planes, screw dislocations, etc., as long as the kinematics offered by the disregistry allows the Burgers vector content to be preserved. Such extensions will be pursued in forthcoming work. A further avenue for future work is the comparison with other models, in particular with atomistics.

Acknowledgments

We would like to thank Franz Roters and Pratheek Shanthraj of the Max Planck Institute for Iron Research for useful discussions on the present work.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research is supported by Tata Steel Europe through the Materials innovation institute (M2i) and Netherlands Organisation for Scientific Research (NWO), under the grant number Stichting voor de Technische Wetenschappen (STW) 13358 and M2i project number S22.2.1349a.

Determination of the damping coefficient

A suitable value for the damping coefficient η in Equation (13(13) $T^i\left(\mathrm{\Delta }\right)={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\dot{\mathrm{\Delta }}$(13) ) can be found through a comparison with DDD. Generally, a dislocation that travels the distance of a Burgers vector b dissipates the work(A1) $W_{\mathrm{d}\mathrm{i}\mathrm{s}}={\int }_{x_j}^{x_j+b}F\left(x\right)\mathrm{d}x$(A1) where $x_j$ and $F\left(x\right)$ denote the initial dislocation position and the force acting on the dislocation, respectively. For DDD the force $F\left(x\right)$ equates the Peach–Koehler force $F\left(x\right)={\tau }^{\mathrm{r}\mathrm{e}\mathrm{s}}b$. The resolved shear stress ${\tau }^{\mathrm{r}\mathrm{e}\mathrm{s}}$ acting on the dislocation is generally related to the dislocation velocity v via a linear drag law(A2) ${\tau }^{\mathrm{r}\mathrm{e}\mathrm{s}}=\frac{B}{b}v$(A2) with the drag coefficient B. Assuming a constant ${\tau }^{\mathrm{r}\mathrm{e}\mathrm{s}}$ over the distance b the dissipated energy equals(A3) $W_{\mathrm{d}\mathrm{i}\mathrm{s}}^{\mathrm{D}\mathrm{D}\mathrm{D}}=Bvb$(A3) Within the Peierls–Nabarro model there is no point force such as the Peach–Koehler force in DDD anymore. The force $F\left(x\right)$ is now calculated by(A4) $F\left(x\right)={\int }_{x_j^{\ast }-\mathrm{\infty }}^{x_j^{\ast }+\mathrm{\infty }}{\tau }^{\mathrm{r}\mathrm{e}\mathrm{s}}\left(x\right)\mathrm{d}x$(A4) where $x_j^{\ast }$ is the current dislocation position. Adopting the shear traction-disregistry law of Equation (13(13) $T^i\left(\mathrm{\Delta }\right)={\hat{T}}^i\sin \left(\frac{2\pi \mathrm{\Delta }}{b^i}\right)+{\eta }^i\dot{\mathrm{\Delta }}$(13) ) for ${\tau }^{\mathrm{r}\mathrm{e}\mathrm{s}}\left(x\right)$ the force $F\left(x\right)$ can be approximated for a single phase model by(A5) $F\left(x\right)={\int }_{x_j^{\ast }-\mathrm{\infty }}^{x_j^{\ast }+\mathrm{\infty }}\eta \dot{\mathrm{\Delta }}\left(x\right)\mathrm{d}x$(A5) Modifying the analytical solution (14(14) $\mathrm{\Delta }\left(x\right)=-\frac{b}{\pi }{\tan }^{-1}\left(\frac{x-x_1}{\zeta }\right)+\frac{b}{2}$(14) ) for an arbitrary dislocation position $x_j$ (A6) $\mathrm{\Delta }\left(x\right)=-\frac{b}{\pi }{\tan }^{-1}\left(\frac{x-x_j^{\ast }}{\zeta }\right)+\frac{b}{2}$(A6) with(A7) $\dot{\mathrm{\Delta }}\left(x\right)=\frac{\mathrm{d}\mathrm{\Delta }}{\mathrm{d}x}v$(A7) one obtains for a constant dislocation velocity v the force(A8) $F\left(x\right)=\eta bv$(A8) and the dissipation(A9) $W_{\mathrm{d}\mathrm{i}\mathrm{s}}^{\mathrm{P}\mathrm{N}}=\eta b^{2}v$(A9) Requiring the dissipated energy in both approaches to be equal gives(A10) $\eta =\frac{B}{b}$(A10) which is an adequate approximation for the damping coefficient.