Comput.m & Structures Vol. 32, No. 314, pp.

563-590, 1989
Printed in Great Britain.

NONLINEAR

ANALYSIS

OF CONCRETE

Co457949/89 $3.00 + 0.00

1989 Maxwell Pqamon Macmillan plc

STRUCTURESt

KLAUS-J~RGENBATH&$@JAN WALCWK, llfi ALAN WELC.!E@~~

and NAGIN MrsrRvl/tt
~Ma~chu~t~ Institute of Technology, Cambridge, MA 02139, U.S.A.
JADINA R & D, Inc., Watertown, MA 02172, U.S.A.
TTaylor Woodrow Services, Ltd, Middlesex, U.K.
Abstract-Solution
capabilities for two- and three-dimensional nonlinear finite element analysis of
concrete structures are presented. The concrete material is modeled including the triaxial nonlinear
stress-strain behavior, tensile cracking, axon-stiff~~g,
compression crushing and str~n-sof~ng.
The
results of various sample analyses are given to demonstrate the constitutive model and the solution
strategies. These analyses include the response predictions of a test beam and of small scale models of
two reactor vessels.

1. INTRODUCTtON
During recent years, interest in nonlinear analysis of
concrete structures has increased steadily, because of
the wide use of plain, reinforced and prestressed
concrete as a structural material, and because of the
development of relatively powerful finite element
procedures [l]. If a realistic nonlinear analysis of a
concrete structure can be carried out, the safety of the
structure is increased and the cost can frequently be
reduced.
Concrete exhibits a complex structural response
with various important nonlinearities; namely, a nonlinear stressstrain
behavior, tensile cracking and
compression crushing material failures and temperature-dependent creep strains [2-91, All these concrete
nonlinearities depend strongly on the triaxial state of
stress, and in addition the nonlinearities introduced
by reinforcing and prestressing steel should in general
be taken into account [l&13].
There are a number of factors that have prevented
the wider acceptability of nonlinear finite element
analysis procedures in the analysis of concrete structures. A first important consideration is that the
~nstitutive
properties of concrete have not as yet
been identified completely, and there is still no generally accepted material law available to model concrete
behavior. A second important factor is that nonlinear
finite element analysis of concrete structures can be
very costly and may require considerable user expertise. The high cost of nonlinear analysis of concrete
structures is largely due to the difficulties encountered
in the stability and accuracy of the solutions. These
difficulties, however, are a direct consequence of the
specitic numerical implemen~tion
of the concrete
nonlinearities.

t Original figures were generated with a colour producing

terminal and submitted in colour.
BCo-author of Part I.
tt Co-author of Part II.

The objective in this paper is to present the threedimensional concrete model available in ADINA and
show some applications. The concrete model was in
many respects already proposed in [2]. However, our
recent improvements on the model have been very
significant and promise to make the model much
more attractive for practical nonlinear concrete
analyses.
~thou~
we could have chosen to present only our
recent improvements on the concrete model, we
believe that it is of more interest to have the complete
model description in this paper with applications.
The basic aim in the development of the model was
to implement in ADINA a model that with the
present constitutive descriptions, numerical methods
and computing equipment available would satisfy the
following two criteria. Firstly, the model should be as
simple as possible, but reproduce the important nonlinear and strength characteristics consistent with
experimental results. Secondly, the model should be
theoretically sound and numerically stable, so that
reliable analysis results are obtained.
The material model is a hypoelastic model based
on a uniaxial stress-strain relation that is generalized
to take biaxial and triaxial stress conditions into
account. Tensile cracking and compression crushing
conditions are identified using failure surfaces, and
strain-softening effects are included in the compression and the tensile regions. A particularly valuable
property of the model is that it is defined by a number
of input parameters that provide considerable
flexibility in its use (and indeed the model may also
be employed to model sand or rock structures).
In the following sections we first present the concrete model in detail. We then give some small-size
applications, merely to demonstrate the details of the
model behavior when it is subjected to various
stress-strain paths. Finally, in Part II of tbe paper we
present the results of the analyses of the Sandia
pressure vessel, an alkali-silica reacted beam and a
PWR prestressed concrete containment. The solution
563

KLAUS-JBRGEN BATHE et al.

564

COMPRESSIVE
-0032

STRAIN. in Iln
-0.001

KUPFER
ET AL
(EXPERIMCNTA

PRESENT

2.6

S
------.

----

-----

Fig. I. Concrete model uniaxial stress-strain

results presented in this section are particularly valuable because they can be compared with and interpreted against laboratory experimental results.
An important issue in nonlinear analysis is the
visualization of the obtained response predictions. In
this paper we use our standard post-processing program ADINA-PLOT to present all of our solution
results.
2. PART I: THE CONCRETE

MATERIAL

MODEL

The model implemented employs three basic features to describe the material behavior, namely (i) a
nonlinear stress-strain
relation including strainsoftening to allow for the weakening of the material
in compression, (ii) failure envelopes that define
cracking in tension and crushing in compression, and
(iii) a strategy to model the post-cracking and crushing behavior of the material. In the solution, the
material can be subjected to cyclic loading conditions,
i.e. the numerical solution allows for unloading and
reloading including deactivation of tensile failures.
However, as will become apparent, the cyclic loading
conditions are only modeled realistically in situations
of essentially proportional loading.
In the following, the material model is described
for infinitesimal displacement conditions using the
engineering stresses triiand engineering strains eii. In
order to analyze problems with large rotation conditions, the total Lagrangian stress and strain variables
must be substituted for the engineering variables [ 11,
and then the model is directly applicable.7

We

use the notation of Refs [I,21 in this paper.

2.1. Stress-strain

I.0

law.

relations

The general multiaxial stress-strain relations are

derived from a uniaxial stress-strain relation 6 versus
e. In this section, we describe the uniaxial and
multiaxial stress-strain relations employed in the
model and concentrate on the behavior prior to
tensile cracking or compression crushing.
In the following discussion, all uniaxial parameters
are identified by a tilde () placed over them (i.e. all
parameters that have been obtained from Fig. I carry
a tilde).
2.1.1. Uniaxial conditions. A typical uniaxial stress,
5, to uniaxial strain, 6, relation (assuming continuous loading of the material) is shown in Fig. 1 (61.
This stress-strain relation shows that there are basically four strain phases; namely, corresponding to
.C,>e>O, &,a~?>~,,, O>e>& and &ce>&,
where c?,is the strain corresponding to the uniaxial
cut-off tensile stress d,, ?,,, is equal to t; . F,, where 5
is an input parameter, CCis the strain corresponding
to the minimum (crushing) stress, d,, that can be
reached, and E?,,is the ultimate compressive strain. If
C,3 % > 0, i.e. the material is in tension, the
stress-strain relation is linear until tensile failure at
the stress d,, and a constant Youngs modulus & is
employed,
5 = &e^
(1)

For 5 < 0, we assume the following relation,

a/a, =

b%/Jw~l~c)
1 + A@/&) + B(P/&)2 + C(9/&,))

(3)

565

Nonlinear analysis of concrete structures

and hence,

mined,
$[l

(E =

!f<f-

- B(e/S$ - 2C(e/&)3]
(4)

[1 + A (6?/frc)+ B(t?/Q + c(Z/ZJr]*

where
A = &l/J% + (P3 - 2P2)&/$

- (2P3 - 3P2 + 1)

(P2 - 2P + l)P
B=(2$/gJ-3)-2A,

C=(2-&/$)+A

and the strength parameters &, d,, Cc,& = 6c/t?c,d,,

&, p = &,I& and & = a./$,, are obtained from uniaxial tests. It may be difficult to determine somewhat
accurate values of d, and CM;,
in which case reasonable
values should be assumed by the analyst.
The stress-strain
relation in eqn (3) assumes
monotonic loading conditions. For unloading conditions and loading back to the stress state from which
unloading occurred, the initial Youngs modulus i$ is
used. For strain states beyond & in compression, it is
usually assumed that the material cannot resist any
load and the stresses are set to zero. However, an
option for a gradual stress unloading beyond 5 is
also included.

- (1 -v)

(6)

where f,, is the maximum value of the loading

function that has been reached during the complete
solution. In unloading, the initial Youngs modulus,
&, is used to form the incremental stress-strain
matrix, both for stiffness and stress calculations.
To obtain the stress-strain relations in loading
conditions, the principal stresses are calculated and
for each principal stress direction a uniaxial tangent
Youngs modulus, &, corresponding to the strain in
the principal stress direction, epi, is evaluated using
eqns (2) and (4). When using eqn (4), the current
strain epiis employed and to account for multiaxial
stress conditions the material variables g,, 6;, & and
I& are replaced by the variables ai, a:, g: and Z:
defined in eqns (13) and (14).
Let lap,, ap2and up3be the principal stresses at time
t, with upsG Iup2< up, and &, , l?p:p2
and i$,, the
corresponding uniaxial Youngs moduli. Hence, the
material is considered as orthotropic with the directions of orthotropy defined by the principal stresses.
The stress-strain matrix corresponding to these directions is, considering three-dimensional stress conditions,

v E,*

v E,,

(1 - v)J$~

v E23

(1 -v)J$3

$1 - 2v) E,*

gl - 2v) E,3

5?,,:pl

1
c=(l+)(l-2)
symmetric

$1 - 2v) E23
(7)
2.1.2. Multiaxial conditions. The behavior of concrete under multiaxial stress conditions is very complex and has not been assessed experimentally in a
complete manner. Various material models with considerable simplifying assumptions have been proposed in the literature. Considering the variability of
concrete materials that need be described in practice,
the objective in this work was to develop an effective
but simple model that provides sufficient flexibility to
the analyst to fit the basic material behavior.
The stress-strain relations are evaluated differently
depending on whether the material is loading or
unloading.
To characterize loading and unloading conditions
we define a loading function J

where v is the Poisson ratio, and the shear modulus

in a coordinate plane is calculated from the weighted
Youngs modulus corresponding to that plane,

G, =

Eij
2(1+V)

= ~

1Upi1$i + 1'Upj1

2(1+V)

(*)

The above stress-strain relations for material loading conditions are only employed in the calculation of
the stiffness matrix at time t. Considering the evaluation of the stress increment, u, from time t to time
t + At, the stress integration is approximated in the
following manner,
u = ee

If=

Epj

IUpiISIUpjI

(9)

(5)

where 9 = (4 s~s~)2, Sij= dii - a,&, a, = l/3 u,,,

and 6, is the Kronecker delta. The material is loading
except when the unloading conditions are deter-

where e is the strain increment. The stress-strain

matrix C is as C in eqn (7) but the uniaxial Youngs
moduli $ are replaced by the moduli &. These are
evaluated using the uniaxial stress-strain relationship

KLAUS-JORGEN
BATHEet al.

566

(a

ENStLE

FAILURE

Fig. 2. Triaxial failure envelopes reduced to two-dimensional conditions.

Fig. 3. Triaxial tensile failure envelope.

in Fig. 1, namely, for GoB 0 we have

qi = $

(10)

and if IcPic 0, we define

r+Ar

f
ap1=

and then

(12)
Here, +ePi and ePiare the strain components at time
t + At and t in the directions of the principal stresses.

INPUT
---KHAN

2.2. Material failure emelopes

To model the failure of the material in tension and
compression in two- and three-dimensional analysis
and to account for multiaxial conditions in the
uniaxial stress-strain law of Fig. 1, failure envelopes
are employed. Figure 2 shows an example of the
failure envelopes used corresponding to two principal stress directions. Regarding the tensile failure it is
noted that considering one principal stress direction
the tensile strength of the material in this direction
does not change with the introduction of a tensile
stress in the other principal stress direction, but a
compressive stress decreases this tensile strength.
The failure envelopes shown in Fig. 2 are a special
case of the triaxial failure envelopes for tensile and
compressive failures shown in Figs 3 and 4.

TO ADINA
AND SlWGY

Fig. 4. Triaxial compression failure envelope.

567

Nonlinear analysis of concrete structures

The triaxial tensile failure envelope is a straightforward generalization
of the envelope for twodimensional
analysis. However, the compressive
failure envelope is more complex and requires more
data input. The envelope can be used to represent a
large number of different envelopes like the biaxial
envelope of Liu et al. [7] and the triaxial failure
surface of Khan and Saugy[8]. The shape of the
compressive failure surface used is largely based on
the experimental results reported by Kupfer et al. [6]
and Launay and Gachon [9].
The compressive failure envelope is input using 24
discrete stress values. Firstly, the values crb,/dc are
input. These values define at what stress magnitudes
a,, the discrete two-dimensional failure envelopes for
additional stresses upzand cprare input. These failure
envelopes are defined by the failure stress values
bp$/c?~(i = 1, . . . ,6;j = 1,2,3) that correspond to the
stress magnitudes up1= tag,, ap2= /I ups(B is a constant) and up2= up,.
The failure envelopes are employed to establish the
uniaxial stress-strain law accounting for multiaxial
stress conditions, and to identify whether tensile or
crushing failure of the material has occurred. Having
calculated the current principal stresses, to establish
the uniaxial stress-strain law it is assumed that up,
and up2 are held constant and the minimum stress
that would have to be reached in the third principal
stress direction to cause crushing of the material is
calculated using the failure envelopes (see Fig. 4). Let
this stress be 8: and y, = a:/~?,, then we also use

5; =y,d,
e: = (C, r: + C,Yl)%

(13)

p:=(C,Y:+C,Yl)&

(14)

where C, and C, are constants defined by the user;

-[I -v2
1

COMPRESSlC+4
STRAIN

Fig. 5. Increase of strength parameters for model under

multiaxial conditions.

state in the failure envelopes. In the following we

consider how one single plane of tensile failure develops and how the material fails in compression crushing.
Tensile failure occurs if the tensile stress in a
principal stress direction exceeds the tensile failure
stress. In this case it is assumed that a plane of failure
develops perpendicular to the principal stress direction. The effect of this material failure is that the
normal and shear stiffnesses across the plane of
failure, and the corresponding normal stress and
shear stresses, are reduced.
Prior to tensile failure the stress-strain material law
is given by eqns (7)-(12). Assuming that upi is larger
than the tensile failure stress, the new material
stress-strain matrix is in the stiffness matrix calculation,

I!&

v$3
E2,

0
0

i.2
0%
2(1 + v)

symm.

normally, C, = 1.4 and C, = -0.4. The constants

8:) Cl, Pi, 5: are employed instead of the unprimed
variables in order to establish, using eqn (4), the
uniaxial stress-strain law under multiaxial conditions
(see Fig. 5).
To identify whether the material has failed, the
principal stresses are used to locate the current stress

where the Epi are the uniaxial Youngs moduli

evaluated in the principal stress directions using eqns
(2) and (4) and the 'E,,
are evaluated using eqn (8).
The constant rl. is a small value, normally lo-, to
reduce the stiffness normal to the tensile failure plane,
and the constant b, normally 0.5, reduces the shear
stiffness in the tensile failure plane.

568

KLAUS-JORGEN
BATHEet al.

For the stress calculation, the following stressstrain matrices are used.
For the tensile stress normal to the tensile failure
plane and the shear stresses in this plane, we use the
total strains to calculate the total stresses with
C=IyIrn.

Gp/, G;]

(6)

where Ef and G$, G{, are evaluated as shown in

Fig. 6. In this figure 5 is a variable, input by the user,
that defines the amount of tension stiffening. Figure
6 shows that G$ is a fraction of the initial shear
modulus, &,/2(1 + v), corresponding to the total
strain normal to the tensile failure plane. Also, Fig. 6
assumes loading from zero stress directly into the
tensile region. If the tensile stress is reached by
unloading from a compressive stress, the strain normal to the tensile failure plane is measured from the
strain at which the stress is zero (see the example of
Fig. 7).
For the remaining stress components, the increment in the stresses is evaluated from the incremental
strains using

symm.

(17)

(1 - v)G,,

where the r$i are the uniaxial Youngs moduli evaluated using eqn (12) and E2, and TG23are evaluated
using eqn (8) but with r&i instead of J!&.
Considering eqns (15), (16) and (17), we note that
by use of eqn (16) the tensile stress normal to the
tensile failure plane and the shear stresses in this
plane are gradually released. Also, the equations
show that plane stress conditions are assumed to exist
at the plane of tensile failure. The factor tf, is not set
exactly equal to zero in order to avoid the possibility
of a singular stiffness matrix. The value to be employed for nS must depend on a number of physical
factors, and in the numerical solution it is at this
time best to leave nJ as a variable that is input by the
analyst.
For visualization, the plane of tensile failure is
referred to as a crack, but we should interpret this
terminology judiciously because a physical crack does
not actually develop at the element integration point.
Instead, the material has failed in one principal stress
direction.
Equations (15), (16) and (17) describe the solution
when tensile failure occurs. To identify compression
failure, the largest principal stress IuP, is employed to
establish from Fig. 4 by interpolation, the biaxial
failure envelope on up2 and ups. The material has
crushed if the stress state corresponding to rr,,*and
a,, lies on or outside the biaxial failure envelope.

STRAIN NWIMAL
RILURE

/
St

t0

PLANE

1
I

FAILURE PLANE

Fig. 6.

Material moduli for stress calculation under tensile

failure. (a) Calculation of Youngs modulus E,, normal to
tensile failure plane. (b) Calculation of shear modulus in
tensile failure plane.
failure is still active. The failure is considered to be
inactive provided the normal strain across the plane
becomes negative and less than the strain at which the
last failure occurred and is active otherwise. Therefore, a tensile failure plane may repeatedly be active
and inactive.
If a tensile failure plane has developed, which may
or may not be active, the material stress-strain
relations are always established as described above
but corresponding to the principal stress directions in
the failure plane and the direction perpendicular to
this plane. Hence, instead of using the principal
stresses and corresponding directions as done for the
unfailed material, the stress conditions along and
normal to the material tensile failure plane are used
to evaluate the stress-strain matrix. Also, when a
failure plane is or was active, a subsequent failure
plane is assumed to form perpendicular to the direction of the one that developed first, once a normal
stress along the original failure plane has reached the
tensile failure stress. It follows that at any integration
point, the direction of the third tensile failure plane
is fixed once failure has occurred in two directions.
It may also happen that after tensile failure of the
material (in one or two directions) the material fails
in compression crushing, which is identified, as usual,
by entering the compression crushing envelope in
Fig. 4 with the principal stress(es) that act(s) along
the tensile failure plane(s).
If the material has crushed in compression, it is
assumed that the material strain-softens into all
directions.
2.4. Post-compression

crushing,

strain-softening

2.3. Post-tensile cracking behavior

behavior

Once a tensile plane of failure has formed, it is

checked in each subsequent solution step whether the

Consider first uniaxial stress conditions. As shown

in Fig. 1, for a uniaxial strain smaller than &, the

Nonlinear analysis of concrete structures

569

2.7. Example solutions

In the following, we present some simple, small-size
applications to demonstrate the details of the model
behavior when it is subjected to various load conditions. In all these examples a constant Poisson ratio
v is used.
Figure 7 summarizes the response obtained when
subjecting a single plane stress element to biaxial
strain-controlled loading. This solution includes initiation, closing and reopening of a crack; loading,
unloading and reloading with large load steps; and
crushing of the material with strain-softening behavior. In this figure, the numbers indicate load steps
where epj and eps are the strain component and in the analysis. Note that the material model can
incremental strain component at time t measured in be used with relatively large load step increments.
Figure 8 shows the material failure for different load
the direction of the principal stress rrp,.
Note that, as shown in Fig. 1, when epsbecomes steps indicated in Fig. 7.
In the next example solution a fully constrained
equal to or less than P:, the stresses are set to zero.
However, an option is available for a gradual stress concrete specimen modeled by four 4-node plane
unloading in which case the constant Youngs mod- stress elements is subjected to a uniform temperature
ulus & = -5:/c: is employed, i.e. in eqn (18) we then loading 8. The initial temperature is zero degrees and
use IS=&.
the temperature is raised to 450C leading to material
crushing, see Fig. 9. Note the increase in the crushing
If unloading of the crushed material in the strainsoftening region occurs characterized by e,r > 0, the strength of the material is due to the biaxial stress
initial Youngs modulus $ is used. In all these conditions.
Next, in Fig. 10, a simple plane stress concrete
calculations the principal stresses are checked individually whether a positive value has been reached, and element subjected to biaxial compression loading is
if so, the stress in the corresponding direction is set considered. After crushing, the material relation of
Sec. 2.4 is used, i.e. all Youngs moduli are set to &.
to zero.
In a further example solution, a single 2D plane
stress element is loaded to trace the post-cracking
2.5. Poissons ratio in the compressive region
behavior of the concrete material model under biaxial
It has been observed in experiments that the ratio
loading in comparison with uniaxial loading (see Fig.
of lateral strain to principal compressive strain re- 11). The shift in the stress-strain 1D response demains constant until approximately 80% of the max- pends on the crp,, up2relation (positive/positive, posiimum compressive stress d,.
tive/negative) as well as on the magnitude of the
Usually we assume in our analyses that the Poisson
Poisson ratio. It should be noted that strain-softening
ratio is constant; however, as an option a value, vr, post-cracking models (r > 1.0) are to be used with
can be used when the material dilates under comprescaution since such models may lead to a nonunique
sion. The value of vr is given by [lo],
solution [14] (see also Sec. 3.2).
Extensive cracking and complex post-cracking beQ
havior can be observed in the analysis of concrete
vs = v whenY~=:QY.
pressure vessels. In this part of the paper we present
detailed solutions of two reduced models. Analyses
of the complete pressure vessel structures are prev~=v,-(v~-v)/~
when y2>yo
sented in Part II of the paper.
Figure 12 shows a simple finite element model of
(19)
a part of the Sandia pressure vessel. This model reprewhere v is the initial Poissons ratio, v, is the maxi- sents a slice of the vessel at the elevation 5.969 m. The
mum value of Poissons ratio at failure and y. is slice is subjected to boundary and loading conditions
corresponding to the whole structural behavior when
usually set to 0.8.
the vessel is under internal pressure. Three eight-node
axisymmetric elements were used to model the con2.6. Temperature eflects
crete of the pressure vessel. The l-node axisymmetric
In some analyses temperature strains, e, need to and 2-node truss elements were used to model the
be included. These strains are taken into account by steel reinforcements. Figure 13 shows crack distribureplacing the total incremental strain, e, in the gov- tions for different load levels. In Fig. 14 a comparison
erning incremental equations by the strain e - elh, with experimental data is presented. Note that in
where e is the thermal incremental strain and is Fig. 13 cracks are plotted at the integration points,
calculated from the temperature conditions.
the size of a cracking symbol depends on the con-

material has crushed and softens with increasing

compressive strain, i.e. 1 is negative.
Under multiaxial stress conditions the compression
crushing is identified using the multiaxial failure
envelope, and once the material has crushed isotropic
using
the uniaxial
are assumed
conditions
stress-strain law with the constants a:, Pi and so on
(see Fig. 5) corresponding to the multiaxial conditions at crushing. The Youngs modulus for the
isotropic material law is evaluated using

KLAUS-JORGENBATHE et al.

s
zz-SS3MlS

II-SS3YlS

zz-SS3LllS

AA-SS3YIS

Nonlinear analysis of concrete structures

571

C.A.S. 32,3.4-F

512

KLAUS-JURGEN
BATHEet al.

513

Nonlinear analysis of concrete structures

.I0

.I5

.I3

TMREE S-NODE
AXISVMTRIC
CONCRETE
ELEMENTS

.I@

SEVEN

2-NODE

(LINES)

AND
ELEVEN
I-NODE
PRESTRESSING
ELEMENTS

(RINGS)
TRUSS

Fig. 12. A finite element model for a slice of the Sandia pressure vessel at elevation 5.969 m (see also Fig.
19). Top: 2D axisymmetric concrete elements. Bottom: prestressing truss elements.

tributing volume at the integration point, circles

represent hoop cracks and two parallel lines a short
distance apart represent cracks in the Y-Z plane. A
smaller distance between the parallel lines is used
when the crack is closed.
A simple model of a part of the PWR pressure
vessel (this vessel is also considered in Part II of the
paper, see Fig. 30) is shown in Fig. 15. As for the
Sandia vessel, a slice of the PWR vessel is considered.
This slice now corresponds to the height 22.30 m.
Two 8-node axisymmetric finite elements were used
to model the concrete of the vessel. The steel hoop,
longitudinal and radial reinforcement is modeled
using l-node axisymmetric and 2-node truss elements. The model is subjected to boundary and
loading conditions that correspond to the response of
the complete structure. The crack distributions for
different solution times are shown in Fig. 16. Note
that three cracks have been formed at some of the
integration points at the final load level. In Fig. 17 the
solution results are compared with the results of an
analysis using a simple force-balancing method.
2.8. Concluding remarks
The objective of Part I of this paper was to
summarize the current status of our concrete material

model. The model can be employed in a number of

important concrete analyses, notably for collapse
analyses of two- and three-dimensional
concrete
structures. One important aspect is the generality of
the model with respect to the stress-strain law and
failure surfaces used. Also, the model is deemed
computationally efficient because relatively large incremental load steps can be employed (see Figs 14
and 17).
Although a powerful model already, further improvements to the concrete model are surely desirable. Some features that would appropriately gain
further research attention are the modeling of timedependent behavior, a more accurate representation
of the post-crushing response and further refinements
in modeling the tensile-stiffening behavior, the effects
of Poissons ratio and cyclic loading conditions.
3.

PART n: APPLICATIONS OF ADINA 5.0 TO THREE

TEST STRUCFURE5

Taylor Woodrow plc (TW), a worldwide engineering group, has been involved in the design and
construction of nuclear power structures for over 30
years. During this period much knowledge has been
gained on the short- and long-term behavior of

I<,1

I<,
POINTS

CRACKlNG
AN0 STRESSES
AT INlECRATlON
POINTS
FOR PRESSURE
6.6261+05
Nltl*+Z

STRESSES
AT INTECRATlON
FOR PRESSURE
3.17E*OS
N/II**2

I 1

STRESSES
POlNTS

CRACKING
AN0 STRESSES
AT INTEGRATION
POINTS
FOR PRESSURE
I .S.BEGS
N/t+--2

FOR PRESSURE
3.66E*05
N/H**2

CRACKlNG
AN0
AT INTEGRATION

LY

Fig. 13. ~elopment

of cracks in the Sandia model of Fig. 12. In 3D analysis a crack at an integration point is represented by two disks a small distance apart,
the crack lying in the plane between the disks. Hence, in this figure parallel lines represent cracks in the Y-Z plane; circles represent hoop cracks. The size of the
two disks corresponds to the contributing volume at the integration point.

4 I t* I

(bl

Nonlinear analysis of concrete structures

STEP
12
OFF GRAPH
ELEVATION

lS.Q6QM

0 ADINA.
I2
+ EXPERImNlAL

STEP

SOLUTION
RESULTS

lIO5

PRESSURE

lN/M**2)

Fig. 14. Displacement-pressure response of the Sandia pressure vessel.

concrete and a comprehensive degree of expertise has

been accumulated in the design, analysis and testing
of a wide range of concrete structures. TW has
extensive in-house testing facilities which at first were
used in support of the nuclear industry and later have
been involved in a wider range of concrete structures
including
the offshore, civil and conventional
building areas. This test work has supplemented
the analytical work which itself has augmented the
testing.
To increase its analytical capabilities it was decided
to acquire a nonlinear finite element system which
included a concrete material model. After reviewing
a number of alternatives TW obtained the ADINA 8 1
code in 1983. In using the package on a wide range
of reinforced concrete structures a number of potential improvements to the concrete model could be
envisaged, and TW decided to add in extra features
and change existing ones (see [lo, IS]). The form of
these changes and their influence were shared with
ADINA R & D, Inc. and have been incorporated

within some of the improvements to the concrete

material model that are found in the latest ADINA
5.0 code.
The work presented in this section was originally
carried out using the modified ADINA 8 1 system but
has subsequently been rerun using the latest ADINA
5.0 version, these later results being presented now.
Three applications are given of which two represent
situations where the original analysis was carried out
prior to the test results being available. The third
analysis supplemented the design but will eventually
be used in conjunction with a physical model at
present under construction in the Taylor Woodrow
Laboratories.
3.1. One-sixth scale model of a reinforced concrete
containmeht
The U.S. Nuclear Regulatory Commission (NRC)
has initiated the construction and testing of a series
of scale models of containment
buildings. This
formed part of a containment integrity program, one

i.

Fig. IS. A finite element model for a slice of the PWR pressure vesselat elevation 22.30 m (seealso Fig. 32).

KLAUS-JGRGEN
BATHEet ai.

576

Nonlinear analysis of concrete structures

Fig. 17. Loadclisplacement

FORCE-BALANCING

511

METHOD

response for the PWR model of Fig. 15.

_i

WING

LINL

._

Fig. 18. General arrangement of the Sandia pressure vessel. It is a 1/6th scale model.

578

KLAUS-J@RGEN
BATHEet

IDINA

ORIGINAL
05700

XVMN

al.

0.000

XWIAX 3.810
Ywll% 0.000
1VAX t t *x3

Fig. 19. Finite element mesh of the Sandia pressure vessel.

Fig. 20. Compressive failure envelope used for analysis of the Sandia pressure vesse1. Comparison of
experimental values, ADINA input values used in the analysis and ADINA default values.

Nonlinear analysis of concrete structures

item of which was the construction of a one-sixth
scale model of a reinforced concrete containment.
This structure was to be analyzed prior to testing by
a number of participating organizations, four from
the U.S. and six from Europe. One of the groups
submitting an analysis was sponsored by the Central
Electricity Generating Board (CEGB) of the United
Kingdom. Work was carried out for the CEGB by the
Nuclear Design Associates (NDA), a joint venture of
Taylor Woodrow and Sir Robert McAlpine.
The analysis formed part of a Round Robin
pretest series of analyses, these being documented in
a NRC document [16] and formed part of the NDA
validation work for the use of the ADINA code [ 17)
Details of the test structure are shown in Fig. 18. It
can be seen that whilst in general being an axisymmetric structure, some degree of divergence from this
occurs at major penetrations, and with the disposition of reinforcement. The ADINA analysis carried
out was assumed to be axisymmetric, the local threcdimensional effects being investigated by separate
linear analyses.
The finite element mesh produced for the ADINA
analysis is shown in Fig. 19. The concrete properties
used in the analysis were as follows,
initial E value ($) = 24,800 MN/m2
crushing stress limit (8,) = -41 MN/m2
crushing strain (&) = -0.00186
ultimate stress (6) = -30 MN/m2
ultimate strain (&) = -0.0032
tensile limit (5,) = 4.1 MN/m2
Poissons ratio (v) = 0.2.
Compressive failure envelopes were based on published results [18], these being shown in Fig. 20 in
order to indicate how they deviate from the ADINA
default values. In general TW has used these failure
curves for most analyses as it is felt they fit a wider
range of test results than the default values, although
the deviation is only of concern in the higher triaxial
stress states. As an aside, it is TWs policy when using
ADINA to try to have a standard approach wherever
possible, therefore when dealing with conventional
reinforced concrete most input parameters are unchanged from one analysis to the next apart from the
compressive limit. This even covers the choice of
tensile limit where, unless otherwise directed by circumstances or direct evidence, the value is assumed
to be 10% of the compressive limit.
The reinforcement was treated as a multilinear
elastic model. The main vertical and hoop reinforcement could be modeled easily but towards the bottom
of the wall there existed diagonal seismic reinforcement. It was considered sufficiently accurate to model
this as an equivalent extra set of vertical and hoop
bars. Orthogonal bars in the base were likewise
replaced by equivalent replacement hoop and radial
bars.

519

The physical model contained a steel liner.

ADINA does not as yet contain an axisymmetric
shell element, therefore the liner was modeled using
hoop and meridional truss elements with due allowance for the Poissons ratio effect that would be
in the plate but which is missing in the equivalent set
of truss elements.
Finally the soil was modeled as a set of spring
elements. Due to the relatively large thickness of the
bottom cap the choice of spring stiffnesses was not
expected to significantly influence the final ultimate
load.
The model was constructed by Sandia National
Laboratories and tested in July 1987, by which time
all analysis work had been carried out, supplied to the
NRC and disseminated in a report of May 1987 [ 161.
The test consisted basically of an initial low pressure
cycle up to 1.15 times design pressure. With design
pressure at 0.317 MN/m2 this represents an initial
loading of 0.365 MN/m2. After returning to zero
pressure a high pressure test was performed up to the
ultimate pressure sustainable by the structure. In the
analysis the initial cycle was not included, the loading
consisting of a gradual application of pressure until
failure was indicated.
The following commentary will help to describe
both the ADINA analysis findings and the complex
behavior predicted.
0.138 MN/m-Cracking
occurs at the wall base
junction due to bending. This represents less than
44% of design pressure, but subsequent cracking in
this area progresses across the section in a controlled
manner. The outcome is that a plastic hinge forms
here.
0.345 MN/m2-Vertical
cracking (that is cracks
due to hoop stresses) commence at the mid-height
and quickly spread to reach the dome spring line at
0.365 MN/m2.
0.400 MN/m2-Horizontal
cracking (due to bending and axial tension) appears on the outside face of
the wall and quickly spreads through the thickness
and along the wall towards the dome and wall base
junction. By about 0.55 MN/m2 the dome and wall
are almost completely cracked through, the pressure
being resisted by the liner and reinforcement.
0.552 MN/m2-Cracking
occurs at the middle of
the underside of the bottom cap. By 0.855 MN/m2
these cracks have extended across the entire underside.
0.565 MN/m2-Initial
yield of the liner at midheight in hoop direction, again rapidly extending up
and down the structure.
0.855 MN/m*-Inside
hoop reinforcement at midheight reaches first yield.
0.890 MN/m*-All layers of hoop reinforcement at
mid-height have yielded resulting in rapid expansion
of the wall.
1.000 MN/m2-By
this pressure most reinforcement has yielded.

KLAUS-J~GEN BATHEet al.

580

3 ] 1 i .

!^! T

5. 969M
/

*
2.

.
4.

6.

a.

10.

12.

#lo5
PRESSURE

(N/Mrtx21

Fig. 21. Analysis and experimental displacement response at the mid-height of the radial movement.

1.100 MN/m*-Crushing
of outer wall at wall-base
junction followed by complete failure at this location.
It should be pointed out that in most work to date
TW has preferred not to use iterations but to progress
the load at a controlled and sufficiently slow rate with
holds of pressure at chosen pressure levels. This,
particularly with our earlier work, seemed to allow
the analysis to progress further but care was necessary especially at points of rapid change and near
ultimate.
A few figures have been produced to illustrate the
behavior of the containment. Figure 21 shows the
comparison between analysis and experiment of the
mid-height radial movement. This is the dominant
movement in the structure and so a good choice for
comparison purposes. As can be seen the actual test
reached 1.0 MN/m2 at which stage excessive leakage
occurred, probably through a large liner tear near to
a major penetration. The correlation between the two
curves is very high, especially as the small deviation
over the range O-O.4 MN/m* can be accounted for by
the initial low pressure test carried out. Reference to
the table of events indicates that much cracking could
have occurred, therefore on repressurization these
regions would have zero tensile strength.
Two other plots (Figs 22 and 23) show global views
of the displaced shape and the extent of cracking.
Figure 22 shows global views of the displaced shape
at various pressure levels, all displacements being
drawn to the same scale. Figure 23 shows the extent
of cracking at various pressure levels.
3.2. Alkali-silica reacted beam
During the period 1984-1989 the Road Directorate
of Denmark initiated a project to study the load
carrying capacity of structural members subjected to
alkali-silica reaction (ASR). One item of this involved studying the shear strength of concrete beams
subject to ASR.

At a conference in Denmark in 1988, attendees

were invited to predict the failure load of one of the
beams that had been tested, but for which the results
had not yet been disclosed. The results of this test
were revealed in a subsequent report [19] together
with photographs of crack patterns on similar tested
beams.
The analysis was carried out on the test beam
shown in Fig. 24. The beam deviates from a conventional test by virtue of the ASR, and from the lack
of shear reinforcement in the zone of maximum shear.
The mesh used is shown in Fig. 25. The relevant
properties chosen for the concrete were based on
cores removed from the beam. Both vertical and
horizontal cores were taken, the former being 40%
stronger. Being unable to incorporate this nonisotropic behavior a mean compressive strength of
33.0 MN/m* was used together with a tensile strength
of 3.3 MN/m*.
The analysis reached an ultimate load of 132 kN at
each load position. In the experiment the specimen
first failed with a shear failure on one side at
127.4 kN.
Figure 26 shows the load-deflection curve at midspan and it can be seen that the final failure was fairly
abrupt even though a great deal of cracking had
occurred in both the flexure and shear zones. This
could be interpreted as basically indicating a shear
failure. Figure 27 shows the extent of cracking at
various load levels whilst Fig. 28, which plots the
largest strains, gives a clear picture of the gradual
flexural damage, initiating on the mid-section bottom
face very early on, and progressing slowly through
the depth and outwards towards the supports. The
subsequent shear damage is seen to be. quite sudden.
Comparisons with a similar beam can be made in
Fig. 29 which shows the same form of progressive
damage with a rapid change on the left hand side at
ultimate load.
On studying the Danish report it was seen that the

581

Nonlinear analysis of concrete structures

1
-...::..+
7 ..:.>..
.>.
*+.

y..
:..

PRESSURE*O.345

RESSURE-

NPA

________
_-_._._

: __________--.-

1__________
_...

_______e-m-e....

0.69

WA

J
i,

Fig. 22. Global view of the displaced shapes of the Sandia pressure vessel. Displacements magnified by
a factor of 167.

KLAUS-JORGENBATHEet al.

582

i,

PRESSURE.0.345

WA

PRESSURE-O.366

WA

-i

PRESSURE-O.

5175

WA

PRESSURE.0.8549

MPA

c,

Fig. 23. Global view of the extent of cracking in the Sandia model. In the original, yellow-l
blue-2 cracks, red-3 cracks.

crack,

583

Nonlinear analysis of concrete structures

LOAD DIAGRAM

nchoraqe

Zone

5hear

Zone

Bendlnq

Shear

Zone

Zone

Anchoraw

Heasure

I"mm

Zon

Fig. 24. Alkali-silica reacted beam. Geometry, reinforcement and load diagram.

Fig. 25. Finite element model, 8-node plane stress elements for the beam of Fig. 24.

DISPLACEHENT

MID-SPAN

Fig. 26. Load-displacement

HEASURED

DISPLACEHENT

AT

POINT

IHETRES)

response as predicted by ADINA for the alkali-silica reacted beam.

584

KLAUS-J~GEN BATHEet al.

585

Nonlinear analysis of concrete structures

% Load of
Utimate

75

88

100

100

Fig. 29. Experimental results of developments of cracks during loading of an alkali-silica reacted beam.

ASR beams were judged against unaffected specimens. On testing these it was found that very early
shear failure initiated in the shear zone resulted in
lower strengths even though the unaffected concrete
had a higher strength. It was judged that the ASR
caused a great deal of micro-cracking of the concrete
making the beam more ductile, hence redistributing
the stresses more efficiently in the shear zone.
To investigate this it was decided to rerun the
analysis with a higher concrete strength of
58.2 MN/m2, this being the quoted cylinder strength
of the undamaged concrete. Just making this change
resulted, as expected, in an increase in ultimate load,
contrary to the test results. On reflection, the shear
zone is basically unreinforced concrete, therefore
cracking should be concentrated in fewer cracks and
should propagate more rapidly. This is very similar to
the situation of no tensile stiffening after failure and
a reduced shear reduction factor and hence a run was
made to use a tensile stiffening factor of 1.2 instead
of 8 and a shear reduction factor of 0.01 instead of
0.5. Whilst the problem showed much more distress
in the shear zone it reached higher than expected
loads.

Study of the crack dispositions and strains still

indicates a smearing of cracks instead of the expected
dominant single crack. Whilst of not too much
significance in conventionally reinforced structures,
this phenomenon could be more important in analyzing problems involving mass concrete.
3.3. Size&l
ment

B PWR prestressed concrete contain-

The main difference between this containment and

the Sandia test structure is that the Sizewell structure
is prestressed by longitudinal, radial and hoop tendons. Nonlinear analyses have been applied to the
full scale structure now under construction
at
Sizewell, United Kingdom and to a tenth scale model
which, at the time of writing, is under construction at
the Taylor Woodrow Laboratories.
Figure 30 shows the form of the structure which is
in essence axisymmetric except for the major wall
penetration (not indicated on this diagram), the
tunnel in the base and the disposition of prestress.
The latter is provided by hoop tendons that cover
overlapping 240 segments of the wall and the lower

.T

Fig. 30. PWR prestmssed concrete containment section through reactor building.

Fig. 31. PWR containment,

tenth scale model: some structural details.

587

Nonlinear analysis of concrete structures

DINA

ORIGINAL
0.5126

XVMlN
XVMAX
YVHIN
YVMAX

0.000
2.512
0.000
6.981

i,

Fig. 32. Finite element mesh for the PWR containment model.

portion of the dome together with longitudinal tendons that cross the dome orthogonally in plane. The
hoop tendons are anchored on three buttresses and
the longitudinal tendons at the prestressing gallery
below the base.
The model details are shown in Fig. 31 which
follows closely the full size structure except for details
in the bottom cap. The mesh used for the tenth scale
model analysis is shown in Fig. 32. A similar mesh
was employed on the full size structure, which was
also being subjected to a more detailed analysis in
which the element size was roughly halved throughout. A comparison of the two sets of results showed
that for this structure the use of two elements through
the wall was sufficient.
Figures 33 and 34 show progressive deflection
behavior and crack build up at various pressure
levels. The accuracy of these predictions will only be
apparent after the forthcoming pressure tests but
during the design process certain hand checks were
made. One of these involved the radial displacement
of the wall at mid-height which is sufficiently removed
from the wall-base junction and the dome spring line
as to make the behavior amenable to solution as a
long cylinder.
C.A.S.
32,M-G

During the pressurization process there are well

defined stages at which the stress distribution could
be calculated, these being during the early elastic
period when the load is shared by the concrete,
reinforcement and prestress, following the total
cracking of the concrete with transfer of the load fully
onto the reinforcement and prestress and finally after
yielding of the reinforcement when the load is carried
solely by the prestress. These conditions result in spot
checks and are shown in Fig. 35 which plots the
predicted radial deflection of the actual vessel at
mid-height against the pressure. Close agreement is
indicated with the ADINA analysis, the latter showing the effect of the post-tension stiffening during the
period in which the load is transferred from the
concrete onto the steel.
As in most analyses, there are areas of approximation and concern where known behavior has by
necessity been ignored. In this analysis one such
approximation is centered around the prestressing
system. The three-dimensionality of the prestressing
system both from the hoop and particularly in the
meridional prestressing could be addressed by means
of three-dimensional analyses although due to the
various cyclic symmetries involved this may need to

588

KLAUS-JORGEN

PRESSURE-.345

MPA

BATHEet al.

PRESSURE-O.

Fig. 33. Progressive deflection behavior at various pressure levels for the PWR containment model as
predicted by ADINA. Disphcements magnified by a factor of 66.

Nonlinear analysis of concrete structures

PRESSURE.0.345

PRESSURE-O.69

HPA

MPA

PRESSURE-O.592

PRESSURE-O.76

HPA

WA

Fig. 34. Progressive crack build up at various pressure levels for the PWR containment model. In the
original, yellow-l crack, blue-2 cracks, red-3 cracks.

KLAUS-J~JRGEN
BATHEet al.

590

RAOIRL

DISP.

FIT

22.30M

consideration are solution strategies, especially with

new automatic loading methods becoming available,
bond between reinforcement and concrete and nonorthogonal cracking.
Concrete in all its many forms still remains a very
complex material to model numerically but given the
open exchange of ideas and mutual respect between
code developers and structural analysts each successful application helps to meet tomorrows challenges
with an extra degree of confidence.
REFERENCES
I. K. J. Bathe, Finite Element Procedures in Engineering
Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982).

RD1NQS.C
0

FORCE-BRLANCING

METHOD
I--r-r-r-

150.
120.
0 ISPLACEMENT
(MM1
Fig. 35. Radial deflection of the PWR containment at
mid-height as a function of applied pressure.
be at the best a 180 analysis and may even require
a full 360 investigation. The test structure includes
all the features of the main structure that could
influence this and the need for a full threedimensional analysis should be evident after the
pressure test.
A more complicating feature in this instance is the
fact that the prestress is unbonded. This is partly
overcome for the vertical portion of the prestress as
only end loads are considered but for the curved
tendons the interaction between displacements and
prestress forces is, whilst not being impossible, a very
complex calculation. ADINA contains means of carrying this out as loads can be made deflection dependent and this has been demonstrated to be applicable
to unbonded prestress, but the algorithm to cover the
general situation is rather formidable.
3.4. Discussion and conclusions

2. K. J. Bathe and S. Ramaswamy, On three-dimensional

nonlinear analysis of concrete structures. Nucl. Engng
Des. 52, 385-409 (1979).
3. W. F. Chen, Plasticity in Reinforced Concrete. McGraw-Hill (1982).
4. Z. P. Bazant and J. C. Chem, Concrete creep at variable
humidity: constitutive law and mechanisms. Mater.
Struct. RILEM. 18(103), 1-19 (1985).
5. S. Gallegos-Cazares and W. C. Schnobrich, Effects of
creep on the behavior of reinforced concrete gable roof
hyperbolic-naraboloids.
Civil Ennineerina Studies.

Siructural Besearch Series No. 543,Unive&ty of Illi:

6.

7.

8.

9.
10.
11.
12.

13.

14.
15.

The three structures detailed in the previous sec-

tions have helped to build up a high degree of

confidence in the ADINA system and the concrete
model. There is still the need of expertise in its
successful application and there are areas that need
further investigation. In particular the effect of more
discrete cracks as would occur in mass concrete needs
more study and from work (not yet published) involving explosive loadings, there is a need to model the
concrete past the ultimate strain state as a commaterial.
This may simply
(Taylor
paction
Woodrows word, not ADINA R & Ds) involve
with
the
concrete
model
marrying
the
Drucker-Praeer cao model. Other tonics for future

16.

17.
18.

19.

nois at Urbana-Champaign (1988).

H. Kupfer, H. K. Hilsdorf and H. Rusch, Behavior of
concrete under biaxial stresses. J. Am. Concr. Inst. 66,
656666 (1969).
T. C. Liu, A. H. Nilson and F. 0. Slate, Biaxial
stress-strain relations for concrete. J. Struct. Div.,
ASCE 98, 1025-1034 (1972).
M. H. Khan and B. Saugy, Concrete for nuclear
reactors. Am. Concr. Inst. Special Publication 34
(1972).
P. Launay and H. Gachon, Concrete for nuclear mattars. Am. Concr. Inst. Special Publication 34 (1972).
J. S. Gedling, N. S. Mistry and A. K. Welch, Evaluation
of material models for reinforced concrete structures.
Comput. Struct. 24, 225-232 (1986).
J. C. Pyle, Analysis of deep composite steel and concrete
slabs. Comput. Struct. 13, 661-671 (1981).
M. Keuser, B. Kepp, G. Mehlhom and R. Rostasy,
Nonlinear static analysis of end-fittings for GFRP
prestressing rods. Comput. Struct. 17, 719-730 (1983).
J. W. Tedesco and W. G. McDougal, Nonlinear dynamic analysis of concrete armor units. Comput. Struct.
21, 189-201 (1985).
Z. P. Bazant, T. B. Belytschko and Ta-Peng Chang,
Continuum theory for strain-softening. J. Engig Me&.
Div.. AXE 110. 16661692 (1984).
N. S. Mistry, Nonlinear analysis of reinforced concrete
structures. Finite Element News 4 (1985).
D. B. Clauss, Round-robin pretest analyses of a 1:6
scale reinforced concrete containment model subject to
static internal pressurisations. U.S. Nuclear Regulatory
Commission NUREG/CR-4913, SAND87-0891 (1987).
S. Kloosterman, Validation report for the ADINA-TW
finite element code. Report PWR/88/83, PMT Document SXB-IC-096064 (19881.
Task Committee on Finite Element Analysis of Reinforced Concrete Structures. A. H. Nilson, Chairman,
State-of-the-art Report on Finite Element Analysis of
Reinforced Concrete, ASCE Special Publication (1982).
Vejdirektoratet (Road Directorate, Ministry of Transport, Denmark), Load carrying capacity of bridges
subjected to alkali-silica reactions. The shear strength
of concrete beams subjected to ASR Interim Report
No. 1 (1986).