Computers and Structures 151 (2015) 73–85

Contents lists available at ScienceDirect

Computers and Structures

journal homepage: www.elsevier.com/locate/compstruc

Application of damage–plasticity models in finite element analysis

of punching shear
Adam Wosatko a,⇑, Jerzy Pamin a, Maria Anna Polak b,1
a
Institute for Computational Civil Engineering, Faculty of Civil Engineering, Cracow University of Technology, Warszawska 24, 31-155 Cracow, Poland
b
Department of Civil and Environmental Engineering, University of Waterloo, 200 University Avenue West, Waterloo, Ontario N2L 3G1, Canada

a r t i c l e i n f o a b s t r a c t

Article history: The paper presents numerical simulations of punching shear in a reinforced concrete slab-column con-
Received 23 December 2013 figuration formerly tested in the laboratory. A brief description of the test program at the University of
Accepted 11 January 2015 Waterloo is reported. For the simulation, a symmetric quarter of the test configuration is employed. Full
three-dimensional finite element discretized geometry is considered together with elastic–plastic
reinforcement embedded as truss elements in concrete. Two regularized numerical models of concrete,
Keywords: formulated within elastic–plastic-damage theories, are applied. The first one, called gradient damage,
FEM simulation
is refined by an additional averaging equation where gradient enhancement involves an internal length
Concrete
Punching shear
scale. In the second one, called rate-dependent damaged plasticity model, a viscoplastic strain rate is
Damage introduced. The results for the first model implemented in FEAP and the second model from ABAQUS
Plasticity are discussed in detail. Different issues of numerical simulation to properly predict punching shear
Regularization behaviour in the slab-column configuration are presented.
Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction description of this phenomenon, some experimental results and

analytical models can be found, in e.g. [1,2].
In multi-storey buildings, designers commonly utilize flat slabs Punching failure of reinforced concrete slabs without shear
on column systems instead of slab-girder-column combinations. reinforcement is brittle and thus dangerous, possibly leading to a
The advantages of flat floors are well-known; they include the progressive collapse of the structure. There are several methods
reduced height of buildings, and more economical design and per- to avoid such failures in flat concrete slabs with shear reinforce-
formance. However, the problem of punching shear failure occurs ment, as described e.g., by [3,4]. Using transverse steel reinforcing
in flat slabs. elements is one of the most effective methods to prevent punching.
When a typical internal slab-column connection is considered, These elements transfer tensile stresses and avoid opening of the
visible cracks usually form first due to bending of a slab, on its sur- punching shear crack. Properly placed and anchored punching
face under tensile stresses, in radial directions from the column. On shear reinforcing elements increase the strength and ductility of
the other hand, punching shear cracks are caused by three-dimen- the slab-column connections.
sional stress conditions with the concentration of shear stresses in Testing of flat concrete slabs for punching shear has been per-
the vicinity of the column. These cracks are formed inside the slab formed in several structural laboratories with many tests done in
and are not visible. As the load (and stresses) increases these cracks the 1950s [5–7]. These tests provided the basis for the calibration
open especially if no transverse reinforcing elements are present; of the ACI design provisions that appeared in the 1963 version of
these internal cracks will eventually reach the surface, join the the code. Since then, many more tests, e.g. [8–10], were done on
flexural cracks, and form a conical failure surface. At failure, crack- slabs in punching shear but little changes were implemented in
ing along circular lines is observed directly around the column on the design formulas. The debate continues on the methods of
the compression surface of the slab and, at a certain distance accounting for concrete contribution to shear carrying capacity,
around the column, on the tensile surface of the slab. The effects of unbalanced moments at the column and of longitudinal
reinforcement on punching strength, as well as the methods and
effectiveness of shear reinforcement in slabs. Although more test-
⇑ Corresponding author.
ing is still needed, not all of these issues can be rationally and
E-mail addresses: awosatko@L5.pk.edu.pl (A. Wosatko), jpamin@L5.pk.edu.pl
(J. Pamin), polak@uwaterloo.ca (M.A. Polak). objectively addressed through testing which is time, space and
1
Principal corresponding author. money consuming. Properly calibrated nonlinear finite element

http://dx.doi.org/10.1016/j.compstruc.2015.01.008
0045-7949/Ó 2015 Elsevier Ltd. All rights reserved.
74 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85

studies can be a useful source of information on slabs behaviour, The specimen in the test frame is shown in Fig. 1(a). The flexural
complementing experimental investigations. Moreover, such mod- reinforcement was formed by 10 M bars of nominal cross-section
els can be used for extensive parametric studies addressing differ- area Ar ¼ 100 mm2. The bars in the tension mat had the spacing
ent aspects of punching shear effects in flat RC slabs. of 100 mm and 90 mm for the upper and lower orthogonal layers
This paper outlines finite element studies on punching shear in in order to produce almost identical bending capacities in the
slabs using damage–plasticity finite element (FE) models. The two directions. The reinforcement bars on the compression side
presented simulations are performed for the slabs tested by formed a grid with 200 mm spacing. The column segments were
[11] with the goal to provide adequate predictions of the slabs reinforced with four 20 M bars.
behaviour, calibrate the models and include rational discussion The full test program consisted of six specimens, two of which
on the material parameters influencing the FE predictions. The had openings constructed near the column. The control slab SB1
problem of punching shear is not trivial to simulate, since the had no shear reinforcement and failed by punching shear in brittle
interaction between flexural and shear failure, while localized manner. All other specimens were strengthened using an increas-
fracture zones in a slab evolve, needs to be addressed. Some ing number of 9.5 mm diameter bolts symmetrically placed in con-
numerical simulations of punching shear are presented in [12– centric rows around the column. They showed substantial ductility
16]. The simulations done herein utilize two constitutive formula- and failed in flexure or combination of flexure and shear., cf. [11].
tions for concrete; namely gradient enhanced damage coupled to However, for all specimens, the flexural cracks initiated at column
plasticity and so-called damaged plasticity, which is available in corners and propagated radially towards the slab edges. In this
commercial FE program ABAQUS [17] and allows an extension paper only slab SB1, without shear reinforcement, is used for the
to viscosity-dependence. For the enhanced damage–plasticity, FE studies. The concrete compressive strength was 44 MPa and
the code developed by [18] is applied. This damage model which the measured tensile strength was 2.13 MPa. Flexural reinforce-
had first been proposed in [19] was developed and implemented ment yield strength was 455 MPa.
into the open-source code FEAP [20]. Although the problem of The results in terms of deflections, strains and crack widths
proper simulation of the localized failure seems to be not severe were monitored. The experimental load–deflection diagram for
in reinforced structures, it has been decided to employ a regular- SB1 is reproduced further in figures together with numerical
ized continuum description to minimize the effects of pathologi- results. The final experimental crack pattern, as seen from the ten-
cal mesh sensitivity and numerical instabilities. Continuum sion side, is presented in Fig. 1(b).
models contain various regularization methods, which are also
called localization limiters [21]. Such a limiter can be included
in many ways, for example: a certain variable averaged by a gra- 3. Review of employed constitutive models
dient operator is incorporated in the formulation [19,22–27]; or
an additional rate-dependent term like viscosity is enclosed in Two inelastic constitutive models have been used in the sim-
the constitutive equation [28–33]. ulations. They both assume that concrete is initially an isotropic
In the presented analysis linear kinematics is assumed. elastic material, described by elastic constants: Young modulus E
and Poisson ratio m. The oldest approach to fracture modelling
called smeared cracking originated in [34] and was further
2. Test program at the University of Waterloo developed for concrete fracture simulations in the eighties, for
an overview of those efforts the reader is referred e.g. to [35–37].
The slab specimens analyzed in this work were tested by [11]. In this paper the authors decided to employ two versions of
These were full-scale models representing interior slab-column coupled inelastic models involving damage and plasticity. They
connections with the column stub of 150  150 mm cross section are both based on the concept of effective stress acting on the
and simply supported along the edges with corners restricted from undamaged skeleton of the material and involve permanent strains
lifting. The overall dimensions of slab-specimens were when a yield limit is reached.
1800  1800  120 mm and they were supported along the The first approach is the gradient-enhanced damage–plasticity
1500  1500 square perimeter placed on neoprene pads. The spec- model described in [38,39] and implemented by the authors in
imens represented portions of a slab-column continuous system, FEAP. This model involves a simplified representation of cracked
bounded by the lines of contraflexure around the column. concrete based on a single scalar damage parameter, but it is non-
The specimens were loaded downwards through the column local and removes the pathological discretization sensitivity com-
until failure. Note that the experimental configuration was in an monly encountered when standard continuum models are used
upside down position in comparison with the real structural case. to simulate softening related to damage. The model incorporates

(a) Experimental setup. (b) Crack pattern for SB1.

Fig. 1. RC slab-column connection. Experimental setup and final crack patterns for slabs SB1 visible from the tension side. Photos quoted from [11].
 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85 75

an internal length parameter and therefore avoids the ellipticity r^ ¼ De e ð4Þ

loss due to material instability. e
where D is the elastic stiffness operator. Combining Eqs. (1) and (4)
The second approach, used in ABAQUS, is originally called
we obtain the relation:
’’damaged plasticity model for concrete and other quasi-brittle
materials’’. It is based on work by [40,41] and designed especially r ¼ ð1  xÞDe e ð5Þ
to reproduce the response of concrete to cyclic loads (which is
in which the damage measure x accounts for the degradation of the
not used in this paper). It provides a good representation of various
elastic stiffness. Damage does not grow during unloading, so that
features of concrete behaviour and a set of model parameters
for a pure damage model the constant (secant) stiffness ð1  xÞDe
allowing the user to adjust it in a rather flexible manner. It is
equipped with a viscoplastic enhancement, which is a natural reg- results in unloading to the origin, i.e. no residual strains remain
after unloading, which means that micro-cracks and micro-voids
ularization technique in dynamics, to a certain extent applicable
also in statics. Among alternative approaches one could cite the are assumed to close completely. The motivation for coupling the
model with plasticity is thus to incorporate irreversible strains.
Hofmann plasticity model [33], which also allows for the descrip-
tion of different properties of concrete in tension and compression The plastic component of the coupled model is formulated in
the effective stress space. The yield condition is as follows
and employs one yield surface within viscoplastic formulation.
In the presented analyses reinforcement is modelled as embed- Fp ¼ r ^ Þ  ry ðjÞ ¼ 0
~ ðr ð6Þ
ded bars of elastic–plastic material with small isotropic hardening,
which means displacement continuity between concrete and steel where r ^; p
~ ðq ^Þ is an equivalent measure of effective stress assumed
pffiffiffiffiffiffiffi
is assumed (bond-slip is neglected for deformed bars). In fact, the in the Burzyński–Drucker–Prager (BDP) form, q ¼ 3J 2 is the clas-
interaction between concrete and reinforcement after cracking is sical deviatoric stress measure, p ¼  13 I1 is the hydrostatic pressure,
incorporated in a simplified manner via the tension stiffening/soft- ry is the yield strength and j is the plastic strain measure (isotropic
ening description of concrete. hardening is assumed).
The plastic multiplier determines the magnitude of plastic
3.1. Gradient-enhanced damage–plasticity (FEAP) strains according to the classical flow rule, which may be non-asso-
ciated, i.e. derived from a potential different than F p . Assuming the
We start the short discussion of the model by selecting a form of standard additive decomposition of strain rates into an elastic and
coupling between plasticity and damage theories. The plasticity a plastic part
theory formulated in stress space is coupled with the scalar dam-
age theory formulated in strain space (intrinsic anisotropy is
_ ¼ _ e þ _ p ð7Þ
neglected). and invoking the local plastic consistency condition F_ p ¼ 0 in order
Considering the damage evolution we distinguish the actual to compute k_ leads to the classical tangential relation:
body with strains  and stresses r and its fictitious undamaged
counterpart with stresses r ^ and strains  ^ . The fictitious counter- r^_ ¼ Dep _ ð8Þ
part represents the undamaged ’’skeleton’’ of the body, and the
In fact, the two theories are coupled by relation (5), in which the
stresses r
^ acting on it are called effective. We adopt the postulate
elastic strain e is computed from the plastic part of the model. The
that the strains observed in the actual body and in its undamaged
stress rate during the evolution of damage and plasticity can be
representation are equal  ¼  ^ [42,43] and that the stresses are
computed by differentiating Eq. (5):
related by means of a scalar damage measure x:
r ¼ ð1  xÞr^ ð1Þ r_ ¼ ð1  xÞ r^_  x_ r^ ð9Þ

The damage x is a function of a damage history parameter jd where r^_ follows from Eq. (8) and the rate of damage during its
and grows from zero to one as jd grows from damage threshold ) is computed as:
evolution (jd ¼ ~
jd0 to some ultimate value. @ x @ ~
Next, we define a damage function which limits the elastic (or x_ ¼ _ ð10Þ
@ jd @ 
elasto–plastic) behaviour of the material in the strain space:
During unloading we have x _ ¼ 0. This leads to the linearized
F d ¼ ~  jd ¼ 0 ð2Þ constitutive relation for the coupled model:
 
where ~  is an equivalent strain measure. During the damage evolu- dx @ ~
  r_ ¼ ð1  xÞ Dep  r
^ _ ð11Þ
tion the history parameter jd ¼ max jd0 ;  ~ . djd @ 
The equivalent strain measure  ~ is defined in a form suitable for
Now, the theory is made gradient-dependent [38] in order to
a quasi-brittle material, taking into account the difference between ensure that numerical simulations of strain localization give mean-
its tensile and compressive strength, cf. [44,45]. Here, the so-called ingful results. In the gradient-enhanced formulation the plasticity
modified von Mises definition of the equivalent strain measure theory remains standard and the damage theory is made nonlocal.
[45] is employed: It is emphasized that in this combination the unstable material
 
~ J 2 ; I1 ; m; r ct ¼ 0 ð3Þ behaviour is caused by damage (the plasticity model in the effec-
tive stress space is hardening).
where J2 and I1 are standard invariants of the (elastic) strain tensor, Following [19], the damage evolution is governed by the follow-
m is Poisson’s ratio, rct is a ratio of uniaxial compressive strength f c ing damage loading function:
and uniaxial tensile strength f t . Alternatively, Mazars definition
[44,46] can be used. F d ¼   jd ¼ 0 ð12Þ
Although damage and plastic processes can be coupled, plastic where the nonlocal strain measure  satisfies the averaging
flow occurs in the undamaged ’’skeleton’’ of the body, so we can equation:
write the elastic constitutive relation between the effective
stresses and elastic strains as follows:   cr2  ¼ ~ ð13Þ
76 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85

with the homogeneous natural boundary condition ðr ÞT n ¼ 0. The F p ðq ^; r

^; p ^ max ; jt ; jc Þ ¼ 0 ð16Þ
parameter c > 0 is related to an internal length scale l and has a unit
where r ^ max is the maximum eigenvalue of the effective stress ten-
of length squared (in one-dimensional situation this relation is
2 sor. Different one-dimensional relations are adopted for the evolv-
c ¼ 12 l , [47]). The internal length is an additional material parameter
ing material strengths in tension and compression, respectively:
which determines the range of non-local interactions of microstruc-
tural deformation carriers. It defines the width of the localization rt ¼ rt ðjt ; j_ t Þ; rc ¼ rc ðjc ; j_c Þ ð17Þ
band. In granular materials it is associated with the grain size, for
The former function is obviously softening and the latter first
concrete it is often related to the maximum aggregate size.
hardening and then softening. Non-associated flow rule is used
It is assumed that the equivalent strain measure is a function of
with a hyperbolic type of BDP potential involving 3 material
elastic strains only, ~ ðe Þ, which means that the damage growth is parameters: dilatancy angle w, so-called eccentricity (deviation
driven by the effective stresses, see Eq. (4), and that the two theo-
from a straight line in p–q space) and initial uniaxial tensile
ries are weakly coupled. If a full coupling exists, as for instance for
strength f t .
plastically induced damage in metals, plastic strains should also
Moreover, two damage variables are defined:
contribute to  ~.
Moreover, the damage growth function xðjd Þ is exponential xt ¼ xt ðjt Þ; xc ¼ xc ðjc Þ ð18Þ
[44] and represents a tensile uniaxial stress–strain relationship
based on experimental observations. It involves the damage to allow for different stiffness degradation in uniaxial tension and
threshold and two more model parameters, i.e. ductility parameter compression. The scalar degradation variable x is derived as a func-
g and numerical parameter a imposing the stress lower bound in tion of xt ; xc and the effective stress.
order to prevent a complete loss of material stiffness which would The model can be equipped with viscoplastic regularization
lead to singularity and premature divergence. according to a generalization of the Devaut–Lions approach, in
The finite element implementation of the model is based on which a viscous upgrade of the plastic strain tensor is derived
the weak-form equilibrium equation and a weak-form of Eq. using an additional viscosity parameter, called the relaxation time
(13). In the ensuing two-field formulation the average strain mea- l, c.f. [17]. The viscoplastic strain rate _ v p which substitutes the
plastic strain rate _ p is computed based on the inviscid state as
sure must be discretized in addition to the displacements. C 0 -con-
follows
tinuity suffices for all shape functions. The coupling to plasticity
influences only the equilibrium equation, while the weak form 1
of the averaging equation is exactly the same as for pure gradient _ v p ¼ ðp  v p Þ ð19Þ
l
damage. Further details of the numerical implementation, consis-
tent linearization and other algorithmic aspects can be found in where p is the plastic strain calculated for the viscosity-indepen-
[18,38]. It should be noted that gradient damage can be useful dent ’’backbone’’ yield surface. The evolution of internal hardening
not only to analyze the macroscale behaviour of concrete. This variables can be treated in a similar way, cf. [52]. The rate depen-
type of modelling can be applied to include effects connected dence related to viscosity provides an additional ductility in the
with the lower levels of observation, e.g. meso- and microscale. model and is an efficient method to overcome the problems of con-
This could include for instance glass-fiber-reinforced composites vergence of cracking simulation algorithm. The inviscid case is
[48], metals [49], biomaterials [50], and damage progress in nano- retrieved when l ¼ 0.
composites [51]. In summary, the parameters which need to be specified for this
In summary, the parameters which need to be specified for this model are: elastic ones, parameters of yield function, dilatancy
model are: the elasticity constants, equivalent strain measure, the angle to define plastic potential, parameters of compression hard-
damage growth function (including the damage threshold), the ening and tension stiffening, specification of tension and compres-
yield function and hardening function, and the internal length scale. sion damage evolution, and finally the viscosity parameter.

3.2. Damaged plasticity model (ABAQUS) 3.3. Short comparison of employed models

The damaged plasticity model in the commercial FE program It is stressed that although the two employed models use
ABAQUS incorporates the influence of moderate confining pressure similar concepts of coupled plasticity and damage, cf. e.g.
and irreversible damage, focusing on failure mechanisms charac- [38,40,43,53], they approach concrete behaviour from different
teristic for quasi-brittle materials, in particular concrete. Certainly, viewpoints. As a result, they represent concrete cracking and fail-
different response is predicted in tension and compression with ure in a different way.
respect to yield strengths (in tension understood as elasticity The former model, implemented in FEAP, is focused on model-
limit), hardening/softening and stiffness degradation. Stiffness ling of regularized material damage, while the plastic counterpart
evolution during cyclic loading and rate sensitivity are included is an upgrade allowing for a better representation of the elastic
in the model. It is based on Eqs. (5) and (7), but the isotropic dam- stiffness degradation which otherwise would be exaggerated. The
age parameter x is here the function of the effective stress and two latter model, available in ABAQUS, is based on the concept of con-
hardening parameters for tension and compression, c.f. [40,41]: crete plasticity and hence yield and plastic potential functions are
  used to represent material failure, while damage variables are
jt merely an upgrade of the model, relevant especially for cyclic load-
x ¼ xðr^ ; jÞ; j ¼ ð14Þ
jc ing to represent properly the elastic stiffness degradation upon
unloading. Therefore, similar model parameters can have different
which are derived from the rate hardening law in the form: values in the two approaches (an example is the dilatancy angle).
_ r Another difference between the two models is the manner in
j_ ¼ khð ^ ; jÞ ð15Þ
which regularization is introduced to prevent the loss of well-
The yield function in the effective stress space is now an posedness of the mathematical model caused by strain softening
upgraded version of the Burzyński–Drucker–Prager (BDP) surface and to avoid pathological discretization sensitivity of numerical
[17]: simulation results.
 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85 77

The gradient-enhanced damage model, based on averaged Table 1

equivalent strain measure according to Eq. (13), belongs to the Numerical simulation data.

strongly non-local models [54], appropriate for both statics and Elastic constants Young’s modulus (GPa) Poisson’s ratio
dynamics. As explained in Section 3.1, an internal length scale Concrete 34.4 0.2
must be specified in addition to the standard material model Steel 205 0.3
parameters. It is noted that the determination of the internal
length parameter for a material is an open issue. It can sometimes Concrete strength (MPa)
be computed from quantities directly measured at the microscopic Tensile, f t 2.13 Compressive, f c 20  f t
level, see e.g. [25]. Otherwise, it must be determined from experi-
Reinforcement Spacing Cover (to axis) Section area Yield strength
ments via an inverse analysis or derived theoretically from
(mm) (mm) (mm2) (MPa)
micromechanics.
In slab:
The extension of the ABAQUS damaged plasticity model for con-
Tensile 100 24 100 455
crete towards viscosity-dependence is a natural upgrade in the Compressive 200 24 100 455
context of dynamic simulations. It assures then that the equations
In column:
of motion remain hyperbolic in the presence of softening, cf. Bars 25 300 455
[29,55]. However, viscosity also provides regularization in the Ties 50 455
static limit as shown for instance by [52,56].

4. Numerical test results

Modified von Mises definition, according to Eq. (3), is employed
and it is treated as the basic case and the ratio r ct ¼ f c =f t is equal
4.1. Simplifications and data
to 20 (this value of r ct approximates the material data). In the sim-
ulation this model can be combined with the Burzyński–Drucker–
A symmetric slab-column geometry and loading are selected for
Prager hardening plasticity. The parameters for the plastic part of
the simulation, hence one quarter of the configuration is analyzed.
the model are: yield strength ry ¼ 2:13 MPa, hardening modulus
The slab and column are discretized with brick elements and linear
for cohesion h ¼ 17:2 GPa, friction and dilatancy coefficients
interpolation is used.
sin u ¼ sin w ¼ 0:5. The behaviour of concrete elements in the col-
A top view of the coarse mesh 1 is presented in Fig. 2(a). Alter-
native meshes are also used in the simulation, in particular with a umn is constrained to be elastic by introducing jo ¼ 6:2  103 . The
double number of elements along the height of the slab (mesh 2). proper value of the internal length scale is examined and discussed
For damaged plasticity model from ABAQUS computations are per- in the next subsection. Note that the internal length influences the
formed with a further densification of elements in each direction damage process according to Eqs. (3) and (13).
(mesh 3). We intend to show that regularization is necessary to The basic parameters for damaged plasticity are as follows:
simulate punching shear behaviour in slab-column connection in dilation angle w ¼ 5:0 (in degrees); strength in uniaxial tension
a proper way. This is the reason why we use different mesh f t ¼ 2:13 MPa, fracture energy Gf ;t ¼ 106:5 N/m or corresponding
densities. stress–strain relationship with residual strength rres t equal to
Perfect bond between concrete and reinforcement is assumed. 0:1 f t or 0:2 f t for the tensile strain which is equal to or greater than
In the numerical simulations the configuration is limited by t1 ¼ 0:02; yield strength for compression f c ¼ 20  f t MPa,
support lines, i.e. a slab segment with dimensions 1500 compressive yield stress ry;c1 ¼ 50:0 MPa for crushing strain
1500  120 mm is computed. Clamps are simulated by supporting c1 ¼ 0:0015. In tension regime a linear softening-constant
nodes near the corner of the slab. The distance between the tensile stress–strain relationship is adopted to represent the residual
reinforcement truss elements is equal to 100 mm in both direc- strength of cracking concrete and to avoid numerical instabilities
tions. The basic data for the computations are specified in Table 1. resulting from zero strength. In particular stress rt decreases line-
The location of the reinforcement is depicted in Fig. 2(b). Static arly from f t to 0:1 f t or 0:2 f t for strain t1 . At strains larger than t1
loading shown in Fig. 2(a) is applied under displacement control, the residual stress remains constant.
i.e. the vertical movement of the upper plane of the column seg- The influence of viscosity parameter l is verified and sensitivity
ment is imposed monotonically. analysis is shown in Section 4.3. When stiffness degradation is
The basic parameters for the gradient damage model are as fol- additionally activated (computations only for mesh 3) then dam-
lows: threshold jo ¼ 6:2  105 , exponential softening with resid- age x is assumed to grow linearly from 0 to 0.9. In uniaxial tension
ual stress parameter a ¼ 0:94 and ductility parameter g ¼ 400. regime this value is achieved for strain t1 ¼ 0:02, while in

(a) Finite elements – mesh 1. (b) Reinforcement in simulated specimen.

Fig. 2. Simulated specimen.

78 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85

compression linear damage growth is assumed until c2 ¼ 0:2. In the whole equilibrium path is not representative. It should be
both cases, when damage reaches 0.9 it is subsequently held con- noted that larger values of internal length result in a larger damage
stant to preserve a residual stiffness and strength. The other model zone. The smoothing effect for l ¼ 20 mm is quite strong and an
parameters are defaults suggested in ABAQUS documentation [17]. unrealistic crack pattern is predicted in a large part of the slab.
The equilibrium paths for the two considered meshes are also plot-
ted for the cases with smaller l (l ¼ 4 mm or 8 mm) in Fig. 3. For the
4.2. Results for gradient-enhanced damage–plasticity case with l ¼ 4 mm and mesh 2 the computations diverge when
the load reaches about 250 kN, similarly to the experimental load
The presentation of the results obtained with the gradient dam- carrying capacity. For the other cases computations diverge earlier.
age model begins with the diagrams of applied load denoted by P If l ¼ 8 mm is adopted, in the final stage a change of stiffness is
versus the maximum deflection of the slab denoted by w (in fact noticed into small softening.
it is the vertical displacement used for load control, but since the The deformation of the slab at the final stage is depicted in Fig. 4
column is very stiff, it is approximately equal to the slab deflection for l ¼ 4 mm and the two adopted meshes. It is visible that for both
at the center). In Fig. 3 the diagrams of the load–displacement meshes the slab is not cut directly at the column and the deformed
curves for gradient damage and three values of internal length elements do not form a clear shear crack as is seen in experiments.
are shown together with the experimental equilibrium path. In Figs. 5 and 6 illustrate the simulated averaged strain patterns
the initial elastic phase the numerical diagram exhibits a too large for the case l ¼ 4 mm and both meshes. The contour plots are made
slope. Most probably the original initial stiffness of the specimen is for the bottom face of the slab and for the cross section along the
affected by precracking due to shrinkage. Moreover, in the numer- vertical symmetry plane (called shortly ’’side view’’). Side views
ical simulation it is implicitly assumed that the support system of of averaged strain distributions show that the simulated slab is
the experimental setup is infinitely stiff while in reality it is not. cracked near the connection, although the shear cone is repro-
After cracking, for deflection w > 1 mm, the slope of the dia- duced in a smeared way. However, for the small value of the inter-
grams is similar to the one observed in the experiment. The peak nal length parameter, many flexural cracks are initiated at column
for pure gradient damage and l ¼ 20 mm occurs too early and corners and propagate radially. We can conclude that the averaged
strain distributions poorly simulate circular cracks at some dis-
tance from the column, hence the punching shear fracture along
inclined cone-shaped surface is rather not clearly reproduced.
The best case is for l ¼ 4 mm and mesh 2, where contour plots
can be interpreted as flexural and punching failure. Fig. 7 presents
the bottom and side views for l ¼ 8 mm. It was expected that the
slightly more regularized model would make the results more cor-
rect and that circular cracks would be reproduced properly
together with radial ones. It turns out that introducing a little
stronger non-locality, which corresponds to about two elements
in the mesh, does not improve the results. Of course the localiza-
tion zone is wider, but l ¼ 8 mm results in dominating averaged
strains closer to the column than it is for l ¼ 4 mm.
In Fig. 8 the diagrams for 3 versions of the gradient damage
Fig. 3. Load–displacement diagrams, influence of internal length scale, gradient model and mesh 2 are presented. The equilibrium path obtained
damage (without plasticity). for damage–plasticity (modified von Mises and coupling with

Fig. 4. Deformation (magnification factor is 10) of slab for two different meshes, gradient damage, l ¼ 4 mm.

(a) Bottom view. (b) Side view.

Fig. 5. Contour plots of averaged strain for mesh 1, gradient damage, l ¼ 4 mm, deflection w  5:5 mm.
 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85 79

(a) Bottom view. (b) Side view.

Fig. 6. Contour plots of averaged strain for mesh 2, gradient damage, l ¼ 4 mm, deflection w  9:5 mm.

(a) Bottom view. (b) Side view.

Fig. 7. Contour plots of averaged strain for mesh 2, gradient damage, l ¼ 8 mm, deflection w  6:75 mm.

plasticity – dashed line) is stiffer than for pure damage (modified

von Mises – bold solid line), since damage is driven by the elastic 250
part of strain and grows slower, but the solution diverges when
the load approaches 200 kN. The results for the alternative loading 200
ForceP [kN]

function according to Mazars’ definition (dotted line), which seems Experiment

150
more accurate for concrete in tension–compression regime, com-
pletely fail, cf. Fig. 8.
100
1
4.3. Results for damaged plasticity model from ABAQUS 50

We start the presentation of the results for this model with dia- 0
0 2 4 6 8 10 12
grams obtained for the coarse mesh 1 compared with the experi-
mental response. In Fig. 9 four options are considered. The Deflection w [mm]
diagrams denoted as Gf ;t and 2Gf ;t use the parameters described Fig. 9. Load–displacement diagrams, influence of tensile softening parameters for
in Section 4.1 (with either Gf ;t or 2Gf ;t as fracture energy value the model without viscosity, mesh 1.
for tension). For these analyses the divergence occurs for deflection

250
250

200
200
ForceP [kN]
ForceP [kN]

Experiment Experiment
150 150
modified von Mises 3
100 and coupling with plasticity 100
modified von Mises 2 2G f t
50 50
Gf t
Mazars definition
0 0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Deflection w [mm] Deflection w [mm]

Fig. 8. Load–displacement diagrams, mesh 2, l ¼ 4 mm, different versions of Fig. 10. Load–displacement diagrams, influence for viscous enhancement,
gradient damage. l ¼ 0:001, mesh 1.
80 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85

of w  1:6 mm. It should be observed that the numerical response diagrams in Fig. 9 are for the cases when tensile response after
does not change significantly when the additional ductility is intro- cracking is modelled by stress–strain diagrams with the residual
duced via artificially increased fracture energy 2Gf ;t . The other two strength rres
t ¼ 0:1 f t or 0:2 f t . In these two cases, the numerical

250
0 01

200 = 0.005

ForceP [kN]
150
5
100

Experiment 0 002
50
0 001
0
0 2 4 6 8 10 12
Deflection w [mm]

Fig. 11. Load–displacement diagrams, influence of viscosity parameter l for Fig. 12. Load–displacement diagrams, influence of viscosity parameter l for
residual tensile strength rres
t ¼ 0:2f t , mesh 1. residual tensile strength rres
t ¼ 0:2f t , mesh 3.

(a) Deformation. (b) Side view.

(c) Bottom view. (d) Diagonal section.

Fig. 13. Equivalent tensile plastic strain and deformation for mesh 1, deflection w  1:5 mm, for point 1 in Fig. 9, case with single value of Gf ;t .

(a) Deformation. (b) Side view.

(c) Bottom view. (d) Diagonal section.

Fig. 14. Equivalent tensile plastic strain and deformation for mesh 1, w  7:5 mm, for point 2 in Fig. 10, case with single value of Gf ;t and l ¼ 0:001.
 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85 81

response is more realistic with the maximum deflection w reaching Fig. 12 presents a similar comparison for mesh 3. The response
almost 6 mm and the maximum load P  180 kN. This type of of the numerical model is slightly more ductile than for the coarser
modelling approach of postcracking behaviour is a starting point meshes, the predicted load carrying capacity grows with the
in the following calculations. increase of l, but the results for l ¼ 0:01 now underestimate it.
The diagrams in Fig. 10 show the influence of the viscous The next series of figures presents the insight into the behaviour
enhancement on the previous computations. For l ¼ 0:001 a con- of the numerical model for selected states marked in load–deflec-
siderable increase of ductility is observed, involving moderate soft- tion diagrams (cf. Figs. 9–12). The results for the coarse mesh 1 are
ening. The diagrams obtained for rres t ¼ 0:1 f t or 0:2 f t approximate shown first. The distribution of the equivalent tensile plastic strain
much better the experimental response. (PEEQT) is monitored in the vertical symmetry plane (called side
In Fig. 11 the sensitivity of the simulation to the value of time view), in the vertical diagonal cross-section and at the bottom sur-
relaxation parameter (viscosity) is verified. The best prediction of face of the slab. Moreover, the displacement mode is shown.
load-carrying capacity is obtained for l ¼ 0:01, while the lack of According to [17] cracking is monitored at sampling points where
viscous regularization results in premature divergence. On the the equivalent tensile plastic strain is larger than zero.
other hand, for l ¼ 1:0 the solution is artificially stiff and unac- The results for the basic set of material parameters are shown in
ceptable. Analogical behaviour as for gradient damage is also Fig. 13 for deflection w  1:5 mm, just before the analysis diverges,
observed – the lack of regularization does not permit one to i.e. for state marked 1 in Fig. 9. The flexural fracture mode is visible
simulate punching shear while too strong non-locality leads to in the slab next to the column with some damage zones propagat-
completely diffuse crack pattern. ing at the bottom of the slab. The distributions of PEEQT for the

(a) Deformation. (b) Side view.

(c) Bottom view. (d) Diagonal section.

Fig. 15. Equivalent tensile plastic strain and deformation for mesh 1, w  12:5 mm, for point 3 in Fig. 10, case with residual strength rres
t ¼ 0:1f t and l ¼ 0:001.

(a) Deformation. (b) Side view.

(c) Bottom view. (d) Diagonal section.

Fig. 16. Equivalent tensile plastic strain and deformation for mesh 1, w  12:5 mm, for point 4 in Fig. 11, case with residual strength rres
t ¼ 0:2f t and l ¼ 0:01.
82 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85

ductility of the model increased, i.e. for other options shown in increased to l ¼ 0:01. The experimental load-carrying capacity is
Fig. 9, monitored at the final analysis step, are quite similar. then approximated very well. However, as can be seen in Fig. 16
In Fig. 14 the monitored fields are shown for the regularized (point 4 in Fig. 11), there is no significant change in the distribution
model with l ¼ 0:001. The analysis can then be continued – the of equivalent tensile plastic strain. This an example of situation in
plots are made for w  7:5 mm (point 2 in Fig. 10). The deforma- which an acceptable simulation of the experimental load–displace-
tion mode and the PEEQT distribution exhibit a distributed dam- ment diagram does not guarantee proper results in terms of
age zone, triangular in cross-section, with the largest values in predicted failure mode.
the middle of the slab, at the column. However, no clear punching From the results presented so far we have concluded that too
cone is formed and the experimental load-carrying capacity is not large values of the relaxation parameter should be avoided. We
reproduced. With the additional ductility, where rres t ¼ 0:1f t , the also suspected that the adopted discretization is too coarse, espe-
limit load is predicted much better, together with a softening cially along the depth of the slab. Using mesh 2 with 10 elements
stage, but for w  12:5 mm (point 3) the inelastic deformation along the depth and intermediate value of l ¼ 0:002 we have
mode does not seem to represent punching (see Fig. 15), it is obtained the results in Fig. 17. The failure mode then exhibits an
rather localized at the bottom surface of the slab next to the inclined fracture zone. However, the flexural crack still seems
column. dominant, the inclination of the crack band at the column is too
The plots of PEEQT for the larger residual strength rres
t ¼ 0:2f t small and the band seems to follow the main reinforcement plane
and w ¼ 12:5 mm are almost the same. As shown in Fig. 11 the as it propagates into the slab, suggesting a mesh bias caused by the
softening final stage is delayed when the viscosity parameter is proportions of element dimension.

(a) Deformation. (b) Side view.

(c) Bottom view. (d) Diagonal section.

Fig. 17. Equivalent tensile plastic strain and deformation for mesh 2, w  12:5 mm, residual strength rres
t ¼ 0:2f t and l ¼ 0:002.

(a) Deformation. (b) Side view.

(c) Bottom view. (d) Diagonal section.

Fig. 18. Equivalent tensile plastic strain and deformation for mesh 3, w  12:5 mm, for point 5 in Fig. 12, case with residual strength rres
t ¼ 0:2f t and l ¼ 0:002.
 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85 83

250

200

Force P [kN]
150
without damage
100

50 Experiment with damage

0
0 2 4 6 8 10 12
Deflection w [mm]
(b) Deformation.
(a) Load-displacement diagrams.

(c) Side view, PEEQT. (d) Side view, damage t.

(e) Bottom view, PEEQT. (f) Bottom view, damage t.

(g) Diagonal section, PEEQT. (h) Diagonal section, damage t.

Fig. 19. Results for mesh 3, case with residual strength rres
t ¼ 0:2f t , viscosity parameter l ¼ 0:002 and active damage component. Comparison of diagrams for deflection
w  12:5 mm: deformation, distributions of equivalent tensile plastic strain (PEEQT), damage for tension xt .

With these results still not satisfactory, we densify the mesh illustrated in Fig. 18. We can observe that the contour plots of
also in plane and the load–displacement curves for this series of the equivalent tensile plastic strain (PEEQT) and damage for ten-
computations are shown in Fig. 12. The equivalent tensile plastic sion xt look almost the same. The multiple diagonal cracks visible
strain distribution in Fig. 18 finally seem to represent a realistic for ’’side views’’ seem to remind punching shear failure, even
fracture mode, although the respective load–displacement diagram though the flexural crack is not excluded.
underestimate the experimental load-carrying capacity (softening
occurs about w ¼ 10 mm). To improve this, we can increase the 5. Final remarks and future work
viscosity parameter, although there is a risk involved of distribut-
ing the process into one triangular fracture zone. The crack at the In the paper two regularized numerical models for concrete,
’’side view’’ (Fig. 18(b)) looks similar to the results obtained for formulated within elastic–plastic-damage theories, have been
gradient damage (Figs. 5(b) and 6(b)). employed in the three-dimensional FE analysis of punching shear
In Fig. 19 the results for mesh 3 are again presented, but for this in the slab-column configuration formerly examined in the labora-
case the damage (stiffness degradation) is additionally activated. tory [11]. The interest has been limited to the simulation of static
The load–displacement diagrams are compared in Fig. 19(a). response of monotonically increasing imposed displacement of the
Although the equilibrium paths are similar, with damage influence column. In particular, the gradient damage model implemented by
the deformations and the internal variable distributions show dif- the authors in FEAP and the rate-dependent damaged plasticity
ferent crack patterns in comparison to the case with no damage model from ABAQUS have been used.
84 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85

The following conclusions can be drawn: [7] Kinnunen S, Nylander H. Punching of concrete slabs without shear
reinforcement. Tech rep No. 158. Stockholm, Sweden: Transactions of the
Royal Institute of Technology; 1960.
 Properly calibrated damage–plasticity models can be used for [8] Dilger WH, Ghali A. Shear reinforcement of concrete slabs. ASCE J Struct Div
predicting shear behaviour and failure in RC slabs. In the cali- 1981;107(ST12):2403–20.
bration it is important to consider load–deflection diagrams, [9] Yamada T, Nanni A, Endo K. Punching shear resistance of flat slabs: Influence of
reinforcement type and ratio. ACI Struct J 1992;89(4):555–63.
displacement patterns and cracking representation. [10] El-Salakawy EF, Polak MA, Soliman MH. Reinforced concrete slab-column edge
 As expected, the FE analysis is very sensitive to the adopted rep- connections with openings. ACI Struct J 1999;96(1):79–87.
resentation of the tensile concrete behaviour. Punching failure [11] Adetifa B, Polak MA. Retrofit of interior slab-column connections for punching
using shear bolts. ACI Struct J 2005;102(2):268–74.
is initiated by tensile cracking, however aggregate interlock, [12] Polak MA. Modelling punching shear of reinforce concrete slabs using layered
shear friction (and reinforcement dowel action) allow the con- finite elements. ACI Struct J 1998;95(1):71–80.
figuration to carry a substantial amount of load after the initial [13] Menétrey P, Willam KJ. Punching shear in reinforced concrete: localized
process. In: deBorst R et al., editors. Proc EURO-C 1998 int conf computational
cracking. This phenomenon must be reflected by the FEM modelling of concrete structures, vol. 2. Rotterdam/Brookfield: A.A. Balkema;
model, otherwise, as shown in Fig. 9, the numerical predictions 1998. p. 867–76.
suggest failure immediately after first cracking. [14] Xiao RY, O’Flaherty T. Finite-element analysis of tested concrete connections.
Comput Struct 2000;78:247–55.
 Even with proper modelling of tension, premature failure is pre- [15] Lewiński PM, Wie˛ch PP. FEM analysis of punching failure of RC slabs versus
dicted due to localized deformation related to flexural and shear results of experimental research. Theor Found Civ Eng 2007;XV:433–40 [in
cracking. A method to avoid such numerical response is to use Polish].
[16] Negele A, Eligehausen R, Ožbolt J, Polak MA. Finite-element simulations on
regularization. The two types of localization limiters used in
punching tests of shear-retrofitted slab-column connections. In: Carpinteri A,
this work both give similar results, they spread damage (crack- Gambarova P, editors. Fracture mechanics of and concrete structures.
ing, inelasticity) effect among surrounding elements. FRAMCOS 6. Catania/Leiden: A.A. Balkema; 2007. p. 911–8.
 Without regularization, only the initial stage of cracking is sim- [17] SIMULIA Abaqus, Abaqus theory manual (6.10). Tech rep, Dassault Systemes,
Providence, RI, USA; 2010.
ulated, then local mesh distortions in damaged zones and at [18] Wosatko A. Finite-element analysis of cracking in concrete using gradient
concrete–reinforcement interface affect negatively the conver- damage–plasticity, Ph.D. dissertation. Cracow: Cracow University of
gence. Possibly element removal algorithm should be used for Technology; 2008.
[19] Peerlings RHJ, de Borst R, Brekelmans WAM, de Vree JHP. Gradient-enhanced
fully damaged elements, but this is impossible in the codes used damage for quasi-brittle materials. Int J Numer Meth Eng 1996;39:3391–403.
here for simulations. Moreover, a more realistic representation [20] Taylor RL. FEAP – a finite element analysis program, version 7.4, user manual.
of bond-slip could reduce the problems, but involves significant Tech rep. Berkeley: University of California at Berkeley; 2001.
[21] Belytschko T, Lasry D. A study of localization limiters for strain-softening in
complications of the numerical model. statics and dynamics. Comput Struct 1989;33:707–15.
 With regularization damage zones representing cracks are not [22] Aifantis EC. On the microstructural origin of certain inelastic models. J Eng
localized, which can make it difficult for the punching cone to Mater Technol 1984;106:326–30.
[23] Mühlhaus H-B, Aifantis EC. The influence of microstructure-induced gradients
form after the initial stage of flexural fracture. This means the on the localization of deformation in viscoplastic materials. Acta Mech
parameters governing the regularization (internal length l or 1991;89:217–31.
viscosity parameter l) should be as small as possible and as [24] de Borst R, Sluys LJ, Mühlhaus H-B, Pamin J. Fundamental issues in finite
element analyses of localization of deformation. Eng Comput 1993;10:99–121.
large as required to obtain proper results (stable and consistent
[25] Geers MGD. Experimental analysis and computational modelling of damage
with experimental reality). and fracture, Ph.D. dissertation. Eindhoven: Eindhoven University of
Technology; 1997.
In future research the numerical simulation will be extended to [26] Aifantis EC. Update on a class of gradient theories. Mech Mater 2003;35:
259–80.
includes the behaviour of slabs strengthened against punching [27] Askes H, Aifantis EA. Gradient elasticity in statics and dynamics: an overview
using shear reinforcement. Also, the issue of the influence of unbal- of formulations, length scale identification procedures, finite element
anced moments and seismic forces on punching shear phenome- implementations and new results. Int J Solids Struct 2011;48:1962–90.
[28] Needleman A. Material rate dependence and mesh sensitivity in localization
non should be studied. This work will have practical applications problems. Comput Methods Appl Mech Eng 1988;67:69–86.
in code calibrations and verifications. [29] Sluys LJ. Wave propagation, localization and dispersion in softening solids,
Ph.D. dissertation. Delft: Delft University of Technology; 1992.
[30] Dubé J-F, Pijaudier-Cabot G, Borderie CL. Rate dependent damage model for
Acknowledgements concrete in dynamics. J Eng Mech 1996;122(10):939–47.
[31] Wang WM, Sluys LJ, de Borst R. Viscoplasticity for instabilities due to strain
softening and strain-rate softening. Int J Numer Meth Eng
The research was supported by the European Union through the 1997;40(20):3839–64.
European Social Fund within project Cracow University of Technol- [32] Glema A, Łodygowski T, Perzyna P. Interaction of deformation waves and
localization phenomena in inelastic solids. Comput Methods Appl Mech Eng
ogy development program – top quality teaching for the prospective 2000;183(1–2):123–40.
Polish engineers University of the 21st century (Contract No. [33] Winnicki A, Pearce CJ, Bićanić N. Viscoplastic Hoffman consistency model for
UDA-POKL.04.01.01-00-029/10-00) and Natural Sciences and concrete. Comput Struct 2001;79:7–19.
[34] Ngo D, Scordelis AC. Finite element analysis of reinforced concrete beams. ACI J
Engineering Research Council (NSERC) at Canada. The support is 1967;64:152–63.
gratefully acknowledged. [35] Rots JG. Computational modeling of concrete fracture, Ph.D. dissertation. Delft:
Delft University of Technology; 1988.
[36] Feenstra PH. Computational aspects of biaxial stress in plain and reinforced
References concrete, Ph.D. dissertation. Delft: Delft University of Technology; 1993.
[37] Mang HA, Hofstetter G. Computational mechanics of reinforced concrete
structures. Braunschweig, Wiesbaden: Vieweg Verlag; 1995.
[1] Menétrey P. Synthesis of punching failure in reinforced concrete. Cem Concr
[38] de Borst R, Pamin J, Geers MGD. On coupled gradient-dependent plasticity and
Compos 2002;24:497–507.
damage theories with a view to localization analysis. Eur J Mech A/Solids
[2] Urban TS. Punching in concrete. Selected problems, monograph 959. Łódź:
1999;18(6):939–62.
Łódź Technical University; 2005. [in Polish].
[39] Pamin J, de Borst R. Stiffness degradation in gradient-dependent coupled
[3] Polak MA, El-Salakawy E, Hammill N. Shear reinforcement for concrete slabs.
damage–plasticity. Arch Mech 1999;51(3–4):407–33.
ACI SP 232 Punching Shear Reinf Concr Slabs 2005:75–96.
[40] Lubliner J, Oliver J, Oller S, Oñate E. A plastic-damage model for concrete. Int J
[4] Megally S, Ghali A. Cautionary note on shear capitals. In: Concrete
Solids Struct 1989;25(3):299–326.
international; 2002. p. 75–82.
[41] Lee J, Fenves G. Plastic-damage model for cyclic loading of concrete structures.
[5] Moe J. Shearing strength of reinforced concrete slabs and footings under
ASCE J Eng Mech 1998;124(8):892–900.
concentrated loads. Tech rep V. D47. Illinois: Portland Cement Association;
[42] Simo JC, Ju JW. Strain- and stress-based continuum damage models - I.
1961.
Formulation, II. Computational aspects. Int J Solids Struct 1987;23(7):
[6] Elstner RC, Hognestad E. Shearing strength of reinforced concrete slabs. J ACI
821–69.
1956;53(2):29–58.
 A. Wosatko et al. / Computers and Structures 151 (2015) 73–85 85

[43] Ju JW. On energy-based coupled elastoplastic damage theories: constitutive [50] Waffenschmidt T, Polindara C, Menzel A, Sergio B. A gradient-enhanced large-
modeling and computational aspects. Int J Solids Struct 1989;25(7):803–33. deformation continuum damage model for fibre-reinforced materials. Comput
[44] Mazars J, Pijaudier-Cabot G. Continuum damage theory – application to Methods Appl Mech Eng 2014;268:801–42.
concrete. ASCE J Eng Mech 1989;115:345–65. [51] Dimitrijevic BJ, Aifantis KE, Hackl K. The influence of particle size and spacing
[45] de Vree JHP, Brekelmans WAM, van Gils MAJ. Comparison of nonlocal on the fragmentation of nanocomposite anodes for Li batteries. J Power
approaches in continuum damage mechanics. Comput Struct 1995;55(4): Sources 2012;206:343–8.
581–8. [52] Simo JC, Hughes TJR. Computational inelasticity. Interdisciplinary applied
[46] Mazars J. Application de la mécanique de l’edommagement au comportement mathematics, vol. 7. New York: Springer-Verlag; 1998.
non linéaire et à la rupture du béton de structure, Ph.D. dissertation, Université [53] Hansen NR, Schreyer HL. A thermodynamically consistent framework for
Paris 6, Paris; 1984. theories of elastoplasticity coupled with damage. Int J Solids Struct
[47] Askes H, Pamin J, de Borst R. Dispersion analysis and element-free Galerkin 1994;31(3):359–89.
solutions of second- and fourth-order gradient-enhanced damage models. Int J [54] Peerlings RHJ, Geers MGD, de Borst R, Brekelmans WAM. A critical comparison
Numer Meth Eng 2000;49:811–32. of nonlocal and gradient-enhanced softening continua. Int J Solids Struct
[48] Geers MGD, de Borst RD, Brekelmans WAM, Peerlings RHJ. Validation and 2001;38(44–45):7723–46.
internal length scale determination for a gradient damage model: application [55] Łodygowski T. Theoretical and numerical aspects of plastic strain localization.
to short glass-fibre-reinforced polypropylene. Int J Solids Struct Monograph 312. Poznań: Poznań University of Technology; 1996.
1999;36(17):2557–83. [56] Winnicki A, Viscoplastic and internal discontinuity models in analysis of
[49] Peerlings RHJ, Brekelmans WAM, de Borst R, Geers MGD. Gradient-enhanced structural concrete, Series civil engineering. Cracow: Cracow University of
damage modelling of high-cycle fatigue. Int J Numer Meth Eng 2000;49(12): Technology; 2007.
1547–69.