0045-7949(94)0050 I -x

COMPARISON OF NONLOCAL APPROACHES IN

CONTINUUM DAMAGE MECHANICS

J. H. P. de Vree, W. A. M. Brekelmans and M. A. J. van Gils

Faculty of Mechanical Engineering, University of Technology. P.O. Box 513. 5600 MB Eindhoven,
The Netherlands

(Received 4 December 1993)

Abstract-A local continuum damage theory and two distinct nonlocal variants are applied to model the
failure behaviour of a construction made of macroscopically brittle material. In the nonlocal formulations
a material characteristic length parameter is introduced associated with the width of the microstructural
damage zone. The numerical implementation of the approaches has been performed in a finite element
code. Simulation results calculated for a plane stress configuration are compared. The local approach
solutions show severe lack of mesh objectivity, whereas both the nonlocal solutions converged after mesh
refinement. By adequate tuning of the nonlocal descriptions mutually similar responses can be obtained,
although intrinsic differences are present in the resulting damage distributions.

1. INTRODUCTION . micropolar (Cosserat) continuum theory [4];

. viscous regularization [ 11, 121;
Continuum damage approaches offer the possibility
. local manipulation of material properties depend-
to simulate the mechanical behaviour of construc-
ing on the element size [I I].
tions of history-dependent material, irreversibly de-
Gradient dependent models as well as micropolar
generating under mechanical loads. The damage
continuum theory require radical modifications of
process is characterized by the development, growth
ordinary finite element codes. It is obvious that
and coalescence of microdefects leading to the for-
viscous regularization is inappropriate for typically
mation and propagation of macrocracks and eventu-
brittle material behaviour. The adaptation of ma-
ally to rupture. Microstructural observations are
terial properties is disputable from a physical point of
useful for the identification and quantification of
view: global responses may reasonably be predicted,
history parameters associated with damage [I]. The
however, calculated local strain and damage distri-
introduction of damage parameters is adopted from
butions will certainly not be in agreement with reality.
Kachanov 121, who first applied the damage concept
In this paper two different nonlocal theories are
to model fracture in creep. Although damage par-
examined and the major structure of the finite el-
ameters are conceptually related to the observable
ement implementation will be explained. The numeri-
local deterioration of the material, generally these
cal results of a characteristic example are reported
parameters are phenomenologically incorporated
and discussed in conclusion.
into the constitutive description. In this paper, brittle
material behaviour will be considered: elastic moduli
are influenced by damage development. The damage 2. THE LOCAL DAMAGE APPROACH
evolution is assumed to be isotropic; explicit time
dependence and thermal effects are excluded from the
considerations.
The increase of damage generally leads to local The influence of isotropic damage on elastic ma-
softening behaviour: the tangential stiffness becomes terial behaviour is described by the so-called damage
negative. It has been shown[3,4] that a local ap- parameter D. This monotonously increasing scalar
proach to softening phenomena may lead to a phys- quantity, 0 < D < 1. expresses the level of material
ically unacceptable localization of the deformation; in degradation [13-l 61. Undamaged material is charac-
a finite element context serious mesh sensitivity terized by D = 0, the complete loss of stiffness and
occurs [5-71. To overcome the deficiencies of the coherence by D = 1. According to the classical con-
classical local modelling a number of basically differ- tinuum damage mechanics approach the preposition
ent approaches can be distinguished: of the principle of effective stress [13, 141 leads to the
l nonlocal modelling of the constitutive behaviour constitutive relationship
14, 8, 91;
l gradient dependent material description [4. lo]; (T = (1 - D)Ht (1)

581
582 J. H. P. de Vree cl al.

with Q the column with Cauchy stress components, t For strain tensors mapped within this surface the
the column with strain components and with H damage rate is zero, while for increasing damage the
Hooke’s elasticity matrix. The product (I-D)H strain tensors should be mapped on the damage
can be considered as the effective elasticity matrix. surface. The damage development for brittle ma-
Implicitly it is assumed that Poisson’s ratio v is not terials is governed by a damage evolution law which
affected by damage. relates the value of the actual damage D to the actual
The damage evolution is determined by a scalar threshold parameter K. In Section 4 a particular
measure of the strain components: the damage equiv- specification of D = D(K) will be used for further
alent strain tcq = C,(E) 2 0. There are several alterna- elaborations.
tives [I 7, 181 to define ccq, weighting the strain
components differently. A suitable selection has to be 2.2. Aspects of the jinite element implementution
made in accordance with material features. In this An arbitrary structural problem is described by the
paper a modification (adopted from Ref. [19]) of the equilibrium requirements, constitutive relationships
Von Mises equivalent strain will be applied: and boundary conditions. This general problem for-
mulation usually only admits a solution by approxi-
mation. The finite element method is applied to
generate the system of discretized nodal equilibrium
equations. The nonlinearity, a consequence of the
1 incorporation of damage, necessitates an incremental
iterative solution procedure. The most recent estimate
+5
for the column with nodal displacements ‘IL,, , in
increment k, iteration i, is determined from
with y denoting the ratio of tensile and compressive
strength and with J, and J; strain invariants defined = ‘u, + “Au, =
kU,+, ‘u, + (“K,) “r! (7)
by
with kK, the global tangential stiffness matrix and ‘r,
J, = L,, + c,, + t,, the column with residual forces, resulting from the
previous iteration step. To calculate the matrix ‘K, the
local tangential stiffness matrix

- L,,C,, + 3(c $ + Et_ + Ei,),. (3)

(8)
The subscripts X, 18,z in eqn (3) indicate the (strain)
components in an orthonormal co-ordinate system.
Increase of damage is only possible if actually teq should be evaluated for the actual iterative approxi-
equals an evolving threshold value K >, 0 and in- mation at every intergration point. Based on the
creases. The damage remains constant if ccq is smaller constitutive modelling in Section 2.1 the matrix kS,
than K or if 6eq is not increasing. These statements are can be decomposed into
mathematically formulated by
“S, = “V, + “w, (9)
b < 0 if teq = K and if & > 0
with “V the elastic and ‘W, the softening contri-
i>=O if teq<ti or if C,,<O. (4) bution, kespectively,

Consistency requires that

“V, = (I - “D,)H, iw,= -
li-= (t’,,) if tcq = K

r;- = 0 if ccq < K (5)

To achieve a consistent iteration scheme the matrix

where the so-called McAuley brackets (.) are defined “W, should be determined according to the following
by elaboration:

(.X)=f(X+IXI). (6)

The nondecreasing threshold strain quantity K is

comparable with the yield stress in the theory of
elastoplasticity. Strain tensors satisfying c,~ = K rep- x h[“c,,, -h ‘K,] (11)
resent the actual damage surface in strain space.
 Comparison of nonlocal approaches in continuum damage mechanics 583

with h[ ] the Heaviside function defined by Analogous to eqn (4) the nonlocal damage quantity
D is assumed to be expressible in a weighted version

1
0 ifx,<O & of the equivalent strain:
h [x] = (12)
1 ifx>O
6 > 0 if Eeq= K and if & > 0
and with ‘K, the actual threshold value given by
6 = 0 if Fees< K or if & 6 0, (17)
‘K, = max[ht,q,, A ‘Kf] (13)
while the threshold K satisfies
where the subscript f refers to the final converged
state in the previous increment (k-l). It is remarked
that for kS, an almost zero matrix is substituted when
llc, approaches the critical value K,. associated with
D(K,) = I, to avoid numerical problems.
The column with residual forces ‘ri in the right to maintain a consistent formulation.
hand side of eqn (7) is equal to the difference of the 3.2. Aspects of the finite element implementation
external load “fe,, and the internal forces “A,,,, to be
determined from the actual stress state ka,: In a finite element context the weighted damage bp
in a particular integration point P can be calculated
from the nonlocal equivalent strain .$ (and the
‘r, = “fc,,,- “J,,, (ka, ). (14)
history) defined by
The iteration process is terminated when the above “IP
1
unbalance becomes negligibly small. -P -
c w’PAV%:; (19)
teq - nip
1 ,,,‘PA V’P ‘p= ’

<p= I
3. NONLOCAL DAMAGE FORMULATIONS
with ip the index of an integration point in a set of
3.1. General concepts P-surrounding integration points, nip the total num-
Two different nonlocal formulations are con- ber of integration points in that set, AVp the inte-
sidered with respect to the applicability in brittle gration volume associated with integration point ip
and 0 the weighting function to be evaluated at the
fracture simulations, particularly the regularization
distance r from point P to the integration point.
of deformation localization. Both procedures will be
Based on eqn (15), the iterative variation of the stress
explained in detail in Sections 3.3 and 3.4. Prelimi-
nary to that the common features are presented. column, necessary to determine the stiffness matrix,
The essential idea of the nonlocal approach is the can be formulated as
negligence of the generally accepted principle of local
6a = (1 - D)HGt - 66Hc
action, stating that, in the absence of temperature
effects, the stress in a material point is completely
determined by the deformation and the deformation
history at that point. The actual internal damage
parameter D figuring in the constitutive eqn (1) is
while SD can be expressed in St,, by
now assumed to be dependent on the strain (and
the strain history) in a limited finite area enclosing the
particular material point where the stress has to be (21)
evaluated. Consequently local strain peaks will
always have a certain transfer to the environment and
By incorporation of the relationships eqns (18) and
thus prevent the localization of the damage, which in
(19) the iterative variation of the nonlocal damage in
turn has a suppressing effect on the progressive
integration point P can be written as
growth of the deformation. The nonlocal constitutive
relationship for an arbitrary material point reads:

o=(l -b)He (15)

where b depends on the distribution of the strain and

the strain history in the vicinity of the material point. The partial derivative of the nonlocal equivalent
In a formal expression the quantity d at time t in a strain in point P to the local equivalent strain in point
point x in the surrounding volume V(.Y) can be ip reads:
written as:
5x4 .I. H. P. de Vree (‘I rrl

, Element , Cell @*V incrcmcnt k in iteration step i. comparable to the

relationships in eqn (10). arc found:

‘V, = (I -“6,)H,

Fig. I. Element mesh and cell grid

expressed in local quantities only (the superscript P.
denoting the particular integration point, is omitted
with in these equations).
U’,’

VP = 1 ,,.‘PA,,‘V’
(24)
3,’= I Among the variety of methods of averaging [6. IO],
in this section a strategy is pursued inspired by the
In cyn (22) the variation of damage in one inte- observation that, at microstructural level, many ma-
gration point P is coupled to the variation of the terials have a characteristic volume where the damage
strain in a certain number (IZ@) of integration points. distribution is almost uniform [21]. The dimension of
reflecting the nonlocal character of the modelling this characteristic volume is related to the size
procedure. The amount for this interaction in a of material inhomogeneities as the constituing
standard finite element schcmc. to achieve a conslst- particles or aggregates.
ent iterative solution strategy. requires scverc adup- The grid method to be presented here. is merely
tations in the computer code. Thercforc the solution adopted from Hall and Hayhurst [Zl]. However, it
technique has been chosen mconsistcntly. neglecting will be applied in a modified manner. A regular cell
the communication between difkrent integration grid overlay is placed over an clement mesh. Under
points. as far as the stiffness matrix is concerned. mesh refinement the grid remains unchanged. The
This implies that cqn (21) is replaced by the measures of the cells are equal to the characteristic
approximation: dimension of the material, indicated by the length
parameter i.. In Fig. I a grid with typical cells is
shown. Essential for the method is the assumption
that in every cell of the grid the nonlocal equivalent
strain and consequently the nonlocal damage is con-
sidercd as constant. The local equivalent strains of all
individual integration points within the cell are aver-
aged with uniform weighting. The assignment of
The following stiffness matrix contributions for elcmcnt integration points to the cells is performed
 Comparison of nonlocal approaches in continuum damage mechanics

Mesh 2: 107 elements

% +.C 5
Fig. 3. Stress-strain relationship.

only once at the beginning of the analysis and does

not change during load incrementation. As a conse-
quence of the procedure proposed, the calculated
Mesh 3: 163 elements Mesh 4: 328 elements
nonlocal damage will be discontinuous over the
boundaries of adjacent cells.

3.4. The continuous average strain method

The nonlocal equivalent strain in an arbitrary
integration point P can also be evaluated as the
weighted average of equivalent strains in all material
points within a certain radius. This method has been
described by Bazant et al. [8,22] and by Saourides
and Mazars [9]. This method principally leads to a
continuous damage field. The nonlocal equivalent Mesh 5: 504 elements

strain in a particular point P is determined with eqn Fig. 5. Meshes for the quarter plate
(I 9) by a summation over the (nip) integration points
with distance r to P not exceeding the material dimensional configurations; however, the back-
characteristic length i (with nip a known number, ground for this proposal is not quite clear.
depending on the local size of the elements in the
relevant vicinity of point P), while the weighting
function w is supposed to be a function of the 4. TEST PROBLEM AND MATERIAL MODEL

distance. For further elaborations w(r) has been

To examine the applicability of the methods out-
selected as a Gaussian weighting function according
lined above a relatively simple configuration model
to Bazant [8]:
with an inhomogeneous strain field is considered: a
w(r) = e-l”,“”
square plate (200 x 200 mm, thickness 1 mm) with a
(27) central circular hole (radius 50 mm). The plate is
supported at one side and loaded with a prescribed
where c( is a control parameter. Bazant and Pijaudier-
displacement at the opposite side. Because of sym-
Cabot [8] recommended LY= 2 in the case of two
metry the analysis is performed for the upper right
quarter of the plate, using boundary conditions as
indicated in Fig. 2. To ensure physical relevance the

=hkdl,

IN1 Mcdl2

=Medl3

=Meatld

Mcdh 5

Fig. 4. The damage D as a function of threshold K/K, Fig. 6. Force-displacement diagram. local modelling
586 J. H. P. de Vree rl ul

0.4 - 0.6

0.6 - 0.8

0.8 - 1.0

Mesh 2 Mesh 5
Fig. 7. Damage distribution, local modelling

material description has been adopted from related 5. RESULTS

research in Ref. [23]. The one dimensional
5.1 Locul modelling
stress-strain relationship (tension) is bilinear as illus-
trated in Fig. 3. This bilinear material behaviour The response results of the analysis with the local
matches with the local constitutive modelling in damage concept are presented in Fig. 6. The variable
Section 2 if D = D(K) is specified by along the horizontal axis is the prescribed displace-
ment of the upper edge of the quarter plate model and
D(ti)=, _lio ‘K’-ri) along the vertical axis the associated total tensile
(253) force. Typical differences in the distributions of the
K (K, -K,,)
damage, determined with a coarse and with a fine
with tiO and K, the initial and the critical values of K, mesh, arc visualized in Fig. 7 for a tensile force equal
respectively. Equation (28) is visualized in Fig. 4. The to 40 N in the post-peak softening regime. The objec-
material data [23] used for the numerical calculations tionable consequence of local modelling is obvious.
are: E =6000MPa. 1’ = 0.3. ;I = IO. tic,= 2 x 10m4, Mesh refinement leads to a vanishing energy dissipa-
K, = 2.7 x IO-’ and L - 10 mm. The nonlocal re- tion and the deformation localizes in an area with
lationship b = B(K), required for elaboration of decreasing volume.
the approaches in Section 3 is taken as equivalent to
eqn (28). 5.2. Nortlocrrl modelling
The quarter of the plate to be analyscd is subdi- Figure X displays the responses calculated with the
vided into isoparametric four-node plane stress el- grid method and the continuous average strain
ements. In order to study the mesh dependence, five method using Gaussian weighting with r = 2, see eqn
different meshes (Fig. 5) are applied to discrctizc the (27). It can be observed that for both methods the
displacement field. The clement size varies from results converge with refinement of the mesh. How-
large to small in the region of the plate where ever. these two nonlocal regularization procedures
high gradients of the strain and consequently of the predict quite different responses. In Figs 9 and 10 the
damage are to be expected. damage distributions are presented, again for a

Grid method Average strain method

INI WI
-

, -

Mesh 1; Mesh 2; Mesh 3; Mesh 4; - Mesh 5;

Fig. X. Force displacement diagrams, nonlocal modellmp.
 Comparison of nonlocal approaches in continuum damage mechanics 587

n
0.0 - 0.2

0.2 - 0.4

0.4 - 0.6

0.6 - 0.8

0.8 - 1.0

Mesh 2 Mesh 5
Fig. 9. Damage distribution, grid method.

0.0 - 0.2

0.2 - 0.4

0.4 - 0.6

0.6 - 0.8

0.8 - 1.0

Mesh 2 Mesh 5
Fig. IO. Damage distribution, average strain method.

tensile force of 40 N on the softening path. Especially a = 4. In Fig. 12 the damage distribution resulting
in Fig. 10 it can be observed that for a course mesh from the average strain method for 51= 4 is displayed.
the damage distribution already closely resembles the
distribution determined with the finest mesh. Based
6. CONCLUSIONS AND DISCUSSION
on the considerations of dissipated energy, the results
of the grid method can be accessed as rather satisfy- Continuum damage theories lead to powerful
ing; therefore these results are accepted as a reference methods for the numerical analysis of the behaviour
and it is examined for which value of the weighting of softening materials, provided that adequate pre-
parameter c( in eqn (27) the average strain method cautions are taken to limit the localization of the
produces a similar response to the grid method. deformation and consequently of the damage. The
Figure 11 gives the responses calculated with the problems of the mesh objectivity of finite element
average strain method for variations of the parameter solutions and correct representation of size effects are
a. Obviously the responses from the grid method and intrinsically related. With the application of dis-
the average strain method approximately coincide for tinguishable nonlocal approaches in continuum dam-
age mechanics, mesh objective mesh results were
70 7 ’
obtained in this paper for brittle material behaviour.

0.2 - 0.4

0.4 - 0.6

0.6 - 0.8

[mm1 0.8 - 1.0

LI =,; .---L1 =2. ~- -II =J. --~-* =4_
, -0 =5;

Fig. 11. Force-displacement diagrams for different values Fig. 12. Damage distribution with 1 = 4, average strain
of z. average strain method. method.
588 J. H. P. de Vree el u/

For the grid method and the weighted average strain analysis
of reinforced concrete beams and slabs. Proc,.
method convergence occurred with mesh refinement. Inst. Civil
Engr.r. 73, 135 I45 ( 1982).
4. R. de Borst and H. B. Miihlhaus. Continuum models
These methods can easily be implemented into a for discontinuous media. Proc,. ht. RILEM/ESIS
classical finite element computer code; the increase Con/: on Fracture Processes in Concrete, Rwk und
of CPU-time needed for the nonlocal calculations Cerumics (Edited by J. G. M. van Mier. J. G. Rots and
proposed is small. Manipulations in the continuous A. Bakker). pp. 601 -618. Noordwijk (1991).
5. S. T. Pietruszczak and 2. Mroz, Finite element analysis
weighted average strain method to produce appro-
of the deformation of strain-softening materials. It?/.J.
priate results, leads to a value of approximately 4 for numrr. Meth. Engng 17, 327-334 (19kl).
the weighting parameter CL This value is not in 6. Z. P. Bazant. F. B. Lin and G. Piiaudier-Cabot. Yield
accordance with the recommendation in the litera- limit degradation: nonlocal contindum model with local
strain. Proc. In/. Corzf: on Compututionul Plu.c-ticitj
ture [S], originating from mathematical consider-
(Edited bv E. Onate. R. Owen and E. Hinton).
ations. This deviation, here actually established for pp. 17%i779. Barcelona (1987).
a particular example, should be examined in a 7. A. Needleman, Material rate dependence and mesh
more general context in the continuation of this sensitivity in localization problems. Cun~pur. Mrth.
research. Appl. Mech. Engng 67, 69-85 (1988).
8. Z. P. Bazant and G. Pijaudier-Cabot, Nonlocal contin-
Apart from the availability of sufficient material
uum damage, localization instability and convergence.
data to describe the local softening phenomena, a J. appl. Mech. 55, 287-293 (1988).
number of suitable choices have to be made to 9. C. Saouridis and J. Mazars. A multiscale approach to
achieve a proper simulation of damaging configur- distributed damage and its usefulness for capturing
structural size effects. Proc. Frunce -C;.S. Workshop on
ations. Using the theory presented in this paper. the
Strain Localizution und Six Effects due to Cracking
following items have to be addressed carefully: and Damage (Edited by J. Mazars and Z. P. Bazant). pp.
. the definition of the equivalent strain; every 391-403, Cachan (1988).
material has its own characteristic internal disinte- IO. N. Triantafylidis and E. C. Aifantis. A gradient ap-
gration mechanism associated with some combi- preach to l&alization of deformation. I. kyperelasiic
materials. J. Elast. 16, 225-237 (1986).
nation of strain components. Therefore a variety
1 I. J. C. Simo. Strain softening and ‘dissipation: a unifica-
of equivalent strain definitions is known in the tion of approaches. Proc. France b’.S. Workshop on
literature leading to different simulation results. Sfruin Localization and Six Efects due to Crucking nnd
The correct choice to simulate experimental obser- Damage (Edited by J. Mazars and Z. P. Bazant),
pp. 440&461, Cachan (1988).
vation can only be accessed by hybrid numerical
12. R. de Borst, Gradient-dependent plasticity: formulation
and experimental research; and algorithmic aspects. Int. J. numer. Merh. Engng 35,
. the particular quantity in the nonlocal approach 521-539 (1992).
that will be averaged; in this paper averaging of 13. J. L. Chaboche, Continuum damage mechanics. Parts I
the equivalent strain has been elaborated. In the and 2. J. appt. Mech. 55, 59- 72 (1988).
14. J. C. Simo and J. W. Ju. Strain- and stress-based
literature also other choices have been made [21];
continuum damage models. Parts I and 2. In/. J. Solids
. the averaging method; two essentially different Struct. 23, 821-869 (1987).
averaging procedures have been considered in this 15. J. Lemaitre, Local approach of fracture. &gng Fruct.
paper, leading to a discontinuous and a continuous Mech. 25, 523-537 (1986).
16. L. M. Kachanov. Introduction to Continuum Dumugu
description of damage evolution. Microstructural
Mechanics. Martinus Nijhoff, Dordrecht (1986).
investigations should give evidence of the physical 17. L. J. M. G. Dortmans and G. de With. Weakest-link
relevance of one of these approaches for a particu- failure predictions for ceramics. Part IV: Application of
lar material. The weighted average strain method mixed mode fracture criteria for multi-axial loading.
offers a number of possibilities for tuning the J. eur. Ceram. Sot. 10, 109-l I4 (1992).
18. W. A. M. Brekelmans, P. J. G. Schreurs and J. H. P. de
modelling description on physical experience by
Vree, Continuum damage mechanics for softening of
adjusting the weighting function [24]. brittle materials. Actu Me&. 93, 133 -143 (1992).
19. J. G. Williams. Stre.~s Ana/[Link] of PoIymrrs. Longman,
Acknoll,ledRement-The results presented in this paper have London (1973).
been obtained by implementation of the theoretical develop- 20. 2. P. Bazant, Mechanics of distributed cracking. .4ppl.
ments into the DIANA finite element code (DIANA Analy- Mech. Rec. 39, 6655705 (1986).
sis by, Delft, the Netherlands). The authors gratefully 21. F. R. Hall and D. R. Hayhurst, Modelling of grain size
acknowledge D. G. Roddeman (from T.N.O. Building and effects in creep crack growth usmg a nonlocal contin-
Construction Research, Delft. The Netherlands) for his uum damage approach. Proc. R. Sot. Land. A 433,
support. 405-421 (1991).
22. Z. P. Bazant and F. B. Lin. Non-local yield limit
degradation. ht. J. numer. Me/h. Engng 26, 1805 1823
REFERENCES (1988).
23. J. Carmeliet and H. Hens, Frdcturc of a fabric re-
I. W. Becker and D. Gross, A one dimensional microme- inforced mortar based on a stochastic approach to
chanical model of elastic-microplastic damage initial damage. Proc. Con/: on Localizrd Domagr II:
evolution. Actu Mech. 70, 221-233 (1987). Fatigue und Frac,ture Mechunics (Edited by M. H.
2. L. M. Kachanov, Time of the rupture process under Aliabadi), pp. 283 -298 (1992).
creep conditions. TX). Akad. Nauk. SSR. Otd. Tekh. 8, 24. Z. P. Bazdnt and T. P. Chang. Instability of nonlocal
26-31 (1958). continuum and stram averaging. J. Engng Mwh. Dir.
3. M. A. Crisfield, Local instabilities in the non linear AXE 110, 1441 1450 (1984).