International Journal of Engineering Science 78 (2014) 124–133
Contents lists available at ScienceDirect
International Journal of Engineering Science
journal homepage: www.elsevier.com/locate/ijengsci
Interpolation–free discrete modeling with gradient matrix:
Case study of edge dislocation in linearly elastic crystal
Alexander Zisman ⇑
Institute of Applied Mathematics and Mechanics, St.-Petersburg State Polytechnic University, 29 Polytechnicheskaya Str., 195251 St.-Petersburg, Russia
a r t i c l e
i n f o
Article history:
Received 19 November 2013
Received in revised form 11 February 2014
Accepted 12 February 2014
Keywords:
Linear elasticity
Crystals
Dislocations
Gradient matrix
Modeling
a b s t r a c t
A meshless and interpolation-free (MIF) method of numerical modeling in solid mechanics
has been formulated with a gradient matrix extending the gradient operation to discrete
data. Nodal strains and, consequently, stresses in this method are expressed immediately
in terms of nodal displacements and the stress divergence in terms of nodal stresses that
makes it possible to get the stress balance equation in a truly discrete form. A trial MIF
model where nodal points correspond to atom positions is employed for a rectilinear edge
dislocation in a linearly elastic crystal. Both the resulting stress level at the dislocation core,
close to the theoretical crystal strength, and respective core dimensions prove to be realistic physically whereas calculated long-range stresses asymptotically approach the related
continuous ﬁelds known in an analytical form for the dislocation in linearly elastic isotropic continuum.
Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction
Interpolations are inherent in numerical models of continuum mechanics aimed to obtain approximate continuous solutions of constitutive equations. In particular, kinematics in the widely used ﬁnite element method (FEM) is formulated with a
priori shape functions expressing a velocity ﬁeld inside each ﬁnite element (FE) in terms of its nodal velocities (Gallager,
1975; Zienkievicz, 1970; Zienkievicz, Taylor, & Zhu, 2005). Such interpolations associated with the meshing intrinsic in
FEM result in some generic drawbacks of this method. On the one hand, the corresponding approximate ﬁelds have artiﬁcial
singularities at FE nodes, edges or faces reducing the solution accuracy. On the other hand, when physically prescribed discontinuities are considered FE faces (edges in 2D) should be placed on respective interfaces whereas nodal points must be
duplicated there. Even with such complications, however, the singularities of arbitrary shape functions will remain irrelevant
to actual discontinuities, which arise between material structural elements of rather complex shape. For example, it would
be hard to use shape functions ascribed to whole grains, which are polyhedrons with more than 20 apices. That is why, in
particular, FEM models of polycrystal deformation usually subdivide each grain into a lot of more primitive elements
(Kanjarla, Van Houtte, & Delannay, 2010; Roters et al., 2010) and become highly consumptive.
Apart from the above mentioned limitations, FEM models become too much awkward in application to large deformations, moving boundaries or propagating cracks, where the mesh should be repeatedly updated. In order to avoid re-meshing
troubles, meshless extensions of FEM have been developed where shape functions are associated with scattered nodal points
rather than certain domains between them and such points may be readily shifted, added or removed in simulation. Detailed
⇑ Tel.: +7 812 5527660; fax: +7 812 5577778.
E-mail address: [email protected]
http://dx.doi.org/10.1016/j.ijengsci.2014.02.015
0020-7225/Ó 2014 Elsevier Ltd. All rights reserved.
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
125
reviews of meshless models (Belyschko, Krongauz, Organ, Fleming, & Krysl, 1996; Chen, Lee, & Escandarian, 2006; Fries &
Matthies, 2004; Ullah & Augarde, 2013) are available and the author in no way claims to survey an extremely rich variety
of related works. However some geometric aspects of them signiﬁcant for the present paper should be addressed. The
moving least square (MLS) approach (Belytschko, Gu, & Lu, 1994; Nayroles, Touzot, & Villon, 1992) and the reproducing
kernel particle method (RKPM) (Liu, Jun, & Zhang, 1995; Liu, Li, & Belitschko, 1997), different in their mathematical origins
Lancaster and Salkauskas (1981) and Lucy (1977), respectively, and rather similar in implementation are most often
employed to derive approximate displacement ﬁelds. Such ﬁelds are usually treated with appropriately smooth shape
functions related to individual nodes, weight functions reﬂecting the local inﬂuence of these nodes, and partly overlaying
supports of the latter, which may be considered as domains of their inﬂuence. Along with the principal re-meshing problem
thus eliminated, other serious troubles inherent in FEM are also reduced, in particular artiﬁcial singularities between
neighboring FEs. Besides, properly modiﬁed meshless methods become superior to FEM in treating physically prescribed
discontinuities (Fries & Matthies, 2004). To this end nodal supports (normally disks in 2D and balls in 3D) are limited by
considered physical interfaces (Belytschko, Lu, & Gu, 1994) or discontinuous enrichment functions are added to continuous
approximations (Belytschko, Ventura, & Xu, 2002).
Despite the above considered advantages, meshless models still are not widely applied in the engineering routine because
they are much more computationally expensive than those of FEM owing to highly complicated treatment of employed
shape and weight functions (Belyschko et al., 1996; Chen et al., 2006; Fries & Matthies, 2004; Pan, Sze, & Zhang, 2004; Ullah
& Augarde, 2013). Furthermore, related approximations for displacement ﬁelds usually are not interpolations in a rigorous
sense because at nodal points they may differ from respective nodal variables. Such deviations become undesirable when
imposing essential boundary conditions and necessitate special treatments (Fernandez-Mendez & Huerta, 2004; Ullah &
Augarde, 2013). Another computational complication of meshless modeling is due to error estimation procedures (Chung
& Belytschko, 1998; Lee, Chung, & Choi, 2003) in the adaptive discretization essential in this numerical method. Meanwhile,
since both the latter and FEM have speciﬁc limitations, efforts have been undertaken to ﬁnd a reasonable compromise between them. Thus in a coupled FE/meshless method (Ullah, Combs, & Augarde, 2013) the whole model is ﬁrst treated with
FEM and then some domains are converted into meshless zones if a predeﬁned error measure is violated. In the extended
FEM (XFEM) reviewed in Belytschko, Gracie, and Ventura (2009), where certain formal grounds of the meshless method
are used and known properties of physical discontinuities are embodied in special enrichment functions, cracks, dislocations
and boundaries considered with a ﬁxed mesh may cross FE volumes. Smoothed FEM (SFEM), where conventional meshbased shape functions have been smoothed within FE (Liu, Nguyen, Dai, & Lam, 2007), at their nodes (Nguyen-Thoi, Vu-Do,
Rabczuk, & Nguyen-Xuan, 2010), edges (Liu, Nguyen-Thoi, & Lam, 2009) or faces (Nguyen-Thoi, Liu, Lam, & Zhang, 2009), are
also cited as competitors to the meshless method, although they retain the re-meshing problem.
Computational drawbacks related to both the mesh-based and meshless shape functions could be avoided if numerical
methods were truly discrete that is not only meshless but also interpolation-free, and the present paper is an effort in this
direction. The author has limited himself to a rather brief overview of the conventional methods, mentioned above, just because of the indicated principal distinction of the proposed approach. An essential prerequisite to develop the meshless and
interpolation-free (MIF) approach in solid mechanics is the increasing computational power that makes it possible to process
very high densities of nodal points and, hence, to image and to consider simulation results in terms of nodal variables only.
Besides, MIF models would appear to be particularly applicable to physically discrete (e.g. crystalline) structures having no
stiff matter in interatomic space. Indeed, any assumption of strain distribution over this space would be completely speculative, whereas geometrically small deviations of similar structures from equilibrium states still should be linearly elastic.
Therefore a mechanistic approach remains reasonable since the strain and, hence, stress tensor can be ascribed to a rather
small number of adjacent atoms (not less than four in 3D or three in 2D) as if they were nodal points of a virtual continuum
fragment. At the same time, although the density of discretization is rapidly increasing, the current mesh-based and meshless methods could not plainly get rid of continuous approximations because the latter are used in calculating stiffness terms
and, eventually, nodal forces and stresses. Therefore a challenging question arises of how to allow for the material stiffness
with no suggestion for continuous ﬁelds between scattered nodal points.
In the present paper, aimed to develop a MIF method of numerical modeling, an elastic response with the material elastic
modulus is considered exclusively at nodal points, whereas an underlying local strain at each node is derived depending on
displacements of neighboring nodes. The predetermined number of thus involved neighbors depends on computational and
physical reasons. In order to implement such a method, one should express constitutive equations of solid mechanics in a
truly discrete form that is ﬁrst of all to extend the conventional gradient operator directly to nodal variables. An appropriate
extension proposed by Zisman and Ermakova (2006) has been called gradient matrix and applied to the strain mapping
(Zisman, Ivanov, Lomov, & Verpoest, 2006) and to the imaging of crystal curvature (Zisman, Van Boxel, Seefeldt, & Van
Houtte, 2008) later on. This matrix is uniquely expressed in terms of nodal coordinates contained in whatever differentiation
domain and derives from the corresponding data sampling the uniform gradient component, thus extracting only a linear
part of an underlying continuous ﬁeld and neglecting unknown residual errors. Similar residuals, like inaccuracy of arbitrary
shape functions, are inherent in discrete models. In any case, the denser situated nodal points, the smaller such errors. In the
present paper the gradient matrix will be employed in order to express nodal strains and hence stresses immediately in
terms of nodal displacements. Furthermore, to properly extend the stress balance equation and thus to enable the discrete
MIF modeling in solid mechanics, the divergence operation will also be expressed in a truly discrete form.
126
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
Another effort to introduce a discrete gradient in numerical methods of solid mechanics has recently been undertaken by
Lu, Qiun, and Han (2008). In this formulation a Voronoi cell is ascribed to each nodal point in order to get a truly meshless
model. However, the Voronoi tessellation remains in a sense a kind of meshing that inherits artiﬁcial limitations of FEM.
Thus, following Lu et al. (2008), one has to assume uniformity of elastic strain inside each cell in calculating stiffness terms,
whereas the strictly prescribed selection of neighboring nodal points disables a potentially useful option to vary the shape
and dimensions of local differentiation domains.
For a case study the present paper considers the stress ﬁeld of a rectilinear edge dislocation situated in a crystal of ﬁnite
cross section. Such a speciﬁc application is quite appropriate to test potentiality of the discrete MIF method in simulating
highly singular continuous ﬁelds. Besides, it would not be practicable to apply the proposed new method to usual boundary
value problems readily tractable with traditional numerical methods. Meanwhile the latter are less adapted to simulation of
internal defects (dislocations, inclusions, etc.), which are inherent in real solids, and it is also advisable to estimate accuracy
of mechanical MIF models in application to physically discrete microscopic structures. Accordingly, nodal points in the considered MIF model have been associated with atom positions in the dislocated crystal. It should be underlined that the
author in no way claims to discover new properties of the lattice dislocation; on the contrary, such a thoroughly examined
defect has been wittingly selected for a trial application in order to verify the proposed method of discrete modeling.
The paper is organized as follows. First, the gradient matrix is deﬁned and its application to nodal variables (displacements, strains, and stresses) is described. Then a procedure of MIF modeling is explained and a trial model of the dislocated
crystal is introduced. Next, the calculation results are reported and compared to the continuous long-range stress ﬁeld of the
dislocation, known in an analytical form, and to the relevant physical data on the dislocation core. In conclusion, potentiality
of the MIF approach, untapped in this case study, is discussed that may apply to both crystal defects and continuous models.
2. Theory
2.1. Gradient matrix
Let N arbitrary vectors u(k) be given over a cluster of respective nodal points r(k) (k = 1, 2, . . ., N) contained in a differentiation domain. The sought tensor (3 3 matrix)
L ¼ $u
ð1Þ
is presumed to be uniform over this domain and to generate a linear ﬁeld u(r) that ﬁts the given discrete data in the best
way. It is convenient to represent the considered variables and points in a Cartesian coordinate system by 3 N matrices
0
ð1Þ
ux
B
U ¼ @ uð1Þ
y
ð1Þ
uz
0
ð2Þ
ux
ðNÞ
::: ux
1
uy
ðNÞ C
::: uy A;
ð2Þ
uz
::: uz
ð2Þ
ð2Þ
ðNÞ
xð1Þ
B ð1Þ
R ¼ @y
xð2Þ
yð2Þ
1
::: xðNÞ
C
::: yðNÞ A
zð1Þ
zð2Þ
:::
ð3Þ
zðNÞ
or, in case of 2D formulation, by their 2 N counterparts. In these terms the desired ideal ﬁt would mean that U is expressed
as
U ¼ U o þ LT R;
ð4Þ
where each column of matrix Uo represents a unique vector uo, related to the coordinate origin, and the superscript T indicates the matrix transposition. However, such an equality is not always possible, whichever L. In 3D formulation Eq. (4) has a
unique solution only in case of four nodal points (N = 4) including the coordinate origin; at N < 4 L becomes uncertain and at
N > 4 no solution exists. Similarly, in 2D case, where matrices U and R have 2 N dimensions and matrix L has 2 2 dimensions, a unique solution exists only at N = 3. Meanwhile notably greater N is often required in order to smooth over discrete
data or to reduce noise of them. Similar difﬁculties are inevitable in discrete models of continuum mechanics where the
number of immediate neighbors around each nodal point is normally more than four in 3D and more than three in 2D. Thus,
the problem becomes overdetermined and can be resolved only with some approximation. Apparently, the best approximation for L should provide the corresponding nodal variables U⁄ least deviated from U in terms of the mean square deviation.
A simple algebraic method to ﬁnd the desired L with no straightforward minimization is based on the gradient matrix
(Zisman & Ermakova, 2006; Zisman et al., 2006, 2008).
½$ ¼ RT ðRRT Þ
1
ð5Þ
that is a sort of Moore–Penrose pseudo-inverse (Moore, 1920; Penrose, 1955) as applied to our coordinate matrix R. The
gradient matrix has N 3 or N 2 dimensions in 3D or 2D formulations, respectively, and is readily recalculated if
nodal points involved in the differentiation domain are re-arranged. As shown in Zisman and Ermakova (2006), existence
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
127
of (RRT)1 and, hence, of [$] in 3D or 2D is ensured only if these points are not all coplanar or collinear, respectively. Besides,
in order to ﬁlter off a uniform component of treated variables, the gradient matrix should be determined relative to the coordinate origin situated in the geometric center of all involved data points so that
N
X
rðiÞ ¼ 0:
ð6Þ
i¼1
Introduced as a discrete counterpart of the conventional gradient, [$] suggests a matrix expression
L ¼ ðU½$ÞT
ð7Þ
that resembles Eq. (1) and, as shown in Zisman and Ermakova (2006), Zisman et al. (2006), provides the best ﬁt of the
product
U ¼ LT R
ð8Þ
to the given U that is
X
T
jjU U jj2 ¼ Sp ðU U ÞðU U Þ ¼
ðuk uk Þ2 ¼ min
ð9Þ
k
with whatever uniform part of U. Vectors uk in statement (9) correspond to respective columns of matrix U⁄ like uk of U.
When comparing the MIF concept to the current meshless method (Belyschko et al., 1996; Chen et al., 2006; Fries & Matthies,
2004; Ullah & Augarde, 2013) where continuous ﬁelds between nodal points are considered, it is notable that the discrete
gradient matrix implies a jumping rather than a moving least square (MLS) procedure.
It is interesting that expression (5) for the gradient matrix that has both needful matrix dimensions and an appropriate
physical dimension (inverse length) was found by intuition ﬁrst, whereas uniqueness and the principal property of respective U⁄ to satisfy statement (9) have been proved afterwards with the help of algebraic properties of orthogonal projections
(Zisman & Ermakova, 2006; Zisman et al., 2006). At the same time, although [$] was particularly introduced in order to extend the spatial gradient to discrete data, it represents a special case of the general pseudo-inverse formulation by Moore
(1920) and Penrose (1955).
The divergence operation, based on the gradient matrix and applicable to discrete data, should also be introduced in order
to properly extend the stress balance equation in what follows. With the identity
div u ¼ $u ¼ Spð$uÞ;
ð10Þ
where the right part is the trace of tensor ru, and Eqs. (1) and (7) kept in mind, the divergence of vector ﬁeld u(r) is expressed in terms of respective nodal variables as
div u ¼ SpðU½$Þ:
ð11Þ
2.2. Nodal strains, stresses and forces by gradient matrix
Let R and [$] be ascribed to a domain containing the considered nodal point together with a number of neighboring ones
and U be the corresponding matrix of geometrically small elastic displacements. Hence the elastic strain, as follows from Eqs.
(1) and (7), is expressed as
e ¼ ð$u þ u$Þ=2 ¼ ððU½$ÞT þ U½$Þ=2;
ð12Þ
and the related nodal stress may be found with Hooke’s law.
The stress balance equation
$r þ f ¼ 0
ð13Þ
allows for both the stress ﬁeld r and the body forces f, although purely mechanical applications often presume f = 0 inside considered solids, as is the case in the present paper. On the other hand, computational solid mechanics deals with
approximate and, hence, somewhat unbalanced stress ﬁelds. Therefore to restore the balance according to Eq. (13) is to
introduce ﬁctitious f – 0 depending on inaccuracy of the calculated stress ﬁeld. Consequently, in order to get rid of ﬁctitious
terms in the ﬁnal solution, one has to compensate f by opposite virtual forces
f u ¼ f
ð14Þ
and then to ﬁnd the stress ﬁeld correction due to fu. A similar approach, employed long ago to introduce defects in continuum (Eshelby, 1957; Hirth & Lothe, 1968; Read, 1953), could be used with discrete MIF models as well. However, Eq. (11),
previously formulated for vectors, cannot immediately apply to r that is a tensor of the second order. To get over the difﬁculty, let us express the stress tensor in a dyadic form
r ¼ S1 i þ S2 j þ S3 k
ð15Þ
128
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
where i, j, and k are unit vectors of a macroscopic Cartesian basis. Then, with expression (11) kept in mind, the stress
divergence
$ r ¼ ið$ S 1 Þ þ jð$ S 2 Þ þ kð$ S 3 Þ
ð16Þ
takes on form
$r ¼ iSpð½S 1 ½$Þ þ jSpð½S 2 ½$Þ þ kSpð½S 3 ½$Þ;
ð17Þ
U and contain respective stress vectors S1(k),
F(m), whereas matrices [S1](m), [S2](m), [S3](m)
where matrices [S1], [S2], and [S3] have the same dimensions as
and S3(k)
(k = 1, 2, . . ., N). Let a ﬁctitious force at nodal point number m be
and [$](m) be
derived over the corresponding differentiation domain. Then the desired extension of Eq. (13) is expressed with Eq. (17)
in terms of nodal variables only:
iSpð½S 1 ðmÞ ½$ðmÞ Þ þ jSpð½S 2 ðmÞ ½$ðmÞ Þ þ kSpð½S 3 ðmÞ ½$ðmÞ Þ þ F ðmÞ ¼ 0;
Fu(m)
S2(k)
ð18Þ
(m)
and the corresponding
= F
is calculated according to Eq. (14). Note that Eqs. (16)–(18) are written in 3D formulation; when applied to 2D models, similar equations would contain neither unit vector k nor subscript 3.
The kinematics of MIF method, as formulated above, may apply to materials of any rheological properties. However, for
deﬁniteness and simplicity sake, let us address Hooke’s law for elastically isotropic matter in the plane strain case:
rxx =G ¼
2ð1 mÞ
2m
exx þ
eyy
1 2m
1 2m
ryy =G ¼
2m
2ð1 mÞ
exx þ
eyy
1 2m
1 2m
ð19Þ
rxy =G ¼ ryx =G ¼ 2exy ¼ 2eyx ;
where G is the shear modulus and m is the Poisson ratio. Therefore all terms containing k in Eqs. (16)–(18) can be omitted and
the normalized vectors, corresponding to stress expression (15), are respectively equal to
S 1 =G ¼
2ð1 mÞ
2m
exx þ
eyy i þ 2eyx j;
1 2m
1 2m
ð20Þ
S 2 =G ¼
2m
2ð1 mÞ
exx þ
eyy j þ 2exy i:
1 2m
1 2m
Hence the matrices for [S1] and [S2] in Eq. (17) are expressed as
½S 1 ¼
2G
1 2m
2G
½S 2 ¼
1 2m
ð1Þ
ð1Þ
ð1 mÞexx þ meyy
ðNÞ
ðNÞ
. . . ð1 mÞexx þ meyy
ð1Þ
...
ð1 2mÞeyx
ð1Þ
...
ð1 2mÞexy
ð1 2mÞeyx
ð1 2mÞexy
!
;
ðNÞ
ðNÞ
ð1Þ
ðNÞ
ðNÞ
með1Þ
. . . mexx þ ð1 mÞeyy
xx þ ð1 mÞeyy
ð21Þ
!
:
2.3. Intermediate conﬁguration for internal stress source in MIF model
Internal stress sources (defects) cannot be treated in the context of a usual boundary value problem where a continuous
ﬁeld of displacements may gradually transform a reference stress-free conﬁguration {R0} into the sought ﬁnal conﬁguration
{RE}, corresponding to external boundary conditions and satisfying the stress balance equation. In particular, to generate a
dislocation together with related elastic ﬁelds in the model interiors is, by deﬁnition, to introduce a characteristic displacement discontinuity (Burgers vector) at any contour drawn around the defect line. That is why in the continuum theory of defects the dislocation formation is treated as a particular sequence of virtual operations (Hirth & Lothe, 1968; Read, 1953): (a)
continuum slitting terminated at the given line, (b) introduction of certain displacement jumps across the cut surface, (c)
restoration of continuity, and (d) elastic relaxation. The last but one operation, as shown by Eshelby (1957), is equivalent
to an appropriate distribution of body forces on the ‘‘cut surface’’ where the deformation compatibility has been violated.
Following this approach unlike traditional numerical methods, we will consider not only the reference conﬁguration {R0}
of the solid model but also its intermediate conﬁguration {RI} with the virtual nodal forces corresponding to the defect insertion. Such forces, in turn, should be reduced with the model relaxation leading to the sought equilibrium conﬁguration {RE}.
Let a differentiation domain and the corresponding coordinate matrix R0 be determined for each nodal point of {R0}. Then
{RI}, roughly similar to the expected ﬁnal {RE} and generally containing virtual nodal forces, should be intuitively selected in
order to introduce the considered defect. Comparison of {RI} to {R0} is not trivial, however, because the latter cannot be entirely converted into the former by a compatible deformation ﬁeld and, hence, there is no point-to-point correspondence
between these conﬁgurations. Nevertheless it remains possible and even apparent on a local scale if two respective domains
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
129
have the same number and similar arrangements of nodal points. In particular, any deformed cluster of N neighboring atoms
outside the dislocation core may be readily related to an appropriate selection of N nodes of the perfect crystal lattice. With
thus determined local matrices RI and respective R0 the corresponding displacement matrices are expressed as
U I ¼ RI R0
ð22Þ
and can be differentiated with the gradient matrices
1
½$0 ¼ RT0 R0 RT0
ð23Þ
found as [$] in expression (5). The nodal strains and stresses are then derived following Eq. (12) and Hooke’s law, respectively. When calculating the related stress divergence with Eqs. (17) and (21) or its 3D counterpart, an updated gradient
matrix
½$I ¼ RTI ðRI RTI Þ
1
ð24Þ
F ðmÞ
u
must be used instead of [$]0 since the previously obtained stresses are considered over {RI}. All virtual nodal forces
in
{RI} are then calculated with Eqs. (14) and (18) and should be employed to relax this intermediate conﬁguration as described
in the next section. It is worth noting that displacements UI in a priori selected intermediate conﬁgurations are not always
geometrically small and hence the linear elastic response assumed above may be a somewhat rough approximation.
However, as considered in the discussion section, the following iteration procedure should reduce this inaccuracy.
2.4. Iteration procedure in MIF modeling
An elastic stiffness matrix predetermined over a whole numerical model routinely reduces the boundary value problem to
simultaneous linear equations (Gallager, 1975; Zienkievicz, 1970; Zienkievicz et al., 2005). In principle, when a MIF model is
defect-free and hence no intermediate conﬁguration is needful, a similar approach is possible. Indeed, the stiffness matrix mentioned above can be found, as evident from Section 2.2, with calculation of nodal forces due to small trial displacements in {Ro}
applied to all nodal points by turns. However, the MIF modeling will cost undertaken efforts only in application to problems,
which could not be treated by conventional methods. This is the case, in particular, if any stress source (defect) is introduced
into the considered matter with some intermediate conﬁguration {RI} that violates the stress balance condition. Then an appropriate procedure should be employed to implement a transition from {RI} to the ﬁnal equilibrium conﬁguration {RE}.
Let RJ, UJ and FJ be the nodal coordinates, displacements and virtual forces, respectively, in iteration number J, whereas
numeration of nodal points is omitted here for brevity. In order to gradually restore the stress balance, nodal displacement
corrections
dU Jþ1 ¼ aF J
ð25Þ
are repeatedly calculated, where J = 0 corresponds to the intermediate conﬁguration, and a has a proper physical dimension
[u/f] and a magnitude that keeps the procedure converging. Following Section 2.2, the corresponding corrections drJ+1 of nodal stresses are found with dUJ+1 and the gradient matrices
1
½$J ¼ RTJ ðRJ RTJ Þ ;
ð26Þ
which are then updated together with respective coordinate matrices:
RJþ1 ¼ RJ þ dU Jþ1 ;
½$Jþ1 ¼ RTJþ1 ðRJþ1 RTJþ1 Þ
1
ð27Þ
Thereafter, with [$]J+1 and the previously found drJ+1, the corrections dFJ+1 are derived like nodal forces in Section 2.2 and the
resulting forces are calculated for iteration number J + 1:
F Jþ1 ¼ F J þ dF Jþ1
ð28Þ
Similar iterations are repeated until magnitudes of virtual nodal forces become negligible with a stated tolerance. The stress
at any nodal point in iteration number M, in particular the last iteration, can be found with the respective stress r0 calculated
in the intermediate conﬁguration {RI} and the stress corrections drJ calculated in all following iterations (J = 1, 2, . . ., M):
rM ¼ r0 þ
M
X
drJ :
ð29Þ
J¼1
3. Trial discrete model: Linearly elastic dislocated crystal
3.1. Intermediate conﬁguration and reduction of nodal forces
Under consideration is a rectilinear edge dislocation in the primitive cubic lattice. The model crystal is inﬁnite along the
dislocation line, while its upper and lower widths are 101b and 100b, respectively, where b is the Burgers vector magnitude,
130
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
and the crystal height is 102b. A representative part of the intermediate conﬁguration {RI}, i.e. an initial arrangement of
atoms treated for nodal points, is schematically shown in Fig. 1a. Although there is no stiff matter between depicted nodes
of the crystal lattice, they in no way lose sight of elasticity of the corresponding continuum model since its stiffness is not
something other than stiffness of elementary cells of the crystal. Following Section 2.2, the isotropic elasticity is assumed
here, whereas all geometric terms correspond to an orthotropic cubic lattice. Such a simpliﬁcation, often used in the materials science, is not inherent in the proposed approach but is employed only in order to facilitate the trial MIF simulation. It is
worth noting, as well, that connecting lines in Fig. 1 do not implicate any meshing in the simulation but are merely
illustrative.
Local coordinate matrices RI in {RI} and the related gradient matrices [$]I employed to differentiate discrete data with
expression (7), involve only the nearest neighbors of each node and take into account local distortions relative to the perfect
lattice (stress-free conﬁguration {R0}), as illustrated in Fig. 1b. For brevity sake, respective layouts at four crystal corners,
each consisting of three nodes, are omitted here. The nodal forces FI have been calculated with [$]I following Section 2.3.
These forces, corresponding to the simple shearing of elementary cells, prove to localize in the two planes crossing the dislocation core, as illustrated in Fig. 1a. The local coordinate matrices, related gradient matrices and nodal forces for the next
iterations are then repeatedly re-calculated, following Section 2.4.
Distribution of virtual nodal forces in the model crystal, calculated with m = 1/3 by the MIF method according to Sections
2.3 and 2.4, are mapped in Fig. 2a and b for the starting ‘‘intermediate’’ state and after 500 iterations, respectively. As one can
see, during the iteration procedure such forces spread over the crystal and diminished by three orders. Considering such a
relaxation degree satisfactory, let us address in the next section the corresponding elastic ﬁelds, generated by the dislocation.
3.2. Calculation results
A map of the model shear stresses rxy shown in Fig. 3a satisfactorily complies with the long-range stress ﬁeld of a straight
edge dislocation in linearly elastic isotropic continuum derived with the known expression (Hirth & Lothe, 1968; Read,
1953). Besides, the corresponding elastic displacements allow one to visualize the lateral proﬁle of the model crystal (right
side) containing the dislocation, as shown in Fig. 3b. When neglecting minor irregularities at the middle plane due to somewhat insufﬁcient relaxation, this result also proves to be quite realistic.
Fig. 1. (a) Intermediate conﬁguration of the dislocated crystal with virtual nodal forces and (b) layouts of nodal points contained in employed
differentiation domains.
Fig. 2. Map of virtual nodal forces in the dislocated crystal: (a) the intermediate conﬁguration to be relaxed and (b) the conﬁguration relaxed by 500
iterations of MIF procedure.
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
131
Fig. 3. (a) Map of calculated shear stresses over the dislocated crystal, and (b) a lateral proﬁle of the crystal in the intermediate conﬁguration (light circles)
and after 500 iterations of MIF procedure (dark circles).
Fig. 4. Shear (a) and tensile (b) stresses in the horizontal and vertical planes, respectively, crossing the dislocation core. The dotted and solid lines
correspond to the MIF modeling (500 iterations) of the dislocated crystal and to the respective classical expressions for inﬁnite continuum.
In order to analyze the modeling results in the most singular area around the dislocation core, it is expedient to consider
the stress distribution in crystal planes crossing the core. The shear stress rxy in the horizontal crystal plane touching the
core top (y/b = 1) is plotted in Fig. 4a. The plot suggests the core width of 4b and the crystal shear strength of about 0.1G.
The tensile stress rxx in the left vertical half plane under the core (x/b = 0.5, y/b 6 0) and in the middle vertical half plane
over the core (x/b = 0, y/b > 0), plotted in Fig. 4b, corresponds to the core height of 2b and the crystal tensile strength of about
0.22G. All these magnitudes satisfactorily comply with the relevant physical estimates (Dieter, 1986; Hirth & Lothe, 1968).
It is worth noting that the calculated rxx component is slightly less in magnitude than its continuum counterpart (Fig. 4b).
This is not at all surprising since the considered crystal is ﬁnite in plane XY and, hence, its local bending due to the dislocation
is less constrained in comparison to the case of inﬁnite matter. Apparently, stronger constraints (bigger crystal) would augment magnitudes of both compressive stresses above the dislocation (Y > 0) and tensile stresses under the latter (Y < 0).
4. Discussion
In the present paper a truly discrete expression (18) for the stress balance equation has been based on the simplest form
of gradient matrix assuming uniformity of the sought value of ru over a differentiation domain. When the latter covers only
one coordination shell of nodes surrounding the central one and almost equidistant from the latter, as is the case in this
work, such an assumption is relevant to discretization; otherwise possible non-linearity of an underlying ﬁeld u(r) will result
in notably more inaccuracy. On the other hand, in case of highly non-uniform distribution of nodal points, inherent in certain
problems, local gradient calculations may involve various quantities of nodal variables. Hence some reﬁnement of the present MIF formulation would be advisable, speciﬁcally allowing for nodes differently remote from a center of their cluster considered. In this regard it is expedient to partition a unique matrix expressions (7) into individual nodal contributions
(i = 1, 2, . . ., N) according to (Zisman & Ermakova, 2006; Zisman et al., 2006):
132
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
N
X
$u ¼
g i ui ;
N
X
g i ¼ ri ri ri
i¼1
!1
ð30Þ
i¼1
where gi is a partial gradient term related to node number i and ui is displacement of the latter. Therefore a weighted evaluation of ru may be expressed as
ð$uÞx ¼
N
X
ðg i ui Þxðr i Þ;
ð31Þ
i¼1
where ri = (riri)1/2 and x(r) is an appropriate weight function decreasing with distance r. At the same time, unlike the current
meshless method (Belyschko et al., 1996; Chen et al., 2006; Fries & Matthies, 2004; Ullah & Augarde, 2013), truly discrete MIF
models should use only nodal values xi = x(ri), whichever forms and properties of employed weight functions. The discussed
reﬁnement, however, is a subject of further works as this would complicate the present introductory paper.
Meanwhile a successful application of MIF modeling in the case study of edge dislocation evidences for a sound potential
of this method in solid mechanics, particularly at its frontier related to microscopic structures. However, a question may
arise of why the assumed linear elasticity properly works on a microscopic scale around the dislocation core. The reason
is that most strong crystal distortions in the selected intermediate conﬁguration (middle layer in Fig. 1a) do matter only during some number of starting iterations and are gradually reduced as the model approaches the balance. Therefore, if the procedure eventually converges as is the case here, elastic strains in the ﬁnal conﬁguration become much smaller and a linear
response to them becomes satisfactory approximation, except for ‘‘interiors’’ of the microscopic core, of course. As to the
assumption of elastic isotropy, it was employed in this trial MIF model for simplicity sake only, whereas Hooke’s law in
its general (anisotropic) form may be used with the same geometrical formalism.
Estimates for the dislocation core dimensions and for the related crystal strength, although they have been made with a
mechanistic model, are in a surprisingly good correspondence to relevant physical data. Apart from the considered modeling
method, such results should be attributed to the physically motivated selection of nodal points. Indeed, the considered lattice dislocation was treated on the scale of its characteristic Burgers vector, determining related elastic ﬁelds. This scale also
corresponds to the least structural length where mechanical response of crystalline materials can be rationally considered.
Based on similar reasons, other numerical models, also essentially mechanistic, have been applied to microscopic structures
of polymers (Wang et al., 2006) and to structured interfaces in living tissues (Bertoldi, Bigoni, & Drugan, 2007). In the present
work the minimum number of adjacent nodal points covered by a differentiation domain around each considered node
seemingly implicates from the physical point of view the minimum cutoff radius of interatomic interaction. Hence it would
be interesting in further MIF works to vary the domain dimensions and, accordingly, the number of involved neighbors in
order to estimate effects of crystal type and of atomic bonding on the shape and dimensions of the dislocation core. Furthermore, similar efforts would concern foundations of non-local elasticity approach in solid mechanics (Eringen, 2002; Kroner,
1967).
As evident from the reported case study, most relevant applications of MIF models should reﬂect material structures. In
particular, a known way to eliminate nonphysical strain singularities, while treating crystal defects with continuous models,
is the gradient elasticity approach assuming certain ‘‘internal lengths’’ and introducing appropriate constitutive equations
(Gutkin & Aifantis, 1999). To justify both these assumptions and equations, the crystal structure should be taken into account
more explicitly, and the discrete MIF modeling is a good expedient for that. Besides, when considering damage processes in
speciﬁcally complicated structures like textile composites, it is not always possible to indicate distinctive cracks (fracture
surface) on local scales. Hence appropriate damage characteristics, expressed in terms of strain or stress, should be ascribed
to representative material domains (Lomov et al., 2008), and the gradient matrix ﬂexibility in variation of the treated domain
dimensions appears to be very favorable in this connection. The same option may facilitate numerical analysis of strong
stress singularities arising in deformed heterogeneous materials around corners and apices of faceted hard particles, etc.
5. Conclusions
The gradient matrix used in this paper is a ﬂexible form of truly discrete derivative operator that applies immediately to
nodal variables with no assumption for continuous ﬁelds between scattered nodal points. This matrix allows one to simply
extend constitutive differential equations to discrete data and hence enables the meshless and interpolation-free (MIF)
numerical modeling.
When applied to displacements, strains and stresses, the MIF modeling with an appropriate iteration procedure makes it
possible to solve the stress balance equation in terms of nodal variables only. For a trial application of this approach in solid
mechanics, a rectilinear edge dislocation in the primitive cubic lattice has been simulated and the obtained results prove to
be in a good correspondence to relevant physical data on the defect core and to the long-range stress ﬁeld of respective dislocation in continuum, known in an analytical form.
Acknowledgement
The author gratefully acknowledges stimulating discussions of this work with C. Teodosiu and P. Van Houtte.
A. Zisman / International Journal of Engineering Science 78 (2014) 124–133
133
References
Belyschko, T., Krongauz, Y., Organ, D., Fleming, M., & Krysl, P. (1996). Meshless methods: an overview and recent developments. Computer Methods in Applied
Mechanics Engineering, 139, 3–47.
Belytschko, T., Gracie, R., & Ventura, G. (2009). A review of extended/genegalized ﬁnite element method (SFEM). Modelling and Simulations in Materials
Science and Engineering, 17, 24. 043001.
Belytschko, T., Gu, L., & Lu, Y. Y. (1994). Fracture and crack growth by element-free Galerkin method. Modelling and Simulations in Materials Science and
Engineering, 2, 519–534.
Belytschko, T., Lu, Y. Y., & Gu, L. (1994). Element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 37, 229–256.
Belytschko, T., Ventura, G., & Xu, J. X. (2002). New methods for discontinuity and crack modeling in EFG. In M. Griebel & M. A. Shweitzer (Eds.). Meshfree
Methods for Partial Differential Equations (26, pp. 37–50). Berlin: Springer Verlag.
Bertoldi, K., Bigoni, D., & Drugan, W. J. (2007). Structural interfaces in linear elasticity. Part I: nonlocality and gradient approximations. Journal of the
Mechanics and Physics of Solids, 55, 1–34.
Chen, Y., Lee, J., & Escandarian, A. (2006). Meshless methods in solid mechanics. New York: Springer.
Chung, H. J., & Belytschko, T. (1998). An error estimate in the EFG method. Computational Mechanics, 21, 91–100.
Dieter, G. E. (1986). Mechanical metallurgy. New York: McGraw-Hill.
Eringen, A. C. (2002). Nonlocal continuum ﬁeld theories. New York: Springer.
Eshelby, J. D. (1957). The determination of the elastic ﬁeld of an ellipsoidal inclusion and related problems. Proceedings of the Royal Society of London, A241.
476-396.
Fernandez-Mendez, S., & Huerta, A. (2004). Imposing essential boundary condition in mesh-free methods. Computer Methods in Applied Mechanics and
Engineering, 193, 1257–1275.
Fries, T. P., & Matthies, H. G. (2004). Classiﬁcation and overview of meshfree methods. Braunsweig: Technical University Braunsweig.
Gallager, R. H. (1975). Finite element analysis. New Jersey: Prentice-Hall.
Gutkin, M. Yu., & Aifantis, E. C. (1999). Dislocations in the theory of gradient elasticity. Scripta Materialia, 40, 559–566.
Hirth, J. P., & Lothe, J. (1968). Theory of dislocations. New York: McGraw-Hill.
Kanjarla, A. K., Van Houtte, P., & Delannay, L. (2010). Assessment of plastic heterogeneity in grain interaction models using crystal plasticity ﬁnite element
method. International Journal of Plasticity, 26, 1220–1233.
Kroner, E. (1967). Elasticity theory of materials with long-range cohesive forces. International Journal of Solids and Structures, 3, 731–742.
Lancaster, P., & Salkauskas, K. (1981). Surfaces generated by moving least squares methods. Mathematics of Computation, 37, 141–158.
Lee, G. H., Chung, H. J., & Choi, C. K. (2003). Adaptive crack propagation analysis with the elementfree Galerkin method. International Journal for Numerical
Methods in Engineering, 56, 331–350.
Liu, W. K., Jun, S., & Zhang, Y. F. (1995). Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids, 20, 1081–1106.
Liu, W. K., Li, S., & Belitschko, T. (1997). Moving least square reproducing kernel methods: I. methodology and convergence. Computer Methods in Applied
Mechanics and Engineering, 143, 113–154.
Liu, G. R., Nguyen, T. T., Dai, KY., & Lam, K. Y. (2007). Theoretical aspects of the smoothed ﬁnite element method (SFEM). International Journal for Numerical
Methods in Engineering, 71, 902–930.
Liu, G. R., Nguyen-Thoi, T., & Lam, K. Y. (2009). An edge-based smoothed ﬁnite element method (ES-FEM) for static, free and forced vibration analyses of
solids. Journal of Sound and Vibration, 320, 1100–1130.
Lomov, S. V., Ivanov, D. S., Truong, T. C., Verpoest, I., Baudry, F., Van den Bosche, K., & Xie, H. (2008). Experimental methodology of study of damage initiation
and development in textile composites in uniaxial tensile test. Composites Science and Technology, 68, 2340–2349.
Lu, J., Qiun, J., & Han, W. (2008). Discrete gradient method in solid mechanics. International Journal for Numerical Methods in Engineering, 74, 619–641.
Lucy, L. B. (1977). A numerical approach to the testing of the ﬁssion thesis. Astronomical Journal, 82, 1013–1024.
Moore, E. H. (1920). On the reciprocal of the general algebraic matrix. Bulletin of the American Mathematical Society, 26, 394–395.
Nayroles, B., Touzot, G., & Villon, P. (1992). Generalizing the ﬁnite element method: diffuse approximation and diffuse elements. Computational Mechanics,
10, 307–318.
Nguyen-Thoi, T., Liu, G. R., Lam, K. Y., & Zhang, G. Y. (2009). A face-based smoothed ﬁnite element method (FS-FEM) for 3d linear and geometrically nonlinear solid mechanics problems using 4-node tetrahedral elements. International Journal for Numerical Methods in Engineering, 78, 324–353.
Nguyen-Thoi, T., Vu-Do, H., Rabczuk, T., & Nguyen-Xuan, H. (2010). A node-based smoothed ﬁnite element method (NS-FEM) for upper bound solution to
visco-elastoplastic analyses of solids using triangular and tetrahedral meshes. Computer Methods in Applied Mechanics Engineering, 199, 3005–3027.
Pan, X., Sze, K. Y., & Zhang, X. (2004). An assessment of the meshless weighted least-square method. Acta Mechanica Solida Sinica, 17, 270–282.
Penrose, R. (1955). A generalized inverse for matrices. Proceedings of the Cambridge Philosophical Society, 5, 406–413.
Read, W. T. (1953). Dislocations in crystals. New York: McGraw-Hill.
Roters, F., Eisenlohr, P., Hantcherly, L., Tjahjanto, D. D., Bieler, T. R., & Raabe, D. (2010). Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity ﬁnite-element modeling: theory, experiments, applications. Acta Materialia, 58, 1152–1211.
Ullah, Z., & Augarde, C. E. (2013). Finite deformation elasto-plastic modeling using an adaptive meshless method. Computers & Structures, 118, 39–52.
Ullah, Z., Combs, W. M., & Augarde, C. E. (2013). An adaptive ﬁnite element/meshless coupled method based on local maximum entropy shape function for
linear and nonlinear problems. Computer Methods in Applied Mechanics Engineering, 267, 11–132.
Wang, Y., Zhang, C., Zhou, E., Sun, C., Hinkley, J., Gates, T. S., & Su, J. (2006). Atomistic ﬁnite elements applicable to solid polymers. Computational Materials
Science, 36, 292–302.
Zienkievicz, O. C. (1970). The ﬁnite element method: from intuition to generality. Applied Mechanics Reviews, 23, 249–256.
Zienkievicz, O. C., Taylor, R. L., & Zhu, J. Z. (2005). The ﬁnite element method: its basis and fundamentals. New York: Elsevier/Butterworth-Heinemann.
Zisman, A., & Ermakova, N. (2006). Deformation and stiffness of ﬁnite element with no assumed interpolation. Journal of Mechanical Behavior of Materials, 17,
219–234.
Zisman, A. A., Ivanov, D. S., Lomov, S. V., & Verpoest, I. (2006). Processing discrete data by gradient matrix: application to strain mapping of textile
composite. In J. Lamon & A. T. Marques (Eds.), Proceedings of the12th European Conference on Composite Materials. Biarritz: CD ROM.
Zisman, A. A., Van Boxel, S., Seefeldt, M., & Van Houtte, P. (2008). Gradient matrix method to image crystal curvature by processing of EBSD data and trial
recognition of low-angle boundaries in IF steel. Materials Science and Engineering: A, 474, 165–172.