A Computational Scheme for the Interaction between an Edge Dislocation and an Arbitrarily Shaped Inhomogeneity via the Numerical Equivalent Inclusion Method

The interactions between dislocations and inhomogeneity may play an important role in strengthening and hardening of materials. The problem can be solved analytically only for limited cases of simple geometry. By employing the recently developed numerical equivalent inclusion method, this work presents an effective computational scheme for studying the stress field due to an edge dislocation in the vicinity of an arbitrarily shaped inhomogeneity. The inhomogeneity is treated as an equivalent inclusion that is numerically discretized by rectangular elements. The mismatch between the matrix and the inhomogeneity materials are formulated through Dundurs’ parameters for numerical stability and robustness. The proposed method can efficiently and accurately evaluate the elastic field of the equivalent inclusion with the assistance of a fast Fourier transform based algorithm, constituting an essential refinement of the existing approach in the dislocation-inhomogeneity literature. Several benchmark examples are examined to demonstrate the flexibility, efficiency and accuracy of the present method.


INTRODUCTION
Barenblatts seminal works [1,2] proposed a cohesive zone model to study the fracture and failure behavior in brittle materials. Based on the hypothesis of cohesive forces, the singularity at the crack tip is circumvented and the stability of crack propagation is discussed [2]. .rom the materials point of view, crack mechanism usually involves a dislocation pile-up. Since the mobility of the dislocation can be remarkably affected by the internal force arising from the inhomogeneity, the crack propagation may be significantly hindered or promoted in the presence of material inhomogeneities, especially when the inhomogeneities are located in the vicinity of the crack tip [36].
The interaction of dislocations with the inhomogeneities is a crucial mechanism for the strengthening and hardening of materials [79]. In continuum mechanics, a crack problem may be solved by the distributed dislocation technique, which is applicable to the crackinhomogeneity interactions, where the solution of the dislocationinhomogeneity is treated as a Greens function [10,11].
Dundurs and Mura [12] first presented a closedform solution for the elastic field of an edge dislocation located outside a circular inhomogeneity, and some related work on circular inhomogeneity were subsequently conducted [13,14]. Stagni and Lizzio [15] analytically solved the interaction between an edge dislocation and an elliptical inhomogeneity. Their solution was represented by a Laurent series expansion. Hills et al. [16] has documented the available classical dislocation solutions in the literature. It is seen, however, that the analytical works on the interactions between dislocation and inhomogeneity are limited to specific and usually simple geometries.
Hutchinson [17] developed an approximate method to study the shielding effect of profuse microcracking at the tip of a macroscopic crack. In his work, the ef-fects of the microcracks relieving the residual stress are manifest on the macroscopic scale as transformation strains. The method is regarded as the zeroth-order approximation [18] in that the transformation strain is determined by the far field loading only, while the local disturbances are neglected. Hutchison approximate method was later employed by Li and his coworkers [1921] to study the interactions between crack/dislocation and inhomogeneities.
The equivalent inclusion method was originally proposed by Eshelby [22] for solving three-dimensional (3D) ellipsoidal inhomogeneity problem. In the context of micromechanics [10], the subdomain containing eigenstrain (i.e., a generic name for non-elastic strains) but having identical material is called an inclusion, to distinguish the term inhomogeneity, which refers to a subdomain with different material from its surrounding matrix. Jin et al. [23] demonstrated that equivalent inclusion method is also a handy and elegant tool for analyzing 2D inhomogeneity in plane elasticity. Zhou et al. [18] developed a numerical approach based on the equivalent inclusion method for handling arbitrarily shaped inhomogeneities of any material constants. The method, termed as numerical equivalent inclusion method, employs Dundurs parameters for describing any mismatch combinations between the inhomogeneity and the matrix materials, and has been demonstrated to be robust and stable over the entire Dundurs plane. In the implementation of numerical equivalent inclusion method, the equivalent inclusion domain is discretized into a system of elementary inclusions, each being subjected to uniform eigenstrains [10] whose magnitude is to be determined through an iterative scheme. The numerical equivalent inclusion method shows a refinement over the wellknown Hutchisons approximate method, because the latter may be viewed as neglecting a coupling term and hence the iteration revoked in determining the unknown eigenstrains.
Even for homogeneous materials, the classical dislocation solution contains a Cauchy type singularity [16], which usually is difficult to be handled by commercial finite element software. In the presence of inhomogeneities, the disturbance incurred may be viewed analogously from an inclusion with properly chosen eigenstrains, following the renowned Eshelbys equivalent inclusion method. It is therefore envisaged that the interactions between dislocations and inhomogeneities may be effectively solved via the numerical equivalent inclusion method.
The paper is organized as follows. In Sect. 2, a pictorial explanation of the equivalent inclusion method, as well as solution of an edge dislocation in a 2D homogenous medium are reviewed first. The elementary solution of a rectangular inclusion is recorded and is then utilized as building blocks for analyzing any arbitrarily shaped inclusion. .urthermore, the implementation of the numerical equivalent inclusion method is presented, where the iterative scheme for determining the equivalent eigenstrain distribution is discussed in detail. In Sect. 3, benchmark examples are utilized to validate the effectiveness of the present solution and the parametric studies on the effects of mesh size, iterative number as well as the aspect ratio are investigated. The main concluding remarks are summarized in Sect. 4, and the detailed expressions for stresses produced by a rectangular inclusion subjected to uniform eigenstrains are given in the appendix.

.ORMULATION
Consider an infinite plane containing an arbitrarily shaped inhomogeneity. An edge dislocation is located outside the inhomogeneity, which is assumed to be perfectly bonded to the matrix. As shown in .ig. 1, the infinitely extended matrix denoted by region 1 has elastic moduli , ijkl C and the subdomain (inhomogeneity) with elastic moduli * ijkl C is designated as region 2. Hereafter, the subscripts 1 and 2 are also adopted in the context to distinguish whether the quantity is associated with the matrix or inhomogeneity.
In view of the equivalent inclusion method [23], the problem in .ig. 1 may be solved by superposing a homogeneous material solution (.ig. 2a) and a corresponding inclusion solution (.ig. 2b). The former subproblem (.ig. 2a) considers an edge dislocation in the absence of the inhomogeneity, and the classical 2D solution will be recorded in Sect. 2.1. The latter sub- problem (.ig. 2b) conceives an equivalent inclusion for the disturbance due to the existence of the inhomogeneity. Since in the current work the inhomogeneity can be any arbitrary shape, the corresponding inclusion in .ig. 2b will be numerically discretized into a number of elementary inclusions, whose formulation will be discussed in Sect. 2.2. The numerical implementation of the equivalent inclusion method will be elaborated in Sect. 2.3, where an effective iterative scheme is presented for determining the distribution of the eigenstrains in the equivalent inclusion.
It has been shown in [23] that the equivalent inclusion method formulation for plane stress is apparently simpler, because the longitudinal eigenstrain component will not vanish in the case of plane strain, unless Poissons ratios of the inhomogeneity and the matrix are the same value or the inhomogeneity degenerates to a cavity. .or simplicity, only plane stress problem is investigated in the current numerical study. However, it is understood that the solution of plane stress may be readily utilized to interpret the case of plane strain, with the conversion of Kolosovs constant [24].

Edge Dislocation Solution in an Infinite Homogeneous Plane
The classical solution of an edge dislocation in an infinite homogeneous medium (.ig. 2a) has been documented in Ref. [16]. In the Cartesian coordinate system, the edge dislocation is located at 1 2 ( , ) x x ′ ′ and has Burgers vector components 1 2 , . b b The induced stresses 0 ij σ at a field point 1 2 ( , ) x x is given by [16]: 11  111  211  1  0  1  22  122  222  2  1  0  112  212  12 Note that the first subscript of the above components is associated with Burgers vector and the remaining two indices are related to the stress components. It is well-known that the above dislocation influence functions contain a singularity when the core of the dislocation is approached. The singularity may lead to numerical issues when a dislocation problem is solved by the commercial finite element software.

Elementary Solution for a Rectangular Inclusion
Since the inclusion can be any arbitrary shape and in order to take advantage of the fast .ourier transform based algorithms [25], the computational domain that encloses the inclusion is discretized into a number of uniform rectangular elements (.ig. 3). The element size is sufficiently small so that the equivalent eigenstrain is assumed to be uniform inside each rectangular patch. A typical rectangular element having dimensions 1 2 ∆ × ∆ and centered at 10 20 ( , ) x x ′ ′ is shown in .ig. 3. The rectangular inclusion, whose sides are parallel to the coordinate axes, is subjected to uniform eigenstrains * . ij ε The stresses inside and outside the rectangular inclusion may be represented in a unified expression [26]: .ig. 2. Schematic of the equivalent inclusion method. An unperturbed field due to a dislocation (a), a disturbance field caused by the corresponding inclusion (b).
where ijkl T is the elementary influence coefficient and are detailed in the Appendix.
With the superposition principle, the resultant elastic field due to the original arbitrarily shaped inclusion (cf. left of .ig. 3) is evaluated by summing contributions from each element. It should be emphasized that a direct superposition may lead to a huge computational burden, especially for a fine mesh. However, the summation has a convolution nature when the computational domain is discretized into uniform rectangular meshes, and hence can be tremendously accelerated by taking advantage of the fast .ourier transform related algorithms [25].

Implementation of the Equivalent Inclusion Method
The fundamentals of equivalent inclusion method have been detailed in [23]. The essence of equivalent inclusion method is established upon the consistency condition that the elastic field disturbance due to the presence of inhomogeneity can be simulated by the resulting field generated by a corresponding inclusion when the eigenstrain is chosen properly. That is, according to Eshelby [22], the consistency condition is enforced inside the inclusion: ijkl kl klmn mn ijkl kl klmn mn kl The details have been discussed in [23], and the counterpart consistency condition is derived as where kl σ is the eigenstresses produced by the unknown eigenstrains (cf. .ig. 2b). As discussed in Sect. 2.2, the eigenstresses may be solved numerically for a given distribution of eigenstrains (.ig. 3). By defining ijkl M ∆ as the compliance moduli difference between the inhomogeneity and the matrix material, i.e.
The above Eq. (6) is applicable for a general 3D configuration. In the present numerical implementation of the plane stress problem, only the in-plane components are retained, and the compliance moduli difference may be presented in matrix form as µ κ + −µ κ + α = µ κ + +µ κ + µ κ − −µ κ − β = µ κ + +µ κ + It is also noted that in Hutchinsons approximation, the disturbance term (eigenstresses) kl σ at the right hand side of is neglected, hence the equivalent eigenstrain components * ij ε are directly estimated by 0 .
ijkl kl M ∆ ε Since no iteration is involved, the approximation is called zeroth order. However, such simple estimation may result in relative rough approximation.
The present computational scheme may start from Hutchinsons approximation, but with succeeding refinement through iteration, as shown in the flow chart (.ig. 4). .or ease of presentation, the superscript (K) .ig. 3. Schematic of the discretization for an arbitrarily shaped inhomogeneity with uniform rectangular elements (color online). in the flow chart denotes the Kth iteration step. Before the start of iteration, the eigenstresses are assumed to be zero. On each step of iteration, Eq. (6) is used to update the current estimation of the equivalent eigenstrains. Then, the current relative error is estimated. Provided that either the result meets the given accuracy constraints or the maximum allowed number of iteration has achieved, the iteration process is stopped. Otherwise, another iteration is performed. In this case, the eigenstresses are re-evaluated using the current eigenstrains (cf. .ig. 3), and then Eq. (6) is recalled to update the equivalent eigenstrains.

RESULTS AND DISCUSSIONS
The present computational scheme is programmed in standard .ORTRAN language and implemented in double precision. The computations are performed on a personal computer with a 3.4 GHz i7 CPU and 16.0 GB memory.
In order to illustrate the application of the aforementioned method, an edge dislocation interacting with an elliptical inhomogeneity in an infinite elastic matrix is considered (.ig. 5). The edge dislocation is located at 1 (2 , 0) a and has Burgers vectors 1 2 ( , ). b b The inhomogeneity is centered at the origin of Cartesian coordinate systems and has semi-axes 1 a and 2 . a The main components of the computational parameters are listed in table.
The closed-form solution obtained by Dundurs and Mura [12], corresponding to a circular inhomogeneity (aspect ratio 2 1 a a = 1), is utilized to validate the present solution. In the benchmark study, the mesh is chosen as 128 128 × and the maximum allowed number of iteration is set to be 5. Comparisons of the normalized shear stress along the axis 1 x are shown in .ig. 6. It is demonstrated that the present numerical equivalent inclusion method solution exhibits excellent agreement with the exact solution.
In order to better illustrate the disturbance due to the presence of the inhomogeneity, the singular Cauchy term is then removed, and the normalized shear stress component along the axis 1 x is compared with the zeroth-order results by Hutchinson approximation [17]. As illustrated in .ig. 7, the present numerical equivalent inclusion method solution consistently shows better agreement with DundursMura solution [12]. More .ig. 4. .low chart of the numerical equivalent inclusion method.
.ig. 5. An edge dislocation interacting with an elliptical inhomogeneity in an infinite elastic medium. related details are also reported later in .ig. 8, where 0th approximation refers to the results by Hutchinsons estimation.
One of the key issues in numerical equivalent inclusion method is the discretization error along the geometrical boundaries, since uniform rectangular mesh has to be used in order to take advantage of fast .ourier transform acceleration. To address this issue, results of different mesh schemes are investigated. As illustrated in .ig. 9, typical results for meshes of 32 × 32 and 512 × 512 are compared with the previous one using 128 × 128. The present computation shows that a finer discretization usually yields better accuracy. However, even the coarsest mesh of 32 × 32 shows satisfactory agreement with the exact solution.
Numerical tests are also performed extensive to illustrate the CPU time savings by virtue of the fast .ourier transform techniques. x y x y O N N N N multiplication operations. Therefore, when the 2D grids are refined by a factor of 2 in both directions, the time consumed by the direct superposition is expected to increase by a factor of 16, as opposed to 4 when fast .ourier transform algorithms are adopted. The above CPU time estimation roughly agrees with our numerical executions (.ig. 10). Note that for a large number of elements, the memory access may also consume considerable time, which could be the cause of the small amount deviations from the CPU time estimation. σ are the number of the field points, the exact value and the approximate numerical solution, respectively. As demonstrated in .ig. 8, the average error for zeroth order approximation may exceed 20%, and would be significantly alleviated after the first iteration. After the 5th iteration, the average relative errors in most cases drop below 5%, sufficient for engineering applications. .or a fixed iteration number, the numerical error drops rapidly from the mesh of 32 × 32 to 128 × 128, and then varies moderately as the discretization becomes even finer.
The effects of the aspect ratio 2 1 a a for the normalized shear stress along the axis 1 x is also investigated in .ig. 11. .or the cases of 0.1, 0.5 and 1.0, the dimensionless results decline gradually and then tend toward a stable value as 1 1 x a ranges from 1.5 to 2.5, while the normalized disturbance shear stresses for the other cases would rise to a peak and then decrease slightly. It can be seen that the shape effects of the inhomogeneity may have a significant effect on the interaction.

CONCLUDING REMARKS
In reality, the inhomogeneities in the engineering materials are not only of arbitrary shape but also distribute in a random way. The interaction between edge dislocations and inhomogeneities may have a significant influence on the crack growth in heterogeneous materials and strain hardening in metal alloys. However, analytical solutions only exist for limited cases. The involved singularity nature may even cause numerical obstacles when solved by the commercial finite element software. The recently developed numerical equivalent inclusion method may a choice for quantitative analysis of the dislocationinhomogeneity interactions. Satisfactory accuracy is achieved by an iterative scheme when compared with Hutchinsons approximate method.
In this work, the fast .ourier transform based algorithms can be seamlessly incorporated with the numerical equivalent inclusion method. Dundurs parameters are employed in the present equivalent inclusion method formulation to enhance the numerical stability and robustness. The effectiveness of the iterative method is illustrated by comparison with the analytical solution. The parametric studies demonstrate that the present method may considerably improve the computational efficiency and accuracy. Although only the circular or elliptical inhomogeneity is exemplified in this work, the discrete rectangular elements can be treated as building blocks, and therefore it is anticipated that the present numerical equivalent inclusion method would be flexible and versatile for analyzing inhomogeneities of arbitrary shape.