AP04-Bittnar1.vp


1 Introduction
In standard continuum mechanics, a solid body is decom-

posed into a set of idealized, infinitesimal material volumes,
each of which can be described independently as far as the
constitutive behavior is concerned. Of course, this does not
mean that the individual material points are completely iso-
lated, but their interaction can take place only on the level of
balance equations, through the exchange of mass, momen-
tum, energy and entropy.

In reality, however, no material is an ideal continuum.
Both natural and man-made materials have a complicated
internal structure, characterized by microstructural details
whose size typically ranges over many orders of magnitude.
The expression “microstructure” will be used as a generic de-
nomination for any type of internal material structure, not
necessarily on the level of micrometers. Some of these details
can be described explicitly by spatial variation of the material
properties. But this can never be done simultaneously over
the entire range of scales. One reason is that such a model
would be prohibitively expensive for practical applications.
Another, more fundamental, reason is that, on a small enough
scale, the continuum description per se is no longer adequate
and needs to be replaced by a discrete model (or, ultimately,
by interatomic potentials based on quantum mechanics).

Constructing a material model, one must select a certain
resolution level below which the microstructural details are
not explicitly “visible” to the model and need to be taken
into account approximately and indirectly, by an appropriate
definition of “effective” material properties. Also, one should
specify the characteristic wavelength of the imposed defor-
mation fields that can be expected for the given type of
geometry and loading. Here, the notion of characteristic
wavelength has to be understood in a broad sense, not only as
the spatial period of a dynamic phenomenon but also as the
length on which the value of strain changes substantially
in static problems. If the characteristic wavelength of the
deformation field remains above the resolution level of the
material model, a conventional continuum description can be
adequate. On the other hand, if the deformation field is
expected to have important components with wavelengths

below the resolution level, the model needs to be enriched so
as to capture the real processes more adequately. Instead of
refining the explicit resolution level, it is often more effective
to use various forms of the so-called enriched or generalized
continuum formulations.

A systematic overview and detailed discussion of general-
ized continuum theories has been given e.g. in the recent
review paper by Bažant and Jirásek (2002). The aim of the
present study is to demonstrate by specific examples how
the need for enriched continuum formulations arises from
discrepancies between experimental observations and theo-
retical predictions based on the standard theories, and also
how the model performance can be improved by adding care-
fully selected enrichment terms.

The enrichments to be discussed here are in general re-
ferred to as nonlocal, but this adjective must be understood in
the broad sense, covering both strongly nonlocal and weakly
nonlocal formulations. Precise mathematical definitions of
strong and weak nonlocality were given by Rogula (1982) and
are also explained in Bažant and Jirásek (2002). Here we only
note that strongly nonlocal theories are exemplified by inte-
gral-type formulations with weighted spatial averaging or
by implicit gradient models, while weakly nonlocal theories
include for instance explicit gradient models. The meaning of
these expressions will become clear from the examples to fol-
low. We will start from enriched formulations of the theory of
elasticity, and then proceed to elastoplasticity and damage
mechanics.

The paper is organized as follows. Section 2 treats the
dispersion of short elastic waves in heterogeneous or dis-
crete media. It is shown that the standard homogenization
procedure erases the information on dispersive properties.
Dispersion laws are then derived for a host of generalized
continuum models, including strain-gradient elasticity, mod-
els with mixed spatial-temporal derivatives, and integral-type
nonlocal elasticity. Advantages and drawbacks of individual
formulations are discussed, and a general framework of non-
local strain-gradient elasticity is outlined.

Section 3 deals with size effects in microscale elasto-
plasticity, in particular with the size dependence of the appar-
ent hardening modulus. Using the academic example of a

16 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004

Nonlocal Theories in Continuum
Mechanics
M. Jirásek

The purpose of this paper is to explain why the standard continuum theory fails to properly describe certain mechanical phenomena and how
the description can be improved by enrichments that incorporate the influence of gradients or weighted spatial averages of strain or of an
internal variable. Three typical mechanical problems that require such enrichments are presented: (i) dispersion of short elastic waves in
heterogeneous or discrete media, (ii) size effects in microscale elastoplasticity, in particular with the size dependence of the apparent
hardening modulus, and (iii) localization of strain and damage in quasibrittle structures and with the resulting transitional size effect.
Problems covered in the examples encompass static and dynamic phenomena, linear and nonlinear behavior, and three constitutive
frameworks, namely elasticity, plasticity and continuum damage mechanics. This shows that enrichments of the standard continuum theory
can be useful in a wide range of mechanical problems.

Keywords: Damage mechanics, dispersion, elasticity, enriched continuum, gradient models, nonlocal models, plasticity, quasibrittle
materials, size effect, strain localization, wave propagation.



semi-infinite shear layer, it is shown that stiffer behavior of
smaller structures can be reproduced with explicit or implicit
gradient plasticity if appropriate boundary conditions are
enforced. The general trends are discussed and compared to
experimental measurements of size effect in plastic torsion of
thin wires.

Section 4 is concerned with localization of strain and
damage in quasibrittle structures and with the resulting tran-
sitional size effect. Mathematical and numerical difficulties
related to the objective description of strain localization due
to softening are explained using the one-dimensional exam-
ple of a tensile bar. It is shown that if a stress-strain law with
softening is incorporated in the standard continuum theory,
the numerical results suffer from pathological sensitivity to
the discretization parameter such as the size of finite ele-
ments. This can be remedied by special enrichments acting as
localization limiters, e.g. by a nonlocal damage formulation.
The onset of localization is studied analytically and relations
to dispersion analysis are pointed out. It is also shown that the
nonlocal model can correctly reproduce the transitional size
effect on the nominal strength of a quasibrittle structure.

All generalized theories presented here introduce a model
parameter with the dimension of length that reflects the in-
trinsic material length scale. The response of a material point
depends not only on the strain and temperature history at
that point but also on the history of a certain neighborhood
of that point or even of the entire body. For this reason, such
theories are classified as nonlocal in the broad sense.

2 Dispersion of elastic waves

2.1 Continuum versus discrete models
In the standard continuum theory, propagation of waves

in a homogeneous one-dimensional linear elastic medium is
described by the hyperbolic partial differential equation

���u Eu� �� � 0, (1)

where � is the mass density, E is the elastic modulus, u ( x, t) is
the displacement and, as usual, overdots stand for derivatives
with respect to time t and primes for derivatives with respect
to the spatial coordinate x. Since � and E are constant coeffi-
cients, equation (1) admits solutions of the form

u x t i kx t ik x ct( , ) ( ) ( )� �� �e e� , (2)

where i is the imaginary unit, � is the circular frequency, k is
the wave number, and c k� � is the wave velocity. Substi-
tuting (2) into (1), we find the condition

� � ���2 2 0Ek (3)

which implies that the magnitude of the circular frequency is
proportional to the magnitude of the wave number. The signs
of � and k are irrelevant – if both change, the real part of (2)
does not change, and if one of them changes, the wave propa-
gates in the opposite direction but otherwise remains the
same. Therefore, we will restrict ourselves to nonnegative
values and write the solution of (3) as

�
�

� k
E

. (4)

Due to the direct proportionality between � and k, the ve-
locity of a harmonic elastic wave, c k E� �� �, is constant,
independent of the wave number k. Since a wave of a general
shape can be represented by a linear combination of har-
monic waves, even such a general wave propagates at constant
velocity c and its shape remains invariant.

The situation is different in a discrete mechanical sys-
tem, which can be best exemplified by a regular infinite
one-dimensional array of mass points connected by linear
springs (Fig. 1a). The governing equations of motion form a
system of ordinary differential equations
Mu K u u K u u j Zj j j j j�� ( ) ( ) ,� � � � � �� �1 1 0 , (5)
where M is the mass of each mass point, K is the spring stiff-
ness, uj is the displacement of mass point number j initially
located at xj, and Z is the set of all integer numbers. For an
assumed harmonic wave of the form

u tj
i k x tj( )
( )

�
�

e
�

(6)

we obtain the condition

� � � � � � �� �
� �

M K K j Z
ik x x ik x xj j j j�2 1 1 01 1( ) ( ) ,

( ) ( )
e e . (7)

The spacing of mass points is assumed to be regular, i.e.,
x x x x aj j j j� � � �� �1 1 . The circular frequency correspond-
ing to wave number k is thus

� � � � � �
2

1
2

1 2
2

K
M

ika
K

M
ka

K
M

ka
( cosh ) ( cos ) sin . (8)

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 17

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 1: (a) Discrete mass-spring model with nearest-neighbor interaction and (b) the corresponding dispersion diagram with the normal-
ized frequency, � M K , plotted as a function of the normalized wave number, ka �



For this mass-spring system, the relationship between �
and k is nonlinear and the wave velocity depends on the wave
number. In such a case, one must distinguish between the
phase velocity, c kp � � , and the group velocity, c kg � d d� . In
the long wave limit (i.e., for k approaching zero), both the
phase velocity and the group velocity tend to c a K M0 � .

As is clear from the � – k diagram in Fig. 1b, c0 is the upper
bound on the velocity of waves at any wave number. As k
increases from 0 to 2� a, the corresponding phase velocity
decreases from c0 to 0. This means that shorter harmonic
waves propagate at slower velocities. For a general wave,
propagation of different harmonic components at different
velocities leads to changes of the wave shape. This phenome-
non is known as dispersion, the equation relating � and k is
called the dispersion equation, the resulting function �( )k is
the dispersion law, and its graph is the dispersion curve.

It is natural to expect that a discrete mechanical model
consisting of mass points M regularly spaced at distance a and
connected by springs of stiffness K is in some sense equivalent
to a continuum characterized by mass density � � M a and
elastic modulus E Ka� . However, this “equivalence” has its
limits. A standard homogeneous linear elastic continuum has
a linear dispersion curve with slope c E Ka M� �� 2 ,

which coincides with the initial slope c a K M0 � of the non-
linear dispersion curve of the discrete model. Both dispersion
curves are almost identical for long waves, but for wave num-
bers comparable to � a (i.e., for wavelengths comparable to
2a) they differ substantially. The standard continuum can be
considered as a long-wave approximation of the discrete
model (or vice-versa). If the actual physical system is close to
the mass-spring model, the “equivalent” continuum model
does not capture the phenomenon of dispersion of short
elastic waves. On the other hand, if a homogeneous elastic
continuum is discretized in space by finite differences (or by
finite elements with a lumped mass matrix), the resulting set

of equations has the form (5) with M a� � and k E a� , where
a is the grid parameter of the numerical method (e.g., size of
finite elements). If these ordinary differential equations are
integrated exactly in time, the solution captures correctly
long-wave phenomena but introduces an artificial (numeri-
cal) dispersion of short waves with wavelengths comparable to
the element size. Note that the numerical dispersion disap-
pears if equations (5) are integrated in time using the central
difference scheme with time step � t a c� , which is just at the
limit of numerical stability. Consequently, propagation of
shocks and other wide-spectrum phenomena are not repre-
sented accurately.

Discrete mass-spring systems are quite realistic models for
crystalline materials on a scale of observation close to the
atomic spacing. Of course, real crystal lattices are three-
-dimensional and interaction forces arise not only between
immediate neighbors but also at longer distance. Still, a one-
-dimensional lattice is an acceptable model for plane waves
propagating perpendicular to a certain set of crystallographic
planes. In this case, each mass point actually represents one
plane of densely packed atoms instead of one single atom.
Interactions at longer distance can easily be incorporated by
adding springs of stiffness K2, K3, …, KN, connecting pairs
of mass points spaced by 2a, 3a, …, Na. A straightforward
extension of the foregoing dispersion analysis yields the dis-
persion relation

� � � � � ��

�

�

�
� �M K Kn ikna
n

N

n
ikna

n

N
�2

1 1

1 1 0( ) ( )e e (9)

from which

� � � �

� �
� �

2
1

2
2

1

2

1
M

K kna
M

K
kna

n
n

N

n
n

N
( cos ) sin . (10)

18 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 2: (a) Dispersion curve of a mass-spring model with interaction up to distance , (b) dispersion curve of aluminum for longitudinal
waves in the direction (after Yarnel et al., 1965)



In the long-wave limit (k � �), the phase velocity c kp � �

tends to

c
a
M

K nn
n

N

0
2

1

�

�
� . (11)

An example of a dispersion curve constructed with N � 3
and K1 � K2 � K3 � K is shown in Fig. 2a. This curve has a more
general shape than for the nearest-neighbor interaction only
(cf. Fig. 1b), but the wave number 2� a still corresponds to
zero frequency. This is natural, because the displacements uj
generated by a harmonic wave with wavelength a are the same
at all the lattice sites. The associated mode is a uniform trans-
lation of the lattice points, which would not lead to any
vibrations.

A striking property of the dispersion law (10) is its period-
icity. This is closely related to the discrete and periodic nature
of the underlying mechanical model. In fact, wave numbers
that differ by an integer multiple of � a correspond to the
same physical state of the mass-point chain. Consequently,
the dispersion curve of this chain is uniquely characterized by
its initial part for wave numbers between 0 and � a. This
range of wave numbers is in the theory of crystal vibrations
called the first Brillouin zone. Another interesting property is
that the frequency range of harmonic waves is limited. The
dispersion law maps the first Brillouin zone onto a certain in-
terval [0, �max] where �max is a limiting circular frequency.
If certain mass point is externally excited with a circular
frequency larger than �max, the vibration does not propagate
through the entire chain and remains localized to the neigh-
borhood of the excited mass point. The discrete chain there-
fore acts as a low-band filter.

Dispersion of elastic waves in crystals is a real physical phe-
nomenon that can be observed and studied experimentally.
Fig. 2b shows an example of a dispersion curve of aluminum,
measured by Yarnel, Warren and Koenig (1965). This specific
curve corresponds to waves propagating in the crystallo-
graphic direction (110), and it can be well approximated by a
function of the form (10).

2.2. Strain-gradient elasticity
If dispersive phenomena were limited to the atomistic

scale, elastic wave propagation could be described by discrete
atomistic models on that scale and by standard continuum
models on any coarser scale. However, dispersion arises not
only due to the discrete character of the crystal lattice, but in
general due to any type of material heterogeneity. Leaving
aside the ideal case of a perfect monocrystal, the internal
structure of all materials exhibits heterogeneities on various
scales. Some defects in crystals can still be treated by atomistic
models, but in most other cases the material needs to be
considered as a continuum (because the relevant scale is
already above the atomistic one) with a certain heterogeneous
microstructure.

For elastic materials, there exist sophisticated and mathe-
matically well-founded homogenization techniques provid-
ing the effective elastic moduli of a homogeneous material
that is in a certain sense equivalent to the heterogeneous one
and can replace it in large-scale simulations. Again, this equiv-
alence is limited and holds with reasonable accuracy for
long-wave phenomena only. In the present context, waves are

considered as long if the wavelength is much larger than the
characteristic size of major heterogeneous features in the in-
ternal structure of the material. If there is a need to describe
shorter waves in a realistic manner, it is in principle possible to
explicitly resolve the details of the heterogeneous internal
structure, but this approach often leads to extreme demands
on the computational resources. Also, since the particular
microstructure is usually not known exactly but only in terms
of a stochastic description, the method of explicit resolution
would need to exploit a Monte-Carlo type of technique or use
stochastic differential equations, which again complicates the
procedure and makes it computationally expensive. As an ele-
gant and efficient alternative, it is possible to construct en-
richments of the standard continuum theory that reflect the
main features of the microstructure without using fast oscillat-
ing material properties.

In standard continuum elasticity it is assumed that the
density of elastic energy stored per unit volume, w, depends
only on the strain tensor, which is directly related to the defor-
mation gradient, i.e., to the first gradient of the displacement
field. The elastic energy stored by the entire body, W, is then
evaluated as the spatial integral of the elastic energy density.
In the one-dimensional setting, one can write

W w u x x
L

� �� ( ( )) ,d (12)
where � �u u xd d is the strain, further denoted as �, and L is
the interval representing geometrically the one-dimensional
body. In linear elasticity, the elastic energy density

w E( )� ��
1
2

2 (13)

is a quadratic function of strain.

One class of enrichments is based on the incorporation
of higher gradients of the displacement field (Toupin, 1962;
Mindlin 1964, 1965). In general, the elastic energy density
can be assumed to depend on �� ���u u u IV, , , etc. The simplest
strain-gradient theory of elasticity uses enrichment by the
second displacement gradient, ��u , which is equal to the strain
gradient, �� , further denoted as �. If we consider one single
material point only, the strain gradient is locally independent
of the strain value. In the linear case, the enriched elastic
energy density potential is written as

w E C( , )� � � �� �
1
2

1
2

2 2, (14)

where C is a higher-order elastic modulus. The variation of
elastic energy density is given by

�
	

	�
��

	

	�
�� 
�� ���w

w w
� � � � , (15)

where 
 	 	� �� �w E is the (Cauchy) stress and
� 	 	� �� �w C is the so-called double stress.

Based on the extended form of the principle of virtual
work, it is possible to derive the static equilibrium equation

( )
 �� � � � �b 0, (16)

where b is the body force density. In dynamics, b is replaced by
the inertial force density, ����u. Combining this with the con-
stitutive equations 
 �� E and � �� C and with the kinematic
equations � � �u and � � ��u , we obtain the wave equation of
strain-gradient elasticity,

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 19

Acta Polytechnica Vol. 44  No. 5 – 6/2004



���u Eu Cu IV� �� � � 0, (17)

which differs from the standard wave equation (1) by the
presence of a term with the fourth spatial derivative of
displacement.

Substituting the assumed harmonic wave solution (2) into
(17), we obtain the dispersion equation

� � � ���2 2 4 0Ek Ck (18)

from which

�
�

�
�

k
E Ck2

. (19)

The dispersion curve is plotted in Fig. 3a. The phase
velocity

c
k

E Ck
p � �

��

�

2
(20)

is not constant, except for the case when C � 0 and the model
reduces to standard elasticity. In strain-gradient elasticity it is
usually assumed that the higher-order modulus C is positive.
This assumption leads to a convex energy potential and per-
mits to generalize certain uniqueness theorems known from
standard elasticity. However, for C > 0 the phase velocity cp
increases with increasing wave number k. We have seen that
the discrete mass-spring model exhibits the opposite trend,
and this is also confirmed by measurements of dispersion
curves in real crystals. Even for heterogeneous continua, the
dispersion curves (determined experimentally or by analyti-
cal solution of some simple cases) typically have negative
curvature. So the strain-gradient theory gives a reasonable
approximation of the dispersion effect only if the higher-or-
der modulus C is negative. Convexity of the elastic potential is
then lost and uniqueness cannot be guaranteed. Indeed, if
C El� � 2 where l is a model parameter with the dimension of
length, the phase velocity of a harmonic wave with wave
number k lcrit �1 vanishes. This means that the equation of
motion (17) is satisfied by a stationary wave of wavelength
2 2� �k lcrit � . A similar result was found for the discrete
mass-spring model, but in that case the stationary wave in
reality represented a uniform translation, because the values
of the displacements had physical meaning only at discrete

points with spacing equal to the critical wavelength. In
contrast to that, a stationary wave in a continuous elastic
medium is physically inadmissible. The problem is aggra-
vated by the fact that, for wave numbers exceeding kcrit, the
circular frequency � solved from the dispersion equation (18)
becomes imaginary. This means that harmonic modes with
wavelengths shorter than 2�l would spontaneously blow up.
The source of this instability becomes apparent if one realizes
that, for short waves, the negative higher-order part the elas-
tic energy, � ��El u2 2( ) , exceeds in magnitude the positive stan-
dard part, E u( )� 2, and so the total energy density becomes
negative. If the amplitude of the wave grows, energy is re-
leased instead of being consumed.

2.3 Models with mixed spatial-temporal
derivatives

Due to the unstable behavior of short waves, equation (17)
is sometimes called the “bad Boussinesq problem”. This
equation can describe dispersion of waves with “moderate”
wave numbers but leads to instabilities if waves shorter than
the critical wavelength 2�l are involved. If the body of interest
is discretized by finite elements, the minimum wave-
length that can be captured by the numerical approximation
is proportional to the element size. Therefore, for meshes that
are sufficiently coarse with respect to the material length pa-
rameter l, the numerical solution leads to reasonable results.
However, upon mesh refinement, the solution becomes pol-
luted by unstable modes rapidly oscillating in space.

Several modifications of the bad Boussinesq problem were
proposed in the literature. Fish, Chen and Nagai (2002a)
replaced the term with the fourth spatial derivative, u IV , by a
term with a mixed derivative, ����u . Their arguments can be re-
phrased and expanded as follows: For small wave numbers,
the fourth-order term in (17) is negligible with respect to the
second-order terms, so we can write Eu u�� 	 ���. Differentiating
this twice with respect to x, we obtain Eu uIV 	 ����� . Finally,
replacing in (17) u IV by ( ) ��� E u�� and C by �El2 yields a modi-
fied wave equation

� ��� ��u Eu l u� �� � �� �2 0 (21)

20 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 3: Dispersion curves of (a) strain-gradient elasticity and gradient models with mixed derivatives, (b) nonlocal elastic models with dif-
ferent weight functions; the normalized circular frequency is �l c0 in case (a) and �a c0 in case (b), and the normalized wave
number is kl in case (a) and ka � in case (b)



which was called by Fish et al. (2002a) the “good Boussinesq
problem”. This problem can be expected to have similar
solutions to the original bad Boussinesq problem (17) at
low wave numbers but a different asymptotic behavior for
high wave numbers. Indeed, the usual procedure leads to the
dispersion equation

� � � ��� � �2 2 2 2 2 0Ek l k (22)

from which

c
k

E
k l

c

k l
p � �

�
�

�

�

�( )1 1
2 2

0
2 2

. (23)

For this model with enrichment by a mixed derivative,
the phase velocity remains real and positive for all wave
numbers. The dispersion curve, plotted in Fig. 3a, is mono-
tonically increasing, has a negative curvature, and for k � �
approaches a horizontal asymptote at circular frequency
� � c l0 . The model can reasonably reproduce dispersion of
waves at moderate wavelengths and does not give rise to
instabilities for very short waves. Its extension to multiple
dimensions is relatively complicated (Fish, Chen and Nagai,
2002b; Nagai, Fish and Watanabe, 2004).

Metrikine and Askes (2002) used a different line of reason-
ing and arrived at an equation of motion with two enrich-
ments terms, proportional to ����u and u IV . In its most general
form, this equation can be written as

� � ��� ( �� )u Eu l u Eu� �� � � �� �� �2 0 , (24)

where l is an internal length and � is an additional model pa-
rameter. For l � 0 or for � � 1, the model reduces to the stan-
dard elastic continuum, while for � � 0 it reduces to the good
Boussinesq problem of Fish et al. (2002a). Metrikine and
Askes (2002) used certain heuristic arguments to link l and �
to the microstructure and then proposed a parameter identi-
fication procedure based on a reflection-transmission test
(Askes and Metrikine, 2002). Of course, one can also consider
l and � as free model parameters and determine them by
optimal fitting of the dispersion curve for a given material.
The dispersion equation corresponding to (24),

� � � � � ��� �� �2 2 2 2 2 2 0Ek l Ek k( ) (25)

yields the phase velocity

c
k

c
k l

k l
p � �

�

�

� �
0

2 2

2 2
1
1

, (26)

where c E0 � �, as usual. If parameter � is nonnegative, the
phase velocity remains real and positive for all wave numbers.
For 0 < � < 1, the dispersion curve has a negative curvature;
see Fig. 3a. So, with a proper choice of parameters, the model
can reasonably approximate dispersion and does not suffer
by unstable behavior of short waves. Its disadvantage is that
the presence of the fourth derivative u IV requires either a
C1-continuous finite element approximation (which is hard
to construct on general meshes in multiple dimensions) or a
mixed approach with independent approximations of sev-
eral fields (e.g., of the displacement field and the strain
field). Also, nonstandard higher-order boundary conditions
are needed on the physical boundary of the investigated
body.

2.4 Integral-type nonlocal elasticity
Another class of enrichments is based on weighted spatial

averaging. The simplest model of this kind can be derived
from the elastic potential

W E x x x
LL

� ��
1
2

( , ) ( ) ( )
 � � 
 
d d , (27)

where E x( , )
 is a function describing the generalized elastic
modulus. The variation of elastic energy is evaluated as

� 
 �� � 
 


 � �� 

W E x x x

E x x x

LL

L

� �

�

��
1
2
1
2

( , ) ( ) ( )

( , ) ( ) ( )

d d

d��

��

�

� �

d

d d




 
 � 
 
 ��

L

LL
E x E x x x

1
2

[ ( , ) ( , )] ( ) ( ) .

(28)

This can be written in the usual form � 
 ��W x x x
L

� � ( ) ( ) d
if the stress is defined as


 
 � 
 
( ) ( , ) ( )x E xs
L

� � d , (29)

where

E x E x E xs( , ) [ ( , ) ( , )]
 
 
� �
1
2

(30)

is the elastic modulus function symmetrized with respect
to its arguments. The corresponding equilibrium equations
derived from the principle of virtual work then keep their
standard form, � � �
 b 0. Consequently, the wave equation
for this model reads

�
	

	

	

	


	 

	


2

2
0

u x t
t x

E x
u t

s
L

( , )
( , )

( , )
� �� d . (31)

Since function E xs( , )
 reflects the strength of long-dis-
tance interaction between points x and 
, its value can be
expected to be negligible if the distance between x and 
 is
large compared to the internal length of the material (which
corresponds to the characteristic size and spacing of major
heterogeneities). For functions Es with a sufficiently fast decay,
the integrals in (29) and (31) make sense even if the integra-
tion domain L is considered as the entire real axis. If the
body infinite and macroscopically homogeneous, function
E xs( , )
 should depend only on the distance between x and 
.
Bearing in mind these restrictive assumptions, we present the
modulus function in the form

E x E xs( , ) ( )
 � 
� �0 0 , (32)

where E0 is a reference value of the elastic modulus and �0 is
a dimensionless even function, further called the nonlocal
weight function. The second term in the wave equation (31)
can then be transformed as follows:

	

	


	 

	


	� 

	

	 

	


x
E x

u t
E

x
x

u t
s( , )

( , ) ( ) ( , )
d d

�





�
� �

�
0

0







�





�

�

�

� �
�

� �E
x
x

u t
E x

u t
0

0
0 0

2	� 

	

	 

	

 � 

	 
( ) ( , )
( )

( ,
d

)
.

	


2
d

�





�

(33)

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 21

Acta Polytechnica Vol. 44  No. 5 – 6/2004



Substituting the assumed harmonic form of an elastic
wave (2) into the transformed wave equation

�
	

	
� 

	 

	


2

2 0 0

2

2
0

u x t
t

E x
u t( , )

( )
( , )

� � �

�





� d (34)

we obtain the dispersion equation

� � ��� �2 0
2

0 0E k k
* ( ) (35)

in which

� �0 0
* ( ) ( )k r rikr� �

�





� e d (36)

is the Fourier image of the nonlocal weight function �0( )r .
Finally, the phase velocity is evaluated as

c
k

E k
c kp � � �

� �

�
�0 0 0 0

*
*( ) ( ) . (37)

Relation (35) shows that there is a unique correspondence
between the dispersion law and the Fourier image of the
nonlocal weight function. If the dispersion law of a certain
material is given in the form �( ) ( )k kc kp� where c kp ( ) is a
known function, it is possible to construct a nonlocal elasticity
model that exactly reproduces the dispersion properties. For
this, it suffices to set E c0 0

2� � where c cp0 0� ( ), and to evaluate

the weight function by the inverse Fourier transform, taking
into account that the phase velocity does not depend on the
sign of the wave number:

�
�

�
�

0 0
0
2

21
2

1
( ) lim ( ) ( ) cos*r k k

c
c k kr k

k

ikr

k

k

p� �
�


�
� e d d

0




� . (38)

Of course, this is possible only if the dispersion law to be
reproduced has reasonable properties, such that the inverse
Fourier transform exists. For instance, for the dispersion law
corresponding to the bad Boussinesq problem, c kp ( ) is given
by (20) and the integral in (38) does not converge (independ-
ently of the sign of C). On the other side, for the good
Boussinesq problem (21) we have c k c k lp ( ) � �0

2 21 and the

inverse Fourier transform (38) yields

�
�

0 2 2
0

1
1

1
2

( )
cos

r
kr

k l
k

l

r
l�

�
�


 �

� d e . (39)

A nonlocal elasticity model with this particular weight
function gives exactly the same dispersion law as the model
with enrichment by a mixed derivative proposed by Fish et al.
(2002a). On an infinite domain, both models are equivalent.

The transformation of the model of Fish et al. (2002a)
into a integral-type nonlocal model can also be performed
directly. At a fixed time instant, equation (21) can be written as

�� ��u l u
E

u� �� � ��2
�

(40)

and interpreted as an ordinary differential equation for the
unknown acceleration ��u, with the current displacement u
considered as known. Equation (40) has the form of the
so-called Helmholtz equation, and its solution satisfying con-
ditions of boundedness (which play the role of boundary
conditions at plus and minus infinity) can be expressed as

��( , ) ( , ) ( , )u x t
E

G x u t� ��

�





�� 
 
 
d , (41)

where G x( , )
 is the Green function of the Helmholtz equa-
tion, formally obtained as the solution of this equation with
the Dirac distribution � 
( ) on the right-hand side. It turns out
that the Green function is in this case given by

G x
l

x
l( , )



�

� �
1
2

e (42)

and so equation (41) is in fact equivalent with (34) if the
nonlocal weight function �0 is selected according to formula
(39).

2.5 Combination of nonlocal averaging and
strain gradients

For the dispersion law corresponding to the gradient
model proposed by Metrikine and Askes (2002), the integral
in the inverse Fourier transform (38) does not converge. So
this model is not equivalent to any integral-type nonlocal elas-
tic model derived from the potential (27). Still, using the al-
ternative procedure based on the Green function, it is possible
to construct a more general nonlocal model equivalent to the
original enriched gradient formulation. Indeed, rewriting
(24) as

�� �� ( )u l u
E

u l u IV� �� � �� �2 2
�

� (43)

and “solving” for the acceleration, we obtain

��( , ) ( , ) ( , ) ( , ) ( , )u x t
E

G x u t l G x u tIV� �� �

�





�� 
 
 
 � 
 
d d
2 

�





�
�



�
�
�

�

�

�
�
�

. (44)

The result resembles the wave equation of strain-gradient
elasticity (17), but with the derivatives ��u and u IV replaced by
their weighted spatial averages. This observation motivates
the development of a nonlocal strain-gradient model with the
elastic energy potential given by

W E x x x

C x x x

LL

L

� �

�

��

�

1
2
1
2

( , ) ( ) ( )

( , ) ( ) ( )


 � � 
 


 � � 
 

d d

d d
L� .

(45)

Taking the variation and using the same line of reasoning
as for the basic version of nonlocal elasticity, we identify the
constitutive equations


 
 � 
 
( ) ( , ) ( )x E xs
L

� � d (46)

� 
 � 
 
( ) ( , ) ( )x C xs
L

� � d (47)
in which C x C x C xs( , ) [ ( , ) ( , )]
 
 
� � 2 is the symmetrized
higher-order modulus. Substituting this into the equation of
motion of the strain-gradient theory, � 
 ��� ( )u � � � � � 0, and
using the kinematic relations � � �u and � � ��u , we obtain the
wave equation

�
	

	

	

	


	 

	


	

	


2

2

2

2

u x t
t x

E x
u t

x
C x

s
L

s

( , )
( , )

( , )

( , )

� �

�

� d
	 

	


2

2
0

u t
L

( , )
d� �

(48)

22 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004



that generalizes equation (31). For the moduli functions in
the special form

E x E xs( , ) ( )
 � 
� �0 0 (49)

C x E l xs( , ) ( )
 � 
� �0
2

1 (50)

dispersion analysis provides the following expression for the
phase velocity:

c k c k k l kp ( ) ( ) ( )
* *� �0 0

2 2
1� � . (51)

Here, �1
* denotes the Fourier image of function �1.

The nonlocal strain-gradient model just presented is quite
general and covers as special cases all the other models
discussed so far. The special choices of weight functions �0
and �1 are summarized in Table 1. The model permits to re-
produce a very wide class of dispersion laws. Which terms
need to be activated depends on the asymptotic behavior of
the dispersion curve at wave numbers approaching infinity. If
the dispersion law �(k) is bounded, it is sufficient to use a reg-
ular weight function �0 obtained by the inverse Fourier trans-
form of ( )� c k0

2. If �(k) tends to infinity but remains of order
O(k), it is possible to use either a weight function �0 with a
singular part of the Dirac type, or regular functions �0 and �1.
Finally, if �(k) grows superlinearly but remains of order O(k2),
function �1 must have a singular Dirac-type part. Still faster
growth could be reproduced by models with second (Mindlin,
1965) or still higher (Green and Rivlin, 1964) strain gradi-
ents, but it seems that dispersion laws of real materials do not
exhibit such behavior, so this question is purely academic. In
fact, all dispersion laws with superlinear growth at k � � are
suspicious because the phase velocity becomes unbounded
and disturbances can propagate at an arbitrary velocity if the
wavelength is selected as sufficiently short.

2.6 Nonlocal model reproducing dispersion of
discrete lattice

It is interesting that the nonlocal elastic model with a
weight function �0 that linearly decreases from its maximum
value at r � 0 to zero at r � a and vanishes for r > a leads to
exactly the same dispersion law as the simplest mass-spring
model with nearest-neighbor interactions and with spacing
a between neighboring mass points. Table 2 shows several
other nonlocal weight functions and their Fourier images that
can be used to construct the corresponding dispersion curves,

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 23

Acta Polytechnica Vol. 44  No. 5 – 6/2004

model �0( )r �1( )r

standard elasticity �( )r 0

strain-gradient elasticity �( )r �( )r

Fish et al. 1
2l

r
le

�

0

Metrikine and Askes 1
2l

r
le

�
�

2l

r
le

�

mass-spring chain
a r

a

�
2 0

The brackets � denote the positive part operator, defined
by x x x� �( ) 2.

Table 1: Special cases of nonlocal strain-gradient elasticity

model �0( )r �0
* ( )k

standard elasticity �(r) 1

mass-spring (neighbors)
a r

a

�
2

2 1
2 2

( cos )� ka
k a

mass-spring (long interaction)
K na r

a n K

n

n

��
�2 2

2 1
2 2 2

K kna

k a n K

n

n

( cos )��
�

Eringen see eq. (54)
2 1

2 2
( cos )� ka

k a
if k

a
�

�

uniform averaging
1

2a
for r a�

sin ka
ka

Fish et al.
1

2a
r ae�

1
1 2 2� k a

Gauss weight function
1 2 2

a
r a

�
e� e 4�k a

2 2

quartic weight function
15

16
1

2

2

2

a
r
a

�
15

3 35 5
2 2

k a
k a ka ka ka[( ) sin cos ]� �

Table 2: Nonlocal weight functions and their Fourier images



which are plotted in Fig. 3b. The dispersion law of a mass-
-spring model with long-distance interactions can be repro-
duced by a nonlocal model with a piecewise linear weight
function whose characteristics uniquely depend on the spring
stiffnesses. If all stiffnesses used by the discrete model are
positive, the weight function is concave (for nonzero r). The
dispersion law gives real frequencies for all wave numbers, but
for wave numbers that are integer multiples of 2� a the fre-
quency vanishes, i.e., the model admits stationary waves of
wavelength a, a 2, a 3, etc. This is natural for the discrete
model, as already explained, but the same property is shared
by the nonlocal continuum model.

For the simplest weight function, constant for r between
0 and a and vanishing for r > a, the dispersion law gives
real frequencies only for wave numbers in intervals [0, �/a],
[2�/a, 3�/a], [4�/a, 5�/a], etc. Between these bands, the fre-
quency becomes imaginary, which indicates an instability. The
potential appearance of periodic modes that carry no energy
for the nonlocal model with uniform strain averaging over a
finite neighborhood was mentioned by Bažant and Chang
(1984). It is clear that every periodic function with period 2a
is mapped by the nonlocal operator onto a zero function,
and therefore has no influence on the nonlocal average.
Each such function can be decomposed into a sum of har-
monic functions of wavelengths 2a, 2 2a , 2 3a , etc., which
correspond to zero frequencies. Thus the static solution can
be modified by a periodic function with period 2a with-
out disturbing equilibrium. The dispersion analysis indicates
that in dynamics the situation is even worse, because har-
monic modes with wave numbers in the intervals (�/a, 2�/a),
(3�/a, 4�/a), etc., are associated with imaginary frequencies
and would blow up.

To avoid the potential appearance of unstable modes, it is
sufficient to use weight functions with positive values of their
Fourier images for any positive wave number k. This is the
case for instance for the Green function of the Helmholtz
equation, and also for the Gauss-type weight function

�
�

0
1 2 2

( )r
a

r a� �e (52)

that is often used by nonlocal models. On the other hand, the
Fourier image of the truncated quartic polynomial function

�0

2

2

2
15

16
1( )r

a
r
a

� � (53)

is positive only for wave numbers smaller than 5.763/a. Con-
sequently, instabilities could develop for fine meshes with ele-
ment size in the order of a.

The relationship between atomic lattices and the non-
local integral models was studied by Eringen (1972), who
proposed a weight function that should correspond to the
mass-spring model. However, Eringen did not use the com-
plete inverse Fourier transform but integrated (38) only in
the limits � � �� �a k a. Consequently, his nonlocal model
would reproduce the dispersion law of the mass-spring model
only for wavelengths larger than 2a. All smaller wavelengths
would be associated with a zero frequency, i.e., stationary
waves could easily appear. Eringen’s weight function

�
� � �

� � � �

�

0

2

3 2

4 2

2 1
1
2

1
2

2

( )

cos

x
a

x
a

a

�
� �

� �
�

�
�

�

�
�

�

� �

si

si

2

2 2 4 2

1
4 2 2

x
x
a

x
a

x
a

x
a

x
a

�


� � �

� �
�

�
�

�

�
� �

�

�
�

�

�
� �

sin si

si� �
1
4 2 2

� ��
�

�
�

�

�
� �

�

�
�

�

�
�
�

��
x
a

x
a

si

(54)

(where si dx
x

� �� � � �
1

0

sin ) is also much more complicated

than the piecewise linear function reproducing the dispersion
of the discrete model exactly for arbitrary wave numbers, and
does not seem to be very practical.

3 Size effects in microscale plasticity

3.1 Experimental observations
Among phenomena that are hard to model and predict by

standard continuum theories, one finds also various forms
of size effects on the apparent “material” properties. Such
effects can be observed already in the elastic range. For
instance, according to standard elasticity, the torsional stiff-
ness of a prismatic beam with a circular cross section should be
proportional to the shear modulus of elasticity, G, and to the
third power of the sectional diameter, D. However, if the tor-
sional stiffness is evaluated experimentally as the ratio be-
tween the torque and the relative twist angle, it turns out that
the expected proportionality to D3 holds with sufficient ac-
curacy only for diameters larger than a certain threshold
value (Morrison, 1939; Lakes, 1986; Fleck, Muller, Ashby and
Hutchinson, 1994). If the results obtained for thick wires are
extrapolated to thin wires, the actual stiffness is underesti-
mated. In the context of the standard theory, this can be inter-
preted as a dependence of the elastic modulus on the size of
the sample. However, such an explanation is not satisfactory
from the theoretical point of view.

In a more fundamental approach, it is admitted that the
standard continuum elasticity theory provides only a large-
-size approximation to the static torsion problem, just as it
provides only a long-wave approximation to the dynamic
wave dispersion problem. If the size of the structure is com-
parable to a certain internal length scale of the material,
higher-order effects appear and the classical concept of a
homogeneous local continuum needs an adjustment. In elas-
ticity, such an adjustment that accounts for the principal
features of the microstructure is provided by various theories
enriched by higher-order gradients or integral-type nonlocal
terms, already exposed in Section 2.

In general, by size effect we mean a situation when a cer-
tain parameter normally considered as a material property
appears to be dependent on the size of the sample or speci-
men for which it is evaluated. The reason for this unexpected
behavior is usually that the specimens are either too small or
too big and the underlying theory is not adequate on the
extreme scale. Important size effects in microscale plasticity
have been detected in indentation tests (Nix, 1989; Ma and
Clarke, 1995; Poole, Ashby and Fleck, 1996), bending of thin
sheets (Stolken and Evans, 1998), plastic torsion (Fleck et al.,

24 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004



1996), and void growth. From the practical point of view,
proper understanding and modeling of size effects on the
small scale is essential for applications of simulation methods
in the analysis and design of micro- and nano-devices used in
modern technology.

3.2 Illustrative example: shear layer
To illustrate the power of plasticity theories extended by

gradient or integral-type nonlocal terms, we will analyze a
rather academic but instructive problem of a material layer of
thickness L, placed between stiff loading plates and loaded
statically by shear. This problem can be modeled in one spa-
tial dimension, which facilitates its analytical solution. The
choice of the coordinate system and the type of loading are
shown in Fig. 4. It is assumed that the elastic properties of the
material are isotropic and plastic flow is isochoric (preserves
volume), so that volumetric and deviatoric effects can be de-
coupled. If this is the case, the relevant components of stress
and strain are the shear stress 
 �xy � and the engineering
shear strain 2� �xy � . Both of them are considered as inde-
pendent of spatial coordinates y and z. From the equilibrium
equation with vanishing body forces, it follows that � is also
independent of x, i.e., it is uniformly distributed in space
and varies only as a function of the pseudo-time parameter-
izing the loading process.

Classical plasticity with linear isotropic hardening is de-
scribed by the basic equations

� � �� �G p( ) (55)
� � sgn( )� � �p � (56)
� , ( , ) , � ( , )� � � � � �� � �0 0 0f f (57)

where �p is the plastic strain and � is the cumulative plastic
strain. The yield function f is given by

f Y( , ) ( )� � � � �� � (58)
where �Y is the current yield stress in shear, evaluated from
the hardening law

� � � �Y H( ) � �0 (59)

with �0 � initial yield stress in shear and H � hardening
modulus (in shear). During monotonic loading with positive
value of shear stress �, there is no difference between the
plastic shear strain �p and the cumulative plastic strain �.
We will therefore replace �p by � and call it simply the plas-
tic strain.

If the hardening modulus is positive and the loading is
monotonic, the stress � uniquely determines the correspond-
ing strain �. Since the stress is uniform, the strain must be uni-
form as well. Therefore, the strain at any point is equal to the
average strain determined as the ratio between the relative
displacement of the loading plates, �v v L v L� � �( ) ( )2 2 ,
and the layer thickness, L. The stress-strain curve can be di-
rectly determined from the measured dependence between v
and the tangential traction on the boundary, t, and it should
be independent of the layer thickness. Therefore, the stan-
dard plasticity model does not indicate any size effect.

3.3 Explicit gradient plasticity
If the layer thickness is comparable to the characteristic

length of the material microstructure, the assumption of the
local dependence of stress on the history of strain at the same
material point becomes questionable. The reason is that the
hardening process does not take place at each infinitely small
material point separately and independently of the surround-
ing points. This can be taken into account by sophisticated
models that consider the details of the hardening mecha-
nisms, e.g., by discrete dislocation models. As a simpler alter-
native, one can use enriched continuum models that take into
account the micromechanical processes only ”on the aver-
age”, but by terms of a higher order than in the standard con-
tinuum theory. Motivated by certain micromechanical consid-
erations, Aifantis (1984) proposed a family of models with the
yield stress dependent not only on the value of the cumulative
plastic strain (internal variable driving the hardening process)
but also on its first and second gradients. In the one-dimen-
sional setting and for linear hardening, the simplest version
of the Aifantis gradient plasticity model replaces the harden-
ing law (59) by

� � � � � �Y , H Hl( )�� � � � ��0
2 (60)

where l is a parameter with the dimension of length.
The elastic part of the model remains unchanged, and so

the strain is uniform in the elastic range. After the onset of
yielding, the yield condition must be satisfied in the plastic
zone. It is easy to show that in the present case the plastic zone
must extend over the entire layer. The yield condition f � 0
combined with the hardening law (60) then provides the ordi-
nary differential equation

� �
� �

� �� �
�

l
H

2 0 (61)

for the unknown function �. This equation should be supple-
mented by boundary conditions at the layer boundaries. The
choice of the specific form of boundary conditions has a
major influence on the solution and on the resulting size
effect. For homogeneous Neumann boundary conditions,
enforcing a vanishing normal derivative of � at the boundary,
the uniform solution obtained with the classical model would
remain valid even for the gradient model. However, if the
shear layer is fixed to rigid or very stiff loading plates and the
bond between the two materials is perfect, the boundary (or
rather the material interface) acts as an obstacle for disloca-
tion motion, which is the main mechanism of plastic flow in
crystalline materials. This has been confirmed by simulations
based on discrete dislocation models (Shu, Fleck, van der
Giessen and Needleman, 2001) . In the extreme case, plastic
flow is completely prevented and the cumulative plastic strain

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 25

Acta Polytechnica Vol. 44  No. 5 – 6/2004

Fig. 4: Material layer between stiff plates loaded by shear



� must vanish at the boundary. For homogeneous Dirichlet
boundary conditions, �( )� �L 2 0, the solution of (61) is

�
� �

( )
cosh

cosh
x

H
x l

L l
�

�
�

�

�
��

�

�
��

0 1
2

. (62)

The distribution of plastic strain (normalized by the fac-
tor ( )� �� 0 H that would correspond to the value of uniform
plastic strain in the standard theory) across the layer thickness
is plotted in Fig. 5b for different relative sizes, L l. If the
layer thickness is substantially larger than the material length,
the plastic strain is almost uniform, except for narrow bound-
ary layers. With decreasing structural size L, the relative im-
portance of these boundary layers with reduced plastic strain
levels increases, which makes the overall response stiffer.

To characterize the overall response of the layer, the rela-
tive tangential displacement of one loading plate with respect
to the other can be evaluated as

�v x x
G

x x

L
G H

L

L

L

L

� � �
�


�
�

��
�

� �
�

�

�

�

�

� ��
�

�

� � �

( ) ( )d d

2

2

2

2

0 L l
L
l

�
�

�
�

�

�
�2

2
tanh .

(63)

Defining the average shear strain ~� � �v L, we can trans-
form (63) into the average stress-strain law

� � � �� � �0 0
~

(~ )Gep , (64)

where � �0 0� G is the limit elastic shear strain and

~
tanhG

G H
l

L
L
l

ep � � �
�

�
�

�

�
�

�


�
�

��

�
1 1

1
2

2

1
(65)

is the elastoplastic (tangent) shear modulus.
According to (64)–(65), the overall response of the speci-

men after the onset of yielding is linear and is equivalent to
the response of the same specimen made of the standard
elastoplastic material with hardening modulus

~

tanh
H

H
l

L
L
l

�
�1

2
2

. (66)

This “apparent” hardening modulus depends not only
on the material parameters H and l but also on the layer

thickness L, i.e., on the size of the specimen for which it is
evaluated; see Fig. 5a. Therefore, it cannot be considered as
an intrinsic material parameter.

If the actual behavior of the material is close to the gradi-
ent model but the experimental results are interpreted using
the standard elastoplastic model, the value of the hardening
modulus will appear to be size-dependent. The reason is that
the model is oversimplified and the comparison between
theory and experiments uses only global characteristics such
as the measured relation between the loading force and the
relative displacement of one loading plate with respect to
the other. If detailed measurements of the strain field inside
the specimen were available, they would reveal a discrepancy
between the actual strain distribution and the theoretical
solution based on the oversimplified assumptions. This would
give a hint regarding the necessary refinement of the theoreti-
cal model by inclusion of higher-order terms. But even if
such detailed measurements are not possible or not available,
the development of an appropriate enriched continuum the-
ory can be guided by the experimentally detected size effect.
It turns out that the size effect for one specific type of test
performed on a series of geometrically similar specimens
can often be well reproduced by several different types of
enriched models that are not necessarily equivalent in a gen-
eral case.

Ideally, the enriched theory should be verified by several
tests leading to different types of stress and strain distribu-
tions, and also supported by micromechanical arguments and
confirmed by observations of the actual processes in the
microstructure. Only then can the model be assumed to rea-
sonably reproduce the actual material behavior and to have
some predictive power. If only one series of experiments is
fitted, the model can usually serve for reliable interpolation in
the limits that have been covered by the experiments, but
extrapolation to smaller or larger sizes can be dangerous.

The diagram linking the shear stress to the average shear
strain for different thicknesses L of the shear layer is plotted
in Fig. 6a, and the dependence of the apparent hardening
modulus

~
H on the specimen size (layer thickness) is shown in

Fig. 5a. If the layer thickness is much larger than the material
length l, the stress-strain curve is practically size-independent
and the apparent hardening modulus

~
H evaluated from the

26 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 5: Explicit gradient plasticity: (a) dependence of the apparent hardening modulus on the size of the specimen (layer thickness),
(b) distribution of plastic strain across the layer



test is very close to the model parameter H, considered here
as an intrinsic material property. So the standard theory can
be used as a good approximation if the specimen is large. For
layer thicknesses smaller than about 20l, the stress-strain
curve becomes size-dependent and it rises above the basic
stress-strain curve valid in the large-size limit. The appar-
ent hardening modulus increases with decreasing size. Even
though the shear test of a layer is difficult to perform, this
trend can be considered as realistic because stronger harden-
ing for smaller specimens is indeed observed in experiments
with plastic torsion of tiny wires or microbending of thin films.
This is documented in Fig. 6b, adapted from Fleck et al.
(1994). The figure shows the dependence of the normalized
torque, M D3, on the normalized angle of twist, �D, which is
equal to the shear strain on the wire surface. According to any
theory with stress at a point dependent on the history of strain
at that point only, the resulting curve should be independent
of the wire diameter D. In reality, this is true only for suffi-
ciently thick wires. The figure clearly shows that for diameters
in the order of dozens of micrometers, the actual response
after the onset of yielding is stronger than expected. Aifantis
(2003) and Fleck and Hutchinson (2001) have demonstrated
that this size effect can be captured by gradient plasticity
theories similar to the one presented here.

3.4. Implicit gradient plasticity
The model with yield stress dependent on the second

gradient of cumulative plastic strain falls into the category of
explicit gradient models. Recently it turned out that certain
advantages can be gained by using implicit gradient form-
ulations, which are closely related to integral-type nonlocal
models. A prominent example is the so-called ductile damage
model of Geers, Engelen and Ubachs (2001). This model
introduces the dependence of the yield stress on the nonlocal
cumulative plastic strain, �, defined implicitly as the solution
of a Helmholtz-type differential equation with the local cumu-
lative plastic strain on the right-hand side. For a one-dimen-
sional problem, this equation reads

� � �� �� �l2 , (67)

where l is, as usual, a model parameter with the dimension of
length. In the simplest case, one could replace � in (58)–(59)

by �, but such a formulation would have certain deficiencies.
Micromechanical arguments based on the idea of a plastically
hardening matrix weakened by growing voids lead to the fol-
lowing expression for the yield stress:

� � � � � � �Y pH( , ) ( )[ ( )]� � �0 1 . (68)

Here, the term � �0 � H represents the yield stress of the
matrix, which exhibits hardening driven by the local value of
the cumulative plastic strain, while 1 � � �p ( ) is an integrity
factor taking into account the reduction of the effective area
by voids that carry no stress. Void propagation is assumed to
be driven by the nonlocal cumulative plastic strain, �, and is
reflected by the ductile damage function �p that vanishes for
� � 0 and increases later on.

Due to the nonlinear format of the hardening law (68),
analysis of the plastic strain distribution in a general state
would lead to a nonlinear differential equation. To allow
for an analytical treatment, we restrict attention to the initial
distribution of the plastic strain rate. Differentiating (68)
with respect to time and using the consistency condition
�

� �f Y� � �� � 0, we obtain the differential equation

H Hp
p

( ) � ( ) � �1 0� � � �� � � �
�

�
� �

d

d
. (69)

At the onset of yielding, �, � and � p vanish in the entire
layer and the derivative d d� �p evaluated at � � 0 has the
same value �p

* at all points x. Taking all this into account

and substituting the rate form of (67) into (69), we get the
linear differential equation with constant coefficients,

( ) � � �*H Hlp� � �� �� � � � �0
2 . (70)

In the absence of body forces, static equilibrium implies
that the shear stress � must be uniform, and so the right-hand
side in (70) is constant. If H p� �� �0 0

* , the general solution

reads

�

�

cosh sinh�� �
�

�

� �
x

H
C

x
l

C
x
l

� � �
2 1 2

, (71)

where C1 and C2 are integration constants and

a Hp� �1 0� �
*

is a dimensionless parameter introduced for convenience.

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 27

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 6: (a) Stress versus average strain for different thicknesses of the shear layer, (b) normalized torque, M D3, versus surface shear
strain for twisted copper wires of various diameters D (after Fleck et al., 1994)



The particular solution depends on the boundary condi-
tions, which should be formulated in terms of the nonlocal
cumulative plastic strain and its normal derivative. For homo-
geneous Neumann boundary conditions, � � �� ( )L 2 0, the
integration constants vanish and the solution is uniform.
Consequently, no size effect is predicted by the model. This is
similar to what happens for the Aifantis model with the ho-
mogeneous Neumann boundary conditions � � �� ( )L 2 0. A
stiff boundary with a perfect bond inhibits propagation of
voids, which can be expressed by the condition that � vanishes
at the boundary. With this homogeneous Dirichlet condition
applied at both parts of the boundary, the particular solution
of (70) becomes

� ( )
� cosh( )

cosh( )
�

�

�

�

�
x

H

x l
L l

� �
�


�

�

�
�2 1 2

(72)

and the corresponding local cumulative plastic strain rate is

� ( )
�

( )
cosh( )

cosh( )
�

�

�
�

�

�
x

H

x l
L l

� � �
�


�

�

�
�2

21 1
2

. (73)

Following the same procedure as in Section 3, we find that
the shear stress rate is linked to the average shear strain rate
by the linear relation

�

~ ~�� �� Gep , (74)

where

~
~G G H

ep � �
�

�
�

�

�
�

�
1 1

1
(75)

is the elastoplastic shear modulus and

~

( ) tanh
H

H
l
L

L
l

�
� �

�

�
�

�

2

21 1
2

2

(76)

is the size-dependent hardening modulus. For very large sizes
L, the hardening modulus tends to H�2, which is the value
corresponding to the local theory with �p computed from �
instead of �. For values of L comparable to l, the hardening
modulus increases with decreasing size, which means that the
average response becomes stiffer. So the present model pre-
dicts a qualitatively similar trend as the explicit version of gra-
dient plasticity.

The ratio between the size-dependent hardening modu-
lus and its large-size limit is plotted in Fig. 7a as a function of
the relative size, L l, for parameter � � 0 5. . Fig. 7b shows the
distribution of plastic strain rate across the layer for different
ratios L l. Even though the general trends are the same as
for the explicit gradient model, certain differences can be
revealed by comparing Fig. 7 with Fig. 5. For the explicit
model, the hardening modulus tends to infinity as the layer
thickness approaches zero, while for the implicit model it
tends to a finite value. The local plastic strain on the bound-
ary is zero for the explicit model (as dictated by the boundary
condition), while for the implicit model it is positive (because
the boundary condition is formulated in terms of nonlocal
plastic strain). For the explicit model, the profile of plastic
strain distribution across the layer thickness keeps the same
shape and grows proportionally during hardening, while for
the implicit model the analytical solution covers only the ini-
tial distribution of the plastic strain rate and later the shape of
the profile would change.

4 Strain localization due to softening

4.1 Problems with objective description of
softening

In many structures subjected to high levels of solicitation,
the initially smooth distribution of strain at a certain critical
stage abruptly changes into a localized pattern with large
strains concentrated in relatively narrow regions. This phe-
nomenon, called strain localization, can be caused for in-
stance by the softening character of the material response.
The general definition of softening is more involved, but in
the one-dimensional case softening means decreasing stress
at increasing strain. The physical source of softening usually
resides in the growth and coalescence of defects such as
voids or cracks. From the micromechanical point of view, this
means that the internal structure of the material is evolv-
ing and the approximate description of the material as a
macroscopically homogeneous one may become question-
able. Indeed, softening incorporated into standard inelastic
continuum models leads to serious mathematical and numer-
ical problems, and enriched theories are needed to provide
an objective description of the softening process.

28 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 7: Implicit gradient plasticity: (a) dependence of the apparent hardening modulus on the size of the specimen (layer thickness),
(b) initial distribution of plastic strain rate across the layer



The essence of the localization problem will be explained
using a one-dimensional example, which could be inter-
preted as localization of shear strain in a layer of elastoplastic
material between two rigid plates (already studied in the pre-
vious section in the context of size effects). However, to better
show various facets of nonlocal theories and their broad appli-
cation field, we will discuss the closely related problem of a
prismatic bar of length L and cross-sectional area A under
uniaxial tension. The bar is made of a softening material de-
scribed by a continuum damage model rather than by an
elastoplastic model. Damage mechanics is frequently used for
quasibrittle materials such as concrete under predominantly
tensile loading. Many results and conclusions of the present
section could be directly reinterpreted in terms of shear bands
in ductile materials such as metals or confined soils.

In the one-dimensional setting, the stress-strain law used
by the damage model reads


 �� �( )1 � E , (77)

where 
 and � are the (normal) stress and strain, E is Young’s
modulus of elasticity, and � is a scalar damage variable
characterizing the current size and density of defects that
reduce the effective area capable of transmitting stress. In
the “virgin”, undamaged state of material with no defects (or
with very small initial defects that are incorporated in the
elastic modulus), the value of the damage parameter is zero,
and it remains zero throughout the elastic stage of loading.
When the elastic limit is exceeded, damage starts growing
and the elastically computed stress E� is reduced by the integ-
rity factor 1 � �. The limit value � �1 corresponds to a fully
damaged material that can no longer carry any stress. The
growth of damage must be described by an appropriate dam-
age evolution law for the internal variable �. This law could
be postulated in the rate form, but a particularly simple and
practical formulation is obtained with a damage law in the to-
tal form

� � g( )� , (78)

where � corresponds to the maximum level of strain reached
in the previous history of the material. Mathematically, the
internal variable � is defined by the loading-unloading condi-
tions

� , , ( ) �� � � � � �� � � � �0 0 0. (79)

During monotonic loading, � is equal to the strain �, and
so the damage evolution function g that appears in (78) can
easily be identified from the monotonic stress-strain curve.

Now suppose that the stress-strain diagram is linear up to
a certain strain level �0, after which stress decreases as a linear
function of strain until the zero stress level is reached at strain
�f (Fig. 8a). This linear softening model can be considered as
the simplest description of concrete cracking under tension.
Due to the heterogeneous and quasibrittle nature of the mate-
rial, a contiguous stress-free crack across the entire section of a
bar does not form instantaneously but is obtained as the final
result of propagation and coalescence of many smaller cracks.
Consequently, even after the onset of cracking the bar can still
transmit a certain force but its residual strength decreases as
the cracks evolve.

Under static loading and in the absence of body forces, the
stress in the bar must be uniformly distributed. In the elastic
range, strain is a unique function of stress, and so the strain
distribution must be uniform as well. When the peak stress
(tensile strength) ft is attained, the uniqueness of the response
is lost. Stress must still remain uniform and it can only de-
crease, but a given stress level can be reached either by further
stretching the material into the softening regime, or by elastic
unloading. Consequently, there are many different spatial dis-
tributions of the strain increments that lead to the same uni-
form stress decrement and thus represent a valid solution sat-
isfying the equilibrium equation and the constitutive law. The
compatibility equations do not represent any important con-
straint in the one-dimensional case, because the strain field
can always be integrated to yield the displacement field. For
example, the material can be softening in an interval of
length Ls and unloading everywhere else. When the stress is
completely relaxed to zero, the strain in the softening region
is �f and the strain in the unloading region vanishes; the total
elongation of the bar is therefore u L s f� � . The length Ls re-
mains undetermined, and it can have any value between 0
and L. This means that the problem has infinitely many solu-
tions; the corresponding post-peak branches of the load-dis-
placement diagram fill the fan shown in Fig. 8b.

The ambiguity is removed if imperfections are taken into
account. Real material properties and sectional dimensions
cannot be perfectly uniform. Suppose that the strength in a
small region is slightly lower than in the remaining portion of
the bar. When the applied stress reaches the reduced strength,
softening starts and the stress decreases. Consequently, the
material outside the weaker region must unload elastically,
because its strength has not been exhausted. This leads to the
conclusion that the size of the softening region cannot exceed
the size of the region with minimum strength. Such a region
can be arbitrarily small, and the corresponding softening

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 29

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 8: (a) Stress-strain diagram with linear softening, (b) fan of possible post-peak branches of the load-displacement diagram



branch can be arbitrarily close to the elastic branch of the
load-displacement diagram. Thus the standard strain-soften-
ing continuum formulation leads to a solution that has several
pathological features: (i) the softening region is infinitely
small; (ii) the load-displacement diagram always exhibits
snapback, independently of the structural size and of the
material ductility; (iii) the total amount of energy dissipated
during the failure process is zero.

From the mathematical point of view, these annoying
features are related to the so-called loss of ellipticity of the
governing differential equation. In the present one-dimen-
sional setting, loss of ellipticity occurs when the tangent
modulus ceases to be positive, i.e., it coincides with the onset
of softening (this is not always the case in multiple dimen-
sions). The boundary value problem becomes ill-posed, i.e.,
it does not have a unique solution with continuous depend-
ence on the given data. From the numerical point of view,
ill-posedness is manifested by pathological sensitivity of the
results to the size of finite elements. For example, suppose
that the bar is discretized by Ne twonode elements with linear
displacement interpolation. If the numerical algorithm prop-
erly captures the most localized solution, the softening region
extends over one element, and we have L L Ns e� . The slope
of the post-peak branch therefore strongly depends on the
number of elements, and it approaches the initial elastic slope
as the number of elements tends to infinity; see Fig. 9a,
constructed for a linear softening law with � �f 0 20� . The
strain profiles at u L� 2 0� for various mesh refinements are
plotted in Fig. 9b (under the assumption that the weakest spot
is located at the center of the bar). In the limit, the profiles
tend to 2 20L x L� �( )� where � denotes the Dirac distribu-
tion. The limit solution represents a displacement jump at the
center, with zero strain everywhere else.

4.2 Nonlocal formulation serving as
localization limiter

In real materials, inelastic processes typically localize in
narrow bands that initially have a small but nonzero width.
Propagation and coalescence of microdefects in the localiza-
tion band can eventually lead to the formation of a displace-
ment discontinuity, e.g., of a macroscopic stress-free crack or
a sharp slip line. The initial thickness of the localization band
depends on the material microstructure and is usually of
the same order of magnitude as the characteristic material
length, determined by the size or spacing of dominant het-

erogeneities. Therefore, it is natural to expect that enriched
continuum theories can better reflect the actual deformation
and failure processes and restore mathematical well-posed-
ness of the boundary value problem. Indeed, when properly
formulated, nonlocal or gradient enrichments regularize the
problem in the sense that the resulting strain field is highly
concentrated in certain narrow zones but remains continu-
ous. The corresponding numerical solutions converge upon
mesh refinement to a physically meaningful limit, and the
numerical results do not suffer by pathological sensitivity to
the discretization. Enrichments that prevent localization of
strain into arbitrarily small regions are called localization
limiters.

Nonlocal material models of the integral type were first
exploited as localization limiters in the 1980s. After some
preliminary formulations exploiting the concept of an imbri-
cate continuum (Bažant, Belytschko and Chang, 1984), the
nonlocal damage theory emerged (Pijaudier-Cabot and Ba-
žant, 1987). Nonlocal formulations were then developed for a
number of constitutive theories, including softening plastic-
ity, smeared cracking, microplane models, etc. For a list of
references, see e.g. Bažant and Jirásek (2002) or Chapter 26
in Jirásek and Bažant (2002).

Generally speaking, the nonlocal approach consists in
replacing a certain variable by its nonlocal counterpart ob-
tained by weighted averaging over a spatial neighborhood of
each point under consideration. In nonlocal elasticity, the
averaged quantity is usually the strain. Nonlocal elastic mod-
els can correctly reflect the experimentally observed disper-
sion of short elastic waves. However, in typical structural
applications, the strain in the elastic regime remains relatively
smoothly distributed (with the exception of stress concentra-
tions and singularities around specific points, e.g., tips of
pre-existing sharp cracks). Steep strain gradients appear only
after the onset of localization and are accompanied by highly
nonuniform distribution of damage. Therefore, most nonlo-
cal models serving as localization limiters reduce to standard
local elasticity at low strain levels, and the nonlocal effects are
considered only in the inelastic regime. For instance, one
widely used nonlocal damage formulation replaces the strain
in the loading-unloading conditions (79) by its nonlocal aver-
age, �, while strain entering the stress-strain law (77) is still
considered as local. According to the modified loading-un-
loading conditions,

� , , ( ) �� � � � � �� � � � �0 0 0, (80)

30 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 9: Pathological effects of mesh refinement on the numerical results obtained with the local damage model: mesh-dependence of
(a) load-displacement diagrams, (b) strain profiles



the internal variable � has the meaning of the largest previ-
ously reached value of nonlocal strain.

If strain has a tendency to localize in a very narrow region,
e.g., in one single cross section, the nonlocal strain becomes
high not only in this region but also in its close neighborhood.
This leads to damage growth in that neighborhood, and the
local strain at the damaged points must be increased in order
to keep the stress distribution uniform. By this mechanism,
strain and damage are prevented from localizing into a single
cross section. The localized region always has a certain finite
size, controlled by the length parameter that appears in the
definition of the nonlocal weight function.

4.3 Localization analysis
The initial bifurcation from a uniform state can be studied

analytically under the simplified assumption that the strain
keeps increasing at all points both for the fundamental uni-
form solution and the bifurcated nonuniform solution. If the
bifurcated solution is considered as a small perturbation of
the uniform solution, the stress perturbation �
 can be linked
to the strain perturbation �� by the linearized equation

�
 �� � �� � �( )1 � �E E (81)
and the perturbation of the damage field �� is linked to the
strain perturbation by the nonlocal relation

� � 
 �� 
 
�( ) ( , ) ( )*x g x

L

� � d , (82)

where g* is the derivative of the damage function g with
respect to its argument � evaluated for the fundamental
(uniform) solution and therefore independent of the spatial
coordinate.

Even though the fundamental solution is considered as
static, we will analyze the evolution of the perturbation as a
dynamic process. This approach provides more insight into
the localization phenomena. In dynamics, the stress pertur-
bation and the displacement perturbation must satisfy the
equation of motion

�� �
��u � � � 0. (83)
For an infinite body and nonlocal weight function in the

form � 
 � 
( , ) ( )x x� �0 , the assumed solution for �u in the
form of a harmonic wave substituted into (81) – (83) leads to
the dispersion equation

� � � � ��� � �2 2 2 01 0( ) ( )
* *� Ek g E k k , (84)

where � is the circular frequency of the wave (not to be con-
fused with the damage variable �), k is the wave number, and
�0

* is the Fourier image of the weight function. The resulting
dispersion law reads

� ��� � �c k g k0 01 �
* * ( ) . (85)

The Fourier image �0
* has a unit value at k � 0 and smaller

values at positive wave numbers k. So if 1 0� � �� g *� , all
wave numbers have real positive frequencies and a small per-
turbation of an arbitrary shape propagates through the body
but does not grow in magnitude. If 1 0� � �� g *� , there
exists a band of low wave numbers between zero and a positive
limit kcrit for which the frequencies are imaginary. This indi-
cates an instability if the perturbation contains a component
with a sufficiently long wavelength. In statics, a stationary

wave of wavelength 2� kcrit can be superimposed on the fun-
damental uniform solution without violating the equilibrium
condition.

The critical wave number can be determined from the
condition � � 0. In a monotonic loading process, the damage
variable � is a function of strain, � �� g( ), and also the deriva-
tive g g*� d d� is a function of strain. The expression under
the square root in (85) vanishes for the wave number satisfy-
ing the condition

�
�

�
�

�

0
1* ( )

( )
( )

k
g
gcrit

�
�
d

d

. (86)

For a given damage function g and nonlocal weight func-
tion �0, this equation can be solved (either analytically or nu-
merically) and the critical wavelength L kcrit crit� 2� can be
determined as a function of strain.

For instance, for an exponential softening curve (more re-
alistic for concrete than the linear one), the damage function
is given by

g
f

( ) exp�
�

�

� �

� �
� � �

�

�

�

�

�
�

�

�

�
�

1 0 0
0

(87)

where �0 is the limit elastic strain and � f is another parameter
controlling the slope of the softening curve. Substituting (87)
and the Fourier image �0

2 2 4* ( ) exp( )k k a� � of the Gauss-like
weight function (52) into (86), we obtain

L
k

a
crit

crit

f

� �

�
�

�

�

�
�

�

�

�
�

2

1
0

� �

�

� �
ln

. (88)

Of course, this expression is valid only in the inelastic
range, i.e., for � �� 0. For the present damage function (87),
the elastic limit coincides with the peak of the stress-strain
curve. If, for instance, � �f �11 0, the critical wavelength for
the uniform state at peak stress is L a acrit � �� ln . .11 10176 .
This example shows that the critical wavelength is propor-
tional to the model parameter that plays the role of the inter-
nal material length (and is often denoted as l).

In an infinite bar at peak stress, the appearance of a sta-
tionary wave of wavelength Lcrit corresponds to a particular
periodic localization pattern. In reality, there exists an ener-
getically more favorable localization pattern that is not peri-
odic but is concentrated into one single interval of length
slightly below the critical wavelength. The exact size of the
localized zone could be solved from a Fredholm integral
equation of the second kind combined with the complement-
arity conditions, but since the solution is not available in
closed form, we will not elaborate on that. If the bar is finite
but longer than Lcrit , it can be expected that at peak stress the
subsequent increments of strain localize into a band of thick-
ness approximately equal to Lcrit . For shorter bars, localiza-
tion can be delayed and the actual behavior depends on the
specific form of the nonlocal averaging operator around
the boundaries. However, in general it can be expected that
localization occurs when the critical wavelength becomes ap-
proximately equal to the bar length. As shown in Fig. 10a, the
critical wavelength monotonically decreases with increasing
strain and asymptotically tends to zero. This property has im-
portant implications for the evolution of the localized strain

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 31

Acta Polytechnica Vol. 44  No. 5 – 6/2004



profile – it indicates that the active zone in which strains are
increasing tends to shrink during the loading process. Such a
trend has indeed been confirmed by numerical simulations.

4.4 Mesh sensitivity and size effect
Fig. 10b shows the load-displacement diagram for strain

localization in a bar under uniaxial tension, calculated by
the finite element method using a nonlocal damage model
with the weight function (53) and with the exponential dam-
age function (87). As the number of elements Nel increases,
the load-displacement curve rapidly converges to the exact
solution. Convergence of strain and damage profiles gener-
ated by an applied displacement u L� 2 0� is documented in
Fig. 10c, d. In contrast to the local model, the process zone
does not shrink to a single point as the mesh is refined. Its size
is controlled by parameter that sets the internal length scale.
This example demonstrates that the nonlocal formulation
serves as a localization limiter and provides an objective de-
scription of the localization process, with no pathological sen-
sitivity of the numerical results to the discretization.

Another important advantage of nonlocal softening mod-
els is that they can realistically describe the size effect on nom-
inal strength of quasibrittle materials. The nominal strength
is understood as the peak load divided by a characteristic area
of the structure. According to the standard (local) version of

perfect plasticity theory, the nominal strength for a set of geo-
metrically similar structures of different sizes should be the
same, independent of the size. For instance, for a beam of a
rectangular cross section, subjected to three-point bending,
the plastic collapse load is

F
M

L
bD

L
p

0

2
0

4
� �



, (89)

where b is the width and D the depth of the cross section, L is
the span, Mp is the plastic limit moment of the cross section,
and 
0 is the yield stress. The nominal strength defined as the
peak load divided by the cross-sectional area,






N
F
bD

D
L

� �0 0 (90)

is equal to the material property 
0 multiplied by the geomet-
rical factor D L, which depends only on the shape of the
structure but not on its size (assuming proportional scaling of
all structural dimensions for three-dimensional similarity, or
at least of the in-plane dimensions L and D for two-dimen-
sional similarity). In contrast to elastic-perfectly plastic struc-
tures, structures made of quasibrittle materials often exhibit
a strong size dependence of the nominal strength. For certain
specimen geometries with initial notches scaled proportion-
ally to other structural dimensions, the large-size limit is
adequately described by linear elastic fracture mechanics,

32 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

(c) (d)

Fig. 10: Nonlocal damage model: (a) critical wavelength as a function of strain (for reference, the stress-strain curve is plotted by the
dashed line), (b) convergence of load-displacement diagram, (c) convergence of strain profile, (d) convergence of damage pro-
file; Nel � number of elements



which predicts proportionality of the nominal strength to
1 D. The transition between no size effect for small sizes
and strong size effect for large sizes can be approximated by
Bažant’s (1984) formula






N
N

D
D

�

�

0

0
1

, (91)

where D0 and 
N0 are parameters to be determined by fitting
of the experimental results.

The transitional size effect characteristic of quasibrittle
structures is nicely reproduced by nonlocal softening models,
if the model parameters are chosen correctly. As an example,
consider the compact tension test of concrete depicted in
Fig. 11a. The experimental results obtained by Wittmann,
Mihashi and Nomura (1990) are shown in the logarithmic
size effect diagram in Fig. 11b, along with the optimal fit by
formula (91) and the results of numerical simulations. The
role of the characteristic structural size D is played by the
ligament length, and the nominal strength is defined as the
peak load divided by the ligament area, bD, where b is the
out-of-plane thickness of the specimen. The experimentally
measured size effect can be accurately reproduced by a non-
local isotropic damage model with different combinations of
internal length parameter a and parameter �f of the damage
function (86). If one parameter is fixed, the other can be
determined by optimal fitting, but both parameters cannot
be determined simultaneously from the size effect on nominal
strength only. A unique parameter identification requires
additional information such as the distribution of strain or
damage in the process zone (Geers, de Borst, Brekelmans
and Peerlings, 1999) or the size effect on fracture energy
(Jirásek, Rolshoven and Grassl, 2004).

5 Concluding remarks
The common denominator of all examples presented in

the preceding sections is that the characteristic wavelength of
the deformation field becomes comparable to the characteris-
tic size of the internal material structure. Here, the notion of
characteristic wavelength has to be understood in a broad
sense, not only as the spatial period of a dynamic phenome-
non but also as the length on which the value of strain

changes substantially in static problems. Such a more general
definition could be based e.g. on a suitably normalized ratio
between the maximum strain and the maximum strain gradi-
ent (both in absolute values). Thus the characteristic wave-
length is necessarily close to the internal material length if the
size of the specimen is not much larger than the size and spac-
ing of major heterogeneities, or if strain localizes due to
softening.

The enrichment terms introduced by various generalized
continuum theories have a differential or integral character,
but all of them can be considered as nonlocal, at least in the
weak sense. They always introduce a model parameter with
the dimension of length, which reflects the internal length
scale of the material.

Three typical cases covered here encompass static and dy-
namic phenomena, linear and nonlinear behavior, and three
classes of material laws, namely elasticity, plasticity and con-
tinuum damage mechanics. This shows that the nonlocal en-
richments can be useful in a wide range of mechanical prob-
lems. Unfortunately, so far there is no general and universally
accepted theory covering this entire range within one unified
framework. Even though the first nonlocal theories were pio-
neered in the 1960s, there many problems remain open and
many issues unresolved. Some of the most challenging ques-
tions include the correct formulation of boundary conditions,
micromechanical justification of models with nonlocal inter-
nal variables, or identification techniques for the internal
length parameter and its possible evolution.

References

[1] Aifantis E. C.: “On the microstructural origin of certain
inelastic models.” Journal of Engineering Mechanics and
Technology, ASME, Vol. 106 (1984), p. 326–330.

[2] Aifantis E. C.: “Update on a class of gradient theories.”
Mechanics of Materials, Vol. 35 (2003), p. 259–280.

[3] Askes H., Metrikine A. V.: “One-dimensional dynami-
cally consistent gradient elasticity models derived from a
discrete microstructure. Part 2: Static and dynamic re-
sponse.” European Journal of Mechanics A, Vol. 21 (2002),
p. 573–588.

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 33

Acta Polytechnica Vol. 44  No. 5 – 6/2004

(a) (b)

Fig. 11: (a) Compact tension specimen, (b) size effect diagram



[4] Bažant Z. P.: “Size effect in blunt fracture: Concrete,
rock, metal.” Journal of Engineering Mechanics, ASCE,
Vol. 110 (1984), p. 518–535.

[5] Bažant Z. P., Chang, T.-P.: “Instability of nonlocal con-
tinuum and strain averaging.” Journal of Engineering Me-
chanics, ASCE, Vol. 110 (1984), p. 2015–2035.

[6] Bažant Z. P., Jirásek M.: “Nonlocal integral formulations
of plasticity and damage: Survey of progress.” Journal
of Engineering Mechanics, ASCE, Vol. 128 (2002),
p. 1119–1149.

[7] Bažant Z. P., Belytschko T. B., Chang, T.-P.: “Contin-
uum model for strain softening.” Journal of Engineering
Mechanics, ASCE, Vol. 110 (1984), p. 1666–1692.

[8] Eringen A. C.: “Linear theory of nonlocal elasticity and
dispersion of plane waves.” International Journal of Engi-
neering Science, Vol. 10 (1972), p. 425–435.

[9] Fish J., Chen W., Nagai, G.: “Nonlocal dispersive model
for wave propagation in heterogeneous media. Part 1:
One-dimensional case.” International Journal for Numeri-
cal Methods in Engineering, Vol. 54 (2002a), p. 331–346.

[10] Fish J., Chen W., Nagai G.: “Nonlocal dispersive model
for wave propagation in heterogeneous media. Part 2:
Multi-dimensional case.” International Journal for Numeri-
cal Methods in Engineering, Vol. 54 (2002b), p. 347–363.

[11] Fleck N. A., Hutchinson J. W.: “A reformulation of strain
gradient plasticity.” Journal of the Mechanics and Physics of
Solids, Vol. 49 (2001), p. 2245–2271.

[12] Fleck N. A., Muller G. M., Ashby M. F., Hutchinson
J. W.: “Strain gradient plasticity: theory and experi-
ment.” Acta Metallurgica et Materialia, Vol. 42 (1994),
p. 475–487.

[13] Geers M. G. D., de Borst R., Brekelmans W. A. M.,
Peerlings R. H. J.: “Validation and internal length scale
determination for a gradient damage model: Applica-
tion to short glass-fibre-reinforced polypropylene.” In-
ternational Journal of Solids and Structures, Vol. 36 (1999),
p. 2557–2583.

[14] Geers M. G. D., Engelen R. A. B., Ubachs R. J. M.: “On
the numerical modelling of ductile damage with an im-
plicit gradient-enhanced formulation.” Revue européenne
des éléments finis, Vol. 10 (2001), p. 173–191.

[15] Green A. E., Rivlin R. S.: “Simple force and stress multi-
poles.” Archive for Rational Mechanics and Analysis, Vol. 16
(1964), p. 325–353.

[16] Jirásek M., Bažant Z. P.: “Inelastic Analysis of Struc-
tures.” John Wiley and Sons, Chichester, U.K., 2002.

[17] Jirásek M., Rolshoven, Grassl, P.: “Size effect on fracture
energy induced by nonlocality.” International Journal for
Numerical and Analytical Methods in Engineering, Vol. 28
(2004), in press.

[18] Lakes R. S.: “Experimental microelasticity of two po-
rous solids.” International Journal of Solids and Structures,
Vol. 22 (1986), p. 55–63.

[19] Ma Q., Clarke D. R.: “Size dependent hardness in silver
single crystals.” Journal of Materials Research, Vol. 10
(1995), p. 853–863.

[20] Metrikine A. V., Askes H.: “One-dimensional dynami-
cally consistent gradient elasticity models derived from
a discrete microstructure. Part 1: Generic formula-
tion.” European Journal of Mechanics A, Vol. 21 (2002),
p. 555–572.

[21] Mindlin R. D.: “Micro-structure in linear elasticity.” Ar-
chive for Rational Mechanics and Analysis, Vol. 16 (1964),
p. 51–78.

[22] Mindlin R. D.: “Second gradient of strain and surface
tension in linear elasticity.” International Journal of Solids
and Structures, Vol. 1 (1965), p. 417–438.

[23] Morrison J. L. M.: “The yield of mild steel with particu-
lar reference to the effect of size of specimen.” Pro-
ceedings of the Institute of Mechanical Engineers, Vol. 142
(1939), p. 193–223.

[24] Nagai G., Fish J., Watanabe K.: “Stabilized nonlocal
model for wave propagation in heterogeneous media.”
Computational Mechanics, Vol. 33 (2004), p. 144–153.

[25] Nix W. D.: “Mechanical properties of thin films.” Metal-
lurgical Transactions, Vol. 20A (1989), p. 2217–2245.

[26] Pijaudier-Cabot G., Bažant Z. P.: “Nonlocal damage the-
ory.” Journal of Engineering Mechanics, ASCE, Vol. 113
(1987), p. 1512–1533.

[27] Poole W. J., Ashby M. F., Fleck N. A.: “Microhardness
of annealed and work-hardened copper polycrystals.”
Scripta Metallurgica et Materialia, Vol. 34 (1996),
p. 559–564.

[28] Rogula D.: “Introduction to nonlocal theory of material
media”, in D. Rogula (ed.), Nonlocal Theory of Material
Media, No. 268 in CISM Courses and Lectures, Springer
Verlag, Wien and New York, 1982, p. 125–222.

[29] Shu J. Y., Fleck N. A., van der Giessen E., Needleman A.:
“Boundary layers in constrained plastic flow: compari-
son of nonlocal and discrete dislocation platicity.” Jour-
nal of the Mechanics and Physics of Solids, Vol. 49 (2001),
p. 1361–1395.

[30] Stolken J. S., Evans A. G.: “A microbend test method for
measuring the plasticity length scale.” Acta Metallurgica et
Materialia, Vol. 46 (1998), p. 5109–5115.

[31] Toupin R. A.: “Elastic materials with couple-stresses.”
Archive for Rational Mechanics and Analysis, Vol. 11 (1962),
p. 385–414.

[32] Wittmann F. H., Mihashi H., Nomura N.: “Size effect on
fracture energy of concrete.” Engineering Fracture Me-
chanics, Vol. 35 (1990), p. 107–115.

[33] Yarnel J. L., Warren J. L., Koenig S. H.: In Lattice Dy-
namics, R. F. Wallis (ed.), Pergamon Press, 1965, p. 57.

Milan Jirásek

Laboratory of Structural and Continuum Mechanics
Swiss Federal Institute of Technology at Lausanne (EPFL)
1015 Lausanne, Switzerland

34 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 44  No. 5 – 6/2004