Nonlocal Theories in Continuum Mechanics

In standard continuum mechanics, a solid body is decomposed 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 isolated, but their interaction can take place only on the level of balance equations, through the exchange of mass, momentum, 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 denomination 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 deformation 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 generalized 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 theoretical predictions based on the standard theories, and also how the model performance can be improved by adding carefully selected enrichment terms. The enrichments to be discussed here are in general referred 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 integral-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 follow. 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 discrete 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, models with mixed spatial-temporal derivatives, and integral-type nonlocal elasticity. Advantages and drawbacks of individual formulations are discussed, and a general framework of nonlocal strain-gradient elasticity is outlined. Section 3 deals with size effects in microscale elastoplasticity, in particular with the size dependence of the apparent hardening modulus. Using the academic example of a


Introduction
In standard continuum mechanics, a solid body is decomposed 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 isolated, but their interaction can take place only on the level of balance equations, through the exchange of mass, momentum, 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 denomination 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 deformation 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 generalized 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 theoretical predictions based on the standard theories, and also how the model performance can be improved by adding carefully selected enrichment terms.
The enrichments to be discussed here are in general referred 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 integral-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 follow. 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 discrete 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, models with mixed spatial-temporal derivatives, and integral-type nonlocal elasticity. Advantages and drawbacks of individual formulations are discussed, and a general framework of nonlocal strain-gradient elasticity is outlined. Section 3 deals with size effects in microscale elastoplasticity, in particular with the size dependence of the apparent hardening modulus. Using the academic example of a 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 transitional size effect. Mathematical and numerical difficulties related to the objective description of strain localization due to softening are explained using the one-dimensional example 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 elements. 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 intrinsic 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.

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 r&& u Eu -¢¢ = 0, where r 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 r and E are constant coefficients, equation (1) where i is the imaginary unit, w is the circular frequency, k is the wave number, and c k = w is the wave velocity. Substituting (2) into (1), we find the condition - which implies that the magnitude of the circular frequency is proportional to the magnitude of the wave number. The signs of w and k are irrelevant -if both change, the real part of (2) does not change, and if one of them changes, the wave propagates in the opposite direction but otherwise remains the same. Therefore, we will restrict ourselves to nonnegative values and write the solution of (3) as Due to the direct proportionality between w and k, the velocity of a harmonic elastic wave, c k E = = w r , is constant, independent of the wave number k. Since a wave of a general shape can be represented by a linear combination of harmonic waves, even such a general wave propagates at constant velocity c and its shape remains invariant.
The situation is different in a discrete mechanical system, 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 where M is the mass of each mass point, K is the spring stiffness, u j is the displacement of mass point number j initially located at x j , and Z is the set of all integer numbers. For an assumed harmonic wave of the form we obtain the condition The spacing of mass points is assumed to be regular, i.e., . The circular frequency corresponding to wave number k is thus For this mass-spring system, the relationship between w 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 k p = w , and the group velocity, c k g = d d w . In the long wave limit (i.e., for k approaching zero), both the phase velocity and the group velocity tend to c a K M 0 = . As is clear from the w-k diagram in Fig. 1b, c 0 is the upper bound on the velocity of waves at any wave number. As k increases from 0 to 2p a, the corresponding phase velocity decreases from c 0 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 phenomenon is known as dispersion, the equation relating w and k is called the dispersion equation, the resulting function w( ) 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 r = 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 = = r 2 , which coincides with the initial slope c a K M 0 = of the nonlinear dispersion curve of the discrete model. Both dispersion curves are almost identical for long waves, but for wave numbers comparable to p 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 = r 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 (numerical) dispersion of short waves with wavelengths comparable to the element size. Note that the numerical dispersion disappears if equations (5) are integrated in time using the central difference scheme with time step D t a c = , which is just at the limit of numerical stability. Consequently, propagation of shocks and other wide-spectrum phenomena are not represented 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 K 2 , K 3 , …, K N , connecting pairs of mass points spaced by 2a, 3a, …, Na. A straightforward extension of the foregoing dispersion analysis yields the dispersion relation In the long-wave limit (k ® ¥), the phase velocity c An example of a dispersion curve constructed with N = 3 and K 1 = K 2 = K 3 = 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 2p a still corresponds to zero frequency. This is natural, because the displacements u j generated by a harmonic wave with wavelength a are the same at all the lattice sites. The associated mode is a uniform translation of the lattice points, which would not lead to any vibrations.
A striking property of the dispersion law (10) is its periodicity. 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 p 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 p 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 interval [0, w max ] where w max is a limiting circular frequency. If certain mass point is externally excited with a circular frequency larger than w max , the vibration does not propagate through the entire chain and remains localized to the neighborhood of the excited mass point. The discrete chain therefore acts as a low-band filter.
Dispersion of elastic waves in crystals is a real physical phenomenon 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 crystallographic direction (110), and it can be well approximated by a function of the form (10).

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 mathematically well-founded homogenization techniques providing 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 equivalence 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 internal 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 elegant and efficient alternative, it is possible to construct enrichments of the standard continuum theory that reflect the main features of the microstructure without using fast oscillating 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 deformation 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 where ¢ = u u x d d is the strain, further denoted as e, and L is the interval representing geometrically the one-dimensional body. In linear elasticity, the elastic energy density 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, ¢ e , further denoted as h. 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 where C is a higher-order elastic modulus. The variation of elastic energy density is given by where s ¶ ¶e e = = w E is the (Cauchy) stress and c ¶ ¶h h = = 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 where b is the body force density. In dynamics, b is replaced by the inertial force density, -r&& u. Combining this with the constitutive equations s e = E and c h = C and with the kinematic equations e = ¢ u and h = ¢¢ u , we obtain the wave equation of strain-gradient elasticity, 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 from which The dispersion curve is plotted in Fig. 3a. The phase velocity 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 permits to generalize certain uniqueness theorems known from standard elasticity. However, for C > 0 the phase velocity c p 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 analytical 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-order 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 l crit =1 vanishes. This means that the equation of motion (17) is satisfied by a stationary wave of wavelength 2 2 p p k l crit = . 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 aggravated by the fact that, for wave numbers exceeding k crit , the circular frequency w solved from the dispersion equation (18) becomes imaginary. This means that harmonic modes with wavelengths shorter than 2pl would spontaneously blow up. The source of this instability becomes apparent if one realizes that, for short waves, the negative higher-order part the elastic energy, -¢¢ El u 2 2 ( ) , exceeds in magnitude the positive standard part, E u ( ) ¢ 2 , and so the total energy density becomes negative. If the amplitude of the wave grows, energy is released instead of being consumed.

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 2pl are involved. If the body of interest is discretized by finite elements, the minimum wavelength 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 parameter l, the numerical solution leads to reasonable results. However, upon mesh refinement, the solution becomes polluted 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 rephrased 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 ¢¢ » r&&. Differentiating this twice with respect to x, we obtain Eu u IV » ¢¢ r&& . Finally, replacing in (17) u IV by ( )&& r E u¢¢ and C by -El 2 yields a modified wave equation  (17) at low wave numbers but a different asymptotic behavior for high wave numbers. Indeed, the usual procedure leads to the dispersion equation 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 monotonically increasing, has a negative curvature, and for k ® ¥ approaches a horizontal asymptote at circular frequency w = c l 0 . 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 reasoning and arrived at an equation of motion with two enrichments terms, proportional to &&¢¢ u and u IV . In its most general form, this equation can be written as where l is an internal length and b is an additional model parameter. to the microstructure and then proposed a parameter identification procedure based on a reflection-transmission test . Of course, one can also consider l and b as free model parameters and determine them by optimal fitting of the dispersion curve for a given material. The dispersion equation corresponding to (24), yields the phase velocity where c E 0 = r, as usual. If parameter b is nonnegative, the phase velocity remains real and positive for all wave numbers. For 0 < b < 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 C 1 -continuous finite element approximation (which is hard to construct on general meshes in multiple dimensions) or a mixed approach with independent approximations of several 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.

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 where E x ( , ) x is a function describing the generalized elastic modulus. The variation of elastic energy is evaluated as This can be written in the usual form d if the stress is defined as where 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, ¢ + = s b 0. Consequently, the wave equation for this model reads Since function E x s ( , ) x reflects the strength of long-distance interaction between points x and x, its value can be expected to be negligible if the distance between x and x is large compared to the internal length of the material (which corresponds to the characteristic size and spacing of major heterogeneities). For functions E s with a sufficiently fast decay, the integrals in (29) and (31) make sense even if the integration domain L is considered as the entire real axis. If the body infinite and macroscopically homogeneous, function E x s ( , ) x should depend only on the distance between x and x. Bearing in mind these restrictive assumptions, we present the modulus function in the form where E 0 is a reference value of the elastic modulus and a 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: Substituting the assumed harmonic form of an elastic wave (2) into the transformed wave equation we obtain the dispersion equation is the Fourier image of the nonlocal weight function a 0 ( ) r . Finally, the phase velocity is evaluated as 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 where c k p ( ) 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 c 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 k p ( ) is given by (20) and the integral in (38) does not converge (independently of the sign of C). On the other side, for the good Boussinesq problem (21) 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 conditions of boundedness (which play the role of boundary conditions at plus and minus infinity) can be expressed as where G x ( , ) x is the Green function of the Helmholtz equation, formally obtained as the solution of this equation with the Dirac distribution d x ( ) on the right-hand side. It turns out that the Green function is in this case given by and so equation (41) is in fact equivalent with (34) if the nonlocal weight function a 0 is selected according to formula (39).

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 elastic model derived from the potential (27). Still, using the alternative 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 and "solving" for the acceleration, we obtain 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 Taking the variation and using the same line of reasoning as for the basic version of nonlocal elasticity, we identify the constitutive equations 2 is the symmetrized higher-order modulus. Substituting this into the equation of motion of the strain-gradient theory, r s c && ( ) u --¢ ¢ = 0, and using the kinematic relations e = ¢ u and h = ¢¢ u , we obtain the wave equation that generalizes equation (31). For the moduli functions in the special form dispersion analysis provides the following expression for the phase velocity: Here, a 1 * denotes the Fourier image of function a 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 a 0 and a 1 are summarized in Table 1. The model permits to reproduce 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 w(k) is bounded, it is sufficient to use a regular weight function a 0 obtained by the inverse Fourier transform of ( ) w c k 0 2 . If w(k) tends to infinity but remains of order O(k), it is possible to use either a weight function a 0 with a singular part of the Dirac type, or regular functions a 0 and a 1 . Finally, if w(k) grows superlinearly but remains of order O(k 2 ), function a 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 gradients, 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.

Nonlocal model reproducing dispersion of discrete lattice
It is interesting that the nonlocal elastic model with a weight function a 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    which are plotted in Fig. 3b. The dispersion law of a mass--spring model with long-distance interactions can be reproduced 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 2p a the frequency 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, p/a], [2p/a, 3p/a], [4p/a, 5p/a], etc. Between these bands, the frequency 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 . 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 harmonic functions of wavelengths 2a, 2 2 a , 2 3 a , etc., which correspond to zero frequencies. Thus the static solution can be modified by a periodic function with period 2a without disturbing equilibrium. The dispersion analysis indicates that in dynamics the situation is even worse, because harmonic modes with wave numbers in the intervals (p/a, 2p/a), (3p/a, 4p/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 is positive only for wave numbers smaller than 5.763/a. Consequently, instabilities could develop for fine meshes with element size in the order of a.
The relationship between atomic lattices and the nonlocal 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 complete inverse Fourier transform but integrated (38) only in the limits -£ £ p p 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 x a x a x a x a x a é ë ê + - ) 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.

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 stiffness 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 torsional stiffness is evaluated experimentally as the ratio between the torque and the relative twist angle, it turns out that the expected proportionality to D 3 holds with sufficient accuracy 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 underestimated. In the context of the standard theory, this can be interpreted 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 comparable 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 elasticity, 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 certain parameter normally considered as a material property appears to be dependent on the size of the sample or specimen 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., 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.

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 spatial 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 decoupled. If this is the case, the relevant components of stress and strain are the shear stress s t xy º and the engineering shear strain 2e g xy º . Both of them are considered as independent of spatial coordinates y and z. From the equilibrium equation with vanishing body forces, it follows that t is also independent of x, i.e., it is uniformly distributed in space and varies only as a function of the pseudo-time parameterizing the loading process.
Classical plasticity with linear isotropic hardening is described by the basic equations where g p is the plastic strain and k is the cumulative plastic strain. The yield function f is given by f Y ( , ) ( ) t k t t k = -(58) where t Y is the current yield stress in shear, evaluated from the hardening law with t 0 = initial yield stress in shear and H = hardening modulus (in shear). During monotonic loading with positive value of shear stress t, there is no difference between the plastic shear strain g p and the cumulative plastic strain k. We will therefore replace g p by k and call it simply the plastic strain.
If the hardening modulus is positive and the loading is monotonic, the stress t uniquely determines the corresponding strain g. Since the stress is uniform, the strain must be uniform 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, Dv v L v L = --( ) ( ) 2 2, and the layer thickness, L. The stress-strain curve can be directly 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 standard plasticity model does not indicate any size effect.

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 surrounding points. This can be taken into account by sophisticated models that consider the details of the hardening mechanisms, e.g., by discrete dislocation models. As a simpler alternative, one can use enriched continuum models that take into account the micromechanical processes only "on the average", but by terms of a higher order than in the standard continuum theory. Motivated by certain micromechanical considerations, 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-dimensional setting and for linear hardening, the simplest version of the Aifantis gradient plasticity model replaces the hardening law (59) by 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 ordinary differential equation for the unknown function k. This equation should be supplemented 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 k 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 dislocation 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 k must vanish at the boundary. For homogeneous Dirichlet boundary conditions, k( ) ± = L 2 0, the solution of (61) is The distribution of plastic strain (normalized by the factor ( ) t t -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 boundary layers. With decreasing structural size L, the relative importance 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 relative tangential displacement of one loading plate with respect to the other can be evaluated as Defining the average shear straing = Dv L, we can transform (63) into the average stress-strain law where g t 0 0 = G is the limit elastic shear strain and is the elastoplastic (tangent) shear modulus. According to (64)-(65), the overall response of the specimen 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 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 gradient 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 theoretical 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 theory 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 general case.
Ideally, the enriched theory should be verified by several tests leading to different types of stress and strain distributions, 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 reasonably 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 modulusH 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 modulusH evaluated from the 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 apparent 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 hardening 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 D 3 , on the normalized angle of twist, fD, 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 sufficiently 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.
The particular solution depends on the boundary conditions, which should be formulated in terms of the nonlocal cumulative plastic strain and its normal derivative. For homogeneous Neumann boundary conditions, ¢ ± = k ( ) 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 homogeneous Neumann boundary conditions ¢ ± = k ( ) L 2 0. A stiff boundary with a perfect bond inhibits propagation of voids, which can be expressed by the condition that k vanishes at the boundary. With this homogeneous Dirichlet condition applied at both parts of the boundary, the particular solution of (70) becomes and the corresponding local cumulative plastic strain rate is 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 where~G is the elastoplastic shear modulus and is the size-dependent hardening modulus. For very large sizes L, the hardening modulus tends to Ha 2 , which is the value corresponding to the local theory with w p computed from k instead of k. 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 predicts a qualitatively similar trend as the explicit version of gradient plasticity.
The ratio between the size-dependent hardening modulus and its large-size limit is plotted in Fig. 7a as a function of the relative size, L l, for parameter a = 05 . . 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 boundary 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 initial distribution of the plastic strain rate and later the shape of the profile would change.

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 phenomenon, called strain localization, can be caused for instance 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 evolving and the approximate description of the material as a macroscopically homogeneous one may become questionable. Indeed, softening incorporated into standard inelastic continuum models leads to serious mathematical and numerical problems, and enriched theories are needed to provide an objective description of the softening process.
The essence of the localization problem will be explained using a one-dimensional example, which could be interpreted as localization of shear strain in a layer of elastoplastic material between two rigid plates (already studied in the previous section in the context of size effects). However, to better show various facets of nonlocal theories and their broad application 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 described 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 where s and e are the (normal) stress and strain, E is Young's modulus of elasticity, and W 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 Ee is reduced by the integrity factor 1 -W. The limit value W =1 corresponds to a fully damaged material that can no longer carry any stress. The growth of damage must be described by an appropriate damage evolution law for the internal variable W. This law could be postulated in the rate form, but a particularly simple and practical formulation is obtained with a damage law in the total form where k corresponds to the maximum level of strain reached in the previous history of the material. Mathematically, the internal variable k is defined by the loading-unloading conditions During monotonic loading, k is equal to the strain e, 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 e 0 , after which stress decreases as a linear function of strain until the zero stress level is reached at strain e 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 material, 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) f t is attained, the uniqueness of the response is lost. Stress must still remain uniform and it can only decrease, 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 distributions of the strain increments that lead to the same uniform stress decrement and thus represent a valid solution satisfying the equilibrium equation and the constitutive law. The compatibility equations do not represent any important constraint 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 L s and unloading everywhere else. When the stress is completely relaxed to zero, the strain in the softening region is e f and the strain in the unloading region vanishes; the total elongation of the bar is therefore u L s f = e . The length L s remains 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-displacement 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 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 e e f 0 20 = . The strain profiles at u L = 2 0 e 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 2 where d denotes the Dirac distribution. The limit solution represents a displacement jump at the center, with zero strain everywhere else.

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 localization band can eventually lead to the formation of a displacement 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-posedness 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 continuous. 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 imbricate 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 plasticity, 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 obtained by weighted averaging over a spatial neighborhood of each point under consideration. In nonlocal elasticity, the averaged quantity is usually the strain. Nonlocal elastic models can correctly reflect the experimentally observed dispersion of short elastic waves. However, in typical structural applications, the strain in the elastic regime remains relatively smoothly distributed (with the exception of stress concentrations 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 nonlocal 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 average, e, while strain entering the stress-strain law (77) is still considered as local. According to the modified loading-unloading conditions, (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 k has the meaning of the largest previously 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.

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 uniform solution and the bifurcated nonuniform solution. If the bifurcated solution is considered as a small perturbation of the uniform solution, the stress perturbation ds can be linked to the strain perturbation de by the linearized equation and the perturbation of the damage field dW is linked to the strain perturbation by the nonlocal relation where g * is the derivative of the damage function g with respect to its argument k 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 perturbation and the displacement perturbation must satisfy the equation of motion rd ds && u -¢ = 0.
For an infinite body and nonlocal weight function in the form a x a x ( , ) , the assumed solution for du in the form of a harmonic wave substituted into (81)-(83) leads to the dispersion equation where w is the circular frequency of the wave (not to be confused with the damage variable W), k is the wave number, and a 0 * is the Fourier image of the weight function. The resulting dispersion law reads The Fourier image a 0 * has a unit value at k = 0 and smaller values at positive wave numbers k. So if 1 0 --> W g * e , all wave numbers have real positive frequencies and a small perturbation of an arbitrary shape propagates through the body but does not grow in magnitude. If 1 0 --< W g * e , there exists a band of low wave numbers between zero and a positive limit k crit for which the frequencies are imaginary. This indicates an instability if the perturbation contains a component with a sufficiently long wavelength. In statics, a stationary wave of wavelength 2p k crit can be superimposed on the fundamental uniform solution without violating the equilibrium condition.
The critical wave number can be determined from the condition w = 0. In a monotonic loading process, the damage variable W is a function of strain,W e = g( ), and also the derivative g g * = d dw is a function of strain. The expression under the square root in (85) For a given damage function g and nonlocal weight function a 0 , this equation can be solved (either analytically or numerically) and the critical wavelength L k crit crit = 2p can be determined as a function of strain.
For instance, for an exponential softening curve (more realistic for concrete than the linear one), the damage function is given by (88) Of course, this expression is valid only in the inelastic range, i.e., for e e ³ 0 . For the present damage function (87), the elastic limit coincides with the peak of the stress-strain curve. If, for instance, e e f =11 0 , the critical wavelength for the uniform state at peak stress is L a a crit = = p ln . . 11 10176 . This example shows that the critical wavelength is proportional to the model parameter that plays the role of the internal material length (and is often denoted as l).
In an infinite bar at peak stress, the appearance of a stationary wave of wavelength L crit corresponds to a particular periodic localization pattern. In reality, there exists an energetically more favorable localization pattern that is not periodic 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 complementarity 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 L crit , it can be expected that at peak stress the subsequent increments of strain localize into a band of thickness approximately equal to L crit . For shorter bars, localization 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 approximately 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 important implications for the evolution of the localized strain 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. 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 damage function (87). As the number of elements N el increases, the load-displacement curve rapidly converges to the exact solution. Convergence of strain and damage profiles generated by an applied displacement u L = 2 0 e 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 description of the localization process, with no pathological sensitivity of the numerical results to the discretization.

Mesh sensitivity and size effect
Another important advantage of nonlocal softening models is that they can realistically describe the size effect on nominal 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 geometrically 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 where b is the width and D the depth of the cross section, L is the span, M p is the plastic limit moment of the cross section, and s 0 is the yield stress. The nominal strength defined as the peak load divided by the cross-sectional area, is equal to the material property s 0 multiplied by the geometrical 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-dimensional similarity). In contrast to elastic-perfectly plastic structures, structures made of quasibrittle materials often exhibit a strong size dependence of the nominal strength. For certain specimen geometries with initial notches scaled proportionally to other structural dimensions, the large-size limit is adequately described by linear elastic fracture mechanics, where D 0 and s 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 nonlocal isotropic damage model with different combinations of internal length parameter a and parameter e 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).

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 characteristic 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 phenomenon 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 gradient (both in absolute values). Thus the characteristic wavelength is necessarily close to the internal material length if the size of the specimen is not much larger than the size and spacing 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 dynamic phenomena, linear and nonlinear behavior, and three classes of material laws, namely elasticity, plasticity and continuum damage mechanics. This shows that the nonlocal enrichments can be useful in a wide range of mechanical problems. 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 pioneered in the 1960s, there many problems remain open and many issues unresolved. Some of the most challenging questions include the correct formulation of boundary conditions, micromechanical justification of models with nonlocal internal variables, or identification techniques for the internal length parameter and its possible evolution.