ATENA Theory PDF
ATENA Theory PDF
ATENA Theory PDF
Na Hrebenkach 55
150 00 Prague
Czech Republic
Phone: +420 220 610 018
E-mail: cervenka@cervenka.cz
Web: http://www.cervenka.cz
Theory
Written by
Vladimír Červenka, Libor Jendele,
and Jan Červenka
1.11 References 13
2 CONSTITUTIVE MODELS 15
2.1 Constitutive Model SBETA (CCSbetaMaterial) 15
2.1.1 Basic Assumptions 15
2.1.2 Stress-Strain Relations for Concrete 18
2.1.3 Localization Limiters 24
2.1.4 Fracture Process, Crack Width 25
2.1.5 Biaxial Stress Failure Criterion of Concrete 25
2.1.6 Two Models of Smeared Cracks 27
2.1.7 Shear Stress and Stiffness in Cracked Concrete 29
2.1.8 Compressive Strength of Cracked Concrete 29
2.1.9 Tension Stiffening in Cracked Concrete 30
2.1.10 Summary of Stresses in SBETA Constitutive Model 30
ATENA Theory i
2.1.11 Material Stiffness Matrices 31
2.1.12 Analysis of Stresses 32
2.1.13 Parameters of Constitutive Model 33
ii
2.9 Microplane Material Model (CCMicroplane4) 73
2.9.1 Equivalent Localization Element 74
2.10 References 78
3 FINITE ELEMENTS 83
3.1 Introduction 83
3.20 Global and Local Coordinate Systems for Element Load 186
iv
4.4.4 The Crisfield Method. 207
4.4.5 Arc Length Step 208
ATENA Theory v
7.2.2 Adams-Bashforth Integration Scheme 255
7.2.3 Reduction of Oscillations and Convergence Improvement 256
1.1 Introduction
This chapter presents the general governing continuum equations for non-linear analysis. In
general, there exist many variants of non-linear analysis depending on how many non-linear
effects are accounted for. Hence, this chapter first introduces some basic terms and entities
commonly used for structural non-linear analyses, and then it concentrates on formulations that
are implemented in ATENA.
It is important to realize that the whole structure does not have to be analyzed using a full non-
linear formulation. However, a simplified (or even linear) formulation can be used in many
cases. It is a matter of engineering knowledge and practice to assess, whether the inaccuracies
due to a simplified formulation are acceptable, or not.
The simplest formulation, i.e. linear formulation, is characterized by the following assumptions:
The constitutive equation is linear, i.e. the generalized form of Hook's law is used.
The geometric equation is linear that is, the quadratic terms are neglected. It means that during
analysis we neglect change of shape and position of the structure.
Both loading and boundary conditions are conservative, i.e. they are constant throughout the
whole analysis irrespective of the structural deformation, time etc.
Generally linear constitutive equations can be employed for a material, which is far from its
failure point, usually up to 50% of its maximum strength. Of course, this depends on the type of
material, e.g. rubber needs to be considered as a non-linear material earlier. But for usual civil
engineering materials the previous assumption is satisfactory.
Geometric equations can be considered linear, if deflections of a structure are much smaller than
its dimensions. This must be satisfied not only for the whole structure but also for its parts. Then
the geometric equations for the loaded structure can then be written using the original (unloaded)
geometry.
It is also important to realize that a linear solution is permissible only in the case of small strains.
This is closely related to material property because if strains are high, the stresses are usually,
although not necessarily, high as well.
Despite the fact that for the vast majority of structures linear simplifications are quite acceptable,
there are structures when it is necessary to take in account some non-linear behavior. The
resulting governing equations are then much more complicated, and normally they do not have a
closed form solution. Consequently some non-linear iterative solution scheme must be used (see
Chapter Solution of Nonlinear Equations further in this document).
Non-linear analysis can be classified according to a type of non-linear behavior:
Non-linear material behavior only needs to be accounted for. This is the most common case for
ordinary reinforced concrete structures. Because of serviceability limitations, deformations
are relatively small. However, the very low tensile strength of concrete needs to be accounted
for.
Deformations (either displacements only or both displacements and rotations) are large enough
such that the equilibrium equations must use the deformed shape of the structure. However
the relative deformations (strains) are still small. The complete form of the geometric
ATENA Theory 1
equations, including quadratic terms, has to be employed but constitutive equations are
linear. This group of non-linear analysis includes most stability problems.
The last group uses non-linear both material and geometric equations. In addition, it is usually
not possible to suddenly apply the total value of load but it is necessary to integrate in time
increments (or loading increments). This is the most accurate and general approach but
unfortunately is also the most complicated.
There are two basic possibilities for formulating the general structural behavior based on its
deformed shape:
Lagrange formulation:
In this case we are interested in the behavior of infinitesimal particles of volume dV . Their
volume will vary dependent on a loading level applied and, consequently, on the amount of
current deformations. This method is usually used to calculate civil engineering structures.
Euler formulation:
The essential idea of Euler's formulation is to study the "flow" of the structural material through
infinitesimal and fixed volumes of the structure. This is the favorite formulation for fluid
analysis, analysis of gas flow, tribulation etc. where large material flows exist.
For structural analysis, however, Lagrangian formulation is better, and therefore attention will be
restricted to this. Two forms of the Lagrangian formulation are possible. The governing
equations can either be written with respect to the undeformed original configuration at time t =
0 or with respect to the most recent deformed configuration at time t. The former case is called
Total Lagrangian formulation (TL) while the latter one is called the Updated Lagrangian
formulation (UL).
It is difficult to say which formulation is better because both have their advantages and
drawbacks. Usually it depends on a particular structure being analyzed and which one to use is a
matter of engineering judgement. Generally, provided the constitutive equations are adequate,
the results for both methods are identical.
ATENA currently uses Updated Lagrangian formulation, (which is described later in this
chapter) and supports the highest, i.e. 3rd level of non-linear behavior. Soon, it should also
support Total Lagrangian formulation.
2
[0X3, tX3, t+tX3]
M0 M1 M2
Configuration t+t
Configuration 0 Configuration t
0 t t+t
[0X2, tX2, t+tX2]
[ X1, X1, X1]
ATENA Theory 3
In this document all the following derivations will be presented in their displacement form and
consequently the principle of virtual displacements will be used throughout.
The following section deals with the definition of the stress and strain tensors, which are usually
used, in nonlinear analysis. All of them are symmetric.
where
0
is the ratio of density of the material at time 0 and t ,
t
t
mn is the Cauchy stress tensor at time t ,
0
t X i , m is the derivative of coordinates, ref. (1.5).
Using inverse transformation, we can express Cauchy stress tensor in terms of the 2nd Piola-
Kirchhoff stress tensor, i.e.:
t
t
mn
t t t
0 X m ,i 0 S mn 0 X n , j (1.4)
0
The elements 0t X i ,m are usually collected in the so-called Deformation gradient matrix:
X 0 t X T
t T
0 (1.5)
where:
T
0 0
T
, 0 , 0
X1 X 2 X 3
t
X T t X 1 , t X 2 , t X 3
4
0
The ratio can be computed using:
t
0
t det( 0t X ) (1.6)
Expression (1.6) is based on the assumption that the weight of an infinitesimal particle is
constant during the loading process.
Some important properties can be deduced from definition of 2nd Piola-Kirchhoff tensor (1.3):
at time 0, i.e. the undeformed configuration, there is no distinction between 2nd Piola-Kirchhoff
0
and Cauchy stress tensors because 00 X E , i.e. unity matrix and the density ratio = 1.,
t
2nd Piola-Kirchhoff tensor is an objective entity in the sense that it is independent of any
movement of the body provided the loading conditions are frozen. This is a very important
property. The Cauchy stress tensor does not satisfy this because it is sensitive to the rotation
of the body. It is energetically conjugated with Green-Lagrange tensor described later.
They’re some other stress tensor commonly used for structural nonlinear analysis, e.g. Jaumann
stress rate tensor (describes stress rate rather than its final values) etc, however they are not used
in ATENA and therefore not described in this document.
If we calculate the length of an infinitesimal fibber prior and after deformation in the original
coordinates, we get exactly the terms of the Green-Lagrange tensor.
The following equation gives relation between variation of Green-Lagrange and Engineering
strain tensors:
t X m t X n
0
t
0 ij t emn (1.9)
X i 0 X j
These are the strain tensors used in ATENA. From the other strain tensors commonly used in
non-linear analysis we can mention Almansi strain tensor, co-rotated logarithmic strain, strain
rate tensor etc.
ATENA Theory 5
1.5 Constitutive Tensor
Although the whole chapter later in this document is dedicated to the problem of constitutive
equations and to material failure criteria, assume for the moment that stress-strain relation can
be written in the following form:
t
0 Sij 0t Cijrs 0t rs (1.10)
This form is applicable for linear materials or in its incremental form it can be used also for
nonlinear materials. The following important relations apply for transformation from coordinates
to time 0 to coordinates at time t :
t
t
t
Cmnpq x t x tC t x t x (1.11)
t 0
0 m ,i 0 n , j 0 ijrs 0 p ,r 0 q ,s
or in the other direction
0
0
t
Cijrs 0 t 0 0
t xi , m t x j , n t Cmnpq t xr , p t xs , q (1.12)
0 t
Using constitutive tensor (1.11) and Almansi strains tt , we can write for Cauchy stresses (with
respect to coordinates at time t ):
t
ij tt Cijrs tt rs (1.13)
The equation (1.13) is equivalent to the equation (1.10) that was written for original
configuration of the structure. It is very important to know, with respect to which coordinate
system the stress, strain and constitutive tensors are defined, as the actual value can significantly
differ. ATENA currently assumes that all these tensors are defined at coordinates at time t .
Sij
ij dV 0t dt R
t dt t dt
0 0 (1.16)
0
V
6
Sij
ij dV tt dt R
t dt t dt
t t (1.17)
t
V
where 0V , t V denotes the structure volume corresponding to time 0 and t and t dt R is the total
virtual work of the external forces. The symbol denotes variation of the entity. Since energy
must be invariant with respect to the reference coordinate system, (1.16) and (1.17) must lead to
identical results.
Substituting expressions for strain and stress tensors the final governing equation for structure
can be derived. They are summarized in (1.18) through (1.29). Note that the relationships are
expressed with respect to configurations at an arbitrary time t and an iteration ( i ) . Typically, the
time t may by 0 , in which case we have Total Lagrangian formulation or t t (i 1) , in which
case we have Updated Lagrangian formulation, where some terms can be omitted. ATENA also
support “semi” Updated Lagrangian formulation, when t conforms to time at the beginning of
time increment, i.e. the beginning of load step. The following table compares the above-
mentioned formulations:
Table 1.6-1 Comparison of different Lagrangian formulation.
t t
t Sij( i ) t tt ij( i ) t dV t t R (1.18)
t
V
nd
where 2 Piola-Kirchhoff stress and Green Lagrange strain tensor are:
t
t t
Sij( i ) t ( i ) t t ( i )
x
t t i ,m
t
t mn t t x (ji,n) (1.19)
t t t
1 t t ( i )
t tt ij( i )
2
t ui, j t t ( i )
u
t i, j
t t ( i ) t t ( i )
u
t k ,i u
t k, j (1.20)
t ij
t t ( i )
t t ( i 1)
t ij t ij( i )
(1.22)
t ij( i ) t eij(i ) tij( i )
where linear part of the strain increment is calculated by:
ATENA Theory 7
1
e(i )
t ij
2
t ui(,ij) t ui(,ij) t tt uk(i,i1) t uk(i,)j t tt uk(i,j1) t uk(i,i) (1.23)
Note that the term t eij( i ) t eij is constant, i.e. independent of t ui( i ) , hence it is on RHS of
(1.28).
where fbi and fsi are body and surface forces, t S and t V denotes integration with respect to the
surface with the prescribed boundary forces and volume of the structure (at time and t ).
The 1st integral in (1.29) accounts for external work on surface (e.g. external forces), the second
one for work done by body forces (e.g. weight) and the last one accounts for work done by
inertia forces, which are applicable only for dynamic analysis problems).
At this point, all the relationships for incremental analysis have been presented. In order to
proceed further, the problem must be discretized and solved by iterations (described in Chapter
Solution of Nonlinear Equations).
8
1.8 Problem Discretisation Using Finite Element Method
Spatial discretisation consists of discretising primary variable, (i.e. deformation in case of
ATENA) over domain of the structure. It is done in ATENA by Finite Element Method. The
domain is decomposed into many finite elements and at each of these elements the deformation
field is approximated by
t
ui h j t uij (1.30)
where
j is index for finite element node, j 1...n ,
n is number of element nodes,
h j are interpolation function usually grouped in matrix H j h1 ( r , s, t ), h2 ( r , s, t ).....hn ( r , s, t ) ,
r , s, t are local element coordinates.
The interpolation functions h j are usually created in the way that h j 1 at node j and h j 0 at
any other element nodes.
Combining (1.30) and equation for strain definition (1.8) it can be derived:
t t
t ( i ) t BL 0 t tt BL( i11) t tt BNL U
( i 1) t t ( i )
(1.31)
where
t
t t ( i )
is vector of Green-Lagrange strains,
t t
U ( i ) is vector of displacements,
t t ( i 1) t t ( i 1) st
t BL 0 , t BL1 , t BNL are linear strain-displacements transformation matrices (the 1 two of
them) and nonlinear strain-displacements transformation matrix (the last one).
Similar equation can be written also for stress tensor.
t t
t S (i ) t t
t C (i ) t t
t (i ) (1.32)
where:
t t ( i )
tS is vector of 2nd Piola-Kirchhoff stress tensor and
t t
t C ( i ) is incremental stress-strain material properties matrix.
Applying the above discretisation for each finite element of the structure and assembling the
results the continuum based governing equations in (1.28) can be re-written in the following
form:
t t
t
M U ( i ) ( t K L t tt K NL
( i 1)
) t tU ( i ) t t R t t F ( i 1) (1.33)
t 2
where
t 2
U is the vector of nodal accelerations,
ATENA Theory 9
t t
R , t t F ( i 1) is the vector of applied external forces and internal forces,
(i )
, ( i1) superscripts indicate iteration numbers.
Note that (1.33) contains also inertial term needed only for dynamic analysis. Finite element
matrices in (1.33) and corresponding analytical expressions are summarized:
Cijrs t ers( i ) t eij( i ) dV
t
K L U ( i ) t BLT t C t BL dV U ( i )
t t t
V t
V
Sij( i 1) tij( i ) dV
t
t t ( i 1)
K NL U ( i ) t t ( i 1) T
BNL t t
Sij( i 1) t t ( i 1)
BNL dV U ( i ) t t
t t t t t t
V t
V
H H
t t T t t T t t
R t f A dV dA t f B dV t t R
t t
A V
10
technique is described briefly. All details and derivations can be found e.g. (ZIENKIEWICZ,
TAYLOR 1989) and ČERVENKA et. al. 1993.
Let us define a vector of stresses xx at element nodes i such as xx xx ,1 , xx ,2 ,....xx ,n ,
T
where the 2nd index indicates element node number. Let us also define a vector
Pxx Pxx ,1 , Pxx ,2 ,....Pxx ,n , whose component are calculated
T
Pxx ,i hi xx d e (1.35)
e
The nodal value xx (with values of xx at nodes i =1..n ) is then calculated as follows:
xx M Pxx
inv
(1.36)
where:
M ij hi h j d e (1.37)
e
M ij hi ij d (1.39)
e
The above values are output as nodal element stress/strain values. It follows to calculate
averaged stress/stain value i xx , yy ,..... xz i in a global node i that is participated by all
elements k with incidence at the global node i .
i ek
i k
(1.40)
k
ek
ATENA Theory 11
where is vector of stresses i xx , yy ,.....xz i at a node i , ek is volume of element k that
has incidence of global node i . It should be noted that in ATENA the same extrapolation
techniques is used for other integration point quantities as well such as: fracturing strains, plastic
strains and others.
u
u N (1.41)
uD
The problem governing equations can then be written:
K NN K ND uN RN
K (1.42)
DN K DD uD RD
ATENA software supports that any constrained degree of freedom can be a linear combination
of other degrees of freedom plus some constant term:
uDi uDi ,0 k uNk (1.43)
k
where uDi ,0 is the constant term and k are coefficients of the linear combination. Of course, the
equation (1.43) can include also the term u
l
l
l
D , however it is transformed into the constant
term.
The free degree of freedom are then solved by
uN K NN
1
R N K ND RD (1.44)
RD K DN uN K DD uD (1.45)
The ATENA simple support boundary conditions mean that the boundary conditions use only
constant terms are uDi ,0 , (i.e. k 0 ). The complex support boundary conditions use the full
form of (1.43).
The boundary conditions as described above allow to specify for one degree of freedom either
Dirichlet, or von Neumann boundary condition, but not both of them the same time. It comes
from the nature of finite element method. However, ATENA can deal also this case of more
12
complex boundary conditions by introducing Lagrange multipliers. The derivation of theory
behind this kind of boundary conditions is beyond the scope of this manual. Details can be found
elsewhere, e.g. in (Bathe 1982). To apply this type of boundary conditions in ATENA, specify
for those degree of freedom both simple load and complex support boundary condition, the latter
one with the keyword “RELAX” keyword in its definition.
Nice feature about ATENA is that at any time it stores in RAM only K NN and all the elimination
with the remaining blocks of K is done at element level at the process of assembling the
structural stiffness matrix.
A special type of complex boundary conditions of Dirichlet type are so-called master-slave
boundary conditions. Such a boundary condition specifies that all (available) degrees of one
finite node, (i.e. slave node) are equal to degrees of freedom of another node (i.e. master node).
If more master nodes are specified, then these master nodes are assumed to form a finite element
and degrees of freedom of the slave node is assumed to be a node within that element. Its (slave)
degrees of freedom are approximated by element nodal (i.e. master) degrees of freedom in the
same way as displacements approximation within a finite element. The coefficients k in (1.43)
are thus calculated automatically. This type of boundary conditions is used for example for
fixing discrete reinforcement bars to the surounding solid element .
1.11 References
BATHE, K.J. (1982), Finite Element Procedures In Engineering Analysis, Prentice-Hall, Inc.,
Englewood Cliffs, New Jersey 07632, ISBN 0-13-317305-4.
ČERVENKA, J., KEATING, S.C., AND FELIPPA, C.A. (1993), “Comparison of strain
recovery techniques for the mixed iterative method”, Communications in Numerical Methods in
Engineering, Vol. 9, 925-932.
ZIENKIEWICZ, O.C., TAYLOR, R.L., (1989), The Finite Element Method, Volume 1,
McGraw-Hill Book Company, ISBN 0-07-084174-8.
ATENA Theory 13
14
2 CONSTITUTIVE MODELS
s De, s x , y , xy , e x , y , xy
T T
(2.1)
where s, D and e are a stress vector, a material stiffness matrix and a strain vector, respectively.
The stress and strain vectors are composed of the stress components of the plane stress state
x , y , xy , Fig. 2-1, and the strain components x , y , xy , Fig. 2-2, where xy is the engineering
shear strain. The strains are common for all materials. The stress vector s and the material matrix
D can be decomposed into the material components due to concrete and reinforcement as:
s s c s s , D Dc D s (2.2)
The stress vector s and both component stress vectors sc , s s are related to the total cross section
area. The concrete stress s c is acting on the material area of concrete A c , which is approximately
set equal to the cross section of the composite material A c A (the area of concrete occupied by
reinforcement is not subtracted).
The matrix D has a form of the Hooke's law for either isotropic or orthotropic material, as will be
shown in Section 2.1.11.
ATENA Theory 15
Fig. 2-2 Components of strain state.
The reinforcement stress vector ss is the sum of stresses of all the smeared reinforcement
components:
n
s s s si (2.3)
i 1
where n is the number of the smeared reinforcement components. For the ith reinforcement, the
global component reinforcement stress ssi is related to the local reinforcement stress si, by the
transformation:
s si T pi si, (2.4)
Asi
where pi is the reinforcing ratio pi , Asi is the reinforcement cross section area. The local
Ac
reinforcement stress si, is acting on the reinforcement area Asi
The stress and strain vectors are transformed according to the equations bellow in the plane
stress state. New axes u, v are rotated from the global x, y axes by the angle The angle is
positive in the counterclockwise direction, as shown in Fig. 2-3.
16
cos( ) 2 sin( ) 2 2 cos( ) sin( )
T sin( ) 2
cos( ) 2
2 cos( ) sin( ) (2.6)
cos( ) sin( ) cos( ) sin( ) cos( )2 sin( ) 2
s ( u ) u , v , uv , s ( x ) x , y , xy
T T
e(u ) u , v , uv , e( x ) x , y , xy .
T T
The angles of principal axes of the stresses and strains, Fig. 2-1, Fig. 2-2, are found from the
equations:
2 xy xy
tan(2 ) , tan(2 ) (2.9)
x y x y
where is the angle of the first principal stress axis and is the angle of the first principal
strain axis.
In case of isotropic material (un-cracked concrete) the principal directions of the stress and
strains are identical; in case of anisotropic material (cracked concrete) they can be different. The
sign convention for the normal stresses, employed within this program, uses the positive values
for the tensile stress (strain) and negative values for the compressive stress (strain). The shear
stress (strain) is positive if acting upwards on the right face of a unit element.
ATENA Theory 17
The material matrix is derived using the nonlinear elastic approach. In this approach the elastic
constants are derived from a stress-strain function called here the equivalent uniaxial law. This
approach is similar to the nonlinear hypoelastic constitutive model, except that different laws are
used here for loading and unloading, causing the dissipation of energy exhausted for the damage
of material. The detailed treatment of the theoretical background of this subject can be found, for
example, in the book CHEN (1982). This approach can be also regarded as an isotropic damage
model, with the unloading modulus (see next section) representing the damage modulus.
The name SBETA comes from the former program, in which this material model was first used.
It means the abbreviation for the analysis of reinforced concrete in German language -
StahlBETonAnalyse.
18
Unloading is a linear function to the origin. An example of the unloading point U is shown in
Fig. 2-4. Thus, the relation between stress c and strain eq is not unique and depends on a load
ef
history. A change from loading to unloading occurs, when the increment of the effective strain
changes the sign. If subsequent reloading occurs the linear unloading path is followed until the
last loading point U is reached again. Then, the loading function is resumed.
The peak values of stress in compression f’cef and in tension f’tef are calculated according to the
biaxial stress state as will be shown in Sec.2.1.5. Thus, the equivalent uniaxial stress-strain law
reflects the biaxial stress state.
The above defined stress-strain relation is used to calculate the elastic modulus for the material
stiffness matrices, Sect. 2.1.11. The secant modulus is calculated as
c
s
Ec (2.11)
eq
It is used in the constitutive equation to calculate stresses for the given strain state, Sect. 2.1.12.
The tangent modulus Ect is used in the material matrix Dc for construction of an element stiffness
matrix for the iterative solution. The tangent modulus is the slope of the stress-strain curve at a
given strain. It is always positive. In cases when the slope of the curve is less then the minimum
value Emint the value of the tangent modulus is set Ect = Emint. This occurs in the softening ranges
and near the compressive peak.
Detail description of the stress-strain law is given in the following subsections.
ATENA Theory 19
(1) Exponential Crack Opening Law
w 3 w w
1 c1 exp c2 1 c13 exp c2 , (2.13)
f t ' ef wc wc wc
Gf
wc 5.14
f t ' ef
where w is the crack opening, wc is the crack opening at the complete release of stress, is the
normal stress in the crack (crack cohesion). Values of the constants are, c1 =3, c2 =6.93. Gf is the
fracture energy needed to create a unit area of stress-free crack, f t ' ef is the effective tensile
strength derived from a failure function, Eq.(2.22). The crack opening displacement w is derived
from strains according to the crack band theory in Eq.(2.18).
(2) Linear Crack Opening Law
c ef ft ' 2G
wc w , wc ' f (2.14)
ft ' ef wc ft
20
(3) Linear Softening Based on Local Strain
ATENA Theory 21
f1 f2
Parameters: c1 ' ef
, c2
ft f t ' ef
Parameters c1 and c2 are relative positions of stress levels, and c3 is the end strain.
22
2.1.2.5 Compression after Peak Stress
The softening law in compression is linearly descending. There are two models of strain
softening in compression, one based on dissipated energy, and other based on local strain
softening.
ATENA Theory 23
2.1.3 Localization Limiters
So-called localization limiter controls localization of deformations in the failure state. It is a
region (band) of material, which represents a discrete failure plane in the finite element analysis.
In tension it is a crack, in compression it is a plane of crushing. In reality these failure regions
have some dimension. However, since according to the experiments, the dimensions of the
failure regions are independent on the structural size, they are assumed as fictitious planes. In
case of tensile cracks, this approach is known as rack the “crack band theory“, BAZANT, OH
(1983). Here is the same concept used also for the compression failure. The purpose of the
failure band is to eliminate two deficiencies, which occur in connection with the application of
the finite element model: element size effect and element orientation effect.
4 noded element
y
crack
direction
2
1
Lc
Lt x
1 ( max 1) , 0; 45 (2.17)
45
24
An angle is the minimal angle ( min 1 , 2 ) between the direction of the normal to the failure
plane and element sides. In case of a general quadrilateral element the element sides directions
are calculated as average side directions for the two opposite edges. The above formula is a
linear interpolation between the factor 1.0 for the direction parallel with element sides, and
max , for the direction inclined at 45o. The recommended (and default) value of max =1.5.
ATENA Theory 25
Fig. 2-14 Biaxial failure function for concrete.
1 3.65a '
f c ' ef f c , a c1 (2.19)
(1 a ) 2
c2
where c1 , c 2 are the principal stresses in concrete and f’c is the uniaxial cylinder strength. In
the biaxial stress state, the strength of concrete is predicted under the assumption of a
proportional stress path.
In the tension-compression state, the failure function continues linearly from the point
c1 0 , c 2 f c' into the tension-compression region with the linearly decreasing strength:
c1
f c ' ef f c' rec , rec (1 5.3278 ), 1.0 rec 0.9 (2.20)
f c'
where rec is the reduction factor of the compressive strength in the principal direction 2 due to
the tensile stress in the principal direction 1.
26
Two predefined shapes of the hyperbola are given by the position of an intermediate point r, x.
Constants K and A define the shape of the hyperbola. The values of the constants for the two
positions of the intermediate point are given in the following table.
ATENA Theory 27
Fig. 2-16 Fixed crack model. Stress and strain state.
The principal stress and strain directions coincide in the uncracked concrete, because of the
assumption of isotropy in the concrete component. After cracking the orthotropy is introduced.
The weak material axis m1 is normal to the crack direction, the strong axis m2 is parallel with the
cracks.
In a general case the principal strain axes and rotate and need not to coincide with the axes
of the orthotropy m1 and m2. This produces a shear stress on the crack face as shown in Fig.
2-16. The stress components c1andc2 denote, respectively, the stresses normal and parallel to
the crack plane and, due to shear stress, they are not the principal stresses. The shear stress and
stiffness in the cracked concrete is described in Section 2.1.7.
28
2.1.7 Shear Stress and Stiffness in Cracked Concrete
In case of the fixed crack model, the shear modulus is reduced according to the law derived by
KOLMAR (1986) after cracking. The shear modulus is reduced with growing strain normal to
the crack, Fig. 2-18 and this represents a reduction of the shear stiffness due to the crack
opening.
For the zero normal strain, there is no strength reduction, and for the large strains, the
strength is asymptotically approaching to the minimum value f c ' ef cf c' .
30
c - shear stress on the crack plane
0
D c H 1 0 ,
L
0 0 G (2.31)
E
1 ,H E1 (1 2 )
E2
ATENA Theory 31
In the above relation E2 must be nonzero. If E2 is zero and E1 is nonzero, then an alternative
1 E
formulation is used with the inverse parameter 2 . In case that both elastic modules are
E1
L
zero, the matrix Dc is set equal to the null matrix.
L
The matrix Dc is transformed into the global coordinate system using the transformation matrix
T from (2.8).
D c TT D cL T (2.32)
The angle is between the global axis x and the 1st material axis m1, which is normal to the
crack, Fig. 2-16.
32
s
where Dc is the secant material stiffness matrix from Section 2.1.11 for the uncracked or
cracked concrete depending on the material state. The stress components are calculated in the
global as well as in the local material coordinates (the principal stresses in the uncracked
concrete and the stresses on the crack planes).
The stress in reinforcement and the associated tension stiffening stress is calculated directly from
the strain in the reinforcement direction.
Parameter: Formula:
Cylinder strength f c' 0.85 f cu'
Tensile strength 2
ft 0.24 f
' ' 3
cu
The SBETA constitutive model of concrete includes 20 material parameters. These parameters
are specified for the problem under consideration by user. In case of the parameters are not
known automatic generation can be done using the default formulas given in the table above. In
such a case, only the cube strength of concrete f’cu (nominal strength) is specified and the
remaining parameters are calculated as functions of the cube strength. The formulas for these
functions are taken from the CEB-FIP Model Code 90 and other research sources.
Used units are MPa.
The parameters not listed in the table have zero default value.
The values of the material parameters can be also influenced by safety considerations. This is
particularly important in cases of a design, where a proper safety margin should be met. For that
reason the choice of material properties depends on the purpose of analysis and the filed of an
application. The typical examples of the application are the design, the simulation of failure and
the research.
ATENA Theory 33
In case of the design application, according to most current standards, the material properties for
calculation of structural resistance (failure load) are considered by minimal values with applied
partial safety factors. The resulting maximum load can be directly compared with the design
loads.
According to some researchers, more appropriate approach would be to consider the average
material properties in nonlinear analysis and to apply a safety factor on the resulting integral
response variable (force, moment). However, this safety format is not yet fully established.
In cases of the simulation of real behavior, the parameters should be chosen as close as possible
to the properties of real materials. The best way is to determine these properties from mechanical
tests on material sample specimens.
2.2.1 Introduction
Fracture-plastic model combines constitutive models for tensile (fracturing) and compressive
(plastic) behavior. The fracture model is based on the classical orthotropic smeared crack
formulation and crack band model. It employs Rankine failure criterion, exponential softening,
and it can be used as rotated or fixed crack model. The hardening/softening plasticity model is
based on Menétrey-Willam failure surface. The model uses return mapping algorithm for the
integration of constitutive equations. Special attention is given to the development of an
algorithm for the combination of the two models. The combined algorithm is based on a
recursive substitution, and it allows for the two models to be developed and formulated
separately. The algorithm can handle cases when failure surfaces of both models are active, but
also when physical changes such as crack closure occur. The model can be used to simulate
concrete cracking, crushing under high confinement, and crack closure due to crushing in other
material directions.
Although many papers have been published on plasticity models for concrete (for instance,
PRAMONO, WILLAM 1989, MENETREY et al 1997, FEENSTRA 1993, 1998 ETSE 1992) or
smeared crack models (RASHID 1968, CERVENKA and GERSTLE 1971, BAZANT and OH
1983, DE BORST 1986, ROTS 1989), there are not many descriptions of their successful
combination in the literature. OWEN et al. (1983) presented a combination of cracking and
visco-plasticity. Comprehensive treatise of the problem was provided also by de BORST (1986),
and recently several works have been published on the combination of damage and plasticity
(SIMO and JU 1987, MESCHKE et al. (1998). The presented model differs from the above
formulations by ability to handle also physical changes like for instance crack closure, and it is
not restricted to any particular shape of hardening/softening laws. Also within the proposed
approach it is possible to formulate the two models (i.e. plastic and fracture) entirely separately,
and their combination can be provided in a different algorithm or model. From programming
point of view such approach is well suited for object oriented programming.
The method of strain decomposition, as introduced by DE BORST (1986), is used to combine
fracture and plasticity models together. Both models are developed within the framework of
34
return mapping algorithm by WILKINS (1964). This approach guarantees the solution for all
magnitudes of strain increment. From an algorithmic point of view the problem is then
transformed into finding an optimal return point on the failure surface.
The combined algorithm must determine the separation of strains into plastic and fracturing
components, while it must preserve the stress equivalence in both models. The proposed
algorithm is based on a recursive iterative scheme. It can be shown that such a recursive
algorithm cannot reach convergence in certain cases such as, for instance, softening and dilating
materials. For this reason the recursive algorithm is extended by a variation of the relaxation
method to stabilize convergence.
where the increments of plastic strain ijp and fracturing strain ijf must be evaluated based on
the used material models.
If the trial stress does not satisfy (2.38), the increment of fracturing strain in direction i can be
computed using the assumption that the final stress state must satisfy (2.40).
Fi f ii n f ti ii t E iikl kl f f ti 0 (2.40)
This equation can be further simplified under the assumption that the increment of fracturing
strain is normal to the failure surface, and that always only one failure surface is being checked.
For failure surface k , the fracturing strain increment has the following form.
Fk f
ij f ik (2.41)
ij
After substitution into (2.40) a formula for the increment of the fracturing multiplier is
recovered.
ATENA Theory 35
kk t f tk kk t f t( wkmax )
and wkmax Lt (ˆ kk f ) (2.42)
E kkkk E kkkk
This equation must be solved by iterations since for softening materials the value of current
tensile strength f t( wkmax ) is a function of the crack opening w , and is based on Hordijk’s formula
(defined in SBETA model).
The crack opening w is computed from the total value of fracturing strain ˆ kk f in direction k ,
plus the current increment of fracturing strain , and this sum is multiplied by the
characteristic length Lt . The characteristic length as a crack band size was introduced by
BAZANT and OH. Various methods were proposed for the crack band size calculation in the
framework of finite element method. FEENSTRA (1993) suggested a method based on
integration point volume, which is not well suited for distorted elements. A consistent and rather
complex approach was proposed by OLIVIER. In the presented work the crack band size Lt is
calculated as a size of the element projected into the crack direction, Fig. 2-20. CERVENKA V.
et al. (1995) showed that this approach is satisfactory for low order linear elements, which are
used throughout this study. They also proposed a modification, which accounts for cracks that
are not aligned with element edges.
36
It is important to distinguish between total fracturing strain ˆij f , which corresponds to the
maximal fracturing strain reached during the loading process, and current fracturing strain ij f ,
which can be smaller due to crack closure, and is computed using (2.44) derived by ROTS and
BLAUWENDRAAD.
where i j , and sF is a shear factor coefficient that defines a relationship between the normal
and shear crack stiffness. The default value of sF is 20.
Shear strength of a cracked concrete is calculated using the Modified Compression Field Theory
of VECHIO and COLLINS (1986).
0.18 f c
ij , i j (2.48)
24 w
0.31
ag 16
Where f c is the compressive strength in MPa, a g is the maximum aggregate size in mm and w
is the maximum crack width in mm at the given location. This model is activated by specifying
the maximum aggregate size a g otherwise the default behavior is used where the shear stress on
a crack surface cannot exceed the tensile strength.
The secant constitutive matrix in the material direction was formulated by ROTS and
BLAUWENDRAAD in the matrix format.
E s E - E(E cr E) -1 E (2.49)
Strain vector transformation matrix T (i.e. global to local strain transformation matrix) can be
used to transform the local secant stiffness matrix to the global coordinate system.
T
E s T Es T (2.50)
It is necessary to handle the special cases before the onset of cracking, when the crack stiffness
approaches infinity. Large penalty numbers are used for crack stiffness in these cases.
ATENA Theory 37
The value of 0 corresponds to unloading to origin (default value for backward compatibility),
fU =1 means unloading direction parallel to the initial elastic stiffness.
The plastic corrector ijp is computed directly from the yield function by return mapping
algorithm.
F p ( ijt ijp ) F p ( ijt lij ) 0 (2.52)
The crucial aspect is the definition of the return direction lij , which can be defined as
G p ( klt ) G p ( ijt )
lij E ijkl then ijp (2.53)
kl ij
where G ( ij ) is the plastic potential function, whose derivative is evaluated at the predictor stress
state ijt to determine the return direction.
The failure surface of MENETREY, WILLAM is used in the current version of the material
model.
2
3P . ' m
F 15
p
' r ( , e ) c 0 (2.54)
fc 6 f c 3 f c'
where
f c'2 f t '2 e 4(1 e 2 ) cos2 (2e 1) 2
m3 , r ( , e) 1
f c' ft ' e 1
2(1 e ) cos (2e 1) 4(1 e ) cos 5e 4e
2 2 2 2
2
38
In the above two formulas the expression f c( eqp ) indicates the hardening/softening law, which is
based on the uniaxial compressive test. The law is shown in Fig. 2-21, where the softening curve
is linear and the elliptical ascending part is given by the following formula:
2
c eqp
f co f c f co 1
(2.57)
c
f'c
f’c0= 2f’t
eq
p
c =f'c/E
p
The law on the ascending branch is based on strains, while the descending branch is based on
displacements to introduce mesh objectivity into the finite element solution, and its shape is
based on the work of VAN MIER. The onset of nonlinear behavior f c'0 is an input parameter as
well as the value of plastic strain at compressive strength cp . The Fig. 2-21 shows typical values
of these parameters. Especially the choice of the parameter f c'0 should be selected with care,
since it is important to ensure that the fracture and plastic surfaces intersect each other in all
material stages. On the descending curve the equivalent plastic strain is transformed into
displacements through the length scale parameter Lc . This parameter is defined by analogy to
the crack band parameter in the fracture model in Sec. 2.2.3, and it corresponds to the projection
of element size into the direction of minimal principal stresses. The square in (2.56) is due to the
quadratic nature of the Menétry-Willam surface.
Return direction is given by the following plastic potential
1
G p ( ij ) I1 2 J 2 (2.58)
3
where determines the return direction. If 0 material is being compacted during crushing,
if 0 material volume is preserved, and if 0 material is dilating. In general the plastic
model is non-associated, since the plastic flow is not perpendicular to the failure surface
The return mapping algorithm for the plastic model is based on predictor-corrector approach as
is shown in Fig. 2-22. During the corrector phase of the algorithm the failure surface moves
along the hydrostatic axis to simulate hardening and softening. The final failure surface has the
apex located at the origin of the Haigh-Vestergaard coordinate system. Secant method based
Algorithm 1 is used to determine the stress on the surface, which satisfies the yield condition and
also the hardening/softening law.
ATENA Theory 39
Fig. 2-22 Plastic predictor-corrector algorithm.
Fig. 2-23. Schematic description of the iterative process (2.73). For clarity shown in two dimensions.
40
G p ( ijt )
Evaluate return direction: mij (2.61)
ij
B A
New plastic multiplier increment: A f Ap (2.65)
f Bp f Ap
G p ( ijt E ( i 1)
mij )
New return direction: (i )
mij (2.66)
ij
Each inequality depends on the output from the other one, therefore the following iterative
scheme is developed.
Algorithm 2:
Step 1: F p ( ( n 1) ij Eijkl ( kl (i 1) klf b ( i 1) klcor ( i ) klp )) 0 solve for (i ) klp
Iterative correction of the strain norm between two subsequent iterations can be expressed as
ATENA Theory 41
(i ) ijcor (1 b) f p ( i ) ijcor (2.74)
f p 1 , then the convergence would be too slow. In this case b can be estimated
f p
as b 1 , in order to increase the convergence rate.
1 f p , then the algorithm is diverging. In this case b should be calculated as
b 1 to stabilize the iterations.
pf
42
where the superscript i denotes values from two subsequent iterations. This will eliminate
problems due to the oscillation of the correction parameter b . Important condition for the
convergence of the above Algorithm 2 is that the failure surfaces of the two models are
intersecting each other in all possible positions even during the hardening or softening.
Additional constraints are used in the iterative algorithm. If the stress state at the end of the first
step violates the Rankine criterion, the order of the first two steps in Algorithm 2 is reversed.
Also in reality concrete crushing in one direction has an effect on the cracking in other
directions. It is assumed that after the plasticity yield criterion is violated, the tensile strength in
all material directions is set to zero.
On the structural level secant matrix is used in order to achieve a robust convergence during the
strain localization process.
The proposed algorithm for the combination of plastic and fracture models is graphically shown
in Fig. 2-23. When both surfaces are activated, the behavior is quite similar to the multi-surface
plasticity (SIMO et al. 1988). Contrary to the multi-surface plasticity algorithm the proposed
method is more general in the sense that it covers all loading regimes including physical changes
such as for instance crack closure. Currently, it is developed only for two interacting models, and
its extension to multiple models is not straightforward.
There are additional interactions between the two models that need to be considered in order to
properly describe the behavior of a concrete material:
(a) After concrete crushing the tensile strength should decrease as well
(b) According to the research work of Collins (VECHIO and COLLINS (1986)) and
coworkers it was established the also compressive strength should decrease when
cracking occurs in the perpendicular direction. This theory is called compression field
theory and it is used to explain the shear failure of concrete beams and walls.
The interaction (a) is resolved by adding the equivalent plastic strain to the maximal fracturing
strain in the fracture model to automatically increase the tensile damage based on the
compressive damage such that the fracturing strains satisfies the following condition:
ft p
ˆkk f eq (2.77)
f c
The compressive strength reduction (b) is based on the following formula based proposed by
Collins:
c rc f c
1
rc , rclim rc 1.0 (2.78)
0.8 170 1
Where 1 is the tensile strain in the crack. In ATENA the largest maximal fracturing strain is
used for 1 and the compressive strength reduction is limited by rclim . If rclim is not specified then
no compression reduction is considered.
ATENA Theory 43
CC3DNonLinCementitious,
CC3DNonLinCementitious2,
CC3DNonLinCementitious2Variable,
CC3DNonLinCementitious2Fatigue (described in section 2.2.10),
CC3DNonLinCementitious2User,
CC3DNONLINCEMENTITIOUS2SHCC (described in section 2.2.11),
and CC3DNonLinCementitious3 (described in section 2.2.12),
with the following differences: CC3DCementitious assumes linear response up to the point when
the failure envelope is reached both in tension and compression. This means that there is no
hardening regime in Fig. 2-21. The material CC3DNonLinCementitious on the contrary assumes
a hardening regime before the compressive strength is reached. The material
CC3DNonLinCementitious2 is equivalent to CC3DNonLinCementitious but purely incremental
formulation is used (in CC3DNonLinCementitious a total formulation is used for the fracturing
part of the model), therefore this material can be used in creep calculations or when it is
necessary to change material properties during the analysis. The material
CC3DNonLinCementitious2Variable is based on the material CC3DNonLinCementitious2 and it
allows to define history evolution laws for selected material parameters. The following material
parameters can be defined using an arbitrary evolution laws: young modulus E , tensile strength
f t ' , compressive strength f c' and f c'0 . It is the responsibility of the user to define the parameters
in a meaningful way. It means that at any time:
ft ' 1
2 f c'0 (2.79)
t/f’t
1.0
t t
(f1 ) Lt/Lch
loc
t
loc ~f1
Fig. 2-24. An example of a user defined tensile behavior for CC3DNonLinCementitious2User
material.
44
c/f’c
1.0
p c
(eqc ) Lc/Lch
loc
c
loc
~ peq
Fig. 2-25. An example of a user defined compressive behavior for CC3DNonLinCementitious2User
material.
G/Gc
1.0
f sh t
(1 ) Lt/Lch
loc
sh
loc
~f1
Fig. 2-26. An example of a user defined shear retention factor for shear stiffness degradation after
cracking.
In the user defined material mode II and III crack stiffness are evaluated with the help of the
shear retention factor rg as:
rg G
cr
Eijij (2.81)
1 rg
where i j , rg min(rgi , rgj ) is the minimum of shear retention factors on cracks in directions
i , j , and G is the elastic shear modulus. Shear retention factor on a crack in direction i is
evaluated from the user specified diagram as shown in Fig. 2-26.
In the above diagrams Lt and Lc represents the crack band size and crush band size respectively
as it is defined Section 2.1.3. Ltch and Lcch represents a size for which the tensile and compression
diagram respectively is valid. For instance it represents the measuring base that was used in an
experiment to determine the strain values in the diagrams above. loc represents the strain value,
after which strain localization can be expected. Usually, this is the strain after which the diagram
ATENA Theory 45
is entering into the softening regime. For instance, the strain value that is used to determine the
tensile strength is calculated based on the following assumptions:
if 1f loc
f
1f 1f
else
It should be realized that the compressive strength of the cracked concrete i.e. (2.83) is a
function of the maximal fracturing strain, i.e. maximal tensile damage at the given point. The
shear strength should be a function of the crack opening. Because of that the shear strength is
specified as a function of the fracturing strain 1f after the localization transformation (2.82).
The shear strength law is specified as a value relative to f t . The compressive strength reduction
is specified as a function relative to f c .
46
t/f’t
1.0
1.0 3/f’c
Fig. 2-27. An example of a user defined tensile strength degradation law due to lateral compressive
stress.
ft
c ts ft
Fig. 2-28: Tension stiffening.
ATENA Theory 47
overestimated. This is the consequence of the fact that the crack band approach assumes that the
crack spacing is larger than a finite element size. In heavily reinforced structures, or if large
finite elements are used, it may occur that the crack spacing will be smaller than finite element
size. This is especially true if shell/plate elements are used. In this case, typically large finite
elements can be used, and they usually contain significant reinforcement. In these cases, it is
useful to provide the crack spacing manually, since otherwise the program will overestimate the
cracking and due to that also larger deflections may be calculated. The program ATENA allows
the user to manually define the crack spacing. This user defined spacing is used as crack band
size Lt in cases when the user defined crack spacing is smaller than the Lt that would be
calculated by formulas presented in Section 2.1.3.
2.2.10 Fatigue
For modelling fatigue behavior of concrete (CEB 1988 and SAE AE-4) under tensile load, a new
material has been implemented in ATENA. The new material
(CC3DNonLinCementitious2Fatigue) is based on the existing three-dimensional fracture plastic
material (CC3DNonLinCementitious2) and uses a stress based model (2.2.10.1). It has an
additional parameter, fatigue , and additional data attributes for base , N , and fatigue , used in the
damage calculation as described in section 2.2.10.2. For details and validation against tests
conducted by KESSLER-KRAMER (2002) see ČERVENKA, PRYL (2007) or PRYL,
CERVENKA, PUKL (2010). Modelling 3-point bending tests with this material is presented in
PRYL, PUKL, CERVENKA (2013) and PRYL, D., MIKOLÁŠKOVÁ, J., PUKL, R. (2014).
48
min
where max is the maximum stress, f is static concrete strength, R , min is minimum
max
stress and is a material constant. The equation (2.86) holds for both compressive and tensile
stresses, however, the value of is not neccessarily the same for tensile and compressive
behavior of a material. The value should be determined from experiments. For example,
=0.052 was used based on the experimental results for load levels 0.7 and 0.9 Fstat when
modelling the test on a probe sealed during curing with a notch from section 3.5.2.4 of
KESSLER-KRAMER (2002) for validation.
Fig. 2-29: Typical S-N line for concrete in compression (KLAUSEN (1978))
The S-N relations mentioned above are mainly obtained by constant amplitude tests. However,
in real structures the stresses are varying. One method which can be of help in this context is the
well-known Palmgren-Miner hypothesis PALMGREN (1924), MINER (1945).
k
ni
N
i 1
1 (2.87)
i
where ni is the number of constant amplitude cycles at stress level i , N i is the number of cycles
to failure at stress level i , and k is the number of stress levels. As a rough tool this hypothesis is
useful, especially concerning steel. It can also be used for concrete although some investigations
have suggested that a value lower than 1 should be used.
ATENA Theory 49
2.2.10.2 Fatigue Damage Calculation
In the implemented model, fatigue damage consists of a contribution based on cyclic stress
(2.2.10.2.1), and an additional contribution from crack opening and closing in each cycle
(2.2.10.2.3). The former is dominant before cracking occurs, the latter in already cracked
regions.
Then, the damage due to fatigue after n cycles is calculated as an increase of the maximum
fracturing strain ˆij f (see section 2.2.3). The maximum fracturing strain in each principal
direction is adjusted by adding
w fatigue n
fatigue , where w fatigue w fail and the failing displacement for the given stress
ElemSize N
w fail invert _ soft _ law( upper ) (see Fig. 2-30).
50
In ATENA 4.0, a single value of fatigue is used to calculate fatigue damage caused by both
tensile and compressive stresses. So far, there is also no special provision implemented for loads
crossing zero, i.e., changing from tension to compression and back in each cycle, which lead to
faster damage according to experimental results presented in CEB 1988 and SAE AE-4. In that
situation, the damage is calculated separately for cyclic loading from 0 to max. compression and
from 0 to max. tension, and then the worse of the both damage values is considered. It should be
also noted that the damage is only introduced in form of maximum fracturing strain, which has
no direct impact on compressive material properties, i.e., the fatigue damage effectively only has
influence on tensile behaviour of the material.
wf1 + ((n_tot - N1) * (wf2 - wf1) / (N2 - N1)) - wf_curr for N1 <= n_tot < N2
wf2 + ((n_tot - N2) * (wfail - wf2) / (N - N2)) - wf_curr for N2 <= n_tot < N
wfail * n_tot / N - wf_curr for N <= n_tot
where
n_tot = n + N_beg, N_curr = N - N_beg,
N_beg = wf_curr * N1 / wf1 for wf_curr < wf1
N1 + (wf_curr - wf1) * (N2 - N1)/(wf2 - wf1) for wf1 <= wf_curr < wf2
N2 + (wf_curr - wf2) * (N - N2)/( wfail - wf2) for wf2 <= wf_curr
0 < N_beg < N
and
wfr_1 = 0.1, Nr_1 = 0.1, wfr_2 = 0.5, Nr_2 = 0.9, N1 = Nr_1 * N, N2 = Nr_2 * N.
1
Available since version 5.3.0
2
In ATENA versions prior to 5.1.3 and 5.3.4: w fatigueCOD n fatigue c fatigueCODload COD
ATENA Theory 51
2.2.10.3 Bringing in Fatigue Damage
It is recommended to introduce the fatigue induced damage into the unloaded structure (i.e., at
the lower stress level). Several other approaches of introducing the damage into the model were
also tested, i.e., introducing the damage at the upper load level or during reloading, but they
usually bring more convergence problems, especially during unloading.
ry c
rack
y crack set
set
r
prima
52
multiple or localized crack mc, lc and association with primary or secondary crack direction ,
)
b) localized cracking regime (softening)
A localized crack forms within a set of multiple cracks if the corresponding normal cracking
strain exceeds the level of mcmb (cracking strain capacity, a material constant).
Opening and sliding displacements of the i , i localized cracks are treated by the crack
lc, lc,
band model (i.e. they are transformed into cracking strains ij , ij by dividing them with
corresponding band width wc or wc).
The overall strain of the RVE is then obtained as a sum of strain of material between cracks
(which may possibly contain nonlinear plastic strain due to compressive yielding), cracking
strains due to multiple cracks, and cracking strains due to localized cracks:
s mc , mc , lc, lc ,
ij ij ij ij ij ij (2.88)
where ijs represents the strain of the continuous material between cracks.
unloading/ unloading/
reloading reloading
mcmb
ATENA Theory 53
Gc
. (2.89)
G
Let us determine stiffness Gc, while considering the most general 2-D case of an element, which
contains two perpendicular sets of multiple cracks and two perpendicular localized cracks. If the
problem is defined in plane , then the total engineering shear strain has only one non-zero
component, which is obtained as:
s mc , mc , lc , lc ,
2 2 2 2 2 , (2.90)
which can be rewritten with use of the shear bridging model (Kabele, 2000) as:
1 1 1 1 1 1
G M mc ,
M mc ,
wc L wc L Gc (2.91)
1 Vf k Gf
L() , for 0
2 2 0 4 k G f
2
1
3 E f d f
L ( ) 0 , for 0 (2.93)
Here Vf is the fiber volume fraction, Gf is the fiber shear modulus, Ef is the fiber Young’s
modulus, df is the fiber diameter, and k is the fiber cross-section shape correction factor. The
quantity and indicates the crack opening in direction and respectively. The parameter
0 represents the limiting value of the crack opening displacement, when no tensile stress can be
transferred across the crack, i.e. the point when the stress-displacement diagram in Fig. 2-32
drops to zero. These parameters are to be supplied by the user with the exception of the
parameter 0 , which is automatically extracted from the provided stress-strain law for tension.
The shear retention factor is then expressed as
1
(2.94)
1 1 1 1
1 G
M ( ) M ( ) wc L( ) wc L( )
mc , mc ,
Note that for an element containing only multiple cracks (before localization) 0 and
mc , mc ,
1/L terms approach zero. For an uncracked element, 0 and 1/M and 1/L
approach zero, giving =1.
54
such as confined reinforced concrete members (columns, bridge piers), nuclear vessels and
triaxial compression tests of plain concrete. A detailed description of the model formulation is
presented in PAPANIKOLAOU and KAPPOS (2007). In this section, only the main differences
between the CC3DNonLinCementitious3 and the CC3DNonLinCementitious2 model are
described, which are mainly focused on the plasticity part of the model (section 2.2.4).
Parameter t in equation (2.99) controls the slope of the softening function and the outmost square
is necessary due to the quadratic nature of the loading surface. The softening function value
starts from unity and complete material decohesion is attained at c = 0. The evolution of both
hardening and softening functions with respect to the hardening/softening parameter is
schematically shown in Fig. 2-33.
ATENA Theory 55
k(κ) / c(κ)
c k
1.0
k
0.8
0.6
c
0.4
0.2
ko
0.0
p
ε v,t κ = εpv
Fig. 2-33: Evolution of hardening (k) and softening (c) functions with respect to the plastic
volumetric strain.
56
Fig. 2-34: Direction (ψ) of the incremental (a) and total (b) plastic strain vectors.
fc (ΜPa) 20 30 40 50 60 70
Εc (MPa) 24377 27530 30011 32089 33893 35497
ν 0.2 0.2 0.2 0.2 0.2 0.2
ft (MPa) 1.917 2.446 2.906 3.323 3.707 4.066
λt 1.043 1.227 1.376 1.505 1.619 1.722
e 0.5281 0.5232 0.5198 0.5172 0.5151 0.5133
fco (MPa) -4.32 -9.16 -15.62 -23.63 -33.14 -44.11
ATENA Theory 57
fc (ΜPa) 80 90 100 110 120
Εc (MPa) 36948 38277 39506 40652 41727
ν 0.2 0.2 0.2 0.2 0.2
ft (MPa) 4.405 4.728 5.036 5.333 5.618
λt 1.816 1.904 1.986 2.063 2.136
e 0.5117 0.5104 0.5092 0.5081 0.5071
fco (MPa) -56.50 -70.30 -85.48 -102.01 -114.00
where J2 denotes the second invariant of stress deviator tensor. The parameter
k p
eq 1
3 y eqp is the maximal shear stress and y is the uniaxial yield stress. This
parameter controls the isotropic hardening of the yield criterion.
N inc
y eqp y H eqp , eqp 2
3 ε p
: ε p (2.102)
i 1
y is the yield stress, H the hardening modulus and eqp is the equivalent plastic strain
calculated as a summation of equivalent plastic strains during the loading history.
In case of von Mises plasticity the plastic potential function is identical with the yield function:
G p ( ij ) F P ( ij ) (2.103)
The associated flow rule is assumed. The background information can be found in (CHEN,
SALEEB 1982, Sec.5.4.2).
The Von Mises model could be used to model cyclic steel behavior including Bauschinger
effect. In this case the yield function is modified as:
58
1
2 σ X : σ X k eqp (r 1)k0 0 (2.104)
where σ is the deviatoric stress, k0 is an initial value of k ( eqp ) according to (2.102), X is the so
called back stress controlling the kinematic hardening:
X 2 3 k1 ε p k2 X eqp (2.105)
In equations (2.104) and (2.105) quantities r , k1 , k2 are material parameters for the cyclic
response. If r is non-zero the cyclic model is activated and it controls the radius of the Von
Mises surface. If r 1 the yielding will start exactly when y is reached. For lower values the
non-linear behavior starts earlier and the slope of the response is mainly affected by parameter
k1 (larger value – higher slope). Parameter k2 on the other hand affects the memory of the cyclic
response. Some examples of various parameter combinations are shown at Fig. 2-35.
ATENA Theory 59
Fig. 2-35: Effect of material parameter choice on cyclic response for E=210 GPa and y = 200 MPa.
60
2.4 Drucker-Prager Plasticity Model
Drucker-Prager plasticity model is based on a general plasticity formulation that is described in
Section 2.2.4. The yield function is defined as:
p
FDP ( ij ) I 1 J 2 k 0 (2.106)
Where and k are parameters defining the shape of the failure surface. They can be estimated
by matching with the Mohr-Coulomb surface. If the two surface are to agree along the
compressive meridian, i.e. 00 , the formulas are:
2sin 6 c cos
, k (2.107)
3 3 sin 3 3 sin
This corresponds to a outer cone to the Mohr-Coulomb surface. The inner cone, which passes
through the tensile meridian where 600 has the constants given by the following expressions:
2sin 6 c cos
, k (2.108)
3 3 sin 3 3 sin
The position of failure surfaces is not fixed but it can move depending on the value of strain
hardening/softening parameter. The strain hardening is based on the equivalent plastic strain,
which is calculated according to the following formula.
eqp min( ijp ) (2.109)
ATENA Theory 61
where determines the return direction. If 0 material is being compacted during crushing,
if 0 material volume is preserved, and if 0 material is dilating. In general the plastic
model is non-associated, since the plastic flow is not perpendicular to the failure surface
The return mapping algorithm for the plastic model is based on predictor-corrector approach as
is shown in Fig. 2-22. During the corrector phase of the algorithm the failure surface moves
along the hydrostatic axis to simulate hardening and softening. The final failure surface has the
apex located at the origin of the Haigh-Vestergaard coordinate system. Secant method based
Algorithm 1 is used to determine the stress on the surface, which satisfies the yield condition and
also the hardening/softening law.
c ,
2
c ft 2
0 1 , c , 0 ft
ft c c 2 ft
2 0
c2
1
ft c
2
0 , ft
62
In tension the failure criterion is replaced by an ellipsoid, which intersect the normal stress axis
at the value of ft with the vertical tangent and the shear axis is intersected at the value of c (i.e.
cohesion) with the tangent equivalent to .
The parameters for the interface model cannot be defined arbitrarily; there is certain dependence
of the some parameters on the others. When defining the interface parameters, the following
rules should be observed:
c
ft , ft c
(2.114)
c 0, f t 0, 0
It is recommended that parameters c, f t , are always greater than zero. In cases when no
cohesion or no tensile strength is required, some very small values should be prescribed.
Trial stress
1 c Initial surface
Final stress
Residual surface
ft
Fig. 2-37: Failure surface for interface elements.
In general three-dimensional case in Fig. 2-37 and equation (2.113) is calculated as:
12 22 (2.115)
ATENA Theory 63
c
min
Ktt
1
Ktt
1
(a) v
ft
min
Knn Knn
1
1
(b) u
Fig. 2-38: Typical interface model behavior in shear (a) and tension (b)
The K nn , K tt denote the initial elastic normal and shear stiffness respectively. Typically for zero
thickness interfaces, the value of these stiffnesses correspond to a high penalty number. It is
recommended not to use extremely high values as this may result in numerical instabilities. It is
recommended to estimate the stiffness value using the following formulas
E G
K nn , K tt (2.116)
t t
where E and G is minimal elastic modulus and shear modulus respectively of the surrounding
material. t is the width of the interface zone. Its value can be selected either on the basis of the
reality. For instance for mortar between masonry bricks the value is typically 10-20 mm.
Alternatively, it can be estimated as a dimension, which can be considered negligible with
respect to the structural size. For instance in case of a dam analysis, where the dam dimensions
are typically in the order of 100 meters, the width of the interface zone can be estimated to be 0.5
64
meters. It is suitable due to numerical reasons if stiffness is about 10 times of the stiffness of
adjacent finite elements.
There are two additional stiffness values that need to be specified in the ATENA input. They are
denoted in Fig. 2-38 as K nnmin and K ttmin . They are used only for numerical purposes after the
failure of the element in order to preserve the positive definiteness of the global system of
equations. Theoretically, after the interface failure the interface stiffness should be zero, which
would mean that the global stiffness will become indefinite. These minimal stiffnesses should be
about 0.001 times of the initial ones.
It is possible to define evolution laws for tensile as well as shear softening by arbitrary
multilinear laws. Examples of such laws are shown in Fig. 2-39. The figure describes bi-linear
softening laws. The break point of this law can be determined for instance by the formula
proposed by Bruehwiler and Wittman (1990).
ft GF
s1 , v1 0.75 (2.117)
4 ft
ft cc 0
I II
GF GF
s1 s1c
u1 u v 1c v
u eqf
u eqf
Fig. 2-39: Example of a softening law for tension and cohesion.
The evolution law depends on the equivalent nonlinear interface relative displacement
Where u f and v fi are the inelastic components of the relative interface displacement on the
basis of their decomposition into elastic and nonlinear, i.e. fracturing part.
u ue u f
(2.119)
vi vi v fi
This approach ensures that the degradation in shear affects also tensile strength and vice versa.
For instance, when the interface is damaged in shear, the tensile strength is reduced as well. The
typical behavior of the interface model with the softening evolution laws is shown in Fig. 2-38
by the dotted lines. The default behavior when no softening law is given is brittle with
immediate drop to zero in tension and to the residual dry friction in shear. The behavior is shown
in Fig. 2-38 by the solid black line.
When user softening laws are defined for the interface material, it is recommended that the
softening law for cohesion is always more ductile then the one for tensile strength, i.e. the
cohesion should be higher than the tensile strength at any time during the softening process.
ATENA Theory 65
Fig. 2-40: Example of a cyclic response of the model in shear under constant normal pre-stress.
2.7.1 Introduction
Reinforcement can be modeled in two distinct forms: discrete and smeared. Discrete
reinforcement is in form of reinforcing bars and is modeled by truss elements. The smeared
reinforcement is a component of composite material and can be considered either as a single
(only one-constituent) material in the element under consideration or as one of the more such
constituents. The former case can be a special mesh element (layer), while the later can be an
element with concrete containing one or more reinforcements. In both cases the state of uniaxial
stress is assumed and the same formulation of stress-strain law is used in all types of
reinforcement. More info about discrete reinforcement is available in Section 10.2.3 Discrete
Reinforcement Embedded in Solid Elements, located near the end of this manual.
66
Fig. 2-41 The bilinear stress-strain law for reinforcement.
The initial elastic part has the elastic modulus of steel Es. The second line represents the
plasticity of the steel with hardening and its slope is the hardening modulus Esh. In case of
perfect plasticity Esh =0. Limit strain L represents limited ductility of steel.
ATENA Theory 67
Fig. 2-43 Smeared reinforcement.
The spacing s of the smeared reinforcement is assumed infinitely small. The stress in the
smeared reinforcement is evaluated in the cracks, therefore it should include also a part of stress
due to tension stiffening (which is acting in concrete between the cracks, section 2.1.9).
scr ' s ' ts (2.120)
where s is the steel stress between the cracks (the steel stress in smeared reinforcement), scr
' '
is the steel stress in a crack. If no tension stiffening is specified ts =0 and scr s . In case of
' '
b
* * 1 b * , *
r
R R0
c1
(2.122)
,
1 0 r c2
1
*R
R
68
where R0 , c1 and c2 are experimentally determined parameters. The Fig. 2-44 shows the
meaning of strain values r , 0 , and stress values r and 0 . These values changes for each
cycle. The values with the subscript r indicate the point where the cycle started, and the
subscript 0 indicates the theorethical yield point that would be reached during the unloading if
the response would not have been modified by the hysteretic behavior. During the calculation of
this point the material stress-strain law is considered (see Sections 2.7.2, 2.7.3)
Nincr .
* f R eq , eq i
eq (2.123)
i 1
Fig. 2-44: Cyclic reinforcement model based on Menegotto and Pinto (1973).
ATENA Theory 69
2.8.1 CEB-FIP 1990 Model Code
b
s
b
max s
, 0 s s1 (2.124)
1
b max
, s1 s s 2 (2.125)
ss
b max f , s 2 s s3
2
(2.126)
s s
max
3 2
b f
, s3 s (2.127)
70
Table 2.8-1: Parameters for defining the mean bond strength-slip relationship for ribbed bars.
2 3 4 5
Value Unconfined concrete* Confined concrete**
Bond conditions Bond conditions
Good All other cases Good All other cases
S1 0.6 mm 0.6 mm 1.0 mm
S2 0.6 mm 0.6 mm 3.0 mm
S3 1.0 mm 2.5 mm clear rib spacing
0.4 0.4
max
2.0 f 1.0 f 2.5 f 1.25 f
C C C C
f
0.15 max 0.40 max
* Failure by splitting of the concrete
**Failure by shearing of the concrete between the ribs
Table 2.8-2: Parameters for defining the bond strength-slip relationship for smooth bars.
ATENA Theory 71
3 4
Table 2.8-3: Parameters for defining the bond strength-slip relationship for smooth bars.
72
2.8.3 Memory Bond Material
The Memory Bond material is an improvement to better capture the response during cyclic
loading and unloading in general. It can be used with any of the above mentioned bond stregth –
bond slip envelope functions. The response only differes after the bond stress sign changes.
Instead of following the same envelope as during loading, the maximum bond stress is
determined by the additional paramater 1 , see Fig. 2-47. Admissible values are res ≤ 1 ≤ max ,
where res is the residual bond stress (last value from the bond stregth – bond slip function) and
max the maximum bond stress (max. value from the bond strength – bond slip function).
In the figure, s is the current slip value, smax the maximum of the absolute slip value ever reached
(damage variable), f ( s ) is the bond strength function.
s 0 1
s 0 1
1 1
N ni n j ij , M
2
mi n j m j ni ij , L li n j l j ni ij
2
(2.128)
where mi and li are chosen orthogonal vectors lying in the microplane and defining the shear
strain components. The constitutive relations for the microplane strains and stresses can be
generally stated as:
N (t ) Ft0 N ( ), L ( ), M ( )
M (t ) Gt 0 N ( ), L ( ), M ( ) (2.129)
L (t ) Gt 0 N ( ), L ( ), M ( )
where F and G are functionals of the history of the microplane strains in time t. For a detailed
derivation of these functionals a reader is referred to BAZANT et al 2000 and CANER and
BAZANT 2000. The macroscopic stress tensor is obtained by the principle of virtual work that is
applied to a unit hemisphere . After the integration, the following expression for the
macroscopic stress tensor is recovered (BAZANT 1984):
3 Nm
ij sij d 6 w sij( ) , where sij N ni n j M mi n j m j ni L li n j l j ni (2.130)
2 1 2 2
where the integral is approximated by an optimal Gaussian integration formula for a spherical
surface; numbers label the points of the integration formula and w are the corresponding
optimal weights.
74
Inside each finite element at each integration point, an equivalent localization element is
assumed. The localization element is a serial arrangement of the localization zone, which is
loading, and an elastic zone (spring), which is unloading. The total length of the element is
equivalent to the crack band size L (width), and can be determined using the same methods as
described in Section 2.1.3 (see Fig. 2-12). The width of the localization zone is given either by
the characteristic length of the material or by the size of the test specimen for which the adopted
material model has been calibrated.
The three-dimensional equivalent element is constructed by three serial arrangements of the
elastic zone (spring) and localization band. The spring-band systems are perpendicular to each
other, and they are arranged parallel to the principal strain directions (Fig. 2-48). The simplified
two-dimensional version is shown in Fig. 2-49. In this arrangement of spring-band systems it is
possible to identify the following unknown stresses and strains:
ijb , 1 iju , 2 iju , 3 iju and ijb , 1 iju , 2 iju , 3 iju
m
where superscript b denotes the quantities in the localization band and the symbol x u with
superscripts u and m defines the quantities in the elastic spring in the direction m .
ij ,ij
u u
ij ,ij
b b
ij ,ij
u u
3
L
ij ,ij
3 u u
h
1
h 2
h
1
L 2
L
Fig. 2-48: The arrangement of the three-dimensional equivalent localization element.
ATENA Theory 75
Finite element Localization
element
1
Elastic
springs
1
ij
u
2
2
1
h ijb iju
Localization h
2 band
L 2
L
Fig. 2-49: The simplified two-dimensional view of the spring-band arrangement.
Ideally, the chosen directions should be perpendicular to the planes of failure propagation. In
ATENA, it is assumed for them to be aligned with the principal axes of the total macroscopic
strain tensor, which in most cases should approximately correspond to the above requirement.
Altogether there are 48 unknown variables. In the subsequent derivations, it is assumed that
these stresses and strains are defined in the principal frame of the total macroscopic strain tensor.
The set of equations available for determining these variables starts with the constitutive
formulae for the band and the elastic springs:
ijb F ( ijb ) (2.131)
The first formula (2.131) represents the evaluation of the non-linear material model, which in our
case is the microplane model for concrete. The second equation (2.132) is a set of three elastic
constitutive formulations for the three linear zones (springs) that are involved in the arrangement
at Fig. 2-48. This provides the first 24 equations, which can be used for the calculation of
unknown strains and stresses.
The second set of equations is provided by the kinematic constrains on the strain tensors.
76
1 b1 1 u 1
11
L1
11 h 11 L 1h
1
22 2 22b 2h 2 22u 2 L 2h
L
1 b 3
33
L 3
33 h 3 33u 3 L 3h
1 1 (2.133)
1
2 L
L
12 1 12b 1h 112u 1L 1h 2 12b 2h 212u 2
L 2h
1 1 1
23 2 23b 2h 2 23u 2 L 2h 3 23b 3h 3 23u 3 L 3h
2 L L
1 1 1
13 1 13b 1h 113u 1L 1h 3 13b 3h 313u 3 L 3h
2 L L
These 6 additional equations can be written symbolically as:
11 1
ij i ijb i h i iju i L i h j ijb j h j iju j L j h (2.134)
2 L L
The next set of equations is obtained by enforcing equilibrium in each direction between the
corresponding stress components in the elastic zone and in the localization band. For each
direction m , the following condition must be satisfied:
ijb m e j m iju m e j for m 1...3 (2.135)
where m e j denotes coordinates of a unit direction vector for principal strain direction m . Since
the principal frame of the total macroscopic strain tensor is used the unit vectors have the
following coordinates:
1
e j 1, 0, 0 , 2 e j 0,1, 0 , 3e j 0, 0,1 (2.136)
The remaining equations are obtained by enforcing equilibrium between tractions on the other
surfaces of the band and the elastic zone (layer) imagined as a spring:
ijb m e j n iju m e j where m 1..3, n 1...3, m n (2.137)
The equation (2.137) is equivalent to a static constraint on the remaining stress and strain
components of the elastic springs. Formulas (2.135) and (2.137) together with the assumption of
stress tensor symmetry represent the remaining 18 equations that are needed for the solution of
the three-dimensional equivalent localization element. These 18 equations can be written as:
ijb m iju for m 1...3 (2.138)
This means that the macroscopic stress must be equal to ijb , i.e., the stress in the localization
element, and that the stresses in all the three elastic zones must be equal to each other and to the
microplane stress ijb . This implies also the equivalence of all the three elastic strain tensors.
Based on the foregoing derivations, it is possible to formulate an algorithm for the calculation of
unknown quantities in the three-dimensional equivalent localization element.
Input: ij , ij , ijb , iju (2.139)
ATENA Theory 77
(i )
i
L jh j L ih
Step 1: d iju i j
Cijkl rkl(i 1) (2.141)
2L L
(i ) ( i 1) (i )
Step 2: iju iju d iju (2.142)
(i ) 2iL jL 2 i L j L i L jh j L ih u
Step 3: ijb ij ij (2.143)
i
L jh j L ih i
L j h j L ih
(i ) (i )
Step 4: rij(i ) ijb iju (2.144)
where Cijlk is the compliance tensor. The above iterative process is controlled by the following
convergence criteria;
d iju (i ) rij( i ) rij( i ) T d iju ( i )
e , e , e (2.145)
ij ijb ijb ij
The macroscopic stress is then equal to the stress in the localization band ijb . More details about
the derivations of the above algorithm as well as various examples of application can be obtained
from the original reference CERVENKA et al. 2004. It should be noted that the described
equivalent localization element is used only if the calculated crack band size L (see Section
2.1.3) in each principal strain direction is larger than the prescribed localization band size h . For
smaller element sizes the equivalent localization approach is not used and mesh-dependent
results may be obtained.
2.10 References
BASQUIN, H.O. (1910), The exponential law of endurance tests, Proc. ASTM, 10 (II).
BAZANT, Z.P, OH, B.H (1983) - Crack Band Theory for Fracture of Concrete, Materials and
Structures, RILEM, Vol. 16, 155-177.
BAŽANT, Z.P., (1984), ‘Microplane model for strain controlled inelastic behavior’, Chapter 3 in
Mechanics of Engineering Materials (Proc., Conf. held at U. of Arizona, Tucson, Jan. 1984),
C.S. Desai and R.H. Gallagher, eds., J. Willey, London, 45-59.
BAŽANT, Z.P., CANER, F.C., CAROL, I., ADLEY, M.D., AND AKERS, S.A., (2000),
‘Microplane Model M4 for Concrete: I. Formulation with Work-Conjugate Deviatoric Stress’, J.
of Engrg. Mechanics ASCE, 126 (9), 944-961.
BIGAL, A.J (1999) - Structural Dependence of Rotation Capacity of Plastic Hinges in RC
Beams and Slabs, PhD Thesis, Delft University of Technogy, ISBN 90-407-1926-8.
BRUEHWILER, E., and WITTMAN, F.H. (1990), “The Wedge Splitting Test, A New Method
of Performing Stable Fracture-Mechanics Tests”, Engineering Fracture Mechanics, Vol. 35, No.
1-3, pp. 117-125.
CANER, F.C., AND BAŽANT, Z.P., (2000) ‘Microplane Model M4 for Concrete: II. algorithm
and calibration.", J. of Engrg. Mechanics ASCE, 129 (9), 954-961.
CEB-FIP Model Code 1990, First Draft, Comittee Euro-International du Beton, Bulletin
d'information No. 195,196, Mars.
CEB 1988, Bulletin D’Information No 188, Fatigue of concrete structures, State of the art report.
78
CERVENKA, V., GERSTLE, K. (1972) - Inelastic Analysis of Reinforced Concrete Panels: (1)
Theory, (2) Experimental Verification and application, Publications IABSE, Zürich, V.31-00,
1971, pp.32-45, and V.32-II,1972, pp.26-39.
CERVENKA, V. (1985) - Constitutive Model for Cracked Reinforced Concrete, Journal ACI,
Proc. V.82, Nov-Dec., No.6,pp.877-882.
CERVENKA, V., PUKL, R., ELIGEHAUSEN, R. (1991) - Fracture Analysis of Concrete Plane
Stress Pull-out Tests, Proceedings, Fracture process in Brittle Disordered Materials, Noordwijk,
Holland, June 19-21.
CERVENKA, V., PUKL, R., OZBOLT, J., ELIGEHAUSEN, R. (1995), Mesh Sensitivity
Effects in Smeared Finite Element Analysis of Concrete Structures, Proc. FRAMCOS 2, 1995,
pp 1387-1396.
CERVENKA, V., PUKL, R. (1992) - Computer Models of Concrete Structures, Structural
Engineering International, Vol.2, No.2, May 1992. IABSE Zürich, Switzerland, ISSN 1016-
8664, pp.103-107.
CERVENKA, V., PUKL, R., OZBOLT, J., ELIGEHAUSEN, R. (1995) - Mesh Sensitivity
Effects in Smeared Finite Element Analysis of Concrete Fracture, Proceedings of FRAMCOS2,
Zurich, Aedificatio.
CERVENKA, V., CERVENKA, J. (1996) - Computer Simulation as a Design Tool for Concrete
Structures, ICCE-96, proceedings of The second International Conference in Civil Engineering
on Computer Applications Research and Practice, 6-8 April, Bahrain.
CERVENKA, J, CERVENKA, V., ELIGEHAUSEN, R. (1998), Fracture-Plastic Material Model
for Concrete, Application to Analysis of Powder Actuated Anchors, Proc. FRAMCOS 3, 1998,
pp 1107-1116.
ČERVENKA, J., BAŽANT Z.P., WIERER, M., (2004), `Equivalent Localization Element for
Crack Band Approach to Mesh Sensitivity in Microplane Model’, submitted for publication, Int.
J. for Num. Methods in Engineering.
ČERVENKA, J., PRYL, D., (2007), `Fatigue Modelling of Crack Growth by Finite Element
Method and Smeared Crack Approach’, Internal Report 2007-08-03-2002-DP, Cervenka
Consulting.
CRISFIELD, M.A., WILLS, J. (1989)- The Analysis of Reinforced Concrete Panels Using
Different Concrete Models, Jour. of Engng. Mech., ASCE, Vol 115, No 3, March, pp.578-597.
CRISFIELD, M.A. (1983) - An Arc-Length Method Including Line Search and Accelerations,
International Journal for Numerical Methods in Engineering, Vol.19,pp.1269-1289.
CHEN, W.F, SALEEB, A.F. (1982) - Constitutive Equations For Engineering Materials, John
Willey \& Sons, ISBN 0-471-09149-9.
DARWIN, D., PECKNOLD, D.A.W. (1974) - Inelastic Model for Cyclic Biaxial Loading of
Reinforced Concrete, Civil Engineering Studies, University of Illinois, July.
DE BORST, R. (1986), Non-linear analysis of frictional materials, Ph.D. Thesis, Delft
University of Technology, 1986.
DRUCKER, D.C., PRAGER, W., Soil Mechanics and Plastic Analysis or Limit Design, Q.
Appl. Math., 1952, 10(2), pp 157-165.
DYNGELAND, T. (1989) - Behavior of Reinforced Concrete Panels, Dissertation, Trondheim
University, Norway, BK-report 1989:1
ATENA Theory 79
FEENSTRA, P.H., Computational Aspects of Bi-axial Stress in Plain and Reinforced Concrete.
Ph.D. Thesis, Delft University of Technology, 1993.
FEENSTRA, P.H., ROTS, J.G., AMESEN, A., TEIGEN, J.G., HOISETH, K.V., A 3D
Constitutive Model for Concrete Based on Co-rotational concept. Proc. EURO-C 1998, 1, pp.
13-22.
ETSE, G., Theoretische und numerische Untersuchung zum diffusen und lokalisierten Versagen
in Beton, Ph.D. Thesis, University of Karlsruhe 1992.
FELIPPA, C. (1966) - Refined Finite Element Analysis of Linear and Nonlinear Two-
Dimensional Structures, Ph.D. Dissertation, University of California, Engineering, pp.41-50.
GRASSL, P., LUNDGREN, K., and GYLLTOFT, K. (2002) “Concrete in compression : A
plasticity theory with a novel hardening law”, International Journal of Solids and Structures,
39(20), 5205-5223.
VAN GYSEL, A., and TAERWE, L. (1996) “Analytical formulation of the complete stress-
strain curve for high strength concrete”, Materials and Structures, RILEM, 29(193), 529-533.
HARTL, G. (1977) “Die Arbeitlinie Eingebetete Staehle bei erst und kurz=Belastung”,
Dissrtation, Univbersitaet Innsbruck
HORDIJK, D.A. (1991) - Local Approach to Fatigue of Concrete, Doctor dissertation, Delft
University of Technology, The Netherlands, ISBN 90/9004519-8.
KABELE, P. (2002) - Equivalent Continuum Model of Multiple Cracking, Engineering
Mechanics 2002, 9 (1/2), pp.75-90, Assoc.for Engineering Mechanics, Czech Republic
KESSLER-KRAMER, CH., (2002) “Zugverhalten von Beton unter Ermüdungsbeanspruchung”,
Schriftenreihe des Instituts für Massivbau und Baustofftechnologie, Heft 49, Karlsruhe.
KLAUSEN, D. (1978), Festigkeit und Schadigung von Beton bei haufig wiederholter
Beanschpruchung, PhD Thesis, University of Technology Darmstadt, 85 pp.
KOLLEGGER, J. - MEHLHORN, G. (1988) - Experimentelle und Analytische Untersuchungen
zur Aufstellung eines Materialmodels für Gerissene Stahbetonscheiben, Nr.6 Forschungsbericht,
Massivbau, Gesamthochschule Kassel.
KOLMAR, W. (1986) - Beschreibung der Kraftuebertragung über Risse in nichtlinearen Finite-
Element-Berechnungen von Stahlbetontragwerken", Dissertation, T.H. Darmstadt, p. 94.
KUPFER, H., HILSDORF, H.K., RÜSCH, H. (1969) - Behavior of Concrete under Biaxial
Stress, Journal ACI, Proc. V.66,No.8, Aug., pp.656-666.
MARGOLDOVA, J., CERVENKA V., PUKL R. (1998), Applied Brittle Analysis, Concrete
Eng. International, November/December 1998.
MENETREY, P., WILLAM, K.J. (1995), Triaxial failure criterion for concrete and its
generalization. ACI, Structural Journal, 1995, 92(3), pp 311-318.
MENETREY, Ph., WALTHER, R., ZIMMERMAN, Th., WILLAM, K.J., REGAN, P.E.
Simulation of punching failure in reinforced concrete structures. Journal of Structural
Engineering, 1997, 123(5), pp 652-659.
MIER J.G.M van (1986) - Multiaxial Strain-softening of Concrete, Part I: fracture, Materials
and Structures, RILEM, Vol. 19, No.111.
MINER M.A. (1945), Cumulative damage in fatigue. Transactions of the American Society of
Mechanical Engineering, 67:A159-A164.
80
OLIVIER, J., A Consistent Characteristic Length For Smeared Cracking Models, Int. J. Num.
Meth. Eng., 1989, 28, pp 461-474.
OWEN, J.M., FIGUEIRAS, J.A., DAMJANIC, F., Finite Element Analysis of Reinforced and
Pre-stressed concrete structures including thermal loading, Comp. Meth. Appl. Mech. Eng., 1983,
41, pp 323-366.
PALMGREN, A. (1924), Die Lebensdauer von Kugellagern. Zeitschrift Verein Deutscher
Ingenieure, 68(14):339-341.
PAPANIKOLAOU, V.K., and KAPPOS, A.J. (2007) “Confinement-sensitive plasticity
constitutive model for concrete in triaxial compression”, International Journal of Solids and
Structures, 44(21), 7021-7048.
PRAMONO, E, WILLAM, K.J., Fracture Energy-Based Plasticity Formulation of Plain
Concrete, ASCE-JEM, 1989, 115, pp 1183-1204.
PRYL, D., CERVENKA, J., and PUKL, R. (2010) “Material model for finite element modelling
of fatigue crack growth in concrete”, Procedia Engineering, 2 (2010) 203–212.
PRYL, D., PUKL, R., and CERVENKA, J. (2013) “Modelling high-cycle fatigue of concrete
specimens in three point bending”, Life-Cycle and Sustainability of Civil Infrastructure Systems
(Eds. Strauss, Frangopol & Bergmeister) 1303–1306.
PRYL, D., MIKOLÁŠKOVÁ, J., PUKL, R. (2014) “Modeling Fatigue Damage of Concrete”,
Key Engineering Materials, 2014 Vols. 577-578, pp. 385-388, ISSN: 1662-9795
RAMM, E. (1981) - Strategies for Tracing Non- linear Responses Near Limit Points, Non- linear
Finite Element Analysis in Structural Mechanics, (Eds. W.Wunderlich, E.Stein, K.J.Bathe)
RASHID, Y.R. (1968), Ultimate Strength Analysis of Pre-stressed Concrete Pressure Vessels,
Nuclear Engineering and Design,1968, 7, pp 334-344.
ROTS, J.G., BLAAUWENDRAAD, J., Crack models for concrete: discrete or smeared? Fixed,
multi-directional or rotating? HERON 1989, 34(1).
SAE, AE-4, Fatigue Design Handbook
SIMO, J.C., JU, J.W., Strain and Stress-based Continuum Damage Models-I. Formulations, II-
Computational Aspects, Int. J. Solids Structures, 1987, 23(7), pp 821-869.
SIMO, J.C., KENNEDY, J.G., GOVINDJEE, S., (1988), Non-smooth Multi-surface Plasticity
and Visco-plasticity. Loading/unloading Conditions and Numerical Algorithms, Int. J. Num.
Meth. Eng., 26, pp 2161–2185.
TAYLOR, G.I., (1938), ‘Plastic strain in metal’, J. Inst. Metals, 62, 307-324.
VAN MIER J.G.M. (1986), Multi-axial Strain-softening of Concrete, Part I: fracture, Materials
and Structures, RILEM, 1986, 19(111).
VECCHIO, F.J., COLLINS, M.P (1986)- Modified Compression-Field Theory for Reinforced
Concrete Beams Subjected to Shear, ACI Journal, Proc. V.83, No.2, Mar-Apr., pp 219-231.
VOS, E. (1983) - Influence of Loading Rate and Radial Pressure on Bond in Reinforced
Concrete, Dissertation, Delft University, pp.219-220.
WILKINS, M.L., Calculation of Elastic-Plastic Flow, Methods of Computational Physics, 3,
Academic Press, New York, 1964.
ATENA Theory 81
3 FINITE ELEMENTS
3.1 Introduction
The preceding chapters dealt with the general formulation of the problem, geometric and
constitutive equations. All expressions were derived independently of the structural shape, the
finite elements used etc. Here, an information about finite elements currently implemented in
ATENA is given.
t
h2
3
7 6 h1
4 9 2
5 s
8
r 1
t t
3 3 h6
h8 7 6
7 6
2 4 9 2
4 9
5 5 s
8 s 8
1 r 1
r
t
h9
3
7 6
4 9 2
5 s
8
r 1
ATENA Theory 83
is not because of its superior properties, but due to the fact that it is a versatile and general
approach with no hidden difficulties and, also very important, these elements are easy to
understand. This is very important particularly in nonlinear analysis. For example it is highly
undesirable to add element-related problems to problems related to e.g. material modeling.
Big advantage of ATENA isoparametric elements is that their interpolation functions hi ( r, s, t )
are constructed in hierarchical manner. Take an example of plane quadrilateral element. Some of
its interpolation functions are depicted in Fig. 3-1. The 1st four functions, i.e. functions h1 ( r, s, t )
to h4 ( r, s, t ) has to be always present in the interpolation set, (to ensure bilinear approximation).
Then, any additional function h6 ( r, s, t ) through h9 ( r, s, t ) can be added independently. This
would involve adding the new function itself and also amendments to the already present
interpolation functions. This approach (and use of C++ templates) makes possible that one
element formulation generates quadrilateral elements with nodes (1,2,3,4), (1,2,3,4,5),
(1,2,3,4,6), ... (1,2,3,4,8), (1,2,3,4,9), (1,2,3,4,5,6), (1,2,3,4,5, 7), ... (1,2,3,8,9), ...
(1,2,3,4,5,6,7,8,9). Additional mid-side points are particularly useful for changing mesh density,
(i.e. element size), see Fig. 3-2, as they allow change of mesh density without need triangular
elements.
Although the concept of hierarchical elements was described for plane quadrilateral elements, in
ATENA it applies for plane triangular elements, 3D bricks, tetrahedral and wedge elements, too.
Always there is a set of basic interpolation function that can be extended by any “higher”
interpolation function.
Apart of interpolation functions finite element properties depend strongly on numerical
integration scheme used to integrate element stiffness matrix, element nodal forces etc. In Atena,
majority of elements are integrated by Gauss integration scheme that ensure n ( n 1) order
accuracy, where n is degree of the polynomial used to approximate the integrated function.
9 10 11 9 10 11
6 7 8 7
6 8
3 4 5 3 4 5
1 2 1 2
84
3.2 Truss 2D and 3D Element
2D and 3D truss elements in ATENA are coded in group of elements CCIsoTruss<xx> ...
CCIsoTruss<xxx>. The string in < > describes present element nodes, (see Atena Input File
Format document for more information). These are isoparametric elements integrated by Gauss
integration at 1 or 2 integration points for the case of linear or quadratic interpolation, i.e. for
elements with 2 or 3 element nodes, respectively. They are suitable for plane 2D as well as 3D
analysis problems. Geometry, interpolation functions and integration points of the elements are
given in Fig. 3-3, Table 3.2-1 to Table 3.2-3.
s
y
2 r
3 CCIsoTruss<xx>
1
CCIsoTruxx<xxx>
x
1 1 1
(1 r ) h3
2 2
2 1 1
(1 r ) h3
2 2
3 (1 r 2 )
Table 3.2-2 Sample points for Gauss integration of 1 node CCIsoTruss<xx> element.
1 0. 2.
ATENA Theory 85
Table 3.2-3 Sample points for Gauss integration of 2 and 3 nodes CCIsotruss<xxx> elements.
1 0.577350269189626 1.
2 -0.577350269189626 1.
The element vectors and matrices for Total Lagrangian formulation , configuration at time t
and iteration (i) are as follows. Note that they are equally applicable for Updated Lagrangian
formulation upon applying changes related to the element reference coordinate system
(undeformed vs. deformed element axis.). The formulation is present for 3-nodes element option.
The 2-nodes variant is obtained by simply neglecting the terms for the element mid-point.
An arbitrary point on the truss element has at reference time t coordinates t X [t x1 , t x1 , t x1 ] :
t
x1 t x11 h1 t x12 h2 t x13 h3
t
x2 t x21 h1 t x22 h2 t x23 h3 (3.1)
t
x3 t x31 h1 t x32 h2 t x33 h3
t t
At time t t ( i 1) the same point has coordinates X ( i 1) :
t t
x1( i 1) ( t x11 t u11( i 1) ) h1 ( t x12 t u12( i 1) )h2 ( t x13 t u13( i 1) )h3
t t
x2( i 1) ( t x21 t u21( i 1) )h1 ( t x22 t u22( i 1) ) h2 ( t x23 t u23( i 1) )h3 (3.2)
t t
x3( i 1) ( t x31 t u31( i 1) )h1 ( t x32 t u32( i 1) ) h2 ( t x33 t u33( i 1) )h3
t t
and at time t t ( i ) coordinates X (i )
t t
x1( i ) ( t x11 t u11( i ) ) h1 ( t x12 t u12( i ) )h2 ( t x13 t u13( i ) )h3
t t
x2( i ) ( t x21 t u21( i ) )h1 ( t x22 t u22( i ) ) h2 ( t x23 t u23( i ) )h3 (3.3)
t t
x3( i ) ( t x31 t u31( i ) )h1 ( t x32 t u32( i ) ) h2 ( t x33 t u33( i ) )h3
t t ( i ) t t
Increment of Green Lagrange strain t 11( i )
t 11 t tt11( i 1) (at time t t , iteration (i )
with
to configuration at time t ) is calculated:
86
t t ( i ) 2 t t ( i 1) 2
l l
r r
(i ) 1
(3.4)
t 11 2
2 tl
r
where truss length differentials are
2 2 2 2
t l t x1 t x2 t x3
r r r r
2 2 2 2
t t l ( i 1) t t x1( i 1) t t x2( i 1) t t x3( i 1)
(3.5)
r r r r
2 2 2 2
t t l ( i ) t t x1( i ) t t x2( i ) t t x3( i )
r r r r
Substituting (3.5), (3.3) into (3.4) after some math manipulation it can be derived:
h1 h1 t 1 h1 h2 h1 h3 t 3
r r x1 r r x12
t
x1
r r
h2 h1 t x1 h2 h2 x1
t 2 h2 h3 t 3
x1
r r 1 r r r r
h h h h h h
3 1 t x11 3 2 x1 3 3 t x13
t 2
r r r r r r
h1 h1 t 1 h1 h2 h h
x2 x2 1 3 t x23
t 2
r r r r r r
t t 1 h2 h1 t 1 h2 h2 h h
t BL 0 x2 x2 2 3 t x23
t 2
2
t l r r r r r r
r h3 h1 t 1 h3 h2 h3 h3 t 3
x x2
t 2
x2
r r 2 r r r r
h h h h h h
1 1 t x31 1 2 x3 1 3 t x33
t 2
r r r r r r
h2 h1 t 1 h2 h2 h h
x3 x3 2 3 t x33
t 2
r r r r r r
h3 h1 t 1 h3 h2 h h
r r x3 r r x3 3 3 t x33 (3.6)
t 2
r r
ATENA Theory 87
h1 h1 t t 1( i 1) h1 h2 t t h1 h3 t t 3( i 1)
r r u1 u12( i 1) u1
r r r r
2 1 t t u1( i 1) h2 h2
h h t t 2( i 1)
u1
h2 h3 t t 3( i 1)
u1
r r 1
r r r r
h h h h h h
3 1 t t u11( i 1) 3 2 t t 2( i 1)
u1 3 3 t t u13( i 1)
r r r r r r
h1 h1 t t 1( i 1) h1 h2 t t 2( i 1) h h
u2 u2 1 3 t t u23( i 1)
r r r r r r
t t 1 h2 h1 t t 1( i 1) h2 h2 h2 h3 t t 3( i 1)
t BL( i11) 2 u2 t t 2( i 1)
u2 u2
t l r r r r r r
r h3 h1 t t 1( i 1) h3 h2 h h
u2 t t 2( i 1)
u2 3 3 t t u23( i 1)
r r r r r r
h h h h h1 h3 t t 3( i 1)
1 1 t t u31( i 1) 1 2 t t 2( i 1)
u3 u3
r r r r r r
h2 h1 t t 1( i 1) h2 h2 t t 2( i 1) h h
u3 u3 2 3 t t u33( i 1)
r r r r r r
h3 h1 t t 1( i 1) h3 h2 t t 2( i 1) h3 h3 t t 3( i 1)
r r u3 u3 u3 (3.7)
r r r r
and
h1 h2 h3
r 0 0 0 0 0 0
r r
1 h1 h2 h3
t t ( n 1)
BNL t 0 0 0 0 0 0 (3.8)
l
t
r r r
r h1 h2 h3
0 0 0 0 0 0
r r r
r
Note that 2-nodes truss element has constant strains along its length and thus the increment of
Green Lagrange strain can be calculated directly, (i.e. not using differentials truss length as it
was the case of (3.4) ):
88
1 l l
t t ( i ) 2 t t ( i 1) 2
t 11( i ) (3.11)
2 t 2
l
This yields a bit simpler element formulation (with the same results). However, for the sake of
preserving unified approach to all truss elements, ATENA uses even in this case the equation
(3.4).
y 2 s
5
1 CCIsoQuad<xxxx>
6
CCIsoQuad<xxxxx>
9 CCIsoQuad<xxxx_x>
8 r ....
3 CCIsoQuad<xxxx_x_x_>
7 ....
4 x CCIsoQuad<xxxxxxxxx>
ATENA Theory 89
Table 3.3-1: Interpolation functions of CCIsoQuad<...> elements.
Table 3.3-2: Sample points for Gauss integration of 4 nodes CCIsoQuad<...> element.
1 0.577350269189626 0.577350269189626 1.
2 0.577350269189626 -0.577350269189626 1.
3 -0.577350269189626 0.577350269189626 1.
4 -0.577350269189626 -0.577350269189626 1.
90
Table 3.3-3: Sample points for Gauss integration 5 to 9 nodes CCIsoQuad<...> elements.
(i )
t 11 t u1,1
(i )
t tt u1,1
( i 1)
t u1,1
(i ) t t ( i 1)
t u2,1 t u2,1
(i ) 1
2
ut 1,1 u
(i ) 2
t
(i ) 2
2,1
t 22( i ) t u2,2
(i )
t tt u1,2
( i 1)
t u1,2
(i ) t t ( i 1)
t u2,2 t u2,2
(i ) 1
2
ut 1,2 u
(i ) 2
t
(i ) 2
2,2
1
(i )
t 12
2
t u1,2 t u2,1
(i ) (i )
2
u(i ) t t
u1( i ) u1( i ) 1 u (i )
(i )
t1 t1
x
t 33 2
x1 t 2 x1 (3.12)
1
ATENA Theory 91
Displacement derivatives:
t t ui( i ) t t ui( i 1)
t
(i )
u
i, j
t x j
(3.13)
t t ( i 1) t t ui( i 1)
t i, j u
t x j
Strains and matrices to calculate them:
t (i ) t t
t BL( i 1) U ( i )
hi
t hi , j
t x j
(3.17)
ui( i ) t t ui( i ) t t ui( i 1)
n
t
x1 hk t x1k
k 1
92
Linear strain-displacement matrix – non-constant part:
t t ( i 1)
l11 t h1,1 t t ( i 1)
l
21 t 1,1 h t t ( i 1)
l
11 t 2,1 h t t ( i 1)
l h
21 t 2,1
t t ( i 1) t t ( i 1) t t ( i 1) t t ( i 1)
l11 t h1,2 l
21 t 1,2 h l
11 t 2,2 h l h
21 t 2,2
t t
t BL( i11) l11 t h1,2 t t l11( i 1) t h1,1
t t ( i 1) t t ( i 1)
l
21 t 1,2 h t t ( i 1)
21 t 1,1 l h t t ( i 1)
l
11 t 2,2 h t t ( i 1)
11 t 2,1l h t t ( i 1)
l h
21 t 2,2 t t ( i 1)
l h
21 t 2,1
t t ( i 1) h1 t t ( i 1) h
l33 t 0 l 2
0
x1 x
33 t
1
x
33 t
1
where
n
t t ( i 1)
11l t hk ,1 t tt u1k ( i 1)
k 1
n
t t ( i 1)
12l t hk ,2 t tt u1k ( i 1)
k 1
n
t t ( i 1)
21l t hk ,1 t tt u2k ( i 1) (3.19)
k 1
n
t t ( i 1)
22l t hk ,2 t tt u2k ( i 1)
k 1
1 n
t
t t ( i 1) t t k ( i 1)
33l hk t 1 u
x1 k 1
t t ( i 1)
0 t h1,1 0 t h2,1 ... 0 t hn ,1
t BNL (3.20)
0 t h1,2 0 t h2,2 ... 0 t hn ,2
h1 h2 h
t 0 t
0 ... t n 0
x1 x1 x1
ATENA Theory 93
2nd Piola-Kirchhoff stress tensor and vector
t tt S11(i 1) t t ( i 1)
t 12 S 0 0 0
t t (i 1) t t ( i 1)
t S 21 t 22 S 0 0 0
t S ( i 1)
0 0 t t ( i 1)
t 11 S t t ( i 1)
t 12 S 0
t t ( i 1) t t ( i 1)
0 0 t 21 S t 22 S 0
0 0 0 0 t t ( i 1)
t S33
(3.21)
t S ( i ) t tt S11(i 1) t t
t
( i 1)
S 22 t t
t
( i 1)
S 21 t t
t S33(i 1)
In case of the simplified 3D analysis, i.e. elements CCIsoQuad2_5<...>, the equations are further
extended as follows:
All element matrices and vectors are computed with respect to element local coordinate
system xlocal ,1 , xlocal ,2 using equations in (3.12) through (3.21). They are transformed into
3D global coordinate system by means of simple transformation:
M global T M local T T , vglobal T vlocal (3.22)
where
M global , M local , vglobal , vlocal are global and local finite element matrices and vectors,
The local element coordinate system (see Fig. 3-5) is defined by local xlocal ,1 , xlocal ,2 , xlocal ,3
coordinates. All of them pass through origin of the global (reference) coordinate system. The
axes xlocal ,1 and xlocal ,2 constitute a local coordinates element plane that is parallel to the element..
The axis xlocal ,3 is perpendicular to the element and the axis xlocal ,1 is defined as a projection of
global x1 axis to the local coordinate element plane. An exception to that is, when the element is
normal to the global x1 . In this case the local xlocal ,1 coincides with the global x2 axis.
The present definition of local element coordinate system depends on plane of the finite element
but it does not depend on its shape itself. This is very important property, as ATENA supports
use of local (instead of global) nodal degrees of freedom and, (of course) these degrees of
freedom must refer to a coordinate system common to all elements of the plane, in which they
lie.
94
x3
3
xl o c a l , 3
4 O
X
X’ 2
xl o c a l , 2
1 x2
x1 xl o c a l , 1
ATENA Theory 95
2D, axisymmetric and 3D problems. Geometry, interpolation functions and integration points of
3 s
CCIsoTriangle<xxx>
...
6 CCIsoTriangle<xxxxxx>
y
5
1
4 r
2
3 s
CCIsoTriangle<xxx>
...
6 CCIsoTriangle<xxxxxx>
y
5
1
4 r
2
x
Fig. 3-6: Geometry of CCIsoTriangle<...> elements.
96
4 4 r (1 r s )
5 4rs
6 4 s (1 r s )
Table 3-2: Sample point for Gauss integration of 3 nodes CCIsoTriangle<...> elements.
ATENA Theory 97
t
4
z CCIsoTetra<xxxx>
9 .............
10 s
3 CCIsoTetra<xxxxxxxxxx>
8
7 6
1
5 r
2 y
x
t
Fig. 3-7 Geometry of CCIsoTetra<...> elements.
t
3
11 10
4 2
12 9 CCIsoBrick<xxxxxxxx>
1
..................
z20 19
CCIsoBrick<xxxxxxxxx...x>
7 18
15 s
8 17 14
6
r 16
13
5 y
x
t
9 3
1
8
7
z 2 15 CCIsoWedge<xxxxxx>
.............
13 14 s
6 CCIsoWedge<xxxxxxxxxx xxxxx>
12 11
4
10 r
5 y
x
98
Table 3.5-1 Interpolation functions of CCIsoTetra<...> elements.
6 4 rs (1 t )
7 4 s (1 r s t )
8 4 rt (1 s )
9 4 s t (1- r )
10 4t (1- r - s - t )
Table 3.5-2 Sample point for Gauss integration of 4 nodes CCIsoTetra<...> element.
ATENA Theory 99
Table 3.5-3 Sample points for Gauss integration of 5 to 10 nodes CCIsoTetra<...> elements.
i=9 i = 10 i = 11 i = 12 i = 13 i = 14 i = 15 i = 16 i = 17 i = 18 i = 19 i = 20
1 1 1
(1 r )(1 s)(1 t ) h 9
1
h12
1
h17
8 2 2 2
2 1
(1 r )(1 s )(1 t )
1 1
h9 h10
1
h18
8 2 2 2
3 1
(1 r )(1 s )(1 t )
1 1
h10 h11
1
h19
8 2 2 2
4 1
(1 r )(1 s )(1 t )
1 1
h11 h12
1
h20
8 2 2 2
5 1
(1 r )(1 s )(1 t )
1
h13
1 1
h16 h17
8 2 2 2
6 1
(1 r )(1 s )(1 t )
1 1
h13 h14
1
h18
8 2 2 2
7 1
(1 r )(1 s )(1 t )
1 1
h14 h15
1
h19
8 2 2 2
8 1
(1 r )(1 s )(1 t )
1 1
h15 h16
1
h20
8 2 2 2
100
9 1
(1 r 2 )(1 s )(1 t )
4
10 1
(1 r )(1 s 2 )(1 t )
4
11 1
(1 r 2 )(1 s )(1 t )
4
12 1
(1 r )(1 s 2 )(1 t )
4
13 1
(1 r 2 )(1 s )(1 t )
4
14 1
(1 r )(1 s 2 )(1 t )
4
15 1
(1 r 2 )(1 s )(1 t )
4
16 1
(1 r )(1 s 2 )(1 t )
4
17 1
(1 r )(1 s )(1 t 2 )
4
18 1
(1 r )(1 s )(1 t 2 )
8
19 1
(1 r )(1 s )(1 t 2 )
4
20 1
(1 r )(1 s )(1 t 2 )
4
Table 3.5-5 Sample points for Gauss integration of 9 to 20 nodes CCIsoBrick<...> element.
102
4 0.7745966692414 0. 0.774596669241483 0.2743484225
83
5 0.7745966692414 0. 0. 0.4389574760
83
6 0.7745966692414 0. - 0.2743484225
83 0.774596669241483
7 0.7745966692414 - 0.774596669241483 0.1714677641
83 0.7745966692414
83
8 0.7745966692414 - 0. 0.2743484225
83 0.7745966692414
83
10 0. 0.7745966692414 0.774596669241483 0.2743484225
83
11 0. 0.7745966692414 0. 0.4389574760
83
12 0. 0.7745966692414 - 0.2743484225
83 0.774596669241483
13 0. 0. 0.774596669241483 0.4389574760
14 0. 0. 0. 0.7023319616
15 0. 0. - 0.4389574760
0.774596669241483
16 0. - 0.774596669241483 0.2743484225
0.7745966692414
83
17 0. - 0. 0.4389574760
0.7745966692414
83
18 0. - - 0.2743484225
0.7745966692414 0.774596669241483
83
19 - 0.7745966692414 0.774596669241483 0.1714677641
0.7745966692414 83
83
20 - 0.7745966692414 0. 0.2743484225
0.7745966692414 83
83
104
Function Include only if node i is defined
Node I
1 hh1 hv1 1
h7
1
h9
1
h13
2 2 2
2 hh2 hv1 1
h7
1
h8
1
h14
2 2 2
3 hh3 hv1 1
h8
1
h9
1
h15
2 2 2
4 hh1 hv2 1
h10
1
h12
1
h13
2 2 2
5 hh2 hv2 1
h10
1
h11
1
h14
2 2 2
6 hh3 hv2 1
h11
1
h12
1
h15
2 2 2
7 hh4 hv1
8 hh5 hv1
9 hh6 hv1
10 hh4 hv2
11 hh5 hv2
12 hh6 hv2
13 hh1 hv3
14 hh2 hv3
15 hh3 hv3
Table 3.5-7 Sample points for Gauss integration of 6 nodes CCIsoWedge<...> element.
2
t uk ,i t uk , j t uk , j t uk ,i
2
t uk ,i t uk , j
(i ) (i )
(3.25)
t x j
(3.26)
t t ( i 1)
t t ( i 1) u
u
t i, j i
xj t
106
t h1,1 0 0 t h2,1 0 0 ... t hn ,1 0 0
0 h 0 0 h2,2 0 ... 0 hn ,2 0
t 1,2 t t
t t
0 0 h
t 1,3 0 0 t h2,3 ... 0 0 t n ,3
h
B
t L0 (3.29)
t h1,2 h
t 1,1 0 t h2,2 t h2,1 0 ... t hn ,2 t hn ,1 0
0 t h1,3 t h1,2 0 t h2,3 t h2,2 ... 0 t hn ,3 t hn ,2
t h1,3 0 h
t 1,1 t h2,3 0 t h2,1 ... t hn ,3 0 t hn ,1
where
hi
t hi , j
t x j (3.30)
t t
t t ( i 1)
l t h1,3 t t ( i 1)
l h t t ( i 1)
l h t t ( i 1)
l h
t BL( i11) t t ( i 1) 13 t t ( i 1) t t ( i 1)
23 t 1,3
t t ( i 1) t t ( i 1)
33 t 1,3
t t ( i 1) t t ( i 1)
13 t 2,3
t t ( i 1)
l11 t h1,2 l11 t h1,1 l21 t 1,2 h 21 t 1,1l h l
31 t 1,2h l h
31 t 1,1 l 11 t 2,2h l h
12 t 2,1
t t l12( i 1) t h1,3 t t l13( i 1) h1,2 t t ( i 1)
l22 t 1,3 h t t ( i 1)
23 t 1,2l h t t ( i 1)
l
32 t 1,3h t t ( i 1)
l h
33 t 1,2
t t ( i 1)
l 12 t 2,3h t t ( i 1)
l h
13 t 2,2
t t ( i 1) t t ( i 1)
t
t t ( i 1) t t ( i 1) t t ( i 1) t t ( i 1) t t ( i 1) t t ( i 1)
l11 t 1,3 h l13 t h1,1 l21 t 1,3 h 23 t 1,1l h l
31 t 1,3h l h
33 t 1,1 l11 t 2,3h l h
13 t 2,1
(3.31)
where
n
t t ( i 1)
l
ij t hk , j t t k ( i 1)
t i u (3.32)
k 1
108
CCLineSpring and CCPlaneSpring elements were created to enable convenient definition of
„uniform“ spring-like conditions along the boundaries. The boundary force at a node i of the
spring element is calculated:
ui kA
Ri (3.35)
n direction
where
k is spring material stiffness parameter set by &MATERIAL SPRING command,
(parameter k has character of multi-linear Young modulus),
ui is displacement at spring element node i ,
A is the area of CCPlaneSpring element or length of CCLineSpring multiplied by
thickness (which defaults to 1 if not specified in element geometry) or the area defined in
element geometry for CCSpring (similarly, with a default of 1 if not specified) for the
respective element,
n is number element nodes, i.e. 1, 2 or 3 for CCSpring, CCLineSpring or CCPlaneSpring
element respectively,
direction is Euclidean norm (i.e. length) of the direction vector, see above.
CCSpring CCLineSpring
CCPlaneSpring
CCSpring
area A
110
Using the quadratic interpolation function, the displacement components u( i ), v( i ) is written
in the terms of triangular coordinates i and nodal displacement vectors :
u ( i ) F ( i )T u, v ( i ) F ( i )T v (3.37)
The displacement vectors u, v contain six components of the nodal displacements and the vector
F ( i ) contains the quadratic interpolation functions in triangular coordinates:
u u1 u2 u3 u4 u5 u6 , v v1 v2 v3 v4 v5 v6
T T
(3.38)
u
( i ) F T (3.42)
v
The stiffness matrix:
K F DF dV
T
(3.43)
V
The matrix F contains partial derivatives of the interpolation function F and the integral in the
last equation is made over the element volume V. The details of the derivation can be found in
FELIPPA 1966 and here only the final matrix equations are presented.
Fig. 3-13 Quadrilateral element (b) composed from two triangular elements (a).
The quadrilateral finite element is composed from two 4-node triangular elements, as shown in
Fig. 3-13. Two degrees of freedom in a node are the horizontal and vertical displacements. The
triangular element is derived from the 6-node triangle by imposing kinematic constraints on two
mid-side nodes. The resulting strain-displacement matrix relation for the 4-node triangle is:
e x x1 x 2 x 3 , e y y1 y 2 y 3 , g x1 x 2 x 3
T T T
(3.45)
u u1 u2 u3 u4 , v v1 v2 v3 v4
T T
(3.46)
The strain interpolation function in the element is linear and is uniquely specified by three nodal
values in the corners of the triangular element, while the displacement interpolation function is
quadratic and is specified by three corners and one mid-side nodal displacement. The
components ui, vi are the horizontal and vertical displacements, respectively, in the node i. The
indexes 1, 2 and 3 denote the corner nodes of a sub-triangle and the index 4 is for the mid-side
node, see Fig. 3-13 (a). The strain-displacement sub-matrices in (3.44) are
3b1 2b3 b2 b3 4b2
1
U b1 3b2 b3 b3 4b1
2S
b1 b2 b3 .
a1 x3 x2 b1 y2 y3
a2 x1 x3 b2 y3 y1
a3 x2 x1 b3 y1 y2
2S a3 b2 a2 b3
where xi, yi are the global Cartesian coordinates of the node i in a sub-triangle, S is the area of the
sub-triangle.
The element stiffness matrix for the 4-node sub-triangle is
K Kuv
K uu (3.48)
K vu K vv
The stiffness matrix K has an order 8 and is so partitioned that the upper four rows correspond to
the horizontal displacement components (index u) and the lower four rows correspond to the
vertical displacement components (index v). The integration of the stiffness coefficients is made
exactly and the resulting sub-matrices are:
K uu St d11 A d13 ( H H T ) d 33 C
K vv St d 22 C d 23 ( H HT ) d33 A
112
where t is the thickness of the element, dij are the coefficients of the material stiffness matrix D,
(3.40). The integration in (3.43) is done explicitly by the following matrix multiplication:
A UT QU, H UT QV, C VT QV (3.50)
Where the area integration matrix Q is:
2 1 1
1
Q 1 2 1 (3.51)
12
1 1 2
The element stiffness matrix of the 5-node quadrilateral, Fig. 3-13(b), is composed of the two 4-
node sub-triangles by summing the stiffness coefficients of the appropriate nodes. The resulting
matrix of the 5-node quadrilateral K10 has the order 10. The coefficients of the matrix can be
rearranged according to the external (index e) and internal (index i) degrees of freedom:
K K ei
K10 ee (3.52)
K ie K ii
The sub-matrices corresponding to two internal degrees of freedom are eliminated by
condensation procedure and the final element stiffness matrix K of the order 8 is obtained:
1
K K ee K ei K ii K ie (3.53)
114
3.8 External Cable
External pre-stressing cables are reinforcing bars, which are not connected with the most of the
concrete body, except of limited number of points, so called deviators, as shown in Fig. 3-15.
This element type is denoted in ATENA as CCExternalCable.
F2 F1 e i p Q f ( ) f r (r ) (3.55)
Fig. 3-17 Forces and displacements in the cable element (cable section).
A section of the cable between the deviators is considered as the uniaxial bar element, Fig. 3-17.
The force F in the cable element depends on: the pre-stressing force P, the displacements of ends
u1, u2 due to structural deformation and the cable slips 1, 2 in the deviators. The slip is
introduced as an additional variable for the external cables. The equilibrium equation of the
cable section is:
F P k (u2 u1 2 1 ) (3.57)
The element stiffness k = Es A/L, where A, L are the cable cross section and the length,
respectively, and Es is the actual secant or tangent modulus derived in the same way as in case of
other reinforcement using bilinear or multi-linear law.
The cable forces F1, F2,…, are determined by applying the above equations in all cable deviators
and by simultaneous solution for slips i .
Introduction of pre-stressing is accomplished by applying an initial slip (cable pull-out) at the
anchor end until a prescribed pre-stressing force is reached. This procedure reflects a real
process of pre-stressing and takes into account the loss of pre-stressing due to friction deviators
and deformation of the structure.
116
segment of the bar. The original length l0 will change to l due to displacement u of the
surrounding body and bar slips .
c c
1 2 i-1 i m
undeformed truss i
deformed truss i
u
i i i+1
u
i+1
lo
l
Discretized form of (3.60) for node i, (considering elements i 1, i ), reads (the bars are of
constant strain type):
Fi L Ai 1 i 1
Fi R Ai i
li li 1 (3.61)
for ( Fi R Fi L ) : Fi R Fi L p c
2
l l
for ( Fi R Fi L ) : Fi R Fi L i i 1 p c
2
If this element acts also as the external cable, see the previous section, then
for ( Fi R Fi L ) : Fi L Fi R dia dib
Fi R Fi L Fi R Fi R dia dib
Fi R Fi L Fi R 1 dia dib )
(3.62)
for ( Fi Fi ) :
R L
Fi Fi d d
R L
i
a
i
b
Fi R Fi L Fi L Fi L dia dib
Fi R Fi L Fi L 1 dia dib
Assembling (3.61) and (3.62) yields final (in)equations for force difference at node i:
li li 1
for ( Fi R Fi L ) :
Fi R Fi L Fi R 1 dia dib
2
p c
(3.63)
l l
for ( Fi R Fi L ) :
Fi R Fi L Fi L 1 dia dib i i 1 p c
2
The above set of (in)equations is calculated in iterative manner. Assume we know the forces at
iteration (k 1) , then the forces at iteration k are:
118
EAi
Fi R , k 1 (u i1 u i ik11 ik 1 )
li
EAi 1
Fi L , k 1 (u i u i1 ik 1 ik11 )
li 1
EAi
Fi R , k Fi R , k 1 ( ik1 ik )
li
EAi 1
Fi L , k Fi L , k 1 ( ik ik1 ) (3.64)
li 1
k k 1 k
i 1 i 1 i 1
k
i i
k 1
i k
k k 1 k
i 1 i 1 i 1
and
for ( Fi R Fi L ) :
EAi EA l l
Fi R
li li 1
2
( ik1 ik ) Fi L i 1 ( ik ik1 ) Fi R 1 dia dib i i 1 p c
for ( Fi R Fi L ) :
EAi EA l l
Fi R
li li 1
2
( ik1 ik ) Fi L i 1 ( ik ik1 ) Fi L 1 dia dib i i 1 p c
EAi 1 EAi EAi 1 k EAi l l
i1
k
i
ik1 Fi R Fi L Fi L 1 dia dib i i 1 p c
li 1 li li 1 li 2
(3.65)
If the above equation is written for all nodes on the bar, we obtain a set of inequalities. It has to
be solved in iterative manner (within each iteration of the main solution loop).
Atena also support so called CCBarWithMemoryBond 2D and 3D elements. They differ from
their original formulation, (i.e. elements CCBarWithBond), in that they have different function
f ( ) for "loading" and f , unload ( ) for "unloading" regime. This means ( min , max ) in the
former and ( min , max ) in the latter case.
In order to obtain more realistic shape, the resulting cohesion stresses are prior their output
smoothed. The smoothing operation for node i is expressed as follows:
i 1li 1 i li
right
li 1 li
i li i 1li 1
left (3.66)
li li 1
120
CCIsoGap<xxxxxx> for 2D and CCIsoGap<xxxxxxxxxxxx> or
CCIsoGap<xxxxxxxxxxxxxxxx> for 3D. The string in < > describes present element nodes, (see
Atena Input File Format document for more information). The elements are derived from the
corresponding isoparametric elements (described in sections 3.3 and 3.4), i.e. they use the same
geometry and nodal ids etc. Geometry of the supported gap elements is depicted in Fig. 3-21.
v 1
2D v 1
r,u(r) u(r)
r,
2 2 5
4 4
6
3 3
CCIsoGap<xxxx> CCIsoGap<xxxxx_x>
w 3D 10 w
1 1
2 7
2 9
r,u(r,s) r,u(r,s) 8
4 3 4
3
5 5 12
11
6 6
s,v(r,s) s,v(r,s)
CCIsoGap<xxxxxx> CCIsoGap<xxxxxxxxxxxx>
13 w
w 1
1 9
2 2
12 16
5 5
4 r,u(r,s) 6 10 4
r,u(r,s) 6 11
3 3
14
8 8
7 15
7
s,v(r,s) s,v(r,s)
CCIsoGap<xxxxxxxx>
CCIsoGap<xxxxxxxxxxxxxxxx>
The interface is defined by a pair of lines, (or surfaces in 3D) each located on the opposite side
of interface. In the original (i.e. undeformed) geometry, the interface lines/surfaces can share the
same position, or they can be separated by a small distance. In this case we speak about the
interface with nonzero thickness.
In the following, the interface behavior is explained on a simple 2-dimensional case, see section
2.6 for a full description of the interface material.
The interface element has two states:
Open state: There is no interaction of the contact sides.
Closed state: There is full interaction of the contact sides. In addition, friction sliding of
the interface is possible in case of interface element with a friction model.
in which u , v are the relative displacements of the interface sides (sliding and opening
displacements of the interface) in the local coordinate system r , s and K tt , K nn are the shear and
normal stiffness, respectively. This coefficient can be regarded as stiffness of one material layer
(real, or fictious) having a finite thickness. It should be understood that the layer is only a
numerical tool to handle the gap opening and closing. F , F are forces at the interface, (again
at the local coordinate system).
The actual derivation of gap elements is now demonstrated for the case of linear 2D gap element
CCIsoGap<xxxx>, see Fig. 3-21. The other elements are constructed in a similar way.
The element has two degrees of freedom defined in the local coordinate system, which is aligned
with the gap direction. They are relative displacements v, u and are defined as follows:
1 1
h1 (1 r ), h2 (1 r )
2 2
u h1u1,4 h2 u2,3
u
v h1v1,4 h2 v2,3
u1
v
1
u2
h 0 h2 0 h2 0 h1 0 v2
u 1 Bu
0 h1 0 h2 0 h2 0 h1 u3
v3 (3.68)
u
4
v4
The rest of the element derivation is the same as in case of any other elements, i.e. the stiffness
matrix K ΒT DBdV , vector of internal forces Q ΒT FdV etc. A numerical integration in
two Gauss points is used to integrate the interface element stiffness matrix. The matrix K a nd
the vector Q are in local coordinate system and therefore before they are assembled in the
problem governing equations they must be transformer in global coordinates.
The stiffness coefficients depend on the gap state. The interface is considered open, if the normal
force F >Rti (Rti is the interface tensile strength force) and the corresponding constitutive law is
(stress free interface):
F 0
(3.69)
F 0
122
The stiffness coefficients are set to small, but nonzero values K ttop , K nnop .
The interface element is considered closed if F Rti. The stiffness coefficients are set to large
values K ttcl , K nncl . It should be noted that the stiffness coefficients are defined only for the purpose
of the numerical iterative solution. (Hint: The values of coefficients in the closed state (the large
values) are based on thickness comparable to the size of neighbor quadrilateral elements. The
minimum values in the open state can be about 1000 times smaller. )
The interface thickness in the out-of plane direction is normally provided as an input parameter.
In the case of axi-symmetric analysis it is however calculated using the formula:
t 2 x (3.70)
where x is the distance from the axis of symmetry.
There are two special options for processing the gap elements:
Initial gap opening
It is possible to "open" gap at a particular load step, typically the first step of the analysis, i.e. we
can introduce to the gaps something like initial element strains in case of ordinary finite
elements. This is achieved by LOAD INITIAL GAP ... INIT_STEP_ID step_id command. Upon
that, during calculation of the (gap) element at the step step_id an artificial opening of the
interface is introduced. Its value is the distance between upper and lower element surfaces/lines
(with reference to undeformed structural shape).
The GAP element load is typically used as follows: we have a structure with a base and upper
blocks. The upper block falls down towards the base block that is typically fixed. The structure is
solved by introducing a layer of gap elements between the base and upper blocks and applying
the GAP element load (for these gaps elements) in the 1st step. As a results, in the first steps the
gaps will open to the distance between the blocks. It involves some tensional forces, but as the
interface material usually sustains only compression forces, they can be neglected. In next steps
the upper block gradually is falling down to the base block until it hits it. At this moment
interface gaps get fully closed, they change their regime form tension to compression and the
upper block gets fully supported by the base block.
Moving gaps 3
Suppose we have a structure has a base block and an upper block sitting on the base block. The
base block is fixed, the upper block is dragged on the upper surface of the base block. The blocks
are not mutually interconnected, only some friction and cohesion forces exist between them.
Such problems can be modelled by the RESET_DISPLS n flag for the CC2DInterface /
CC3DInterface. If this flag is input, then the upper and bottom surface/lines for all corresponding
elements are realigned at the end of each step as shown for 2D elements in the following picture.
The 3D gaps element are realigned in the same way.
Of course, the boundary surface/lines projection of the gap interface (and thus its "moving" can
be used in more complex situation, but the essence of the described technique remains the same.
The layer of interface elements are typically connected to the bottom/ upper block of structure by
MASTER SLAVE NODAL LISTS boundary conditions, where we must not forget to use the
flag PROCESS_FLAG USE_CURRENT_COORDS. It will assure that after realigning the
interface gets properly connected to the rest of the (deformed) structure.
3
Available starting from ATENA version 4.3.1.
y
y
2
1
1
CCCircumferentialTruss CCCircumferentialTruss2
x x
124
In the following structural vectors and matrices for the CCIsoTruss element are derived.
Development of the CCIsoTruss2 is much the same. In fact, it is CCIsoTruss acting at the centre-
point of the CCIsoTruss2 element with its cross-sectional area calculated as explained above.
The element vectors and matrices for Total Lagrangian formulation (TL), configuration at time t
and iteration (i) are as follows. Note that they are equally applicable for Updated Lagrangian
formulation (UL) upon applying changes related to the element reference co-ordinate system
(undeformed vs. deformed element axis.).
The truss element center has at reference time t and t t ( i 1) co-ordinates t X [t x1 , t x1 ] and
t t ( i 1) ( i 1) ( i 1)
X [t t x1 , t t x1 ] , respectively. The element length (at respective time) is its length is
t
l 2 t x1 and t t ( i 1)
l 2 t ( t x11 t u11( i 1) ) .
t t ( i ) t t
Increment of Green Lagrange strain t 11( i )
t 11 t tt11( i 1) (at time t t , iteration (i )
with
to configuration at time t ) is calculated:
1 l l
t t ( i ) 2 t t ( i 1) 2
(i )
(3.71)
tl
t 11
2 2
where truss length t t ( i )
l 2 t ( t x11 t u11( i 1) t u11( i ) ) . Note that t u11( i ) is co-ordinate increment
(i) ( i 1)
( t t x1 t t x1 ) . Substituting expressions for element length into (3.71) yields:
(i )
4 2 x u
t 1
1
t 1( i 1)
1 t u11( i ) t x11 t u11( i 1)
2 2
x
t 11
t 1 2
1
(3.72)
1( i 1) 1( i ) 2
u1( i )
u u t
1 u 1( i )
t 1
1 t 1
t 1
x1 t 1 2
t 1 t 1
x 1 2 x 1
t t 1
t BL 0 t
(3.73)
x11
t t ( i 1)
t
u11( i 1)
B (3.74)
x
t L1
t 1 2
1
and
t t ( n 1) 1
BNL (3.75)
x
t
t 1 2
1
x11 t u11( i )
t 2
t t
t X (i )
1,1 1 e(i )
t 11 t
(3.77)
x11
2 t 1 2
4 x1 u1 4 x1 x11 t u11( i ) t x11
2 t 1 t 1( i ) 2 t 2
t t ( i )
l
e(i )
t 11 (3.78)
4 2 t x11
t t 1
l 2 x1
126
The essential point in the element’s derivation is that displacements and rotations fields are
approximated "independently", (see e.g. (Jendele 1981), where similar approach is used for
plates). This means that they are handled separately. Unlike in true Mindlin theory our
formulation matches geometric equations automatically. However, a special technique is used to
improve the element’s shear behaviour (Hinton and Owen 1984).
The first formulation of this element proposed by Ahmad was linear but since that time many
improvements have been achieved. The most important is the application of reduced or selective
integration scheme that reduces or totally removes locking of the element. Also, many authors
extended the original formulation to geometrically and later also materially nonlinear analysis.
One such an advanced form of the element is the formulation implemented in ATENA.
On input, the Ahmad element uses the same geometry as 20 nodes isoparametric brick element,
i.e. CCIsoBrick<xxxxxxxxxxxxxxxxxxxx>, see Fig. 3-27. This is needed, in order to be able to
use the same pre- and postprocessors’ support for the shell and native 3D brick (i.e. hexahedron)
elements. After the 1st step of the analysis, the input geometry will automatically change to the
external geometry from Fig. 3-27. As nodes 17 and 18 contain only so called bubble function,
the element is post-processed in the same way is it would be the element
CCIsoBrick<xxxxxxxxxxxxxxxx>. Internally, all element’s vectors and matrices are derived
based on the internal geometry as depicted also in Fig. 3-27.
With shell elements, the best connection at edges is to cut both at 45 degrees, or a different
corresponding angle if the thicknesses are not the same, or if connected at other than right angle,
see Fig. 3-24 (a). Another option is to use a volume brick element at the corner, which is the only
feasible way when more than two shells are connected, see Fig. 3-24 (b). The nodes on the
surface connected to the volume element have to be listed in the INTERFACE subcommand in
the shell geometry definition for correct behavior. Connecting like in Fig. 3-25 is not
recommended, as the master-slave relations induced by the fixed thickness of the shell may
cause numerical problems.
Shell1
Brick Shell2
Shell3
(a) (b)
Fig. 3-24: Shell - recommended connection (a) 2 shells (b) 3 shells
128
Fig. 3-26 Ahmad element coordinate systems
t t
The vectors V1k , V2k , V3k are defined as follows: Firstly, two auxiliary vectors V1 , V3 are
t t t
t
calculated. Vector V3 at a point is defined as a line joining bottom and top coordinates at the
t
node k (prior any deformations, i.e. at reference configuration). The second vector V1 is normal
t
to V3 and is parallel to plane of global X 1G and
0 0
X 3G . Hence:
(3.79)
V1 V1 1 , V1 2 , V1 3 V3 3 , 0, V3 1
t t t t t t
t t
If V3 is parallel to X 2G (i.e. V3 1 V3 3 0 ), V1 is defined by
0 t t
V1 V3 2 , 0, 0
t t
(3.80)
t t t t
After that, the coordinate system V1k , V2k , V3k itself is defined. The vector V3k is constructed in
t
the same way as was the vector V3 , however, current, i.e. deformed configuration is used. The
remaining two vectors are defined as vector product:
t t
V2k V3k V1
t
(3.81)
t t t
V1k V2k t V3k (3.82)
t x1 t x1
r s
t t
t x x
X 3L 2 2 (3.83)
r s
t t
x3 x3
r s
t t t
X 2L X 3L V1k
t t t
X 1L X 2L X 3L (3.84)
t t t t
As the nodal coordinate system V1k , V2k , V3k can rotate along V3k , the local coordinate system
t t t t
would X 1L , X 2L , X 3L rotate simultaneously along X 3L . This definition allows for user defined
shell local coordinate system that is common for all shell elements, irrespective of their
t t
incidences. Note that unlike V3k the vector X 3L is always normal to the element mid-plane
surface.
Curvilinear coordinate system
This system is used to calculate derivatives and integration in element integration points. Its
coordinates are r , s for in-plane direction and t in direction of element thickness, see Fig. 3-26.
The in-plane displacements are approximated by Lagrange, Hetherosis or Serendipity
approximation similar 2D isoparametric elements. For the 3rd direction, i.e. through the depth of
the element. linear approximation is used within the frame of the shell layer concept.
130
t3 Input geometry
11 10
4 2
12 9
1
20 19
7 18
15 s
8 17 14
6
r 16
13
5
t3
11 10 External geometry
4 2
12 9
117
7
15 s
8 14
18 6
r 16
13
5
t
Internal geometry
1
2 8
3 9 7
s
4 6
5
r
132
3.12.2 Geometry Approximation
The coordinates of the top and bottom element surface is used to define the element geometry:
t x1 t x1k ,top t x1k ,bot
t t N
1 t t k ,top 1 t t k ,bot
x2 x hk 2 x2 2 x2 (3.85)
t x3 k 1 t x3k ,top t x3k ,bot
where N=8 is number of nodes per element, (geometry is always interpolated by 8-nodes
Serendipity interpolation, irrespective of displacement interpolation), h(r,s) is k-th interpolation
t x1k ,top t x1k ,bot
function, r,s,t are isoparametric coordinates (see Fig. 3-27), t x2k ,top and t x2k ,bot are vector of
t x3k ,top t x3k ,bot
top and bottom coordinates of point k, see Fig. 3-29.
t
X k ,top [ t x1k ,top , t x2k ,top , t x3k ,top ]
node k t
X k ,mid [t x1k , mid , t x2k ,mid , t x3k ,mid ]
t
X k ,bot [t x1k ,bot , t x2k ,bot , t x3k ,bot ]
rotations with respect to vectors v1k and v 2k respectively. These degree of freedoms (DOFs) are
used throughout the whole element’s development. However, in order to improve compatibility
of the present shell element with other 3D elements implemented in ATENA, externally the
T
element uses t u1top ,t u2top ,t u3top ,t u1bot ,t u2bot DOFs, i.e. displacements at the top and bottom of the
element. The 6th displacement, i.e. u3bot is eliminated due to application of shell theory that
assumes 33 0 .
Approximation of the original three "displacement" and two rotation degrees of freedom is
independent. Nevertheless, the curvatures used in governing element equations use all of them in
the sense dictated by geometric equations. This approach enables to satisfy not only equilibrium
equations for membrane stresses and in-plane shear (in mid-surface) as it is the case of popular
Kirchhoff hypothesis, but also to satisfy equilibrium condition for transversal shears (normal to
mid-surface).
Note that in the following derivation of the element we will deal with the original set of
element’s DOFs , see (10). Every point thus has five degree of freedom,
T
t u1mid ,t u2mid ,t u3mid ,t ,t . Displacement vector is calculated by:
t k , mid t V2k t
V1k 1
t u1k u1 1 t k
t k t t k , mid t t k
N
V k
u2 u hk u2 2 thick k
t
V1 2 t k (3.88)
22
t u3k k 1
t u3k , mid t k t k
V2 3 V1 3
T
The original displacement vector at point k has the form t u1mid ,t u2mid ,t u3mid ,t ,t . Unlike in the
case of geometry approximation, were N=8, displacements approximation accounts also for
displacement in the element mid-point, i.e. N=9. The ninth function h is so called bubble
function.
u k ,mid [ t u1k ,mid , t u2k ,mid , t u3k ,mid ]
t
V3k
node k
k
t t
V2k
t
k
t
V1k
134
3.12.4 Strain and Stresses Definition.
The 2nd Piolla Kirchhoff tensor and Green Lagrange strain tensor is used. They are calculated
and printed in the local coordinate system t x '1 , t x ' 2 and t x '3 .
Green - Lagrange tensor.
The general definition for Green-Lagrange strain tensor has the form (see eq. (1.8)):
1 t
t
0 ij
2
0 ui, j 0t u j ,i 0t uk ,i 0t uk , j (3.89)
Using the above equation and applying the Von-Karman assumption, Eqn. (3.89) can be written
as:
t u1
2
x1
t
1 u3
t
t u2 2 t x1
0 11
t
2
x2
t
t 1 t u3
0 22 t u t u t
2 0t12 t 1 t 2 2 x2 0t L 0t NL (3.90)
t x2 x1 t u t u
2 0 13 t t 3 t 3
u t
u
2 t 1 3 x2 x1
0 23 t x t x
3 1 0
t
u2 u3
t
0
t x t x
3 2
information about their calculation is beyond the scope of this publication. It is available e.g. in
(Jendele 1992).
2nd Piolla Kirchhoff tensor.
Energetically conjugated with Green - Lagrange tensor is 2nd Piolla Kirchhoff tensor and this
tensor is used by the present shell element. Remind that we account for all stresses with
exclusion of normal stress which is perpendicular to shell mid surface (as it is usual practice in
shell analysis). This is the reason, why we introduced local coordinate system and all expression
are derived with respect to it.
Obviously the local coordinate system varies depending on element deformation and thus it is
necessary to re-compute (each iteration) the transformation matrix T (that relates local and
global coordinate systems).
Note that that it is possible to abbreviate full 3 by 3 element tensor to the above vector, because
of adopting Von Karmann simplifying assumption.
where hi are values of interpolation function at point (0,0), ai are corresponding node values,
a9 is departure in the center (i.e. computed value corresponding to degree of freedom at center)
and a9 is total value in center.
Depending on combination how many nodes and integration points are used, we distinguish the
Serendipity, Lagrange and Heterosis degenerated element variants, see Fig. 3-31.
136
Serendipity element.
This element was used in the original Ahmad work. It comprises eight nodal points (center point
corresponding to bubble function is omitted).
Gauss integration scheme is used for integration . It can be integrated by full, reduced or
selective integration procedure. Using full integration ,i.e. at three by three sample points,
element exhibits shear locking for thin and even moderately thick element. If reduced integration
is used. the problem of locking is significantly improved without creating spurious energy modes
on structure level. However, in case of thin element there are two non communicable spurious
energy mode on element level.
It should be noted that there were reported some difficulties, if some unfavorable constraints are
applied. Nevertheless the element is popular. If reduced integration is used the provided results
are relatively good.
Fig. 3-31 Node notation for element variants of the Ahmad shell element
Nine node Lagrangian element.
The nine point Lagrangian element is still considered to be the best variant of the degenerated
element. This is especially because of its versatility. For full integration scheme there is no
problem with membrane and shear locking for very thin plate and shell application. If element is
moderately thick, shear locking can be improved by reduced integration scheme. However, in
that case the element exhibits rank deficiency.
138
Fig. 3-33 Integration schemes and sampling point notation
The steps during selective integration of shear can been explained on example of integration
arbitrary function f ( r , s ) :
1/ First we calculate the value of f at sampling points that corresponds to two by two integration
rule, i.e
f (-0.5773, -0.5773), f (-0.5773, 0.5773), f (0.5773, -0.5773), f (0.5773, 0.5773)
2/ Using bilinear approximation we calculate the values of f at points that correspond to three by
three integration rule. There are two possibility to that.
The first one is based on original approximation area and the main idea is that we calculate the
value of function f at "corners" of isoparametric element (i.e. r 1., s 1. ):
k 1
4
f (0.5773, 0.5773) f i hi' (0.5773, 0.5773)
k 1
where fi are element nodal values of function f and hi' are interpolation functions corresponding
to two by two interpolation and a node i.
140
Fig. 3-35 The layer model
tV k t
V2k 1
t
V3k 1 0 0 0
11 t k , N ,top
u1 t k
t k ,top
t t u1 t u1k , N ,top
t k ,top V1 2 V2k 2 V3k 2 0 0 0 t k , N ,top t k , N ,top
u
u2 t k t t 2k , N ,top u2
t u3k ,top V1 3 V2k 3 V3k 3 0 0 0 u t k , N ,top
u3
t k ,bot T2 t k , N ,bot
3
t k t k , N ,bot
(3.96)
u1 0 V3 1 u1 u1
t t
0 0 V1k 1 V2k 1
t u k ,bot t k , N ,bot
t k
t k , N , bot
u2 u
t 2k ,bot 0 V3 2 t k , N ,bot t 2k , N ,bot
t t
0 0 V1k 2 V2k 2
u3
t t k
u3 u3
0 0 0 V1k 3 V2k 3 V3 3
Complete transformation of the original DOFs to the new element formulation DOFs
The final transformation from the original to the new element DOFs at a node k is obtain by
substituting (3.95) into (3.96). Thus we can write
t u1k ,top
t k ,top t u1k ,mid t u1k ,mid
u2 t k ,mid t k ,mid
t u3k ,top u2 u2
t k ,bot T2 T1 u3 T t u3k ,mid
t k , mid
(3.97)
u1 t k t k
u
t k , bot
t k ,bot
2
t k tk
u3
where T
142
thick k t k thick k t k
t k 2 t 2 t 2 t t t t t t t t t t t t
V1 1 V2k 1 V3k 1 V1k 1 V1k 2 V2k 1 V2k 2 V3k 1 V3k 3 V1k 1 V1k 3 V2k 1 V2k 3 V3k 1 V3k 3 V1 1 V2 1
2 2
t k t k
V V V thick t k
k k
t k t k t k t k t 2 t 2 t 2 t t t t t t thick t k
V1 1 V1 2 V2 1 V2 2 V3 1 V3 3 V1k 2 V1k 3 V2k 2 V2k 3 V3k 2 V3k 3
k k k
1 2 2 2 3 2 V1 2 V2 2
2 2
t k t k
V V V thick t k
k k
t k t k t k t k t t t t t t t 2 t 2 t 2
thick t k
V1 1 V1 3 V2 1 V2 3 V3 1 V3 3 V1k 2 V1k 3 V2k 2 V2k 3 V3k 2 V3k 3
k k k
1 3 2 3 3 3 V1 3 V2 3
2 2
thick k t k thick k t k
t k 2 t 2 t 2 t t t t t t t t t t t t
V1 1 V2k 1 V3k 1 V1k 1 V1k 2 V2k 1 V2k 2 V3k 1 V3k 3 V1k 1 V1k 3 V2k 1 V2k 3 V3k 1 V3k 3 V1 1 V2 1
2 2
V V V
k
t k t k t k t k t k t k t 2 t 2 t 2 t t t t t t thick t k thick k t k
V1 1 V1 2 V2 1 V2 2 V3 1 V3 3 V1k 2 V1k 3 V2k 2 V2k 3 V3k 2 V3k 3
k k k
1 2 2 2 3 2 V1 2 V2 2
2 2
t k t k
V V V thick k t k thick k t k
t t t t t t t t t t t 2 t 2 t 2
V1 1 V1 3 V2k 1 V2k 3 V3k 1 V3k 3 V1k 2 V1k 3 V2k 2 V2k 3 V3k 2 V3k 3 k k k
V1 3 V2 3
1 3 2 3 3 3
2 2
In a very similar way we can define inverse transformation, i.e. from the new DOFs to original
one. Without any derivation the matrix reads:
t u1k ,top
t k , mid
1 u t k ,top
t k , mid
u
u2
2 t u k ,top
t k , mid
3 u T' t 3k ,bot (3.98)
u1
t k
t u k ,bot
t
k t 2k ,bot
u3
u3k ,bot c1top t u1k ,top c2top t u2k ,top c3top t u3k ,top c1bot t u1k ,bot c2bot t u2k ,bot (3.100)
where:
t t t t
V1k 3 V1k 1 + V2k 3 V2k 1
c top
V V 1
1 2 2
t k t k
1 3 2 3
t t t t
V1k 3 V1k 2 + V2k 3 V2k 2
c2top
V V 1
t 2 t 2
k k
1 3 2 3
1 V V yy
t 2 t 2
k k 2
1 3 2 3 z
c top
V V 1
3 2 2
t k t k
1 3 2 3
t t t t
V1k 3 V1k 1 + V2k 2 V2k 1
c bot
V V 1
1 2 2
t k t k
1 3 2 3
t t t t
V1k 3 V1k 2 + V2k 3 V2k 2
c bot
V V 1
2 2 2
t k t k
1 3 2 3
(3.101)
The DOFs u1k ,bot or u2k ,bot can be eliminated in the same way. During the execution of the
t
element, it is recommended to constrain one of u1k ,bot ,t u2k ,bot ,t u3k ,bot based on which solution is
the most stable, (i.e. maximum denominator in (3.101)).
Constraining DOFs at the centre of Hetherosis element
A special attention needs to be paid to the 9th mid-plane node of Hetherosis element, when we
have to additionally constrain t u1k , mid , t u2k ,mid , t u3k , mid . Thus, of the 6 DOFs we need to constrain 4
of them.
t
For example, suppose we want to keep free u2k ,top and t
u3k ,top and we need to
fix t u1k ,top ,t u1k ,bot ,t u2k ,bot ,t u3k ,bot . Equation (3.99) from the previous paragraph needs to be added by
three more equations. These are:
144
t u1k ,top
t k ,top
u
t u1k ,mid T11' T12' T13' T14' T15' T16 t 2k ,top 0
'
t k ,mid ' u
u2 T21 T22 T23 T24 T25
' ' ' '
T26' t 3k ,bot 0 (3.102)
u
t u3k ,mid T31' T32' T33' T34' T35'
T36' t 1k ,bot 0
u
t 2k ,bot
u3
t
Equations (3.99) and (3.102) are then solved for u1k ,top ,t u1k ,bot ,t u2k ,bot ,t u3k ,bot as a linear
combination of t u2k ,top and t u3k ,top .
1
t u1k ,top T11' T14' T15' T16' T12' t u2k ,top T13' t u3k ,top
k ,bot ' ' t k ,top ' t k ,top
u1 T21 T24' T25' T26' T22 u2 T23 u3 (3.103)
u2k ,bot T31' T34' T35' T35' T32' t u2k ,top T33' t u3k ,top
t k ,bot top top t k ,top top t k ,top
u3 c1 1 c1bot c2bot c2 u2 c3 u3
Again, there are several alternatives regarding of which of the 6 DOFs to keep and which to
eliminate. The best option is chosen the same way as described in Section 0.
146
3.13.1 Geometry and displacements
2
t t ( i ) t
xi hk t t X ik ,(i ) ak t tVi nk ,(i )
2
where i=1,2,3 is index relating to global axes x1 , x2 , x3 , (i.e. x,y,z), k 1...nG , nG = number of
the element's nodes used to approximate geometry, typically 8 or 9. Note that due to Shell
theory the shell thickness at node k 0 ak t ak t t ak(i 1) t t ak( i ) ak . The symbol t t Vi nk ,( i ) is
ith coordinate, ( i 1, 2, 3 for coordinate x, y , z ), of the vector V n at node k at time t t ,
iteration (i ) . The vector V n is normal to the shell. Later we will use also vectors V 1 ,V 2 ,
(V 1 V 2 V n ) . They will constitute base vectors for shell's bending rotations , .
t t
t t
ui(i 1) hk t t X ik ,( i 1) ak Vi nk ,(i 1) t X ik ak tVi nk
t t
2 2
(3.106)
t
2
hk t t X ik ,(i 1) t X ik ak t t Vi nk ,(i 1) t Vi nk
Note that in this case k 1...n , n is number of nodes to approximate displacements. Current
implementation of the shell elements assumes ng n , (which differs for Ahmads elements).
t
ui hk
t t
X ik ,(i ) t t X ik ,(i 1) ak
2
t t
Vi nk ,(i ) t t Vi nk ,(i 1)
(3.107)
t
hk U ik ak
2
Vi Vi nk ,(i 1)
t t nk ,( i ) t t
At each node the element has 5 DOFs: 3 displacements U ik and two rotations k , k described
below:
T
For the next derivation let us assume a general vector vL v1L , v2 L , v3 L with unit length that is
subject to rotations
[ L , L , L ]T , (where the subscript L indicates that both the vector and the rotations are defined
with respect to the local coordinate system (defined by t t Vi1k ,(i 1) , t tVi 2k ,(i 1) , Vi nk ,(i 1) ). The
t t
rotations of the vector will produce displacements, (all in the local CS)
u1L 0 v3 L v2 L L v3 L v2 L L
u v 0
v1L L v3 L L v1L L (3.109)
2 L 3L
u3 L v2 L v1L 0 L v2 L L v1L L
Now assume the same behaviour for a vector normal to the shell's surface (again in the local CS
and unit lenght), i.e. vL t t ViLnk ,(i 1) 0, 0,1 . When this vector gets rotated, it produces
displacements, (see (3.110):
u1 V1 k
t t 1 ,( i 1)
L t t V12k ,(i 1) L
u t t V 1k ,(i 1) t t V 2k ,(i 1)
2 2 L 2 L (3.111)
u3 V3
t t 1k ,( i 1)
L V3
t t 2 k ,( i 1)
L
148
Substituting now L t t k ,(i 1) , L t t k ,(i 1) and uk t t Vi nk ,( i 1) , k 1..3 we can write
final equations for displacements due to rotations, (for iteration (i-1) and (i) and the difference):
t t
Vi nk ,( i 1) t t
k ,(i 1) t t
Vi 2k ,(i 1) t t k ,(i 1) t t
Vi1k ,(i 1)
t t
Vi nk ,( i ) t t
k ,( i ) t t
Vi 2k ,( i 1) t t k ,( i ) t t
Vi1k ,(i 1)
Vi nk ,( i ) t t Vi nk ,( i 1) t t k ,( i ) t t k ,( i 1) Vi 2k ,( i 1) t t k ,(i ) t t k ,(i 1)
t t t t
t t
Vi1k ,(i 1)
k t t
Vi 2k ,(i 1) k t t
Vi1k ,(i 1)
(3.112)
Hence, they represent rotation along two user defined vectors t t
Vi nk ,(i 1) . It is important to note
that the vectors t t
Vi nk ,(i 1) moves as the structure deforms.
t
ui hk U ik ak k
2
t t
Vi 2k ,(i 1) k t t
Vi1k ,(i 1)
(3.113)
If we need to connect a node k of the 2D shell element to an ambient solid element, for each
shell node we proceed as follows:
The shell's displacement at the bottom uibot , at the top uitop and their difference is:
a
uibot hk U ik k k
2
t t
Vi 2k ,(i 1) k t t
Vi1k ,( i 1)
a
uitop hk U ik k k
2
t t
Vi 2k ,(i 1) k t t
Vi1k ,(i 1)
(3.114)
The ambient solid's displacements at the same location is, (i.e. at the top and bottom of the shell
at node k)
where hhlbot and hhltop are the solid's interpolation functions at location top and bottom of the
shell's at node i, UU ibot ,l ,UU itop ,l are corresponding nodal displacements of the solid element.
The displacements' dofs U ik are connected as usually, i.e. using master-slave boundary
conditions. The rotational dofs k , k are computed by
Eqn. (3.117) sets 3 equations, i.e. uitop bot uuitop bot , i 1..3 for only 2 variables k and k . This
contradiction is resolved by transforming (3.117) into the shell's local coordinate system and
neglecting displacements perpendicular to the shell's surface, (i.e. in direction of thickness of the
shell). This conforms with the shell theory that states constant thickness of shells. The final
bot bot
equations read uitop
L
uuitop
L
, iL 1..2 , where the index L denotes local coordinate system of
the shell. Eqn. (3.117) then changes to
(3.118)
hhltopUU1top ,l hhlbotUU1bot ,l
ak k
V12k ,(i 1) k
t t t t
V11k ,( i 1)
TG 2 L hhltopUU 2top ,l hhlbotUU 2bot ,l TG 2 L ak k
t t 2k ,( i 1)
V2
k t tV21k ,( i 1)
hhltopUU 3top ,l hhlbotUU 3bot ,l
a k
k t t 2k ,( i 1)
V3
k t tV31k ,( i 1)
If we need to connect a node k of the 2D shell element to a node l of an ambient beam element,
i.e. to connect the shell's three global displacements and two local rotations to three global
displacements and rotations used by the beam CCIsoBeamBar, for each shell node we proceed as
follows:
150
Throughout this section the notation from the previous section is used. The shell's displacements
U ik are connected to beam's displacements UU il as usually, i.e. using master-slave boundary
conditions. As for rotations, the shell uses local rotations k , k along the vectors
t t
Vi 1k ,( i 1)
, t t
Vi 2k ,( i 1)
, and they must be transformed into three rotations l , l , l along
global axes x1 , x 2 , x 3 . Hence, for each slave node k of the shell we can write:
UU1l
U 11
k
0 0 0 0 0
UU 2l
U 0 1 0 0 0 0
k
2 UU l
U 0
k
0 1 0 0 0 l3 (3.119)
t t 1k ,( i 1)
3
0
k
0 0 t t
V11k ,(i 1) t t
V21k ,( i 1) V3 l
k t t 2k ,( i 1)
0
t t 2k ,( i 1) t t 2k ,( i 1)
0 0 V1 V2 V3 l
Note that it is assumed that the node k of the shell is located at position of the beam node l.
The elements are derived using Green-Lagrange strains and 2nd Piola Kirchhoff stresses. Green-
Lagrange strains at (ith iteration), i,j-axis x,y,z are calculated as follows :
t t ( i )
t ij
2
1 t t ( i ) ( i ) t t ( i ) (i ) t t (i ) t t (i )
t ui , j t u j ,i t u m ,i t um , j
1
( t tt ui(,i j1) t ui , j ) ( t tt u (ji,i1) t u j ,i ) ( t tt um(i,i1) t um ,i )( t tt um(i,j1) t um , j ) (3.120)
2
t t ( i 1)
t eij tij
t ij
where:
1 t t (i 1) t t ( i 1) t t (i 1)
t t ( i 1)
t ij
2
t ui, j t u j ,i t uk ,i t t ( i 1)
u
t k, j
1
e
t ij
2
t ui , j t u j ,i
t t ( i 1)
t u m ,i t u m , j t um , j t um ,i
t t ( i 1)
t eij0 t eij1
1 (3.121)
t eij
0
2
t ui, j t u j ,i
1 t t (i 1)
t eij
1
2
t u m ,i t u m , j t u m , j t um ,i
t t ( i 1)
1
t ij
2
t u m ,i t u m , j
f x y z f f
r r r
r x x
f x y z f f
J (3.122)
s s s
s y y
f x y z f f
t t t t z z
f f
x r
f J 1 f (3.123)
y s
f f
z t
Derivatives of coordinates at t with respect to r,s,t to calculate J:
t xi hk t k t t t k
X i ai Vni
r r 2
t xi hi t k t t t k
X i ai Vni (3.124)
s s 2
t xi 1
hi t ai tVnki
t 2
152
t t ui(i 1) hk t t k ,(i 1) t
Ui ak t t
dVi nk ,(i 1)
r r 2
t t ( i 1)
ui h t
k t t U ik ,(i 1) ak t t
dVi nk ,(i 1)
s s 2
t t ui(i 1) 1 (3.126)
hk ak t t
dVi nk ,(i 1)
t 2
t t
U ik ,(i 1) t t X ik ,( i 1) t X ik
t t
dVi nk ,( i 1) t tVi nk ,(i 1) t Vi nk
t t ui(i 1) t inv ,k t t ui(i 1) t inv ,k t t ui(i 1) t inv , k t t ui(i 1)
J j1 J j2 J j3 (3.128)
t xj r s t
ui , k hk t
xj
t
t J inv
j1
r
U i ak
k
2
k
Vi 2k ,(i 1) k
t t t t
Vi1k ,(i 1)
hk k t
t J inv
j2
,k
s
U i ak
2
k t t 2k ,( i 1)
Vi k t tVi1k ,( i 1)
(3.131)
1
t J inv
j3
,k
2
hk ak k t tVi 2k ,(i 1) k t tVi1k ,(i 1)
After some rearrangement Eqn. (3.154) yields:
154
ui h h ,k
U ik k t J inv
j1
,k
k t J inv
j2
xj
t
r s
ak h h ,k t
k Vi 2k ,( i 1) t k t J inv
t t
j1
,k
k t J inv
j2 J inv
j3
,k
hk
2 r s
a h h ,k
k k t tVi1k ,( i 1) t k t J inv
j1
,k
k t J inv
j2 t J inv
j3
,k
hk
2 r s
U ik t h kj k t t
gi1k ,(i 1) t G kj k t t
gi2k ,( i 1) t G kj
(3.132)
ak
t t
gi1k ,(i 1) t t
Vi 2k ,( i 1)
2
ak
t t
gi2k ,(i 1) t t
Vi1k ,( i 1)
2
At this place, we can derive final expression to compute linear and nonlinear strans increments.
Linear strains t eij( i ) are calculated as follows:
t e11( i )
(i ) u1( i )
t e22
t e33 (i ) ...
t e (i ) ( i ) B1L0 B1L1 ... B kL0 B kL1 ... B n B n u (ki )
L0 L1
(3.133)
2 t e12
2 e(i ) ...
t 23 u ( i )
(i ) n
2 t e13
t h1k 0 0 t t
g11k ,(i 1) t G1k t t
g12k ,( i 1) t G1k
0
t
h2k 0 g 12k ,(i 1) t G2k
t t
g 22k ,( i 1) t G2k
t t
0 0 t k
h3 t t 1k ,( i 1) t k
g3 G3 t t 2k ,( i 1) t k
g3 G3
B kL0 t k t t 1k ,( i 1) t k t k
(3.134)
h2
t k
h1 0 g1 G2 t t g 12k ,(i 1) t G1k t t 2k ,( i 1) t k
g1 G2 t t g 22k ,(i 1) G1
0 t
h3k t
h2k t t
g 12k ,( i 1) t G3k t t g31k ,(i 1) t G2k t t
g 22k ,( i 1) t G3k t t g32k ,(i 1) t k
G2
t k t t 1k ,( i 1) t k
h3 0 t k
h1 g1 G3 t t g31k ,(i 1) t G1k t t 2k ,( i 1) t k
g1 G3 t t g32k ,(i 1) t k
G1
where
T
u (ki ) U1k ,(i ) , U 2k ,(i ) , U 3k ,( i ) , k ,(i ) , k ,( i ) at node k.
The second part of B kL , i.e B kL1 , is derived from 2 t eij1 B kL1 u (ki ) , at node k:
1 t t (i 1) k t k
2
t u m ,i
Um hj k t t
g 1mk ,(i 1) t G kj k t t
g m2k ,(i 1) t G kj
1
2
t tt um( i,j1) U mk t hik k t t
g 1mk ,( i 1) t Gik k t t
g m2k ,(i 1) t Gik
1 k t t (i 1) t k
2
U m t um,i h j k t tt um(i,i1) t t g 1mk ,(i 1) t G kj k t tt um(i,i1) t t g m2k ,(i 1) t G kj
1
U mk t tt um(i,j1) t hik k t tt um(i,j1) t t g 1mk ,(i 1) t Gik k t tt um(i,j1) t t g m2k ,( i 1) t Gik
2
Introducing
t t um(i 1)
lmj(i 1) t t ( i 1)
u
t m, j
t x j
1k j,(i 1) t t g 1mk ,(i 1) lmj(i 1) (3.135)
m
we can write
1 k (i 1) t k
e
1
t ij
2
U m lmi h j k 1ki,(i 1) t G kj k 2k i,(i 1) t G kj
(3.136)
1
U mk lmj(i 1) t hik k 1k j,(i 1) t Gik k 2k ,(j i 1) t Gik
2
l11(i 1) t h1k l21(i 1) t h1k l31(i 1) t h1k k ,( i 1) t k
11 G1 k ,( i 1) t k
21 G1
l12(i 1) t h2k l22(i 1) t h2k l32(i 1) t h2k k ,( i 1) t k
12 G2 k ,( i 1) t k
22 G2
l13(i 1) t h3k l23(i 1) t h3k l33(i 1) t h3k 13k ,( i 1) t k
G3 k ,( i 1) t k
23 G3
B kL1 (i 1) t k
l11 h2 l12(i 1) t h1k
t k
l21(i 1) h2 l22(i 1) t h1k
t k
l31(i 1) h2 l32(i 1) t h1k
t k k ,( i 1) t k
11 G2 12 k ,( i 1) t k
G1 k ,( i 1) t k
21 G2 22 k ,( i 1)
G1
l (i 1) t
h3k l13(i 1) t h2k l22(i 1) t h3k l23(i 1) t h2k l32(i 1) t h3k l33(i 1) t h2k k ,( i 1) t k
12 k ,( i 1) t k
G3 13 G2 k ,( i 1) t k
22 k ,( i 1)
G3 23 t k
G2
12(i 1)
l11
t
h3k l13(i 1) t h1k l21(i 1) t h3k l23(i 1) t h1k l31(i 1) t h3k l33(i 1) t h1k k ,( i 1) t k
11 k ,( i 1) t k
G3 13 G1 k ,( i 1) t k
21 k ,( i 1)
G3 23 t k
G1
(3.137)
156
S11
0 S11
0 0 S11 SYM
S12 0 0 S 22
t t
t S ( i 1)
0 S12 0 0 S 22 (3.138)
0 0 S12 0 0 S 22
S 0 0 S 23 0 0 S33
13
0 S13 0 0 S 23 0 0 S33
0 0 S13 0 0 S 23 0 0 S33
Then matrix B NL B1NL ... B kNL ... BnNL is composed so that (at a node k)
ij sij ij u (ki ) B kNL
T T t t
t S ( i 1) B kNL u (ki ) (3.139)
where states for variation of the following entity. It can be shown that the matrix B NL can be
set in the following shape:
t h1k 0 0 t t
g11k ,(i 1) t G1k t t
g12k ,( i 1) t G1k
t t 2k ,( i 1) t k
0
t
h1k 0 t t
g 12k ,(i 1) t G1k g2 G1
0 t t 2k ,( i 1) t k
0 t
h1k t t
g31k ,(i 1) t G1k g3 G1
t k t t 2k ,( i 1) t k
h2 0 0 t t
g11k ,(i 1) t G2k g1 G2
t t 2k ,( i 1) t k
B kNL 0 t
h2k 0 t t
g 12k ,(i 1) t G2k g2 G2 (3.140)
t k t t 1k ,( i 1) t k t t 2k ,( i 1) t k
0 0 h2 g 3 G2 g3 G2
t hk t t 1k ,( i 1) t k t t 2k ,( i 1) t k
0 0 g G g1 G3
3 1 3
t t 1k ,( i 1) t t 2k ,( i 1) t k
0 t
hk
3 0 g 2
t
Gk
3 g2 G3
t k t t 1k ,( i 1) t k t t 2k ,( i 1) t k
0 0 h3 g 3 G3 g3 G3
Having the matrices (3.134), (3.137), (3.138) and (3.140) these are used to compute the element's
stiffness matrix, mass matrix, element loads etc. in exactly the same way as it is done for other
Atena's element.
A family of 3D isoparametric shell elements is presented, see the figure below. Their properties
lie between degenerated Ahmad shell elements from Section 3.12 and full 3D brick elements
from Section 3.5.
Shape and kinematic behaviour resembles that of the shell's element. All points through the
shell's thickness remain located on a line passing thru the corresponding top and bottom nodes of
the shell, however unlike in the classical shell theory, their distance can change. As for degrees
of freedom, (DOFS), a typical 3D isoparametric shell element has 9 nodes at the top and nine
nodes at the botom surface, each of them having 3 DOFS, (i.e. 3 displacements). A similar 2D
shell element would feature 9 nodes located at the shell's midplane, each of them having 5
DOFS, (3 displacements plus 2 rotations).
The new elements use full 3D static equations. i.e. the elements consider all 6 components of 3D
stress and strain vector. Geometrical and material nonlinearity is supported. The governing
equations are calculated and integrated in material points. Gauss integration is used in shell's
158
plane direction, whilst layered concept is employed throughout the thickness of the shells, (i.e.
rectangur quadrature). As each layer can use different material model, some layers can be
employed for modelling of embedded reinforcement. The elements typically use 3 x 3 x
number_of_layers integration (i.e. material) points.
The elements are suitable for both shallow and deep shells and are extremely simple for use,
because they can be input and output as usual 3D solid hexahedral elements with 8, 20 or 27
nodes. Hence, these shells can be hadled with most 3D pre- and post-processors. They also use
standard 3D material models, element loads and other boundary conditions designed for
hexahedral elements.
The presented shell elements are particularly useful for structures that combine solid 3D
elements and shell elements, because they do not imply any additional shell kinematic constraint
that would harm an anjancent 3D solid elements. (Typical shell elements assume t 0 that
enforces the same displacements of the corresponding top and bottom nodes in direction of their
connecting line). They are designed for bent shells and to analyze these structures (with the same
accuracy) they require far less finite elements compared to a similar analysis using standard
hexahedral elements. On the other hand, the 3D behaviour of these elements involves a small
overhead, so that standard 2D shell elements (with only 5 stress/strain components per material
point) can perform in some cases slightly better. Nevertheless, the overhead is well paid off by
easy of use of the presented elements, their nice 3D visualization, simple connection to adjacent
3D solid parts of the structure etc. In addition, the hiearchical isoparametric space interpolation
(used for the presented 3D shell elements) ensures that finer and coarser meshes are easy to
connect. Of coarse, this feature must be supported by pre- and postprocessor being used.
The shell’s geometry at the configuration time t and t dt , (iteration (i-1) and (i)), is defined by:
1 t t k ,top 1 t t k ,bot
t
xi hk Xi Xi
2 2
t t ( i 1) 1 t t t k ,top (i 1) 1 t t t (3.141)
xi hk Xi X ik ,bot ( i 1)
2 2
t t (i) 1 t t t 1 t t t k ,bot (i )
xi hk X ik ,top (i ) Xi
2 2
Displacements at time t t (i 1) , i=1,2,3 for global axes x,y,z, at iteration (i 1) reads :
t t
ui( i 1) t t xi( i 1) t xi (3.142)
Substituting (3.141) into (3.142), i=1,2,3 for global axes x,y,z, we can derive
160
t t 1 t t t k ,top (i 1) 1 t t t k ,bot (i 1) 1 t t k ,top 1 t t k ,bot
ui(i 1) hk Xi Xi hk Xi Xi
2 2 2 2
1 t t t k ,top (i 1) t k ,top 1 t t t k ,bot (i 1) 1 t t k ,bot
hk Xi Xi
2
Xi Xi
2 2
1 t t t 1 t
hk U ik ,top (i 1) t t
U ik ,bot (i 1)
2 2
(3.143)
where
t t
U ik ,top (i 1) t t X ik ,top (i 1) t X ik , t t
U ik ,bot (i 1) t t X ik ,bot (i 1) t X ik
t t
Displacement increments within i-th iteration are calculated as ui(i ) t t xi( i ) t xi(i 1) :
1 t t t 1 t t t
ui(i ) hk U ik ,top (i ) U ik ,bot (i ) (3.144)
2 2
t t
where U ik ,top (i ) t t X ik ,top (i ) t X ik (i 1) , t t
U ik ,bot ( i ) t t X ik ,bot ( i ) t X ik ( i 1) . In the above X ik ,top
t t
and X ik ,bot is top and bottom nodal coordinate coordinate of node i. Similarly U ik ,top ( i 1) ,
t t
U ik ,bot (i 1) denotes displacements at the same node.
The elements are derived using Green-Lagrange strains and 2nd Piola Kirchhoff stresses. Total
Lagrangian formulation is employed, but after each load step we transform the analysed model
(and its stress and other tensors) to the coordinate system defined by the current shape of the
model. (The standard Total Lagrangian formulation calculates all with respect to the original
coordinate system without any transformation; Updated Lagrangian formulation carries all the
transformation each transformation, BATHE(1982). )
The shell's total strains at time t t , i-th iteration, are calculated: ( i,j=1..3 for axis x,y,z)
1 t t (i ) t t (i ) t t (i ) t t ( i )
t ij
t t ( i )
2
t ui, j t u j ,i t uk ,i t uk , j
1
t tt ui(,i j1) t ui(,i j) t tt u (ji,i 1) t u (ji,i) t tt uk( i,i 1) t uk( i,)i t tt uk( i,j1) t uk( i,)j
2
(3.145)
1
e
(i )
t ij
2
t ui(,ij) t u (ji,i) t tt uk(i,i 1) t uk(i,)j t tt uk(i,j1) t uk(i,)i
(3.146)
1
t (i )
ij t uk(i,)i t uk(i,)j
2
162
t xi hk 1 t t k ,top 1 t t k ,bot
Xi Xi
r r 2 2
xi hk 1 t t k ,top 1 t t k ,bot
t
Xi Xi (3.150)
s s 2 2
t xi hk
t
2
t
X ik ,top t X ik ,bot
t t
The above expressions are employed to obtain derivatives of (total) displacements ui(i 1) with
respect to r,s,t. They are needed to calculate strains (3.146).
where
hk 1 t t inv hk 1 t t inv hk bot t inv hk 1 t t inv hk 1 t t inv hk
hktop, j t J inv
j1 J j2 J j3 , hk , j J j1 J j2 J j3
r 2 s 2 2 r 2 s 2 2
At this place, we can derive final expression to compute linear and nonlinear strans increments.
Linear strains t eij( i ) are calculated as follows , see (3.146):
t e11( i )
(i ) u1( i )
t e22
t e33 (i ) ...
t e (i ) ( i ) B1L0 B1L1 ... B kL0 B kL1 ... B n B n u (ki )
L0 L1
(3.156)
2 t e12
2 e(i ) ...
t 23 u ( i )
(i ) n
2 t e13
t hktop,1 0 0 t hkbot,1 0 0
top top
0 t h k ,2 0 0 t hk ,2 0
0 0 h top
0 0 top
t hk ,3
B kL0 top top
t k ,3
top top (3.157)
t hk ,2 th k ,1 0 t hk ,2 t hk ,1 0
0 top top top top
th t hk ,2 0 t hk ,3 t hk ,2
top k ,3
t hk ,3 0 t hktop,1 top
t hk ,3 0 h top
t k ,1
T
where u (ki ) U1k ,top ( i ) , U 2k ,top ( i ) , U 3k ,top ( i ) , U1k ,bot ( i ) , U 2k ,bot ( i ) , U 3k ,bot ( i ) at node k.
Introducing
t t ui(i 1)
lij( i 1) t t ( i 1)
u
t i, j (3.158)
t x j
we can write
t t ( i 1)
u
t m ,i t um( i,) j t tt um(i,j1) t um(i,)i
lmi(i 1) hktop, j t t
U mk ,top ( i ) hkbot, j t t
U mk ,bot ( i ) lmj(i 1) hktop,i t t
U mk ,top (i ) hkbot,i U mk ,bot (i ) (3.159)
t t
t t
U mk ,top ( i ) lmi( i 1) hktop, j lmj( i 1) hktop,i t t U mk ,bot ( i ) lmi(i 1) hkbot, j lmj(i 1) hkbot,i
164
Finally, matrix B L1 yields
B kL1
l11(i 1)t hktop,1 ( i 1) top
l21 t hk ,1 l31( i 1) t hktop,1 l11( i 1) t hkbot,1 l21(i 1) t hkbot,1 l31( i 1)t hkbot,1
l12( i 1) hktop,2 ( i 1) top
l22 hk ,2 l32(i 1) hktop,2 l12( i 1) hkbot,2 l22(i 1) hkbot,2 l32( i 1) hkbot,2
l13( i 1) hktop,3 l23(i 1) hktop,3 l33(i 1) hktop,3 l13( i 1) hkbot,3 l23(i 1) hkbot,3 l33(i 1) hkbot,3
( i 1) top (i 1) top ( i 1) bot
l11 hk ,2 l12 hk ,1 l21(i 1) hktop,2 l22 ( i 1) top
hk ,1 l31(i 1) hktop,2 l32( i 1) hktop,1 l11(i 1) hkbot,2 l12(i 1) hkbot,1 ( i 1) bot
l21 hk ,2 l22(i 1) hkbot,1 ( i 1) bot
l31 hk ,2 l32 hk ,1
l ( i 1) htop l (i 1) htop l22(i 1) hktop,3 l23
( i 1) top
hk ,2 l32(i 1) hktop,3 l33( i 1) hktop,2 l12(i 1) hkbot,3 l13( i 1) hkbot,2 ( i 1) bot
l22 hk ,3 l23(i 1) hkbot,2 l32(i 1) hkbot,3 l33( i 1) hkbot,2
12( i 1) ktop,3 13(i 1) ktop,2 ( i 1) top ( i 1) top
l11 hk ,3 l13 hk ,1 l21 hk ,3 l23 hk ,1 l31(i 1) hktop,3 l33( i 1) hktop,1 l11(i 1) hkbot,3 l13(i 1) hkbot,1 l21 hk ,3 l23(i 1) hkbot,1
( i 1) bot
l31(i 1) hkbot,3 l33( i 1) hkbot,1
(3.160)
t t
Assembling stresses at time t t , iteration (i-1) into matrix t S (i 1) , participation of nonlinear
strains tij( i ) is, see (3.146)
1
t t
t Sij(i 1) tij(i ) t t
t Sij(i 1) tij(i ) t t
t Sij(i 1) t uk( i,)i t uk( i,)j
2
1 (i )
t tt Sij(i 1) t uk(i,)i t u k , j t u k ,i t u k , j
(i ) (i )
2
t tt Sij( i 1) t uk(i,)i t uk(i,)j
T (3.161)
u1(i ) u1(i )
... ...
u (ki ) B
T t t ( i 1)
NL
tS B u (ki )
NL
... ...
u ( i ) u (i )
n n
t hktop,1 0 0 t hkbot,1 0 0
top bot
0 t h k ,1 0 0 t h k ,1 0
0 0 top
t hk ,1 0 0 bot
t hk ,1
top
t hk ,2 0 0 t hkbot,2 0 0
B NL
0 hktop,2 0 0 hkbot,2 0
k
t t
0 0 t hktop,2 0 0 bot
t hk ,2
htop 0 0 bot
0 (3.162)
t k ,3 t hk ,3 0
0 t h top
k ,3 0 0 t h bot
k ,3 0
top bot
0 0 t h k ,3 0 0 t hk ,3
(3.163)
Using (3.147) and (3.148) it follows to present final expression for computation of space
derivatives of f ( x1 , x2 , x3 ) :
f h k
J inv
ji h k h k (r , s, t ) hˆ k (r , s ) h k (t )
k rj
F ;
xi
hˆ k k
k
f
J inv
ji kF
ˆk
h ( r , s ) h
(t )
J ji Fk
inv
rj
h k hˆ k
h
xi rj
r j
k k k
hˆ k k k h
hˆ k k k h
hˆ k k k h
J1i Fk
inv
h hˆ
J inv
F h hˆ J inv
F h ˆ
h
r k s
k t
2i 3i
r s t
k
hˆ k
ˆ
h k
h
J1inv k
inv k
inv ˆk
i Fk r h J 2i k F s h J 3i Fk h t
(3.164)
Having all the matrices and relationships above, the rest of derivation of the presented
isoparametric shell elements is straightforward. Simply use the matrices B L0 , B L1 , B NL and
t t ( i 1)
tS to calculate structural stiffness matrices t K L , t tt K NL
( i 1)
, vectors of nodal forces
t t
F (i 1) , and loads t t R as described in the Section Problem Discretisation Using Finite
Element Method earlier in this document.
166
3.16 Curvilinear Nonlinear 3D Isoparametric Layered Shell Wedge
Elements
This section describes wedge shell finite elements. Their properties and their derivation are much
the same as that for hexahedral shell finite elements CCIsoShellBrick<xxxxxxxx> ...
CCIsoShellBrick<xxxxxxxxxxxxxxxxxx> described in the previous chapter. The only diference
in that they feature wedge shape. Their geometry is depicted in the figure below
Depending on number of element nodes these finite elements call CCIsoShellWedge<xxxxxx>
... CCIsoShellWedge<xxxxxxxxxxxx>.
168
t
a2 Geom etry
t
b2
t 1
Vt
t
t
b3 a3
6
13
t
a1 1
Vs 14 5 Bric k nodes
7 16
t
s
b1 18 15 8
19
r 17 6
20 5 Beam 3D nodes
7
10 8
2
9
1 11
10
12 9
3 11 Beam 1D nodes
4 14
12
2
1 15 t
3
4
w
z
13 z 2
2 s
2
y r
x
x
y v
u Isoparam etric shape
Global c oord. system and elem ent dofs
Fig. 3-40 CCBeamNL element
vectors t Vt , tVs depicted in Fig. 3-40, in a cross section i, at time t, which define local
coordinate axis s,t. The symbols t ai , t bi refers to dimensions of the cross section i, time t; see
the figure, too.
Geometry of the beam at time t dt is defined in a similar way:
t dt t s
x hi t dt X i t ai t dtVi tx t bi t dtVi sx
2 2
t dt t t s s
y hi t dt Yi t ai t dtVi y t bi t dtVi y (3.166)
2 2
t dt t s
z hi t dt Z i t ai t dtVi tz t bi t dtVi sz
2 2
The element’s displacements at time t dt is calculated as follows:
t dt
u t dt x t x
t dt
v t dt y t y (3.167)
t dt
w t dt z t z
and displacement increments within a iteration:
t s
u hi U i t ai Vi tx t bi Vi sx
2 2
t t s s
v hi Vi t ai Vi y t bi Vi y (3.168)
2 2
t s
w hi Wi t ai Vi t z t bi Vi sz
2 2
In the above equation the vectors Vi t , Vi s are Vi t t dt Vi t t Vi t and Vi s t dt Vi s t Vi s and they are
approximated by
170
Vi tx t dt
Vi tz iy t dt
Vi y iz
t
t
Vi y t dt
Vi t x iz Vi tz ix
t dt
Vi tz Vi y ix
t dt t
Vi t x iy
t dt
(3.169)
Vi sx t dt
Vi sz iy t dt
Vi y iz
s
s
Vi y t dt
Vi sx iz t dt
Vi sz ix
Vi sz t dt
Vi y ix
s t dt
Vi sx iy
The parameters ix , iy , iz are rotations around the global axis, with respect to beginning of the
current load step. Note that (3.169) is valid only approximately.
172
t dt
The matrix t BNL . :
It is constructed in the way that
U 1
u V
x 1
W1
u x
y 1
1y
u z
z 1
U
v 2
x V2
v W
t dt
t BNL
2 t dt
BNL U
y 2x t
v y
2
z 2z
w
U 3
x V
w 3
y W3
x
w 3y (3.172)
z 3
z
3
The detailed expressions for calculating t dtt BNL are given in (3.175) and (3.176). The equations
are important because they present the way, how spatial derivatives of all the displacements are
calculated. The entries in t dtt BNL are thus used to setup also the matrix t dtt BL 0 and t dtt BL1 .
These matrices are computed as follows:
t dt t dt
t BL 0(1,i ) t BNL (1,i )
t dt t dt
t BL 0(2,i ) t BNL (5,i )
t dt t dt
t BL 0(3,i ) t BNL (9,i )
t dt t dt
(3.173)
t BL 0(4,i ) t BNL (4,i ) t dtt BNL (2,i )
t dt t dt
t BL 0(5,i ) t BNL (6,i ) t dtt BNL (8,i )
t dt t dt
t BL 0(6,i ) t BNL (7,i ) t dtt BNL (3,i )
174
t dt 1 hi
t BNL (1,1) J1,1
r
t dt
t BNL (1,2) 0
t dt
t BNL (1,3) 0
t dt
t BNL (1,4) 0
t dt 1hi t t t dt tz s t t dt sz 1 1 t t dt sz 1 1 hi t t dt tz
t BNL (1,5) J1,1 ai Vi bi Vi J1,2 hi bi Vi J1,3 ai Vi
r 2 2 2 2 r
t dt 1 hi t t t dt t y s t t dt s y 1 1 t t dt s y 1 1 hi t t dt t y
t BNL (1,6) J1,1 ai Vi bi Vi J1,2 hi bi Vi J1,3 ai Vi
r 2 2 2 2 r
t dt 1 hi
t BNL (2,1) J 2,1
r
t dt
B
t NL (2,2) 0
t dt
t BNL (2,3) 0
t dt
t BNL (2,4) 0
t dt 1hi t t t dt t z s t t dt sz 1 1 t t dt sz 1 1 hi t t dt tz
t BNL (2,5) J 2,1 ai Vi bi Vi J 2,2 hi bi Vi J 2,3 ai Vi
r 2 2 2 2 r
t dt 1 hi t t t dt t y s t t dt s y 1 1 t t dt s y 1 1 hi t t dt t y
t BNL (2,6) J 2,1 ai Vi bi Vi J 2,2 hi bi Vi J 2,3 ai Vi
r 2 2 2 2 r
t dt 1 hi
t BNL (3,1) J 3,1
r
t dt
t BNL (3,2) 0
t dt
t BNL (3,3) 0
t dt
t BNL (3,4) 0
t dt 1 hi t t t dt tz s t t dt sz 1 1 t t dt sz 1 1 hi t t dt tz
t BNL (3,5) J 3,1 ai Vi bi Vi J 3,2 hi bi Vi J 3,3 ai Vi
r 2 2 2 2 r
t dt 1 hi t t t dt t y s t t dt s y 1 1 t t dt s y 1 1 hi t t dt t y
t BNL (3,6) J 3,1 ai Vi bi Vi J 3,2 hi bi Vi J 3,3 ai Vi
r 2 2 2 2 r
t dt
t BNL (4,1) 0
t dt 1 hi
t BNL (4,2) J1,1
r
t dt
t BNL (4,3) 0
t dt 1 hi t t t dt tz s t t dt sx 1 1 t t dt sz 1 1 hi t t dt tz
t BNL (4,4) J1,1 ai Vi bi Vi J1,2 hi bi Vi J1,3 ai Vi
r 2 2 2 2 r
t dt
t BNL (4,5) 0
t dt 1 hi t t t dt tx s t t dt sx 1 1 t t dt sx 1 1 hi t t dt tx
t BNL (4,6) J1,1 ai Vi bi Vi J1,2 hi bi Vi J1,3 ai Vi
r 2 2 2 2 r
t dt
t BNL (5,1) 0
t dt 1 hi
t BNL (5,2) J 2,1
r
t dt
t BNL (5,3) 0
t dt 1 hi t t t dt tz s t t dt sx 1 1 t t dt sz 1 1 hi t t dt tz
t BNL (5,4) J 2,1 ai Vi bi Vi J 2,2 hi bi Vi J 2,3 ai Vi
r 2 2 2 2 r
t dt
t BNL (5,5) 0 (3.175)
t dt 1 hi t t t dt tx s t t dt sx 1 1 t t dt sx 1 1 hi t t dt t x
t BNL (5,6) J1,1 ai Vi bi Vi J1,2 hi bi Vi J1,3 ai Vi
r 2 2 2 2 r
176
t t
The stress matrix t Sij from (1.34) has he form:
t tt xx t t
t xy t t
t xz
t t
t yy t t
t yz
t t
t zz
t t
t xx t t
t xy
t t
t xz
t t
t Sij
t t
t t
(3.177)
t yy t yz
t t
t zz
xx
t xz
t t t t t t
sym.
t t xy
t t
t yy
t t
t yz
t zz
t t
As already mentioned, stress-strain relations are calculated in r,s,t coordinate system, hence we
need equations for their transformations from global x,y,z coordinate system to the isoparametric
system with r,s,t coordinates and vice versa.
Let us denote t dt T , t dtT transformation matrices for strain and stress transformation from
global to isoparametric coordinate system, so that:
t tt xx
t t
t rr
t t t yy
t t t dt t zz
t t
t rs T t t
t tt rt t xy
t t
t tt yz
t xz
(3.178)
t t
t xx
t t
yy
t tt rr t
t t t dt zz
t t
t rs T
t
t tt rt
t t
t xy
t t
yz
t
t t
t xz
Then the transformation matrices are calculated by:
t dt V rx 2
V rz
ry 2 2 ry ry t dt
t dt t dt
V 2t dt V rx t dt
V 2t dt V V rz 2t dt V rx t dt
V rz
t dt ry t dt sy s ry t dt ry t dt sy
T t dt V rx t dtV sx t dt
V V t dt
V rz t dt
V sz t dt
V rx t dt
V y t dt V V sx t dt
V V sz t dt V rz t dt
V t dt rx t dt s z
V V t dt V rz t dtV sx
t dt V rx t dtV tx t dt
V
ry t dt
V
ty t dt
V rz t dt
V sz t dt
V rx t dt ty
V t dt
V
ry t dt
V tx t dt
V
ry t dt
V
tz t dt
V rz t dt
V
ty t dt rx t dt t z
V V Vt dt rz t dt t x
V
(3.179)
1 0.577350269189626 1.
2
2 -0.577350269189626 1.
1 0.774596669241483 0.555555555555556
3 2 0. 0.888888888888889
3 -0.774596669241483 0.555555555555556
1 0.861136311594053 0.347854845137454
2 0.339981043584856 0.652145154862546
4
3 -0.339981043584856 0.652145154862546
4 0.861136311594053 0.347854845137454
5 1 0.906179845938664 0.236926885056189
2 0.538469310105683 0.478628670499366
178
3 0. 0.568888888888889
4 -0.538469310105683 0.478628670499366
5 -0.906179845938664 0.236926885056189
1 0.932469514203152 0.171324492379170
2 0.661209386466265 0.360761573048139
3 0.238619186083197 0.467913934572691
6
4 -0.238619186083197 0.467913934572691
5 -0.661209386466265 0.360761573048139
6 -0.932469514203152 0.171324492379170
Although the 2-nodes integration may be sufficient, (with respect to the quadratic displacement
transformation), a higher order integration scheme will yield better result in a case of high
material nonlinearity and/or in a case of a very curved beam geometry.
As for integration within the cross-section, i.e. in s,t coordinates, trapezoidal quadrature is used.
The element cross-section is subdivided into ns , nt “strips” as depicted in the following figure.
t
dtnt
2 s
dt2
dt1
ds1 dsns ind ividual “weight“
ds2 2 and m aterial
Geometry
t
Vt
1
6
13
Vs
1
14 5 Brick nodes
7 16
s
18 15 8
19
r 17 6
20 5 Beam 3D nodes
7
10 8
2
9
1 11
10
12 9
3 11 4
12
2
1 t
3
4
w
z
z 2
2 s
2
y r
x
x
y v
u Isoparametric shape
Global coord. system and element dofs
Fig. 3-42 CCIsoBeamBrick12_3D and CCIsoBeamBrick8_3D elements
Shape of cross section can be any quadrilateral, i.e. it need not be only a rectangle as depicted
above. The elements are particularly useful for analyses of structures, where beam elements must
be combined with 3D solid and/or shell elements.
180
Derivation of the element is much the same as that for CCIsoShell element, i.e. Equations
(3.142) and (3.144) thru (3.164) remain valid. Geometry and displacement approximation
(3.143) is replaced by:
1 s t k , front 1 s t k ,back 1 t t k ,top 1 t t k ,bot
t
xi hk Xi Xi Xi Xi
2 2 2 2
t t ( i 1) 1 s t t k , front ( i 1) 1 s t t k ,back (i 1) 1 t t t k ,top (i 1) 1 t t t
xi hk Xi Xi Xi X ik ,bot (i 1) (3.183)
2 2 2 2
t t (i) 1 ds t t 1 s t t 1 t t t 1 t t t
xi hk X ik , front (i ) X ik ,back (i ) X ik ,top (i ) X ik ,bot (i )
2 2 2 2
hk hk (r ) are 1D interpolation functions, see the interpolation function for CCIsoTruss
elements. The same notation is used for CCIsoShell Elements.
The element is calculated in integration points, (i.e. material points) that are located similar to
CCBeamNL_3D elements, refer to Fig. 3-41. The element can use any 3D material model.
Different materials can be specified for each material points, (or points in cross section). Some
of them can be used for modelling of embedded reinforcement. (Btw. discrete reinforcement can
be employed, too). The elements support both material and geometric nonlinearity.
3 t
w
z
1 z 2
2 s
2
y r
x
x
y v
u Isoparametric shape
If we need to connect a node k of the 1D beam element (having 3 displacements and 3 rotations
degrees of freedom) to an ambient solid element (having 3 displacement dofs per node), for each
beam node we proceed as follows:
Using (3.168) and (3.169) the beam's displacement are,
t
u hk U k t ak
2
t dt
Vktz ky t dt t
Vk y kz
st
2
bk t dt
Vksz ky t dt s
Vk y kz
t st
v hk Vk t ak
2
Vktx kz
t dt t dt
Vktz kx
2
bk t dt
Vksx kz t dt
Vksz kx
(3.184)
t
w hk Wk t ak
2
t dt
Vk y kx
t t dt
Vktx ky
st
2
bk t dt
Vk y kx
s t dt
Vksx ky
Writing (3.184) for local coordinate system and dt 0 , and noting that
t
V [0, 0,1], V
tL t sL
[0,1, 0], k 1... we derive
t s
u L hk U kL t ak kyL t bk kzL
2 2
t
v L hk VkL t ak kxL (3.185)
2
s
wL hk WkL t bk kxL
2
182
Displacement difference of nodes top-bottom, i.e. s 0, t 1 , hk 1
yL 1 wL ,top bot
k t a 0 0 0 0 0 L ,right left (3.189)
kzL k u
1 v L ,right left
0 0 0 t 0 0 L ,right left
bk w
Transforming o global coordinate system, ( TG 2 L is transformation matrix from global to local
CS):
1 1 u top bot
0 t
0 0 0 top bot
2 ak 2t b v
k
x
y 1 T 0 wtop bot
k TG 2 L t a 0 G2L
T
0 0 0 0 (3.190)
k 0 TG 2 L u right left
kz
1 v right left
0 0 0 t 0 0 right left
bk w
and
u right left hhlrightUU right
l
hhlleft
v right bot hhlright VVright
l
hhlleftVVleft
l
(3.192)
wright bot hhlrightWWright
l
hlleftWWleft
l
where hhltop , hhlbot , hhlright , hhlleft are the solid's interpolation functions at locations corresponding
to point at beam's node i cross section at [ s, t ] equal to [0,1], [0, 1], [1, 0], [1, 0] , respectively.
l l l l
UU top ,UU bot ,UU right ,UU left are corresponding nodal displacements of the solid element.
The beam displacement dofs U ik are connected into solid elements as usually, i.e. using master-
slave boundary conditions. For rotational dofs kx , ky , kz at node k we can write
UU1top
top
VV1
WW1top
u top hh1top 0 0 hhNtop 0 0
top
v 0 hh1top 0 ... 0 hhNtop 0 ...
(3.193)
wtop 0 0 hh1top 0 0 hhNtop
UU top
N
VVN
top
top
WWN
utop hhtop UU top
Similarly
ubot hhbot UU bot
u right hh right UU right (3.194)
uleft hhleft UU left
UUtop
u u hh
top bot top
hh bot
0
0 UU bot
(3.195)
right left
u u 0 0 hh right hhleft UU right
left
UU
Finally, combining (3.190) and (3.195) yields the resulting equations for connecting the beam's
rotation and the ambient solid's displacements [UU ,VV , WW ]T
184
3.19.2 Integrated forces and moments for shells
Integrated forces for shells are computed as follows:
t/2
Nx' x ' x ' dz '
t / 2
t/2
Ny' y ' y ' dz '
t / 2
t/2
Nz ' z ' z ' dz '
t / 2
t/2
Qx ' z ' dz '
t / 2 x ' z '
t/2
Qx ' y ' x ' y ' dz ' (3.196)
t / 2
t/2
Qy ' z ' y ' z ' dz '
t / 2
t/2
M y' x ' x ' z dz '
t / 2
t/2
M x' y ' y ' ( z ) dz '
t / 2
t /2
Kx' y ' x ' y ' ( z ) dz '
t / 2
The actual values of the forces and moments are calculated by extrapolation of stresses from IPs
into finite element nodes, (please refer to Section "Extrapolation of Stress and Strain to Element
Nodes" in Chapter CONTINUUM GOVERNING EQUATIONS. The process is as follows:
Let us take an example of N x ' that is calculated by integration of x ' x ' thru element's thickness.
The stress x ' x ' at element nodes is extrapolated from stresses in IPs ˆ x ' x ' by
M ij hi h j dVe
Ve
where Ve stands for element volume. Using (3.196) and writing (3.197) for extrapolation within
shell mid-plane e , (i.e. integration trhu e instead of Ve ) we can write
PPxx , i
e t/2
t / 2
hi ˆ x ' x 'dz d e hi ˆ x ' x ' dVe
Ve
(3.198)
1
MM ij hi hj d e hi hj dVe
e Ve t
where t t (r , s ) is element thickness at r,s. The integration for extrapolation is carried out over
, because the forces and moments are the same trhu shell thickness. Note that h h (r , s ) is
e k k
Ve
hi hj dVe
e t/2
t / 2 i
h (r , s ) hj (r , s )dt d e
h (r , s ) h (r , s ) dt d
t/2
i j e
e t / 2
(3.199)
hi (r , s ) hj (r , s ) t (r , s) d e hi hj t d e
e e
1
hi h j dVe hi hj d e MM ij
Ve t e
The forces and moments act on the plane (x'y'). They are calculated similar way to (3.198),
1
however, MM ij hi h j dr h h dVe , where b h is area of the beam's cross section and
le Ve b h i j
le is element length.
186
load. Nevertheless, some elements are internally defined in a local coordinate system and it can
be employed for an element load definition, too. Location of such a local system, (if it exists) has
been described together with description of the associated finite element. For example, local
coordinate systems are defined for plane 3D isoparametric elements, shell and beam elements
etc. On the other hand, elements such as tetrahedrons, bricks and others are defined in directly in
global coordinate system and therefore a local element load is treated as if it were input as a
global element load.
An exception to the above are truss elements. Although they are defined in global coordinate
system, they do support local element load. Their local coordinate system (for element loading
only) is defined as follows:
local X axis points in direction of the truss element,
local Y axis is normal to local X axis and lies in the global XY plane,
its positive orientation is chosen so that the local X and local Y forms a right-hand (2D)
coordinate system in the plane defined by these local axes,
local Z axis is vector product of the local X and local Y axes, (for 3D case only).
if the truss is parallel to global z, then local X points in direction of global Z, local Y
coincides with global Y and local Z has opposite direction of the global X, (for 3D case
only).
2D 3D
ZG
YG N2
YG
YL
N1
N1 XL YL
YG
XG
XL N2
XG
Fig. 3-44 Local and global coordinate systems for truss element N1-N2, (e.g. loaded element edge)
Specification of a boundary load deserves slightly more attention. Firstly, it is applied only to an
element’s edge or an element’s surface, (see also the note below), as opposed to e.g. an element
body load that is for the whole element. Local coordinate system is thus defined by location of
the loaded edge or surface. Secondly, a boundary load definition must include a reference to a
selection, which contains nodes to be loaded. Their order in the list is irrelevant, as what really
matters is the order in which they appear in the element incidences. When processing a boundary
load, ATENA loops thru all element’s surfaces and edges, (in the order specified in the table
below) and checks appropriate incidental nodes. If the tested node is present in the list of loaded
boundary nodes, it is picked up and put into incidences of a new planar or line element. This
element is later used to process the boundary load. It is its local coordinate system, that is
(possibly) used to deal with local/global load transformations.
Note that only one surface or one edge of each element can be loaded in a single boundary load
specification. If more element’s surfaces or edges are to be loaded, use more boundary load
definitions. Violation of this rule causes an error report and skipping of the offending boundary
load.
188
3D edge load
2D edge load planar element
planar element
n1 ZG n1 XL ZL
YG XL
YL
YL n2 n3
n2
n3
YG
XG
XG
3D edge load
solid element 3D surface load
solid element
n1 X L Z L n1
ZG ZL
ZG
YL
YL
n3 XL
n2 n2 n3
n6 n6
n5 YG YG
XG n5
XG
Fig. 3-45 Examples of positioning local coordinate system used by surface and element load for 2D
and 3D elements
Transport analysis does not distinguish between local and global element loads. Hence, a local
element “load” is treated as being a global load. The actual load value is always scalar, (unlike
vectors in statics) and it is assumed positive for flow out of the element.
3.21 References
AHMAD, S., B. M. IRONS, ET AL. (1970). "Analysis of Thick and Thin Shell Structures by
Curved Finite Elements." International Journal of Numerical Methods in Engineering 2:
419-451.
BATHE, K.J.(1982), Finite Element Procedures In Engineering Analysis, Prentice-Hall, Inc.,
Englewood Cliffs, New Jersey 07632, ISBN 0-13-317305-4.
190
4 SOLUTION OF NONLINEAR EQUATIONS
The main objective of this chapter is to review methods for solution of a set of nonlinear
equations. Several methods, which are implemented in ATENA are described later in this
Chapter. However, all of them need to solve a set of linear algebraic equation in the form
Ax b (4.1)
where A, x , b stands for a global structural matrix and vectors of unknown variables and rhs of
the problem, respectively. Hence, this problem is discussed first.
d a11 a77
T
a22 a33 a44 a55 a66
u a12 a67
T
a13 a23 a13 a24 a34 a15 a25 a35 a45 a46 a56
l a21 a31 a32 a76
T
a13 a43 a43 a51 a52 a53 a54 a64 a65 (4.3)
p 0 1 3 5 9 11 12
T
For each column i of the matrix A the vector p stores location of ai (i 1) within the array u ,
resp. l . If A is symmetric, then u l and only l is stored. Note the a direct solver we have to
The vector a stores for each column of A first diagonal element, followed by all non-zero
elements, from the top to the bottom of the column. The vector c stores row index of each entry
in the vector a . r stores location of all diagonal elements aii within a appended by an artificial
pointer to an 1n 1 , where n dim( A) .
Both of the above equations are computed easily, because the involved matrices have triangular
pattern. Hence, the solution of (4.7) represents back substitution only. If A is symmetric, (which
is usually the case), then
192
U LT (4.8)
194
DBCG I ssds sbcg S,NS 4*(11)+8*(1+8*n)
LUBC I ssilus sbcg S,NS 4*(13+4*n+nl+nu)+8*(1
+8*n+nu+nl)
DCGS I ssds scgs S,NS 4*(11)+8*(1+8*n)
LUCS I ssilus scgs S,NS 4*(13+4*n+nl+nu)+8*(1
+8*n+nu+nl)
DOMN I ssds somn S,NS 4*(11)+8*(1+4*n+nsave
+3*n*(nsave+1))
LUOM I ssilus somn S,NS 4*(13+4*n+nu+nl)+8*(1
+nl+nu+4*n+nsave+3*n
*(nsave+1))
DGMR I ssds sgmres S,NS 4*(31)+8*(2+n+n*(nsav
e+6)+nsave*(nsave+3))
LUGM I ssilus sgmres S,NS 4*(33+4*n+nl+nu)+8*(2
+n+nu+nl+n*(nsave+6)+
nsave*(nsave+3))
In the above:
n is number of degree of freedom of the problem. nel is the number of nonzeros in the lower
triangle of the problem matrix (including the diagonal). nl and nu is the number of nonzeros in
the lower resp. upper triangle of the matrix (excluding the diagonal).
ssds Diagonal Scaling Preconditioner SLAP Set Up. Routine to compute the
inverse of the diagonal of a matrix stored in the SLAP Column format.
ssilus Incomplete LU Decomposition Preconditioner SLAP Set Up.Routine to
generate the incomplete LDU decomposition of a matrix. The unit lower
triangular factor L is stored by rows and the unit upper triangular factor U is
stored by columns. The inverse of the diagonal matrix D is stored. No fill in
is allowed.
ssics Incompl Cholesky Decomposition Preconditioner SLAP Set Up. Routine to
generate the Incomplete Cholesky decomposition, L*D*L-trans, of a
symmetric positive definite matrix, A, which is stored in SLAP Column
format. The unit lower triangular matrix L is stored by rows, and the inverse
of the diagonal matrix D is stored.
ssd2s Diagonal Scaling Preconditioner SLAP Normal Eqns Set Up. Routine to
compute the inverse of the diagonal of the matrix A*A'. Where A is stored in
SLAP-Column format.
As for the solution procedure, i.e. the latter of the two solution phases, the most commonly used
method is Conjugate gradient method (with incomplete Cholesky preconditioner) (Rektorys
1995). The flow of execution is as follows:
196
r 1 b A x1
z 1 M 1 r 1
r i zi
i 1 i 1
i
r z
p i z i i p i 1
(4.11)
r i zi
i
i
p Ap i
x i 1 x i i p i
r i 1 r i i Ap i
z i 1 M 1 r i 1
i i 1
This solution procedure is implemented in scg routine.
The iterative solvers in ATENA are based on SLAP package (Seager and Greenbaum 1988) that
where modified to fit into ATENA framework. The authors of the package refer to (Hageman
and Young 1981), where all of the implemented solution techniques are fully described.
4
Available starting from ATENA version 5.
198
Conjugate-Gradients method is applied. You can select the method using only one input
parameter.
where:
q is the vector of total applied joint loads,
The R.H.S. of (4.12) represents out-of-balance forces during a load increment, i.e. the total load
level after applying the loading increment minus internal forces at the end of the previous load
step. Generally, the stiffness matrix is deformation dependent, i.e. a function of p , but this is
usually neglected within a load increment in order to preserve linearity. In this case the stiffness
matrix is calculated based on the value of p pertaining to the level prior to the load increment.
The set of equations (4.12) is nonlinear because of the non-linear properties of the internal
forces:
f (kp ) kf ( p ) (4.13)
All the quantities for the (i-1)-th iteration have already been calculated during previous solution
steps. Now we solve for p i at load level q using:
pi pi 1 pi (4.16)
As pointed out earlier, equation(4.15) is nonlinear and therefore it is necessary to iterate until
some convergence criterion is satisfied. The following possibilities are supported in ATENA
( k marks k -th component of the specified vector):
piT pi
rel . disp
piT pi
( q f ( pi 1 ))T ( q f ( pi 1 ))
rel . force
f ( pi )T f ( pi )
(4.17)
piT ( q f ( pi 1 ))
rel .energy
piT f ( pi )
The first one checks the norm of deformation changes during the last iteration whereas the
second one checks the norm of the out-of-balance forces. The third one checks out-of-balance
energy and the fourth conditions checks out-of-balanced forces in terms of maximum
components (rather then Euclid norms). The values of the convergence limits are set by
default to 0.01 or can be changed by input command SET.
The concept of solution nonlinear equation set by Full Newton-Raphson method is depicted in
Fig. 4-1:
200
Loading
q
Loading increment
p p p Deformation
0 1 2
The modified Newton-Raphson method is shown in Fig. 4-2. Comparing Fig. 4-1 and Fig. 4-2 it
is apparent that the Modified Newton-Raphson method converges more slowly than the original
Full Newton-Raphson method. On the other hand a single iteration costs less computing time,
because it is necessary to assemble and eliminate the stiffness matrix only once. In practice a
careful balance of the two methods is usually adopted in order to produce the best performance
for a particular case. Usually, it is recommended to start a solution with the original Newton-
Raphson method and later, i.e. near extreme points, switch to the modified procedure to avoid
divergence.
Loading
q
Loading increment
p0 p1 p2 Deformation
p3 p4
i i 1 i 1 (4.23)
The notation is explained in Fig. 4-3. The matrix K can be recomputed for every iteration
(similar to Full Newton-Raphson method) or it can be fixed based on the 1st iteration for all
subsequent iterations (Modified Newton Raphson method). The vector q does not mean in this
case the total loading at the end of the step but only a reference loading "type". The actual
loading level is a multiple of this.
202
The scalar is an additional variable introduced by the Line-search method, which will be
discussed later. The scalar is used to accelerate solutions in cases of well-behaved load-
deformation relationships or to damp possible oscillations, if some convergence problems arose,
e.g. near bifurcation and extreme points.
q 0
0 q
q 1
1 q
q 2 g
q 3 1
R2
g0
Load increment
Loading R1
T 0
0
0
q start
p 1
p 2
p 3
00 11 2 2
p0 p1 p2 Deformation
Substituting (4.21) through (4.25) into (4.20), the deformation increment i 1 can be calculated
from:
K i 1 RHS i 1 i 1q gi 1 (4.26)
Hence:
It remains only to set the additional constraint for i 1 and i 1 and the whole algorithm is
defined. Thus compared to the Newton-Raphson methods in which we solve n dimensional non-
linear problem, the Arc-length method need to solve a (n + 2) dimensional problem, where the
first n unknowns correspond to deformations and the last two are i 1 and i 1 .
ni 1 i 1 i 1 (4.30)
where:
is scalar that relates dimensions of to size of deformation space,
start is a (n+1) dimensional vector similar to i 1 , however its (n + 1)-th coordinate equal
to start .
n1
n2
n3
t1 t3
t2
204
Ri 1 piT1 i 1
i 1 (4.33)
piT1 T 2 (i 1 start )
To obtain i 1 by (4.33) the residual Ri 1 must be defined. In fact, it also define type of Arc-
length constrain being used. The types supported in ATENA are described below.
n1
n2
n3
t1 t3
t2
The step length s and angle are depicted in Fig. 4.3-4. The norm of the vector ti 1 is
calculated using (4.29):
2
ti 1 piT1pi 1 2 ( i 1 start ) 2 (4.35)
Based on the similar triangles (see Fig. 4.4-), the following can be derived:
rl 1 tl 1
(4.37)
t 'l 1 s t 'l
s 2 ( ti ' s )
Ri 1 (4.38)
ti '
The vector ti ' 1 is calculated using (4.35). By substituting the above equations into (4.33) final
expression for i 1 is obtained.
From the above derivation it is clear that in practice we at first employ Normal Update Method
(Chapter 4.4.1) to solve for ti ' and ni'1 and thereafter we correct the i 1 in order to satisfy
the constraint ti 1 ti s .
206
||ri-1|| n’i-1
|| r’i-1|| = ||t’i - s||
||ti-1||
ti-1 ni-1
ti
s
||t’i||
s = step length
p
Fig. 4.4-5 Explicit orthogonal method.
This method is usually utilized to analyze geometrically nonlinear structures, particularly
stability problems. Its main feature is robustness and compared with the "classical" Crisfield
cylinder method (see below) it avoids the problem of the choice of the proper i 1 root (the
condition ti 1 ti s while expressing vector length analytically). As for convergence, the
method is comparable to the method 4.4.3, but has the advantage that it preserves the step length.
a3 2 ( i 1 start ) 2 2 i T1 i 1 s 2
Equation (4.41) has generally two roots i 1 and hence we must decide which of them to use.
There exist several strategies but ATENA chooses that root i 1 , for which cos( ti 1 , ti ) 0 (or
higher of them), i.e. direction of new increment as close as possible to direction of the previous
increment (within the same step).
ni 1
si 4 si 1 (4.44)
n
n
si si 1 (4.45)
ni 1
where
si and si 1 is Arc length step length in the current and the previous load increment,
respectively.
n and ni 1 is desired number of iterations and number of iterations in the previous step.
n is typically 5-6.
Hence:
208
( p ) p T
p p
0 p g ( p) g ( p )T 0 (4.47)
d p 0
g ( p0 ' ) g ( p0 ) g ( p0 ' ) g ( p0 )
g ( p0 ) g ( p0 ) || p0 p0 || g ( p0 )
|| p0 ' p0 || '
(4.48)
and using :
p p0
p (4.49)
The final expression for ' can be derived:
g ( p0 )T
' (4.50)
g ( p0 )T g ( p0 ' )T
Practical experience suggests that the value of parameter should be kept in interval < 0.1 – 5>.
4.6 Parameter
The parameter scales the deformation space p to the loading dimension . If 0 , the
solution for i 1 is searched on an area of a cylindrical shape of radius equal to step length
s (Crisfield method) and the axis normal to the p (deformation) space. The solution is the point
of intersection of this area and the line, defined by the energy gradients of structure and by the
applied load at point p . If 0 , the solution is carried out in the same way on ellipsoidal or
spherical space.
The higher value of , the higher "weight factor" for changes in loading space compared to
displacement increments.
ATENA currently supports the following formulae for setting and optimization of (for current
step j ). They are reviewed below.
req (4.51)
( p )
This value (due to nonlinearities) will not match req . Therefore, for step j we will modify j
as follows:
j req j 1
req
(4.53)
j 1
req req ( j 1 p )
req
j 1 j 1
j 1 j 1 j 1
j 1
j
j 1
( j 1 p )
The above optimisation process is initialized in the first step by assuming that
0 1, 0 1, ( j 1 p ) T , where T is displacement corresponding to master Arc-length
load increment defined earlier in this chapter. Hence
req req req
0 req T (4.54)
j 1 j 1
1
j 1 1
( j 1 p) T
The parameters j in all subsequent steps are calculated using (4.53). If ratio of displacements
changes ( j p ) to load changes ( j ) in the last load step increase, then the equation (4.54)
(4.55) increases in the current step, thereby puts higher „weight factor“ on loads compared to
displacements. Hence, the equation (4.54) tends to keep constant importance of loading space
irrespective of displacements. Note that the equation (4.54) corresponds to
BETA_FORCES_DISPLS_RATIO_CONSTANT.
The second supported strategy is different. In ATENA it is referred to as
BETA_RATIO_CONSTANT method and it tries to keep constant coefficients, whilst
managing the coefficients . Thus, it works in opposite way as compared to the first strategy
described above.
From (4.52) we can write for steps (j-1) and j
210
j 1 ( j 1 p )
j 1
j 1
j ( j p)
j
j
j ( j p) j 1 ( j 1 p)
j j 1
( j 1 p)
(4.56)
j j 1
j 1 ( j p)
j
j 1 j j j
and if we assume , then and the above equation yields
( j 1 p ) ( j p) j 1 j 1
( j 1 p)
j 1
j j 1 (4.57)
( j p)
j
( j 1 p )
j 1
If in subsequent steps changes, the procedure is trying to compensate for that by re-
( j p)
j
( p)
adjusting the coefficients . In other words, this strategy is trying to keep constant, (i.e.
relative importance of load vs. displacement spaces).
212
Fig. 4-5 Optimization of dofs numbering
4.8 References
BATHE, K.J.(1982), Finite Element Procedures In Engineering Analysis, Prentice-Hall, Inc.,
Englewood Cliffs, New Jersey 07632, ISBN 0-13-317305-4.
CRISFIELD, M.A. (1983) - An Arc-Length Method Including Line Search and Accelerations,
International Journal for Numerical Methods in Engineering, Vol.19, pp.1269-1289.
CUTHILL, E. and J. MCKEE (1969). Reducing the Bandwidth of Sparse Symmetric Matrices.
Proc. 24th Nat. Conf. ACM.
DAVIS, T., AMESTOY, P., DUFF, I.S (1995) - An Aproximate Minimum Degree Ordering
Algorithm, Comp. and Information Science Dept., University of Florida, Tech. Report
TR-94-039.
DUFF, I. S. and KOSTER, J. (1999) - The Design and Use of Algorithms for Permuting Large
Entries to the Diagonal of Sparse Matrices. SIAM J. Matrix Analysis and Applications,
20(4):889-901.
FELIPPA, C. (1966) - Refined Finite Element Analysis of Linear and Nonlinear
Two-Dimensional Structures, Ph.D. Dissertation, University of California, Engineering,
pp.41-50.
GIBBS, N. E., W. G. POOLE, et al. (1976). "An Algorithm for Reducing the Bandwidth and
Prole of a Sparse Matrix." SIAM Journal of Numerical Analysis 13(2).
KARYPIS, G. and KUMAR, V. (1998) - A Fast and High Quality Multilevel Scheme for
Partitioning Irregular Graphs. SIAM Journal on Scientific Computing, 20(1):359-392.
LIU, J.W.H. (1985) - Modification of the Minimum-Degree Algorithm by Multiple Elimination.
ACM Transactions on Mathematical Software, 11(2):141-153.
LI, X.S., DEMMEL, J, W. (1999) - A Scalable Sparse Direct Solver Using Static Pivoting. In
Proceeding of the 9th SIAM conference on Parallel Processing for Scientific Computing, San
Antonio, Texas, March 22-34.
214
RAMM, E. (1981) - Strategies for Tracing Non- linear Responses Near Limit Points, Non- linear
Finite Element Analysis in Structural Mechanics, (Eds. W.Wunderlich,E.Stein, K.J.Bathe)
REKTORYS, K. (1995). Přehled užité matematiky. Prague, Prometheus.
SLOAN, S. W. and M. F. RANDOLF (1983). "Automatic Element Reordening for Finite
Element Analysis with Frontal Solution Schemes." Int. Journal for Numerical Methods in Eng.
19: 1153-1181.
SEAGER, M. K. and A. GREENBAUM (1988). A SLAP for the Masses, Lawrence Livermore
National Laboratory
SCHENK, O., GARTNER, K. and FICHTNER, W. (2000) - Efficient Sparse LU Factorization
with Left-right Looking Strategy on Shared Memory Multiprocessors. BIT, 40(1):158-176.
SCHENK, O. and GARTNER, K. (2004) - On Fast Factorization Pivoting Methods for Sparse
Symmetric Indefinite Systems. Technical Report, Department of Computer Science, University
of Basel, submitted.
SONNEVELD, P. (1989) - CGS, a Fast Lanczos-Type Solver for Nonsymmetric Linear
Systems. SIAM Journal on Scientific and Statistical Computing, 10:36-52.
VONDRACEK, R. (2006) - The Use Of The Sparse Direct Solver In The Egineering
Applications Of The Finite Element Method. Theses for Ph.D. Czech Technical University,
Prague.
Fig. 5-1 Decomposition of stress history into stress steps (left) or impulses (right).
The sense of Stieltjes integral is given in the above figure.
Equation (4.58) has to be modified for the case of 2 and 3D analyses for practical analyses. This
is done below. It is important to note that (4.58) applies for any stress and strain history and it is
defined in incremental form. It means that at a particular time t stress at t t depends only on
current material state at time t and stress increment at time t t , i.e. d .
The final form of the above equations reads:
t B( ( ))
(t ) (t , ) B( ( )) ( ) d 0 (t ) (4.59)
t'
where:
(t ) = is stress vector at time t , (note the bar atop of a symbol indicates vector),
218
Notice the way the equation (4.59) is written. Long-term and short-term material behaviour is
separated. The former is encapsulated in compliance function (t , t ') , whereas the short-term
behaviour is comprised in the matrix B( ( )) . This assumption brings significant simplification
of the creep and shrinkage analysis and it is believed that for most practical analysis the induced
inaccuracy is acceptable.
Substituting t t ' t , t 0 into (4.59) and applying load increment (t ') (t ') (i.e.
loading from zero level) at time t ' , it can be derived
The matrix D(t ') corresponds to reciprocal value of the well known secant Young modulus
E (t ') in the case of 1D stress-strain conditions. In the case of plane stress conditions, the matrix
B( ( )) reads (4.62), etc.
1 0
B= 1 0 (4.62)
sym. 2(1 )
n
t t '
1 1
'(t , t ') 1 e
(4.63)
E (t ') 1 E (t ')
where :
= are so called retardation times,
n = number of approximation functions, i.e. this parameter is related to the input parameter
number of retardation times.
E (t ') = instant Young modulus at time t ' ,
(tr ) r r 1 r (4.65)
(tr ) r r 1 r (4.66)
n
1 1 1
1 , r (4.67)
r 1/2
E Er 1/2 1 E ,r 1/ 2
t
r
r 1 e (4.68)
tr
220
i
n
t
r 1 e *r 1 r0
(4.69)
1
(tr ) r r 1 r (4.70)
e
*
t r
*
E
r
r 1/2
r (4.71)
r r 1
Er 1/2
1
Er 1/2 E (tr ) E (tr 1 ) = constant average secant Young modulus at time incremenent
2
tr ,
1 1
E r 1/2
2
2
E (tr ) E (tr 1 ) Er Er 1 = constant average value of Dirichlet
coefficient E at tr ,
1
B r 1/2 B(tr ) B(tr 1 ) = average value of the matrix B at tr .
2
Equation (4.64) thru (4.71) defines all necessary relations to complete the creep and shrinkage
analysis in ATENA. Of course, they are supplemented by relations used by short-term material
constitutive model, i.e. equations for calculating the matrix B.
At each time increment, a typical short-term alike analysis is carried. Difference between the
short-term analysis and the described analysis of one step of the creep and shrinkage is that the
latter one uses especially adjusted Young modulus E r 1/2 and initial strain increments r to
account for creep and shrinkage. After each step these have to be updated. It involves mainly
update of and r . With these values a new E
r
r 1/2 is calculated and the next temporal
analysis step is carried out.
In the above m is number of retardation times per log(t ) unit, m 1 . By default this constant is
in ATENA set to 1. If required, a more detailed approximation is possible, i.e. any value m 1
can be used. In the program this parameter is input as a number of retardation times per time unit
in logarithmic scale. For a typical concrete creep law a certain optimal value can be determined
and it is independent of a structure being analyzed. Note however, that the value depends on the
choice of time units.
Example: If the retardation times parameter is set to 2, the creep law will be approximated by
two approximation points for the time interval between 0 - 1 day, two points for the interval 1 -
10 days, then two points for 10 - 100 days, etc.
Therefore, the proper values will depend on the choice of time units. If the time unit is a day, the
recommended value is 1 - 2.
Start time 1 must be chosen sufficiently low, so that Dirichlet series can account for processes
in very young concrete, right after its loading has been applied. As a default, ATENA uses the
expression 1 0.1 t ' .
t
n (4.73)
2
The above limits are applicable for the case, when the coefficients E (t ') of Dirichlet series in
(4.63) are calculated by Least-square method (Jendele and Phillips 1992).
ATENA also supports alternative way of calculation of the coefficients E (t ') of Dirichlet series
in (4.63). In this case, Inverse Laplace transformation (Bazant and Xi 1995) is used instead. This
method requires 1 0 , typically 1E-3 and
n t (4.74)
Comparing the above two approaches, it can be said that Least-square method yields
approximation of the compliance function at discrete times, whereby Inverse transformation is
based on continuous approach. In some cases Least-square method results in better convergence
behaviour, however it sometimes suffers from numerical problems during calculation due to ill-
222
posed problem for solution of E (t ') . It is left to experience and engineering judgment to decide,
which of the method is more appropriate for a particular solution.
Integration times or sample times tr are calculated in similar way. In this case, the times are
uniformly spread in log(t t ') time scale. They are generated starting from the 1st loading time
t ' . Hence , we can write
r i 1
tr ti 10 l
ti 1 ti , t ' ti (4.75)
where l 2 is number of time increments per unit of log(t t ') and ti 1 t ' 0.1 ti 0.1 days.
Each new major load increment or decrement causes the generation procedure (4.75) must start
again from small time increments. This parameter defines the number of time steps, the program
will use to integrate the structural behavior. Creep or other nonlinear effects will cause a
redistribution of stresses inside the structure. In order to properly capture such processes a
sufficiently small time steps are needed. Its definition depends on the type of the analyzed
structure as well as on the choice of time units. For typical reinforced concrete structures and for
the time unit being a day, it is recommended to set this parameter to 2. This will mean that for
each load interval longer then 1 day, two sub-steps will be added. For a load that is interval
longer then 10 days, 4 sub-steps will be added. For an interval longer than 100 days, it will be 6
sub-steps, etc.
The creep and shrinkage analysis in ATENA requires that the user set number of retardation
times m and number of time increments l per unit of log time, (unless the default values are
OK). He/she also specifies time span, i.e. 1 and n . Then, retardation times are generated, i.e.
an appropriate command is issued. It follows to set stop time of the analysis. Usual input data
describing structural shape, material etc. are given thereafter, however, there are three important
differences from time-independent analysis:
1. Material model for concrete contains data for long‐term as well as for short‐term material
model.
2. Step data must include information about time, at which the step is applied.
3. It is recommended to input data for all intended load time steps prior the steps are executed. It
helps the generation of integration (intermediate) times
Intermediate time steps, i.e. times tr as well retardation times are generated automatically. The
analysis proceeds until the stop time is reached. If no stop time is specified, it is assumed to be
time of the last load step. If the time span for retardation times does not covered step load times,
the solution is aborted, giving an appropriate error message.
The following data summarized input parameters for the supported models. Note that some
models allow improved prediction based on laboratory data. If it is the case, the model input the
corresponding experimentally measured values. Also, some model can account for material point
history of humidity h(t ) and temperature T (t ) . Again, a model supports this feature, if it can
input adequate data.
Table 5.5-1 : List of material parameters for creep and shrinkage prediction – definition and
description
224
Concrete. type Type of concrete according to ACI. Type 1 1
is Portland cement etc. Types 1,3 accepted
for static analysis, types 1-4 accepted for
transport analysis.
226
End of curing Time at beginning of drying, i.e. end of days 7
curing.
a , Autogenous shrinkage at infinity time, - 0
(typically negative!)
t ts
a a , (0.99 min(0.99, ha , ) tanh
a
228
Table 5.5-2: Input parameters needed by individual creep and shrinkage prediction models
Model name B3 B3- BP- CEB ACI CSN BP1 BP2 Gen EN MC
impr KX eral 1992 2010
Model No. 3 4 8 2 1 5 6 7 9 10 11
Concrete. Type x x x x x x x
Cement class x x
Aggregate x x
Thickness x x x x x x x x x x
S /V
Strength f cyl 28 x x x x x x x x x x
Strength x x x
f cyl 0,28
Fracture x x x
energy G f ,28
Strength ft 28 x x x
Young m. E28 x x x x x x x
Ambient x x x x x x x x x x
humid. h
Ratio ac x x x x x x
Ratio wc x x x x x x
Ratio as
Ratio sa x x
Ratio g s x x
Ratio sc x x
Shape factor x x x x x
Slump x
Air content x
Cement mass x
Concr. density x x x x x x
a , x
a x
ts x
ha , x
Current x x x x x x x x x xx x
I/D
time t
Load time t’ x x x x x x x x x xx x
Tot.water loss x x x
w
Water loss x x
w(t)
Shrink. x x x x x x x x xx x
0 (t )
I
Compl. x x x
(t , t ')
Humidity x x x xx x
h(t )
Hi
Temperat. x x x xx x
T (t )
Compl. x
(t , t ')
Shrink. x
0 (t )
Di
Strength x
f cyl (t )
The above parameter “Concrete type” actually referes to a cement type according to the ACI
classification. It used in the creep analysis. The following table brings description of widely
recognized cement types. Note that only types 1,3 are supported in Atena static analysis. The
230
transport analysis in Atena recognizes types 1-4. The remaining types are described just for
information.
Table 5.5-3: Cement types according to ACI classification
ATENA
Concrete Cement type Description
type
5
Air-entraining cements
6
Blended hydraulic cements produced by intimately and uniformly intergrinding or blending two or more types of
fine materials. The primary materials are portland cement, ground granulated blast furnace slag, fly ash, silica fume,
calcined clay, other pozzolans, hydrated lime, and pre-blended combinations of these materials. The letter “X”
stands for the percentage of supplementary cementitious material included in the blended cement. Type IS(X), can
include up to 95% ground granulated blast-furnace slag. Type IP(X) can include up to 40% pozzolans.
7
All portland and blended cements are hydraulic cements. "Hydraulic cement" is merely a broader term. ASTM C
1157, Performance Specification for Hydraulic Cements, is a performance specification that includes portland
cement, modified portland cement, and blended cements. ASTM C 1157 recognizes six types of hydraulic cements.
5.6 References
ACI_COMMITTEE_209 (1978). Prediction of Creep, Shrinkage and Temperature Effects in
Concrete Structures. Detroit, 2nd draft, ACI.
BATHE, K. J. (1982). Finite Element Procedures in Engineering Analysis. Englewood Cliffs,
New Jersey 07632, Prentice Hall, Inc.
BAZANT, Z. AND T. SPENCER (1973). "Dirichlet Series Creep Function for Aging Concrete."
ASCE Journal of Engineering and Mechanical Division: 367-387.
BAZANT, Z. P. (1988). Mathematical Modeling of Creep and Shrinkage of Concrete. New
York, John Wiley & Sons.
BAZANT, Z. P. AND S. BAWEJA, EDS. (1999). Creep and Shrinkage Predicition Model for
Analysis and design of Concrete Structures: Model B3. Creep and Shrinkage of Concrete,
ACI Special Publicatino.
BAZANT, Z. P. AND J. K. KIM (1991). "Improved Prediction Model for Time-Dependent
Deformation of Concrete: Part 1- Shrinkage." Materials and Structures 24: 327-345.
BAZANT, Z. P. AND J. K. KIM (1991). "Improved Prediction Model for Time-Dependent
Deformation of Concrete: Part 2- Basic Creep." Materials and Structures 24: 409-421.
BAZANT, Z. P. AND J. K. Kim (1991). "Improved Prediction Model for Time-Dependent
Deformation of Concrete: Part 3- Creep at Drying." Material and Structures 25: 21-28.
BAZANT, Z. P. AND J. K. KIM (1991). "Improved Prediction Model for Time-Dependent
Deformation of Concrete: Part 4- Temperature Effects." Material and Structures 25: 84-
94.
BAZANT, Z. P. AND L. PANULA (1978). "Practical Prediction of Time-dependent
Deformations of Concrete; Part 1: Shrinkage." Material and Structures 11 (65): 301-316.
BAZANT, Z. P. AND L. PANULA (1978). "Practical Prediction of Time-dependent
Deformations of Concrete; Part 3: Drying Creep." Material and Structures 11 (65): 415-
423.
BAZANT, Z. P. AND L. PANULA (1978). "Practical Prediction of Time-dependent
Deformations of Concrete; Part 4: Temperature Effect on Basic Creep." Material and
Structures 11 (66): 424-434.
BAZANT, Z. P. AND L. PANULA (1978). "Practical Prediction of Time-dependent
Deformations of Time-dependent Deformation of Concrete; Part 2: Basic Creep."
Material and Structures 11 (65): 317-328.
BAZANT, Z. P. AND L. PANULA (1978). Simplified Prediction of Concrete Creep and
Shrinkage from Strength and Mix. Struct. Engng. Report No. 78-10/6405. Evanston,
Illinoins, Northwestern University, Dep. of Civ. Engng.
BAZANT, Z. P. AND F. H. WHITTMAN (1982). Creep and shrinkage in Concrete Structures.
New York, John Wiley & Sons.
BAZANT, Z. P. AND Y. XI (1995). "Continous Retardation Spectrum for Solidification Theory
of Concrete Creep." Journal of Engineering Mechanics 121(2): 281-287.
232
BETON, C. E.-I. D. (1984). CEB Design Manual on Structural Effects on Time Dependent
Behaviour of Concrete. Saint Saphorin, Switzerland, Georgi Publishing Company.
HAGEMAN, L. AND D. YOUNG (1981). Applied Iterative Methods. New York, Academic
Press.
JENDELE, L. AND D. V. PHILLIPS (1992). "Finite Element Software for Creep and Shrinkage
in Concrete." Computer and Structures 45 (1): 113-126.
REKTORYS, K. (1995). Přehled užité matematiky. Prague, Prometheus.
SEAGER, M. K. AND A. GREENBAUM (1988). A SLAP for the Masses, Lawrence Livermore
National Laboratory.
where tc is the construction phase, ti initiation period, tp propagation period and tr post-repair
period.
We aim at predicting the initiation period, without going into propagation or post-repair phases.
Carbonation and chloride ingress are two leading mechanisms contributing to reinforcement
corrosion. Both of them are described further. The initiation phase ends with the beginning of
reinforcement corrosion. Fig. 6-1 brings more detailed description of initiation and propagation
phases and their relationship to concrete events. Prediction of initiation period represents a
preventive measure which is affected above all by concrete cover thickness, concrete
composition, and environment. It makes sense to change design in the beginning rather than
mitigating reinforcement corrosion later. Acceleration of carbonation and chloride ingress on
crack appearance is taken into account.
8
Not available in ATENA version 5.1 and older. Development/testing implementation of CARBONATION,
CHLORIDES, and ASR in version 5.3.
6.1 Carbonation
Carbonation depth of a sound (uncracked) concrete reads (Papadakis and Tsimas 2002)
2 De ,CO 2CO2
xc t A1 t (4.77)
0.218(C kP )
where xc is the carbonation depth, De,CO2 is the effective diffusivity for CO2, C is the Portland
cement content in kgm-3, k<0.3,1.0> is the efficiency factor of supplementary cementitious
material (SCM-slag, silica, fly ash), P is the amount of SCM in kgm-3, CO2 is the volume
fraction of CO2 in the atmosphere taken as 3.6e-4 and t is the time of exposure. The effective
diffusivity in m2s-1 is given by the empirical equation (Papadakis and Tsimas 2002)
3
(W 0.267 C kP ) /1000)
De ,CO 2 6.1106 (1 RH ) 2.2 (4.78)
C kP W
c 1000
where W is the water content in concrete in kgm-3, c is the cement density in kgm-3 assumed as
3150 kgm-3 and RH is the relative humidity of ambient air. Eqs. (4.77)(4.78) allow predicting
either carbonation depth or induction time of uncracked concrete.
Cracked concrete leads to faster carbonation. This acceleration is given in the form (Kwon and
Na 2011)
xc (t ) (2.816 w 1) A1 t (4.79)
where w is the crack width in mm, A1 is the carbonation velocity according to Eq.(4.77).
Eq. (4.79) allows computing carbonation depth or induction time. Note that crack 0.3 mm
increases carbonation depth by a factor of 2.54. This also means that induction time is 6.46 times
shorter compared to a sound concrete.
236
In reality, crack may grow during any service time. Thus, Eq. (4.79) needs to be recast to
incremental form. An increment of carbonation depth in a given time step t is evaluated from
the total derivative by differentiating Eq. (4.79)
xc t
2.816
wi 1 1 A1
t
2.816 A1 t
w (4.80)
2 ti 0.5 2 wi 0.5
where wi+1 is the crack width at the end of the time step, ti+0.5 is the mid-time. It is assumed that
nonzero w at a frozen time t has no effect on carbonation depth, thus the term w can be left
out. Eq. (4.80) allows predicting either carbonation depth or induction time of gradually cracking
concrete.
6.2 Chlorides
Implemented model for chloride ingress is based on (Kwon, Na et al. 2009). Let us consider 1D
transient problem of chloride ingress in concrete with initially free chloride content
x
C x, t CS 1 erf (4.81)
2 Dm t f ( w)t
where CS is the chloride content at surface in kgm-3, Dm is the averaged diffusion coefficient at
time t in mm2 s-1, x is the position from the surface in mm and f(w) gives acceleration by
cracking and equals to one for a crack-free concrete. Cs and C can be related to concrete volume
or to binder volume, however, the units must be kept consistently through the computation.
Diffusion coefficient D(t) is assumed to decrease over time t according to the power law
238
Fig. 6-2. Evolution of actual and mean diffusion coefficients for standard concrete, based on data
from (Kwon, Na et al. 2009).
Let us assume characteristic value Cs 10.3% of chlorides per binder for submerged concrete
without further reductions (Table 8.5 in DuraCrete). The critical level for corrosion is 1.85 % per
binder (Table 8.7 in DuraCrete). The concrete cover is taken as 100 mm. Computed induction
240
time according to Eq. (4.81) is summarized in Table 6.3-1. Crack width is considered since the
beginning of exposure.
Table 6.3-1 Induction time for chloride corrosion of submerged concrete, in dependence on original
crack width.
The corrosion rate for the carbonation depends on the corrosion current density icorr [µA/cm2],
which ranges between 0.1-10 (passive corrosion-high corrosion) and depends on the quality and
the relative humidity of the concrete (Page CL, 1992). This model predicts amount of corroded
steel during the whole propagation period tp. The corrosion rate is based on Faraday’s law
(Rodriguez, 1996), determined as follows:
xcorr (t ) 0.0116icorr t (88)
where xcorr is the average corrosion rate in the radial direction [m/year], icorr is corrosion current
density [µA/cm2] and t is calculated time after the end of induction period [years].
By integration of Eq. (1), it is obtained the corroded depth for 1D propagation:
t
xcorr (t ) 0.0116icorr t Rcorr dt (89)
tini
where xcorr is the total amount of corroded steel in radial direction [mm] and Rcorr is parameter,
depends on the type of corrosion [-]. For uniform corrosion (carbonation) Rcorr = 1, for pitting
corrosion (chlorides) Rcorr = <2; 4> according to (Gonzales at.al., 1995) or Rcorr = <4; 5.5>
according to ( Darmawan &, 2007).
Effective bar diameter for both types of corrosion is obtained from:
d t dini 2 xcorr (t ) (90)
where d(t) is evolution of bar diameter in time t, d ini is initial bar diameter [mm], ψ is
uncertainty factor of the model [-], mean value ψ = 1 and xcorr is the total amount of corroded
steel according to (2).
The corrosion rate for chlorides is more complicated because it is affected by concentration of
chlorides in the concrete. Calculation of corrosion current density was formulated by Liu and
Weyer’s model (Liu, Weyers, 1998):
3006
icorr 0.926 * exp 7.98 0.7771ln1.69Ct 0.000116 RC 2.24t 0.215 (91)
T
where icorr is corrosion current density [µA/cm2], Ct is total chloride content [kg/m3 of concrete]
on reinforcement which is determined from 1D nonstationary transport, T is temperature at the
depth of reinforcement [K] and Rc is ohmic resistance of the cover concrete [Ω] (Liu, 1996) and t
is time after initiation [years]:
RC exp8.03 0.549 ln1 1.69Ct (92)
The average corrosion rate in radial direction is determined further when plugging (93),(94) to
(1). The total amount of corroded steel in radial direction stems from (2) and the effective bar
diameter from (3).
Cracking of concrete cover for both carbonation and chlorides can be estimated from DuraCrete
model which provides realistic results (DuraCrete, 2000). The critical penetration depth of
corroded steel xcorr,cr is formulated as:
C
xcorr ,cr a1 a 2 a3 f t ,ch (95)
d ini
where parameter a1 is equal 7.44e-5 [m], parameter a2 is equal 7.30e-6 [m], a3 is
[-1.74e-5 m/MPa], C is cover thickness of concrete [m], dini initial bar diameter [m], ft,ch is
characteristic splitting tensile strength of concrete [MPa].
The critical penetration depth of corroded steel xcorr,sp for both carbonation and chlorides is
calculated from (DueaCrete, 2000) as:
wd w0
xcorr , sp xcorr , cr (96)
b
where parameter b depends on the position of the bar (for top reinforcement 8.6 µm/µm and
bottom 10.4 µm/ µm), wd is critical crack width for spalling (characteristic value 1 mm), w0 is
width of initial crack (known from previous ATENA computation) and xcorr,cr depth of corroded
steel at the time of cracking [m].
After spalling of concrete cover, corrosion of reinforcement takes place in direct contact with the
environment. To determine the rate of corrosion of reinforcement after spalling, Table 2 gives
rates of reinforcement corrosion (Spec-net, 2015).
242
Corrosivity zone (ISO 9223) Typical environment Corrosion rate for first year (µm/yr)
Category Description Mild steel Zinc
C1 Very low Dry indoors ≤1,3 ≤0,1
C2 Low Arid/Urban inland >1,3 a ≤25 >0,1 a ≤0,7
C3 Medium Coastal and >25 a ≤50 >0,7 a ≤2,1
industrial
C4 High Calm sea-shore >50 a ≤80 >2,1 a ≤4,2
C5 Very High Surf sea-shore >80 a ≤200 >4,2 a ≤8,4
CX Extreme Ocean/Off-shore >200 a ≤700 >8,4 a ≤25
6.5 References
244
7 TRANSPORT ANALYSIS
As pointed out in the previous section, creep material behaviour of concrete strongly depends on
moisture and temperature conditions. Some constitutive models for creep in ATENA can pay
regards to these factors and based on previously computed moisture and temperature histories
within the structure they can predict concrete behaviour more accurately. This section describes
a module called CCStructuresTransport that is used to calculate the histories. A more accurate
creep analysis then typically consists of two steps: firstly execute CCStructuresTransport
module and calculate the moisture and humidity histories of the structure and secondly execute
CCStructuresCreep module to carry out the actual static analysis. Of course, for both analyses
we have to prepare an appropriate model. Export/Import of the results between the modules is
already done by ATENA automatically.
To be exact, both the transport and static analysis should be executed simultaneously but as
moisture and temperature transport does not depend significantly on structural deformations, i.e.
coupling of the analyses is low, the implemented “staggered” solution yields sufficiently
accurate results.
The governing equations for moisture transport read (for representative volume REV] :
w ( we wn )
div( J w ) (5.1)
t t
where:
w is total water content defined as a ratio of weight of water at current time t to weight of
water and cement at time t0 0 in REV, [mass/mass], e.g. [kg/kg]
we , wn = stands for the amounts of free and fixed (i.e. bound) water contents, [mass/mass],
is gradient operator.
Note that in (5.2) only diffusion of water vapor is considered. Moisture advection is negligible.
The equations (5.1) and (5.2) can be also written as being dependent on w or relative moisture h .
A relationship between h and w is given by
w w( h ) (5.3)
Using (5.3) Equation (5.2) can be written as follows
J w Dh h (5.4)
If hydration we want to add heat Qh (t ) , which expr esses amount of hydration heat
t t t t
J
within unit volume i.e Qh , 3 , Equation (5.5) changes to
m
T Qh
t
CT (T Tref ) Qh CT
t
t
div ( J T ) (5.6)
J
Heat flux JT , 2 is calculated by
m s
JT KT grad (T ) (5.7)
and KT stands for heat conductivity, e.g. [J/(day.m.K)].
Note that Equation (5.5) accounts for heat transport due to conduction only. Heat advection is
negligible. In (5.5) we can also neglect hydration heat, because in large times its impact for heat
transfer is small. On the other hand, we cannot neglect concrete moisture consumption due to
hydration process. According to (Bazant and Thonguthai 1978; Bazant 1986) hydration water
content wh can be calculated by:
1
t 3
wn wh 0.21 c e (5.8)
e te
where
e = 23 days, te is equivalent hydration time in water at temperature 25 0C that corresponds to
the same degree of hydration subject to real age, moisture and temperature conditions of the
material. The parameter c relates to the amount of cement and is calculated by(5.53). If
temperature ranges from 0 to 100 0C , te is computed by
246
te h T dt (5.9)
where dt is time increment after the mould has been removed and coefficients T , h are
calculated by
1
h (5.10)
1 (3.5 3.5h) 4
U 1 1
T exp h (5.11)
R T0 T
Uh
In the fraction the symbol U h stands for activation energy of hydration and R is gas
R
Uh
constant. According to (Bazant 1986) 2700 0 K . T, T0 are real and reference concrete
R
0
temperature expressed in K . The reference temperature is given by
T Qh
LHST CT (5.15)
t t
248
RHSh J w J h (5.16)
RHST J T (5.17)
The strip over an entity in the above equations means that the entity is vector. (Scalar entities do
not have the strip). The fluxes J w J h are identical, i.e. the subscript w indicates also moisture
phase. Using the above notation Equations (5.1) and (5.5) can be written as follows
LHSh div( RHSh )
(5.18)
LHST div( RHST )
The LHSh includes time derivative of moisture. It is computed using the following expressions:
wh wh (te )
te
h T (5.19)
t
T
w N T w; w N w (5.24)
T
T NT T ; T N T
where
h , w, T stands for vectors of the corresponding entities. The vectors have dimension n equal
to number of finite nodes of the problem.
N is vector of interpolation, (i.e. shape) functions,
N1 N 2 N n
x ...
x x
N N 2 N n
N 1
T
...
y y y
N1 N 2 N n
z ...
z z
Using (5.24) Equations (5.20) and (5.21) can be written in the form
250
h w T
LHSh chh N T chw N T chT N T ch 0
t t t
(5.25)
h w T
LHST cTh N T cTw N T cTT N T cT 0
t t t
and
(5.26)
RHST kTh N h kTw N w kTT N T kT 0
T T T
The resulting set of equations are solved iteratively using finite element method, see
(Zienkiewicz and Taylor 1989), (weak formulation, in which the shape functions N are used as
weight function):
N LHS
V
h
div( RHSh ) dV 0
(5.27)
N LHS
V
T
div( RHST ) dV 0
where V is volume of the analyzed structure. Each of the above equations represents a set of
equations with dimension equal to number of finite nodes n. Note that div( RHSh ) and
div ( RHST ) are scalars !
In the next derivation the two parts of (5.27) are dealt with separately.
h w T
N LHS dV N c
V
h
V
hh NT
t
chw N T
t
chT N T
t
ch 0 dV
h w
c
V
hh NN T dV chw NN T dV
t V t
... c
V
h0 NdV (5.28)
h w
cchh cchw ...cc h 0
t t
T h T w T T
V N LHST dV V
N cTh N cTw N cTT N cT 0 dV
t t t
(5.29)
h w
ccTh ccTw ...ccT 0
t t
(5.30)
ccTh cTh NN T dV ; ccTw cTw NN T dV ; ... ccT 0 cT 0 NdV
V V V
The second part of (5.27) are calculated using Green theorem (5.36):
N div( RHS ) dV N n
RHS h dS N RHSh dV
T
h s
V S V
N nsT
S
k N h k N w k N T k dS
hh
T
hw
T
hT
T
h0
(5.31)
kk h 0 N kh 0 dV
V (5.32)
...
kk T 0 N kT 0 dV
V
and also
252
J hh N nsT khh N dS
T
...
J h 0 N nsT kh 0 dS
S
JT 0 N nsT kT 0 dS
S
(5.33)
Using (5.28) to (5.33) the original governing equations (5.27) can be written as follows:
h w T
cchh cchw cchT cc h 0 kkhh h kkhw w kkhT T kk h 0
t t t
J hh h J hw w J hT T J h 0
(5.34)
h w T
ccTh ccTw ccTT ccT 0 kkTh h kkTw w kkTT T kk T 0
t t t
J Th h J Tw w J TT T J T 0
After sorting the unknown variables h , T by finite nodes into a single vector , Equation (5.34)
will read
cc kk cc 0 kk 0 J J 0 (5.35)
t
The right-hand side (5.35) is non-zero only for non-zero prescribed boundary conditions and
hence it has character of “load” vector in a static analysis.
In (5.31) we used Green theorem. It states:
(5.36)
where
where index (i )
indicates number of iteration and t t ( i )
is increment of the unknowns for time
t t and iteration ( i ) :
t t ( i ) -1 t t (i ) J
t t (i 1) K (5.40)
254
7.2.1 -parameter Crank Nicholson Scheme
This scheme comprises a number of well established integration procedures. It depends, what
value of the parameter is used. The set of equations (5.38) is solved for time t t ,
whereby the vector of unknown variables is calculated as a linear combination of the
corresponding vectors at time t and t t . Hence
t t
t (1 ) t t (5.41)
Depending on a particular value of the parameter we get the well known Euler implicit
integration (for =1), trapezoidal Crank Nicholson scheme (for =0.5), Galerkin integration
method (for =2/3) or even Euler explicit scheme (for =0), which is only conditionally stable.
Solution predictor:
t
t t
t t (5.42)
t
Solution corrector:
t t 1
t t t (5.43)
t t
J .
Using the above after some mathematical manipulation we derive final expressions for K,
These read:
= K 1 C
K
t
1
J J K t t 1 t C t t t (5.44)
t
1
K J
t t t prev
t
t
t t
2
t
(5.45)
t prev
2 t t prev t
where
index prev indicates that the entity comes from time preceding time t Note that we assume that all
entities from time t are already known and we solve for their values at time t t .
Solution corrector:
t t 2 t
t t t (5.46)
t t t
J K t t t t t t
n
2
n 1 n n 1
2
C 2t t t 2t
t t
n 1 n n 1
2 t
n 1
tn tn 1 tn t t tn
2 2 2
(5.48)
J t t t t
2 2
n n 1 n n 1
1
K J
where is a new damping factor. The factor is typically set to something in range 0.3...1
depending on current convergence behaviour of the problem.
256
Type of cement
w
Water-cement ratio wc
c
As already pointed out, the model does not account for aggregate, i.e. it predict moisture move
only in pores filled by water-cement paste.
w
The main entity of the model is water content w w(h, t , T , ) . It is defined as follows:
c
Gw
w (5.50)
Gw,0 G c
where
kg
Gw is water content in mortar at time t , 3 ,
m of morter
kg
Gw,0 is water content at time zero, 3 ,
m of morter
kg
Gc is amount of cement at time zero, 3 .
m of morter
Mortar here stands for mixture of water and cement. If concrete material is to be considered, then
w can be calculated by
Vconcrete
Gw w
Vmortar G
w (5.51)
V V w,0 G
G c
Gw,0 concrete G c concrete
Vmortar Vmortar
Vconcrete
where is ratio of total volume to (only) volume of mortar (i.e. water and cement) and G
Vmortar
are corresponding amounts of water and cements in concrete, (i.e. not only in
kg
mortar!) 3 .
m of concrete
w
The model itself already accounts for moisture used by hydration process. i.e. 0 . As a
t
result, wh according to (5.19) need not be implemented.
On the other hand, if moisture losses due to hydration are to be computed by the model based on
w
(5.19), it is important to fix 0 and to modify wh , so that it predicts “relative” moisture
t
content w used throughout whole derivation CCStructuresTransport. The original function for
wh was written for absolute weight of water and hence, for “relative” content of water Equations
(5.8) must be rewritten into
c te
3
0.21 G 1 1
e te c
G te 3 Gc te 3
wh 0.21 0.21 (5.52)
w,0 G c w,0 G
c
G G e te Gw,0 Gc e te
and the constant c from (5.8) becomes equal to
Gc c
G
c (5.53)
Gw,0 Gc G w,0 G
c
More detailed description of the model is beyond scope of this document and the reader is
referred to in (Xi, Bazant et al. 1993; Xi, Bazant et al. 1994).
CCTransportMaterial
CCTransport material is a simple constitutive law that allow users to enter laboratorily measured
moisture and heat characteristics. Refering to Equations (5.1) and (5.5) heat and moisture flow
governing equations can be written in the following general form:
Heat :
Q h T w
CTh CTT CTw CTt KTh grad (h) KTT grad (T ) KTw grad ( w) KTgrav
t t t t x
Moisture :
w h T w
Cwh CwT Cww Cwt Dwh grad (h) DwT grad (T ) Dww grad ( w) Dwgrav
t t t t x
(5.54)
The parameters CTh , CTT … K wgrav are calculated as:
258
CTh CTh0 fChTh (h) f CTTh (T ) fCtTh (t )
CTT CTT
0
fChTT (h) f CTTT (T ) fCtTT (t )
CTw CTw
0
f ChTw (h) fCTTw (T ) fCtTw (t )
CTt CTt0 f ChTt (h) f CTTt (T ) f CtTt (t )
Cwh Cwh
0
fChwh (h) fCTwh (T ) fCtwh (t )
CwT CwT
0
f ChwT (h) fCTwT (T ) fCtwT (t )
Cww Cww
0
fChww (h) fCTww (T ) fCtww (t )
Cwt Cwt0 fChwt (h) f CTwt (T ) f Ctwt (T )
CCTransportMaterialLevel7 material
CCTransport materialLevel7 is an extension of the above CCMaterialTransport material in the
way it automatically computes moisture and temperature capacity and conductivity/difussivity
incl. "sink" terms regarding hydration, (i.e. rate of hydration heat and mosture consumption
during connrete hydration). In terms of the above nomenclature this upper material level
calculates CTT , KTT , CTt , Cwh , Dwh , Cwt . As already mentined, the presented material adds on its
bottom level, i.e. CCMaterialTransport. All parameters and characterisctics from the bottom
level, (ie. those from CCMaterialTransport) can still be input and used. They typically serve for a
refinement/addition of parameters generated by the upper material level. The result from the
bottom an upper level are simply added to form final characteristics of the material model
CCTransportMaterialLevel7. Note that default values of CTT , KTT , CTt , Cwh , Dwh , Cwt in the bottom
level are by default set to zero.
Hydration heat and affinity hydration model
Consider hydrating cement under isothermal temperature 25oC a relative humidity h 1 . At this
temperature, the rate of hydration maturity factor , 0...1 can be expressed by chemical affinity
A A ( ) :
25 25
A25 (5.56)
t
where A stands for the chemical affinity, [ s 1 ], The expression already include coefficient
E
exp a . Hence A 25 is not normalized and refers to temperature 25oC. For different
RT
J
temperature it is replaced by A , see (5.60). R is gas constant 8314.41 , T is temperature,
kmol K
[K] and Ea is 40 kJ/mol. It is worthy to note the incorporation of the maturity method into (5.56)
. A characteristic time might be introduced to express an affinity A (Bernard, Ulm et al. 2003).
The affinity property can be obtained experimentally or analytically. Using experimental
approach, heat flow q (t ) that corresponds to the hydration heat Qh Qh (t ) is meassured in an
isothermal calorimetry.
Alternatively, the hydration material parameters are computed by an analytical micro-scale
model that accounts for the majority of underlying chemical reactions as well as topology of
cement grains (with the consequence to hydration kinetics). The solution stems from (Smilauer
and Bittnar 2006) and it employs discrete hydration model CEMHYD3D (Bentz 2005) allowing
to account for particle size distribution of cement, chemical composition of cement, temperature
and moisture history in concrete etc.
Qh
(5.57)
Qh , pot
1 Qh
A25 (5.58)
Qh , pot t t
Qh
where Qh, pot is potential hydration heat, [J/kg]. Hence the normalized heat flow under
Qh , pot
isothermal 25oC equals to chemical affinity A 25 .
Cervera et al. (Cervera, Oliver et al. 1999) developed an analytical form of the affinity which
was refined in (Gawin, Pesavento et al. 2006). A slightly modified formulation is proposed here
260
B Ea
A 25 B1 2 exp exp (5.59)
RT
where B1 ,[ s 1 ], B2 , [-] are coefficients to be calibrated, is the ultimate hydration degree, [-],
and represents microdiffusion of free water through formed hydrates, [-]. The parameters in
(5.59)express isothermal hydration at 25◦C.
When hydration proceeds under varying temperature, maturity principle expressed via Arrhenius
equation scales the a nity to arbitrary temperature T
E 1 1
AT A 25 exp a (5.60)
R 273.15 25 T
For example, simulating isothermal hydration at 35oC means scaling A 25 with a factor of 1.651
at a given time. This means that hydrating concrete for 10 hours at 35oC 35◦C releases the same
amount of heat as concrete hydrating for 16.51 hours under 25◦C. Note that setting Ea 0
ignores the e ect of temperature and proceeds the hydration under 25◦C.
Gawin et al. (Gawin, Pesavento et al. 2006), among others, added the effect of relative humidity.
The extension of (5.58) leads to
1 Qh
AT h (5.61)
Qh , pot t t
1
h (5.62)
1 a ah
4
where h h (h) accounts for the reduction of capillary moisture. h is relative humidity r,
(Bazant and Najjar 1972). a is material parameter, typically a 7.5 . Depending on curing
conditions is calculated as follows:
Sealed curing:
w/c
, 1 (5.63)
0.42
Saturated curing:
w/c
, 1 (5.64)
0.36
w / c is water-cement ratio.
Substituting (5.59) and (5.62) into (5.61) yields final equation to predict development of
hydration heat. As it is difficult to express function analytically (from (5.59), (5.61)), the
above equations are integrated numerically.
Heat capacity
The model assumes the following components of concrete: aggregate, filler, water and cement.
Total mass of concrete in one cubic meter results from individual masses of components:
mconcr maggregate m filler m paste
(5.65)
m paste mcement mwater
where mconcr is mass of concrete per a unit volume. Similarly for mass of aggregate maggregate ,
mass of filler m filler , mass of water mwater and mass of cement mcement . Corresponding volumes
are Vaggregate maggregate / aggregate , V filler m filler / filler etc. i stands for mass density of the phase
i. Having total volume Vconcr Vaggregate V filler Vwater Vcement , we can calculate phase fractions
f aggregate Vaggregate / Vconcr and similarly for the remaining phases.
Heat capacity and its evolution of cement paste (cement+water) was studied in (Bentz 2007) at
230C for w/c between 0.3 and 0.5. The capacity of fresh cement paste yields
where Cconcrete is concrete capacity (per unit volume) and akin for aggregate, filler and cement
paste. The last term, i.e. C paste depends also on degree of hydration and is calculated by
The heat capacity of structural concrete spans the range between 0.8 and 1.17 Jg-1K-1. A former
Czech standard CSN 731208 declares 840 and 870 Jkg-1K-1 for dry and saturated mature concrete,
respectively. Caggregate is approximately 840 Jkg-1K-1 for basalt and limestone, 790 Jkg-1K-1 for
granite, 800 Jkg-1K-1 for sand. Ccement is about 750 Jkg-1K-1 and Cwater is 4180 Jkg-1K-1 .
Heat conductivity
The thermal conductivity of cement paste was found to remain in the range 0.9-1.05 Wm-1K-1 for
arbitrary degree of hydration, for both sealed and saturated curing conditions, and for w/c from
0:3 to 0.4 (Bentz 2007). Water in the capillaries has the thermal conductivity 0.604 Wm-1K-1
(Bentz 2007). The thermal conductivity of hardened concrete varies between 0.85 and 3.5 Wm-
1 -1
K (Neville 1997) p.375, depending strongly on an aggregate type.
Thermal conductivity also depends on the saturation state of concrete. For example, a structural
concrete made from normal-weight aggregate with a unit mass of 2240 kg/m3 yields = 1.696
Wm-1K-1 for protected and 1.904 for weather exposed conditions (Neville 1997) , p. 376.
262
Figure 7-1. Thermal conductivity of concrete according to the Czech code CSN 731209.
Figure 7-1 summarizes thermal conductivities for ordinary concrete depending on concrete unit
mass and saturation conditions, according to (Neville 1997) and a former Czech standard CSN
731208. The latter considers 1.5 for a dry concrete and 1.7 Wm-1K-1 for a water-saturated
concrete.
Faria et al. (Faria, Azenha et al. 2006) applied the evolution of concrete conductivity with
regards to
0 1.0 0.248
where is the conductivity of fully hardened concrete, i.e. in infinite time.
wh wh
Ch , t (5.70)
t t
wh Qw, pot c ,[kg] (5.71)
264
Moisture capacity
Moisture content at unit volume w,[kgm -3 ] is calculated a simple expression
w wf
b 1 h (5.72)
bh
where w f ,[kgm-3 ] is the free water saturation and b is dimensionless approximation factor,
which must always be greater than one. It can be determined from the equilibrium water content
w80 at relative humidity h 0.8 by substituting the corresponding numerical values in equation
(5.72):
h( w f w80 )
b (5.73)
w f h w80
w w f b 1 b
Ch (5.74)
h b h
2
The above expression is applicable for analyses using reference unit volume. If reference unit
weight of the structure is preferred, then we employ moisture capacity C C / , kg/kg , where
is concrete density, kg/m 3 .
Moisture diffusion
The present model accounts for diffusivity mechanism of moisture transport. It is valid for dense
concrete, which has not mutually connected pores and moisture convection thru pores (being
kg
driven by water pressure) can be neglected. Hence, moisture flux qh , 2 is calculated by the
m s
kg
equation qh Dh h , where total moisture diffusivity Dh , is calculated as sum of water
ms
Dhw and watervapour Dhwv diffusivity:
w
Dhw Dww (5.76)
h
where water diffusivity Dww , m 2 / s is
w
1
wf
3.8 A21000
Dww (5.77)
w
2
f
As in the presented model relative humidity h is the primary variable used to analyze moisture
transport, Dpwv must be transformed to Dhwv . This is done by:
p ( psat h)
Dhwv D pwv D pwv D pwv psat (5.80)
h h
Any expression to calculate pressure of saturated water vapour can be used. The presented model
uses
aT
psat 611e T0 T
, Pa (5.81)
266
Table 7.3-1 Parameters of affinity hydration model used for CEM I.
The above table is based on fitting predicted results from CEMHYD3D analysis by (5.59), see
Table 7.3-5 and Figure 7-3. The simulations were carried out on CEMHYD3D’s microstructures
50 × 50 × 50 µm and with the activation energy 38.3 kJ/mol. Saturated curing conditions were
assumed, since sealed conditions will be obtained from coupling with moisture transport. Table
7.3-5 specifies input data for selected Portland cements.
Figure 7-3 Fit of selected cements to the affinity model, w/c = 0.4
The majority of concretes is produced from blended cements (CEM II - CEM V), hence it is
necessary to scale down Q pot by approximately 30 %. This is a common Portland clinker
substitution in the majority of blended cements in Europe.
There are other default parameters, which are not specified here: QW POT= 0.42, TH INIT = 0,
ALPHA INIT = 0, TEMPERATURE INCR MAX =0.1, H80 = 0.8, TEMP0 = 234.18, A WV =
17.08, TEMP0 ICE = 272.44 ,A WV ICE = 22.44
The parameter A ≈ 7.5 expresses hydration slow-down with regards to relative humidity. The
hydration practically stops at ≈ 0.8.
268
Table 7.3-3 Parameters specifying specific heat conductivity for concrete
components
Ready-mix concrete is assumed, which requires rather higher w/c due to workability and
pumping issues. The parameters CEMENT DENSITY, WATER DENSITY, AGGREGATE
DENSITY, FILLER DENSITY are provided in Table 7.3-2 in the units [kg/m 3 ].
Table 7.3-4 Approximate composition for major concrete classes used in EN206-1
270
7.4 Fire Element Boundary Load
When undertaking heat transfer calculations it is important that relevant thermal properties of
materials and heat transfer coefficients at the boundaries are defined for the entire temperature
interval of the load.
The convection and emissivity heat flux on a boundary exposed to fire is calculated as follows:
qn hc (Tg Tb ) r (Tg4 Tb4 ) (5.84)
where
= Stefan-Boltzmann constant [5.67x10-8 W/m2 K4],
Tg = absolute temperature of radiation source [K],
r = resulting emissivity factor of the radiation source and the heated surface [-],
Adiabatic boundary
Adiabatic boundary surface refers to a boundary surface, where no heat can pass in (and/or out)
the structure. Structural symmetry lines and areas are good example of this boundary conditions.
272
Nominal HV fire – Temperature of the heat source is calculated by (5.82) and T1 (unless
it is manually input as temp_g_ref) is set to 1100 [°C].
Modified HC fire – This definition is much the same as the above with the only
difference that default value for T1 is 1300 [°C].
Generic fire, (also refered to as User curve fire) - Temperature of the heat source is
assumed constant and is set value of temp_g_ref . If temp_g_ref is not inputed, then 1100
[°C] is used.
In any case, the generated (or directly inputed) curve for T (t ) can be additionally modified in
time by a user supplied function time_id. The function takes one parameter, which is time of the
fire and it specifies a coffecient, by which the originally generater (or inputed) boundary
conditions should be multiplied. Of course, load variation in space can be modified by coeff_x,
coeff_y coefficients etc. in the same way as for any other generated element load, (for more
details see Atena Input file manual).
The first part qh,1 includes moisture flux driven by gradient between ambient and surface relative
humidity hg and hb . hcw stands for moisture convection coefficient of the concrete-air interface.
The second part qh ,2 , accounts for moisture flux due evaporation driven by gradient of humidity
air ratio at the interface, i.e. xb and xg with the evaporation coefficient . By default
kg
(25 19v), 2 , where v is ambient air velocity, [m/s]. For more information see
m s
http://www.engineeringtoolbox.com/evaporation-water-surface-d_690.html.
The humidity air ratio, [-] is calculated as follows, (i reflects conditions in ambient air, i.e. i=g,
or in surface of the structure i.e. i=b):
9
Available starting from ATENA version 5.
pa p hi pvw.sat (5.88)
Here p stands for total air pressure, (typically normal air pressure p=101325Pa), hi is relative
humidity and pvw.sat is partial pressure of saturated water vapour at Ti, (see
http://en.wikipedia.org/wiki/Density_of_air)
7.5Ti
Ti 237.3
pwv , sat 610.78 10 (5.89)
The third part is moisture flux evaporated from concrete calculated by CEMSTONE, see
http://www.cemstone.com/concrete-evaporation-forecast-engineers.cfm. It yields nearly the
same values as provided by ACPA calculator, see
http://apps.acpa.org/apps/EvaporationCalculator.aspx. TCg , TCb is ambient and surface
temperature in Celsia.
274
kJ
hwe 2270, is assumed. More information available at
kg
http://www.engineeringtoolbox.com/evaporation-water-surface-d_690.html.
Both moisture and heat fluxes are typically computed using only their first or second part.
Therefore, the related ATENA input commands allows to read some boolean flags that specify,
which parts of the above fluxes should by accounted for and which should be skipped. For more
information refer to the ATENA input file manual.
7.6 References
BAZANT, Z. P. (1986). Mathematical Modelling of Moisture Diffusion and Pore Pressure,
Chapter 10. Concrete at High Temperature. Z. P. Bazant: 198-237.
BAZANT, Z. P. and W. THONGUTHAI (1978). Pore Pressure and Drying of Concrete at High
Temperature. Proceedings of the ASCE.
CELIA, M. A. and P. BINNING (1992). "A Mass Conservative Numerical Solution for Two-
Phase Flow in Porous Media with Application to Unsaturated Flow." Water Resour. Res
28(10): 2819-2828.
CELIA, M. A., T. BOULOUTAS, et al. (1990). "A General Mass-Conservative Numerical
Solution for the Unsaturated Flow Equations." Water Resour. Res 27(7): 1438-1496.
DIERSCH, H. J. G. and P. PERROCHET (1998). On the primary variable switching technique
for simulating unsaturated-saturated flows, http://www.wasy.de/eng/prodinfo/flow/
swpool/swpool.htm#fef_manuals.
HUGHES, J. R. (1983). Analysis of Transient Algorithms with Particular Reference to Stability
Behaviour. Computational Methods for Transient Analysis, Elsevier Science Publishers
B.V.
JENDELE, L. (2001). ATENA Pollutant Transport Module - Theory. Prague, Edited PIT, ISBN
80-902722-4-X.
JENDELE, L. and D. V. PHILLIPS (1992). "Finite Element Software for Creep and Shrinkage
in Concrete." Computer and Structures 45 (1): 113-126.
REKTORYS, K. (1995). Přehled užité matematiky. Prague, Prometheus.
SEAGER, M. K. and A. GREENBAUM (1988). A SLAP for the Masses, Lawrence Livermore
National Laboratory.
WOOD., W. L. (1990). Practical-Time Stepping Schemes. Oxford, Clarenton Press.
XI, Y., Z. P. BAZANT, et al. (1993). "Moisture Diffusion in Cementitious Materials, Adsorbtion
Isotherms." Advn. Cem. Bas. Mat. 1: 248-257.
XI, Y., Z. P. BAZANT, et al. (1994). "Moisture Diffusion in Cementitious Materials, Moisture
Capacity and Diffusivity." Advn. Cem. Bas. Mat. 1: 258-266.
276
8 DYNAMIC ANALYSIS
ATENA software support four methods to execute dynamic analyses. These are:
Newmark’s method,
Wilson
Modified Wilson .
Note that Hughes method with 0 reduces to Newmark’s method and Modified
Wilson is just an extension to Wilson .
The governing equations for dynamic analysis read:
Hughes method:
M ut t C 1 α u t t αu t K 1 α u t t αu t 1 α R t t α R t
Newmark method:
M ut t Cu t t Ku t t R t t
u t t 1
2βt 2 t 2 γ 1 2v _ t t 2u 2u γ t
t t
t 1 γ u
t
u
2 tβ 2 tβ
1 2u βt u t 2u t 2u 2u t t
t t 2 t t 2 t
ut t
2 t 2β
Substituting (5.92) into (5.91) and after some mathematical manipulation the requested
displacement increment at iteration l can be calculated:
u K inv
eff Reff (5.93)
where (for using structural damping C M M K K ) effective stiffness and RHS vector are:
K eff M M K K
(5.94)
Reff M ( M1 M2 ) K ( K1 K2 ) 0
The coefficients above are calculated using the following expressions. They are summarized (by
solution method):
278
Hughes method:
1 1
α M γ 1 1 α M γ 1 t
2 2
M α 1 M t 1 ut M u
β 2β β tβ
u t t αu u t t γ M
t t
t 2β βt
1 1
2 α 2 K γ 1 α K γ αu t t u t t γ K
K α 1 K tut K u t
β β βt
1 1
α M γ 1 1 α M γ 1 t
M1
2 2
α 1 M t 1 ut M u
β 2β β tβ
M2
u t t
αγt M γ t M 1
t β
2
1 1
2 α 2 K γ 1 α K γ
K1 α 1 K tut K u t
β β (5.95)
u t t 1 α γ K
K2
βt
0 R t t (1 α) R t α F t t (1+α)+F t α
Newmark method:
M γ 1 γ 1 t
M M t 1 ut M M u
2β 2β β tβ
u t t u γ M
t t
t 2β βt
K
K γ γ
K tut K K u t
u γ K t t
2β β βt
M γ 1 γ 1 t
M1 M t 1 ut M M u
2β 2β β tβ
u t t γt M 1
M2
t 2β
K γ γ
K1 K tut K K u t
2β β
u γ K
t t
K2
β t (5.96)
0 R t t
F t t
M
1
u 2
t
u
2 θ3 2 θ3 θ t 2 θ3
3u θt M 2
t t
M2
θ 2 t 2
1 3θ 2 2θ t K 1 θ3 θ 2 t 2 1 6θ 2 2 K 1 2θ3 2θ t
K
1
u
t
u t
2 θ 3
2 θ 3
2 θ 3
2 θ 3
3u K
t t
K2
θt
Wilson method :
R t t 1 1 t 1
0 2
2 3 R F t t 1 3 F t
θ θ θ θ
Modified Wilson method
(5.97)
1 1 F t t R t t
0 3 F t 3
θ θ θ θ
The parameters , are integration parameters used by Newmark and Hughes method.
Their value is essential for convergence of the this time marching scheme. It can be shown that
1 1 1 1
, corresponds to linear acceleration within the time step. Values ,
2 6 2 4
yields constant acceleration. The integration scheme is unconditionally stable, if
1 1 1 1
, 0.25( ) 2 and it is only conditionally stable for , 0.25( ) 2 provided
2 2 2 2
that the stability limit is fulfilled:
1 1 2
2
2
2 2 2
tcrit (5.98)
2
where is modal damping parameter.
The above defines the condition for time increment t for a linear conditionally stable case:
280
t 1
Tn 2
2 (5.99)
3
0.551
As for Wilson and Modified Wilson method they use parameter. Its value is 1 and
the scheme is unconditionally stable for 1.4 . It essentially specifies time, for which time we
calculate the governing equations (5.91), i.e t t . For 1 Wilson and Modified
Wilson method yield the same solution expressions and equations and these are also the same
1 1
as those for Newmark and Hughes methods with , , 0 .
2 6
Modified Wilson method assemble the governing equations for time t t . As a result, all
Von Neumann boundary conditions must be given for t t , (e.g. concentrated load, load by
MASS_ACCELERATION etc). It does not apply to Dirichlet boundary conditions that are (as
usually) input for t t , (e.q. prescribed displacement, acceleration etc.).
The fact that Modified Wilson method executes for t t also affects output/draw of results
in structural material points. Within iterations, (e.g. for monitors at iterations) they are
printed/drawn for t t . After the iterations process has been completed, they are
printed/drawn for t t as usually. Internal forces are always printed for t t and the same
for external forces.
As described above Modified Wilson method behaves in a bit nonstandard way. Particularly
input of R t t is unpractical. To alleviate these difficulties and inconvenience Atena also offers
Wilson method. Although it still solves the governing equations for time t t , it uses
several extrapolation, (e.g. R t t R t ( R t t R t ), F t t F t ( F t t F t ) ) so that it
suffices with R t t and F t t only. Consequently, it inputs all boundary conditions and
print/draw all result for t t akin to any other solution method for dynamic analysis. On the
other hand it is at price of accuracy because the extrapolation is linear whereby the loading and
internal forces is not!
Remind that for dynamic analysis concentrated forces, element body/boundary load.. is input in
incremental form and it is "cummulated" in the structure. The same applies for prescribed
displacements.
Prescribed velocities, accelerations... must be input as total load. MASS_ACCELERATION
must be also input in total values(and in each step it is also recalculated from scratch.
More details on the methods’ convergency can be found in (Hughes 1983) and (Wood. 1990).
where:
i is i-th structural eigenvector,
Using the fundamental properties of eigenmodes iT M 1, iT K i we can rewrite (5.101):
M K i 2ii (5.102)
Equations (5.100) introduces 2 parameters for damping and thus, if only 2 values of i are to be
used, they are directly substutited in (5.101), (resp. (5.102)) and solved for from this set of
equations.
However, in practice structural damping is more complicated and some sort of compromise must
be done. In this case structural damping properties are typically measured for more eigenmodes
and optimal values of coefficients M , K are calculated by least square method, i.e. we are
M K i 2ii
2
seeking minimum of the expression . It yields the following set of
equations
wi2 wi2i
M K wi2 2
i i2 i i i
(5.103)
M wi2 K wi2i2 2 wi2ii
i i i
There exists other assumptions to account for structural damping, however their use is typically
significantly more complex and more costly in terms of both RAM and CPU.
282
8.2 Spectral analysis
A proper selection of the solution time increment dt is essential for each dynamic analysis. If it
is too large, the computed results will suffer from unacceptable inaccuracies. We will probably
miss some important peaks in the loading history and the analysis as a whole may even diverge.
On the other hand, use of a too small value of dt will yield an analysis that is pointlessly
expensive in terms of execution time and its demands towards CPU/RAM resources. In addition,
its postprocessing is more laborious and prone to errors.
The spectral analysis described in this section is designed to assist the engineer in setting suitable
dt. The main idea of the procedure is based on approximation of the structural loading f (t ) by
Fourier series f FT (t ) , i.e. f (t ) f FT (t ) , refer e.g. to http://en.wikipedia.org/wiki/Fourier_series
. Both f (t ), f FT (t ) have one independent variable, which is structural time t .
a0 N 2 2
f FT (t ) an sin n t bn cos nt (5.104)
2 n 1 T T
where N denotes number of harmonics used for the approximation, n is harmonic-th id and
2 2
sin n t and cos n t are n-th approximation functions, (i.e. n-th harminics). Eqn.
T T
(5.104) is suitable for approximation of a function (e.g. f (t ) ) in interval t 0...T . Its
Fourier coefficient are calculated as follows, see
http://stelweb.asu.cas.cz/~slechta/fourier/fourier.html ,
http://www.mathstools.com/section/main/fourier_series_calculator#.VCFKkhZIpKI:
T
2
a0
f ( ) d
0
2
T
2
an f ( ) sin n d (5.105)
0 T
2
T
2
bn f ( ) cos n d
0 T
Now let us introduce a coefficient cn an2 bn2 and create a spectrum diagram of the loading.
2
For each harmonic from (5.104) plot a point, whose coordinates are n, cn' . Such a point
T
By default the FFT analysis uses full modal spectrum, i.e. n 1..N in (5.104). However, the
modal spectrum can be filtered, e.g. n n1..m1 , n2 ..m2 ,...nk ..mk ,...nL ..m1 . In this case only values n
from within the L intervals are used. This technique can be used to filter out some noise signals,
etc.
Let's take an example: Assume a simplified ElCentro accelerogram loading conditions, whose
acceleration in time are depicted by the green line in the figure below:
Let's approximate this function by Fourier series. In the first case we use 300 harmonics, i,e,
N 300 . The approximated accelerations are shown by the blue line, see the figure above. In the
second case, we use only 50 harmonics and the corresponding approximation function is drawn
by the red line. Plotting the functions in more detail, it can be seen that the approximation with
N 300 is fairly accurate whilst the approximation with N 50 is rather crude, see the figure
below. This conclusion is endorsed by calculated average relative absolute error of the
approximations. These are respectively 0.0314256 and 0.878789
284
The spectrum diagram below shows contribution of individual approximation harmonics. It
detects what harmonics are or are not important. Looking at the diagram we see that the highest
important harmonics is the one with log10(Tn ) 1 , i.e. Tmin 0.1 . Therefore, the recommended
T
solution time increment is dt min 0.01 . This dt should ensure reasonably accurate results
10
from dynamic analysis of a structure that is loaded by the investigated accelerogram.
286
9 EIGENVALUES AND EIGENVECTORS ANALYSIS
This section describes methods used by ATENA software to calculate structural eigenvalues
and eigenvectors. Good knowledge of eigenmodes of a structure is in many cases essential for
understanding its behaviour and selection of a method for its further analysis. It applies to
statics and particularly to dynamic analyses, in which case it helps choosing a proper time
increment in subsequent loading steps. It also help in avoiding or reducing oscillation of the
structure.
Ku 2Mu (8.1)
where
K, M is stiffness and mass matrix of structure,
u is vector of eigenvector’s nodal displacements,
is circular eigenfrequency
We are looking for a non-trivial solution, so that we solve for 2 that comes from
det(K 2 M ) 0 (8.2)
ui Ψci (8.3)
where
Ψ is matrix of base vectors k , k 1..m ,
uiT Kui
(ui ) (8.4)
uiT Mui
It can be proved that (ui ) converges from upper side to the corresponding circular
frequency i2 . The condition of minimum of (ui ) yields:
(ui )
0, k 1..m (8.5)
ci ,k
If we introduce
A Ψ T KΨ , B Ψ T MΨ (8.6)
the condition (8.5), after substituting (8.6), results in
Aci i2 Bci (8.7)
which is an equation of eigenproblem of matrices A,B. This problem has dimension m ,
which is significantly smaller than the original dimension n.
CT AC D (8.8)
288
then the matrices A and D have identical eigenvalues and they are diagonal elements of the
matrix D. The transformation matrix C is calculated in iterative manner
C S1 S 2 .......S k , k 1.. (8.9)
where the individual S k has the following form
1 0 0 0
1
1
Sk 0 cos( ) 0 sin( ) 0 (8.10)
0 1 0
0 sin( ) cos( ) 0
0 0 0 1
The entries cos( ), sin( ) are put in i,j rows and columns and they are constructed in the
way that they will zeroize aij after the transformation. The other diagonal elements are equal
to 1 and the remaining off-diagonal elements are 0.
In the case of general eigenproblem the whole procedure of constructing S k is very similar.
The matrices S k now adopt the shape
1 0 0 0
1
1
Sk 0 1 0 a 0 (8.11)
0 1 0
0 b 1 0
0 0 0 1
Notice that the matrix S k is not orthogonal anymore. The two variables a,b are calculated to
zeroize off-diagonal elements i,j of the both matrices K and M. Eigenmodes of the problem
are then calculated as
aii'
'
i
2
(8.12)
bii
where aii' , bii' are diagonal elements of transformed (and diagonalized) matrices A, B.
ui ,2 K 1 fi ,1 (8.14)
and the iterating is stop, when ui ,k 1 ui ,k . The above described algorithm tends to converge
to the lowest eigenmodes. If any of these are to be skipped, the initial eigenvector ui ,1 must be
orthogonal to the corresponding eigenvectors. In practice, the vector ui ,k must orthogonalized
with respect to the skipped eigenvectors even during the iterating procedure, as the initial
orthogonality may get (due to some round-off errors) lost.
In the above
m is number of projection eigenmodes, (reasonably higher than the number of required
eigenmodes),
U k is matrix of columnwise arranged m einvectors after k- th iteration,
A k 1 , B k 1 are transformed stiffness and mass matrices of the problem, (having dimension
m<<n),
Ck 1 is matrix of eigenvectors of A k 1 , B k 1 , see (8.9)
290
Δ 2 is matrix with eigenmodes (on its diagonal). Notice that eigenmodes for transformed and
the original eigenmode problem are the same.
The steps 1 thru 4 are repeated until the difference between the two subsequent operations is
negligible.
The solution algorithm (8.16) is in ATENA a bit modified in order to reduce CPU time and
RAM resources and is described below:
Step1- Inverse iteration method:
ˆ MU
U k 1 k
U
KU ˆ
k 1 k 1
ˆ MU
U
k 1 k 1
B U MU T T U
U ˆ
k 1 k 1 k 1 k 1 k 1
The advantage of this procedure over the one defined in (8.16) is that now you don’t need to
store the original and factorised form of the matrix K. Only the factorised form is needed
during the iterations.
A special issue in this method is how to setup the initial vectors U1 . This is what we do in
ATENA. The first vector contains the diagonal elements of M. The next vectors are
constructed in the way that they have zeros everywhere except one entry. This entry
m
correspond to maximum ii and is set to 1.
kii
The procedure as it is , (because of Inverse iteration method), cannot solve for zero
eigenmodes. This may be a problem, especially if we want to analyze structural rigid body
motions or spurious energy modes. If this is the case, shift matrix K by an arbitrary value s ,
i.e. solve the associated eigenproblem
(K s M )us s2 Mus (8.18)
The original eigenvalues and eigenvectors are then calculated by
9.2 References
BATHE, K. J. (1982). Finite Element Procedures in Engineering Analysis. Englewood Cliffs,
New Jersey 07632, Prentice Hall, Inc.
JENDELE, L. (1987). The Orthogonalization of Multiple Eigenvectors in Subspace Iteration
Method. IKM - XI. Internationaler Kongress ueber Anwendungen der Mathematik in
der Ingenieurwissenschaften, Weimar.
WOOD., W. L. (1990). Practical-Time Stepping Schemes. Oxford, Clarenton Press.
292
10 GENERAL FORM OF DIRICHLET BOUNDARY CONDITIONS
A unique feature of ATENA software is the way, in which it implements Dirichlet boundary
conditions. It supports to constraint any degree of freedom (DOF) by a linear of any number
of other structural DOFs. The proposed method of applying and processing the boundary
conditions is computationally efficient and memory economical, because all constraint
degrees of freedoms (DOFs) are eliminated already during assembly of structural global
stiffness matrix and load vectors. The adopted concept has wide range of use and several its
possibilities are discussed. At the end of the Section a few samples are given.
K
j 1
ij u j ri , i 1..n (9.1)
where Kij is an element i, j of a predictor matrix K, (i.e. usually structural stiffness matrix),
ri is an external force, (or unbalanced force) applied into i-th structural degree of freedom
(DOF) and finally ui is displacement (or displacement increment) at the same DOF. Such a
set of equations is always accompanied by many boundary conditions (BCs). They can be one
of the following:
Von-Neumann boundary conditions, (also called right-hand side (RHS) BCs). Number and
type of these BCs has no impact on dimension n of the problem (9.1). They are accumulated
in the vector r . This vector is assembled on the per-node basis for concentrated nodal forces
and/or per-element basis for nodal forces being equivalent to element loads.
The second type of boundary conditions are Dirichlet boundary conditions, (also called left-
hand side (LHS) BCs). ATENA implementation of this type of BCs is now described. A
simple form of such BCs reads:
ul 0, l 1, n
(9.2)
ul ul 0 , l 1, n
These kinds of BCs typically represent structural supports with no displacements, (the first
equation) or with prescribed displacements ul 0 , (the second equation). Although most LHS
BCs are of the above form, (and only a few finite element packages offer anything better),
they are cases, when a more general LHS BC is required. Therefore, ATENA software
provides solution for implementing a form of Dirichlet BCs, where each degree of structural
freedom can be a linear combination of any other degrees of freedom. Mathematically, this is
expressed by:
ul ul 0
k1, n
lk uk , l 1, n (9.3)
ul ul 0 lk uk (9.4)
Substituting (9.4) into the Equation (9.1) yields
n n
j 1, j l
Kij u j Kil ul
j 1, j l
Kij u j Kil (ul 0 lk uk ) ri , i 1..n (9.5)
K
j 1
ij K il lk kj u j ri Kil ul 0 , i 1..n (9.6)
The above set of equations could be already used to solve for the unknown displacements (or
displacement increments) u j . kj stands for k. j Kronecker delta tensor. The trouble is,
however, that even though the matrix K might be symmetric, the set of equations (9.6) is not
symmetric anymore. Thus, to preserve the symmetry, add an lk multiple of the row l , i.e.
n
lk K lj K ll lk kj u j lk rl K ll ul 0 (9.7)
j 1
to the row k, i.e.
n
K
j 1
kj K kl lk kj u j rk K kl ul 0 (9.8)
K
n
kj K kl lk kj lk K lj K ll lk kj u j
j 1
n
K
j 1
kj lk K lj K kl lk lk lk K ll kj u j (9.9)
rk K kl ul 0 lk rl lk K ll ul 0
294
Hence, the final form of the governing set of equations will read:
n
K
j 1
ij K il lk kj ik lk K lj ik kj lk2 K ll u j
(9.10)
ri K il ul 0 ik lk rl K ll ul 0
K
j 1
ij u j ri , i 1..n (9.11)
where
K
r
r1 K1l ul 0
...
ri Kil ul 0
... (9.13)
rk K kl ul 0 lk rl K ll ul 0
...
rj K jl ul 0
...
rn K nl ul 0
is now also
Providing the original matrix K is symmetric, i.e. K ij K ji , then the matrix K
symmetric, i.e. Kij K ji .
The displacement ul constrained by Equation (9.4) has a constant part ul 0 and a variable part
lk u k , in which ul depends only on a single u k . A more general form of this BC would be,
if ul depends on more displacements. It corresponds to the following form of the boundary
condition:
In this case, the displacement ul is calculated as a constant part ul 0 plus a linear combination
lk of displacements u k . k can be any displacement, i.e. k 1..n . Replacing BC defined
by Equation (9.4) by the above Equation (9.14), the equation will change to the form
n
K K ij il lk kj ik lk Klj ik kj lk2 Kll u j
j 1 k ,k l (9.15)
ri Kil ul 0 ik lk rl K ll ul 0
k ,k l
The problem is, however, that displacements u k in (9.16) need not be free but fixed by
another BC, k can run also through l, (resulting in a recursive formulation), more BCs can be
specified for the same ul , a particular BC can be specified more times and in more forms etc.
For example, we may have a set boundary equations that contains BCs
u1 u2 , u2 u1 (9.17)
or it can contain
u1 u2 , u2 u1 , u1 0.003 (9.18)
Both of these are correct. Unfortunately, the set can also contain
u1 u2 , u2 0.003, u2 u1 , u1 0.003 (9.19)
which is definitely wrong. Therefore, before any use of such set of BCs it is necessary to
detect and fix all redundant and contradictory multiple BCs present in it. It is easily done in
case of a simple set of BCs like the one above, but in real analyses with thousands of BCs in
the form (9.16), (some of them quite complex, i.e. k runs through many DOFs) the only way
to proceed is to treat (9.16) as a set of equations to be solved prior their use in (9.13).
Redundant BCs are ignored and contradictory BCs are fulfilled after their summation. Let us
suppose that all structural constraints are specified in the set of equation (9.16). This can be
written in matrix form
ul ul 0 A u k (9.20)
The above relationship represents a system of algebraic linear equations. The system is
typically non-symmetric, sparse and has different number of rows (i.e. number of BCs) and
columns, (i.e. number of master and slave DOFs). Moreover, it is often ill-conditioned, with a
296
number of equations being linear combinations of the others, e.g. see the example in (9.17).
At the beginning it is often not known, which DOF is dependent, (i.e. slave) and which is
independent, (i.e. master), (e.g. see also (9.17)).
Based on the above properties the following procedure has been developed to solve the
problem (9.20):
1. Allocate "columns" for all slave and master DOFs, i.e. loop through all BCs in (9.16)
and allocate DOFs for both slave (i.e. LHS) and master (i.e. RHS) displacements ui .
2. Allocate storage for the matrix A and vectors ul , ul 0 in (9.20). The matrix has lr
number of rows (see (9.16)) and lc number of columns. lc is dimension of the DOFs
map created in the point add. 1.
3. Assemble the matrix A and the vectors ul , ul 0 .
4. Detect constant BCs, i.e. ul ul 0 and swap rows of A so that the rows corresponding
to constant BCs are pushed to the bottom.
5. Detect constant fixed DOFs, i.e. those with lk 0 and variable fixed DOFs, i.e. that
are those dependent on other (master) DOFs and having lk 0 .
6. Swap columns of A , so that the former DOFs are pushed to the right and the latter
DOFs to the left. The operations described at the point 5 and 6 are needed to assure
order, in which the constrained DOFs are eliminated. This is important for later
calculation of the structural reactions.
7. Using Gauss method triangulize the set of BC equations. The triangulization is carried
out in the standard way with the following differences.
a. Before eliminating entries of A located in column below akk , check for a non-
zero entry in the row k. If all its entries are zero, then ignore this row and
proceed to the next one. (It is the case of multiple BCs having the same
content).
b. Check, whether the row k specifies BC for constant or variable DOF, (see
explanation in the point 5 above). In the former case push the row k to the
bottom and proceed to the next row.
c. Swap columns k ...lc so that abs ( akk ) becomes maximum.
d. If akk 0 , swap lines k ...lr to find a non-zero entry in akk . Thereafter,
swap columns k ...lc to find maximum akk .
e. Eliminate entries below akk as usually.
As it was already mentioned, the matrix A is typically very sparse. Hence, a special storage
schemes is used that stores only non-zero entries of A. The data are stored by rows. Each row
has a number of data series, i.e. sequences or chunks of consecutive non-zero data (within the
row). The data are in a three-dimensional container. For each such chunk of data we also
need to store its first position and length. This is done in two two-dimensional containers.
a11 0 0 0 0 0 0
0 a22 a23 0 0 a26 a27
0 0 a33 0 0 0 0
A0 a42 0 a44 0 0 0 (9.21)
0 0 0 0 a55 a56 a57
0 a62 0 0 0 a66 0
0 0 0 0 0 0 a77
It is stored as follows ( A.data stores the actual data, A.rowbase stores base indices for non-
zero entries in rows, A.rowlength contains dimension of non-zero data chunks; all arranged
by rows):
A.data (1)(1)(1) a11
A.data (2)(1)(1) a22 , A.data(2)(1)(2) a23 , A.data (2)(2)(1) a26 , A.data (2)(2)(2) a27
A.data (3)(1)(1) a33 ,
A.data (4)(1)(1) a42 , A.data(4)(2)(1) a44
.......
A.rowbase(1)(1) 1
A.rowbase(2)(1) 2, A.rowbase(2)(2) 6
(9.22)
A.rowbase(3)(1) 3
A.rowbase(4)(1) 2, A.rowbase(4)(2) 4
.....
A.rowlength(1)(1) 1
A.rowlength(2)(1) 2, A.rowlength(2)(2) 2
A.rowlength(3)(1) 1
A.rowlength(4)(1) 1, A.rowlength(4)(2) 1
A number of optimisation techniques are used to speed up the process of triangularization of
the matrix A. These are summarized below:
The data are stored by rows and the elimination is also carried out by rows. (Row-based
storage is also more convenient during assembling the A from (9.16)). All the operations
needed for the elimination are carried out only for nonzero data. Their horizontal position is
stored in A.rowbase and A.rowlength , hence it is no problem to skip all zero entries. A
typical total number of columns lc , see (9.16), is of order from thousands to hundred
thousands DOFs. On the other hand a.rowlength is on average only of order of tens. This is
where the CPU savings comes from.
By the way, the same mapping of non-zero entries is also used for columns. This is achieved
by additional arrays A.columnbase and A.columnlength that are also included in the storage
scheme A. (Their construction is similar to A.rowbase and A.rowlength ; instead by rows
298
they are arranged by columns). These two additional arrays make possible to skip all zero
entries during column-base operations. The resulting significant increase of triangularization
speed pays off for a small amount of an extra RAM that is needed to store A.columnbase and
A.columnlength .
The adopted procedure of triangularization many times swaps lines and/or columns of A. In
view of the adopted storage scheme it can be quite expensive procedure. To alleviate this
problem, the storage scheme includes four additional arrays, namely A.rowindex ,
A.rowinverseindex , A.columnindex and A.rowinverseindex . At the beginning
A.rowindex(i) i and similarly A.rowinverseindex (i ) i, i 1...lc . When a row r1 should be
swapped with a row r2 , the data in A.data remains unchanged and we swap only
corresponding row indices in A.rowindex , (and accordingly also entries in the array for
inverse mapping A.rowinverseindex ). The same strategy is used for swapping the columns.
As a result any swapping operation does not require any moving of actual data, (except of
swapping corresponding indices for mapping the rows and columns) and thus it is extremely
fast. On the other hand, in order to access aij we must use ai ' j ' , where i ' rowindex(i) and
j ' columnindex( j ) . The incurred CPU overhead is well acceptable, because the matrix A is
very sparse.
300
Incompatible meshes on the contact between the blocks using CBCs
Fig. 10-2 (cont) Mesh generation from simple blocks
In the above example two blocks are connected to form a structure, where the top (smaller)
block is placed atop of the bottom (larger) block. Position of the top block is arbitrary with
respect to the bottom block. Unless the concept of CBCs is used, the meshes on interface of
the two blocks must be compatible, (see top of Fig. 10-2). On the other hand, the proposed
CBCs allow using of incompatible meshes, (see the bottom of Fig. 10-2). In this case the
mesh in each block is generated independently, which is significantly simpler. After they are
done, the proposed CBCs are applied to connect the interface nodes. (Typically the surface
with the finer mesh is fixed to the surface with the coarse mesh). The latter approach also
demonstrates possibility of a mesh refinement while still using well-structured and transparent
meshes. This is particularly useful in the case of complex numerical models.
bar 1
bar 2
hi ( r , s ), U i are element interpolation function and Ui are nodal displacements for the
underlying solid element, respectively. For displacement at the node n we can write
4
u (rn , sn ) hi (rn , sn ) U i . ( rn , sn ) are coordinates of the node n. Comparing this formula with
i 1
302
Postprocessing of this element and an ordinary nonlinear hexahedral element is the
same. Consequently, this element does not need any extra support for visualisation of
the results. It makes its implementation and use simple.
Derivation of all ij coefficients and ui 0 constants for all nodes 1 to 12 is beyond the scope
of this document. Nevertheless, a similar procedure is used, as it was in the previous example.
ATENA package covers also Ahmad element for curved shell structures, see Chapter 3.12.
The usual 2D shape of shell element is in the same manner replaced by geometry of a 3D
nonlinear hexahedral element. Originally, the shell element has 3 displacements and 2
rotations at each node located at middle thickness of the shell. These 5 DOFs are in by use of
CBCs replaced by 3 displacements at the top and 2 displacements at the bottom at the
respective nodes from the hexahedron, (i.e. brick) geometry. Advantages of this approach are
the same as those in the case of the curvilinear beam above: simpler pre/post processing,
simpler connection to the adjacent 3D elements, no need to support rotational DOFs during
pre/post processing, no need for extra support for geometry of the shell element.
10.3 References
BATHE, K.J. (1982), Finite Element Procedures In Engineering Analysis, Prentice-Hall, Inc.,
Englewood Cliffs, New Jersey 07632, ISBN 0-13-317305-4.
JENDELE, L, CERVENKA, J, CERVENKA V., " Bond in Finite Element Modelling of
Reinforced Concrete", Proceedings Euro-C 2003, Computatinal Modelling of Concrete
Structures, Swets & Zeitlinger, Lisse, The Netherlands, ISBN 90 809 536 3, 793-8036
crack opening law ........................................ 20 cracking .................................... 17, 25, 35, 41, 52
definition.................15, 83, 191, 217, 245, 293 Crank-Nicholson integration .......................... 255
306
CCIsoQuad element..................................... 90 formulation..................................................... 2
CCIsoTetra element ..................................... 99 general............................................................ 2
CCIsoTriangle element ................................ 96
R
CCIsoWedge element ................................ 104
Rayleigh Ritz method ..................................... 288
problem discretisation.................................... 9
Q10/Q10Sbetaelement............................... 111 S
shell/Ahmad element ................................. 136 Serendipity element ........................ 126, 137, 145
truss 2D/3D element .................................... 85 shape function................................................... 83
introduction ........................................................ 1 SHCC................................................................ 52
Inverse Iteration method................................. 289 S-N curve.................................................... 48, 50
Inverse Subspace Iteration method................. 287 solver .............................................................. 191
isoparametric formulation ................................ 83 Arc length................................................... 202
Arc Length step.......................................... 208
J
Consistently Linearized Method ................ 205
Jacobi method................................................. 288
Crisfield Method ........................................ 207
L direct sparse................................................ 193
Lagrangian element .........................126, 137, 145 Explicit Orthogonal Method ...................... 206
Cholesky .................................................... 192
M
iterative ...................................................... 193
moisture.......................................................... 245
linear .......................................................... 191
multiaxial........................................................ 218
Linear search .............................................. 208
Multipoint constraint ...................................... 293
Modified Newton Raphson ........................ 201
T W
time equivalent ............................................... 246
Wöhler curve .................................................... 50
Transport analysis .......................................... 245
triaxial ........................................................ 54, 56 X
Xi-Bazant model............................................. 256
308