PLASTIC-DAMAGE MODEL FOR CYCLIC LOADING OF CONCRETE STRUCTURES - Lee & Fenves 1998

Prévia do material em texto

PLASTIC-DAMAGE MODEL FOR CYCLIC LOADING
OF CONCRETE STRUCTURES
By Jeeho Lee! and Gregory L. Fenves/ Member, ASCE
ABSTRACT: A new plastic-damage model for concrete subjected to cyclic loading is developed using the
concepts of fracture-energy-based damage and stiffness degradation in continuum damage mechanics. Two dam-
age v~ables, .one for te~sile damage and the other for compressive damage, and a yield function with multiple-
hardemn~ vanables are tntroduced .to account for different damage states. The uniaxial strength functions are
factored mto two parts, correspondmg to the effective stress and the degradation of elastic stiffness. The con-
s~itutive relations .for elastopla~tic ~esponses ar~ decoupled from the degradation damage response, which pro-
vIdes advantages m the numencal Implementation. In the present model, the strength function for the effective
stress is used to control the evolution of the yield surface, so that calibration with experimental results is
convenient. A simple and thermodynamically consistent scalar degradation model is introduced to simulate the
effec~ of damage on el~stic stiffness and ~ts recovery during crack opening and closing. The performance of the
plastIc-damage model IS demonstrated WIth several numerical examples of simulating monotonically and cycli-
cally loaded concrete specimens.
INTRODUCTION
The modeling of crack initiation and propagation is one of
the most important aspects in the failure analysis of concrete
structures. The cracking process in concrete is distinguished
from cracking of other materials, such as metal and glass, in
that it is not a sudden onset of new free surfaces but a con-
tinuous forming and connecting of microcracks (Mehta and
Monteiro 1993). The formation of microcracks is represented
macroscopically as softening behavior of the material, which
causes the localization and redistribution of strain in a struc-
ture.
This phenomenological behavior at the macroscopic level
can be modeled by classical plasticity (Pramono and Willam
1989; Chen 1994). On the other hand, the microcracking pro-
cess also causes stiffness degradation, which can be observed
in concrete subjected to cyclic loading (Karsan and Jirsa 1969;
Gopalaratnam and Shah 1985). Modeling the stiffness degra-
dation is difficult to represent with classical plasticity. In con-
tinuum damage mechanics, the degradation can be modeled by
defining the relationship between stresses and effective
stresses. Comprehensive reviews on continuum damage me-
chanics are given by Kachanov (1986) and Lemaitre and Cha-
boche (1990). Several models for concrete based on continuum
damage mechanics have been developed (Mazars 1986; Ma-
zars and Pijaudier-Cabot 1989; Cervera et al. 1995). However,
without inelastic (or plastic) strains the continuum damage me-
chanics theory cannot provide an appropriate dilatancy control,
which is very important for simulating plain and reinforced
concrete structures under multiaxial loading.
In the plastic-damage approach, stiffness degradation is im-
bedded in a plasticity model (Simo and Ju 1987; Ju 1989;
Lubliner et al. 1989; yazdani and Schreyer 1990). In the cou-
pled elastoplastic-damage model (Simo and Ju 1987; Ju 1989),
the effective-stress concept in continuum damage mechanics
is introduced to represent stiffness degradation. The advantage
of this model is that the damage in stiffness degradation, which
'Postdoctoral Res. Fellow, Dept. of Civ. and Envir. Engrg., Univ. of
California, Berkeley, CA 94720.
'Prof., Dept. of Civ. and Envir. Engrg., Univ. of California, Berkeley,
CA.
Note. Associate Editor: George V. Voyiadjis. Discussion open until
January I, 1999. To extend the closing date one month, a written request
must be filed with the ASCE Manager of Journals. The manuscript for
this paper was submitted for review and possible publication on February
3, 1997. This paper is part of the Journal of Engineering Mechanics,
Vol. 124, No.8, August, 1998. ©ASCE, ISSN 0733-9399/98/0008-
0892-0900/$8.00 + $.50 per page. Paper No. 15055.
892/ JOURNAL OF ENGINEERING MECHANICS / AUGUST 1998
originally is coupled in the constitutive relations, can be de-
coupled from the plastic deformation by linearizing the evo-
lution equations. However, as seen in other models based on
continuum damage mechanics, calibration of the parameters
determining the evolution of a yield surface with experimental
data is difficult because most experimental data are based on
the stress.
In a model (referred to as the Barcelona model in this paper)
proposed by Lubliner and coauthors (Lubliner et al' 1989), a
fracture-energy-based scalar damage variable is used to rep-
resent all damage states. In addition to the damage variable,
the model introduces elastic and plastic degradation variables
to simulate the degradation of elastic stiffness. The degradation
variables are coupled with the plastic deformation in the con-
stitutive relations, making it convenient to calibrate the pa-
rameters with experimental results. However, the coupled re-
lations give a complicated and unstable numerical algorithm,
which causes spurious plastic unloading during iteration (Lee
and Fenves 1994).
Because a quasi-brittle material under cyclic loading under-
goes several damage states, such as tensile cracking, com-
pressive failure, and stiffness degradation, an assumption of a
single damage variable is not adequate. To account for differ-
ent damage responses, multiple hardening (or damage) varia-
bles can be used (Murray et al. 1979; Mazars 1986; Ohtani
and Chen 1988; Mazars and Pijaudier-Cabot 1989). Isotropic
continuum damage mechanics models with multiple damage
variables (Mazars 1986; Mazars and Pijaudier-Cabot 1989)
cannot represent the different effects of damage on the tensile
and compressive strengths because the damage variables even-
tually contribute to the same isotropic evolution for both
strengths.
When cyclic loading ranges between the tensile and the
compressive, the recovery of degraded stiffness is observed
during unloading from the tensile region to the compressive
region (Reinhardt 1984). The stiffness recovery is a conse-
quence of closing the previously opened cracks. Several mod-
els for simulating the stiffness recovery have been suggested
in the context of isotropic and anisotropic damage models (Or-
tiz 1985; Ju 1989; Cervera et al. 1995; Hansen and Schreyer
1995).
In the present study, a new plastic-damage model for cyclic
loading is developed using the concept of fracture-energy-
based damage similar to the Barcelona model. Two damage
variables, one for tensile damage and the other for compressive
damage, account for the different damage states. The yield
function proposed by Lubliner and coauthors (Lubliner et al.
1989) is modified using multiple damage (or hardening) var-
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.
FRAMEWORK FOR PLASTIC-DAMAGE MODEL
(6)
(7)
(8)
(lOa,b)
K =).H(O', K)
K =Hi(r(a), a, K)
f, =f,(K,); fc =fc(KC>
it can be seen directly that D represents the degradation of the
elastic stiffness.
The plastic strain rate is evaluated by a flow rule, which is
assumed to be generated from a scalar plastic potential func-
tion 41. For a plastic potential defined in the effective-stress
space, the plastic strain is given by
As will be shown later, the function H can be derived consid-
ering plastic dissipation.
From experimental tests of quasi-brittle materials, it is ob-
served that there is a surface in the stress space, for which a
failure or damaged state can be determined (Chen 1994). The
surface, called the yield surface, evolves with the damage var-
iable, or the hardening variable in the context of classical plas-
ticity, and limits the current admissible stresses. Because the
responses from tensile damage andcompressive damage on
quasi-brittle materials are quite different, it is not realistic to
represent all damage states by a single parameter. To account
for the different damage responses of concrete in tension and
compression, a multihardening (or multisoftening) yield func-
tion, which has more than one variable to describe the evolv-
ing yield surface, is used. If two state variables f, and!c, rep-
resenting the uniaxial tensile strength and compressive
strength of the material, respectively, describe the yield func-
tion, the admissible stress states are constrained by the follow-
ing condition:
where r =a zero-degree homogeneous function of the stress.
Because the scalar degradation damage variable (1 - D) is used
and r(&) = r(x&) for any positive scalar quantity x, H(r(&),
iT, K) =H(r(a), iT, K) or simply
where)' is not the derivative of the Lame's constant, but a
nonnegative function referred to as the plastic consistency pa-
rameter.
Another internal variable set, in addition to the plastic strain,
is needed to represent the damage states. In the present model,
it is assumed that the damage variable K is the only necessary
state variable and its evolution is expressed as
F(O', f" fc) S 0 (9)
It is assumed that F is a first-degree homogeneous function
with respect to all three variables; in other words xF«(J', f,,!c)
= F(x(J', xf" xfc) for any positive scalar quantity x. Also, it is
assumed that F itself is an isotropic function in the stress
space. The yield function remains isotropic even though its
hardening behavior provides a more general evolution that the
simple expansion or contraction of the yield surface.
In the present model, the uniaxial strength functions in (9)
are expressed in terms of two damage parameters K, and Kc
Assuming that both relations in (10) can be factored into the
degradation damage and effective-stress responses in the same
form as (5), it follows that:
(3)
(2)
(4)
(la,b)
O'=D:(J'
From the concept of the continuum damage theory, the
stress is mapped into the effective stress iT by a rank four
tensor D
Plastic-damage and continuum damage models have been
developed using thermodynamical approaches (Ju 1989; Han-
sen and Schreyer 1994). Constitutive relations are derived
from thermodynamic potential functions and, consequently,
are consistent in the thermodynamic context. In the thermo-
dynamical approach, because well-defined potential functions
are required and the associativity between state variables and
their conjugate forces are preferred, it is difficult to construct
the constitutive system based on physical observations, which
is more appropriate for modeling quasi-brittle materials such
as concrete. For instance, applying the associative plastic flow
rule to a concrete model gives a poor result in dilatancy eval-
uation. In the present study, the constitutive equations for a
plastic-damage model are derived in a more direct, but still
thermodynamically consistent, manner. A thermodynamic in-
terpretation of the present model is given in Appendix I.
In the incremental theory of plasticity, the strain tensor £ is
decomposed into the elastic part E' and the plastic part rt',
which for linear elasticity is given by
iables. The uniaxial strength functions are factored into two
parts corresponding to the effective stress and degradation of
the elastic stiffness. The latter is referred to as degradation
damage in this paper. The constitutive relations for elastoplas-
tic responses are decoupled from the degradation damage re-
sponse. The strength function for effective stress is used to
control the evolution of the yield surface, so that calibration
with experimental results is convenient. Finally, a simple and
thermodynamically consistent stiffness recovery scheme is in-
troduced for simulating crack opening and closing. In the pres-
ent approach, although the damaged elastic stiffness remains
isotropic, the plastic-damage model provides the separate ev-
olution of the tensile and compressive strengths and induced
directional damage by plastic strains.
Several numerical examples of concrete modeling are
shown with the experimental counterparts to demonstrate the
effectiveness of the present plastic-damage approach. The de-
velopment of the model is limited to infinitesimal deformation
theory, which is adequate for modeling concrete materials. It
is noted that i and i l respectively, denotes the eigenvalue ma-
trix and the ith eigenvalue of x.
where the elastic stiffness E is a rank four tensor; and (J' =
stress tensor. In plastic-damage theories, the plastic strain rep-
resents all irreversible deformations including those caused by
microcracks. From (1) the stress-strain relation is
Because the effective stress is defined with the undamaged
elastic stiffness from (2), it becomes
where Eo = initial elastic-stiffness tensor. For practical reasons,
isotropic degradation damage is assumed in many cases. If a
scalar degradation damage variable 0 S D < 1 is used to rep-
resent the isotropic damage, then D = 1I(1 - D)I, where I is
the rank four identity tensor, and consequently from (3)
0' = (1 - D)O'; 0'= (1 - D)Eo:(E - £P) (5a,b)
By comparing (2) and (5), it follows that E = (1 - D)Eo, and
f, =[1 - D,(K,)]],(K,); fc =[1 - Dc(Kc)]IC<Kc) (lla,b)
where 0 s D, < 1 and 0 S Dc < 1 = tensile and compressive
degradation damage responses, respectively; and], and Ie =
functions that represent the effective-stress responses in the
uniaxial tensile and compressive stress functions, respectively.
Because D, (or De) is a function of only K, (or Ke), the single
degradation damage variable D in (5) can be used to describe
JOURNAL OF ENGINEERING MECHANICS / AUGUST 1998/893
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.
both tensile and compressive degradation responses if it is de-
fined as
Because P is a first-degree homogeneous function and with
the definition of Din (12), substituting the relations in (5) and
(13) into (9) leads to
F(a,j"lc):50 (14)
By introducing a new vector function C = C(K) =
C(J,(K,),!c(KJ), the yield function in (14) is rewritten as a
function of stress and damage variables F(a, K) = P(a,J.,lc)
(20)
(22)
The uniaxial version of the damage variable, denoted by 0 :5
KIl :5 I, is defined as
L
O'
KI( = ..!. (T1( eP) deP
gl( 0
is assumed to be determined by the plastic strain. In the Bar-
celona the model, the relation between the uniaxial stress, de-
noted by (Til, and the corresponding scalar plastic strain, de-
noted by f!, is assumed as
(TI( =/1(0[(1 + al()exp(-b~) - all exp(-2b~)] (19)
where /1<0 = initial yield stress, defined as the maximum stress
without damage; and all and bll = constants. The function in
(19) is factorized as shown in (II), assuming that the degra-
dation also takes an exponential form
where dll = a constant. The effective stress is then given by
al( =/1(0 [(1 + al()(exp(-bl(eP»I-(dMlb,.J - al(exp(-b~»2-(dM'btrl]
(21)
(15)
(13a,b)
F(i.7, K) :5 0
f, = (I - D)j,; fc = (I - D)lc
D = D(K) = 1 - (I - D,)(I - Dc) (12)
Note that D in (12) satisfies the condition 0 :5 D < 1 and
becomes D, if Dc = 0 (uniaxial tensile case) and Dc if D, = 0
(uniaxial compressive case). Accordingly, without losing the
consistent results in uniaxial responses, (11) can be replaced
by
in which F also is an isotropic, first-degree homogeneous sca-
lar function with multiple-hardening evolution.
With this framework, the elastoplastic responses in the pres-
ent plastic-damage model are described only in terms of a
and K
where
(23)
The total stress is determined by evaluating the degradation
damage from (12)
The present constitutive system, which is constructed by the
relations in (16)-(18), can be implemented effectively for nu-
merical computation, because the solution for the elastoplastic
responseis separate from the degradation damage process.
i.7 = Eo:(e - eP ) E {i.7IF(i.7, K) :5 O} (l6a)
F/ = ~yo;)<I>(i.7) (l6b)
K = ~H(i.7, K) (l6c)
In addition to (16), the following relations from Kuhn-Tucker
complementary conditions and the plastic consistency condi-
tion are used as loading/unloading conditions:
(27)
(24)
(25b)
(25a)
(200)
(26b)
stress in terms
The quantity gil = dissipated energy density during the entire
process of microcracking. Because the capacity for the dissi-
pated energy per unit volume cannot be given as a material
property, gil should be derived from other known material
properties such as the fracture energy. Assuming that the frac-
ture energy in the uniaxial tensile state Gil and its counterpart
in the uniaxial compressive state are given as a material prop-
erties gil is the specific fracture energy normalized by the lo-
calization zone size, also referred to as the characteristic length
Ill, which leads to gil = GJIIl (Lubliner et al. 1989; Oliver
1989). To maintain objective results at the structural level, III
should be an objective value or assumed as a material property.
It is still an open issue whether the counterpart of the fracture
energy in compressive failure is a material property (Chen
1994).
By substitution of (19) into (23), the relationship between
the constants is obtained
/1(0 ( al()gl( = bl( 1+"2
Combining (19) and (22) gives the uniaxial
of KI(
where <l>1l = 1 + all(2 + all)KIl' Similarly, the effective stress
and degradation damage variable are written in terms of KIl
0'1( =ll(KI()
0'1( =/1(0 [CJ (1 + al( - V<!>I(KI(» r-<d,.lbMJ) V<!>I(KI()
DI( = 1 - [CJ (1 + al( - V<!>I(KI(»r,.l~)
To model the different damage states for tensile and com-
pressive loaing, the damage variable K in (8) is defined using
K, and Kco which are independent state variables: K = [K, Kc].
(18)
(17a-c)
a = [1 - D(K)]a
~ ~ 0; ~F = 0; ~F = 0
CONSTITUTIVE RELATIONS
Evolution of Damage
The damage in quasi-brittle materials can be defined by
evaluating the dissipated fracture energy required to generate
microcracks. In the Barcelona model (Lubliner et aI. 1989;
Oller et aI. 1990), one scalar damage variable, which is the
combination of tensile and compressive damage, represents all
damage states. Although a scalar damage variable model can
simulate monotonic loading problems, it is not appropriate for
modeling cyclic behavior of concrete because the material ex-
periences a different damage process under tensile and com-
pressive loading. In the present study, two damage variables,
one for tensile damage and the other for compressive damage,
are defined independently, and each variable is factorized into
the effective-stress response and stiffness degradation response
parts.
Considering a uniaxial tensile or compressive stress state,
the state variable X E {t, c} is introduced. The state is uniaxial
tensile for X = t and uniaxial compressive for X = c. The stress
894/ JOURNAL OF ENGINEERING MECHANICS / AUGUST 1998
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.
where
From (22) and (25), the damage evolution equation for a uni-
axial state is given by
(36)
/
/
/
/
/
""""""""/
/
""(fby·lby)
IX = fhO - feO
2!hO -!eO
where 0 :S IX < 0.5 and (3) the uniaxial initial yield tensile
stress ftO, which generally is assumed to be the same as the
maximum tensile strength, is used with feO to give
1~a(a/l +,j3J2> = Co
FIG. 1. Lubllner (1989) Yield Function In Plane Stress Space
tion of the yield surface in the principal stress space is shown
in Fig. 1.
The parameters in (35) are evaluated so that the yield func-
tion accurately represents the experimental data. In the original
model, three independent material strength data for the plane
stress condition are used to determine the following parame-
ters: (1) The initial cohesion Co is given by assuming the uni-
axial compressive loading state Co = -!co, where!co is the uni-
axial initial yield compressive stress; (2) by considering the
biaxial initial yield compressive stress fbo, with !co the coeffi-
cient IX is given by
I ry
(28)
(31)
(33)
0 if & = 0
(t (~/») (30)
r(6-)= otherwise
(t I~/I)
. I ~ ). PKl( =- Jl«Kl( £
gl(
To extend the uniaxial version of the damage evolution equa-
tion to the multidimensional case, the scalar plastic strain rate
f! in (28) is assumed to be evaluated in the general three-
dimensional case by the relation
if = >"V~cI>(&) (32)
By substituting (32) into (31), the damage evolution equation
is given in a form similar to (8)
where (x) = (I x I + x)/2 denotes the Macaulay bracket function.
Substitution of (29) into (28) gives an evolution equation in
the form of
h(6-, K) = [r(8o)fr(K,)/g, 0 .,0 ]
o 0 (1 - r(a»!c(Kc)/gc
in matrix form. For the eigenvalue matrices & and fp, it is
assumed that the diagonal elements are the eigenvalues in al-
gebr~ic order. The flow rule in (6) can be modified into one
for Ff
£P =8'l(r(8o)£~ + 8Cl«(1 - r(8o»£~n (29)
where 8 = Kronecker delta; and £~ax and £~In = the algebrai-
cally maximum and minimum eigenvalues, respectively, of the
plastic strain tensor eP • The scalar quantity r(6') is a weight
factor, such that 0 :S r :S I, and it is defined as
where:A = h·Vcl-<I>(&-).
Rate-Independent Plasticity
The yield function [see (15)] for frictional materials like
concrete and rock can be rewritten as
where C = cohesion; and F = a scalar function of invariants of
the stress tensor. To make a yield surface more realistic for
concrete materials, many models include J3 = 1/38:8:8 term, or
the Lode angle, as well as /1 = tr a and J2 = 1/28:8 (Chen
1994), where 8 denotes the deviatoric stress tensor. In the Bar-
celona model, the algebraically largest principal stress, denoted
by doh is included instead of the Lode angle. This term gives
a yield function that better represents the two different behav-
iors in the tensile region and the compressive region. The yield
function for the Barcelona model is
(37)
(38)
(39a)
F(a) =C(K)
J3 =!co (IX - 1) - (1 + IX)
fro
The biaxial tensile strength depends on IX and ~, and it is
always a little less (approximately 1.3%) than the uniaxial ten-
sile strength.
In the Barcelona model, the cohesion c is only a variable
representing material strength, and the damage variable only
affects the cohesion. Consequently, (34) can be expressed as
Eq. (38) models isotropic hardening behavior in the classical
plasticity sense, which gives good results for monotonic load-
ing. Eq. (38) is not appropriate for modeling the cylic behavior
of concrete because the evolution of one strength hardly af-
fects the other strength, so it is necessary to use two cohesion
variables. A yield function with two cohesion variables is ob-
tained by defining ~, which is a constant in the Barcelona
model, and the cohesion such that
(34)
(35)
F(a) - C = 0
where IX and ~ = dimensionless constants; and amax = alge-
braically maximum principal stress. A plane stress cross sec-
Cc(K)J3 =-- (1 - IX) - (1 + IX)
C,(K)
(39b)
JOURNAL OF ENGINEERING MECHANICS I AUGUST 1998 I 896
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.
(42)
Because al is not a material property, which can be measured
explicitly, it is convenient to evaluate a, as a function of d,lb,
and D, such that
For the tensile case, a different approach is used because
the stiffness degradation at the maximum tensile stress usually
is assumed to be zero. Assuming that the degradation value in
a uniaxial state, denoted by D" is given at a, =/,c/2, then from
(25)
(48)
(49)
(50)
(51)
~ _ 10g(1 - Dc)
be - log ( 1 ~cae)
[2(1 - Dl"d, -1]
a, =2(1 - D,)b,ld'[(1 - D,)b,ld, - 1]
~~ l+a,+~
V <!>(K,) = 2
Substitution of (49) into (27) gives
d, 10g(1 - D,)
b, =10g[(1 + a,) - ~] - 10g(2a,)
APPLICATIONS TO CONCRETE
(39c)
F(u, K) =_1_ [WI +~ + ~(K)<&max)] - Cc(K) (40)
1 - a
where the invariants are now functions of the effective stress.
The flow rule gives the relation between the plastic flow
direction and the plastic strain rate. In contrast with metals,
the nonassociative flow rule is necessary to control the dila-
tancy in modeling frictional materials (Chen and Han 1988).
Because the yield surface in (40) is a combined geometric
shape from two different Drucker-Prager type functions, a
Drucker-Prager type function is used as the plastic potential
function in the present model
<I> =VU; + all; <I> = IIsll + all (41a,b)
where Iisil =~ denotes the norm of the deviatoric stresses.
The parameter ap is chosen to give proper dilatancy. The plas-
tic strain rate is obtained from (6) and (41)
eP = i.. (11:11 + apI)
where c, = a tensile cohesion; and Cc = a compressive cohesion.
For the assumed relations between the cohesion and damage
variables, given by c, = 1,(K,) and Ce = -lc(Kc), the yield func-
tion in (15) becomes
By substituting (46) into (27)
An efficient numerical algorithm for the present model has
been developed (Lee 1996) and implemented as a part of the
finite-element program FEAP (Taylor 1992). The derived con-
sistent algorithmic tangent stiffness (Lee 1996) is used in place
of the continuum tangent stiffness to maintain quadratic con-
vergence rate in solving the nonlinear equations by the New-
ton-Raphson method (Simo and Taylor 1985, 1986). In this
section, several numerical examples of the present plastic-
damage model for concrete are presented. Two-dimensional,
plane stress quadrilateral isoparametric elements are used. Dis-
placement control is used to apply the loads, and for all cases
the Poisson ratio is 0.18.
Monotonic Uniaxial Loading Test
To compare the present rate-independent plastic-damage
model with the existing experimental results (Karsan and Jirsa
1969; Gopalaratnam and Shah 1985), the following material
properties are used: (1) For the tensile case, Eo = 3.1 X 104
MPa, the maximum tensile strength /: = 3.48 MPa, G, = 40
N/m, i, = 82.6 mm, and up = 0.2, D, = 0.5; and (2) for the
compressive case, Eo = 3.17 X 104 MPa, the maximum com-
pressive strength/; = -27.6 MPa, Ge = 5,690 N/m, ic = 82.6
mm, up =0.2, and Dc =0.4. Fig. 2 shows the simulated tensile
and compressive stress-strain curves, respectively, which agree
well with the experimental data, using a single-element patch
(mesh A; 82.6 X 82.6 mm). It is noted that the averaged strain
over the assumed effective specimen size (82.6 mm) is used
for reproducing the experimental results.
In Fig. 2(a), the experimental data shows a very steep post-
peak curve because the data were measured by a wide strain
gauge over the localized zone. If the averaged stress-strain
relation in a smaller localized zone is considered, the absolute
value of the negative tangential stiffness after the ultimate
stress point will decrease. The numerical model cannot yield
such a large negative stiffness shown in the experimental result
for tensile loading because the assumed stress curves are con-
structed as the difference of two exponential functions [see
(19)].
Quasi-brittle materials, such as soil, rock, and concrete, are
known to undergo a significant volume change (dilatancy),
which is caused by the large amont of the inelastic strain. The
(44)
(46)
(47)
s(6') =So + (1 - so)r(&)
Stiffness Degradation and Calibration with
Experimental Data
From experimental tests of concrete it is clear that the deg-
radation of stiffness from microcracking occurs in tension and
compression and becomes more significant as the strain in-
creases. Under cyclic loading, the mechanism of the stiffness
degradation becomes complicated because of the opening and
closing of microcracks. The modeling of the crack opening!
closing behavior can be implemented by elastic-stiffness re-
covery duirng the elastic unloading process from tensile state
to compressive state. The degradation damage variable defined
in (12) is modified by multiplying a parameter, such that 0 $
S $ 1, to D,
D(K, 6') = 1 - (1 - D e(K»(1 - s(6')D,(K» (43)
The parameter s is chosen to represent the stiffness recovery
as follows:
where 0 $ So $ 1 is a constant to set the minimum value of
s. Note that r{fr) = r{a). By substituting (43) into (5), the total
stress is given by
u = (1 - Dc(K»(1 - SD,(K»Eo: (e. - £P) (45)
The ratio d~/b~, in (26) and (27) can be evaluated by spec-
ifying the degradation values in each uniaxial state. To perform
the calibration for compressive cases, it is assumed that a value
of the elastic-stiffness degradation Dc can be obtained from
experiments at the maximum compressive stress. If K; is the
compressive damage value at which the effective uniaxial
function lc is at its maximum absolute value, it follows from
(25) that
)]
d<lb<
1 + ae+a ---
e 2
which leads to the relation between dclbe and Dc such that
896/ JOURNAL OF ENGINEERING MECHANICS / AUGUST 1998
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.
0.0010.0001o
STRAIN
(a)
-0.0001
MUHB
oL............................L................-'-........L..................--o......L_........................J
-0.001
'r--"':"""":---:-....,...----.-----r-----,
~---
0.0001 0.0002 0.0003 0.0004 0.0001 0.0001
STRAIN £11
(a)
4r------------------.
-0.001o0.001
EXPERIMENTAL
NUMERICAL (<Xp-O.15)
NUMERICAL (<Xp-O.2)
oL--'-........--o........'--.l.-__..............--.........._ .......-o...-I
0.01
STRAIN
(b)
FIG. 3. Numerical Solutions of Biaxial Loading with Uniaxial
Loading: (a) Tensile Case; (b) Compressive Case Compared
with experimental Result (Kupfer et al. 1969)
havior. The hysteresis on reloading cannot be simulated by the
model because of the rate-independent elastic loading/unload-
ing assumption.
In the third case for cyclic uniaxial loading, mesh A (50 X
50 mm) is subjected to cyclic tensile-low compressive load.
The material properties used for this case are different from
the previous case because the test material had very high
Young's modulus (Reinhart 1984): Eo = 3.9 X 104 MPa,f: =
3.2 MPa, f: = -40 MPa, G, =45 N/m, Ge =450 N/m, 1, =50
mm, Ie =50 mm, ctp = 0.2, D, =0.8, and De =0.5. As shown
in Fig. 5, the numerical simulation agrees well with the result
of the experimental test for cyclic tensile-low compressive
loading. The microcrack opening and closing are simulated by
the degradation function. Because there is a limitation in mod-
eling such discrete phenomenon by a continuum model, in-
cluding the present plastic-damage model, the sudden recovery
of stiffness near the final cyclic loading stage, which is ob-
served experimentally, is difficult to simulate precisely. The
overall response of the model for cyclic tensile-compressive
loading is close to the experimental results.
As the last example of the cyclic loading, high tensile-high
compressive full cyclic load is imposed on mesh A (82.6 X
82.6 mm) with the following material properties: Eo =3.1 X
l~ MPa; f: = 3.3 MPa,f: = -30 MPa, G, = 60 N/m, Ge =
6,000 N/m, I, = 82.6 mm, Ie = 82.6 mm, ctp = 0.2, D, = 0.7,
and De = 0.4. The result is illustrated in Fig. 6 with the loading
path O-A-B-C-D-E-F-G. During the tensile unloading-com-
pressive loading path (paths A-B and D-E), the stiffness re-
covery is shown. The compressive strength after tensile load-
ing is not degraded (see path A-B-C), which means that the
effect of compressive damage on the yield surface and the
elastic stiffness is not affected considerably by the tensile dam-
age. On the other hand, the tensile strength depends on the
compressive damage because the stiffness degradationaccu-
JOURNAL OF ENGINEERING MECHANICS / AUGUST 1998/897
-0.001
M!!SHA
-0.001
06----.....L---..L.-""---1---...L.-""----I
0.000
-t5 ·11
~ ·10
oS
·25
dilatancy can be modeled by the numerical model. It is noted
that the dilatancy is controlled by the ctp parameter in the plas-
tic potential function. Considering a numerically admissible
range of ctp, which is described in Lee (1996), and the pa-
rameters of the yield function, ctp = 0.2 is used as a realistic
value in these examples.
-30 r-------------------,
-0.002 -o.ooa
STRAIN £11
(b)
FIG. 2. Numerical Solutions of Monotonic Uniaxial Loading
Compared with Experimental Results: (a) Tensile Case (Gopa-
larantnam and Shah 1985); (b) Compressive Case (Karsen and
Jlrse 1969)
Biaxial Loading Test
A four-element mesh (mesh B; 82.6 X 82.6 mm) is used to
test the model under biaxial loading with same material data
as for the uniaxial loading test. The results for biaxial tensile
loading and compressive loading are shown in Figs. 3(a and
b), respectively, illustrating the typical effects of biaxial load-
ing. The compressive results with two different values of ctp
(0.15 and 0.2) are compared with the experimental result
(Kupfer et al. 1969). The result with ctp = 0.15 shows good
agreement with the experimental one for all components of the
strain, whereas ctp = 0.2 yields a slightly different response for
the out-of-plane strain caused by the larger dilatancy.
Cyclic Loading Test
Mesh A (82.6 X 82.6 mm) is used for a cyclic tensile load-
ing test, and the comparison with the experimental data (Go-
palaratnam and Shah 1985) is illustrated in Fig. 4(a). Mesh A
also is subjected to compressive cyclic loading with the same
material data for the monotonic uniaxial compressive case, ex-
cept Eo = 3.17 X 104 MPa. Fig. 4(b) shows the numerical
result compared with the experimental data (Karsan and Jirsa
1969). For both cases the degradation of stiffness is simulated
at each unloading/reloading cycle as well as the softening be-
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.
-35 ,-----------------__...,
A
·20
10 ,---"""T'----T'"-----r----...,
0.00010.00040.0002 0.0003
STRAIN ell
(a)
0.0001
3
2
4r------------------,
-30L.-_........::a.l'-==:::::..._--I. ---J._---'-_~
-0.003 -0.002 -0.001 0.000 0.001
STRAIN ell
FIG. 6. Numerical Solution for Full Cyclic Loading (Path: O-A-
SoC-D-E-F-G)
The sensitivity of the numerical solution to element size is
examined to check the ability of the characteristic length (Ol-
iver 1989) to provide mesh objectivity. Three (mesh C, D, and
E) with the same specimen dimensions (50.8 X 25.4 mm) but
different element sizes are tested to examine the sensitivity to
mesh size, as illustrated in Fig. 7. The material properties are
the same as those used in the monotonic uniaxial tensile load-
ing case except the characteristic length, which is defined as
the element size for each case (mesh C, 25.4 mm; mesh D,
12.7 mm; and mesh E, 6.35 mm). A slightly smaller tensile
strength was imposed on the left-most element to force strain
localization.
Fig. 7 shows the load-displacement curves measured at the
right end. Although the left-end element is undergoing strain
localization, the other elements unload elastically after the
peak strength is reached. Because the ratio of the softening
bandwidth to the length of elastic unloading zone is different
in each case, the global load-displacement responses cannot
be identical, even though the dissipated energy for all three
cases is equal. Nevertheless, the global responses approach a
limit, as the softening bandwidth becomes smaller.
CONCLUSIONS
Mesh-Sensitivity Tests
The load path in Fig. 6 resembles the hysteresis curve of
the experimental results in Fig. 4(b), because energy is dissi-
pated during the tensile reloading (path C-D-E). This example
shows that the present model can simulate quasi-brittle mate-
rial behavior subjected to cyclic loading for which the hyster-
esis curve in the tensile region is negligible compared with
the compressive counterpart.
MISHA
-0.004 -0.001-0.001
oK-...J.~~~~~~~:::::!=-----.J
0.000
-30
~ ·25
:!
·20
6
~ ·15
l/l~ ·10
-0.002 -0.003
STRAIN e11
(b)
FIG. 4. Numerical Solutions of Cyclic Uniaxial loading C0m-
pared with Experimental Results: (a) Tensile Case (Gopalaratnam
and Shah 1985); (b) Compressive Case (Karsan and Jirsa 1969)
60
l/l
~
2
4,----,---.....---r----r---r----,.--...,
·2
-4L-_......L.--o_J.-........-1.--' J.-........-1.--' .L-__-..l
·2E-GS 0 2E-GS 4E-GlI 1I!-GlI 8E-GS 0.0001 0.00012
STRAINe11
FIG. 5. Numerical Solution of Cyclic Tensile-Compressive
Loading Compared with Experimental Result (Reinhart 1984)
mulated by the compressive damage is not recovered during
tensile unloading (path CoD). This is physically reasonable,
because the compressive failure mainly is caused by dilatancy,
which affects both compressive and tensile strength. It is noted
that the stiffness recovery takes place only when the tensile
cracks are closed.
The constitutive relations for a new plastic-damage model
based on plasticity and continuum damage mechanics have
been developed in this paper. The concept of the fracture-en-
ergy-based multiple-hardening variables is used to represent
tensile and compressive damage independently. The constitu-
tive system for the elastoplastic responses is described com-
pletely by the effective-stress and damage variable, which can
lead to a decoupled algorithm for the effective-stress compu-
tation and degradation evaluation. A thermodynamically con-
sistent scalar model is used to simulate the stiffness degrada-
tion and recovery.
898/ JOURNAL OF ENGINEERING MECHANICS / AUGUST 1998
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.
100 .------------------......,
(54)
APPENDIX I. THERMODYNAMICAL
CONSIDERATIONS
(58)
(61)
(62)
0' = (1 - J) iJ'I'o
1 - w iJ£'
-q, - PTl9 + (l - w)O':e _.!. (q·V6) ~ 0
6
iJ'I'O'K s 0
iJK
The stress defined in (61) becomes identical with (45) by as-
suming that '1'0 for the present model is given by
'I' ='1'(£', K, d, 6); 'I' = (l - d)'I'0(£', K, 6) (59a,b)
Here, '1'0 denotes the total potential energy function of material
without the damage in stiffness degradation and d represents
unrecoverable damage in the elastic stiffness, which is defined
ip the present model as d = 1 - (1 - D,)(l - Dc). Note that
d is a nonnegative function, which implies that the damaging
process is unrecoverable. Assuming an isothermal state and the
equality as Tl = -(l/p)(iJ'I'/iJ6) and then by the chain rule and
(52), (58) reduces to
[(1 - w)O' - (l - a) ~~:] :e' + (1 - w)O':eP
J (
iJ'I'o .) .:.
-(1-) ""K'K +d'l'o~O
u (60)
In (60), (l - w)a: eP ~ 0 b«(cause of the con\;exity of the
plastic potential function and d'l'0 ~ 0 because d and '1'0 are
positive. Hence, to satisfy (60), the following relations are re-
quired as sufficient conditions:
where w =a parameter such that 0 s w < 1.
The Helmholtz free energy per unit volume 'I' can be de-
fined using e,
'I' = e, - PTl6 (55)
where Tl = entropy per unit mass; and 6 = thermodynamic
temperature. Differentiating (55) with respect to time and com-
bining with (53) gives
q, = e - e/ - p,,6 - PTl9 (56)
By substituting (54) into (56), (52) becomes
ph = V·q + p,,6 + PTl9 + q, - (l - w)O':e (57)
Substitution of (57) into the second law of thermodynamics
gives the Clausius-Duhem inequality
In continuum damage mechanics the Helmholtz free energy
for damaged material is written as
0.02
(a)MESHC
(c) MESH E
(b) MESH D
MESHE
MESHC
MESH 0
71
0.001 0.01 0.011
DISPLACEMENT (mm)
FIG. 7. Load versusDisplacement Curves for Different
Meshes
Through the examples the present plastic-damage model is
shown to give numerical simulations that agree well with ex-
perimental data for concrete under cyclic loading as well as
monotonic loading. Stiffness degradation and recovery, caused
by crack opening/closing, are simulated properly. This shows
that the proposed formulation with cooperation of plasticity
provides a simple but effective method. The model also shows
a good control of dilatance, which is very important to sim-
ulating concrete structures under multiaxial loading. Conse-
quently, it can be concluded that the present plastic-damage
model can represent accurately the behavior of concrete struc-
tures under monotonic and cyclic loading.
The equations of a rate-independent constitutive model re-
sult in an ill-posed initial-boundary value problem for material
softening (Hill 1958; Needleman 1988; de Borst et al. 1993;
Neilsen and Schreyer 1993). Consequently, uniqueness of so-
lution for the problem is not guaranteed. One result of the ill-
posedness is mesh sensitivity in a finite-element analysis. The
characteristic length concept is used to maintain the objectivity
of the numerical solution. It should be noted that use of the
nonassociative flow rule gives another source of ill-posedness
and this kind of ill-posedness can be avoided by choosing
appropriate parameters in the plastic potential function (Ru-
nesson and Mroz 1989; Neilsen and Schreyer 1993).
(63)The first law of thermodynamics is expressed as
e = 0': e - V· q + ph (52)
where e =internal energy per unit volume; q =heat flux vec-
tor; P = density; and h = source heat per unit mass. The in-
ternal energy is assumed to be divided into two parts, an ex-
plicit internal energy, denoted bye.. and an implicit internal
energy, denoted bye/
(53a,b)
The explicit part is the internal energy itself in the classical
approach. The implicit part is introduced to explain the stiff-
ness recovery in the thermodynamic sense. It can be assumed
that the implicit internal energy rate balances a portion of the
stress power
1
'l'o='2£':Eo:£' + 'I'~(K) + 'I'~(6)
where 'I'~ =contribution of the damage variable (the harden-
ing variable in the classical plasticity context) to '1'0 and de-
fining w =(l - s)D,I(l- sD,). It is noted that because 0 s
w < 1, e/ is always smaller than the mechanical input power,
the stress power. The dissipated energy during the tensile and
compressive damage process associated with the effective
stress can be used for the definitions of 'I'~
'I'~(K) = g, + gc - (f' C,(K) dK + f< Cc(K) dK) (64)
Because (iJ'I'oIiJK) = (iJ'I'~/iJK), and K, and Kc are monotonically
increasing functions, the condition in (62) becomes satisfied
JOURNAL OF ENGINEERING MECHANICS / AUGUST 1998 / 899
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.
K
a;, a:
r =
s =
s =
a, 13
a p
£, e", eP
aco a, = parameters for uniaxial compressive and tensile
strength functions, respectively;
C, Co = cohesion and initial cohesion, respectively;
Ce, Ct compressive and tensile degradation damage var-
iables, respectively;
D.p = continuum elastoplastic tangent stiffness tensor;
D, Dc, D, = total, compressive, and tensile degradation dam-
age variables;
Dc, D, = measured elastic stiffness degradation for uniaxial
compression and tension, respectively;
E, Eo = elastic and undamaged elastic stiffness tensors;
Eo = Young's modulus;
F = yield function;
fbO' fcQ' ftO = initial yield stresses for biaxial compression, uni-
axial compression, and uniaxial tension, respec-
tively;
!C, f, = uniaxial compressive and tensile strength func-
tion, respectively;
compressive counterpart of fracture energy and
fracture energy, respectively;
8e, 8, = compressive and tensile specific fracture energy
parameters, respectively;
/1 = first invariant of stress tensor;
12 = second invariant of deviatoric stress tensor;
Ie, It = compressive and tensile characteristic length pa-
rameters, respectively;
stress weight factor function;
deviatoric stress tensor;
parameter for stiffness recovery;
coefficients of yield function;
=coefficient of plastic potential function;
= total, elastic, and plastic strain tensors, respec-
tively;
TJ = entropy per unit mass;
6 = thermodynamic temperature;
damage variable vector;
Ke, K t = compressive and tensile damage variables, re-
spectively;
>... = plastic consistency parameter;
p =density of material;
0', a = stress and effective stress tensors;
iT, 6- = eigenvalue matrices for stress and effective stress
tensors, respectively;
maximum compressive and tensile strength val-
ues, respectively;
<I> plastic potential function; and
qr, qro = Helmholtz free energy per unit volume for dam-
aged and undamaged materials, respectively.
Pramono, E., and Willam, K. (1989). "Fracture energy-based plasticity for-
mulation of plain concrete." J. Engrg. Mech., ASCE, 115(6), 1183-
1203.
Reinhardt, H. W. (1984). "Fracture mechanics of an elastic softening
material like concrete." Heron, 29(2), 1-42.
Runesson, K., and Mroz, Z. (1989). "A note on nonassociated plastic
flow rules." Int. J. of Plasticity, 5, 639-658.
Simo, J. C., and Ju, J. W. (1987). "Stress and strain based continuum
damage models-I. Formulation." Int. J. Solids and Struct., 23,
821-840.
Simo, J. C., and Taylor, R. L. (1985). "Consistent tangent operators for
rate-independent elastoplasticity." Computer Methods in Appl. Mech.
and Engrg., 48,101-118.
Simo, J. C., and Taylor, R. L. (1986). "A return mapping algorithm for
plane stress elastoplasticity." Int. J. for Numer. Methods in Engrg., 22,
649-670.
Taylor, R. L. (1992). "FEAP: A finite element analysis program for en-
gineering workstations." Rep. No. UCBISEMM-92 (Draft Version),
Department of Civil Engineering, University of California, Berkeley,
Calif.
Yazdani, S., and Schreyer, H. L. (1990). "Combined plasticity and dam-
age mechanics model for plain concrete." J. Engrg. Mech., ASCE,
116(7), 1435-1450.
APPENDIX III. NOTATION
The following symbols are used in this paper:
(65)
ACKNOWLEDGMENTS
which shows that the present plastic-damage model is ther-
modynamically consistent.
APPENDIX II. REFERENCES
This work has been supported by the U.S. Bureau of Reclamation,
under Agreement No. 1425-5-FC-81-08970. The writers thank Larry Nuss
for his support. The writers also thank Professor Jacob Lubliner for help-
ful discussions.
Cervera, M., Oliver, J., and Faria, R. (1995). "Seismic evaluation of
concrete dams via continuum damage models." Eanhquake Engrg. and
Struct. Dynamics, 24, 1225 -1245.
Chen, W. F. (1994). Constitutive equations for engineering materials, Vol.
2, Elsevier, New York, N.Y.
Chen, W. F., and Han, D. J. (1988). Plasticity for structural engineering.
Springer-Verlag, New York, N.Y.
de Borst, R., Sluys, L. J., Muhlhaus, H.-B., and Pamin, J. (1993). "Fun-
damental issues in finite element analyses of localization of deforma-
tion." Engrg. Computations, 10, 99-121.
Gopalaratnam, V. S., and Shah, S. P. (1985). "Softening response of plain
concrete in direct tension." ACI J., 82(3), 310-323.
Hansen, N. R., and Schreyer, H. L. (1994). "A thermodynamically con-
sistent framework for theories of elastoplasticity coupled with dam-
age." Int. J. Solids and Struct., 31 (3), 359 - 389.
Hansen, N. R., and Schreyer, H. L. (1995). "Damage deactivation." J.
Appl. Mech., 62, 450-458.
Hill, R. (1958). "A general theory of uniqueness and stability in elastic-
plastic solids." J. Mech. Phys. Solids, 6, 236-249.
Ju, J. W. (1989). "On energy-based coupled elastoplastic damage theo-
ries: Constitutive modeling and computational aspects." Int. J. Solids
and Struct., 25(7), 803-833.
Kachanov, L. M. (1986). Introduction to continuum damage mechanics.
Martinus Nijhoff Publishers, Dordrecht, The Netherlands.Karsan, 1. D., and Jirsa, J. O. (1969). "Behavior of concrete under com-
pressive loadings." J. Struct. Div., ASCE, 95(12), 2535-2563.
Kupfer, H., Hilsdorf, H. K., and Rusch, H. (1969). "Behavior of concrete
under biaxial stresses." ACI J. 66(8), 656-666.
Lee, J. (1996). "Theory and implementation of plastic-damage model for
concrete structures under cyclic and dynamic loading," PhD disserta-
tion, Department of Civil and Environmental Engineering, University
of California, Berkeley, Calif.
Lee, J., and Fenves, G. L. (1994). "Numerical implementation of plastic-
damage model for concrete under cyclic loading: Application to con-
crete dam." Rep. No. UCB/SEMM-94/03, Department of Civil Engi-
neering, University of California, Berkeley, Calif.
Lemaitre, J., and Chaboche, J.-L. (1990). Mechanics of solid materials.
Cambridge University Press, New York, N.Y.
Lubliner, J., Oliver, J., Oller, S., and Onate, E. (1989). "Aplastic-damage
model for concrete." Int. J. Solids and Struct., 25(3), 299-326.
Mazars, J. (1986). "A model of a unilateral elastic damageable material
and its application to concrete." Fracture toughness and fracture en-
ergy of concrete, F. H. Wittmann, ed., Elsevier Science Publishers,
Amsterdam, The Netherlands, 61-71.
Mazars, J., and Pijaudier-Cabot, G. (1989). "Continuum damage theory-
application to concrete." J. Engrg. Mech., ASCE, 115, 345-365.
Mehta, P. K., and Monteiro, P. 1. M. (1993). Concrete: Structure, properties,
and materials, 2nd Ed., Prentice-Hall, Inc., Englewood Cliffs, N.J.
Murray, D. w., Chitnuyanondh, L., Rijub-Agha, K. Y., and Wong, C.
(1979). "Concrete plasticity theory for biaxial stress analysis." J.
Engrg. Mech., ASCE, 105, 989-1006.
Needleman, A. (1988). "Material rate dependence and mesh sensitivity
in localization problems." Computer Methods in Appl. Mech. and
Engrg., 67, 69-85.
Neilsen, M. K., and Schreyer, H. L. (1993). "Bifurcation in elasto-plastic
materials." Int. J. Solids and Struct., 30(4), 521-544.
Ohtani, Y., and Chen, W. F. (1988). "Multiple hardening plasticity for
concrete materials." J. Engrg. Mech., ASCE, 114(1), 1890-1910.
Oliver, J. (1989). "A consistent characteristic length of smeared cracking
models." Int. J. Numer. Methods in Engrg., 28,461-474.
Oller, S., Onate, E., Oliver, J., and Lubliner, J. (1990). "Finite element
nonlinear analysis of concrete structures using a plastic-damage
modeL" Engrg. Fracture Mech., 35(1/2/3), 219-231.
Ortiz, M. (1985). "A constitutive theory for the inelastic behavior of
concrete." Mech. of Mat., 4,67-93.
900 I JOURNAL OF ENGINEERING MECHANICS I AUGUST 1998
 J. Eng. Mech., 1998, 124(8): 892-900 
D
ow
nl
oa
de
d 
fr
om
 a
sc
el
ib
ra
ry
.o
rg
 b
y 
Io
w
a 
St
at
e 
U
ni
ve
rs
ity
 o
n 
01
/2
4/
19
. C
op
yr
ig
ht
 A
SC
E
. F
or
 p
er
so
na
l u
se
 o
nl
y;
 a
ll 
ri
gh
ts
 r
es
er
ve
d.