Two-dimensional nonlocal Eshelby’s inclusion theory: eigenstress-driven formulation and applications

The classical Eshelby’s theory, developed based on local linear elasticity, cannot be applied to inclusion problems that involve nonlocal (long range) elastic effects often observed in micromechanical systems. In this study, we introduce the extension of Eshelby’s inclusion theory to nonlocal elasticity. Starting from Eringen’s integral formulation of nonlocal elasticity, an eigenstress-driven nonlocal Eshelby’s inclusion theory is presented. The eigenstress-driven approach is shown to be a valid mathematical extension of the classical eigenstrain-driven approach in the context of nonlocal inclusion problems. Two individual numerical approaches are developed and applied to simulate inclusion problems and numerically extract the corresponding nonlocal Eshelby tensor. The numerical results obtained from both approaches confirm the validity of the derived nonlocal Eshelby tensor and its ability to capture the non-uniform eigenstress distribution within an elliptic inclusion. These results also help reveal the fundamental difference between the mechanical behaviour of the classical local and the nonlocal inclusion problems. The eigenstress-driven nonlocal inclusion theory could provide the necessary theoretical foundation for the development of homogenization methods of nonlocal heterogeneous media.

The classical Eshelby's theory, developed based on local linear elasticity, cannot be applied to inclusion problems that involve nonlocal (long range) elastic effects often observed in micromechanical systems.In this study, we introduce the extension of Eshelby's inclusion theory to nonlocal elasticity.Starting from Eringen's integral formulation of nonlocal elasticity, an eigenstress-driven nonlocal Eshelby's inclusion theory is presented.The eigenstress-driven approach is shown to be a valid mathematical extension of the classical eigenstrain-driven approach in the context of nonlocal inclusion problems.Two individual numerical approaches are developed and applied to simulate inclusion problems and numerically extract the corresponding nonlocal Eshelby tensor.The numerical results obtained from both approaches confirm the validity of the derived nonlocal Eshelby tensor and its ability to capture the non-uniform eigenstress distribution within an elliptic inclusion.These results also help reveal the fundamental difference between the mechanical behaviour of the classical local and the nonlocal inclusion problems.The eigenstress-driven nonlocal inclusion theory could provide the necessary theoretical foundation for the development of homogenization methods of nonlocal heterogeneous media.

Introduction
Eshelby's inclusion theory has long served as a foundational concept in classical linear elasticity of heterogeneous materials.By employing the well-known four-step thought experiment, Eshelby obtained a closed-form expression to quantify the stress and strain fields associated with an elastic inclusion in an elastic matrix [1]; the quantity enabling the quantification of these fields is often known as the Eshelby tensor.Thanks to its theoretical versatility and physical relevance, Eshelby's solution to the embedded inclusion problem has been a cornerstone of solid mechanics since the 1950s; later, this theory was applied to the modelling of, to name a few, composite materials [2,3] as well as material defects and interfaces [1,4,5], and to the homogenization of heterogeneous micromechanical structures [6][7][8].
Despite the extensive body of literature on the original Eshelby's inclusion problem and its micro-mechanical applications, there are still material configurations showing non-negligible discrepancies between the theoretical prediction and the experimental observation [9][10][11].Recent theoretical and experimental investigations have revealed the existence of non-negligible sizedependent effects (also known as nonlocal effects) in micro-mechanical systems.The occurrence of this effect can be related to interatomic forces [12], micropolar effects [13,14], couple stress [15,16], interface stress [4,17] and microcrack interactions [18].The origin of size-dependent effects in these materials can be attributed to the intricate microstructure that can span multiple spatial scales.The presence of these microstructural features introduces localized material heterogeneities, which in turn alter the strain and stress distribution, create strain concentration and cause perturbations of the mechanical response.The interplay between these microstructural features and the associated physical behaviour occurs at various scales and eventually manifest as nonlocal effects at the continuum level [19,20].
The presence of nonlocal effects within the microstructure poses substantial challenges to the applicability of Eshelby's inclusion theory, as well as those micromechanical analysis methods, for example, mean-field theory and homogenization methods [2,8], that depend on the classical Eshelby's approach.Indeed, all these theories are based on the local-reaction principle [21], therefore they are unable to capture any nonlocal behaviour [12,22,23].These local theories, which assume that the response of a material at a point depends only on the state of deformation in an infinitesimally small neighbourhood around that point, are inadequate to account for the presence of long-range forces [18,19].In that regard, classical local linear elasticity, which serves as the theoretical postulate of Eshelby's inclusion problem [24], also fails to capture nonlocal effects and therefore limits the potential application of Eshelby's method to the more general class of nonlocal micromechanical systems.
The limitations of classical local elasticity theories have motivated the development of non-classical continuum theories that provide a promising theoretical foundation to address the incompatibility between Eshelby's theory and the aforementioned nonlocal effects.An examination of existing nonlocal formulations in the field of elasticity highlights two main classes of nonlocal continuum theories: the gradient elasticity theory [25,26] and the integral nonlocal elasticity [27,28].In gradient elasticity theory, the enrichment of nonlocal effects is achieved by incorporating high-order gradient terms [29,30].From a mathematical perspective, the high-order gradient enhancement can be seen as a higher-order Taylor series approximation to the nonlocal fields such that it enables a more accurate characterization of nonlocal perturbations than classical local elasticity theory.Nevertheless, the gradient formulation is still localized due to the local action of differential operators (recall that differential operators are defined by taking the limit with respect to an infinitesimal value, which eventually leads to a localized behaviour), and therefore often lacks a clear physical interpretation of specific microstructural features or material behaviour [31].
Extensive theoretical investigations and a wide range of practical applications have further confirmed the unique advantages offered by integral elasticity in modelling nonlocal effects, when compared with gradient elasticity.Unlike the gradient approaches, the enrichment of nonlocal effects in integral nonlocal theories is achieved by incorporating weighted convolutional integral operators [27].From a mathematical standpoint, the intrinsic nonlocal characteristic of integral operators facilitates the effective inclusion and quantification of nonlocal effects.The weight function (or the nonlocal kernel) measures the strength of the nonlocal interactions between two distant points, while the support of the integral operator defines the extent of the long-range actions and it is commonly referred to as the nonlocal horizon [32].This coherent interpretation of integral formulations allows a better understanding of integral over gradient theories [19,33] and has been applied to solve numerous practical problems that exhibit nonlocal effects [34,35].
To date, only a handful of studies have been conducted to expand Eshelby's theory to the realm of nonlocal elasticity.For gradient nonlocal formulations, Reid et al. [36] first investigated the twodimensional Eshelby's inclusion problem with strain-gradient nonlocal contributions and pointed out that the nonlocal length scale is a good indicator of the extent to which the inclusion/matrix interface is modified by gradient effects.Later, Zhang et al. [37] studied size-dependent inclusion problems by considering a modified strain gradient elasticity with couple stresses.Using the same gradient approach, Maranganti et al. [38] investigated an inclusion problem with electromechanical coupling effects at the nanoscale.Based on these early studies, Gao [9,10,39] has systematically explored the strain gradient solution for a range of Eshelby's inclusion problems.More recent work on the strain gradient reformulation of Eshelby's inclusion problems can be found in [40][41][42].By contrast, investigations into integral-based nonlocal Eshelby's inclusion problems are considerably less developed.To the best of the author's knowledge, only one study [43] has attempted to solve Eshelby's inclusion problem with integral nonlocal formulations.
This study aims at bridging this theoretical gap by extending the classical Eshelby's inclusion theory to nonlocal problems through the incorporation of Eringen's integral nonlocal elasticity [27].Theoretical challenges arise when attempting this task.In the original Eshelby's problem, the constant eigenstrain assumption was adopted to produce an equivalent constant eigenstress field [1]; this assumption allowed the analytical derivation of the Eshelby tensor in local inclusion problems, without introducing any spurious body force (since the divergence of constant eigenstress is always zero).This approach cannot be directly extended to nonlocal cases when considering Eringen's theory, as the stress in Eringen's theory is defined via an integral formulation such that the eigenstress field within an ellipsoidal inclusion becomes fully nonlocal and fail to maintain a constant property due to the interaction between the nonlocal horizon and the inclusion boundary.To reconcile this critical contradiction between Eshelby's inclusion problem and Eringen's nonlocal elasticity, in this study, we propose a different approach by introducing a constant eigenstress assumption (in contrast to the constant eigenstrain assumption initially used by Eshelby).By adopting the proposed assumption, we are able to maintain a uniform nonlocal eigenstress field, which in turn facilitates the analytical derivation of the field tensor representing the stress-strain relationships within the inclusion.This field tensor serves a similar purpose as the classical Eshelby tensor, describing the mechanical behaviour of the inclusion.To underscore this approach, we introduce the terminology eigenstress-driven approach.Correspondingly, the resulting field tensor is renamed as the eigenstress-driven Eshelby tensor.This nomenclature highlights some key differences of the proposed approach and the key role played by the eigenstress in determining the mechanical behaviour of the inclusion within the framework of nonlocal elasticity.
The overall objectives and contributions of this study can be summarized as follows: -Definition of eigenstress-driven nonlocal Eshelby's inclusion problem.We employ Eringen's integral nonlocal elasticity approach theory along with the eigenstress-driven approach to generalize the classical Eshelby's inclusion problem (and the well-known four-step thought experiment) to the nonlocal case.The eigenstress-driven Eshelby's inclusion problem is shown to be well defined and is equivalent to the classical eigenstrain-driven definition when considering the limit case of local linear elasticity.-Solution procedure and analytical derivation of eigenstress-driven nonlocal Eshelby tensor.We derive the closed-form expression of the eigenstress-driven nonlocal Eshelby tensor by leveraging the extended four-step thought experiment.The nonlocal Green's function and the corresponding Fourier transform are presented.We note that this study appears to be the first attempt at employing the Fourier transform approach to solve the integral nonlocal elasticity equation and to derive the nonlocal Eshelby tensor.-Numerical solution and validation.We develop two individual numerical approaches to assess and validate our theoretical findings.The two approaches serve different purposes.The first approach, indicated as direct numerical computation (DNC), is based on numerical integration techniques in order to numerically evaluate the nonlocal Eshelby tensor.The second approach is built upon a multiscale multimesh finite-element method (M 2 -FEM) [44].This approach is specifically tailored to simulate the nonlocal Eshelby's inclusion problem by leveraging the exceptional capabilities of the M 2 -FEM in solving integral nonlocal elasticity problems.
The remainder of this paper is structured as follows.In §2, we briefly revisit Eringen's integral nonlocal elasticity theory as well as the classical Eshelby's inclusion problem.Following the preliminary introduction, in §3, we derive Green's function for a nonlocal elastic continuum by using the Fourier transform technique.By incorporating Greens function approach and the extended four-step thought experiment, in §4, we present the rigorous definition of the eigenstress-driven nonlocal Eshelby's inclusion problem and derive the eigenstress-driven nonlocal Eshelby tensor.Specific discussion is presented to justify the rationale of the eigenstressdriven approach and its equivalence to the classical eigenstrain-driven approach.Numerical approaches and simulation results are presented in § §5 and 6, respectively, to numerically assess and validate the nonlocal Eshelby tensor.Finally, some concluding remarks are presented in §7.

Preliminaries
In this section, we briefly review the theoretical foundation and physical background that form the basis of this study.We begin by reviewing the mathematical formulation of Eringen's integral nonlocal elasticity theory.Building upon the nonlocal elastic formulation, we then revisit the concepts of eigenstrain and eigenstress, originally defined in the context of local elasticity, and extend their definitions and physical implications to the nonlocal case.By employing the proposed nonlocal elastic framework, we establish the nonlocal Eshelby's inclusion problem and outline the associated four-step thought experiment.Before proceeding to the theories, it is necessary to clarify that in this study we focus on solving two-dimensional nonlocal Eshelby's inclusion problem only.This choice is primarily due to the intricate interplay between integral nonlocal elasticity and Eshelby's problem, which makes the numerical validation unapproachable as the problem dimension increases.Meanwhile, it should also be noted that while our analytical derivations are based on two-dimensional problems, they can be readily extended to a threedimensional case [24].

(a) Integral nonlocal elasticity
The origins of strain-driven nonlocal approaches can be attributed to Eringen's nonlocal elasticity theory, where the concept of nonlocality is introduced into continuum mechanics by considering an integral form of constitutive relations [27].According to Eringen's integral approach to nonlocal elasticity [45], the stress-strain constitutive relations at any material point in a threedimensional homogeneous nonlocal solid Ω can be expressed using integral operators in the following form: where σ , C and are the nonlocal stress, the stiffness tensor and the local strain tensor, respectively.Further, x denotes a specific material point on the nonlocal solid Ω; and x ∈ H(x) is a dummy spatial variable denoting the (remaining) material points that interact with x via longrange connections and are contained within a characteristic nonlocal horizon H.The strength of the nonlocal effects is captured by the nonlocal kernel function denoted as K(x, x ) with support on the nonlocal horizon H.Note that, in this study focused on two-dimensional elasticity, both Ω and H are two-dimensional Euclidean spaces.Moreover, due to physical constraints, the horizon should never exceed the material domain, that is H(x) ⊆ Ω is always a subspace of Ω, ∀ x ∈ Ω.The nonlocal definition of the constitutive relations presented in equation (2.1) enables a clear interpretation of nonlocality: the stress σ (x) at point x is due to the combined effect of the strain (x ) at points x contained within the nonlocal horizon H around x.The strength of the point-topoint nonlocal interaction between x and x is then determined by the nonlocal kernel function K(x, x ).
The strong form of the governing equations can be derived by variational minimization of the potential energy functional [46] obtained from the constitutive framework presented in equation (2.1) 2) The governing equations are subject to the following boundary and initial conditions: where u and u 0 denote the displacement field and the prescribed displacement on the boundary, respectively,

(b) Eshelby's inclusion problem and four-step thought experiment
Eshelby's inclusion problem focuses on analysing the distribution of the elastic strain field that arises within a three-dimensional ellipsoidal (or elliptical, in the case of two-dimensional problems) inclusion embedded in an unbounded matrix (where both the inclusion and the matrix are assumed to be homogeneous and have the same material properties).Specifically, the objective in Eshelby's inclusion problem is to determine the distribution of elastic strain that arises from the interaction between the inclusion and its surrounding matrix, under the condition that the inclusion is subjected to a prescribed uniform eigenstrain * . 1shelby brilliantly addressed the problem in an analytical way by employing the renowned four-step thought experiment described below: -Step 1: Remove the inclusion from the matrix.When outside the matrix, the inclusion is subjected to free boundary conditions, and its deformation is governed by the original uniform eigenstrain * .Note that at this step the stress field within the inclusion is zero since there is no body force and boundary force applied.-Step 2: Apply the surface traction t to the inclusion in order to force the inclusion to deform back to its original shape.The surface traction creates a negative strain field that cancels out the initial eigenstrain * , and generates an internal stress field that is the opposite of the eigenstress σ * , which is a fictitious stress produced by the eigenstrain [1,24] where C is the material stiffness tensor.Note that the eigenstress in equation (2.4) is defined through the classical local constitutive relations, as we only introduce the classical Eshelby's inclusion problem at this point.A detailed discussion of the eigenstress and its formulation will be presented in the next section.-Step 3: Put the inclusion back into the matrix.At this step, there is no change in the elastic responses within both the inclusion and the matrix.-Step 4: Remove the surface traction t.Following the removal of the external traction force, the inclusion tends to deform and return back to the initial state it had in Step 1.However, due to the confinement imposed by the matrix, the inclusion cannot deform freely and can only interact with the matrix, eventually resulting in the generation of a constrained elastic strain and stress field.
The four-step thought experiment is schematically illustrated in figure 1.Note that, as clarified in the introduction, the local Eshebly's inclusion problem is concerned with a uniform eigenstrain *2 (not a uniform eigenstress σ * ).In the subsequent sections, we will discuss in detail how this eigenstrain-driven approach is not applicable when considering the nonlocal extension of Eshelby's inclusion problem and, as a result, introduce an eigenstate-driven approach.

Green's function and Fourier transform approach
Prior to introducing the specific problem configuration, we first introduce the mathematical tools that are necessary to solve the nonlocal Eshelby's inclusion problem.Specifically, we extend the traditional Green's function and Fourier transform approach to be applicable to a nonlocal scenario.It should be highlighted that, due to the presence of the non-vanishing nonlocal horizon, the derivation of Green's function to the partial diff-integral governing equations and the associated Fourier transform approach becomes complicated and differs from the derivation in local situations.Given the importance and complexity of these mathematical tools, in this section, we decide to present the full derivation before accessing the nonlocal Eshleby's inclusion problem.

(a) Static Green's function in integral nonlocal linear elasticity
Similarly to the case of local elasticity, the static Green's function in nonlocal elasticity is defined by its physical meaning where u i (x, y, b j ) indicates the ith displacement component seen at x due to a point load b j in the jth direction applied at y (note that both x ∈ R 2 and y ∈ R 2 are within an unbounded twodimensional domain).Substituting the above definition into the nonlocal constitutive relation in equation (2.1), the resulting nonlocal stress is where the subscript ,l x is used to denote the partial derivative with respect to the variable x in the direction l.In analogy with the classical local case, the assumption of isotropic and homogeneous material is also applied in order to derive the nonlocal Green's function.Further substituting into equation (2.2) and ignoring the inertial term, we obtain the governing equation in terms of the static Green's function for the integral nonlocal case where δ(x − y) denotes the Dirac delta function and δmi denotes the second-order identity tensor.
It is important to note that, in contrast to the governing equations for a local system where only partial differential operators are involved [48], the governing equations for a nonlocal system incorporate a integro-differential formulation.For the sake of brevity, in the following, we use tensor T = T kmlj (x, y) to denote the integro-differential formulation of the nonlocal Green's function.
The tensor T kmlj (x, y) can be simplified by applying the high-dimensional Leibniz-Reynolds transport theorem [49] T kmlj (x, y) = where n r (x ) is the normal vector to the surface of the nonlocal horizon H(x) (which is represented by ∂H(x), see the surface integral term), and is the spatial velocity of points (in analogy with the temporal velocity in classical Leibniz-Reynolds transport theorem) ∂H(x), which describes how points on the boundary of the nonlocal horizon H(x) move as x changes its position.The symbol • is used to differentiate the point and the corresponding spatial velocity defined on ∂H(x) from the original point x.Note that equation (3.5) is valid since the nonlocal horizon is set to be homogeneous at any point throughout the whole material domain (recall that we have made the assumption of isotropic and homogeneous material).In other words, the shape of the nonlocal horizon should remain unchanged and can be defined in the following form: where H 0 is the shape of nonlocal horizon when x = 0. We highlight that the variation of x in equation (3.6)only shifts the nonlocal horizon without changing its shape and therefore guarantees translation invariance.Further leveraging the divergence theorem, we have such that the original surface integral form of T  (x, y) the partial derivative is calculated with respect to x.To unify the partial derivative in both terms, we leverage the symmetry property3 of the kernel function K(x, x ) [46] K(x, x ) = K(x , x), (3.8) which implies that the kernel function only depends on the relative distance between x and x .By defining a new variable τ = x − x, the derivative of K(x, x ) with respect to x and x can be recast into the following form: and similarly Substituting equations (3.9) and (3.10) back into equation (3.4), we obtain Notably, unlike the original expression in equation (3.4), the updated expression of T kmlj (x, y) shown above avoids the derivative of the nonlocal integral by switching the order between the partial derivative in the jth direction and the nonlocal integration.In addition, by using the change of variables the domain of nonlocal integration is now transformed into a fixed domain H 0 ; the partial derivative in G km,l τ j τ now is taken with respect to the new variable τ , not x .Using this updated expression and putting it back into equation (3.3), we obtain the simplified governing equations in terms of the nonlocal Green's function between the two points x and y when considering an isotropic, homogeneous and infinitely large nonlocal solid.The same conclusion can be easily shown to be applicable to T kmlj since where the nonlocal integral containing G km,l τ j τ (y − x − τ ) should depend on y − x only (this is obvious since τ is the variable of integration and it vanishes after calculating the integral).
Defining the new variable x = y − x and substituting it back into equation (3.4), T kmlj can be further recast into the following form: We immediately note that the new expression of T kmlj is in fact the convolution of two functions where K(•) is called the truncated kernel function defined by This unique property implies the possibility of using the convolution theorem and the Fourier transform to solve for the nonlocal Green's function.Substituting the convolution form of T kmlj (x) back into equation (3.12) and applying the Fourier transform on both sides of the equation, we obtain note that δ(x − y) vanishes after applying Fourier transform since the Fourier transform of the Dirac delta function is unity.To simplify the algebraic expression, in the following, we do not distinguish different differentiation variables and ignore the superscripts (such as τ and x) over the partial derivative index (for example, see G km,lj in equation (3.14)).
The solution to equation (3.17) requires obtaining the Fourier transform of G km,lj (x), which can be done by leveraging the differentiation property of the Fourier transform.Suppose the Fourier transform of the nonlocal Green's function and the corresponding inverse Fourier transform are defined by where g km (ξ ) is the Fourier transform of the nonlocal Green's function and ξ is the variable in Fourier space, and the coefficient 1/(2π ) 2 is employed in the inverse Fourier transform since we are considering two-dimensional problems.Using the above definition, we obtain which leads to Substituting it back into equation (3.17) and simplifying, we obtain the final governing equation of the nonlocal Green's function in Fourier space as where is the Fourier transform of the truncated kernel function K(x), which is the only term that contains nonlocal factors.This is an important finding because it allows us to mathematically separate and isolate the influence of nonlocality exclusively within k(ξ ).As expected, in the limit to local behaviour one should obtain K(x) = δ(x) and k(ξ ) = 1, such that equation (3.21) becomes the governing equation of the static Green's function in classical local linear elasticity [24].To obtain g km (ξ ) in equation (3.21), we first define the following second-order tensor: where |ξ | is the length of the ξ and ξl = ξ l /|ξ | is the direction vector of ξ (with unit length).Substituting this definition back into equation (3.21) and taking the inverse of w ik on both sides of the equation, we obtain an alternative form of g km (ξ ) as follows: where w −1 km is the inverse of w ik such that they satisfy w ik w −1 km = δim .Further substituting equation (3.24) back into equation (3.19), we obtain In many practical cases, k(ξ ) is often assumed to be independent of the orientation and only a function of the magnitude of ξ .Hence, equation (3.25) can be further written as In fact, this assumption can be easily satisfied by considering a radially symmetric kernel function K(x); for example, a two-dimensional Gaussian function with truncated circular support H 0 .In the subsequent section, we will use the nonlocal Green's function to solve the nonlocal inclusion problem and derive the nonlocal Eshelby tensor.

Derivation of eigenstress-driven nonlocal Eshelby tensor
In this section, we incorporate the four-step thought experiment to formulate the nonlocal inclusion problem.To account for the presence of nonlocal effects, we propose the eigenstressdriven definition of the nonlocal Eshelby tensor.This definition allows extending the classical Eshelby's findings and solution procedure from the local case to the nonlocal case.The eigenstress-driven nonlocal Eshelby's problem is then solved by incorporating the nonlocal Green's function obtained in §3.The eigenstress-driven approach addresses the challenges posed by nonlocality and is proven to be equivalent to the traditional eigenstrain-driven approach in local problems, hence providing a more comprehensive and universal framework for studying nonlocal inclusion problems.

(a) Eigenstress-driven nonlocal inclusion problem
Before discussing the eigenstress-driven nonlocal inclusion problem, we assess the limitations associated with the traditional eigenstrain-driven approach.Following the limitation analysis, we then propose the eigenstress-driven approach as an alternative and utilize it to extend the classical Eshelby's problem to the nonlocal case.
Schematics illustrating the possible interaction of the nonlocal horizon with the boundary of the inclusion.The inclusion V 0 consists of (1) the reduced inclusion V − , and (2) its complement V c − .The nonlocal horizon for points within V − remains intact and always keeps the shape H 0 , while the nonlocal horizon for points within V c − is truncated by the inclusion's boundary.The shape of the reduced inclusion is determined by both the inclusion V 0 and the nonlocal horizon H 0 .V − converges to V 0 as H 0 → 0.

(i) Limitation of eigenstrain-driven definition
Consider an eigenstrain field prescribed inside the inclusion V In the classical local Eshelby's inclusion problem, the eigenstrain field within the inclusion is set to be constant in order to create a constant eigenstress field [24].However, this technique is no longer valid when considering integral nonlocal elasticity because the nonlocal horizon H (or, equivalently, the area of integration in the definition of nonlocal stress, see equation (2.1)) can be physically truncated by the inclusion boundary (see figure 2 for details).Therefore, it is necessary to define a non-uniform eigenstrain distribution to ensure a constant eigenstress field in nonlocal problems.
Based on the above discussion, we define a new eigenstrain field through the following equation: such that now σ * is the constant nonlocal eigenstress defined within the inclusion and H − (x) = H(x) ∩ V 0 is the truncated nonlocal horizon, which determines the eigenstrain distribution.The inclusion V 0 can be divided into two sub-regions depending on the position of x where V − is called the reduced inclusion and V c − is the complement of V − .Note that, unlike the definition of nonlocal stress in a homogeneous medium (see equation (2.1)), the nonlocal eigenstress must be defined by using the nonlocal integral operator with truncated support H − (x).This is because the eigenstress should be determined by eigenstrain only, which is defined exclusively within the inclusion V 0 (not the whole domain).
Figure 2 shows the detailed geometric configuration of the nonlocal Eshelby's inclusion problem.It can be found that the nonlocal horizon H − (x) = H(x) holds when x ∈ V − .This implies that the nonlocal horizon at points within the reduced inclusion will not be truncated and always retain its original shape.Correspondingly, the nonlocal horizon of points within V c − is always truncated by the inclusion boundary and therefore it is mathematically impossible to obtain a uniform solution * (x) that satisfies equation ( 4 Eshelby's problem.Additionally, the classical eigenstrain-based definition of the local Eshelby tensor is not applicable and will fail to reproduce the same findings (such as the constant constrained stress field in the elliptic inclusion) in nonlocal problems.To address this issue, we perform the four-step experiment to derive a nonlocal Eshelby tensor based on the constant eigenstress assumption (or for short, the eigenstress-driven nonlocal Eshelby tensor).

(ii) Eigenstress-driven four-step experiment
We extend the original four-step thought experiment introduced in §2 to a nonlocal problem (where both the matrix and the inclusion are nonlocal) in order to derive the eigenstress-driven nonlocal Eshelby tensor.Specifically, we take the constant eigenstress assumption to define the nonlocal inclusion problem and follow the four-step procedure (see figure 1 for reference) to obtain the constrained nonlocal stress within the inclusion.
Step 1: Remove the inclusion from the matrix.The strain and stress state in both the matrix and the inclusion can be obtained by M = 0, I = * and σ M = 0, σ I = 0 Step 2: Apply the surface traction t to S 0 to force the inclusion return to its original shape M = 0, I = * + el = 0 Steps 3 and 4: Put the inclusion back into the matrix and remove the traction force t.Mechanical interactions arise at the surface between the inclusion and the matrix, and therefore create the constrained strain and the constrained nonlocal stress field (denoted by c and σ c , respectively) in both domains The last two steps are the most critical to derive the nonlocal Eshelby tensor and differ the most compared with the derivation in the local Eshelby's inclusion problem [24].Note that, when the inclusion is considered outside the matrix (i.e.step 1 and 2), the nonlocal horizon needs to be truncated as described previously (figure 2).However, when the inclusion is put back into the matrix, the truncated part of the nonlocal horizon is automatically refilled by the same homogeneous nonlocal material from the matrix.In this regard, the constrained nonlocal stress σ c should be computed using H(x) rather than the truncated H − (x) such that the stress inside the inclusion which recovers the additivity property ( l S is the local Eshelby tensor and I is the identity tensor).
The lack of additivity shown in equation (4.9) raises questions regarding the definition of the nonlocal Eshelby tensor using the constant eigenstrain assumption.As indicated by equation (4.10), the eigenstrain-driven definition of Eshelby tensor S l is effective in local problems to describe the linear-like relationship between the constrained strain c and eigenstrain * .However, it is impractical to reproduce the same linear-like relationship by using the same eigenstrain-driven definition of Eshelby tensor, as the additivity is not preserved in nonlocal problems.
(iii) Remarks on the use of constant eigenstress The constant eigenstress assumption deserves some additional remarks.Broadly speaking, the eigenstress can be physically interpreted as the stress that is required to produce the same deformation described by the eigenstrain.Alternatively, the eigenstress can also be seen as the residual stress that exists when the material deformation is zero.In this regard, the constant eigenstrain and constant eigenstress assumption are equivalent in classical local Eshelby's inclusion problems (as a constant eigenstrain always induces a constant eigenstress, and vice versa).However, this equivalency breaks up immediately after introducing the integral constitutive relations shown in equation (2.1).As previously discussed, in the nonlocal Eshelby's inclusion problem, a constant eigenstrain induces a non-uniform eigenstress field within the inclusion, which introduces a spurious body force since the divergence of non-uniform eigenstress is nonzero.The presence of spurious body forces implies that, in Step 2 of the virtual experiment, it is physically unfeasible to generate a mechanical strain through boundary tractions that counterbalances the prescribed constant eigenstrain.In fact, researchers [50,51] have pointed out that the constant nonlocal eigenstress can be easily generated, for instance, by considering a uniform external thermal load within the inclusion.Therefore, it is recommended to adopt the constant eigenstress assumption in nonlocal Eshelby's inclusion problem, rather than the constant eigenstrain assumption as in classical Eshelby's theory.

(b) Eigenstress-driven nonlocal Eshelby tensor (i) Definition of eigenstress-driven nonlocal Eshelby tensor
We define the nonlocal Eshelby tensor nl S σ based on the eigenstress-driven approach discussed above.Note that to standardize the notation, in the remaining section of this study, we denote nonlocal and local quantities with left superscripts nl and l , respectively.Additionally, we use right superscripts σ and to represent variables related to the eigenstress-driven and eigenstrain-driven approaches, respectively.From equation (4.7) in the four-step procedure, we define the nonlocal eigenstress-driven Eshelby tensor nl S σ as follows: such that the nonlocal stress inside the inclusion can be written as Following the procedure outlined in the thought experiment, the constrained displacement in both the inclusion and the matrix after Step 4 is generated by adding the cancellation force, with the same magnitude but in the opposite direction to the surface traction t.The constrained displacement can be obtained by using the nonlocal Green's function or in index notation By taking its gradient and substituting the result into the infinitesimal strain definition, we obtain the explicit expression of the constrained strain c c ij (x) = Further applying the nonlocal constitutive relation, we obtain where the explicit expression of the eigenstress-driven nonlocal Eshelby tensor nl S σ ijkl is clearly indicated.Note that the superscripts σ and nl indicate that the Eshelby tensor is eigenstressdriven and nonlocal.
It can be easily seen that the eigenstress-driven definition is also applicable to local inclusion problems.Indeed, the relationship between the proposed eigenstress-driven and the traditional eigenstrain-driven Eshelby tensor, when considering both of them to be in the local case, can be established through the following equation: l S σ : C = C : l S , (4.17) which indicates that l S σ and l S are conjugate (the left superscript l means local and the right subscript means it is eigenstrain-driven).We merely note that, while this eigenstress-driven formulation is applicable to both local and nonlocal cases, the strain-based Eshelby tensor cannot be properly defined when considering the nonlocal case (see the previous discussion and equation (4.9)).Therefore, the conjugate property breaks when considering the nonlocal case due to the limitations of the eigenstrain-driven definition.

(ii) Analytical formulation
The nonlocal Eshelby tensor can be further simplified.Applying the divergence theorem and changing the variable of differentiation (see the derivation in equation (3.7)), we obtain such that the problem can be reduced to evaluating the following integral: where D ijkl (x) is often called the auxiliary tensor [24].Therefore, the relationship between the auxiliary tensor and the nonlocal Eshelby tensor can be established by V 0 e −iξ y dV(y) For a two-dimensional elliptic inclusion with semi-axes lengths r in = (a, b), the integral Q(ξ ) can be calculated as [24] Q where J (•) is the Bessel function of the first kind and is the Hadamard product [52].Therefore, the auxiliary tensor can be further written as Rewriting the above integral in the polar coordinate system, we have where the radial symmetry condition (see the discussion at the end of §3) is implicitly imposed in k(|ξ |) (note that k(|ξ |) depends on the magnitude of ξ only, not the direction).It is worth noting that the separation between the radial integral and the angular integration allows us to define an intermediate function H(x, r in , θ) that solely accounts for the contribution from the radial integral.By substituting equation (4.24) back into equation (4.20), we complete the derivation of nonlocal Eshelby tensor.

(c) Integrability of nonlocal Eshelby tensor
The integral form of the auxiliary tensor D ijkl shown in equation (4.24) is integrable and can be analytically solved when k(|ξ |) = 1 (which corresponds to the local elasticity case).It follows that D ijkl (x) is constant when x ∈ V 0 , which leads to the well-known result of the Eshelby inclusion problem [1].The integrability of equation (4.24) should be re-examined when considering nonlocal problems.Note that the Fourier integral with respect to the variable |ξ |, which is indicated by H(x, r in , θ ), is an improper integral.The convergence of H(x, r in , θ) is strongly affected by the behaviour of k(|ξ |).Appropriate conditions should be applied to k(|ξ |) in order to ensure the convergence of the improper integral with respect to |ξ |.A practical condition, although either neither sufficient nor necessary, is that This condition is critical because in many cases, the Fourier transform of the kernel functions k(|ξ |) → 0 as |ξ | → ∞.Therefore, it is highly probable that the limit presented in equation (4.25) diverges to infinity and the integral H(x, r in , θ) does not exist.
, ρ ∈ (0, 1), (4.27) which satisfies a certain regularity condition so that its Fourier transform k(|ξ |) converges to 1 − ρ as |ξ | → ∞, not to zero (note that ρ is the ratio of the nonlocal kernel).This approach can be generally applied to most cases of nonlocal kernel functions (such as the Gaussian kernel, see [53]), not the power-law kernel only.
Although both techniques can prevent the divergence of the improper integral, an analytical solution is still unlikely to be obtained due to the complexity of the nested integral in the formulation, and hence numerical methods are required to evaluate the nonlocal Eshelby tensor.
In the next section, we choose the numerical route in order to provide solutions to this problem.Finally, we merely note that while this study focuses on the general solution to the nonlocal Eshelby's inclusion problem, it may indeed be mathematically feasible to design a specialized nonlocal kernel and Q(ξ ) (for instance, by leveraging the strong symmetry in three-dimensional spherical inclusions), such that the overall integral D ijkl (x) in equation ( 4.24) becomes constant or analytically solvable.In fact, studies [54][55][56] have successfully demonstrated that the closed-form solution to the nonlinear Eshelby tensor can be analytically obtained by considering properly selected configurations.

Numerical assessment
In this section, we employ numerical techniques to compute the eigenstress-driven nonlocal Eshelby tensor and examine its characteristics.The selected benchmark problem includes an elliptic inclusion V 0 with semi-axis length r in = (a, b) embedded within an infinite domain.For any point inside and outside the inclusion, we assume a circular nonlocal horizon H 0 of radius r that forms a radially symmetric support of the nonlocal kernel function which is used in the definition of mixed local-nonlocal kernel function as discussed above to ensure the Fourier integral convergence The nonlocal strength governed by the above definition of the kernel function is determined by (1) the nonlocal horizon H 0 , (2) the nonlocal kernel parameter τ and (3) the nonlocal kernel ratio ρ.
Figure 3 shows the geometric configuration of the benchmark nonlocal Eshelby's inclusion problem.Unlike figure 2, in which the nonlocal horizon can be truncated when computing the nonlocal eigenstress, the constrained nonlocal stress is computed throughout the entire material domain with the nonlocal horizon kept constant (also see the discussion on the difference between I(•) and I − (•)).2) aligning the two coordinate axes with the semi-major and semi-minor axes of the elliptic inclusion V 0 .The circular-shaped nonlocal horizon H(x) with radius r is defined at any point x within the whole infinite space domain.θ is the angle between the point position vector x and the unit Fourier vector ξ (which rotates along the unit circle defined in the Fourier space).
In the following, we develop two different approaches in order to (1) solve Eshelby's nonlocal inclusion problem described above, and (2) evaluate the associated nonlocal Eshelby tensor.Specifically, in the first approach, we leverage the multimesh finite-element method (M 2 -FEM) [44] to compute the nonlocal elastic response.In the second approach, we choose to evaluate the nonlocal Eshelby tensor by approximating equation (4.18) directly.

(a) Solution to nonlocal Eshelby's problem via finite-element methods
Using the M 2 -FEM, the overall nonlocal elastic response and the constrained stress within the inclusion are calculated.More specifically, the FE approach is applied to simulate step 4 of the Eshelby's procedure (see §4).Consider a bounded, square-shaped homogeneous nonlocal solid with an elliptic inclusion located at the centre (figure 4).The effect of the eigenstress σ * , which is initially transformed to the surface load t in the thought experiment, is now considered as the external load applied within the inclusion.The corresponding weak-form formulation for the nonlocal Eshelby's problem, based on the simulation configuration, can be derived as follows [44]: where σ in the weak-form formulation is exactly the constrained stress σ c in Eshelby's inclusion problem.The detailed derivation, algorithm and numerical implementation of M 2 -FEM for nonlocal problems can be found in [44].Both σ and its relationship (which in theory, is governed by the nonlocal Eshelby tensor) with the eigenstress σ * can be determined once the problem is solved by FEM.
In the FE simulations, boundary conditions should be selected appropriately in order to approximate the infinite medium considered in the Eshelby's problem (recall that the matrix domain is unbounded, see the description in figure 1).In the following FE simulations, we consider a finite size domain and employ specifically designed boundary conditions to numerically mimic the conditions of an infinite medium.The boundary conditions should be designed carefully in order to minimize their effect on the overall distribution of nonlocal stress σ .Considering that the normal stress and shear stress respond differently to the imposed boundary conditions, in the subsequent sections, we consider two different sets of boundary conditions, named BC1 and BC2, to accurately simulate the distribution of normal stress (σ xx and σ yy ) and shear stress components (σ xy ).The choice of different boundary conditions is described in figure 4 and numerical tests assessing the sensitivity of the stress components to the different boundary conditions are included in the electronic supplementary material.
where the Fourier transform of the nonlocal kernel function can be obtained by substituting the truncated one Knl (x) defined in equation (5.1) where J 0 (•) is the Bessel function of the first kind.Substituting equations (5.4,5.5)back into the Fourier integral in equation (4.24), we obtain Here, we also use the superscript mix to denote the mixed local-nonlocal definition.The above expression involves improper integrals (integration with respect to |ξ | ∈ [0, ∞)), nested integrals (integral with respect to |x| ∈ [0, r]) and transcendental functions (Bessel function J 0 (•) and J 1 (•)), which make the integral extremely difficult to evaluate.We use the technique of upper integral limit truncation in order to approximate the above improper integral.Note that the Bessel function J 0 ( ( By substituting equation (5.7) into equation (5.6), we obtain the following approximation: where H nl (•) and H l (•) are the nonlocal and local parts defined as and Equations (5.9a) and (5.9b) can be solved by different approaches.The local part, as proved by Eshelby, can be analytically determined as constant when x ∈ V 0 [24].The nonlocal part, although it includes nested integrals and transcendental functions, can be approximated by using numerical integration methods.The auxiliary tensor D and the nonlocal Eshelby tensor nl S σ can be evaluated using equation (4.24) once H mix is obtained.

Simulation results and analysis
We consider two numerical cases, one involving a circular inclusion and another involving an elliptic inclusion.Numerical results from both the FEM and DNC are presented and compared to substantiate our findings.Additional parametric studies and analyses are conducted to further investigate the convergence behaviour of the Eshelby tensor (from nonlocal to local) by varying nonlocal parameters.Material elastic properties in all the following simulations are fixed to be μ = 1[Pa] and λ = 1[Pa], and the eigenstress: 0.04 0.02 0.02 0.04 , ( is used in all the following computations of constrained stress and FE simulation.Note that the choice of using arbitrary values for both the material properties and the eigenstress does not hinder the generality of the results.Other simulation setups, including the reference local value of constrained stress, the eigenstrain-driven local Eshelby tensor and the eigenstress-driven local Eshelby tensor are also listed in the electronic supplementary material.

(a) Circular inclusion
In this numerical case, we consider a circular inclusion with radius r in = 0.

(i) FEM results
As discussed in the previous section, when using the FE approach the normal and shear stress response is obtained for the nonlocal inclusion problem by considering two different boundary conditions BC1 and BC2 specified in figure 4. The prescribed constant eigenstress inside the circular inclusion, which is treated as the external stress load, is given in equation (6.1). Figure 5 presents the nonlocal stress distribution within the whole square-shaped domain.A first glance at the M 2 -FEM simulation results indicate that all three nonlocal stress components exhibit a visible dependence on the presence of the inclusion.Specifically, the stress distribution  inside the inclusion tends to be more uniform compared to the distribution outside.This characteristic is aligned with the elastic response in local problems, where the stress/strain distribution is found to be homogeneous within the entire inclusion (see the simulation results obtained from COMSOL in the electronic supplementary material).The stress components inside the inclusion, as predicted by the M 2 -FEM, closely match the reference values obtained from the local inclusion problems as well.As highlighted in figure 5 (see the blue cross sign and the blue value displayed on the colour bar), all three values of stress components at the centre of inclusion are in close agreement with the reference value provided in the electronic supplementary material.
The primary difference between local and nonlocal simulation results is the perturbation of the stress distribution located near the boundary of the inclusion (see the shaded region in figure 5).As previously discussed in §4, the presence of the nonlocal horizon separates the original inclusion into two distinct sub-regions: the reduced inclusion V − and its complement V c − .The stress distribution, as illustrated in the zoom-in inset, becomes increasingly heterogeneous as points move from V − into V c − , and eventually exhibits a discontinuity when crossing the inclusion boundary.The distinct behaviour observed in the nonlocal simulation is a result of the interaction between the nonlocal horizon and the inclusion boundary.For any point x ∈ V c − , the corresponding nonlocal horizon H(x) is not a subset of the inclusion V 0 ; namely, there exists a separate part of H(x) with varying shape and area (depending on the position of x) that is located outside of V 0 .This separate portion of H(x) captures the discontinuous stress distribution across the inclusion boundary and eventually leads to the perturbation of stress for points inside V c − .During practical simulations, the impact brought by stress discontinuity is less pronounced when nonlocal effects are not significant enough and the stress distribution tends to be less heterogeneous inside V c − (see the ratio between the shaded region area and the total area of V c − between two dashed lines in figure 5).The impact of nonlocal effects and stress distribution will be discussed in depth in the convergence study.
(ii) DNC results from the reduced horizon and eventually becomes heterogeneous and anisotropic when it exceeds the inclusion boundary (highlighted by the red dashed line).Again, we remark that the strong dependence on the inclusion boundary observed in both FEM and DNC results is due to the non-resolvable nonlocal horizon and its interaction with the inclusion boundary.
The constrained nonlocal stress distribution exhibits a similar inclusion-boundary dependency.It is worth noticing that, compared with the nonlocal Eshelby tensor, the distribution of constrained nonlocal stress shown in figure 7 is more sensitive to the inclusion boundary: the banded area that shows the inhomogeneous distribution of nl σ c near the inclusion boundary is much wider than the band shown in figure 6.This is because nl σ c is computed by taking the nonlocal integral of nl S σ (see the definition in equation (4.11)) over the horizon which covers points outside the inclusion boundary (figure 2).In this regard, the constrained nonlocal stress depends not only on nl S σ inside V 0 (which is only slightly affected by the inclusion boundary), but also on nl S σ outside V 0 (which exhibits strong heterogeneity and anisotropy, see figure 6).A direct comparison between figures 5 and 7 reveals that the distribution of the constrained nonlocal stress predicted by both FEM and DNC are in good agreement.This consistency confirms the validity of the derived analytical expression of the nonlocal Eshelby tensor in equation (4.20) and its effectiveness in predicting the constrained nonlocal stress.Furthermore, the proximity of nl σ c within the reduced inclusion (predicted by both methods) to the reference value in the local case confirms the postulate proposed earlier in §4 (that is, nonlocality should affect the Eshelby tensor less at points that are closer to the centre of the inclusion, as the stress perturbation at those points is less affected by the inclusion-matrix interface).

(b) Elliptic inclusion
In this numerical case, we simulate a more general nonlocal inclusion problem consisting of an elliptic inclusion.The semi-axis length is selected as (a in , b in ) = (0.8, 0.6)[m].Boundary conditions and nonlocal configurations (including horizon size, nonlocal parameter and nonlocal ratio) are all set to be the same as in the circular inclusion case (see the specification presented at the beginning of §6).

(i) FEM results
Figure 8 presents the FE simulation results for the elliptic inclusion problem.It can be observed that the nonlocal stress distribution in the elliptic inclusion case exhibits similar inclusiondependent behaviour to that identified in the circular inclusion case.The nonlocal stress inside the inclusion, particularly when inside the reduced inclusion, is uniformly distributed and its value is very close to that found for the local inclusion problem (see electronic supplementary material).Notably, the consistency observed between the circular and elliptic cases for the nonlocal problem aligns well with Eshelby's findings in the local inclusion problem, hence substantiating the effectiveness of the stress-driven nonlocal Eshelby tensor and its efficacy in extending the Eshelby tensor formalism to nonlocal cases.
The anisotropic behaviour of the nonlocal stress distribution is also well predicted by the M 2 -FEM.Note that, unlike the circular inclusion which shows rotational symmetry up to an infinite order, the elliptic inclusion only possesses a twofold rotational symmetry.(ii) DNC results   (

c) Convergence study
We also explore numerically the convergence behaviour of the nonlocal Eshelby tensors.Specifically, we investigate the variation of the constrained stress field and the nonlocal Eshelby tensor distribution by varying the nonlocal kernel parameter τ .These numerical simulations provide insights into the sensitivity of the results to changes in these parameters, allowing for a better understanding of the behaviour of the nonlocal Eshelby tensors under different conditions.To streamline the analysis, we focus on circular inclusion problems, keeping both the geometric configurations and the material properties consistent with those discussed in §5a.This simplification allows us to isolate and study the convergence behaviour of the nonlocal Eshelby tensors more effectively.xx | in figure 11 to better visualize the stress variation with respect to τ and study the convergence behaviour.It can be observed that, as τ increases from 0.002 to 0.01, the peak value of the stress difference (see the bright yellow colour) located near the inclusion boundary becomes smaller and smoother, but the area where stress oscillations occur around the inclusion boundary becomes larger (see the increasing rim thickness in each subfigure from a to c).These observations suggest that the presence of nonlocality tends to smoothen the discontinuity of stress across the inclusion boundary.Additionally, as the strength of nonlocality increases, the smoothing effect becomes more significant and affects a larger area.Correspondingly, the smoothing effect reduces as τ decreases and eventually vanishes when the problem transits from nonlocal to fully local.
The convergence behaviour can be further identified from the results of the nonlocal Eshelby tensor presented in figure 12 ).Notably, as τ decreases from 0.008 to 0.002, the nonlocal effect also reduces so that the curve of the nonlocal Eshelby tensor is less distorted, which eventually converges to the reference local value.

Summary and outlook
This study focused on the conceptualization, definition, analytical solution and numerical assessment of a nonlocal-enriched formulation of Eshelby's inclusion problem.The main objective was to extend the classical Eshelby's theory to incorporate nonlocal effects and provide a comprehensive understanding of the effects of inclusions in nonlocal media.The conceptualization involved introducing the eigenstress-driven approach, which adopts the constant eigenstress assumption, as an alternative to the traditional eigenstrain-driven approach, in nonlocal elasticity.The definition of the nonlocal Eshelby's inclusion problem started from formulating the problem via Eringen's integral nonlocal elasticity and then continued by incorporating the eigenstress-driven four-step thought experiment.Analytical solutions to the nonlocal inclusion problem, which are mathematically represented by the eigenstress-driven nonlocal Eshelby tensor, were derived using mathematical tools including the nonlocal Green's function and the Fourier transform approaches.The resulting eigenstress-driven nonlocal Eshelby tensor provided a concise and effective representation of the relationship between the nonlocal constrained stress and the nonlocal eigenstress field within the inclusion.The validity of the proposed theory and the nonlocal Eshelby tensor was assessed through numerical simulations by using both DNC and finite-element methods.The excellent correspondence between the two approaches provided evidence of the validity and effectiveness of the eigenstress-driven nonlocal Eshelby tensor formulation.In conclusion, the proposed nonlocal inclusion theory provides a powerful framework to study and analyse the mechanical behaviour of nonlocal media, and serves as the theoretical foundation for the future development of advanced homogenization methods for nonlocal heterogeneous media.Future work will focus on extending the eigenstressdriven approach to more general nonlocal inclusion problems including three-dimensional and nonlinear problems.
Data accessibility.The additional data in the form of equations and discussions supporting this article have been uploaded as part of the electronic supplementary material.The code can be accessed here: https://zenodo.org/records/10827725 [57].Supplementary material is available online [58].writing.After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the publication.

Figure 1 .
Figure 1.Four-step thought experiment as proposed by Eshelby.
kmlj (x, y) is now transformed into the volume integral form similar to T (a) kmlj (x, y).Note that the key and only difference between the updated form of T (b) kmlj (x, y) and T (a) kmlj (x, y) shown in equation (3.4) is the partial derivative in the jth direction: in T (b) kmlj (x, y) the partial derivative is calculated with respect to x , while in T (a) kmlj

. 5 )
where I − (•) represents the modified nonlocal integral operator with truncated support H − (see the definition in equation (4.2)).This step demonstrates the importance of making the constant eigenstress assumption in a nonlocal Eshelby's inclusion problem: the nonlocal stress inside the inclusion σ I is constant only when the eigenstress σ * is constant.This also confirms that the constant eigenstress assumption is more convenient to use than the constant eigenstrain assumption when analysing the nonlocal Eshelby's inclusion problem.The surface traction is obtained by using the definition of Cauchy stress t = −σ * • n.(4.6) ) contains two different nonlocal integral operators I(•) and I − (•) (for the sake of simplicity, we use I(•) to denote the original integral operator with non-truncated support H 0 ).Apparently, the additivity between I(•) and I − (•) is not applicable given that the two integral operators are inconsistent as they are defined on different supports.Clearly, in the limit to the local case (i.e. when the nonlocality vanishes), both integral operators become identity operators, and equation (4.9) becomes σ I = C : ( c − * ) = C : ( l S − I) * , (4.10)

Figure 3 .
Figure 3. Geometric configuration of the nonlocal Eshelby's benchmark problem.The Euclidean coordinate system in the space domain is established by (1) setting the origin point at the centre, and (2) aligning the two coordinate axes with the semi-major and semi-minor axes of the elliptic inclusion V 0 .The circular-shaped nonlocal horizon H(x) with radius r is defined at any point x within the whole infinite space domain.θ is the angle between the point position vector x and the unit Fourier vector ξ (which rotates along the unit circle defined in the Fourier space).

Figure 4 .
Figure 4.The two different sets of boundary conditions used in FE simulations to approximate an infinite domain.(a) Combination of fixed and free boundary conditions (BC1) that are applied to study the behaviour of the normal stress components σ xx and σ yy in Eshelby's nonlocal inclusion problem.(b) Combination of fixed and free boundary conditions (BC2) applied to study shear stress component σ xy .Note that fixed boundary conditions (u x = u y = 0[m]) are marked by slanted lines; free boundary conditions are indicated by the free edge (no slanted lines).Note that, in each configuration, the length of each boundary is the same on each side of the domain.The size of the square-shaped matrix medium is set to be 4[m] × 4[m].
7 [m].The size of the nonlocal horizon at each point, either inside the inclusion or the matrix, is set to be r = 0.2 [m].The two nonlocal and non-dimensional parameters that define the mixed nonlocal kernel function (see equation (5.2)) are given as τ = 0.002 and ρ = 0.8, respectively.

Figure 5 .
Figure 5. FE simulation of the nonlocal stress distribution when considering a circular inclusion: (a) distribution of normal stress σ xx , (b) distribution of normal stress σ yy and (c) distribution of shear stress σ xy .The circular inclusion V 0 and the reduced inclusion V − are outlined by the dotted red line and the dotted grey line, respectively.Specific nonlocal stress values at the centre of the inclusion (labelled by the blue cross sign) are marked on the colour bar with blue colour.The shaded red area in the zoom-in inset indicates the perturbed region induced by nonlocality.

Figure 8 .
Figure 8. Nonlocal stress distribution obtained via FE analysis for the elliptic inclusion case: (a) distribution of normal stress σ xx , (b) distribution of normal stress σ yy and (c) distribution of shear stress σ xy .The circular inclusion V 0 and the reduced inclusion V − are outlined by the dotted red line and dotted grey line, respectively.Specific nonlocal stress values at the centre of the inclusion (labelled by the blue cross sign) are marked on the colour bar with blue colour.The shaded red area in the zoom-in axis indicates the perturbed region of stress induced by nonlocality.

Figure 9 .Figure 10 .
Figure 9. Distribution of nonzero components of the nonlocal Eshelby tensor in the extended elliptic horizon V + .Note that the nonlocal Eshelby tensor preserves the minor symmetry only in the reduced elliptic inclusion V − .

Figures 9 and 10
Figures 9 and 10 present the DNC results for the elliptic inclusion case.Again, similar to the DNC results in the circular inclusion case, both the nonlocal Eshelby tensor ( nl S σ ) and the constrained nonlocal stress ( nl S σ ) in the elliptic inclusion are found to be almost uniform inside the reduced inclusion and heterogeneous when approaching the inclusion boundary.As discussed above, the major difference in the nonlocal Eshelby tensor between the circular and elliptic inclusion cases is their symmetry.Specifically, due to the shape of the inclusion, the nonlocal Eshelby tensor in elliptic inclusion problems loses major symmetry within the whole inclusion.This can be verified by comparing the distribution of S 1122 and S 2211 (see sub-figures (b1) and (b2) in figure9).The constrained nonlocal stress distribution obtained by DNC and FEM shows excellent agreement also in the elliptic case.Specific values of each stress component located within the reduced inclusion are all in alignment with the reference values in the corresponding local inclusion problems.Combining the simulation results from both circular and elliptic inclusion cases, we conclude that the stress-driven nonlocal Eshelby tensor exhibits good mathematical properties and serves as a viable extension of the classical local Eshelby tensor.

Figure 11 .Figure 12 .
Figure 11.Distribution of the absolute difference between nonlocal stress and local stress log |σ xx − σ 0 xx | with different nonlocal kernel parameter τ : (a) τ = 0.002, (b) τ = 0.006 and (c) τ = 0.01.Simulation results show that the perturbed area of stress near the inclusion boundary becomes larger and also smoother as τ increases.Note that all the data are on a logarithmic scale for visualization.
Figures 11 and 12  present the variation of the constrained nonlocal stress σ xx (computed using FEM) and nonlocal Eshelby tensor (computed using DNC) when applying different nonlocal kernel parameters, respectively.