- Email: [email protected]
A fast stress integration algorithm for reinforced concrete sections with axial loads and biaxial bending
Computers and Structures 82 (2004) 213–225 www.elsevier.com/locate/compstruc
A fast stress integration algorithm for reinforced concrete sections with axial loads and biaxial bending J.L. Bonet a, M.L. Romero a
b,*
, P.F. Miguel a, M.A. Fernandez
a
Department of Civil Engineering, Polytechnic University of Valencia, Valencia 46022, Spain Department of Technology, Universitat Jaume I, Campus Riu Sec, Castell on 12071, Spain
b
Received 27 November 2002; accepted 10 October 2003
Abstract This paper deals with the implementation of two integration procedures for arbitrary geometry reinforced concrete cross-sections with axial forces and biaxial bending, from service load to ultimate load. Two methods are proposed: the first method is suitable for a general stress field, and the second method is applicable in sections in which the stress field is uniform at least in one direction. Both methods decompose the integration area into thick layers parallel to the most tensile stressed fiber, whose definition depends on the constitutive equation of the concrete. The integration of the stress field of each thick layer is transformed into a double integral or a path integral over the perimeter of this layer. They are evaluated by a Gauss quadrature. The methods have been tested in different section types by modifying the depth and the inclination angle of the neutral axis. The results obtained for the different alternatives are compared in accuracy and in speed with the results obtained from the classical fiber decomposition method and also with those from an alternative method proposed by a different author. 2003 Elsevier Ltd. All rights reserved. Keywords: Non-linear analysis; Reinforced concrete; Stress integration; Cross-section analysis; Biaxial bending
1. Introduction In general, reinforced concrete structures are subjected to biaxial bending and axial forces due to the geometry, the shape of the cross-section or the nature of external forces. In edification or civil engineering many examples follow this behaviour such as the columns in the corners of buildings and elements affected by seismic and wind forces. With the aim of analysing this type of structures, non-linear analysis programs and, particularly, those based on the finite element method, require the computation of the constitutive equation of the cross-section
*
Corresponding author. Tel.: +34-9647-28134; fax: +349647-28106. E-mail address: [email protected] (M.L. Romero).
and the evaluation of the internal forces through the integration of the stress fields. These programs perform this operation many times; hence its optimisation leads to an important reduction in the computing time consumed. The runtime is mainly consumed during the evaluation of the constitutive stiffness matrix and the internal equilibrium forces of the reinforced concrete element, as demonstrated by Bonet [1] and Romero et al. [2]. According to the type of analysis required, different cross-section integration methods are adopted. If the stress–strain relation is linear, the resolution of the stress integration can be achieved analytically, Br€ ondumNielsen [3] and Rodrıguez-Gutierrez and Dario Aristizabal-Ochoa [4]. Otherwise, the integration procedure is performed numerically by dividing the section into layers for symmetric bending or into cells (fibers) if the cross-section is subjected to biaxial bending moments, Mari [5].
0045-7949/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2003.10.009
214
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
The problem of section estimation in a generalized finite element context has been treated by a lot of authors as Lai et al. [6] who developed an analytical model to simulate the hysteretic and stiffness degrading behavior of reinforced concrete members subjected to axial load and biaxial bending interaction. Kaba and Mahin [7] proposed the first fiber beam element, which uses force shape functions to compute the element flexibility and the displacement shape functions to compute the element-resisting forces. The element had some convergence problems for concrete due to the softening descending branch. Later on, Zeris and Mahin [8,9] improved this element through a finite-element model for the analysis of the non-linear behavior of reinforced concrete columns under uniaxial excitation, but their procedure to compute the element stiffness matrix and the element-resisting forces was not entirely general and required an ad hoc approach to solve softening problems. Bousias et al. [10] proposed a model for the incremental force-deformation behaviour of reinforced concrete sections and members, under generalized load or deformation histories, including cyclic loading, up to ultimate deformation. At the section level the model was of the bounding surface type and accounted for the coupling between the two directions of bending and between them and the axial direction. For the construction of the member tangent flexibility matrix on the basis of the section tangent flexibility matrix, a piecewise-linear variation along the member was assumed for the nine terms of the tangent section flexibility matrix. Spacone et al. [11] presented a fibre beam-column element for the non-linear static and dynamic analysis of reinforced concrete frames. The element formulation was flexibility-based and relied on force interpolation functions that strictly satisfy the equilibrium of bending moments and axial force along the element. They proposed an iterative algorithm for the determination of the resisting forces during the element state determination, which was accurate and stable, even in the presence of strength loss, tracing very well the highly non-linear behaviour of R/C members under cyclic load combinations of bending moment and axial force. The technique based on the decomposition of a dense grid of cells is not numerically efficient, due to the large amount of information needed to characterise the section and the high number of operations required to implement the stress integration with an allowable error level, Miguel et al. [12]. These problems are more important in the structural analysis, during the computation of the global stiffness matrix and the evaluation of the internal forces. Moreover, this method can produce errors in the exact adaptation of the grid of cells to the section stress integration. If the neutral axis is unknown during the numerical procedure, the equilibrium may be achieved by iterating its location. In this case, the convergence
procedure is not assured, because the variation of the internal forces of the cross-section is discrete and depends on the size of the grid. To solve this problem, Rasheed and Dinno [13] proposed an algorithm which uses the secant stiffness approach for the direct prediction of the complete monotonic response of RC rectangular sections under symmetric bending moments. For biaxial bending moments and axial loads many studies have been carried out on the interaction surface of the cross-sections. Worth particular mention is the research by Brondum [14,15] and Yen [16], who presented methods for the analysis of arbitrary cross-sections, but using a rectangular stress block. In addition, Kawakami et al. [17] presented a complete algorithm for the analysis of RC members, including the effect of both regular and prestressed reinforcements, but their method was very difficult to implement in a computer program. Subsequently, Rodrıguez-Gutierrez and Dario Aristizabal-Ochoa [4] presented a method to determine the biaxial interaction diagrams for any orientation of the neutral axis of RC short columns with any geometry, obtaining analytical closed form expressions for equilibrium forces. The same authors [18,19] have latterly presented a method to determine (for arbitrary sections) an analytical model of moment–prestress–curvature (M– P–/) and also for slender columns. Recently, Fafitis [20] developed a method for the computation of the interaction surface of reinforced concrete sections subjected to axial loads and biaxial bending based on Green’s theorem. The integration of the stresses over the compressed region of the concrete section was transformed into line integrals along the compressive perimeter of this region. This method is suitable when the integrated function is uniform over one direction (as is the case where the section is under an increasing monotonic load), where the points located in parallel lines to the neutral axis are under the same stress level. This method is therefore not applicable to sections that do not meet this condition. When a reinforced concrete structure is analysed using numerical methods, the integrating functions in the cross-sections might not be constant along the lines parallel to the neutral axis. This occurs when constitutive equations that take into account the different behaviour of concrete in loadunload processes are implemented. The cyclic loads or the stress distribution processes in the structure may cause some regions of the section experience permanent strains, which are responsible for a non-uniform stress distribution in the lines parallel to the neutral axis. In contrast to the classical cells or layers methods, Miguel et al. [12] developed a numerical algorithm based on the decomposition of the non-cracked part of the section in quadrilaterals. In each of them, the integration of the section was performed by a Gauss–Legendre quadrature.
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
In addition, the same authors presented an alternative integration procedure to improve the accuracy for the same runtime, Bonet et al. [21]. Hence, in the case that the stress function to integrate does not approximate to a low order polynomial or is defined by steps, it would be advisable to subdivide the non-cracked zone by parallel ‘‘thick layers’’ to the neutral axis (Fig. 1). Each one of these layers will configure a polygon that will be decomposed into quadrilaterals by the previously defined algorithm. This integration method will henceforth be referred to as thick layers integration (TLI). It is presented in the next section. Table 1 shows a summary of the integration methods outlined above. This paper presents the TLI method (for non-cylindrical stress fields) and a modified thick layer integration (MTLI) method for the analysis of reinforced concrete
215
sections with any geometric shape (including holes). The algorithms are implemented for combined biaxial bending moments with axial forces, from any load level: service or ultimate stress. The second method is suitable only for sections in which the stress field is uniform over one direction, and it is based on the Gauss–Legendre quadrature. The fundamentals of the methods are explained in the paper and various possibilities for this purpose are analysed. In the proposed formulation of interaction strength estimations plain sections remain plane and perfect bond and lack of shear postulates are implicit. The obtained results for the different alternatives will be compared in accuracy and in speed with the results obtained from the fiber decomposition models and the method proposed by Fafitis [20].
Quadrilateral 1
N
Quadrilateral
Quadrilateral 2 Quadrilateral 3
Type
1
Hollow
2
Solid
3
Solid
εc Reference line.
Strain (a)
Stress
σc
(b)
Fig. 1. Integration with ‘‘thick layers’’ (TLI): (a) decomposition in ‘‘thick layers’’ and (b) decomposition in quadrilaterals.
Table 1 Summary of the integrations methods Classical integration methods of the stress field Method
Procedure
Layers decomposition method
• • • •
Cells decomposition method
The section is decomposed in ‘‘nb’’ layers The stress and Et is evaluated in the center of each layer The section is decomposed in a mesh of ‘‘n m’’ cells/fibers The stress and the tangent elasticity modulus (Et) are evaluated in the centre of each cell
Methods based on the Gauss–Lengendre quadrature Fixed integration points
The Gauss points are ‘‘fixed’’ and independent of the neutral axis position
Variable integration points Non-cylindrical stress fields Without layers decomposition (Miguel et al.) With layers decomposition Cylindrical stress fields Without layers decomposition (Fafitis)
With layers decomposition
Steps: • Decomposition in ‘‘Nc ’’ quadrilaterals • Double integration of each quadrilateral Thick layer integration method (TLI) Steps: • Decomposition in ‘‘nL’’ quadrilaterals • Double integral fi Path integral Modified thick layer integration (MTLI)
216
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
2.1. Decomposition of the integration area
The algorithm has limitations in some specific cases, where the problem is non-monotonic even under monotonic bending increase and constant axial load, due to e.g.:
• Interaction surface estimation. • Cross-sectional analysis. • Moment–curvature diagrams computation.
The integration area is decomposed into quadrilaterals to integrate stress fields in each region by a trivial quadrature method. Once the polygonal contour of the integration area has been determined, a reference line is chosen, Fig. 2a, and the orthogonal projections points from the vertexes of the contour are obtained (Fig. 2b). The section is then decomposed into quadrilaterals and triangles from the two consecutive vertexes and its projections, Fig. 2c. In this way, quadrilaterals and triangles with ‘‘positive’’ or ‘‘negative’’ areas are obtained, Fig. 2d.
At present, the authors are working in the inclusion of a new monitored variable to avoid previous limitations.
2.2. Internal force and constitutive matrix of the quadrilaterals
• Concrete degradation of reinforcement buckling. • Stress re-distribution in structures. But the proposed method performs adequately in:
The stress integration in each quadrilateral is performed using a Gauss quadrature. The quadrilaterals can be transformed into a square of 2 · 2 units by side, Zienkiewicz [22], using a coordinate transformation. The internal forces and tangent constitutive matrix of the quadrilaterals are obtained trivially. The sign of the Jacobian states whether the quadrilateral is solid or hollow. To this end the vertex is enumerated in a clockwise direction (positive) or anti-clockwise (negative). The proposed algorithm implements this enumeration automatically. The stresses are obtained by the
2. Thick layers integration method The contour where the stresses must be integrated can be considered as an n-sides polygon. If the section is made up of curved contours, it will also be reduced by polygons. The integration procedure of the section can be divided in three phases or steps: decomposition of the integration area in quadrilaterals or triangles, stress integration of each quadrilateral and layer subdivision.
Quadrilateral- 1 Integration Area
Quadrilateral- 2 2
3
2
4
4 - (4’)
Reference Line
Reference Line
(3’)
(2’)
1 — (1’)
1
Quadrilateral- 3
3
(a)
(c)
(b) QUADRLATERAL NUMBER
VERTEX
TYPE
2 SOLID DEGENERATED QUADRILATERAL
1 3
1
4 2
1
3
2 4 3
1 2
4 3
SOLID NON DEGENERATED QUADRILATERAL EMPTY DEGENERATED QUADRILATERAL
(d) Fig. 2. (a) Integration area and reference line, (b) vertex projection, (c) quadrilaterals determination and (d) vertex enumeration.
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
addition of each stress that actuates over the quadrilaterals. In addition, the constitutive matrix of the sections is obtained by composition of the quadrilaterals. 2.3. Thick layers decomposition Certain considerations must be taken into account in the proposed method. First, the utilisation of non-linear structural analysis programs for reinforced concrete structures establishes the use of iterative algorithms. Hence, during the convergence process its is observed that the neutral axis and the stress field are different for all the cross-sections into which the beams of the structure are divided. Furthermore, the cracking phenomenon in the concrete characterises the behaviour of this type of structure. In the cracked region of the cross-section, the concrete fibers do not add stiffness to the structure. Moreover, the integration area may be made up of all the cross-section or only of the regions with a non-zero stress field. In the first case, the location of the Gauss points will be fixed during the whole convergence process. Hence, the decomposition of the integration area will be performed once during the procedure. Otherwise, if the stress field influences the integration area, it will be redefined in each iteration. Hence, the allocation of the integration points will change during the process. The accuracy of the integrals depends for each quadrilateral on the number of Gauss points used and the polynomial form of the stress–strain relation. When the stress field or the elasticity modulus functions are not adjusted to a small order polynomial, a substantial improvement in the integration definition is brought about by decomposing the area into non-zero stress regions, by ‘‘thick layers’’ parallel to the neutral axis (Fig. 1b). In this situation, each layer defines an n-sides polygon where the proposed method can be applied.
217
The smaller deformation point of the layer may define the reference line. In consequence, the internal forces or the constitutive matrix are obtained by summing the terms regarded with the quadrilaterals of each layer into which the integration area is divided.
3. Modified thick layers integration method The second method is suitable only for sections in which the stress field is uniform over one direction, what can be called ‘‘cylindrical stress field’’, and in this way the stresses are integrated for arbitrary geometry RC cross-section with biaxial moments and axial forces. Again, the contour where the stresses must be integrated can be considered as an n-sides polygon. The proposed procedure can be divided into three phases: decomposition into ‘‘thick layers’’ of the specific area of integration (non-zero stress field of the concrete fibers), transformation of the stress integral for each thick layer into a ‘‘path’’ integral over its perimeter, and finally the evaluation of the ‘‘line’’ integral of stresses with a Gauss quadrature. If the integration function is defined by a curve that approximates adequately in its total range to a polynomial function, the subdivision of the integration area into thick layers is unnecessary. 3.1. Transformation of the stress integral into a ‘‘path integral’’ over the perimeter of the concrete integration area If the stress field is uniform in a particular direction (rc ðy; zÞ ¼ rc ðzÞ; Ect ðy; zÞ ¼ Ect ðzÞ), the stress integral of the integration area can be transformed into a path integral over the perimeter of the integration area (Fig. 3). The internal forces can be obtained by the expressions:
z
z
2 Integration Area
L1
L2
y dz (y/2-y0)
3
(z-z0)
dz
1 dy y
(y0, z0)
L3
Reference Line
y L5
4
L4
Fig. 3. Transformation of the stress integral.
5
Reference Line
218
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
Nc ¼ ¼
Z Z Z ZAc
rc ðy; zÞ dy dz
3.2. Evaluation of the stress integral by the Gauss quadrature
rc ðzÞ dy dz
As the integration area contour is approximated by a polygon, the stress field integral over the perimeter L, can be obtained by decomposing this integral side by side through the perimeter (Fig. 4).
Ac
I
y rc ðzÞ dz ZL Z ðz z0 Þ rc ðzÞ dy dz Mcy ¼ Ac I ¼ ðz z0 Þ y rc ðzÞ dz ZL Z ðy y0 Þ rc ðzÞ dy dz Mcz ¼ Ac I ¼ ðy=2 y0 Þ y rc ðzÞ dz ¼
ð1Þ
Z Z
¼ f ðy; zÞ dy dz ¼ Ac
¼
Z ZAc I
SLi ðzÞ
ð3Þ
i¼1
SLi ðzÞ ¼
Z
ziþ1
hðzÞ dz zi
In order to perform the stress integral of a determined side of the contour of the integration area (Li ), the wellknown Gauss–Legendre method is used. The next coordinate transformation must be implemented:
Ect ðy; zÞ dy dz Ect ðzÞ dy dz
Ac
y Ect ðzÞ dz ZL Z ðz z0 Þ Ect ðzÞ dy dz Dct ð1; 2Þ ¼ I Ac ¼ ðz z0 Þ y Ect ðzÞ dz ZL Z ðy y0 Þ Ect ðzÞ dy dz Dct ð1; 3Þ ¼ I Ac ¼ ðy=2 y0 Þ y Ect ðzÞ dz ZL Z ðz z0 Þ2 Ect ðzÞ dy dz Dct ð2; 2Þ ¼ Ac I ¼ ðz z0 Þ y Ect ðzÞ dz ZL Z ðz z0 Þ ðy y0 Þ Ect ðzÞ dy dz Dct ð2; 3Þ ¼ Ac I ¼ ðz z0 Þ ðy=2 y0 Þ Ect ðzÞ dz ZL Z ðy y0 Þ2 Ect ðzÞ dy dz Dct ð3; 3Þ ¼ Ac I ¼ ½y 2 =12 þ ðy=2 y0 Þ2 y Ect ðzÞ dz ¼
nL X
where nL is the number of sides that forms the integration area:
The elements of the tangent stiffness matrix of the section (Dct ) can also be obtained: Z Z
hðzÞ dz ¼
L
L
Dct ð1; 1Þ ¼
I
z¼
ðziþ1 zi Þ ðziþ1 zi Þ ðn þ 1Þ þ zi ) dz ¼ dn 2 2
ð4Þ
The stress field integral of each side of the contour is
SLi ðzÞ ¼
Z
ziþ1
hðzÞ dz ¼ zi
ð2Þ
L
where Ac is the integration area, L is the contour of this area, rc is the concrete stress at each fiber (y; z), Ect is the tangent elasticity modulus of deformation of the fiber (y; z); and (y0 ; z0 ) are the coordinates of the centre of reference of forces of the section. The ‘‘y’’ axis is situated parallel to the direction in which the stresses are uniform.
¼
ðziþ1 zi Þ 2
Z
þ1
hðnÞ dn
1
n ðziþ1 zi Þ X ki hðni Þ 2 i¼1
ð5Þ
where ni is the value of the curvilinear coordinate of the Gauss point, ki is the weight associated to the Gauss point and ‘‘n’’ is the number of Gauss points.
z (yi+1,zi+1) y
Li dz
(y/2-y0) (z-z ) 0
(y0, z0)
y = yi +
(yi -yi) (z -z ) i (zi+1-zi)
(yi-zi)
y Reference Line
Fig. 4. Evaluation of the stress integral.
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
40 5
ð6Þ The proposed algorithm allows us to integrate the stresses in a hollow section or in a solid section. Hence, the vertexes of the quadrilateral must be enumerated again in counter-clockwise (positive sign) or clockwise (negative sign).
20
Lj
5
L
" # nL n X ðziþ1 zi Þ X ¼ ki hðni Þ 2 j¼1 i¼1
90
Ac
20
In this way, the integral is I Z Z f ðy; zÞ dy dz ¼ hðzÞ dz
219
3.3. Stresses integration through ‘‘thick layers’’ The accuracy of the integrals depends for each one of the sides on the number of Gauss points used and the polynomial form of the stress–strain relation. When the function to integrate does not follow a low order polynomial, or is defined as a steps function, a large number of Gauss points has to be used to obtain an acceptable accuracy. In these cases it is better to subdivide the integration area into thick layers parallel to the neutral axis (Fig. 1a). Each one of the layers define an n-sides polygon where the proposed integration method can be performed. Subsequently, the application of this method to the classical stress–strain relationship of the concrete (Fig. 1a) may be performed using five layers. Hence, the internal forces and the constitutive matrix of the section can be obtained as the sum of the terms of each side that composes each perimeter of each ‘‘thick layer’’ into which the integration area has been decomposed. This method will be hereafter referred to as modified thick layers integration (MTLI).
4. Numerical examples The methods described above were tested in three different sections, in which the strength of concrete, and the depth and the inclination angle of the neutral axis were modified with the purpose of integrating the stress field. Four methods were implemented: the classical fiber decomposition approach, the integration method proposed by Fafitis [20], the TLI and the MTLI proposed in this paper. The comparison between the different methods was performed in terms of the time and the accuracy obtained. The sections tested were: 20 · 20 cm rectangular; hollow 20 · 20 cm rectangular with 5 cm thick walls and a section in ‘‘T ’’ (Fig. 5). The ‘‘exact’’ solution is considered as the solution solved with the MTLI method with 48 Gauss points for
15
5
20
5
15
60 Fig. 5. Section in ‘‘T ’’ (cm).
each side (Li). Two strengths of concrete were considered for all the cases, 25 and 80 MPa. The constitutive equation used corresponds to the Model Code 1990 [23], and also includes the tensile behaviour and the tension stiffening effect through a model of gradual descending unload, CEB-1995, [24]. In this case the concrete has a cylinder strength distribution. The tensile strain descending branch shape depends on the stress of the most tensile stressed reinforced bar of the section. In the numerical implementation of this study the steel contribution has been ignored because it does not influence the accuracy and the runtime of the integration process. Due to this, the equation adopted for this branch will be a straight line from the maximum tensile stress (fct ) and the ultimate strain deformation for the concrete, which corresponds to the plastic steel deformation (ecend ¼ fy =Es ), Fig. 6. For each case, the stress distribution was integrated assuming that the ultimate compressed concrete fiber had a deformation of 0.0035 and the most tensile stressed fiber of )0.0025 (Fig. 7).
σc
Eci 1
fc
Ec1
ε c,end
1
εc εc1
εc,lim
-fct Fig. 6. Concrete stress–strain relationship adopted.
220
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
For each angle, three depths of the most tensile stressed fiber were computed, corresponding to 1/3, 2/3 and 3/3 of the maximum height (hmax ) orthogonal to the neutral axis (Fig. 7).
hmáx 3.5
Fiber under tension
εc
-2.5
Strain
σc
Stress
Fig. 7. Stress distribution.
To perform the integration, the ‘‘TLI’’ method was used with a mesh of points of 2 · 2, 3 · 3, 4 · 4, 6 · 6, 8 · 8, 10 · 10, 12 · 12, 15 · 15, 20 · 20 and 32 · 32. In each one of the quadrilaterals where the integration area was divided, and for the Fafifis method and the ‘‘MTLI’’ method, the following set of points was used: 2, 3, 4, 6, 8, 10, 12, 15, 20 and 32, for each side of the perimeter of the total integration area or the perimeter of each thick layer that forms the integration area, respectively. For the fiber integration approach, a variable mesh was used from 20 · 20 to 1200 · 1800 cells, depending on the type of section, in order to obtain a better adaptation of the mesh to the geometry of the section (Table 2). Different inclination angles from the neutral axis were tested, depending on the type of section (Table 3). Table 2 Cells mesh Section type
Cell mesh
Rectangular
20 · 20, 40 · 40, 100 · 100, 200 · 200, 400 · 400 20 · 20, 40 · 40, 100 · 100, 200 · 200, 400 · 400 60 · 90, 120 · 180, 300 · 450, 600 · 900, 1200 · 1800
Hollow rectangular Section in ‘‘T ’’
Table 3 Position of the last fiber under tension Section type
Inclination of the neutral axis
Depth
Rectangular
0, 30, 45
Hollow rectangular
0, 30, 45
Section in ‘‘T ’’
0
hmax , 1/3 hmax , 2/3 hmax hmax , 1/3 hmax , 2/3 hmax hmax , 1/3 hmax , 2/3 hmax
5. Results In this section the previously computed numerical examples corresponding to the three section types are analysed. In each example, the error and the runtime are calculated. It was observed that the computing integration time using quadrature was extremely small and lower than computer clock accuracy. Hence, the runtime was measured for the computation of 10 test series. Each series repeats the analysis of the section 50 times. This procedure was repeated for the four algorithms. The error obtained in each method in relative percentage to the exact value is Fi Fexac 100 Errorð%Þ ¼ Fexac
ð7Þ
where Fi is the exact result obtained in case ‘‘i’’ for force F (axial or bending) and Fexac is the exact value for force F (axial or bending). Table 4 shows, for all the types of sections and number of Gauss points tested, the runtime and the maximum errors obtained for all forces, corresponding to the method proposed by Fafitis. It should be noted that the errors for a small number of Gauss points (2, 3, 4 and 6) are not allowable, principally for the 80 MPa concrete. This result appears because the numerical integration using the Gauss– Legendre quadrature of any function f ðxÞ is based on the interpolation of the function with an n-order polynomial, in such a way that the Gauss–Legendre quadrature of n-order, exactly integrates a polynomial of order 2n 1 or lower. Because the stress field is defined by a step function and there is no continuity in the derivative, the polynomial interpolation can produce important integration errors. On the other hand, it can be observed that the errors increase when the strength of the concrete increases. This behaviour comes from the form of the constitutive equation of the concrete where the deformation of the maximum stress is independent of the concrete strength (Fig. 8). However, since the stress function takes on a narrower shape as the strength increases, a higher number of Gauss points are required for sections with 80 MPa concrete than for sections with 25 MPa concrete to obtain the same degree of accuracy.
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
221
Table 4 Results of the Fafitis method Section type
fc (MPa)
Number of Gauss points 2
Rectangular
25 80 25 80 25 80
Hollow rectangular Section in ‘‘T ’’
3
4
6
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
1.32 1.32 2.57 2.64 4.17 4.29
12.47 65.88 12.47 65.98 6.43 81.09
1.65 1.65 3.13 3.07 5.21 5.21
10.50 19.3 10.50 19.3 3.70 32.95
1.92 1.92 3.79 3.85 6.32 6.26
3.57 13.39 3.57 13.39 1.63 29.03
2.63 2.59 4.84 4.99 7.96 8.07
0.22 7.09 0.22 7.09 0.16 21.66
Fig. 8. Uniaxial behaviour of the concrete in compression.
If the integration area is decomposed into thick layers, and each thick layer is also subdivided into quadrilaterals (TLI method), for a small number of Gauss points (2 · 2, 3 · 3, 4 · 4 and 6 · 6) the error is reduced considerably (Table 5). In this way, it can once again be observed that the errors of the TLI method rise when the concrete strength increases. However, the computation time increases notably compared with that of the Fafitis method. If the stress field of the section is uniform in a fixed direction, the ‘‘TLI’’ method does not optimise the
location of the Gauss points, because the integral in this direction can be obtained analytically (Eq. (6)). For this reason, the runtime increases in comparison with the Fafitis method. Through the proposed ‘‘MTLI’’ method, the location of the Gauss points is optimised, and hence, the runtime decreases considerably (Table 6), obtaining an accuracy level similar to the ‘‘TLI’’ method. Table 7 compares the degree of accuracy and the computation time of the three Gaussian methods versus the number of Gauss points for the 80 MPa concrete strength. It can be observed that although, in general, a higher number of Gauss points allows us to obtain a higher accuracy, the results obtained through the Fafitis method do not follow this behaviour due to the irregularity that exists in the stress field of the section. On the other hand, it can be observed that the errors in the integration methods ‘‘TLI’’ and ‘‘MTLI’’ are similar. However, the runtime of the proposed method is notably lower than that of the ‘‘TLI’’ method, and this behaviour is relatively important as the number of Gauss points increases in the integration. Table 8 shows the errors and the computation time for the rectangular section (both hollow and solid), corresponding to the fiber decomposition method. In Table 9 the results of the section in ‘‘T ’’ are presented. It can be observed that although the relative errors obtained are small, the computation time is higher
Table 5 Results of the ‘‘TLI’’ method Section type
Rectangular Hollow rectangular Section in ‘‘T ’’
fc (MPa)
25 80 25 80 25 80
Number of Gauss points 2·2
3·3
4·4
6·6
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
11.48 12.25 20.43 23.06 12.57 13.61
0.10 0.41 0.13 0.45 0.12 0.80
21.91 22.84 37.40 42.91 19.89 21.42
0.0014 0.045 0.00172 0.0485 0.0016 0.39
34.77 37.02 60.37 69.43 29.98 32.78
2.16e)5 0.0022 2.16e)5 0.0018 2.13e)5 0.145
72.02 76.73 124.62 142.52 58.06 63.88
2.89e)9 0.0052 3.6e)9 0.00069 3.5e)9 0.016
222
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
Table 6 Results of the proposed ‘‘MTLI’’ method Section type
fc (MPa)
Number of Gauss points 2
Rectangular Hollow rectangular Section in ‘‘T ’’
25 80 25 80 25 80
3
4
6
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
4.39 6.21 5.66 9.12 12.19 13.34
0.54 0.25 0.54 0.34 0.12 0.80
5.94 8.02 7.25 11.76 15.32 16.75
0.0068 0.1223 0.0068 0.244 0.0016 0.39
6.97 10 8.95 14.28 18.62 20.38
9.58e)5 0.0448 9.58e)5 0.10 2.43e)5 0.145
9.45 13.89 11.92 19.27 24.71 27.24
1.4e)6 0.00494 1.4e)6 0.012 3.5e)9 0.016
Table 7 Results (rectangular section; neutral axis angle: 45; depth: hmax ; fc ¼ 80 MPa) Gauss points
Integration method Fafitis
2 3 4 6 8 10 12 15 20 32
‘‘TLI’’
‘‘MTLI’’
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
2.03 2.69 3.41 4.67 5.77 7.2 8.51 10.39 13.56 21.26
17.08 6.88 0.73 3.54 0.42 1.56 0.52 0.77 0.25 0.07
13.84 26.03 42.4 87.72 150.22 229.32 326.14 502.29 880.51 2219.1
0.38 0.04 0.0015 0.00053 6.12e)05 5.74e)06 0.52e)6 5.8e)09 0 0
5.38 7.04 8.63 11.92 15.16 18.35 21.86 26.47 34.55 54.05
0.148 0.158 0.068 0.0082 0.00088 9.06e)05 9.23e)06 2.88e)07 1e)09 0
Table 8 Results of the fiber decomposition method Section type
fc (MPa)
Cell mesh 10 · 10
20 · 20
100 · 100
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
Rectangular
25 80
22 22
2.19 4.63
82.5 82
0.43 2.30
519 513.5
0.035 0.51
Hollow rectangular
25 80
19 16.5
2.19 4.64
66 69
0.43 2.30
406 412
0.035 0.51
Table 9 Results of the fiber decomposition method Section type
fc (MPa)
Cell mesh 60 · 90
Section in ‘‘T ’’
25 80
120 · 180
360 · 540
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
192 195
0.056 0.52
758 774.5
0.016 0.023
4737.5 4847
0.0037 0.008
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
223
25 MPa 100 Error (%)
10
Cells
1
Fafitis
0.1
TLI
0.01
MTLI
0.001 0.0001 1E-05 1E-06 1E-07 1E-08 1E-09 1E-10 1
10
100
1000
10000
100000
Time (sg) 80 MPa Error (%)
100 10
Cells
1
Fafitis TLI
0.1 0.01
MTLI
0.001 0.0001 1E-05 1E-06 1E-07 1E-08 1E-09 1E-10 1
10
100
1000
10000
100000
Time(sg)
Fig. 9. Efficiency of the fours methods for 25 and 80 MPa for a rectangular section with a ¼ 45 and depth ¼ hmax .
than that obtained with the ‘‘TLI’’ and the ‘‘MTLI’’ method. In general, as the mesh density increases, the error decreases (Tables 8 and 9). However, the accuracy of this method depends on how the cell mesh adapts to the integration area, in that a mesh with a higher number of cells, but distributed in worse conditions respective to the geometry of this area will produce higher errors. In order to compare the four methods, Figs. 9 and 10 are presented. These figures present the percentage error on the Y -axis, and the integration time on the X -axis. The objective of a good algorithm is to achieve lower time consumption with fewer errors. Fig. 9 shows one numerical example in particular: rectangular section with a ¼ 45 and the depth of the neutral axis to the maximum. When the three Gaussian methods are compared with the cell method, it can be noted that the ‘‘MTLI’’ method is the most efficient. This is followed by the ‘‘TLI’’ method, and then the Fafitis method. The Fafitis method is faster than the MTLI but presents inadmis-
sable errors. MTLI is slower because it spends time decomposing the integration area into layers. For a nonlinear numerical procedure, the error must be fixed to arrive at the convergence, which is what reduces the efficiency of the Fafitis method. For example, for the case of Fig. 9 with 25 MPa, and for an admissible error of 0.001%, the ‘‘MTLI’’ method is 2.7 times faster than the Fafitis method, and 1.3 times faster than the ‘‘TLI’’ method. The difference between the TLI and MTLI methods and the Fafitis method is that the former decompose the stress function by trams and in consequence the interpolation polynomial adapts better to the stress function, optimising the properties of the Gaussian integration. It can also be noted that the time integration for the numerical examples of the 80 MPa examples are slower than for the 25 MPa examples. The reason is that the stress–strain relationship of the 80 MPa presents a narrower behaviour and it needs a higher order polynomial integration and in consequence more integration points and more time for the same accuracy.
224
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225 25 MPa 100
Error (%)
10 1 0.1 0.01 0.001
CELLS
0.0001 1E-05
FAFITIS
1E-06 1E-07
MTLI
TLI
1E-08 1E-09 1E-10 1
10
100
1000
10000
100000
Time (sg) Fafitis
Cells
MTLI
TLI
80 MPa 100 10 1 0.1
Error (%)
0.01
CELLS
0.001 0.0001 1E-05 1E-06
TLI
MTLI
FAFITIS
1E-07 1E-08 1E-09 1E-10 1
10
100
1000
10000
100000
Time (sg) Fafitis
Cells
MTLI
TLI
Fig. 10. Efficiency of the fours methods for 25 and 80 MPa for all the examples computed.
6. Conclusions This paper deals with the implementation of two integration procedures for arbitrary geometry reinforced concrete cross-sections with axial forces and biaxial bending, from service load to ultimate load. Two methods are proposed: the TLI method is suitable for a general stress field, and the MTLI method is applicable in sections in which the stress field is uniform at least in one direction. Both methods decompose the integration area into thick layers parallel to the most tensile stressed fiber, whose definition depends on the constitutive equation of the concrete. The integration of the stress field of each thick layer is transformed into a double (TLI) or a path integral (MTLI) over the perimeter of this layer, and they are evaluated by a Gauss quadrature. In the study of the four methods analysed (fibers, Fafitis [20], ‘‘TLI’’ and ‘‘MTLI’’), applied to different concrete section types, the following conclusions are reached:
(a) In all the cases, for an equivalent computation time, the accuracy of the Gaussian integration methods is notably higher than for the fiber method. (b) The accuracy of the Gaussian integration methods increases in relation to the number of Gauss points used, except in the method proposed by Fafitis [20]. The accuracy level of this method depends on the form of the function used to represent the concrete stress–strain relationship. (c) The accuracy level of the fiber decomposition method does not depend only on the mesh density but also on its adaptation to the geometry of the section. (d) Due to the form adopted by the concrete constitutive equation defined by the Model-Code 90 (which depends on the strength), the accuracy level is reduced as the strength of concrete is increased for all the integration methods. (e) Both the thick layer integration (TLI) and modified thick layer integration (MTLI) method show an excellent accuracy level (lower than 0.02%) for a 6 · 6 Gauss mesh points per quadrilateral or 6 Gauss
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
points per side of the perimeter that forms the thick layer. (f) When the accuracy level and the runtime for the three Gaussian methods are compared with the fiber decomposition method, it is observed that the ‘‘MTLI’’ method is more efficient, followed by the ‘‘TLI’’ method and finally the Fafitis method. Finally, the proposed ‘‘modified thick layers’’ method (MTLI), with regard to the accuracy, efficiency, and continuity in the stress field integration is recommended for implementation in non-linear reinforced concrete frameworks programs.
Acknowledgements This paper was supported by the Spanish ‘‘Ministerio de Ciencia y Tecnologia’’ under project MAT200202461.
References [1] Bonet JL. Metodo simplificado de calculo de soportes esbeltos de hormig on armado de secci on rectangular sometidos a compresi on y flexi on biaxial, PhD thesis, Civil Engineering Department, Technical University of Valencia, July 2001 [in Spanish]. [2] Romero ML, Miguel PF, Cano JJ. A parallel procedure for nonlinear analysis of reinforced concrete three dimensional frames. Comput Struct 2002;80(16–17):1337–50. [3] Br€ ondum-Nielsen T. Serviceability analysis of concrete sections under biaxial bending. ASCE, J Struct Eng 1997;123(1):117–9. [4] Rodrıguez-Gutierrez JA, Dario Aristizabal-Ochoa J. Biaxial interaction diagrams for short RC columns of any cross section. J Struct Eng, ASCE 1999;125(6):672–83. [5] Mari AR. Nonlinear geometric, material and time dependent analysis of three dimensional reinforced and prestressed concrete frames. Report No. USB/SESM-84/12, Department of Civil Engineering, University of California, Berkley, California, USA, June 1984. [6] Lai SS, Will GT, Otani S. Model for inelastic biaxial bending of concrete members. J Struct Eng-ASCE 1984;110(11):2563–84. [7] Kaba S, Mahin SA. Refined modeling of reinforced concrete columns for seismic analysis. Report EERC 8403, Earthquake Engineering Research Center, University of California, Berkeley, California, 1984.
225
[8] Zeris CA, Mahin SA. Analysis of reinforced-concrete beam-columns under uniaxial excitation. J Struct EngASCE 1988;114(4):804–20. [9] Zeris CA, Mahin SA. Behavior of reinforced-concrete structures subjected to uniaxial excitation. J Struct EngASCE 1991;117(9):2640–56. [10] Bousias SN, Panagiotakos TB, Fardis MN. Modelling of RC members under cyclic biaxial flexure and axial force. J Earthquake Eng 2002;6(2):213–38. [11] Spacone E, Filippou FC, Taucer FF. Fibre beam-column model for non-linear analysis of R/C frames. 1. Formulation. Earthquake Eng Struct 1996;25(7):711–25. [12] Miguel PF, Bonet JL, Fernandez MA. Integraci on de tensiones en secciones de hormig on sometidas a flexocompresi on esviada. Rev Int Metodos Numer para el calculo dise~ no ingenierıa 2000;16(2):209–25 [in Spanish]. [13] Rasheed HAS, Dinno KS. An efficient nonlinear analysis of RC sections. Comput Struct 1994;53(3):613–23. [14] Brondum T. Ultimate limit states of cracked arbitrary cross sections under axial loads and biaxial bending. ACI Concr Int 1982;4(11):51–5. [15] Brondum T. Ultimate flexural capacity of cracked polygonal concrete sections with circular holes under biaxial bending. ACI Concr Int 1987;84(3):212–5. [16] Yen JR. Quasi-Newton method for reinforced concrete column analysis and design. J Struct Eng, ASCE 1991; 117(3):657–66. [17] Kawakami MT, Kagaya M, Hirata M. Limit states of cracking and ultimates strength of arbitrary cross sections under biaxial loading. ACI J 1985;82(1):203–12. [18] Rodriguez JA, Aristazabal-Ochoa JD. M–P–/ diagrams for reinforced, partially, and fully prestressed concrete sections under biaxial bending and axial load. J Struct Eng, ASCE 2001;127(7):763–73. [19] Rodriguez JA, Aristazabal-Ochoa JD. Reinforced, partially, and fully prestressed slender concrete columns under biaxial bending and axial load. J Struct Eng, ASCE 2001;127(7):774–83. [20] Fafitis A. Interaction surfaces of reinforced-concrete sections in biaxial bending. J Struct Eng, ASCE 2001; 127(7):840–6. [21] Bonet JL, Miguel PF, Fernandez MA, Romero ML. Efficient procedure for stress integration in concrete sections using a Gauss–Legendre quadrature. In: Topping BHV, editor. Proceedings of the Eighth International Conference on Civil and Structural Engineering Computing. Stirling, United Kingdom: Civil-Comp Press; 2001. p. 53. [22] Zienkiewicz OC, Taylor RL. In: The finite element method, vol. 1. CIMNE, Mc-Graw-Hill; 1995. [23] Comite Euro-internacional du beton: CEB-FIB Model Code 1990. CEB Bull 1991;203–5. [24] Comite Euro-internacional du beton: new developments in non-linear analysis methods. CEB Bull 1995;229:94.
A fast stress integration algorithm for reinforced concrete sections with axial loads and biaxial bending J.L. Bonet a, M.L. Romero a
b,*
, P.F. Miguel a, M.A. Fernandez
a
Department of Civil Engineering, Polytechnic University of Valencia, Valencia 46022, Spain Department of Technology, Universitat Jaume I, Campus Riu Sec, Castell on 12071, Spain
b
Received 27 November 2002; accepted 10 October 2003
Abstract This paper deals with the implementation of two integration procedures for arbitrary geometry reinforced concrete cross-sections with axial forces and biaxial bending, from service load to ultimate load. Two methods are proposed: the first method is suitable for a general stress field, and the second method is applicable in sections in which the stress field is uniform at least in one direction. Both methods decompose the integration area into thick layers parallel to the most tensile stressed fiber, whose definition depends on the constitutive equation of the concrete. The integration of the stress field of each thick layer is transformed into a double integral or a path integral over the perimeter of this layer. They are evaluated by a Gauss quadrature. The methods have been tested in different section types by modifying the depth and the inclination angle of the neutral axis. The results obtained for the different alternatives are compared in accuracy and in speed with the results obtained from the classical fiber decomposition method and also with those from an alternative method proposed by a different author. 2003 Elsevier Ltd. All rights reserved. Keywords: Non-linear analysis; Reinforced concrete; Stress integration; Cross-section analysis; Biaxial bending
1. Introduction In general, reinforced concrete structures are subjected to biaxial bending and axial forces due to the geometry, the shape of the cross-section or the nature of external forces. In edification or civil engineering many examples follow this behaviour such as the columns in the corners of buildings and elements affected by seismic and wind forces. With the aim of analysing this type of structures, non-linear analysis programs and, particularly, those based on the finite element method, require the computation of the constitutive equation of the cross-section
*
Corresponding author. Tel.: +34-9647-28134; fax: +349647-28106. E-mail address: [email protected] (M.L. Romero).
and the evaluation of the internal forces through the integration of the stress fields. These programs perform this operation many times; hence its optimisation leads to an important reduction in the computing time consumed. The runtime is mainly consumed during the evaluation of the constitutive stiffness matrix and the internal equilibrium forces of the reinforced concrete element, as demonstrated by Bonet [1] and Romero et al. [2]. According to the type of analysis required, different cross-section integration methods are adopted. If the stress–strain relation is linear, the resolution of the stress integration can be achieved analytically, Br€ ondumNielsen [3] and Rodrıguez-Gutierrez and Dario Aristizabal-Ochoa [4]. Otherwise, the integration procedure is performed numerically by dividing the section into layers for symmetric bending or into cells (fibers) if the cross-section is subjected to biaxial bending moments, Mari [5].
0045-7949/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2003.10.009
214
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
The problem of section estimation in a generalized finite element context has been treated by a lot of authors as Lai et al. [6] who developed an analytical model to simulate the hysteretic and stiffness degrading behavior of reinforced concrete members subjected to axial load and biaxial bending interaction. Kaba and Mahin [7] proposed the first fiber beam element, which uses force shape functions to compute the element flexibility and the displacement shape functions to compute the element-resisting forces. The element had some convergence problems for concrete due to the softening descending branch. Later on, Zeris and Mahin [8,9] improved this element through a finite-element model for the analysis of the non-linear behavior of reinforced concrete columns under uniaxial excitation, but their procedure to compute the element stiffness matrix and the element-resisting forces was not entirely general and required an ad hoc approach to solve softening problems. Bousias et al. [10] proposed a model for the incremental force-deformation behaviour of reinforced concrete sections and members, under generalized load or deformation histories, including cyclic loading, up to ultimate deformation. At the section level the model was of the bounding surface type and accounted for the coupling between the two directions of bending and between them and the axial direction. For the construction of the member tangent flexibility matrix on the basis of the section tangent flexibility matrix, a piecewise-linear variation along the member was assumed for the nine terms of the tangent section flexibility matrix. Spacone et al. [11] presented a fibre beam-column element for the non-linear static and dynamic analysis of reinforced concrete frames. The element formulation was flexibility-based and relied on force interpolation functions that strictly satisfy the equilibrium of bending moments and axial force along the element. They proposed an iterative algorithm for the determination of the resisting forces during the element state determination, which was accurate and stable, even in the presence of strength loss, tracing very well the highly non-linear behaviour of R/C members under cyclic load combinations of bending moment and axial force. The technique based on the decomposition of a dense grid of cells is not numerically efficient, due to the large amount of information needed to characterise the section and the high number of operations required to implement the stress integration with an allowable error level, Miguel et al. [12]. These problems are more important in the structural analysis, during the computation of the global stiffness matrix and the evaluation of the internal forces. Moreover, this method can produce errors in the exact adaptation of the grid of cells to the section stress integration. If the neutral axis is unknown during the numerical procedure, the equilibrium may be achieved by iterating its location. In this case, the convergence
procedure is not assured, because the variation of the internal forces of the cross-section is discrete and depends on the size of the grid. To solve this problem, Rasheed and Dinno [13] proposed an algorithm which uses the secant stiffness approach for the direct prediction of the complete monotonic response of RC rectangular sections under symmetric bending moments. For biaxial bending moments and axial loads many studies have been carried out on the interaction surface of the cross-sections. Worth particular mention is the research by Brondum [14,15] and Yen [16], who presented methods for the analysis of arbitrary cross-sections, but using a rectangular stress block. In addition, Kawakami et al. [17] presented a complete algorithm for the analysis of RC members, including the effect of both regular and prestressed reinforcements, but their method was very difficult to implement in a computer program. Subsequently, Rodrıguez-Gutierrez and Dario Aristizabal-Ochoa [4] presented a method to determine the biaxial interaction diagrams for any orientation of the neutral axis of RC short columns with any geometry, obtaining analytical closed form expressions for equilibrium forces. The same authors [18,19] have latterly presented a method to determine (for arbitrary sections) an analytical model of moment–prestress–curvature (M– P–/) and also for slender columns. Recently, Fafitis [20] developed a method for the computation of the interaction surface of reinforced concrete sections subjected to axial loads and biaxial bending based on Green’s theorem. The integration of the stresses over the compressed region of the concrete section was transformed into line integrals along the compressive perimeter of this region. This method is suitable when the integrated function is uniform over one direction (as is the case where the section is under an increasing monotonic load), where the points located in parallel lines to the neutral axis are under the same stress level. This method is therefore not applicable to sections that do not meet this condition. When a reinforced concrete structure is analysed using numerical methods, the integrating functions in the cross-sections might not be constant along the lines parallel to the neutral axis. This occurs when constitutive equations that take into account the different behaviour of concrete in loadunload processes are implemented. The cyclic loads or the stress distribution processes in the structure may cause some regions of the section experience permanent strains, which are responsible for a non-uniform stress distribution in the lines parallel to the neutral axis. In contrast to the classical cells or layers methods, Miguel et al. [12] developed a numerical algorithm based on the decomposition of the non-cracked part of the section in quadrilaterals. In each of them, the integration of the section was performed by a Gauss–Legendre quadrature.
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
In addition, the same authors presented an alternative integration procedure to improve the accuracy for the same runtime, Bonet et al. [21]. Hence, in the case that the stress function to integrate does not approximate to a low order polynomial or is defined by steps, it would be advisable to subdivide the non-cracked zone by parallel ‘‘thick layers’’ to the neutral axis (Fig. 1). Each one of these layers will configure a polygon that will be decomposed into quadrilaterals by the previously defined algorithm. This integration method will henceforth be referred to as thick layers integration (TLI). It is presented in the next section. Table 1 shows a summary of the integration methods outlined above. This paper presents the TLI method (for non-cylindrical stress fields) and a modified thick layer integration (MTLI) method for the analysis of reinforced concrete
215
sections with any geometric shape (including holes). The algorithms are implemented for combined biaxial bending moments with axial forces, from any load level: service or ultimate stress. The second method is suitable only for sections in which the stress field is uniform over one direction, and it is based on the Gauss–Legendre quadrature. The fundamentals of the methods are explained in the paper and various possibilities for this purpose are analysed. In the proposed formulation of interaction strength estimations plain sections remain plane and perfect bond and lack of shear postulates are implicit. The obtained results for the different alternatives will be compared in accuracy and in speed with the results obtained from the fiber decomposition models and the method proposed by Fafitis [20].
Quadrilateral 1
N
Quadrilateral
Quadrilateral 2 Quadrilateral 3
Type
1
Hollow
2
Solid
3
Solid
εc Reference line.
Strain (a)
Stress
σc
(b)
Fig. 1. Integration with ‘‘thick layers’’ (TLI): (a) decomposition in ‘‘thick layers’’ and (b) decomposition in quadrilaterals.
Table 1 Summary of the integrations methods Classical integration methods of the stress field Method
Procedure
Layers decomposition method
• • • •
Cells decomposition method
The section is decomposed in ‘‘nb’’ layers The stress and Et is evaluated in the center of each layer The section is decomposed in a mesh of ‘‘n m’’ cells/fibers The stress and the tangent elasticity modulus (Et) are evaluated in the centre of each cell
Methods based on the Gauss–Lengendre quadrature Fixed integration points
The Gauss points are ‘‘fixed’’ and independent of the neutral axis position
Variable integration points Non-cylindrical stress fields Without layers decomposition (Miguel et al.) With layers decomposition Cylindrical stress fields Without layers decomposition (Fafitis)
With layers decomposition
Steps: • Decomposition in ‘‘Nc ’’ quadrilaterals • Double integration of each quadrilateral Thick layer integration method (TLI) Steps: • Decomposition in ‘‘nL’’ quadrilaterals • Double integral fi Path integral Modified thick layer integration (MTLI)
216
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
2.1. Decomposition of the integration area
The algorithm has limitations in some specific cases, where the problem is non-monotonic even under monotonic bending increase and constant axial load, due to e.g.:
• Interaction surface estimation. • Cross-sectional analysis. • Moment–curvature diagrams computation.
The integration area is decomposed into quadrilaterals to integrate stress fields in each region by a trivial quadrature method. Once the polygonal contour of the integration area has been determined, a reference line is chosen, Fig. 2a, and the orthogonal projections points from the vertexes of the contour are obtained (Fig. 2b). The section is then decomposed into quadrilaterals and triangles from the two consecutive vertexes and its projections, Fig. 2c. In this way, quadrilaterals and triangles with ‘‘positive’’ or ‘‘negative’’ areas are obtained, Fig. 2d.
At present, the authors are working in the inclusion of a new monitored variable to avoid previous limitations.
2.2. Internal force and constitutive matrix of the quadrilaterals
• Concrete degradation of reinforcement buckling. • Stress re-distribution in structures. But the proposed method performs adequately in:
The stress integration in each quadrilateral is performed using a Gauss quadrature. The quadrilaterals can be transformed into a square of 2 · 2 units by side, Zienkiewicz [22], using a coordinate transformation. The internal forces and tangent constitutive matrix of the quadrilaterals are obtained trivially. The sign of the Jacobian states whether the quadrilateral is solid or hollow. To this end the vertex is enumerated in a clockwise direction (positive) or anti-clockwise (negative). The proposed algorithm implements this enumeration automatically. The stresses are obtained by the
2. Thick layers integration method The contour where the stresses must be integrated can be considered as an n-sides polygon. If the section is made up of curved contours, it will also be reduced by polygons. The integration procedure of the section can be divided in three phases or steps: decomposition of the integration area in quadrilaterals or triangles, stress integration of each quadrilateral and layer subdivision.
Quadrilateral- 1 Integration Area
Quadrilateral- 2 2
3
2
4
4 - (4’)
Reference Line
Reference Line
(3’)
(2’)
1 — (1’)
1
Quadrilateral- 3
3
(a)
(c)
(b) QUADRLATERAL NUMBER
VERTEX
TYPE
2 SOLID DEGENERATED QUADRILATERAL
1 3
1
4 2
1
3
2 4 3
1 2
4 3
SOLID NON DEGENERATED QUADRILATERAL EMPTY DEGENERATED QUADRILATERAL
(d) Fig. 2. (a) Integration area and reference line, (b) vertex projection, (c) quadrilaterals determination and (d) vertex enumeration.
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
addition of each stress that actuates over the quadrilaterals. In addition, the constitutive matrix of the sections is obtained by composition of the quadrilaterals. 2.3. Thick layers decomposition Certain considerations must be taken into account in the proposed method. First, the utilisation of non-linear structural analysis programs for reinforced concrete structures establishes the use of iterative algorithms. Hence, during the convergence process its is observed that the neutral axis and the stress field are different for all the cross-sections into which the beams of the structure are divided. Furthermore, the cracking phenomenon in the concrete characterises the behaviour of this type of structure. In the cracked region of the cross-section, the concrete fibers do not add stiffness to the structure. Moreover, the integration area may be made up of all the cross-section or only of the regions with a non-zero stress field. In the first case, the location of the Gauss points will be fixed during the whole convergence process. Hence, the decomposition of the integration area will be performed once during the procedure. Otherwise, if the stress field influences the integration area, it will be redefined in each iteration. Hence, the allocation of the integration points will change during the process. The accuracy of the integrals depends for each quadrilateral on the number of Gauss points used and the polynomial form of the stress–strain relation. When the stress field or the elasticity modulus functions are not adjusted to a small order polynomial, a substantial improvement in the integration definition is brought about by decomposing the area into non-zero stress regions, by ‘‘thick layers’’ parallel to the neutral axis (Fig. 1b). In this situation, each layer defines an n-sides polygon where the proposed method can be applied.
217
The smaller deformation point of the layer may define the reference line. In consequence, the internal forces or the constitutive matrix are obtained by summing the terms regarded with the quadrilaterals of each layer into which the integration area is divided.
3. Modified thick layers integration method The second method is suitable only for sections in which the stress field is uniform over one direction, what can be called ‘‘cylindrical stress field’’, and in this way the stresses are integrated for arbitrary geometry RC cross-section with biaxial moments and axial forces. Again, the contour where the stresses must be integrated can be considered as an n-sides polygon. The proposed procedure can be divided into three phases: decomposition into ‘‘thick layers’’ of the specific area of integration (non-zero stress field of the concrete fibers), transformation of the stress integral for each thick layer into a ‘‘path’’ integral over its perimeter, and finally the evaluation of the ‘‘line’’ integral of stresses with a Gauss quadrature. If the integration function is defined by a curve that approximates adequately in its total range to a polynomial function, the subdivision of the integration area into thick layers is unnecessary. 3.1. Transformation of the stress integral into a ‘‘path integral’’ over the perimeter of the concrete integration area If the stress field is uniform in a particular direction (rc ðy; zÞ ¼ rc ðzÞ; Ect ðy; zÞ ¼ Ect ðzÞ), the stress integral of the integration area can be transformed into a path integral over the perimeter of the integration area (Fig. 3). The internal forces can be obtained by the expressions:
z
z
2 Integration Area
L1
L2
y dz (y/2-y0)
3
(z-z0)
dz
1 dy y
(y0, z0)
L3
Reference Line
y L5
4
L4
Fig. 3. Transformation of the stress integral.
5
Reference Line
218
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
Nc ¼ ¼
Z Z Z ZAc
rc ðy; zÞ dy dz
3.2. Evaluation of the stress integral by the Gauss quadrature
rc ðzÞ dy dz
As the integration area contour is approximated by a polygon, the stress field integral over the perimeter L, can be obtained by decomposing this integral side by side through the perimeter (Fig. 4).
Ac
I
y rc ðzÞ dz ZL Z ðz z0 Þ rc ðzÞ dy dz Mcy ¼ Ac I ¼ ðz z0 Þ y rc ðzÞ dz ZL Z ðy y0 Þ rc ðzÞ dy dz Mcz ¼ Ac I ¼ ðy=2 y0 Þ y rc ðzÞ dz ¼
ð1Þ
Z Z
¼ f ðy; zÞ dy dz ¼ Ac
¼
Z ZAc I
SLi ðzÞ
ð3Þ
i¼1
SLi ðzÞ ¼
Z
ziþ1
hðzÞ dz zi
In order to perform the stress integral of a determined side of the contour of the integration area (Li ), the wellknown Gauss–Legendre method is used. The next coordinate transformation must be implemented:
Ect ðy; zÞ dy dz Ect ðzÞ dy dz
Ac
y Ect ðzÞ dz ZL Z ðz z0 Þ Ect ðzÞ dy dz Dct ð1; 2Þ ¼ I Ac ¼ ðz z0 Þ y Ect ðzÞ dz ZL Z ðy y0 Þ Ect ðzÞ dy dz Dct ð1; 3Þ ¼ I Ac ¼ ðy=2 y0 Þ y Ect ðzÞ dz ZL Z ðz z0 Þ2 Ect ðzÞ dy dz Dct ð2; 2Þ ¼ Ac I ¼ ðz z0 Þ y Ect ðzÞ dz ZL Z ðz z0 Þ ðy y0 Þ Ect ðzÞ dy dz Dct ð2; 3Þ ¼ Ac I ¼ ðz z0 Þ ðy=2 y0 Þ Ect ðzÞ dz ZL Z ðy y0 Þ2 Ect ðzÞ dy dz Dct ð3; 3Þ ¼ Ac I ¼ ½y 2 =12 þ ðy=2 y0 Þ2 y Ect ðzÞ dz ¼
nL X
where nL is the number of sides that forms the integration area:
The elements of the tangent stiffness matrix of the section (Dct ) can also be obtained: Z Z
hðzÞ dz ¼
L
L
Dct ð1; 1Þ ¼
I
z¼
ðziþ1 zi Þ ðziþ1 zi Þ ðn þ 1Þ þ zi ) dz ¼ dn 2 2
ð4Þ
The stress field integral of each side of the contour is
SLi ðzÞ ¼
Z
ziþ1
hðzÞ dz ¼ zi
ð2Þ
L
where Ac is the integration area, L is the contour of this area, rc is the concrete stress at each fiber (y; z), Ect is the tangent elasticity modulus of deformation of the fiber (y; z); and (y0 ; z0 ) are the coordinates of the centre of reference of forces of the section. The ‘‘y’’ axis is situated parallel to the direction in which the stresses are uniform.
¼
ðziþ1 zi Þ 2
Z
þ1
hðnÞ dn
1
n ðziþ1 zi Þ X ki hðni Þ 2 i¼1
ð5Þ
where ni is the value of the curvilinear coordinate of the Gauss point, ki is the weight associated to the Gauss point and ‘‘n’’ is the number of Gauss points.
z (yi+1,zi+1) y
Li dz
(y/2-y0) (z-z ) 0
(y0, z0)
y = yi +
(yi -yi) (z -z ) i (zi+1-zi)
(yi-zi)
y Reference Line
Fig. 4. Evaluation of the stress integral.
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
40 5
ð6Þ The proposed algorithm allows us to integrate the stresses in a hollow section or in a solid section. Hence, the vertexes of the quadrilateral must be enumerated again in counter-clockwise (positive sign) or clockwise (negative sign).
20
Lj
5
L
" # nL n X ðziþ1 zi Þ X ¼ ki hðni Þ 2 j¼1 i¼1
90
Ac
20
In this way, the integral is I Z Z f ðy; zÞ dy dz ¼ hðzÞ dz
219
3.3. Stresses integration through ‘‘thick layers’’ The accuracy of the integrals depends for each one of the sides on the number of Gauss points used and the polynomial form of the stress–strain relation. When the function to integrate does not follow a low order polynomial, or is defined as a steps function, a large number of Gauss points has to be used to obtain an acceptable accuracy. In these cases it is better to subdivide the integration area into thick layers parallel to the neutral axis (Fig. 1a). Each one of the layers define an n-sides polygon where the proposed integration method can be performed. Subsequently, the application of this method to the classical stress–strain relationship of the concrete (Fig. 1a) may be performed using five layers. Hence, the internal forces and the constitutive matrix of the section can be obtained as the sum of the terms of each side that composes each perimeter of each ‘‘thick layer’’ into which the integration area has been decomposed. This method will be hereafter referred to as modified thick layers integration (MTLI).
4. Numerical examples The methods described above were tested in three different sections, in which the strength of concrete, and the depth and the inclination angle of the neutral axis were modified with the purpose of integrating the stress field. Four methods were implemented: the classical fiber decomposition approach, the integration method proposed by Fafitis [20], the TLI and the MTLI proposed in this paper. The comparison between the different methods was performed in terms of the time and the accuracy obtained. The sections tested were: 20 · 20 cm rectangular; hollow 20 · 20 cm rectangular with 5 cm thick walls and a section in ‘‘T ’’ (Fig. 5). The ‘‘exact’’ solution is considered as the solution solved with the MTLI method with 48 Gauss points for
15
5
20
5
15
60 Fig. 5. Section in ‘‘T ’’ (cm).
each side (Li). Two strengths of concrete were considered for all the cases, 25 and 80 MPa. The constitutive equation used corresponds to the Model Code 1990 [23], and also includes the tensile behaviour and the tension stiffening effect through a model of gradual descending unload, CEB-1995, [24]. In this case the concrete has a cylinder strength distribution. The tensile strain descending branch shape depends on the stress of the most tensile stressed reinforced bar of the section. In the numerical implementation of this study the steel contribution has been ignored because it does not influence the accuracy and the runtime of the integration process. Due to this, the equation adopted for this branch will be a straight line from the maximum tensile stress (fct ) and the ultimate strain deformation for the concrete, which corresponds to the plastic steel deformation (ecend ¼ fy =Es ), Fig. 6. For each case, the stress distribution was integrated assuming that the ultimate compressed concrete fiber had a deformation of 0.0035 and the most tensile stressed fiber of )0.0025 (Fig. 7).
σc
Eci 1
fc
Ec1
ε c,end
1
εc εc1
εc,lim
-fct Fig. 6. Concrete stress–strain relationship adopted.
220
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
For each angle, three depths of the most tensile stressed fiber were computed, corresponding to 1/3, 2/3 and 3/3 of the maximum height (hmax ) orthogonal to the neutral axis (Fig. 7).
hmáx 3.5
Fiber under tension
εc
-2.5
Strain
σc
Stress
Fig. 7. Stress distribution.
To perform the integration, the ‘‘TLI’’ method was used with a mesh of points of 2 · 2, 3 · 3, 4 · 4, 6 · 6, 8 · 8, 10 · 10, 12 · 12, 15 · 15, 20 · 20 and 32 · 32. In each one of the quadrilaterals where the integration area was divided, and for the Fafifis method and the ‘‘MTLI’’ method, the following set of points was used: 2, 3, 4, 6, 8, 10, 12, 15, 20 and 32, for each side of the perimeter of the total integration area or the perimeter of each thick layer that forms the integration area, respectively. For the fiber integration approach, a variable mesh was used from 20 · 20 to 1200 · 1800 cells, depending on the type of section, in order to obtain a better adaptation of the mesh to the geometry of the section (Table 2). Different inclination angles from the neutral axis were tested, depending on the type of section (Table 3). Table 2 Cells mesh Section type
Cell mesh
Rectangular
20 · 20, 40 · 40, 100 · 100, 200 · 200, 400 · 400 20 · 20, 40 · 40, 100 · 100, 200 · 200, 400 · 400 60 · 90, 120 · 180, 300 · 450, 600 · 900, 1200 · 1800
Hollow rectangular Section in ‘‘T ’’
Table 3 Position of the last fiber under tension Section type
Inclination of the neutral axis
Depth
Rectangular
0, 30, 45
Hollow rectangular
0, 30, 45
Section in ‘‘T ’’
0
hmax , 1/3 hmax , 2/3 hmax hmax , 1/3 hmax , 2/3 hmax hmax , 1/3 hmax , 2/3 hmax
5. Results In this section the previously computed numerical examples corresponding to the three section types are analysed. In each example, the error and the runtime are calculated. It was observed that the computing integration time using quadrature was extremely small and lower than computer clock accuracy. Hence, the runtime was measured for the computation of 10 test series. Each series repeats the analysis of the section 50 times. This procedure was repeated for the four algorithms. The error obtained in each method in relative percentage to the exact value is Fi Fexac 100 Errorð%Þ ¼ Fexac
ð7Þ
where Fi is the exact result obtained in case ‘‘i’’ for force F (axial or bending) and Fexac is the exact value for force F (axial or bending). Table 4 shows, for all the types of sections and number of Gauss points tested, the runtime and the maximum errors obtained for all forces, corresponding to the method proposed by Fafitis. It should be noted that the errors for a small number of Gauss points (2, 3, 4 and 6) are not allowable, principally for the 80 MPa concrete. This result appears because the numerical integration using the Gauss– Legendre quadrature of any function f ðxÞ is based on the interpolation of the function with an n-order polynomial, in such a way that the Gauss–Legendre quadrature of n-order, exactly integrates a polynomial of order 2n 1 or lower. Because the stress field is defined by a step function and there is no continuity in the derivative, the polynomial interpolation can produce important integration errors. On the other hand, it can be observed that the errors increase when the strength of the concrete increases. This behaviour comes from the form of the constitutive equation of the concrete where the deformation of the maximum stress is independent of the concrete strength (Fig. 8). However, since the stress function takes on a narrower shape as the strength increases, a higher number of Gauss points are required for sections with 80 MPa concrete than for sections with 25 MPa concrete to obtain the same degree of accuracy.
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
221
Table 4 Results of the Fafitis method Section type
fc (MPa)
Number of Gauss points 2
Rectangular
25 80 25 80 25 80
Hollow rectangular Section in ‘‘T ’’
3
4
6
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
1.32 1.32 2.57 2.64 4.17 4.29
12.47 65.88 12.47 65.98 6.43 81.09
1.65 1.65 3.13 3.07 5.21 5.21
10.50 19.3 10.50 19.3 3.70 32.95
1.92 1.92 3.79 3.85 6.32 6.26
3.57 13.39 3.57 13.39 1.63 29.03
2.63 2.59 4.84 4.99 7.96 8.07
0.22 7.09 0.22 7.09 0.16 21.66
Fig. 8. Uniaxial behaviour of the concrete in compression.
If the integration area is decomposed into thick layers, and each thick layer is also subdivided into quadrilaterals (TLI method), for a small number of Gauss points (2 · 2, 3 · 3, 4 · 4 and 6 · 6) the error is reduced considerably (Table 5). In this way, it can once again be observed that the errors of the TLI method rise when the concrete strength increases. However, the computation time increases notably compared with that of the Fafitis method. If the stress field of the section is uniform in a fixed direction, the ‘‘TLI’’ method does not optimise the
location of the Gauss points, because the integral in this direction can be obtained analytically (Eq. (6)). For this reason, the runtime increases in comparison with the Fafitis method. Through the proposed ‘‘MTLI’’ method, the location of the Gauss points is optimised, and hence, the runtime decreases considerably (Table 6), obtaining an accuracy level similar to the ‘‘TLI’’ method. Table 7 compares the degree of accuracy and the computation time of the three Gaussian methods versus the number of Gauss points for the 80 MPa concrete strength. It can be observed that although, in general, a higher number of Gauss points allows us to obtain a higher accuracy, the results obtained through the Fafitis method do not follow this behaviour due to the irregularity that exists in the stress field of the section. On the other hand, it can be observed that the errors in the integration methods ‘‘TLI’’ and ‘‘MTLI’’ are similar. However, the runtime of the proposed method is notably lower than that of the ‘‘TLI’’ method, and this behaviour is relatively important as the number of Gauss points increases in the integration. Table 8 shows the errors and the computation time for the rectangular section (both hollow and solid), corresponding to the fiber decomposition method. In Table 9 the results of the section in ‘‘T ’’ are presented. It can be observed that although the relative errors obtained are small, the computation time is higher
Table 5 Results of the ‘‘TLI’’ method Section type
Rectangular Hollow rectangular Section in ‘‘T ’’
fc (MPa)
25 80 25 80 25 80
Number of Gauss points 2·2
3·3
4·4
6·6
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
11.48 12.25 20.43 23.06 12.57 13.61
0.10 0.41 0.13 0.45 0.12 0.80
21.91 22.84 37.40 42.91 19.89 21.42
0.0014 0.045 0.00172 0.0485 0.0016 0.39
34.77 37.02 60.37 69.43 29.98 32.78
2.16e)5 0.0022 2.16e)5 0.0018 2.13e)5 0.145
72.02 76.73 124.62 142.52 58.06 63.88
2.89e)9 0.0052 3.6e)9 0.00069 3.5e)9 0.016
222
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
Table 6 Results of the proposed ‘‘MTLI’’ method Section type
fc (MPa)
Number of Gauss points 2
Rectangular Hollow rectangular Section in ‘‘T ’’
25 80 25 80 25 80
3
4
6
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
4.39 6.21 5.66 9.12 12.19 13.34
0.54 0.25 0.54 0.34 0.12 0.80
5.94 8.02 7.25 11.76 15.32 16.75
0.0068 0.1223 0.0068 0.244 0.0016 0.39
6.97 10 8.95 14.28 18.62 20.38
9.58e)5 0.0448 9.58e)5 0.10 2.43e)5 0.145
9.45 13.89 11.92 19.27 24.71 27.24
1.4e)6 0.00494 1.4e)6 0.012 3.5e)9 0.016
Table 7 Results (rectangular section; neutral axis angle: 45; depth: hmax ; fc ¼ 80 MPa) Gauss points
Integration method Fafitis
2 3 4 6 8 10 12 15 20 32
‘‘TLI’’
‘‘MTLI’’
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
2.03 2.69 3.41 4.67 5.77 7.2 8.51 10.39 13.56 21.26
17.08 6.88 0.73 3.54 0.42 1.56 0.52 0.77 0.25 0.07
13.84 26.03 42.4 87.72 150.22 229.32 326.14 502.29 880.51 2219.1
0.38 0.04 0.0015 0.00053 6.12e)05 5.74e)06 0.52e)6 5.8e)09 0 0
5.38 7.04 8.63 11.92 15.16 18.35 21.86 26.47 34.55 54.05
0.148 0.158 0.068 0.0082 0.00088 9.06e)05 9.23e)06 2.88e)07 1e)09 0
Table 8 Results of the fiber decomposition method Section type
fc (MPa)
Cell mesh 10 · 10
20 · 20
100 · 100
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
Rectangular
25 80
22 22
2.19 4.63
82.5 82
0.43 2.30
519 513.5
0.035 0.51
Hollow rectangular
25 80
19 16.5
2.19 4.64
66 69
0.43 2.30
406 412
0.035 0.51
Table 9 Results of the fiber decomposition method Section type
fc (MPa)
Cell mesh 60 · 90
Section in ‘‘T ’’
25 80
120 · 180
360 · 540
t (sg)
Error %
t (sg)
Error %
t (sg)
Error %
192 195
0.056 0.52
758 774.5
0.016 0.023
4737.5 4847
0.0037 0.008
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
223
25 MPa 100 Error (%)
10
Cells
1
Fafitis
0.1
TLI
0.01
MTLI
0.001 0.0001 1E-05 1E-06 1E-07 1E-08 1E-09 1E-10 1
10
100
1000
10000
100000
Time (sg) 80 MPa Error (%)
100 10
Cells
1
Fafitis TLI
0.1 0.01
MTLI
0.001 0.0001 1E-05 1E-06 1E-07 1E-08 1E-09 1E-10 1
10
100
1000
10000
100000
Time(sg)
Fig. 9. Efficiency of the fours methods for 25 and 80 MPa for a rectangular section with a ¼ 45 and depth ¼ hmax .
than that obtained with the ‘‘TLI’’ and the ‘‘MTLI’’ method. In general, as the mesh density increases, the error decreases (Tables 8 and 9). However, the accuracy of this method depends on how the cell mesh adapts to the integration area, in that a mesh with a higher number of cells, but distributed in worse conditions respective to the geometry of this area will produce higher errors. In order to compare the four methods, Figs. 9 and 10 are presented. These figures present the percentage error on the Y -axis, and the integration time on the X -axis. The objective of a good algorithm is to achieve lower time consumption with fewer errors. Fig. 9 shows one numerical example in particular: rectangular section with a ¼ 45 and the depth of the neutral axis to the maximum. When the three Gaussian methods are compared with the cell method, it can be noted that the ‘‘MTLI’’ method is the most efficient. This is followed by the ‘‘TLI’’ method, and then the Fafitis method. The Fafitis method is faster than the MTLI but presents inadmis-
sable errors. MTLI is slower because it spends time decomposing the integration area into layers. For a nonlinear numerical procedure, the error must be fixed to arrive at the convergence, which is what reduces the efficiency of the Fafitis method. For example, for the case of Fig. 9 with 25 MPa, and for an admissible error of 0.001%, the ‘‘MTLI’’ method is 2.7 times faster than the Fafitis method, and 1.3 times faster than the ‘‘TLI’’ method. The difference between the TLI and MTLI methods and the Fafitis method is that the former decompose the stress function by trams and in consequence the interpolation polynomial adapts better to the stress function, optimising the properties of the Gaussian integration. It can also be noted that the time integration for the numerical examples of the 80 MPa examples are slower than for the 25 MPa examples. The reason is that the stress–strain relationship of the 80 MPa presents a narrower behaviour and it needs a higher order polynomial integration and in consequence more integration points and more time for the same accuracy.
224
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225 25 MPa 100
Error (%)
10 1 0.1 0.01 0.001
CELLS
0.0001 1E-05
FAFITIS
1E-06 1E-07
MTLI
TLI
1E-08 1E-09 1E-10 1
10
100
1000
10000
100000
Time (sg) Fafitis
Cells
MTLI
TLI
80 MPa 100 10 1 0.1
Error (%)
0.01
CELLS
0.001 0.0001 1E-05 1E-06
TLI
MTLI
FAFITIS
1E-07 1E-08 1E-09 1E-10 1
10
100
1000
10000
100000
Time (sg) Fafitis
Cells
MTLI
TLI
Fig. 10. Efficiency of the fours methods for 25 and 80 MPa for all the examples computed.
6. Conclusions This paper deals with the implementation of two integration procedures for arbitrary geometry reinforced concrete cross-sections with axial forces and biaxial bending, from service load to ultimate load. Two methods are proposed: the TLI method is suitable for a general stress field, and the MTLI method is applicable in sections in which the stress field is uniform at least in one direction. Both methods decompose the integration area into thick layers parallel to the most tensile stressed fiber, whose definition depends on the constitutive equation of the concrete. The integration of the stress field of each thick layer is transformed into a double (TLI) or a path integral (MTLI) over the perimeter of this layer, and they are evaluated by a Gauss quadrature. In the study of the four methods analysed (fibers, Fafitis [20], ‘‘TLI’’ and ‘‘MTLI’’), applied to different concrete section types, the following conclusions are reached:
(a) In all the cases, for an equivalent computation time, the accuracy of the Gaussian integration methods is notably higher than for the fiber method. (b) The accuracy of the Gaussian integration methods increases in relation to the number of Gauss points used, except in the method proposed by Fafitis [20]. The accuracy level of this method depends on the form of the function used to represent the concrete stress–strain relationship. (c) The accuracy level of the fiber decomposition method does not depend only on the mesh density but also on its adaptation to the geometry of the section. (d) Due to the form adopted by the concrete constitutive equation defined by the Model-Code 90 (which depends on the strength), the accuracy level is reduced as the strength of concrete is increased for all the integration methods. (e) Both the thick layer integration (TLI) and modified thick layer integration (MTLI) method show an excellent accuracy level (lower than 0.02%) for a 6 · 6 Gauss mesh points per quadrilateral or 6 Gauss
J.L. Bonet et al. / Computers and Structures 82 (2004) 213–225
points per side of the perimeter that forms the thick layer. (f) When the accuracy level and the runtime for the three Gaussian methods are compared with the fiber decomposition method, it is observed that the ‘‘MTLI’’ method is more efficient, followed by the ‘‘TLI’’ method and finally the Fafitis method. Finally, the proposed ‘‘modified thick layers’’ method (MTLI), with regard to the accuracy, efficiency, and continuity in the stress field integration is recommended for implementation in non-linear reinforced concrete frameworks programs.
Acknowledgements This paper was supported by the Spanish ‘‘Ministerio de Ciencia y Tecnologia’’ under project MAT200202461.
References [1] Bonet JL. Metodo simplificado de calculo de soportes esbeltos de hormig on armado de secci on rectangular sometidos a compresi on y flexi on biaxial, PhD thesis, Civil Engineering Department, Technical University of Valencia, July 2001 [in Spanish]. [2] Romero ML, Miguel PF, Cano JJ. A parallel procedure for nonlinear analysis of reinforced concrete three dimensional frames. Comput Struct 2002;80(16–17):1337–50. [3] Br€ ondum-Nielsen T. Serviceability analysis of concrete sections under biaxial bending. ASCE, J Struct Eng 1997;123(1):117–9. [4] Rodrıguez-Gutierrez JA, Dario Aristizabal-Ochoa J. Biaxial interaction diagrams for short RC columns of any cross section. J Struct Eng, ASCE 1999;125(6):672–83. [5] Mari AR. Nonlinear geometric, material and time dependent analysis of three dimensional reinforced and prestressed concrete frames. Report No. USB/SESM-84/12, Department of Civil Engineering, University of California, Berkley, California, USA, June 1984. [6] Lai SS, Will GT, Otani S. Model for inelastic biaxial bending of concrete members. J Struct Eng-ASCE 1984;110(11):2563–84. [7] Kaba S, Mahin SA. Refined modeling of reinforced concrete columns for seismic analysis. Report EERC 8403, Earthquake Engineering Research Center, University of California, Berkeley, California, 1984.
225
[8] Zeris CA, Mahin SA. Analysis of reinforced-concrete beam-columns under uniaxial excitation. J Struct EngASCE 1988;114(4):804–20. [9] Zeris CA, Mahin SA. Behavior of reinforced-concrete structures subjected to uniaxial excitation. J Struct EngASCE 1991;117(9):2640–56. [10] Bousias SN, Panagiotakos TB, Fardis MN. Modelling of RC members under cyclic biaxial flexure and axial force. J Earthquake Eng 2002;6(2):213–38. [11] Spacone E, Filippou FC, Taucer FF. Fibre beam-column model for non-linear analysis of R/C frames. 1. Formulation. Earthquake Eng Struct 1996;25(7):711–25. [12] Miguel PF, Bonet JL, Fernandez MA. Integraci on de tensiones en secciones de hormig on sometidas a flexocompresi on esviada. Rev Int Metodos Numer para el calculo dise~ no ingenierıa 2000;16(2):209–25 [in Spanish]. [13] Rasheed HAS, Dinno KS. An efficient nonlinear analysis of RC sections. Comput Struct 1994;53(3):613–23. [14] Brondum T. Ultimate limit states of cracked arbitrary cross sections under axial loads and biaxial bending. ACI Concr Int 1982;4(11):51–5. [15] Brondum T. Ultimate flexural capacity of cracked polygonal concrete sections with circular holes under biaxial bending. ACI Concr Int 1987;84(3):212–5. [16] Yen JR. Quasi-Newton method for reinforced concrete column analysis and design. J Struct Eng, ASCE 1991; 117(3):657–66. [17] Kawakami MT, Kagaya M, Hirata M. Limit states of cracking and ultimates strength of arbitrary cross sections under biaxial loading. ACI J 1985;82(1):203–12. [18] Rodriguez JA, Aristazabal-Ochoa JD. M–P–/ diagrams for reinforced, partially, and fully prestressed concrete sections under biaxial bending and axial load. J Struct Eng, ASCE 2001;127(7):763–73. [19] Rodriguez JA, Aristazabal-Ochoa JD. Reinforced, partially, and fully prestressed slender concrete columns under biaxial bending and axial load. J Struct Eng, ASCE 2001;127(7):774–83. [20] Fafitis A. Interaction surfaces of reinforced-concrete sections in biaxial bending. J Struct Eng, ASCE 2001; 127(7):840–6. [21] Bonet JL, Miguel PF, Fernandez MA, Romero ML. Efficient procedure for stress integration in concrete sections using a Gauss–Legendre quadrature. In: Topping BHV, editor. Proceedings of the Eighth International Conference on Civil and Structural Engineering Computing. Stirling, United Kingdom: Civil-Comp Press; 2001. p. 53. [22] Zienkiewicz OC, Taylor RL. In: The finite element method, vol. 1. CIMNE, Mc-Graw-Hill; 1995. [23] Comite Euro-internacional du beton: CEB-FIB Model Code 1990. CEB Bull 1991;203–5. [24] Comite Euro-internacional du beton: new developments in non-linear analysis methods. CEB Bull 1995;229:94.