INTERNATIONALJ O U R N A LFOR N U M E R I C A L METHODS I N ENGINEERING,VOL.

38, 361 1-3633 (1995)

ANISOTROPIC MODELLING AND NUMERICAL

SIMULATION OF BRITTLE DAMAGE IN CONCRETE

SANJAY GOVINDJEE

Department of Civil Engineering, Division of Structural Engineering and Structural Mechanics, University of California,
Berkeley, C A 94720, U.S.A.

GREGORY J. KAY
University of California, Lawrence Liuermore National Laboratory, P.O. Box 808, Livermore, C A 94550, U.S.A.

JUAN C. SIMO*
Division of Applied Mechanics, Department of Mechanical Engineering, Stanford University, Stanjord, C A 94305, U.S.A.

SUMMARY
A framework for damage mechanics of brittle solids is presented and exploited in the design and numerical
implementation of an anisotropic model for the tensile failure of concrete. The key feature exploited in the
analysis is the hypothesis of maximum dissipation, which specifies a unique damage rule for the elastic
moduli of the solid once a failure surface is specified. A complete algorithmic treatment of the resulting
model is given which renders a useful tool for large-scale inelastic finite element calculations. A rather simple
three-surface failure model for concrete, containing essentially no adjustable parameters, is shown to
produce results in remarkably good agreement with sample experimental data.

KEY WORDS: brittle damage; softening; anisotropic damage; characteristic length; concrete modeling

1. INTRODUCTION
Degradation of the stiffness properties under loading (i.e. damage) is a complex phenomenon
displayed by many materials. Often knowledge of the actual micromechanical processes involved
is quite limited. Typically, experimental measurements are available that define an envelope of
admissible states in a given range of stress or strain space. This is referred to as the ‘failure surface’.
Points lying interior to the domain enclosed by the failure surface define states that do not lead to
damage evolution in the material. In some situations, it is also possible to characterize the
changes in shape experienced by the failure surface under sustained loading. This limited
information is often all that exists about the damage in the material.
For brittle materials such as ceramics, glasses, rocks or plain concrete, the stress response can
be characterized via linear elastic stress-strain relations. Damage in these materials is reflected at
the macroscopic level in the progressive degradation experienced by the elastic moduli under
sustained tensile loading, at states lying on (or above) the failure surface. The additional

*Deceased

CCC 0029-5981/95/213611-23 Received 14 December 1993

0 1995 by John Wiley & Sons, Ltd. Revised 13 December 1994
3612 S.GOVINDJEE. G. J. KAY AND J. C. SlMO

information needed to construct a complete constitutive model is an evolution law that character-
izes this degradation process; i.e., a damage rule. The pioneering work of Kachanov' furnishes the
first example of a damage rule, restricted in principle to either one-dimensional or isotropic
theories, via the concept of effective stress. Proposals aimed at generalizing this original idea to
account for anisotropy often introduce a fourth-order tensor that maps between the stress and
effective stress spaces, see e.g., work of Lemaitre' for a review of the extensive literature on this
subject. The precise identification of this tensor in actual materials, plain concrete in particular,
appears not to be readily accomplishable.
The purpose of this paper is to address, in detail, continuum and computational issues involved
in the modelling of brittle materials, specifically in plain concrete, by exploiting an alternative
framework that yields a fully anisotropic damage rule free from adjustable parameters for a given
failure surface. This approach employs an internal variable formalism that relies on two phenom-
enological assumptions which complement our limited information on the micromechanics
leading to damage evolution:
(i) The choice of the elastic moduli as the set of internal (damage) variables: The idea of
adopting the elastic moduli themselves or, equivalently, the compliance tensor as internal
variables first appears in works of Ortiz3 and Simo and J u . ~
(ii) The principle of maximum (damage) dissipation: This single hypothesis/assumption auto-
matically renders a fully anisotropic damage evolution law for the tensor of elastic moduli
which, as alluded to above, is free from any adjustable parameters once the failure surface
has been defined. The use of the principle of maximum dissipation as a means of uniquely
defining the damage rule first appears in Simo and J u . ~
It should be noted that the second hypothesis is also the basis for the classical model of
associative plasticity. In this context, assumption (ii) is often credited to von Mises' and is known
as the principle of maximum plastic work. The continuum model so constructed, therefore,
furnishes the counterpart in continuum damage mechanics of associative plasticity; and impor-
tantly it is completely defined once a failure surface is specified. Additionally, it is pointed out that
assumption (ii) is also known to be valid in the disparate area of the modelling of damage due to
Mullin's effect.6 The validity of the assumption can, of course, only be checked by examining the
predictions of so constructed models.
As in the case of softening plasticity, damage degradation of the elastic moduli leads to strain
softening and, therefore, leads inevitably to the appearance of strongly discontinuous solutions;
see Reference 7. The underlying mechanism can be traced to a local loss of Hadamard's strong
ellipticity condition; see e.g. Reference 8 for one of the many derivations of this well-known fact.
As a result, rate-independent models exhibiting softening render ill-posed the boundary value
problem for the quasi-static case (or the initial boundary value problem in the dynamic case). In
simulations, lack of well-posedness manifests itself in numerical solutions exhibiting a strong
mesh-dependence. Motivated in part by these undesirable features, a number of authors (e.g.
Reference 9) have questioned the common point of view which regards loss of stiffness as a true
material pqoperty of a continuum.
A number of techniques have been proposed to circumvent the preceding difficulties, from the
early approaches of Pietruszczak and Mroz'' or Baiant and Oh" based on the use of a charac-
teristic length to more recent methodologies that introduce a 'regularization parameter'. Repre-
sentative examples of the latter approach are the viscous regularization, advocated in Reference
12 and others, and regularization techniques based on the use of generalized continuum, as in
Reference 13 or Reference 14. An alternative approach that retains the rate-independent charac-
ter of the softening model and restores well-posedness, without resorting to a characteristic
 BRITTLE DAMAGE IN CONCRETE 3613

length, is described in the work of Simo et a1.15 Since the emphasis in this paper is placed on
modelling issues, amongst the many proposed techniques, only two methods are considered in
detail the viscous regularization and the characteristic length method in the form given in
Reference 16. The goal here is to illustrate how existing regularization approaches fit within the
proposed framework as opposed to a detailed discussion of their merits.
An outline of the remainder of this paper is as follows: In Section 2 we summarize the
continuum framework for the design of continuum damage models based on assumptions (1) and
(ii) (adapting from the presentation in Reference 17.) In Section 3 we consider in some detail
computational issues involved in the implementation of the general model. These include
numerical integration, inclusion of rate effects and an alternative implementation of the concept
of characteristic length. Section 4 describes the formulation of an anisotropic continuum damage
model for the tensile and shear cracking of concrete. A set of three simple damage surfaces are
postulated for concrete, which together with the damage evolution formulation of Section 2
renders a set of anisotropic evolution laws for the rank 4 stiffness tensor of the material and a set
of softening variables. These evolution laws are then integrated using the algorithmic formulation
of Section 3 and placed within a finite element code. Example calculations are shown to
demonstrate various features of the model and to establish comparisons with some available
experimental data. These examples demonstrate the full three-dimensional anisotropic features
inherent in the model.

2. SUMMARY OF THE GENERAL MODEL

Consider a continuum solid occupying a region B c R3. Since the interest here is in brittle
materials, the deformation of the solid is well described by the linear strain theory and the
assumption is made that no permanent (plastic) strains remain in the solid upon unloading. It
follows that the stress state is characterized via the stress-strain relations:
Q = C:E (1)
where Q is the Cauchy stress tensor, E is the infinitesimal strain tensor, and C is the rank 4 stiffness
tensor of the material. To account for damage degradation in the material, the tensor C at time t is
considered to be a function of the deformation history up to time t. According to the discussion in
the introduction, it is further assumed that the elastic domain (or admissible set of stress states) is
the set of stresses defined by M (given) damage functions 4kas

In this description, S stands for the set of symmetric rank 2 tensors and q is viewed as an internal
variable that determines the evolution in time of the shape of the admissible set with progressive
damage. An example of a damage criterion suitable for brittle materials such as plain concrete is
the Rankine condition.
Without any further information on the detailed micromechanics governing damage evolution,
a damage rule can be specified by introducing the additional assumption of maximum damage
dissipation, as follows: First, regard the elastic moduli themselves as internal variables, leading to
a free energy function for the material of the form
“ ( E , C, a) = ~ E:C+
: ES(a) (3)
where a is the internal variable conjugate to q, in the sense that q = - a,S(a). Use of the Legendre
transformation yields the dual of (3,involving the stress field Q and the internal variable q. which
is the starting point in or ti^.^ Equation (3), on the other hand, is taken as a point of departure in
3614 S. GOVINDJEE, G.J. K A Y AND J. C . SIMO

Simo and J u . ~The first term in (3) represents the elastic free energy in the system while the second
term is the free energy associated with the progressive degradation of the material.
The second step in the analysis involves the enforcement of maximum damage dissipation.
Recall that in the absence of thermal effects, the internal dissipation 9 is given byl8
9 : =- Y + a:E 2 0 (4)
where superposed dots indicate ordinary time differentiation. Using (1)-(3) and q = - a,S(ct)
yields
9 = ;cT:oi:cT + qi 2 0 (5)
where D = @ - ' denotes the material compliance tensor. The hypothesis of maximum dissipation
states that among all admissible states lying in E, the actual state of stress is the one which
maximizes (5). For materials that exhibit softening, this hypothesis is not strictly enforced.
Instead, the hypothesis is reduced to finding the critical point of dissipation since no such
maximum exists. To perform this constrained optimization problem, one first constructs the
associated Lagrangian

where Y k 2 0 are known as generalized Lagrange multipliers or consistency parameters. The

critical point of this Lagrangian is characterized by the classical optimality conditions (see e.g.
Reference 19, Section 10.8):

a,2 = 0, a q 9= o (74

Denoting by Q the tensor (or outer) product symbol, conditions (7) result in the damage flow
rules

supplemented by the Kuhn-Tucker unilateral constraint conditions (7a), i.e.

M
(Pk< 0, Yk 2 0 and Yk$k=O (9)
k= 1

Relations (9) give precise conditions for assessing active or inactive damage evolution in the
material. Their significance is identical to the loading/unloading conditions for multi-surface
plasticity in the form first stated in Koiter." In a strict sense, the damage rules and the derivation
sketched above are valid if the constraints are of the form 4 k ( ( T , q ) =fk(a)- hk(q), with &(a)
homogeneous of degree one. This is a situation often encountered in practice and, in particular, is
the case of interest discussed in Section 4.
In formal analogy with classical plasticity, observe that equations (8) and (1) imply the
following expression for the time rate of change of stress:
 BRITTLE DAMAGE IN CONCRETE 361 5

Table I. Rate independent, anisotropic damage model

Free energy [with elastic moduli C and compliance tensor D = @-‘I:
Y(&,c,a) = f E : c : E + S(a)
Stress-strain and internal variable relations:
a=@:& and q = -d,S
Damage rule for compliance tensor and internal variable:

Kuhn-Tucker damagelno damage conditions:

M
$k(a, 9)< 0, ?’k 2 0 and 1 4k(a, 4)Y t = 0
k=l

Consistency condition:
M

k=l
1 dk(0, 4)Y k =0

The structure of this expression is identical to that arising in the classical model of associative
plasticity, although its physical interpretation is quite different. In (lo), the ‘added’ strain
contribution arises from the loss of stiffness in the material whereas the analogous term in
plasticity theory is due to the change in permanent (plastic) strains at constant stiffness.
For completeness, the rate-independent constitutive equations are summarized in Table I. As
these equations stand, however, the resulting model will yield an ill-posed initial boundary value
problem if the material exhibits softening. A brief discussion of the issues involved is given later in
Section 2.2.

2.1. Rate effects

For materials that display rate effects, the model summarized in Table I can be easily extended
to incorporate first-order rate dependencies. For single surface models (i.e. M = l), the extension
is patterned after the Prager-Perzyna model of viscoplasticity, and amounts to replacing the
Kuhn-Tucker optimality loading/unloading conditions by the consitutive equations:

where ( . ) = +[(.) + 1. I] is known as the Macauley bracket and rj is a scalar governing the
viscosity of the relaxation process. For hardening materials, it is well-known that as q -+ 0 the
rate-independent solution is recovered. Whether the same property holds for softening materials
is an open question. The widespread use of the viscoplastic regularization as a means of
circumventing ill-posedness of the rate-independent problem rests on the belief that this is not the
case.

Remark 2.1. Difficulties arise in the generalization of (11) to the case of an elastic domain
defined by M > 1 intersecting surfaces, unless further restrictions are placed on the gradients d,,#k
of the constraints; see Simo et a1.” for additional details. In particular, the conditions set forth in
this latter reference show that the inviscid limit is well-behaved (for hardening materials) if the
3616 S. GOVINDJEE. G. J. KAY AND J. C. SIMO

gradients are stiffness-orthogonal, i.e. if

for j # k. This observation will prove useful in the rate-dependent extension of the model
described in Section 4.

2.2. Dissipation and characteristic lengths

The equations in Table I are known to produce ill-posed initial boundary value problems
when the elastic domain defined in equation (2) contracts with the evolution of the internal
variable q. The severe difficulties associated with lack of well-posedness are quite different from
those arising from global lack of uniqueness of the initial boundary value problem (see DeBorst”
for an illustration of these latter difficulties). Read and Hegemier,’ for instance, have advanced
a physical mechanism that correlates ‘ill-posedness’ of conventional softening models with
attempts to model a structural response as a continuum response, thus suggesting the introduc-
tion of a length scale in the problem. From a purely mechanical standpoint, a source of difficulty
in the model in Table I lies in the expression for the dissipation, which is intrinsically defined per
unit volume and given by equation (5). Typically, however, for brittle materials the dissipation
mechanism is the generation of cracks. Such a mechanism dissipates energy on a per unit area
basis. Hence the model derived above and the physical mechanism are incompatible in their basic
length scales.
As pointed out in the introduction, a number of approaches can be found in the literature that
introduce a length scale into the model in an attempt to resolve the inconsistency between the
model and the physical mechanism. A review of the main methods currently employed can be
found in Sluys;* these include higher-order gradient theories, rate effect theories, Cosserat
theories and direct characteristic length theories. The inclusion of rate effects as described above
provides for one rather simple method of introducing an implicit length scale into the formula-
tion. This method, however, is lacking in several respects. One, in the limit of the rate-independent
model, the equations may once again produce ill-posed initial boundary value problems. Two, the
viscosity parameter governs both the rate characteristics of the material as well as the character-
istic length of the material see Sluys8 These two phenomena should, however, be governed by
separate parameters as they are separate physical mechanisms not intimately associated with
each other in general. Note also that in an algorithmic setting there may be time step restrictions
on otherwise ‘unconditionally stable’ integration algorithms in order to maintain the ‘well-
posedness’ of the problem.23 Finally, no conclusive proof currently exists on the effectiveness of
this approach other than experience gained in numerical simulations. Direct characteristic length
methods, on the other hand, require a reasonable understanding of the physical cause of the
material’s degradation; in the case of the brittle materials, this information is often available.
An example of a direct characteristic length theory is that of Oliver.16 In this theory,
a dissipation constraint is used to enforce a desired amount of dissipation in a volume of material
with a single crack growing through it. The appeal of this construction lies in its theoretical and
computational simplicity over, perhaps, more correct but also more complex approaches such as
the non-local damage theory of Pijaudier-Cabot and B a ~ a n or t ~ Cosserat
~ theories as in-
l ~ basic idea is that, if the fracture energy per unit area crack advance Gf is a known
D e B 0 r ~ t . The
material property, then the energy dissipated per unit volume gr in the continuum theory is
related through the expression

9f = G f P
 BRITTLE DAMAGE IN CONCRETE 3617

where I* is known as the characteristic length and is given as the reciprocal of the directional
derivative of the crack indicator function in the direction of the crack plane normal; see Oliver16
for more details. Equation (13) provides a practical restriction on the softening properties of the
model. This is seen to be the case since
r
gf= J 9dT
r
where the path of integration r in (14) is from a state of no damage to a state of total failure.
Equation (14) involves the softening properties of the model through the dissipation rate 9 and
hence together with (13) provides a restriction on the softening parameters of the material. Note
that since (14) is a path-dependent integral the path chosen must correspond to that used to
experimentally determine Gf.An application of this method is shown in Section 4.
A completely different approach based on the mathematical analysis of the discontinuous
solutions that arise in the presence of softening is described in Simo et al." The 'regularized'
softening model in this latter reference strictly remains rate independent and the initial boundary
value problem becomes well-posed, leading to numerical solutions exhibiting no mesh-depend-
ence whatsoever. The notion of a characteristic length plays no role in this approach at the
continuum level. The results of this analysis provide a partial justification for the preceding
methodology.

3. COMPUTATIONAL ISSUES
Because of the non-linear nature of the equations in Table I, numerical methods are usually
employed when solving initial boundary value problems. The discussion that follows is restricted
to their use in the standard strain driven framework that exists in most finite element
codes- whether globally performing load or displacement controlled calculations. This frame-
work assumes that the strain E", the stiffness @", and the internal variables a" are all known at time
t". In addition, it is assumed that the strain E " + ~ at time t " + ' > t" is known and that the stress and
values of the stiffness and internal variable at time t " + l are desired. In what follows, superscript
n's and n + 1's will always denote quantities evaluated at times t" and t"", respectively.

3.1. Return mapping algorithm

The method of solving for the desired unknowns parallels very closely the so-called return
mapping algorithms of classical plasticity (see the comprehensive review in Simo3'). The only
substantial difference lies in the details of the consistency parameter calculation when softening is
involved, as will be seen in Section 4. To begin, the flow rules for the compliance tensor and the
internal variables are integrated using the implicit backward Euler rule. This results in the
following two expressions:

k=l

where Ayk 2 0 is a discrete consistency parameter. Note that equation (I 5 ) can be inverted using
the Sherman-Morrison-Woodbury formula (see for example Reference 25, Section 8.3) to give
the stiffness tensor at time t"".
3618 S. GOVINDJEE. G. J. K A Y AND J. C. SIMO

Equation (15) can now be combined with the stress-strain relation (1) at time t"" to give an
implicit expression for the current value of the stress tensor as
r M 1

The expression for the current value of the internal stress is simply

q"+' = - &S(oI,+') (18)

We now require the Kuhn-Tucker damage/no-damage conditions (9) to be satisfied at time t"+
i.e. for ( k = 1,. . , , M ) :
M
4;" < 0, AYk b 0 and 4:'' AYk = 0 (19)
k= 1

Thus at the end of every time step in a calculation the constraints on the admissible stress states
are always satisfied. Further examination of equations (15)-( 19) results in the observation that,
just as in the time continuous case, if 4;' < 0, then Ayk = 0 and no damage evolution takes place
with respect to the kth damage surface. Algorithmically, this motivates the notion of a trial
'elastic' predictor state defined by Ay, = 0 for all k. The trial stress is given by
e".+l - cn:E"+l
trial - (20)
and the trial internal stress is given by

4::; = 4" (21)

A predictor-corrector-type algorithm can now be defined as:
1. For each k check whether &(fit:;, d 0. If so, then assume the surface is inactive and
Set AYk = 0.
2. If &(G:;~;, > 0, then assume the surface is active and that Ay, > 0.
3. If all the surfaces are inactive, then
1 = G ntrial
+ 1
9 q"+ = qYrTa:, a n + ' = a", rpf' = D" (22)
fully defining the state of the system at time t"+'.
4. If any of the surfaces are active, then use the consistency conditions that 4;'' = 0 for
k E {active surfaces} and equations (15)-(18) to simultaneously solve for the unknowns AYk,
e n + ' and
, qntl, and then update the remaining state variables and EDn+'.

Remark 3.1. When the orthogonality condition (12) is not met, a fifth step must be added
to algorithm to insure that the proper set of active surfaces was chosen. For further details see
Simo et a/.''

Remark 3.2. In general, these equations form an intractable highly non-linear system of
equations that must be simultaneously solved for the unknowns; however, for a certain class of
popular damage surfaces this non-linearity is tractable from the standpoint of practical calcu-
lations. Consider the case of a single active damage surface of the form
4 b 4 ) =.fW - Of + 4 = 0 (23)
 BRITTLE DAMAGE IN CONCRETE 3619

where cf is a scalar constant. In this case, equations (17)-( 18) reduce to

By substituting equation (15) into (23), one obtains an equation in terms of a"+' and Ay. This
equation together with equation (24) forms a complete set of equations that may solved for the
unknowns o n + and A?. In the case where dJis a constant, the system of equations reduces even
further into a single (in general non-linear) equation for the scalar Ay. These arguments also hold
true for the case of multiple damage surfaces but yield M equations in Ayk.

3.2. Rate-dependent model

The development of the algorithm for the rate-dependent case exactly mimics the rate-
independent case. The viscoplastic flow rules are first integrated with a backward Euler difference
scheme to give

where At = t"+' - t". The stress at time t"" can now be written as

The internal stress expression (18) with anil now given by equation (27) also holds.
In the rate-dependent problem there are no Kuhn-Tucker conditions, however, an algorithm
similar to the rate-independent algorithm can be used:
For each k check whether 4k(a[i:, & ):; d 0. If so, then assume the surface is inactive and set
<CP;+')At/q = 0.
If (bk(orrL:,q:2a:) > 0, then assume the surface is active and that (&+') At/q > 0.
Same as in the previous algorithm.
If any of the surfaces are active, use the constitutive relations (1 1) for 4;" (k E {active
surfaces}) and equations (26)-(28) and (18) to simultaneously solve for (4;' ' ) At/q, a n +and
'
q"' ',and then update the remaining state variables.

Renzark 3.3. For the general case, the main difference between the rate-dependent and rate-
independent case is that Ayk has been replaced by (&+')At/?. The other difference lies in the fact
that 4;" # 0 for k E {active surfaces} in the rate-dependent case. However, since 6:'' > 0 for
k E (active surfaces}, one can write the analogue to the rate-independent equation 4;' = 0 as

Note that in equation (29) the quantity (4;' I ) At/q is considered as one variable to be solved for;
therefore, in the limit as r] + 0 the rate-independent case is recovered (assuming the orthogonality
condition (12) is satisfied). This observation makes it possible to use the same code for both cases
by always using (29) in the iteration process.
3620 S. GOVINDJEE, G. J. KAY AND J. C . SIMO

3.3 Characteristic length

In conjunction with the algorithms outlined above, an algorithm for calculating the character-
istic length is needed. The discussion that follows pertains to the method proposed by Oliver l 6
and is restricted to the context of isoparameteric Co finite elements. Recall, first, Oliver's
algorithm for calculating the characteristic length, I*:

where N A ( are the element shape functions, Nnodes

a ) is the number of element nodes, zA are the
nodal values of the crack indicator function 7, and n(x) is the normal to the smeared crack plane at
a point x. According to Oliver's definition, the crack indicator function should take the value of
1 in front of the crack field and 0 behind the crack field. Front and back are determined by the
following relation:
1 (xA - x,).n(x) > 0
0 (xA - x,).n(x) < 0

where x represents the point in the body at which one wishes to evaluate the characteristic length,
xA are the nodal co-ordinates, and x, represents the centre of the element in which x lies. In the
context of the 2-D problems examined by Oliver, equation (31) is adequate. However, for 3-D
problems using the common 8-node Co brick element equation (31) leads to a discontinuous
characteristic length function (30).This can best be seen by considering an example: Consider the
8-node Co brick in the isoparametric domain where points are denoted by 4 = (tl, t2,t3).In this
case,

In (32), n = (a& )n/ (1 (8,S)n I( which by construction is of unit norm; hence n can be completely
defined by the spherical angles (0, cp) in the isoparametric co-ordinate frame; 8 is the polar angle
measured from the t1axis and cp is the azimuthal angle measured from the t3 axis. Doing so
allows one to write

I coscp J
Equation (32) is now evaluated for two cases that should produce characteristic length values
arbitrarily close to each other. For case 1, (0, cp) = (x/2,x/4 + E ) where E -4 1. This implies by (31),
that z3 = z4 = z7 = z8 = 1 and z1 = z2 = z 5 = z6 = 0 when using a standard counter-clockwise
node numbering scheme. The second case corresponds to (0, cp) = ((42) + 6, (n/4) + E ) for arbit-
rarily small but positive 6. In this case, z4 = z5 = z7 = za = 1 and z1 = T~ = z3 = 76 = 0.
If (32) is to be continuous with respect to the orientation of the crack field, then cases 1 and
2 should produce the same result within a small tolerance that is a function of 6. For case 1,
equation (32) gives
 BRITTLE DAMAGE IN CONCRETE 362 1

For case 2, equation (32) gives

8,z-n = II(8,S)nIIa(sin8sincp + cosq)(1 + (35)
If one assumes a full quadrature rule, the reciprocal of the characteristic length at the Gauss
points ( & 1/&, k l/& k 1/J?) reduces to

where the sign in (1 f $1 depends upon the Gauss point chosen. Clearly (36) does not in general
tend to (34) in the limit of vanishing 6. Only when t1t3= 0, does (35)reduce to (34).
The difficulty that arises is due to the discontinuous nature of the definition of the nodal values
for the crack indicator function. To circumvent this difficulty, a simple continuous definition for
T~ can be made as follows:

(xA- x,).n(x) - zmin

?A = ,pax - Tmin (37)

where

The redefinition of (31) as (37) removes the continuity problems noted above and yet remains
faithful to the original idea proposed by Oliver.

4. APPLICATION: CONCRETE’
In this section, an application is made of the developments of Sections 2 and 3. The example
problem to be considered is the brittle tensile failure of concrete. Failure in the material is
assumed to initiate when the 1st principal stress in the material exceeds some threshold value. In
the discussion that follows, the principal direction associated with this principal stress is assumed
known and given by the eigenvector n. This vector defines the normal to a ‘smeared crack field
that is locally fixed in the material after initiation. Across this smeared crack field, tensile tractions
are limited by some critical value that decreases exponentially with increasing deformation.
Likewise, shear tractions across the smeared crack field are limited by some critical value that
decreases with increasing deformation. The progressive degradation of strength as employed here
is a generalization of the original smeared crack model of RashidZ7 and represents the more
modern viewpoint that the degradation is progressive; see e.g. DeBorst and Nauta.28

4. I . Damage surfaces
To begin, a set of damage surfaces must be defined. Based on the description above, three
coupled surfaces are postulated-the first to control the tensile tractions across the smeared

‘A more detailed exposition of the derivations in this section may be found in Govindjee et
3622 S. GOVINDJEE, G. J. K A Y AND J. C. SIMO

crack field and the second two to control the shear tractions across the smeared crack field:
41 = S1 -fn + k,q < 0 (40)
4 2 = 1Sz:al -f, + k,q < 0 (41)
43 = IS3:al -fs + ksq < 0 (42)
In the above, S1 = n @n, Sz = i ( n @ m + m @ n), and S3 = i ( n @ I + I @ n); m and 1 denote unit
vectors that together with n form an orthonormal basis; f, is the critical tensile traction that can
be held across the smeared crack field;f, is the critical shear traction that can be held across the
smeared crack field; q is an internal scalar softening stress tending tof,; and k, and k, are coupling
constants that are chosen so that the tensile and shear tractions across the smeared crack field
reach their asymptotic values simultaneously.
Theoretically, one would like to choose k, = 1; however, for computational reasons it is better
to choose k, = (1 - fl,) where fl, is a small number. This gives the material a small amount of
residual tensile strength. The constant k, is similarly set to (fS/f,) (1 - fl,), where fl, represents the
residual shear strength of the material. Ideally, this value should be pressure dependent; however,
for simplicity this added complexity is avoided here. Additionally, the softening stress is assumed
to be given by the following relationship:
q =f,(l - exp[ - Ha])= - d,S (43)
where H is a softening decay constant that is determined by equations (13) and (14). The
motivation for this expression comes from the experimental observation of exponential-type
softening behaviour.

4.2. Flow rules and consistency parameters

Using (40) and (42) in (8) gives the damage flow rules:

Remark 4.2. Flow rule (44) indicates that if the material is initially isotropic then the ortho-
gonality condition (12) is satisfied. This follows since S1:S2 = S1:S3 = S2:S3= 0 and
S1:Co:S2 = S1:Co:S3 = Sz:Co:S3= 0 where Co is the initial undamaged isotropic stiffness of
the material.

Remark 4.2. The consistency parameters may be easily calculated using standard techniques;
see Govindjee et a1.26for further details.

4.3. Dissipation constraint

For brittle materials like concrete, the fracture toughness (energy per unit crack advance) is
usually measured under pure Mode-I conditions. Thus the appropriate expression for the
dissipation (14) is
 BRITTLE DAMAGE IN CONCRETE 3623

where all the shear terms have been assumed zero. Simplifying (46) with the aid of (40) and
substituting into (14) gives
P ,

The path of integration is from a state of no-damage, a = 0, to a state of complete damage a,. The
state of complete damage is defined by q(aJ =fn. However, in the present case there is no such
state for finite a; see equation (43). Therefore, a, is set as the value where q reaches a predefined
fraction of its asymptotic limit. Equally, one could define a, through the relation Ha,= 1,where
for instance a value of i= 4 would correspond to a value of q within 2 per cent of its asymptotic
limit fn. Doing so in (47) implies that

Thus, the softening decay constant is constrained by (13) to be

4.4 Algorithmic issues

The algorithmic approximation of the above equations follows directly from Section 3. In
particular note that Remark 3.2 applies to the postulated damage surfaces. When the algorithm
given in Section 3 is used to advance the solution from time t" to time t"", the primary issue is
Step 4 in the algorithm; how does one solve for the unknowns? Because Remark 3.2 applies to the
proposed model, the problem reduces to a non-linear root-finding problem in several variables.
To make this root-finding problem transparent, first consider the special case where only the
tensile surface is active. The case of multiple active surfaces is considered afterwards.
When only the tensile surface is active, the equation to be solved is

4, = S 1 : f 9 + l- f n + kn4n+l = 0 (50)
With the use of the time integrated flow rules, this equation reduces to a single non-linear
equation with only one unknown, Ayl:

qL(Ayl) = Sl:aLk: - AyISI:@":S1 -f. + k, f. - k,f,exp( - Hk,Ay,)exp( - Ha")= 0 (51)

This equation can be solved using Newton's method or any other suitable root finder.

Remark 4.3. Unlike conventional metal plasticity, equation (51) is not convex in Ay,; rather, it
is concave. This follows by considering the second derivative of +1 with respect to Ay,:

--
a241
- -fnk:Hexp(- Hk,Ay,)exp( - Ha") < 0

Thus a qualitative plot of 41reveals that the starting iterate for a Newton-like method must lie to
the right of the curve's apex. The apex location is given by the condition that a41/aAyl = 0.
A more simple way of selecting a starting value is to choose a value to the right of the desired root
3624 S. GOVINDJEE, G. J. KAY AND J. C. SIMO

by setting the terms giving exp( - Hk,a"+l) to zero in (51). This gives a starting value of

Algorithmically, this is convenient since all the required quantities are used elsewhere in the
algorithm. If the single surface chosen had been 42, the initial starting value for a Newton
iteration method can be obtained in a similar manner. This leads to

An analogous result holds for &.

The case where all three surfaces are active requires the solution of a set of three non-linear
equations in three unknowns:

41(Ayk) = S1 + k,q"+'(Ayk) = O
- A71 S1 :@":S1-fn
I -~Ay2
42(AYk) = 1 S ~ : n ~ ~s2:@":&-fs + k,q"+'(Ayk) = 0 (55)
-fs + k,q"+'(Ay,) = 0
43(Ayk)= (S3:a::/ 1 - Ay, S3:@":S3
This set of equations can be solved, as in the single surface case, using Newton's method. The
initial guess for the iteration scheme can be chosen by applying the same technique as was used in
the single surface case; i.e. by setting exp( - Ha"' ') = 0 and solving for (Ay?, A?!, A$). Note this
produces the same starting value formulae as before, because the procedure of setting the
exponential to zero uncouples the three questions.

Remark 4.4. To use the expressions given above, a knowledge of the material stiffness is
required. As mentioned before, the Sherman-Morrison-Woodbury formula may be used to
invert the compliance expression to give the stiffness:
3
@"+I = @" + 1 rk(@":Sk)@(Cn::Sk)
k=l

where

Remark 4.5. Crack closure is easily incorporated into this scheme by setting the material
stiffness to its undamaged value whenever S1:o < 0. Upon re-encountering a state where
S1:n 3 0 the damaged stiffness tensor can be 'recovered'.

4.5. Algorithmic tangents

When solving initial boundary value problems it is often desirable to have the algorithmic
tangent modulus of the material model. Note that this modulus differs from the continuum
tangent modulus as was pointed out by Simo and Taylor.29 By following the well-established
~ ' tangent modulus for the case of only the tensile surface active
procedure presented in S i m ~ , the
 BRITTLE DAMAGE IN CONCRETE 3625

is given by
(S1:C") @ ( S , :C")
@$' = C" -
S1:C":Sl- k,2fnHexp( - Ha"")
In the three-surface case, the identical procedure leads to the following expression for the tangent

4.6. Rate effects

As noted before, the rate-dependent case can be easily treated once the rate-independent case
has been worked out. In particular, the system of three equations that must be solved in the rate-
dependent case are derived from (55) by: (1) replacing Ayk by ( & + ' ) A t / q and (2) replacing
instances of s k : c " : s k by S k : @ " : s k + q/At. In a Newton iteration method ( & + ' ) A t / q are treated
as the unknowns. The starting values for the iteration scheme are obtained from the rate-
independent formulas via the two substitutions given above. The identical substitutions likewise
convert the rate-independent algorithmic moduli (58) and (59) to the rate dependent case.

Remark 4.6. The regularizing effect of viscosity" is now transparently seen through its effect
on the algorithmic tangent moduli. Consider the single-surface case given in equation (58) which
converts to
(S1:C")@(S,:@")
- (62)
Cg2visc = @"
(Sl :Cn:S1 + 5)
- k, f. Hexp( - Ha"+')

Stability of a Newton-Raphson finite element global solver will depend strongly on how close the
denominator of the second term on the right is to zero. The effect of the viscosity is easily seen to
make the denominator more positive and to lessen any numerical ill-conditioning caused by
a small denominator. Note however that one may be required to decrease the time step At
appreciably to take advantage of this effect-depending, of course, on the material properties.
Thus viscous behaviour is seen to provide an added stabilizing effect on numerical calculations
for sufficiently small time steps-in addition to providing an implicit characteristic length scale
for the model.
3626 S.GOVINDJEE, G. J. KAY A N D 1. C . SlMO

4.7. Example calculations

To close this section, a series of numerical examples are provided to demonstrate the effec-
tiveness of the model. The first example is concerned with plain concrete and examines the effect
of the shear surfaces & and &. The second example is the analysis of a notched plain concrete
beam, the third is the analysis of a lightly reinforced, and the last is the analysis of an overly
reinforced beam, all under three-point bending.
Calculations were made using either the finite element code FEAP from the University of
California at Berkeley or the finite element code NIKE3D from the lawrence Livermore National
Laboratory. All calculations were done with 8 node brick elements with a standard Galerkin
formulation and the characteristic length method described above.

Tapered wedge. In this example a tapered wedge is analysed. The geometry, material properties
and loading for this example are shown in Figure 1; the loading is displacement controlled across
the entire top face of the wedge. The three hidden faces in Figure 1 are symmetry planes. Shown in
Figure 2 are two load-deflection plots for the response of the wedge on its top face. The dashed line
represents the analysis when only surface (no shear surfaces considered) is used in the calculation
and the solid line represents the analysis when all three damage surfaces are employed. The two
curves remain identical up to a displacement of 0.0013 in where they diverge with the solid line
giving a more realistic result. This is due to the 3-D nature of the stress field near the centreline of
the wedge where the smeared cracks are not prefectly parallel to the horizontal plane. The result is
that the imposed loading produces shear stresses across these crack planes. In the case of the solid
line, the damage surfaces 42and $3 become activated, and control and degrade ths transmission of
stress-resulting in the total separation of the block. In the case of the dashed line, these shear
stresses are not controlled and eventually become the dominant mechanism for the transfer of load
across the smeared crack planes-resulting, eventually, in the load-displacement plot having
a non-physical positive slope at large extensions. The control of shear stresses (tractions) across
smeared crack planes is clearly required for problems of this nature.

Notched plain concrere beam. In this example a notched plain concrete beam in three-point
bending is analysed. The geometry, material properties and loading for this example are shown in

, Displacement Driven Face

Figure 1. Tapered wedge shown in one-eighth symmetry with 27 8 node brick elements. Material properties:
E = 4 x lo6 psi, v = 0.21, .f. = 400 psi, = 2000 psi. fi, = 0.01, G , = 1.72 Ib/in, and 7 = 0 psi sec
 BRITTLE DAMAGE IN CONCRETE 3627

350 I I I 1

300

Shear surfaces on -
250

-2
* 200 Shear surfaces off -----
v

-2
3 150

100

50

0-
0 1 2 3 4 5
E-3 Displacement (in)

Figure 2. Load versus displacement response for tapered wedge example. The solid line gives the response when both
tensile and shear tractions are taken into account. The dashed line gives the response when only tensile tractions are
accounted for

Figure 3, where the rear and left faces of the geometry shown are symmetry planes. The material
properties given are the mean values reported from experiments on this beam geometry; see
Malvar and Warren31 and Malvar and F ~ u r n e yThe . ~ ~critical shear value was obtained from
a simple Mohr’s circle argument as half of the reported compressive strength, fel = 4206 psi; the
tensile retention value was chosen for numerical stability close to machine zero and the shear
retention value was arbitrarily set to a moderately small number. In the range chosen the actual
value of the shear retention parameter does not affect the beam’s overall response appreciably.
The calculation itself was displacement controlled. As the deflection is increased, the calculation
results predict a narrow smeared crack originating at the notch root and propagating up the
centre of the span. The resulting load-deflection curve under the load point is plotted in Figure 4;
the solid line is the calculation and the triangles represent the mean of the data from several
experimental runs. The overall agreement is seen to be reasonably good and refinement of the
mesh and the properties would likely improve the agreement. Note that agreement beyond
a deflection of 0018 in. deteriorates rapidly. This is due to the fact that the response of the beam is
3628 S.GOVINDJEE. G. I. KAY AND J. C.SlMO

Figure 3. Finite element discretization of the notched plain concrete beam. The beam is shown in quarter symmetry.
Material properties reported: E = 3.15 x 10' psi, fm = 449.6 psi and GI= 0.436 Ib/in. Material properties estimated:
v = 0.2, fr = 2103 psi, 3/, = 003, q = 0 psi sec

Figure 4. Load versus displacement curves under the load point for the notched plain concrete beam. The solid line gives
the calculation and the triangles give the experimental data Malvar and Fourney3'

effectively governed by only a few 8 node brick elements in the top fibres of the beam; these
elements are known to respond poorly in pure bending.

Moderately reinforced concrete beam. In this example, a moderately reinforced beam under
three-point displacement driven bending is analysed. The geometry, material properties and
loading for this example are shown in Figure 5, where the left and the rear faces are symmetry
planes. This problem corresponds to an experiment conduct by Burns Sei~s~~-beam J4.Note,
however, that for simplicity the reinforcing pattern for this beam was altered by removing the
stirrups as in Kwak and F i l i p p ~ u The
. ~ ~ material properties are those reported by Burns and
S e i ~ s where
, ~ ~ indicated. The tensile strength of the concrete was estimated using the relation
f. = 5 fi wheref,' was reported to be 4820 psi; the shear strength was estimated to be five times
the tensile strength. The Young's modulus of the concrete was taken from Kwak and F i l i p p o ~ . ~ '
The remaining concrete properties were arbitrarily set within the spread of data for concrete. For
the rebar, the hardening modulus was estimated from the experimental data.
 BRIlTLE DAMAGE IN CONCRETE 3629

Full Cross Section 20 in

t I 9
5 in
Figure 5. Finite element discretization of the moderately reinforced concrete beam. The beam is shown in quarter
symmetry with the inset showing the reinforcing pattern. Material properties reported: (concrete) E = 3.5 x lo6 psi, (rebar)
E = 29.5 x lo6 psi and u, = 44.9 x lo3 psi. Material properties estimated: (concrete) Y = 0.2, f ,= 350 psi, 1; = 1750 psi,
8, = 005, Gf= 1.0 Ib/in, q = 0 psi sec, (rebar) v = 0.3, and H = 27.0 x lo3psi

Displacement (in)

Figure 6. Load versus displacement curves under the load point for the moderately reinforced concrete beam. The solid
line gives the calculation and the triangles give the experimental data Burns and S e i s P

In the calculation, as the loading is applied to the beam, a diffuse crack field propagates up
from the lower fibres of the beam with a dominant crack in the Centre of the span. The cracks
furthest from the centreline of the beam also curve in towards the loading. This description
corresponds well to the observed crack pattern. Figure 6 shows the resulting load-deflection
curve underneath the load point; the solid line is the calculation and the triangles are the
data. Up to a load of about 10 kips the beam behaves elastically. At this point the concrete begins
to crack and load is transferred into the reinforcing rods (which were discretely modelled with
3630 S. GOVINDJEE, G.J. K A Y AND J. C. SIMO

12 in
41-1

Full Cross

t 1-1
7 in
Figure 7. Finite element discretization of the heavily reinforced concrete beam. The beam is shown in quarter symmetry
with the inset showing the reinforcing pattern. Material properties of concrete (all estimated): E = 3.4 x lo6 psi, v = 0.167,
f. = 350 psi, & = 1750 psi, PI= 0.05, GI = 4.0 Ib/in, and 9 = 0 psi sec. Material properties of lower four rebar ( # 9’s):
(reported) E = 31.6 x lo6 psi, (estimated) I’ = 0.3, (reported) rry = 80.5 x lo3 psi, and (estimated) H = 31.6 x lo3 psi. Ma-
terial properties of upper two rebar ( # 4‘s): (reported) E = 29.2 x lo6 psi, (estimated) \’ = 0.3 (reported) uY= 50.1 x lo3 psi,
and (estimated) H = 29.2 x lo3 psi. Material properties of stirrups rebar ( # 2’s): (reported) E = 273 x 10” psi, (estimated)
1’ = 0.3, (reported) op= 47.2 x lo3 psi, and (estimated) H = 27.5 x lo3 psi

elasto-plastic beam elements). When the reinforcing rods reach their yield point, the beam loses
almost all of its tangent stiffness, as indicated by both the data and the calculation results in the
figure. Note that in this calculation the effects of the shear damage surfaces are not appreciable
until the reinforcing rods begin to yield and, that after yielding, the slope of the load-deflection
curve is strongly governed by both the hardening modulus of the reinforcing rods and the shear
degradation properties of the concrete. As noted in the original report and later by Kwak and
’ the level of loading shown, the effect of reinforcing rod slippage is minimal. It is
F i l i p p ~ u , ~for
further remarked that even though the structure has not softened globally, locally the concrete
softens appreciably.

Heavily reinforced concrete beam. In this example, an heavily reinforced beam (1.53 per cent
volume of rebar) under three-point bending is analysed. The geometry, material properties and
loading for this example are shown in Figure 7, where the left and rear faces are symmetry planes;
(the stirrup spacing is 8-25 in). The beam corresponds to an experiment reported in Bresler and
S~ordelis~~-beam A l . The tensile strength of the concrete was estimated to be roughly one-tenth
of the reported compressive strength of 3.49 ksi and the shear strength of the concrete was
estimated to be roughly half of this value. The Young’s modulus and Poisson ratio were taken
from Kwak and Filippou3’ and the remaining concrete properties were estimated within the
acceptable range for concrete. The hardening moduli for the reinforcing rods were estimated from
the experimental data-note however that there is no rebar yield in this calculation.
The heavy reinforcing pattern makes the inclusion of bar slip in this example important. For
the calculation, the bar slip behaviour was modelled with a simplistic stick-slip law which was
implemented in the 1-D slideline logic in NIKE3D by Maker and L a ~ r s e n A. ~slip ~ value of
150 psi for a wetted area load was chosen by trial and error.
In the calculation (displacement controlled), a very diffuse crack field propagates from the
lower fibres of the beam towards the top. This crack field eventually forms a series of prominent
 BRITTLE DAMAGE IN CONCRETE 363 1

Figure 8. Smeared crack pattern projected onto the side of the beam shown in half symmetry. The intensity of the shading
indicates the intensity of crack opening (damage)

100 -

80 -

60 -

40 -

n 0.1 0.2 0.3 0.4 0.5

Displacement (in)

Figure 9. Dad versus displacement curves under the load point for the heavily reinforced concrete beam. The solid line
gives the calculation and the triangles give the experimental data Bresler and S c o r d e l i ~ ~ ~

cracks spaced out along the length of the beam and the overall crack field tends to point towards
the load point. Figure 8 shows the smeared crack field projected onto the outerface of the beam
half way through the loading; the intensity of the crack opening is indicated by the intensity of the
shading. Note, in particular, that horizontal cracks underneath the load point are predicted by
the model. The overall crack field qualitatively corresponds well to the experimental observa-
tions. Figure 9 shows the resulting load-deflection curve underneath the load point; the solid line
represents the calculation and the triangles represent the data. The agreement is seen to be good.
If slippage of the rebar had not been included in the calculation, the predicted response of the
beam would have diverged sharply from the actual data at a deflection of roughly 0.3 in.

5. CLOSURE
An alternative framework for continuum damage mechanics has been presented and exploited in
the design and numerical implementation of a constitutive model for the tensile failure of plain
concrete. For a given failure surface, often specified via experimental testing, the model automati-
cally furnishes a generally fully anisotropic damage evolution law under the single assumption of
3632 S. GOVINDJEE, G.J. K A Y AND J. C. SIMO

maximum (damage) dissipation. Once the failure surface is specified, no adjustable parameters are
present in the model. A detailed algorithmic treatment has been presented which addresses issues
related to time integration as well as the inclusion of regularization approaches leading to
mesh-independent calculations. Although no spurious mesh-dependencies arise, the global solu-
tion to the initial boundary value problem generally retains the characteristic lack of uniqueness
inherent to this class of problems. In spite of recent advances, the design of robust procedures for
dealing with this issue remains unsettled.
The proposed three failure surface model for plain concrete exhibits reasonable agreement with
limited experimental results given the gross nature of the assumptions, thus lending support to
the hypothesis of maximum dissipation that dictates the form of the damage rule. It is emphasized
that the proposed damage model for the brittle failure of concrete contains a minimal set of
material constants, which can be determined from standard tests.
ACKNOWLEDGEMENTS

The authors would like to thank Dr. G. L. Goudreau of the Lawrence Livermore National
Laboratory for sponsoring this work, Prof. F. Filippou of the University of California at Berkeley
for many helpful discussions on the behaviour of concrete, Drs. B. Maker and R. Whirley of the
Lawrence Livermore National Laboratory for their help in implementing versions of the model in
the LLNL finite element codes DYNA3D and NIKE3D, and Prof. R. L. Taylor of the University
of California at Berkeley for many helpful discussions and for graciously providing a version of
his finite element code FEAP for the development of the model.

REFERENCES
1. L. M. Kachanov, ‘On creep rupture time’, IZV Akad, Nauk SSSR, Otd. Techn. Nauk, 8, 26-31 (1958).
2. I. L. Lemaitre, A Course on Damage Mechanics, Springer, Berlin, 1992.
3. M. Ortiz, ‘A constitutive theory for the inelastic behavior of concrete’, Mech. Mat., 4, 67-93 (1985).
4. I. C. Simo and J. W. Ju, ‘Strain- and stress-based continuum damage models-I. Formulation’, Int. J . Solids Struct.,
23, 821-840 (1987).
5. R. Hill, The Mathematical Theory of Plasticity, Oxford University Press, Oxford, 1983.
6. S. Govindjee and J. C. Simo, ‘Transition from micro-mechanics to computationally efficient phenomenology: Carbon
black filled rubbers incorporating Mullins’ effect’, J . Mech. Phys. Solids, 40, 213-233 (1992).
I. Z. P. Baiant and T. Eklytschko, ‘Wave propagation in a strain softening bar; Exact solution’, J. Eng. Merh., 111,
381-389 (1985).
8. L. J. Sluys, ‘Wave propagation, localization and dispersion in softening solids’, Doctoral Thesis, Delft University of
Technology, 1992.
9. H. E. Read and G. A. Hegemier, ‘Strain softening of rock, soil and concrete-A review article’, Mech. Mat., 3,271-294
( 1984).
10. S . T. Pietruszczak and 2. Mroz, ‘Finite element analysis of deformation of strain-softening materials’, Int. j . numer.
methods eng., 17, 327-334 (1981).
11. Z. P. Baiant and B. H. Oh, ‘Crack band theory for fracture of concrete’, Muter. Constr., 16(93), 155-177 (1983).
12. A. Needleman, ‘Material rate dependence and mesh sensitivity in localization problems’, Comput. Methods Appl.
Mech. Eng., 67, 68-85 (1987).
13. H. L. Schreyer, ‘Analytical solutions for nonlinear strain-gradient softening and localization’, J. Appl. Mech., 57,
522-528 (1990)’
14. R. DeBorst, ‘Simulation of strain localization: A reappraisal of the Cosserat continuum’, Eng. Comp., 8, 317-332
(1 99 I).
15. J. C. Simo, J. Oliver and F. Armero, ‘An analysis of strong discontinuities induced by strain-softening in rate-
independent inelastic solids’, Comp. Mech., 12, 277-296 (1993).
16. J. Oliver, ‘A consistent characteristic length of smeared cracking models’, Int. j . numer. methods. eng., 28, 461-474
(1989).
17. J. C. Simo, ‘Some aspects of continuum damage mechanics and strain softening’, Contract report to Aptek Inc. for
contract No. DNA-001-85-C-0264, 1988.
18. C. Truesdell and W.Noll, ‘The non-linear field theories of mechanics’, in S. Fliigge (ed.),Handbuch der Physik, Band
11l/3, Springer, Berlin, 1965.
19. D. G. Luenberger, Linear and Nonlinear Programming, Addison-Wesley, 1984.
 BRITTLE DAMAGE IN CONCRETE 3633

20. W. T. Koiter, ‘Stress-strain relations, uniqueness and variational theorems for elastic-plastic materials with a singular
yield surface’, Quart. Appl. Math., 11, 350-354 (1953).
21. J. C. Simo, J. G. Kennedy and S. Govindjee, “on-smooth multisurface plasticity and viscoplasticity. Load-
ing/unloading conditions and numerical algorithms’, Inr. j . numer. methods eng., 26, 2161-2185 (1988).
22. R. DeBorst, ‘Computation of post-bifurcation and post-failure behavior of strain-softening solids’, Comp. Struct., 25,
21 1-224 (1987).
23. J. C. Simo, ‘Strain softening and dissipation: A unification of approaches’, in Z. P. Baiant and J. Mazars, (eds.),
Cracking and Dumage, Elsevier Applied Science, New York 1988, pp. 440-461.
24. G . Pijaudier-Cabot and Z. P. Baiant, ‘Nonlocal damage theory’, J . Eng. Mech., 113, 1512-1533 (1987).
25. J. E. Dennis, and R. B. Schnabel, Numerical Methods f o r Unconstrained Optimization and Nonlinear Equations,
Prentice-Hall, Englewood Cliffs, N.J., 1983.
26. S. Govindjee, G. J. Kay and J. C. Simo, ‘Anisotropic modeling and numerical simulation of brittle damage in
concrete’, Report N o . UCB/SEMM-94/18, University of California at Berkeley, Department of Civil Engineering,
1994.
27. Y. R. Rashid, ‘Ultimate strength analysis of prestressed concrete pressure vessels’, Nucl. Eng. Des., 7 , 334-344 (1968).
28. R. DeBorst and P. Nauta, “on-orthogonal cracks in a smeared finite element model’, Eng. Comp., 2,35-46 (1985).
29. J. C. Simo and R. L. Taylor, ‘Consistent tangent operators for rate independent elastoplasticity’, Comput. Methods
Appl. Mech. Eng., 48, 101-1 18 (1986).
30. J. C. Simo, ‘Topics in the numerical analysis and simulation of classical plasticity’. in P.G.
Ciarlet and J. L. Lions (eds.),
Handhook fbr Numerical Analysis, Vol. 111, North-Holland, Amsterdam (in press).
31. L. J. Malvar and G. E. Warren, ‘Fracture energy of three-point bend tests on single-edge-notched beams’, Exper.
Mechanics, 28,266-272 (1988).
32. L. J. Malver and M. E. Fourney, ‘A three dimensional application of the smeared crack approach., J . Eng. Frac.
Mech., 35, 251-260 (1990).
33. N. H. Bruns and C. P. Seiss, ‘Load-deformation characteristics of beam-column connections in reinforced concrete’,
Civil Eng. Studies, SRS No. 234 (1962).
34. 9. Bresler and A. C. Scordelis, ‘Shear strength of reinforced concrete beams’, J . Am. Conc. Inst., 60, 51-72 (1963).
35. H. G . Kwak, and F. C. Filippou, ‘Finite element analysis of reinforced concrete structures under monotonic loads’,
Report No. LICEISEM M-90/14, University of California at Berkeley, Department of Civil Engineering, 1990.
36. 9. N. Maker and T. A. Laursen, ‘A finite element formulation for rod/continuum interactions: the one dimensional
slideline’, Int. j. numer. methods eny., 37, 1-18 (1993).