‘Explicit’ and ‘implicit’ non-local continuous descriptions for a plate with circular inclusion in tension

Increasing application of composite structures in engineering field inherently speed up the studies focusing on the investigation of non-homogeneous bodies. Due to their capability on capturing the size effects, and offering solutions independent of spatial discretization, enriched non-classical continuum theories are often more preferable with respect to the classical ones. In the present study, the sample problem of a plate with a circular inclusion subjected to a uniform tensile stress is investigated in terms of both ‘implicit’/‘weak’ and ‘explicit’/‘strong’ non-local descriptions: Cosserat (micropolar) and Eringen theories, by employing the finite element method. The material parameters of ‘implicit’ model is assumed to be known, while the nonlocality of ‘explicit’ model is optimized according to stress concentration factors reported for infinite Cosserat plates. The advantages/disadvantages, and correspondence/non-correspondence between both non-local models are highlighted and discussed apparently for the first time, by comparing the stress field provided for reference benchmark problem under various scale ratios, and material parameter combinations for matrix-inclusion pair. The results reveal the analogous character of both non-local models in case of geometric singularities, which may pave the way for further studies considering problems with noticeable scale effects and load singularities.

Non-locality, by definition, implies the presence of internal lengths and spatial dispersion properties in wave propagation [22]. Within this framework, as suggested in [13,14,23,25], we can distinguish between 'implicit'/'weak' and 'explicit'/'strong' nonlocal formulations. The former is referred to continua with extra degrees of freedom of various kind (i.e. [12,14,15]), in which there exist additional equations of motion containing non-standard primal fields, named generalised continua, microcontinua, or even multifield continua. While in the latter (i.e. [16,20,21]) primal fields of classical elasticity are preserved, yet the equations of motion contain integral, integro-differential or finite-difference operators in the spatial fields. In the end, both approaches avoid physical inadequacies, and mesh-dependency in numerical solutions [1,2].
Among 'implicit' non-local models, many papers showed the advantages of micropolar theory, which has been widely exploited for describing heterogeneous materials with microstructures (e.g. masonrylike structures, granular, blocked and layered materials, rock assemblages, reinforced composites, and so on [26][27][28][29][30][31][32][33]). The micropolar model was first developed by Cosserat brothers [17], and improved by many researchers over the years such as, Eringen [18], Nowacki [19], etc. The theory accounts for the rotation of individual material points by integrating size effects through an additional kinematic and work-conjugated dynamic descriptors representing the material microstructure. Respectively, the stress tensor is expressed in terms of skew-symmetric stress and couple-stress parts. Due to its limited non-local character, the theory is classified as 'implicit' or 'weak' non-local model. Meanwhile, Eringen's nonlocal theory, which was first introduced by Kröner [34], Krumhansl [35], Kunin [22], and Eringen [36] and were further improved by Eringen and Edelen [20], and Eringen [16,37,38], is concerned with the physics of materials bodies whose behaviour at a material point is influenced by the state of all points of the body, hence extensively employed to investigate the mechanical behaviour of nano/micro sized structures (e.g. carbon nanotubes, graphene sheets, nonhomogenous elastic continua, etc.) [39][40][41]. The theory possesses 'explicit' or 'strong' non-locality, since the small-scale parameter is integrated directly into the convolution type constitutive equation to capture the size effects. Accordingly, stress at each point is related to the strain of entire domain through an attenuation-type kernel function that contains information about internal length. As the internal lengths may have many different interpretations depending on the structure [16,23], 'explicit' nonlocal continua may also be used represent the material microstructure like 'implicit' one, even if it has often not.
As an inevitable outcome of widespread application of composite materials, the studies that have focused on quantifying the effect of heterogeneities through investigation of interface between inclusion and matrix materials gained speed. This problem has been widely studied considering both classical [42][43][44] and non-classical theories via employing different solution techniques. As examples of early studies focused on non-local models, Lubarda [45] provided the analytical solutions in the context of anti-plane strain couplestress elasticity, while Zhang and Sharma [46] employed strain gradient elasticity with couple stresses. From a different point of view, Dong, Wang and Rubin [47], treated the interface as a finite thickness Cosserat shell/membrane medium and compared their results with explicit solutions. More recently, Atroshchenko et al. [48], modelled the matrix-inclusion structure as Cosserat media considering various material properties and showed the effect of imperfect and perfect interfaces via employing boundary element method (BEM), while Fantuzzi et al. [32] investigated the classical matrix-inclusion problem considering micropolar theory using both strong and standard forms of the finite element method.
Due to its rather complex character in the implementation, integral form of Eringen's non-local theory based finite element (FE) formulation, originally derived by Polizzotto [49], has been only conducted by a limited number of researchers. Pisano, Sofi and Fuschi [50] are the first ones, who implemented the FE formulation in a two-dimensional (2D) nonlocal elasticity context for solving boundary-value problems. Inherently, the inclusion problem was first addressed by the same authors on the basis of a socalled ''enhanced'' strain-difference based Eringentype nonlocal model [39] via modelling a piecewisehomogenous domain [51,52]. Although the effect of the non-local material parameters on strain field was widely and clearly discussed, the studies are limited with square inclusion and the stress evaluation of mentioned problem is only recently examined by Pisano and Fuschi [53], who obtained spurious stress oscillations around the inclusion when full integration is considered.
After the descriptions of 'implicit' and 'explicit' non-locality had been suggested [13,14,23], Trovalusci [25] has focused on their definitions in a comparative way, for the first time. Yet, these different descriptions of non-locality have not been employed to study a common problem to look for possible correspondences and differences. With this motivation, in the present study, the effect of a circular inclusion as in particle and/or fibre reinforced composite materials is investigated on the basis of both Cosserat ('implicit') and the integral form of Eringen's two-phase local/nonlocal ('explicit') theories under uniaxial loading by employing the finite element method. Hence, Eringen's theory is not here used to model solids at atomistic scale, but rather the focus is on structures at a larger scale dominated by other kinds of discrete nature (e.g. heterogeneity) yielding an apparent non-local mechanical behaviour. The inclusion is incorporated into the model as a heterogeneity possessing different material properties than the surrounding matrix. The intrinsic length scale parameter of Cosserat model is assumed to be known while the non-locality of Eringen model is optimized using an evolutionary algorithm to minimize the difference between stress concentration factors obtained for infinite 'implicit' and 'explicit' non-local plates. After constitutive parameter detection of 'explicit' model, the plates are investigated under various material combinations for matrix-inclusion pair. To mimic the effect of both finite and infinite plates, wide range of edge lengths for a fixed radius of inclusion is considered. To the best of our knowledge, this represents the first attempt to investigate the possible analogous/dissimilar characters between conducted 'implicit' and 'explicit' non-local continuum theories in the simplified context of a two dimensional inclusion-matrix problem under plane strain assumption.
The paper is organized as follows. In Sect. 2, the governing equations are outlined, and corresponding finite element formulations are derived in the context of both 'implicit' and 'explicit' non-local theories. Section 3 is devoted to numerical examination of matrix-inclusion problem, alongside with a brief explanation about optimization process. The comparison of 'explicit' and 'implicit' type non-local theories and the influence of small-scale parameters are also discussed focusing on the normal stress field. Finally, some concluding remarks about capabilities, advantages/disadvantages of both non-local models are presented in Sect. 4.

Materials and methods
Continuum theories considered in the present study are briefly explained, and the displacement-based finite element formulations are derived regarding twodimensional plane-strain problems within the linearized kinematical framework. Matrix and inclusion are assumed as linear, elastic and isotropic materials possessing different mechanical properties. FE formulation of micropolar continuum model is implemented within the environment of the software COMSOL MultiphysicsÓ [54], while for Eringen's nonlocal model an in-house code is developed using Wolfram Mathematica [55]. The superscripts E and M, that refer to Eringen's and micropolar model, respectively, are used to distinguish the parameters appearing in both theories, and possessing different interpretations.
The body under investigation can be regarded as a set of material points in 3D Euclidean space, occupying a domain X 0 , and enclosed with a boundary, C 0 which is subjected to traction vector, t, and coupletraction vector, m(considered for Cosserat media only). A Cartesian coordinate system xyz is used with a suitable origin for the parameterization of positions of material points, x. In the two-dimensional framework (2D) it is assumed that X 0 is formed with a uniform and symmetric thickening of a 2D region X 2 x; y f g, enclosed with boundary C, along z axis by an amount of h.

Overview
In micropolar theory, material particles that constitute continuum are described in terms of both their positions and orientations leading to following linearized kinematic relations: where / k stands for the components of micro-rotation vector, and e M ij and v kj denote the components of strain and curvature tensors with e ijk being the usual thirdorder permutation tensor. One may recognize that, if the micro-rotations are constrained to follow the local rigid rotations, (macro-rotations): the classical kinematic relations, i.e., e ij ¼ u i;j þ u j;i À Á 2, will be recovered and the continuum becomes a couple-stress continuum [26,56].
From balance considerations, each component of surface traction and surface couple-traction, denoted by t M i and m k , respectively, are described as: where r M ij and l kj are the components of the nonsymmetric stress and couple-stress tensors, respectively, and n j being the components of unit outward normal vector of the continuum boundary. Hence, if body forces and body couples are neglected, the equilibrium equations take the following form.
Considering linear elasticity, the stress-strain relation of an isotropic micropolar continua can be represented as: which requires six elastic material constants for the complete description, i.e., Lamé's constants k and l, and four additional parameters, a, b, c and v related to micropolar theory [57]. It is important to note that, in the present work, second Lamé's parameter, l, is not equal to shear modulus, G, as in the classical convention, but holds following relation: Accordingly, Poisson's ratio, m, can be described as follows: In Cosserat model, the size effects and relative rotations are incorporated through an internal characteristic length, l c , and a coupling number, N, in such a way that: If these parameters (i.e. l c , N) are taken small enough. the Cosserat effects become negligible and the body behaves as Cauchy continua. Note that this finding is valid only in the case of materials with at least the orthotetragonal symmetry, or materials belonging to more restricted symmetry classes, as of course isotropic materials [26,58].

Finite element formulation
For FE modelling, the field variables within an element e (i.e. u e , / e ) are approximated considering a natural coordinate system f; g and using related interpolation function matrices: where the over tilde symbol is used to indicate the nodal values. Since the present study restrained with 2D media only, the degrees of freedom (DOFs) per each point is reduced to two in-plane displacement components along x and y directions, i.e., u x and u y , and one out-of-plane micro-rotation component along z direction, i.e., / z . N M u is employed for interpolation of nodal in-plane displacements and consists of quadratic interpolation functions, while N / is employed for nodal out-ofplane micro-rotations, and includes linear type shape functions. As it was mentioned in Fantuzzi et al. [32], and Fantuzzi et al., [33], and clearly showed in the Fig. 1a, all the nodes of the nine-node Lagrange element possess displacement-type DOFs, whereas micro-rotation DOFs are attached only to the four corner nodes. As a matter of fact, it has been demonstrated that the use of a lower interpolation order for micro-rotations allows an optimal compromise between numerical accuracy and computational efficiency to be achieved. In other words, the adoption of the same interpolation functions for both displacements and micro-rotations usually leads to negligible accuracy improvements against a notable increase in the number of DOFs. Hence corresponding shape function matrices take the following form: fields, which are obtained using differential matrix operator L M , permutation vector M, and gradient operator r: are derived as follows where d M e indicates the unknown vector of nodal displacements, and the superscripts appeared in Eq. 12 3 (and also in Eq. 10) refer to node numbers, while the matrices B M ee and B M ev include the derivation of interpolation functions. Since N M u and N / are taken as functions of natural coordinates, the derivations with respect to physical coordinate systems given in Eq. 11 1 and Eq. 11 3 should be performed using chain rule, and inverse Jacobian matrix of element e, J À1 e . Inherently, constitutive relations given in Eq. 5 are transformed to: where r M e and l e vectors, described as denote non-symmetric stress field and couple stress field of an element, while constitutive matrices, D M ee and D M ev are defined as follows:  Fig. 1 The a 9-node quadrangular, b 4-node quadrilateral element illustrated in A standard displacement-based FE formulation is provided using the principle of minimum total potential energy (i.e.
where for the static analysis, total potential energy functional,P M , which is written in terms of total elastic strain energy, U M and external work potential W M , must be minimum for equilibrium.
in the absence of body forces and couples ð Þ with N total , A m and d M being the total number of elements, area of mth element, and global nodal displacement vector. Due to 'weak' non-local character of constitutive equation, element formulation of an element m is derived as follows, by including the derivative of external work potential, f M m : By performing proper assemblage operations, the well-known linear equation system is achieved: where K M refers to global stiffness matrix and f M is the global nodal force vector of discretized 2D micropolar media. For integration operations in Eq. 17, a standard Gauss integration technique (GQ) with 3 9 3 Gauss points for second order element, is adopted.

Overview
The assumptions of linear elasticity provide the following kinematic relations for Eringen's nonlocal model: where e E ij refer to components of the strain tensor. For the continuum to be in balance, interactions between material points, that are characterized through traction forces only, t E , can be described as follows: where r E ij refers to components of the symmetric stress tensor. In the absence of body forces, the following equilibrium equation is carried out, similar to classical (local) elasticity theory: As a matter of fact, the localisation of global (integral) balance laws (e.g. mass, momentum, moment of momentum, energy) would introduce non-local residuals (localisation residuals), due to the long-range effects [20,21,37]. This contradicts with global balance laws being valid for each infinitesimal volume element isolated from the body [16,20,21,37]. However, by a new interpretation of stress tensor and internal energy density, the local form of balance laws is recovered only with additional jump conditions [16]. Moreover, axiom of objectivity provides the vanishing of non-local body force and couple residuals. As the only remaining residual associated with non-locality, strain-rate dependent non-local energy residual [39,59], dies out under the quasi-static condition, following constitutive equation for linear elastic isotropic solids is recovered [16,37].
To avoid any conflict with its previous interpretation in micropolar model, constitutive equation is expressed using shear modulus, G, instead of second Lamé's parameter, l, while first Lamé's parameter, k, is retained. The convolution type constitutive equation of the integral form of Eringen's two phase local/ nonlocal model implies that the stress at a point is linked to the strain of the entire domain through a kernel function, s, which is assumed to be bi-exponential in the present study, in accordance with literature [60,61], and represented as follows: where r (i.e. r ¼ x À x j j) refers to the Euclidian distance between investigated point of interest, x, and its neighbour points, x, [49,51,52], and d being the dimension of the structure, while j stands for nonlocal parameter. As one can see, two additional parameters, n and j, both of which attribute to non-local character of structure, are introduced with Eringen's non-local theory. The former symbol stands for fraction coefficient, that ranges from 0 to 1, and is used to regulate the weights of local and non-local parts, while the nonlocal parameter that depends on both internal length scale, a 0 , and material constants e 0 : defines the characteristic of the kernel function. Hence, assuming n equals to 1 yields full local model, and 0 corresponds to full non-local one, while the increasing values of j produces more pronounced non-locality.

Finite element formulation
For FE modelling of 2D non-local media, the domain is discretized using four-node quadrilateral elements (see Fig. 1b), and displacement field within an element e (i.e. u e ) is approximated using interpolation matrix including linear shape functions, N u , and nodal displacement vector, d E e as in the form: where L E and B E e refer to differential operator, and strain matrices. As previously mentioned, the derivations with respect to physical coordinate systems given in Eq. 26 2 are performed by the chain rule, and employing the inverse of Jacobian matrix of element e, J À1 e . The constitutive relation given in Eq. 22 is then transformed to: while stress tensor and constitutive matrix are described as below: Note that the over bar in Eq. 27, denotes that the related matrix is written in terms of f; g as a requirement of convolution form of Eringen's constitutive relation, hence, should not be confused with overbars in prescribed surface traction, t, and surface couple-traction, m.
Finally, the weak form of displacement based FE formulation is derived based on minimum total potential energy principle (i.e. oP E d E Â Ã od E i ¼ 0; i ¼ 1; 2; . . .; N total ð Þ ) with following definitions for elastic strain energy U E and external work potential W E .
Due to 'strong' non-local character of Eringen's nonlocal model, the element formulation is provided as follows: The subscripts m and n stand for element number. The term with coefficient n represents the local part of two-phase model, while the terms with coefficient 1 À n correspond to the non-local part. In the non-local part, the first term stands for the contribution of the mth element to its own energy, while second and third terms account for the influence exerted on the mth element by the remaining elements, and the influence exerted by the mth element to other elements, respectively. Second and third terms in the non-local part, which expand element stiffness matrix to overall dimension are equal to each other only for homogenous material properties: One can see that, Eq. 30 has the same nature with the corresponding ones reported in literature [47], and to be accordance with the literature, k E mn and k E nm are called as cross-stiffness matrices from now on. At last, the following linear equation system is obtained after proper assemblage operations.
The integration operations given in Eq. 30 are performed using GQ method. For the local part (i.e. k m ), 2 9 2 Gauss sampling points (GP) are sufficient for first order elements, while for the nonlocal part, the number of GP varies depending on nonlocal parameter, element length and calculated part of stiffness matrix (i.e. k mm or k mn ). To reduce the computational burden, only the elements inside the influence zone (see Fig. 2) are assumed to contribute to crossstiffness matrices k mn and k nm , while the effect of others is negligible due to the decaying nature of attenuation function [50,62]. It should also be noted that, this simplification inherently leads to some sort of banding of the global stiffness matrix K E .

Numerical simulations
This section is devoted to the numerical examination of linear elastic plates with a circular inclusion under uniform tensile stress. As loading condition, a constant tensile stress of 100 MPa is applied to left and right vertical boundary edges, while boundary conditions are imposed by limiting the horizontal and vertical displacements of vertical and horizontal axes located in the middle, respectively (Fig. 3).
Despite the symmetric character of the problem (along both x and y axes), analyses are performed regarding the full domain, since in the context of Eringen's nonlocal theory, imposing symmetry is much more complicated than only considering the symmetric part of the model [63]. It is a result of the missing neighbour elements that should have contributed to cross-stiffness matrices due to long-range interactions. The simulations are repeated considering two different mesh configurations of different refinement, both characterized by different scale ratios: for all models, respectively named Model 1, 2, 3, the radius of inclusion, a, is assumed constant, 0.05 [m], while the width of the plate (structural/macro length), L, is modified in size: L/a = 3.0, 10.0, 20.0 to consider both finite and infinite domains (see Fig. 4). As the main purpose of the present study, the effect of different material configurations (see Table 1) is studied in the context of two distinct non-classical theories, Cosserat (micropolar) and Eringen' nonlocal models, as described below: MAT 1: Both the matrix and the inclusion are modelled as Cauchy.
MAT 2: The inclusion is Cauchy and the matrix is Cosserat/Eringen.

MAT 3: The inclusion is Cosserat/Eringen and the matrix is Cauchy.
MAT 4: Both the matrix and inclusion are modelled as Cosserat/Eringen.
In this way, it is aimed to examine the incorporation of size effects through 'implicit' and 'explicit' nonlocal, scale dependent continua, specifically focusing on stress concentration factor (SCF) as one of the main design parameters for the problem under investigation:   where r 0 and r max correspond to nominal stress and maximum stress (e.g. obtained at interface point A), respectively. Note that, for Eringen's nonlocal theory, the inhomogeneity of the structure is modelled only by different values of fraction coefficient and Young's modulus, and spatial distance is used to evaluate the attenuation function regardless of the material properties. However, for nonhomogeneous nonlocal models, more refined models are available. For instance, [39] considers the inhomogeneity as a source of additional attenuation effects, and incorporates it by introducing 'equivalent distance' and 'reference internal length', rather than only incorporating the nonlocal stress-strain law with different material constants.
For the transition between Cauchy (classical) and Cosserat/Eringen (non-classical) models, different strategies are developed based on the conducted continuum theory. For micropolar theory, it is achieved via tailoring the value of internal characteristic length, l c and coupling number, N, as given in Table 1. For small values of corresponding material parameters, the size effect becomes negligible, and the model behaves as a Cauchy continuum, while an opposite trend is obtained for larger values. On the other hand, when Eringen's nonlocal theory is considered, the locality and/or non-locality of inclusion and matrix can simply be arranged by changing the fraction coefficient, n, and/or the non-local parameter, j, appearing in the constitutive equation. For a fully local media, fraction coefficient should be equal to unity (i.e. n = 1), while any smaller value (i.e. n \ 1) lets the incorporation of size effects. To modify the non-locality of Eringen's model in accordance with micropolar theory (regarding the sample problem), an evolutionary algorithm called Differential Evolution Method [64], which has been widely employed as a global optimization tool in engineering problems, is utilized. Details of the identification procedure can be found [65], in which all the steps of optimization algorithm are presented in detail for determination of scale parameter of carbon nanotubes and atomic chains. In this study, it is intended to reduce the relative difference between stress concentration factors obtained through Eringen's nonlocal theory based FE models (for L/a = 20.0) and explicit expressions of stress concentration factors reported for infinite plates employing Cosserat theory [48]. Consideration of SCF to construct the objective function is very advantageous computationally; however, other strategies considering global behaviour are possible, such as minimizing the differences of strain energy or stress distribution over the domain. Then again, focusing on a single parameter as it is done herein turns out to be suggestive enough to look for possible hints on equivalencies between two different nonlocal descriptions of solids, which is the main aim of the study. As it is less expensive computationally, the parameter to be optimized is selected as n (0:0 n 1:0), while j is taken as 0.2a m for all the models, and the radius of influence zone i.e. R I , is assumed as 0.05 [m] for given j. It should be noted that, although an arbitrary value in accordance with literature [52,53] is assigned for j, the attention was paid to keep it reasonable in comparison with radius of inclusion. Eventually, the objective function to be minimized, takes the following form.
where the double bracket refers to norm operation, the subscripts E and C stand for Eringen and Cosserat models, respectively, and the sub-numbers in parenthesis correspond to the material configurations (i.e. MAT 2, MAT 3, MAT 4). The numerical values of parameters used for different material configurations can be found in Table 2. Before any further progress, the procedure that is followed to attain nodal stresses should be explained. After obtaining the normal stress field of an element, it is a straightforward task to calculate the stresses at each node of corresponding element just by substituting their coordinates (i.e. direct evaluation method). Later, single nodal stress values are calculated by averaging the stress contribution of all the elements sharing the corresponding node. Note that, this approach is valid only if the nodes are surrounded by the elements possessing same kind of material properties. Hence, to obtain the correct stress values for interfacial nodes (e.g. point A), the contributions of matrix and inclusion elements are handled separately resulting in two distinct values.
It is assumed that, the optimization is achieved either the total number of stall iterations exceeds 50, or the objective function is minimized up to 0.1%. For this problem, the former condition is fulfilled and optimum fraction coefficient is obtained as: n opt ¼ 0:597873 with a convergence rate illustrated in Fig. 5. Once the nonlocal material parameters of Eringen's model are fixed, the results of finite element plate models employing either Eringen or Cosserat theories are analysed considering all scale ratios, material configurations, and domain discretization. Corresponding SCFs of both approaches are listed in Table 3, while the effect of discretization is studied through calculating the relative differences between both mesh configuration: Even though second-order elements clearly suppress the first-order ones in problems with high gradients (e.g. stress concentration), the capability of 'explicit' non-local model is studied only in terms of four-node elements. This restriction is a result of excessive computational burden of formation of crossstiffness matrices (i.e. k mn , k nm ) that emerges due to different Jacobians of quadrilateral elements. Hence, for comparison purposes, linear formulation is also implemented to 'implicit' non-local model. According to results listed for MAT 1, it is evident that the solutions are independent from the program used (i.e. COMSOL MultiphysicsÓ and Wolfram Mathematica). To see the global behaviour, the zoom-in and zoom-out looks of contour plots of normal stresses are illustrated at Figs. 6 and 7, together with their variation along vertical axis y (see Fig. 8). It is worth underlining that since FE method enforces continuity of displacement field only, the discontinuity of contour maps of the normal stresses among linear elements were required to be averaged by performing linear interpolation to achieve smooth and physically meaningful stress fields.
For all scale ratios following inferences can be made independent of the conducted non-local continuum theory:  • As the scale ratio decreases, SCF increases regardless of size effects, • Having a matrix with size effects leads to decrease in the SCFs (e.g. approximately 14% for MAT 2, MAT 4), which highlights the necessity of using a non-classical type theory in such a case, • The presence of internal length in inclusion seems to have negligible effect on the stress field of matrix region, • For matrix with size effects, the dependence of SCFs on spatial discretization is weakened compared to Cauchy type matrix. This situation is more evident for second order elements, which may be due to the following reason: Since the error of stresses are known to be minimum at GQ points, moving away from GQ points would introduce more numerical errors in quadratic elements compared to linear elements of same dimensions due to higher estimation of stress gradients.
The above-mentioned results are in accordance with the existing literature [32,48]. In the following subsections, the comments about each model with different scale ratios will be detailed with focusing on the conducted non-local theory.
3.1 Plate Model 1 (L/a = 3.0) As the first step, a plate with scale ratio of L/a = 3.0 is examined. In this particular case, the dimension of inclusion is comparable with the length of the plate which indicates the presence of a more definite nonlocality with respect to other scale ratios. According to Figs. 6, 7 and 8, two major differences are encountered between the results of 'explicit' and 'implicit' non-local models: • An increase of stress close to interface is attained for non-local type inclusion (i.e. Mat (3) and Mat (4)), • A decrease in stress close to boundaries of the domain is attained for non-local type matrix (i.e. Mat (2) and Mat (4)) when Eringen's nonlocal theory is employed.
Both phenomenon can be explained as a result of 'strong' non-local character of the Eringen model. For the former, the high strain field of matrix around interface, which falls into the influence zone of the inclusion part, results in a rise in stress field of inclusion; while for micropolar theory, such behaviour is not observed. In fact, for micropolar theory, the stress field of inclusion changes depending on the presence of size effects in matrix. One must bear in mind that if an imperfect interface is present, such an interaction between the two regions may weaken or even vanish for both non-local theories. For the latter, the reason of discrepancy seems to be the well-known boundary effect, which occurs due to missing contribution of non-existing neighbour elements. It is important to notice that boundary effects may be of importance for domains that have dimensions comparable with the influence zone of nonlocal parameter: In order to compensate the decreasing stresses, close to boundaries, the stress at remaining regions needs to be higher due to balance requirement. Indeed, for this particular scale ratio, SCFs of 'explicit' non-local model turns out to be higher than expected for MAT 2 and MAT 4. It is evident if one compares the difference between SCFs of 'implicit' and 'explicit' non-local models regarding all scale ratios: As illustrated in Fig. 8 In this section, plates with scale ratios of L/a = 10.0 and L/a = 20.0 are commented together since no significant difference in terms of stress fields are attained for Models 2 and 3. What is said for L/a = 3.0 about 'explicit' and 'implicit' non-local models holds with a slight adjustment: although boundary effects still exist, it has negligible effect around inclusion.
As for L/a = 20.0, the plate can be assumed as an infinite medium with a small inclusion, SCFs tabulated in Table 3  reported in the literature (Atroshchenko et al. [48]), by calculating their relative difference (i.e. at point A): Comparative results determined via Eq. 36 are listed in Table 4. The solutions of both non-local models appear to be in a good agreement with existing literature except for MAT 3: The relative error with respect to the analytical solution tends to increase for both non-local models with linear elements as the mesh is refined. This drawback is suppressed with nine-node elements, and the inherent advantage of using a second-order element formulation in case of strong material discontinuity is pointed out. Moreover, very similar situation is encountered for MAT 1, which suggests this drawback is stemmed from purely computational issues. Contrarily, for material configurations possessing a non-classical type matrix as in MAT 2 and MAT 4, the results tend to be independent from the spatial discretization (see Fig. 8) and tend to converge to explicit solutions for 'implicit' non-local model. Although this seems a bit different for 'explicit' non-local model, the slight increase in error is most likely due to the simplifications made for numerical calculations: number of Gauss sampling  (2), c Mat (3), d Mat (4) points and range of influence zone is selected such that the numerical integration with GQ technique provides less than %2 error with respect to numerical integration schemes embedded in Mathematica.
Consideration of larger number of GQ points and wider range of influence zone will reduce the errors, but will increase the computational burden dramatically. Note that 'implicit' model formulation does not require any more effort for numerical integration than local elasticity formulation.

Conclusion
In the present study, the elasticity problem of a plate with a circular inclusion under uniaxial tensile loading is investigated in terms of both 'implicit'/'weak' and 'explicit'/'strong' non-classical continuum models, Cosserat (micropolar) and Eringen's nonlocal theory [4,13,14], respectively, by employing the finite element method. For both theories, inclusion is incorporated as a heterogeneity possessing different material properties than matrix. The internal characteristic length of Cosserat model is assumed to be known while the non-locality of Eringen's model is tuned via optimizing the value of fraction coefficient in order to be in accordance with the analytical expressions of SCF reported for infinite Cosserat plates. Hence, Eringen's theory is not here used to model solids at atomistic scale, but rather the focus is on structures dominated by other kinds of discrete nature (e.g. heterogeneity) yielding an apparent nonlocal mechanical behaviour. In this way, apparently for the first time, a comparison between micropolar and Eringen's nonlocal theory of elasticity is performed, both in local and global manners, in the simplified context of two-dimensional inclusion-matrix problem that possess geometric singularities. According to the numerical simulations performed for different scale ratios and material configurations, optimization of fraction coefficient of 'explicit' description seems to provide a reasonable approximation to stress concentration factors obtained by 'implicit' theory, in which the peak stresses at the interface are reduced approximately about 14% for non-local type matrix configuration. Such significant decrease in SCFs highlights the necessity of utilizing a non-local type theory in the presence of size effects. On the other hand, major discrepancies are observed especially at the boundaries of non-local type inclusion and non-local type matrix. This discrepancy is attributed to 'strong' non-local character of Eringen model, that captures the interaction of particles within a neighbourhood which, also, paves the way to modelling the well-known boundary effect in atomic simulations. These observations provide an insight to possible advantages and drawbacks of both non-local descriptions, depending on the application. For instance, if this problem is a sub-model extracted from inner parts of a structure, boundary effects are unlikely to play a role, at least as far as the stress concentration factor is concerned, which justifies the use of 'implicit' description. On the other hand, to examine possible inclusions close to boundaries by sub-modelling, it is apparent that boundary effects may alter the stress concentration, and, in the end, the ultimate strength of the structure. In that case, 'explicit' model seems to be the better option. However, the computational burden introduced by 'explicit' modelling, both due to crossstiffness matrices and numerical integration, should not be underestimated. This naturally brings into consideration using a coupled non-local modelling in which the interior of the domain is modelled by 'implicit' approximation while the boundaries are modelled by 'explicit' description when necessary, by suitable displacement and traction conditions [66]. This may be the next step once a correspondence between two different non-local descriptions are put into evidence.
The numerical examples presented herein can be considered as a first evolution for matrix-inclusion problem, while more rigorous treatment may be required for a better approximation to the actual physical phenomena. For instance, in further studies inclusion and matrix parts can be considered as separate regions, with enforcing some contact constrains for interfacial relation [48] and by considering inhomogeneity as a source of additional attenuation effects [39]. In that case, the Euclidian distance appeared in kernel function of Eringen's nonlocal model cannot be taken equal to spatial distance, which will yield different results from those currently obtained. In any case, a connection between two different descriptions of non-locality has been demonstrated quantitatively, which is expected to pave the way to more enhanced modelling of the structures with size effects by unifying the superior aspects of these approaches both theoretically and computationally. To look for deeper correspondences, the study has been extended to problems with localised loads by focusing on displacement and stress fields [67] while further examination on crack problems is an ongoing project.