Computers and Structures 98-99 (2012) 33–54
Contents lists available at SciVerse ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
Analysis of arbitrary composite sections in biaxial bending and axial load Vassilis K. Papanikolaou
⇑
Laboratory of Reinforced Concrete and Masonry Structures, Civil Engineering Department, Aristotle University of Thessaloniki, P.O. Box 482, Thessaloniki 54124, Greece
a r t i c l e
i n f o
Article history: Received 1 November 2011 Accepted 13 February 2012 Available online 11 March 2012 Keywords: Reinforced concrete Composite sections Biaxial bending Stress integration Ultimate strength Moment–curvature
a b s t r a c t
A new methodology is presented for the ultimate strength and moment–curvature analysis of arbitrary composite sections under biaxial bending and axial load. The definition of section geometry and material properties can be unconditionally complex, based on an object-oriented implementation. Stress integration is performed using a Green path integral, with an adaptive strain-mapped Gaussian sampling. Derivative-free solution strategies for the calculation of incremental and ultimate response are applied. Results are presented in the form of moment–curvature curves, ultimate strength interaction curves and 3D failure surfaces. The performance of the methodology is demonstrated through various case studies, comparisons and benchmarks. 2012 Elsevier Ltd. All rights reserved.
1. Introduction
Reinforced concrete and composite structural elements, widely used in buildings and bridges, are generally subjected to a combined action of biaxial bending and axial load. This is a result of their their section section geometry geometry and material material compositi composition, on, their their position position and orientation in the structure and, most importantly, the nature of external loading. Columns at corners or under two-way slabs, bridge piers and composite decks, subjected to wind or earthquake excitat excitation, ion, are represen representat tative ive cases. cases. Assessing Assessing the adequacy adequacy of these sections at ultimate limit state (usually via interaction diagrams) or providing information on their inelastic response gradually up to failure (in the form of moment–curvature curves), is a computationally intense task, mainly due to material nonlinearities and geometrical complexities, particularly for composite sections. tions. Theref Therefore ore,, addre addressi ssing ng this this proble problem m requi require ress efficie efficient nt solution algorithms, characterized characterized not only by the requisite robustness but also by execution speed, since section analysis is often a repetitive task [1] within the broader computational frameworks of nonlinear structural analysis, design, assessment and retrofit. The problem of arbitrary arbitrary section analysis has received attention in the literature since the early 1960’s (e.g. [2–4] [2–4]), ), however this attention has been intensified during the last decade along with the advent of inexpensive personal computers. Various analytical, numerical and mixed methodologies for analyzing analyzing sections of varying complexity complexity in terms terms of geometry geometry and material material composition composition have been suggested, with varying degrees of reported efficiency ⇑
Tel.: +30 2310995662. E-mail address:
[email protected]
0045-7949/$ - see front matter 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2012.02.004 doi:10.1016/j.compstruc.2012.02.004
in terms of convergence stability and speed. Since the literature is extensive, a proper categorization is attempted hereinafter, in order to assess whether the research field under discussion is open for further novelties and improvements. In this direction, five key characteristics characteristics of section analysis methodologies have been identified, fied, namely namely (a) the type type of arbit arbitrar rary y sectio section n addre addresse ssed, d, i.e. i.e. whether whether it is reinforce reinforced d concrete concrete (R/C) (R/C) or generall generally y composite composite,, (b) the form of material constitutive laws, e.g. polynomial or arbitrary functions, (c) whether section subdivisions are imposed prior to stress integration, (d) the stress integration scheme and (e) the solution strategies for applying force equilibrium conditions. A non-exhaustive yet representative representative directory of previous studies on the subject is summarized in Table 1, 1 , providing specific information on the aforementioned five key topics for each study. A critical review on this categorized literature summary leads to the following remarks: A significant number of suggested algorithms are limited to R/C sections, consisting of a single concrete surface and a group of individual fibers for reinforcement bars (i.e. two distinct materials). In order to bypass such limitations, the first prerequisite of the present methodology will be the unconditional section complexity, i.e. unlimited number of section components (surfaces, faces, fiber fiber groups groups and and – newly newly introduc introduced ed – lines) lines),, each each assigned to a different material constitutive law. Most studies impose restrictions on material stress–strain constitutive laws, usually in the form of specific Code directives (e.g. [5,6] (e.g. [5,6]), ), or more elaborate piecewise polynomial laws (e.g. [7,8]). [7,8] ). For this reason, the second prerequisite of the present meth method odol olog ogy y will will be the the use use of full fully y arbi arbitr trar ary y mate materi rial al
34
V.K. Papanikolaou Papanikolaou / Computers and Structures 98-99 (2012) 33–54
Table 1
Features of selected previous studies on the analysis of arbitrary R/C and composite sections.
Authors
Section type Material constitutive law Sect Sectio ion n subd subdiv ivis isio ion n for surfaces
Stre Stress ss inte integr grat atio ion n
Solu Soluti tion on stra strate tegi gies es
Charalampakis and Koumousis [8] Koumousis [8]
Compos Composite ite
Chen et al. [16] al. [16]
Composite
Chiorean [18] Chiorean [18]
Composite Composite
Rosati et al. [5] al. [5]
Comp Compos osit ite e
Sfakianakis [15] Sfakianakis [15]
Composite Composite
Sousa and Muniz [7] Muniz [7]
Compos Composite ite
Brøndum–Nielsen [25] Brøndum–Nielsen [25],, Dundar and Sahin [22] Sahin [22],, Yen [21] Yen [21] Bonet et al. [1] al. [1]
Reinforced concrete Reinforced concrete Reinforced concrete Reinforced concrete Reinforced concrete Reinforced concrete Reinforced concrete Reinforced concrete
Fafitis [17] Fafitis [17]
Pallarés et al. [6] al. [6]
Penelis [3] Penelis [3],, De Vivo and Rosati [26] Rosati [26],, Alfano et al. [27] al. [27] Rodriguez and Ochoa [20] Werner [4] Werner [4] Yau et al. [23] al. [23]
Sug Suggest gested e d met met ho hodolo dolog gy
Comp Compos osit ite e
Piecew Piecewise ise up to cubic cubic polynomials Parabolic-linear
Curvil Curviline inear ar trapez trapezoid oidss
Closed Closed-fo -form rm per trapez trapezoid oid
Deriva Derivativ tive-f e-free ree (Brent (Brent))
No subdivisions
Closed-form
Parabolic-li Parabolic-linear near with softening Para Parabo boli licc-li line near ar
No subd subdiv ivis isio ions ns
Green Green/G /Gau auss ss-L -Lob obat atto to with with adaptive bisection Closed Closed-fo -form rm per polygo polygon n
Derivative-free (regula-falsi) With derivatives (arc-length) With With deriva derivativ tives es (NR) (NR)
F ib iber integration
Incremental search
Clos Closed ed-f -for orm m per per poly polygo gon n
–
One One poly polygo gon n per per r e part No subdivisions
Parabolic-li Parabolic-linear near with softening Piecew Piecewise ise up to cubic cubic polynomials Constant (rectangular)
One polygon per r e part No subdivisions
Piecewise nonpolynomial (arbitrary) Parabolic-linear
Polygons (thick layers) per r e part No subdivisions
Non-polynomial – linear No subdivisions (EC2) Para Parabo boli licc-li line near ar One One poly polygo gon n per per r e part Parabolic-linear with Trap Trapez ezoi oids ds softening Para Parabo boli licc-li line near ar One One poly polygo gon n per per r e part Constant (rectangular) No No subdivisions
Piec Piecew ewis ise e nonnonpolynomial (arbitrary)
No subd subdiv ivis isio ions ns
constitutive laws in piecewise form, as reported only in Bonet et al. [1] [1].. This enables the analysis of sections with materials of non-polyn non-polynomia omiall stress–st stress–strain rain relation relationship shipss (e.g. (e.g. [9] for high-strength high-strength concrete and nonlinear analysis, [10] analysis, [10] for for confined concrete and [11,12] and [11,12] for for plain concrete). The most critical element that affects both the accuracy and speed of a section analysis algorithm is the stress integration scheme scheme.. There There are are three three major major paths paths that that can be follow followed, ed, namely (a) fiber integration (e.g. [13–15]), [13–15]), (b) analytical integration using closed-form functions (e.g. [8,16] (e.g. [8,16])) and (c) numerical integration, usually in a form of Gaussian sampling on a Green Green path integral integral (e.g. (e.g. [17,18]). [17,18]). The first approach is both approximate and slow, requiring an increased fiber mesh density and a proporti proportional onal number of arithmet arithmetic ic operatio operations ns to reach an acceptable level of accuracy. Therefore, this method is now obsolete for ultimate strength section analysis and is used only for non-cylindrical stress fields, e.g. in cases of cyclic loading and/or load path dependency [1] [1].. As far as the second approach is concerned, its main advantage is that it yields exact and quick results (especially for lower order functions), it is however totally restricted to a specific stress–strain expression, which is in conflict with the aforementioned second prerequisite. Consequently, the sole path to provide a generalized solution tion for for arbi arbitr trar ary y mate materi rial al cons consti titu tuti tive ve laws laws is the the implement implementatio ation n of a suitable suitable numerical integration integration scheme. However, it has been reported that, in certain cases, numerical integration may be more expensive than analytical methods, while low order numerical integration may yield unacceptably large errors [1,19] errors [1,19] because because the required order of numerical integration is not known a priori [8] [8].. These issues will be confirmed and event eventual ually ly addres addressed sed in the presen presentt method methodolo ology gy,, by applying applying an improved improved adaptive strain-mapped Gaussian sampling on a Green path integral, demonstrating fast execution
Closed-form
With derivatives as finite differences (NR) 2D Gauss or Green/Gauss Green/Gauss per – polygon Green/Gauss – Closed-form
With derivatives (NR)
Closed Closed-fo -form rm per polygo polygon n
With With deriva derivativ tives es (NR) (NR)
Clos Closed ed-f -for orm m per per trap trapez ezoi oid d
With With deri deriva vati tive vess (NR) (NR)
Closed Closed-fo -form rm per polygo polygon n
Nested Nested iterat iteration ionss
Closed-form
Derivative-fr ee ee (regula-falsi)
Green Green/G /Gau auss ss with with adap adapti tive ve strain-mapping
Derivative-free (Brent)
with customizable accuracy, which expresses the third prerequisite of the present method. A closely related issue to the aforementioned stress integration efficiency is whether section subdivisions are required in the employed employed integration integration scheme. There There are many studies (e.g. (e.g. [1,7,8,20])) where this geometric manipulation is imposed prior [1,7,8,20] to stress integration, naturally leading to a reduction in execution speed. For this reason, the suggested scheme will be formulated without need for any subdivisions (e.g. [6,18] [6,18]). ). The majority of solution strategies presented in the literature, for applying force equilibrium conditions, are based on secant or tangent schemes (e.g. Newton-Raphson), which require the prior calculation of derivative measures (e.g. section stiffness). However, the disadvantages of (a) the additional computational computational cost for calculat calculating ing derivati derivatives ves (e.g. (e.g. using using finite differences differences [21,22])) and (b) the inherent non-convergence issues related [21,22] to secant/tangent methods, has led some researchers to adopt simpler and more straightforward derivative-free derivative-free procedures [8,16,23],, which, under certain conditions, demonstrate stable [8,16,23] and fast performa performance nce.. Following Following the third third prerequi prerequisite site,, the present methodology will apply advanced derivative-free derivative-free methods [8] ods [8],, providing, where necessary, various techniques to guarantee convergence. convergence. In order to satisfy the aforementioned aforementioned three prerequisites of the suggested methodology, namely (a) unrestrained unrestrained section complexity (b) arbitrary material constitutive laws and (c) fastest possible execution with customizable accuracy, a new feature set is suggested in the last row of Table of Table 1. 1. In the subsequent chapters, the mathematical formulation of the methodology will be deployed, together with validation tests, case studies, extended comparisons with previous studies and benchmarks. The ultimate goal of the present study is to suggest a new robust and fast procedure for
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
the ultimate and moment–curvature analysis of arbitrary composite sections under biaxial bending and axial load, which may provide in turn, a solid basis for further developments in the fields of structural design, assessment and optimization [24]. 2. Mathematical formulation 2.1. Geometry and material definitions
Following the first prerequisite of the present methodology, the section under consideration can be consisted of an unlimited number of individual section components, namely (a) surfaces (S 1, S 2, . . ., S n), describing surface-distributed materials such as concrete, masonry or structural steel, (b) multi-segment lines (L1, L2, . . ., Ln), usually simulating linearly distributed reinforcement or fiber-reinforced polymer (FRP) strips and (c) fibergroups (FG1, FG2, . . ., FGn), i.e. dimensionless fibers representing individual
35
reinforcement bars or tendons (Fig. 1). The above three component types can contribute either positively or negatively to the section response (i.e. the corresponding internal actions will be added or subtracted respectively), allowing the efficient embedment of voids or the use of multi-nested materials. This feature not only provides an explicit description of section openings, without resorting to complicated ‘ fictitious cuts’ [3,4,17,28], but also the means to avoid double-contribution errors, i.e. by subtracting the concrete material under structural steel or reinforcement fibers [17], similarly to the ‘ background material’ concept [8]. Surfaces are defined as a counterclockwise series of vertex coordinates ( x j, y j), together with the included angle value of the succeeding circular arc segment x j (ignored for straight segments). It is also noted that a negative angle value leads to the subtraction of the corresponding circular sector from the surface. Therefore:
S i ! ½ð x1 ; y1 ; x1 Þ; ð x2 ; y2 ; x2 Þ . . . ð xn ; yn ; xn Þ
ð1Þ
Lines are defined as a series of vertex coordinates (multi-segment), including the contribution area of the succeeding segment ( A j):
Li ! ½ð x1 ; y1 ; A1 Þ; ð x2 ; y2 ; A2 Þ . . . ð xn1 ; yn1 ; An1 Þ; ð xn ; yn Þ
ð2Þ
Similarly, fibergroups are defined as a set of fiber coordinates and corresponding contribution areas:
FGi ! ½ð x1 ; y1 ; A1 Þ; ð x2 ; y2 ; A2 Þ . . . ð xn ; yn ; An Þ
ð3Þ
Fig. 2 depicts the definition of the three aforementioned section component types. After the definition of individual section components, an adaptive linearization process is performed on each surface component (S i) which contains circular arc segments, transforming them to closed polygons, in order to maintain compatibility with the Green path integral transformation concept [17,28], discussed later. This algorithm is an improved version of the constant divisions approach by Kwan and Liauw [28], outlined as follows: all arcs of equal length are bisected in consecutive steps (sweeps), until the surface area of the running sweep differs from the previous one, less than a requested tolerance:
jðs Atot s1 Atot Þ=s1 Atot j < A TOL
Fig. 1. Composition of arbitrary composite section.
ð4Þ
where A tot is the total surface area, s is running sweep and A TOL is the tolerance value. Suggested values of A TOL may be in the range of 1–2%. In the most unfavorable case of a single circle, this leads to a 64-segment polygon for A TOL = 1% (default), with an area accuracy of 0.16%. For A TOL = 2%, it leads to a 32-segment polygon, with
Fig. 2. Definition of individual section components.
36
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
(e.g. exponential [12]), which will be handled by adaptive sampling. Moreover, three distinct strain threshold values should be declared for each material, which will be used in a later stage for the direct definition of the ultimate strain profile: e cu is the strain corresponding to the ultimate compressive strength. e co is the strain corresponding to the compressive strength under pure compression. e tu is the strain corresponding to the ultimate tensile strength.
Fig. 3. Adaptive linearization procedure.
an area accuracy of 0.64%. In addition, it has been reported in Charalampakis and Koumousis [8] that ultimate response discrepancies between exact circle integration and equivalent polygonal representation are negligible for more than 48 vertices. Fig. 3 shows a worked example of the adaptive linearization process. The advantage of this improved linearization method is that smaller arc segments (e.g. rounded edges of steel sections) introduce fewer additional vertices to the corresponding section surface, thus improving the efficiency of the integration procedure. The second prerequisite of the present method allows fully arbitrary constitutive laws for the assigned materials. For this reason, the material definition is consisted of a series of arbitrary stress– strain (r e) function parts ( f 1, f 2, . . ., f n) in a piecewise form (Fig. 4a). These function parts may preferably be polynomials, for exact Gaussian integration; however they can be of any other form,
Examples of specific material definitions are shown in Fig. 4b and 4c for concrete (parabolic-linear) and steel (bilinear) respectively, according to EN1992 [9]. It is possible that any of these criteria can be ignored when defining the ultimate profile (e.g. for concrete tensile strength). Throughout the present formulation, various geometrical and loading parameters of the considered section are defined ( Fig. 5). Resultant internal actions ( N , M X , M Y ) are calculated with respect to the reference coordinate system XCY, where C is the section origin, that can be either the geometrical, elastic or plastic center [29]. It has been previously reported that the plastic center yields better convergence of the solution process, especially for axial loads approaching the compressive capacity (because it is contained in the isoload contour [16]). Thus, it will be adopted throughout this study, and will be further improved. The exact location of the neutral axis is determined by its orientation ( h) and depth ( d), while the current strain profile by the curvature ( u) and the strain-atorigin (eo = ud) [8]. The neutral axis depth and the strain profile are always defined on the rotated (by angle h ) coordinate system xCy, while a positive curvature is assumed to induce compression at the top section fiber. The strain profile is considered linear, following the Bernoulli–Euler assumption, i.e. plane sections remain plane and perpendicular to the axis of the element. Finally, moment vectors M X and M Y are pointing to the positive direction of X and Y axes respectively, and the inclination of their resultant vector is denoted by angle a . It is pointed out that the orientation angle of the neutral axis ( h) and the inclination angle of the resultant
Fig. 4. (a) Generalized material definition, applications on (b) concrete and (c) steel.
37
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
to the rotated coordinate system xCy (see Fig. 5). Since the section is composed of a specific number of components, i.e. (nS) surfaces, (nL) lines and (nFG) fibergroups, the resultant actions are the summation of individual component contributions as follows: nS
R ¼
nL
X
½signðS i ÞRS ;i þ
i¼1
nFG
X
½signðLi ÞRL;i þ
i¼1
X
½signðFGi ÞRFG;i
ð5Þ
i¼1
where R is the resultant action (N , M x, M y, compression negative), n is the component count, indices S , L and FG correspond to surface, line and fibergroup components respectively and sign() denotes whether the component contribution is positive (actual material) or negative (void). Following the integration process, moments M x and M y are transformed back to the reference system XCY as follows:
M X ¼ M x cosðhÞ þ M y sinðhÞ M Y ¼ M x sinðhÞ þ M y cosðhÞ
ð6Þ
The integration algorithms for the individual component types (S i, Li, FG i) are described below. Fig. 5. Definition of section loading parameters.
moment (a) are generally not coincident ( h – a), except for very special cases. It will be shown later that this discrepancy is the root of various divergence problems, which will be treated accordingly. The geometry and material definitions described in this chapter are efficiently structured in computer code, using a pure object-oriented approach (Fig. 6). Initially, two basic classes are introduced: (a) Material, including material properties and respective stress calculation methods and (b) Points, including a set of coordinates ( x j, y j) and geometry manipulation methods (e.g. rotate, pan etc.). Subsequently, three descendant classes from ancestor class Points are formed, one for each section component type, namely Surface, Line and FiberGroup. Each class includes additional properties (e.g. included angles or contributing areas), an instance of the Material class (object) as the assigned material and a method to perform stress integration. Moreover, the Surface class includes the method to perform adaptive linearization. The fundamental Section class contains an unlimited number of section component objects (Surface, Line, FiberGroup), together with methods to define centroids, determine the ultimate limit state and summate the internal actions from the individually integrated component contributions. Finally, a descendant of Section class, named SectionSolver , is adding the necessary methods to apply force equilibrium conditions and calculate the incremental and ultimate section response (solution strategies). The advantage of this object-oriented approach (as opposed to traditional procedural programming) is that any section under consideration can be simply created as a SectionSolver object, encapsulating all the necessary properties and numerical procedures (methods) to perform section analysis (analysis engine). These procedures will be described in detail in the subsequent paragraphs.
2.2.1. Surface stress integration For the stress integration of individual surfaces ( S i) it has been already justified that, in order to handle arbitrary stress–strain constitutive laws, the most appropriate solution is to apply a numerical integration scheme. In this study, the Gauss–Legendre quadrature method, originally introduced in Fafitis [17], will be adopted. Initially, the resultant actions ( RS ,i) are represented by area integrals of general form:
RS ;i ¼
ZZ
xr ys rð yÞdxdy
ð7Þ
S i
where S i represents the surface under consideration, ( x, y) are the surface planar coordinates, r ( y) is the stress function expressed in terms of the coordinate y, and (r , s) are exponents depending on the requested resultant action, i.e. (r , s) = (0,0) for axial N S ,i, (0,1) for moment M xS ,i and (1,0) for moment M yS ,i. By applying Green’s theorem, the area integral is transformed into a line integral along the closed boundary L i that encloses the area S i as follows:
1 RS ;i ¼ r þ 1 ¼
1 r þ 1
I X
x
1 y rð yÞdy ¼ r þ 1
r þ1 s
Li
n ‘i
n‘i
X "I
r þ1 s
x y rð yÞdy
j¼1
‘ j
#
I j
where the surface boundary Li is represented by a closed polygon of n ‘i segments. Therefore, the integration problem breaks down to the integration of each individual polygon segment ‘ j , denoted as the elementary integral I j. The equation of this line segment may be expressed as:
2.2. Stress integration
Following the geometry and material definitions, perhaps the most demanding numerical task that greatly influences both the accuracy and speed of the analysis, is the employed stress integration scheme. The fundamental integration problem is defined as follows: calculate the resultant internal actions of the composite section (axial N and moments M X , M Y ) for a prescribed neutral axis orientation (h) and strain profile ( eo, u), denoted as I(eo, u, h). The section origin (C ) is initially set to the plastic center [16,29] and the reference coordinate system XCY is rotated by angle h, resulting
ð8Þ
j¼1
Fig. 6. Object-oriented formulation of the analysis engine.
38
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
‘ j ! x ¼ a j y þ b j
ð9Þ
where a j and b j is the slope and intercept of line ‘ j , respectively. Introducing (9) in (8), the elementary integral is further simplified into a line integral of a single-variable nonlinear function F j( y):
I j ¼
I
r þ1 s
ða j y þ b j Þ y rð yÞdy ¼
‘ j
I
F j ð yÞdy
ð10Þ
‘ j
At this stage, the Gauss–Legendre integration scheme for the elementary integral I j is applied. This scheme is theoretically yielding exact solutions for polynomial functions F j( y) of order n, when a minimum number of nG ¼ n n 2 þ 1 Gauss points are utilized (the symbol ‘n’ denotes integer division). In the original study of Fafitis [17], a constant number of Gauss points nG = 3 was selected, claiming exact results for Code provided stress–strain laws for concrete, e.g. stress-block or parabolic-linear [9,30]. Nonetheless, it has been observed that, even for these piecewise-polynomial functions (e.g. parabolic-linear), the originally suggested application of Gaussian quadrature per segment ð‘ j Þ, may lead to non-exact results, because it is highly probable that multiple function parts will be mapped on the same polygon segment, exhibiting non-continuity in the derivative and virtually resulting to a non-polynomial form. This observation is clearly depicted in Fig. 7, for the case presented in the original study [17] (circular points). It is therefore suggested herein, that the sampling process should be applied per function part within each segment (rectangular points), in order to guarantee exact solutions for the widely used piecewise-polynomial stress– strain relationships. Table 2 shows a comparison between the original and the improved integration scheme (described in detail hereinafter) for the example shown in Fig. 7, where it is demonstrated that the latter yields absolutely exact results for all considered actions.
A schematic view of the suggested strain-mapped Gaussian sampling procedure is shown in Fig. 8. Initially, for each function part f k of the constitutive law (Fig. 4), the values of the (predefined) strain limits elim A,k and elim B,k are mapped on the segment ‘ j under consideration by a simple strain-to-coordinate transformation, which leads to the corresponding coordinate limits y lim A,k and y limB,k:
yðeÞ ¼ ðe eo Þ=u
ð11Þ
where u and eo are the current strain profile parameters. If any coordinate limit falls outside a segment limit ( y A, j or yB, j), it is replaced by the segment limit itself. If both coordinate limits fall outside the same segment limit, the respective function part f k is not contributing to the integral and is subsequently ignored. Finally, Gaussian sampling is applied for all contributing function parts of the considered segment ‘ j , resulting to the elementary integral I j, following (10):
I j ¼
I
: F j ð yÞdy¼
‘ j
X" nf
k¼1
nG
k 1 ð ylim B;k ylim A;k Þ ½wm F j ð ym Þ 2 m¼1
X
#
ð12Þ
where nf is the number of stress–strain function parts, nG k is the number of utilized Gauss points i.e. the quadrature order, ym is the coordinate where function F j is sampled and wm is the corresponding weight. The coordinate ym for the Gauss the point number m is derived from the corresponding Gauss index km = [1. . .1] as follows:
ym ¼
1 km ð y þ ylim A;k Þ þ ð ylim B;k ylim A;k Þ 2 lim B;k 2
ð13Þ
Finally, the function F j is sampled at coordinate ym, using the function part f k, in stress–strain terms, as follows:
F j ð ym Þ ¼ ða j ym þ b j Þr þ1 ysm f k ðeo u ym Þ
Fig. 7. Improvement of the Gaussian integration scheme suggested by Fafitis [17].
Table 2
Comparison between Fafitis [17] and present integration scheme. h = 45, d = 9.972 in
N (kip)
M x (kip-in)
M y (kip-in)
Exact (100 Gauss points) Fafitis (reported values) Fafitis (simulated) Fafitis – error% Present strain-mapped scheme
1999.66 1997.29 1997.14 0.13% 1999.66
7405.90 7411.72 7412.03 0.08% 7405.90
2820.60 2796.55 2797.00 0.84% 2820.60
ð14Þ
39
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
Fig. 8. Adaptive strain-mapped Gaussian integration.
An important advantage of the suggested strain-mapped technique is that each function part ( f k) can acquire a different number of Gauss points (nGk) which, if set equal to ðn þ r þ 1 þ sÞ n 2 þ 1 for polynomial parts of order n, produces exact integration results with an optimum reduction in computational cost, as compared to a constant quadrature order [1,17,19]. Another advantage is that no geometric manipulations, such as section subdivisions [1,4,5,7,8,20,26,27], are imposed prior to stress integration, which obviously enhances the overall analysis performance. More specifically, it is assumed that the present strain-mapped approach yields similar results to the MTLI method suggested by Bonet et al. [1] (where the section is decomposed into thick layers corresponding to each function part and integrated separately), yet without the laborious section decomposition and the restriction of a constant Gauss quadrature order. There are certain cases where a function part ( f k) is of non-polynomial form (e.g. ascending branch in high-strength concrete stress–strain relationship, according to EN1992 [9]) and therefore, Gaussian integration can yield only approximate estimates. For this reason, an adaptive integration technique is introduced herein, which is described as follows: each function part ( f k) data structure is upgraded with memory, which stores the largest strain range integrated for that part Demax,k, throughout the analysis process (initially set to zero): Demax;k ¼ Maxðjuð ylim B;k ylim A;k ÞjÞ
ð15Þ
nGk þ1 m¼1
nGk
½wm F j ð ym Þ
X
nGk
!,X
½wm F j ð ym Þ
m¼1
RL;i ¼
I
X " I
!
½wm F j ð ym Þ
m¼1
# X
n‘i
r s
x y t ð yÞrð yÞdy ¼
Li
n ‘i
r s
x y rð yÞdy ¼
t j
‘ j
j¼1
I j
ð17Þ
j¼1
where the multi-segment line Li is represented by a series of n‘i consecutive segments of thickness t j. Therefore, the integration problem breaks down again to the integration of the elementary integral I j, using Gaussian quadrature (Fig. 8), similarly to Eqs. (10) and (12):
I j ¼ t j
I I
r s
ða j y þ b j Þ y rð yÞdy ¼ t j
‘ j
During the running integration, if the strain range surpasses the respective stored value, an adaptive process is triggered, which sequentially increases the utilized number of Gauss points, until a prescribed integration accuracy is reached:
X
2.2.2. Line integration The possibility to include and directly integrate linearly distributed materials in a composite section is introduced for the first time in the present study. This feature is practical for simulating closely spaced reinforcement (without resorting to a dense array of bar fibers), linear reinforcement distributions of unknown bar spacing and FRP strips [31]. It can also reduce modeling time and potential errors, as the exact choice of rebar diameter and spacing usually emanates from Code regulations rather than strength requirements. The adaptive strain-mapped integration method described in the previous paragraph can be directly adapted to line components (Li), following a few necessary modifications. The resultant actions ( RLi), contributing to the total section response (Eq. (5)) are calculated using the following line integral, similarly to (8):
I j ¼ t j
‘ j
I
F j ð yÞdy
ð18Þ
‘ j
: F j ð yÞdy¼
nGk
X" X nf
k¼1
Ak 2
½wm F j ð ym Þ
m¼1
#
ð19Þ
where A k is the contributing area of the function part f k of segment ‘ j :
< I TOL ð16Þ
where I TOL is a tolerance value (default 5‰). The main advantage of the suggested adaptive integration technique is that is capable to automatically handle fully arbitrary stress–strain functions with customized accuracy and optimized execution speed, following the third prerequisite of the present methodology.
ylim B;k ylim A;k Ak ¼ A j yB; j y A; j
ð20Þ
and A j is the area of segment ‘ j , as defined in (2). All other formulations are identical to the surface integration formulation. However, contrary to surface integration, horizontal segments are now contributing to the total integral and they are specially treated as follows:
40
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
N j ¼ A j rð y A; j Þ M x; j ¼ A j y A; j rð y A; j Þ x A; j þ xB; j M y; j ¼ A j rð y A; j Þ 2
ð21Þ
For the exact integration of polynomial functions, the direct line integration needs even fewer Gauss points than surface integration, specifically ðn þ r þ sÞ n 2 þ 1 for polynomial parts of order n. This renders this approach superior in terms of execution speed than the equivalent fiber discretization of linearly distributed materials, which may require notably more arithmetic operations to reach an equivalent level of accuracy. 2.2.3. Fibergroup integration The integration procedure for each fibergroup component ( FGi) is straightforward and performed analytically as follows: nFGi
RFG;i ¼
X
½ A j x jr y js rðeo u y j Þ
ð22Þ
j¼1
where nFGi is the number of group fibers, A j is the fiber contributing area and ( x j, y j) are the fiber coordinates, as defined in (3). 2.3. Multicriteria ultimate limit states
ecu;i =ð ymax;i dÞ if d < y max;i if d < y min uu;i ¼ u min ¼ Min eco;i =ð yo;i dÞ etu;i =ð ymin;i dÞ if d > y min;i
8>< >:
yo;i ¼ y min þ ð ymax;i ymin Þ
Following the implementation of the stress integration scheme that calculates the resultant internal forces for a given neutral axis orientation (h) and strain profile ( eo, u), the next step is to establish a method to derive the ultimate strain profile ( uu, eou) for a given natural axis depth ( d) and subsequently integrate the ultimate strain profile, as a function of the natural axis depth and orientation: I(eou, uu, h) I(d, h). The resultant actions emanating from this strain profile define the locus of points in the triaxial stress space that correspond to the conventional (or nominal) ultimate strength of the section [8]. It has to be noted though, that the conventional ultimate strength may slightly differ in certain cases from the physical one, which corresponds to the top of the respective moment–curvature diagram [15]. However, in the present study, the former definition will be employed, since it is compatible with Code requirements (e.g. [9,30]) and its derivation is deterministic, leading to fast and accurate solutions, as opposed to the laborious iterative calculation of moment–curvature maxima [8,15]. This derivation will be the basis for the construction of ulti?
mate strength interaction curves and 3D failure surfaces, using special solution strategies presented in the subsequent chapter. Since the section consists of an unlimited number of components and different materials, multiple failure criteria should be taken into account in order to determine the global ultimate strain profile that corresponds to the conventional ultimate strength of the section. This is a necessary improvement of the existing literature, where either a single concrete [16] or a dual concrete/steel [18] strain criterion is taken into account for composite sections. Moreover, the pure compression criterion imposed by EN1992 §6.1(5) [9] and other Codes (referred hereinafter as ‘ C-restriction’, after the pertinent figure 6.1 of EN1992) is generally neglected, except in Rosati et al. [5]. According to the present approach, for a given natural axis depth (d), the ultimate strain profile ( uu,i, e ou,i) is initially defined for each section component ( S i, Li, FGi) (always on the rotated coordinate system xCy), using the three material strain threshold values (ultimate strength criteria) already declared in the material definition (ecu, eco, etu), as follows (Fig. 9):
eou;i ¼ uu;i d
ecu;i eco;i
ð23Þ
ð24Þ ð25Þ
where ymax,i and ymin,i are the component coordinate limits and ymin is the whole section coordinate minimum. The first row in Eq. (23) corresponds to the ultimate compressive strength, the second to the ultimate compressive strength under pure compression (C-restriction) and the third to the ultimate tensile strength. The derived value of uu,i essentially corresponds to the curvature value attained on the first triggering of any of the aforementioned failure criteria. It is also observed in Fig. 9, that the C-restriction criterion is applicable only when the neutral axis falls below the section (d < y min ), i.e. the whole section is under compression. It is pointed out that the Crestriction as described in EN1992 [9] refers to the contribution of concrete in a reinforced concrete section, i.e. an effectively composite section. It is deemed that this restriction should still hold for concrete that participates in fully arbitrary composite sections. Either way, the application of the C-restriction is fully optional
Fig. 9. Ultimate strain profile for individual section components.
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
and can be easily disabled. However, it has been observed (see next chapters) that it leads to stable solutions even for the uncommon cases of high axial compression. Following the calculation of all component ultimate curvature values (uu,i), the resulting ultimate curvature for the whole section is simply derived as the globally minimum value, as follows:
uu ¼ Minðuu;i Þ 8 ðS i ; Li ; FGi Þ eou ¼ u u d
ð26Þ ð27Þ
The above definition of the global ultimate curvature value implies that the section is considered to have reached its ultimate limit state when any of its components first reaches its own ultimate state, as defined by its corresponding material criteria (Eq. (23)). It is also noted that each section component may be optionally excluded from the process of ultimate limit state definition, upon engineering judgment. The advantage of this approach is that it can automatically handle sections of any complexity. Figs. 10 and 11 show schematically the derivation of the multicriteria ultimate strain profile for composite sections under prevailing flexure and pure compression (where C-restriction is also triggered), respectively.
41
2.4. Solution strategies
The last component of the present methodology is the implementation of robust and rapid solution strategies for applying equilibrium conditions between external and internal actions and subsequently derive the requisite moment–curvature curves, ultimate strength interaction curves and 3D failure surfaces. Following the prerequisite of speed optimization, derivative-free procedures will be adopted herein, avoiding the additional computational cost of derivative calculation and any inherent non-convergence issues related to secant/tangent methods.
2.4.1. Interaction curves and failure surfaces The ultimate failure surface is defined as the locus of resultant action points (N , M X , M Y ) in stress space, calculated from the stress integration procedure applied on the ultimate strain profile, for a prescribed position of the neutral axis: I (d, h). Furthermore, planar interaction curves are usually expressed either as (a) M X M Y contours for a constant external axial force N or (b) M N contours for a constant inclination of the resultant moment vector a .
Fig. 10. Ultimate strain profile of composite section under prevailing flexure.
Fig. 11. Ultimate strain profile of composite section under pure compression.
42
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
For the construction of the 3D failure surface, there are three different approaches of increasing computational cost, namely: The d h method where both neutral axis depth ( d) and orientation (h) values are incremented, and resultant points are derived from direct integration of the corresponding ultimate strain profile. The N h method where the natural axis depth d is iteratively calculated for incrementing values of axial load ( N ) and natural axis orientation ( h), producing the so-called ‘ isogonic ’ meridians [8,20]. The N a method, where both the natural axis depth ( d) and orientation ( h) are iteratively calculated for incrementing values of both axial load ( N ) and moment vector inclination ( a) [8,20].
Fig. 12 depicts the aforementioned methods, where it is observed that the d h approach produces an irregular wireframe with non-planar, non-equidistant meridians, as opposed to the more elaborate N h and N a methods. It is apparent, that interaction curves (surface cuts) are produced in the same sense, specifically using the N h (or N a) approach for M X M Y contours (for a constant N ) and the N a approach for N M contours (for a constant a ). As a result, there are two problems to be addressed:
given N ; h ! solution ! d given N ; a ! solution ! d; h
ð28Þ ð29Þ
For the solution of (28) and (29), the root-finding Brent method [32] will be employed, successfully introduced for section analysis in Charalampakis and Koumousis [8]. The Brent method is considered superior than similar derivative-free Regula-Falsi [16,23] or bisection methods, providing guaranteed and fast convergence with fewer iterations, for well bracketed roots. The Brent root-finding method may be symbolically represented as:
vt ¼ B vv AB ½UðvÞ ¼ w t
ð30Þ
where v t is the solution (root), sought in bracket [v A. . .vB], which yields the requested (target) result w t , using the evaluation procedure U(v). It is noted that in the current implementation, the target
value w t may be nonzero and the evaluation procedure U is noninvertible (i.e. cannot be solved with respect to vt ). The convergence criteria for the Brent method are the requested accuracies vTOL and wTOL for v and w variables respectively. By applying definition (30) in problem (28), the solution is expressed as:
dt ¼ B dd AB ½Iðd; hÞ ¼ N t
ð31Þ
where d t is the natural axis depth solution that lies in the [d A. . .dB] bracket and N t is the given (target) axial force. Similarly, the solution to problem (29) is defined as: h
d
ht ¼ B h AB ½Iðd ¼ B dB ½Iðd; hÞ ¼ N ; hÞ ¼ a t
ð32Þ
a ¼ arctanðM Y =M X Þ
ð33Þ
A
where ht is the natural axis orientation solution that lies in the [h A. . .hB] bracket and at is the target moment vector inclination. This solution involves a double-nested application of the Brent method; the axial force (N ) is equilibrated in the inner loop and the inclination of the moment vector (a) is equilibrated in the outer. A crucial element that determines the performance of the aforementioned solution strategy is the effective definition of the rootfinding brackets [ v1. . .v2]. In the trivial case where extended brackets are selected e.g. ½d A . . . dB ¼ R or ½h A . . . hB ¼ ½0 . . . 2p, apart from the apparent computational overhead (i.e. more iterations are needed to reach convergence), it is possible that multiple or no roots exist, resulting to undesirable non-convergence issues [23]. For this reason, a new adaptive bracketing technique is introduced herein, depicted in Fig. 13 and described as follows: starting from an initial value v0, the corresponding w 0 = U(v0) value is first calculated. From the initial w0 and target wt values, the search direction is defined as: positive (+ v) if w t > w0 and negative ( v) otherwise. The value of v is then incremented towards this direction, and corresponding w values are calculated in the form of wi = U(vi), where i is the increment number. The search procedure stops when the target w t value lies between w i 1 and w i and the resulting search bracket for the Brent method is [ v A. . .vB] = [vi1. . .vi].
Fig. 12. Various 3D interaction surface construction methods of increasing computational cost.
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
43
Fig. 13. Preliminary bracket-search procedure before applying derivative-free solution.
The advantage of this preliminary process is that, by introducing a very few additional section integrations, the target N t and at values are always bracketed (in Eqs. (31) and (32)), and hence convergence of the solution is guaranteed with reduced computational cost, compared to the aforementioned trivial bracketing practice. Further details on this technique are the following: it has been measured that, for di values, an exponential incrementing scheme is providing the fastest performance, while linear incrementing is suitable for h i values, specifically:
di ¼ d0 þ ð ymax ymin Þebþi hi ¼ h 0 þ iðp=16Þ
ð34Þ ð35Þ
where b ¼ 2 and d 0 = 0 for solving Eq. (31) and b ¼ 3, h 0 = at and d 0 is equal to the previously converged value of d (inner loop), for solving Eq. (32). Finally, the tolerance values used in the Brent method are:
8>< >:
N TOL ¼ 10 4 ðN max N min Þ dTOL ¼ 10 4 ð ymax ymin Þ hTOL ¼ a TOL ¼ 10
2
ð36Þ
deg
where N max , N min are the section tensile and compressive capacities and y max , y min are the section coordinate limits (for the current natural axis orientation h ). Nonetheless, in the process of constructing the 3D failure surface using the most demanding N a method, there were observed certain instability issues for unsymmetrical sections, which can be explained as follows (Fig. 14): although it has been reported in previous studies [16] that the employment of the plastic center ( C )
provides guaranteed convergence, because the origin of the M X M Y axes always stays within the respective interaction curve (isoload contour), it is observed that, especially for high tensile axial loads, this consideration does not hold. This is because any materials without tensile strength (e.g. concrete) abruptly stop contributing under tensile loads and hence the ‘updated’ plastic center of the still active materials in tension (e.g. steel) deviates significantly from the one originally considered. As a result, the origin of the M X M Y axes actually falls outside the interaction curve, thus there are certain moment vector inclination angles ( a) which are not available. Hence, the requested solution cannot be possibly bracketed (Eq. (32)). For this reason, a new self -correcting centroid is herein suggested, in order to establish convergence, as follows. It has been observed that failure surfaces for symmetrical sections are surfaces of revolution about a vertical axis connecting the apexes of ultimate axial capacity ( N min , N max ). This vertical axis coincides with the locus of M X M Y origins (i.e. a vertical line passing from 0, 0) and these origins are always contained in the respective horizontal contours (for constant N ) . On the contrary, for unsymmetrical sections, the failure shape is skewed with an inclined axis connecting these apexes. It is therefore apparent that there may be regions where the M X M Y origins, still lying on a vertical line, fall outside the respective isoload contours and hence lead to non-convergence. The concept of the self-correcting centroid is to define proper eccentricities for the plastic center ( C ) that, when imposed on the section geometry, will move the respective M X M Y origin inside the isoload contour. This technique is depicted in Fig. 15 and is described as follows: initially, the moment skewness of the failure surface (i.e. the slope of the inclined axis) is defined for each direction as:
Fig. 14. Convergence issues in constructing N a interaction surfaces for unsymmetrical sections.
44
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
Fig. 15. Definition of interaction surface skewness in unsymmetrical sections.
M X ;max M X ;min N max N min M Y ;max M Y ;min M Y ;skew ¼ N max N min M X ;skew ¼
ð37Þ ð38Þ
where (M X ,max, M Y ,max) and (M X ,min, M Ym in ) are the nonzero corresponding moment values for N max and N min (at surface apexes). The plastic center eccentricities ( e X , eY ) for a given axial force N are calculated as:
M Y ;min þ M Y ;skew ðN N min Þ N M X ;min þ M X ;skew ðN N min Þ eY ðN Þ ¼ N e X ðN Þ ¼
ð39Þ ð40Þ
Having calculated the above eccentricities, the section is accordingly shifted before applying the solution procedure (Eq. (32)). The outcome of this procedure is that any potential solution instabilities during the construction of the N a failure surface (or interaction curves) for unsymmetrical sections are eliminated, as shown in the right of Fig. 14. 2.4.2. moment–curvature curves The derivation of uniaxial moment–curvature curves for a given axial force in the context of the present methodology is straightforward; initially, the uniaxial direction i.e. the neutral axis orientation h is selected and a series of section integrations for linearly
incrementing values of curvature u are performed as follows: for each curvature step u, the given external axial load is equilibrated to the respective internal resultant N , by iterating the value of strain-at-origin eo [8]. As a result, the corresponding moment (MX or M Y ) is plotted vs. the curvature value. This iterative scheme is performed exactly in the same way as the ultimate strength solution, using the Brent method as follows:
eot ¼ B eeoBoA ½Iðeo ; u; hÞ ¼ N t
ð41Þ
where e ot is the strain-at-origin solution that lies in the ½ eoA . . . eoB bracket and N t is the target axial force. For adaptively deriving the application bracket, an exponential incrementing scheme was applied (Fig. 13):
eoi ¼ eo0 þ j eo0 je3þi
ð42Þ
where e o0 was set equal to its previously converged value and the tolerance value eoTOL was set equal to 107. It is noted that the above process is not stopped when a predefined failure criterion is reached (e.g. ultimate strain), but is rather terminated when section collapse occurs (zero calculated moments). This provides a complete picture of the section response (including post-peak response), where any failure criteria may be marked a posteriori (from the stress–strain history).
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
2.4.3. Synopsis The suggested methodology for arbitrary section analysis, consisting of geometry and material definitions, the integration scheme, the derivation of ultimate limit states and solution strategies, all bound under a object-oriented implementation, has emerged the following points of novelty:
(a) Adaptive linearization. (b) Adaptive strain-mapped integration. (c) Integration of linearly distributed materials. (d) Derivation of multicriteria limit states. (e) Adaptive bracketing for derivative-free methods. (f) Self-correcting centroid for unsymmetrical sections. In the following chapters, the performance of the suggested methodology will be demonstrated through various case studies, benchmarks and comparisons with previous studies.
3. Validation and case studies
The object-oriented analysis engine described in the previous chapter was integrated into a user-friendly graphical front-end software (Fig. 16), which contains all the necessary features to conduct the forthcoming validation tests, comparisons and benchmarks. Before proceeding, it is herein introduced a testbed consisting of five R/C and composite sections of increasing complexity, for future reference. Section (A) is a simple R/C symmetrical column section, (B) is similar to (A), but with an unsymmetrical placement of reinforcement, (C) is a bilaterally jacketed R/C section with different concrete and steel materials for the core and the jacket, (D) is the well-known composite section example reported in Chen et al. [16] and finally, (E) is a conceptual, complex bridge deck composite section consisting of two different concrete grades, openings, fiber and linearly distributed reinforcement, structural steel sections of two different grades and FRP strips. Detailed geometry and material definitions are provided in Fig. 17.
45
The first validation test investigates the applicability of the adaptive linearization procedure, as it is applied on the plain concrete curved section of Fig. 3. The analysis results are presented in Fig. 18 as a series of moment interaction curves for an increasing number of linearization sweeps. It is observed that section response rapidly converges to a stable solution, for an adaptive linearization process of 1% tolerance, demonstrating the effectiveness of the method. The second test investigates the applicability of the adaptive strain-mapped integration scheme, especially on non-polynomial functions. In Fig. 19, section (A) is analyzed using different concrete stress–strain relationships, namely the parabolic-linear and nonlinear models of EN1992 [9] and the Popovics model [12], with only the former holding a piecewise polynomial form. It is observed that continuous and smooth moment interaction curves are derived for all cases, also demonstrating the expected non-convexity for softening laws (as reported in [18]) for high levels of axial compression ð m ¼ N = Atot ¼ 0:80Þ. Therefore, the robustness of the integration scheme and the accompanying solution algorithm is confirmed. The third test investigates the validity of line integration, by comparing a linearly distributed (perimetric) reinforcement model to three equivalent fiber discretizations of increasing density. From the moment interaction curves depicted in Fig. 20, it is observed that fiber modeling rapidly converges to the linear distribution model for higher rebar density, which verifies the applicability of line integration. Note that using line distributed reinforcement is slightly conservative. In Fig. 21, moment interaction curves for composite sections (C) and (E) under varying axial loads are presented. For the same sections, full moment–curvature curves (up to collapse) for various natural axis orientations ( h) are depicted in Fig. 22. From the above analyses and especially for section (E), it is concluded that the present methodology can successfully handle exceptionally complex composite sections without any instability issues. However, the actual accuracy of the analysis results will be investigated by various comparisons to the existing literature, presented in the subsequent chapter.
Fig. 16. Graphical front-end for the analysis engine.
46
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
Fig. 17. Testbed including five R/C and composite sections of increasing complexity.
Fig. 18. Section response during adaptive linearization.
47
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
150
M Y (kNm)
N = -3000 kN ( ν = -0.80)
EC2 Parabolic-Linear
-25
EC2 Nonlinear (Sargin)
-20
Popovics
-15
EC2 Parabolic -Linear σ (Μ Pa)
-10
100
-5
ε
0 0.0000
-0.0005
-0.0010
-0.0015
-0.0020
-0.0025 -0.0030
50
EC2 Nonli near
-25 -20
-0.0035
σ (Μ Pa)
-15 -10
0 -200
-150
-100
-50
0
50
100
150
200
-5
ε
0 0.0000
-0.0005
-0.0010
-0.0015
-0.0020
-0.0025
-50
-20
-0.0035
Popovics
-25
M X (kNm)
-0.0030
σ (Μ Pa)
-15 -100
-10 -5
Bilinear steel
0
f y = 500 MPa
0.0000
ε -0.0005
-0.0010
-0.0015
-0.0020
-0.0025
-0.0030
-0.0035
-150
Fig. 19. Section response using various concrete constitutive models.
400
4 bars
N = -1000 kN
M Y (kNm)
8 bars 16 bars
300
Linear
200
100
0 -1000
-800
-600
-400
-200
0
200
400
600
800
1000
-100 50
M X (kNm) -200
-900
-650
-300
Section 30×70 cm, cov er 5 cm f c = 20 MPa, EC2 parabo lic-linear 2
f y = 500 MPa, bilinear, A tot = 40 cm -50
-400
Fig. 20. Section response using fiber and linear reinforcement modeling.
4. Comparisons with previous studies
The first comparison refers to the analysis of steel equal angle sections, reported in Charalampakis [33], which is considered a
computationally demanding task due to the section nonconvex and curved geometry. In Fig. 23, moment interaction curves for a steel section L40 40 4 rotated by 45 , for zero and high tensile and compressive axial loads ( m = 0, m = ±0.8) are shown. It is ob-
48
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54 N = 500 kN
320
N = 0 kN N = -1000 kN
M Y (kNm)
Section C
N = -2000 kN N = -3000 kN
240
with respect to plastic center
N = -4000 kN
N = 50000 kN N = 25000 kN N = 0 N = -50000 kN N = -100000 kN N = -150000 kN
100000
M Y (kNm)
Section E
80000
with respect to plastic center 60000
160 40000
80 20000
0 -320
-240
-160
-80
0 0
80
160
240
320
-60000
-50000
-40000
-30000
-20000
-10000
0
10000
20000
30000
40000
-20000 -80
M X (kNm)
M X (kNm)
-40000
-160 -60000
-240
Nmax = 626.3 kN
Nmax = 56419 kN
Nmin = -4232.9 kN
-80000
Nmin = -156348 kN -320
-100000
Fig. 21. Moment interaction curves for sections C and E under various axial loads.
Fig. 22. Moment–curvature response of sections C and E for various natural axis orientations.
served that the results from both studies are absolutely coincident, which confirms the accuracy of the present method for single-material sections of curved shape. Moreover, perfect symmetry in numerical results is observed for equal axial loading of opposite signs, which demonstrates the reliability of the solution. The next comparison refers to reinforced concrete sections (i.e. two materials) of irregular geometry, reported in Rosati et al. [5]. In Fig. 24, moment interaction curves for a G-shaped (left) and multicell (right) R/C sections, under various tensile and compressive axial loads are shown. It is observed that the results from both studies are almost coincident, albeit it is pointed out that in Rosati et al. [5], the concrete area under steel fibers is not subtracted, resulting to minor differences (more pronounced for high compressive axial loads). Fig. 25 shows comparisons of moment interaction curves for a composite section reported in Chen et al. [16] (section D), between the present study and various researchers [8,15,16]
(N = 4120 kN). Similar results are also reported in [18], which was not included in the plots for simplicity. It is again shown that the differences between all studies are marginal, which verifies the accuracy of the present methodology. At this point it should be noted that during the comparison process some discrepancies were also detected. Fig. 26 shows moment (left) and axial-moment (right) interaction curves for a composite section reported in Chen et al. [16]. Initially, this section was analyzed using the present methodology and the one by Charalampakis and Koumousis [8] (using the pertinent software by the same Authors [34]), producing only marginal differences (for N = 3000 kN). These are attributed to the present fiber modeling of steel bars as opposed to the area representation assumed in [8]. However, it is noted that when the above curves are compared to the respective ones in Chen et al. [16] and Chiorean [18], significant differences are observed (Fig. 26, left). On the contrary, when the present method is directly compared to Chiorean [18] in terms of axial-moment capacity curve (Fig. 26, right), only marginal dif-
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
49
Fig. 23. Comparison of moment interaction curves for steel angle section, between present method and Charalampakis [33].
Fig. 24. Comparison of moment interaction curves for G-shaped and multicell section, between present method and Rosati et al. [5].
ferences are observed. This is a controversial issue, since the uniaxial moment capacity range corresponding to a moment interaction curve (here: M X ,max M X ,min = 153.2 kNm for M Y = 0 and N = 3000 kN) directly maps on a horizontal cut of the respective axial-moment interaction curve for the same axial load (N = 3000 kN, red arrows in Fig. 26, right), irrespective of the section origin type considered (plastic or geometric). It is therefore believed, that the above discrepancies are simply attributed to a modeling oversight. Actually, if the structural steel web is removed from the section, all the above differences are vanishing. At this point, the challenging computational task of handling concrete softening in the ultimate section response was investi-
gated. In Fig. 27, moment capacity curves are compared between the present study and Charalampakis and Koumousis [34] software for varying degrees of concrete softening (s), using the aforementioned composite section. It is observed that the analysis results from both methods are coincident, which demonstrates the robustness of the present methodology even for softening functions, where non-convex interaction curves may result under higher compressive axial loads (see also Fig. 19). In Fig. 28, the analysis is taken further producing full 3D failure surfaces with and without considering the softening effect. It is shown that the present method produces smooth surfaces without non-convergence issues, where increasing softening leads to smaller ultimate strength vol-
50
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
Fig. 25. Comparison of moment interaction curves for composite section reported in Chen et al. [16], between present method and various researchers [8,15,16].
Fig. 26. Comparison of moment and axial-moment interaction curves for composite section reported in Chen et al. [16], between present method and various researchers [8,16,18].
umes, with more pronounced pinching under higher compression levels. It should be noted here that the effect of concrete softening would be more pronounced (as in [18]) if the C-restriction criterion, which limits the quasi-concentric compressive strain, was (improperly) ignored. The final case is investigating whether the present method can successfully analyze large-scale geometries and simplistic material laws, which may possibly lead to numerical (e.g. floating point) issues, if not properly handled in computer code. In Fig. 29 (left), a comparison is presented between axial-moment interaction curves for a large reinforced spillway pier section, reported in Stefan and Léger [35], using linear elastic materials for concrete and steel. It is observed that the results between the two studies are coincident.
Moreover, the 3D interaction surface of a powerhouse section (consisting of two non-contiguous concrete regions) is depicted in Fig. 29 (right) with exactly the same shape as the one reported in [35]. From the above comparisons presented in this chapter, it is concluded that the present methodology not only exhibits a robust numerical behavior but also provides accurate results for a wide range of geometry and material characteristics. 5. Benchmarks
Apart from the indispensable numerical stability that should characterize any section analysis methodology, it is deemed that
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54
51
Fig. 27. Comparison of moment interaction curves for composite section reported in Chen et al. [16], between present method and Charalampakis and Koumousis [8], under different levels of concrete softening.
Fig. 28. Effect of concrete softening level on section interaction surfaces.
execution speed is also of great importance, since it is not only a repetitive task inside various structural analysis computational frameworks but also an indication of a successful computer code implementation. However, this issue is usually overlooked in similar studies, rendering any pertinent speed comparison difficult or impossible. In this study, a consistent benchmarking procedure is presented for the first time, including various benchmarking systems for executing various algorithmic tasks. Specifically, four representative personal computers of the last decade were employed, the detailed characteristics and performance indices of which are described in Table 3. On the above systems, three different numerical procedures were assigned, namely (a) one million section integrations I(d, h) with random values of d = [ 2h. . . 2h] and h = [0 . . .2p] ( h is the section height), (b) the construction of a 3D failure surface with 36,000 vertices using the N h method and (c) the same surface using the more demanding N a method. These procedures were executed for all five reference sections depicted in Fig. 17, and exe-
cution times are reported in Table 4. It is pointed out that for each numerical procedure, only the corresponding object-oriented programming (OOP) method was profiled, i.e. the reported times do not contain any processing overhead from graphics, logs etc. It is deemed that the efficiency of the present methodology in terms of speed is satisfactory, performing thousands of section integrations or 3D point evaluations per second. The main scope of this benchmarking report is to establish a solid reference basis for future studies, simply by applying (at least) linear extrapolation based on the above known system performance characteristics (Table 3). To the best of the Author’s knowledge, the only consistent benchmarking report in the literature is presented in Bonet et al. [1], where three bare concrete sections (solid rectangular, hollow rectangular and double-T) are integrated for different neutral axis depths and inclinations ( d, h). Integration is performed using various integration schemes, such as the Fafitis method [17], already described in paragraph 2.2 and the Modified Thick Layer Integra-
52
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
Fig. 29. Comparison of axial-moment interaction curves for reinforced spillway pier section (left) and 3D interaction surface for powerhouse section (right), between present method and Stefan and Léger [35].
Table 3
Characteristics of benchmarking systems.
a b c
Feature/index
Obsolete PC
Older PC
Regular PC
High-end PC
Processor CPU/FSB (MHz) SuperPIa LinXb (Linpack) 7-Zipc
Intel Pentium 4 2800/200 50.78 s 1.79 GFlops 1307 MIPS
AMD Athlon64 4400+ 2640/240 33.59 s 3.60 GFlops 2367 MIPS
Intel Core 2 Duo E6600 3000/333 20.55 s 9.23 GFlops 3011 MIPS
Intel Core i7 920 3800/200 11.43 s 12.48 GFlops 3823 MIPS
v1.5 mod XS, 1M digits (www.overclock.net/downloads/138140-super-pi.html). v0.64, 32 bit, 1 thread, size 5000 (www.xtremesystems.org/forums/showthread.php?201670-LinX-A-simple-Linpack-interface). v9.20, 32 bit, 32M dictionary, 1 thread, total rating (www.7-zip.org). All tools extracted on 7th September 2011.
tion (MTLI) method [1]. Since both the present and the MTLI methods are improved integration schemes originating from the Gauss/ Green concept introduced by Fafitis [17], it is deemed interesting to attempt a comparison between them in terms of execution speed. Of course, a direct comparison in terms of running time is not evenhanded, because system details in [1] are not available. For reference only, a direct comparison shows that the present method is at least 50 times faster than MTLI (using the ‘Older PC’ system of Table 3). However, a just alternative is to perform an indirect comparison procedure as follows: for both present and MTLI methods, a ‘speed index’ is defined as the time ratio between the method under consideration and the reference Fafits scheme (using three Gauss points), i.e. the percentage increase in computational time. This index will be evaluated against the reported accuracy for each method, which is expressed as the maximum normalized error of all resultant actions with respect to exact integration (using 48 or 100 Gauss points respectively). In Table 5, the results of this indirect comparison are presented. For the present method, a series of 21 section depth/orientation combinations (described in Table 3 in [1]) with 100,000 integrations per combination was performed and compared against the respective results for the MTLI method for 2, 3 and 4 Gauss points, reported in [1]. This procedure was repeated for both normal strength ( f c = 20 MPa) and high strength ( f c = 80 MPa) concrete, featuring a non-polynomial stress–strain law according to Model Code 90 [11], including concrete tensile response.
Table 4
Benchmarks of stress integration and interaction surface construction on various systems.
Section
Obsolete PC
Older PC
Regular PC
High-end PC
6.89 6.55 20.51 83.38 175.65
4.55 4.54 14.08 59.05 126.85
2.00 2.11 5.11 20.70 52.21
1.13 1.16 2.97 12.72 32.49
7.62 9.20 14.67 64.59 284.75
5.00 6.00 9.55 42.23 183.68
1 million integrations with random d, h (s)
A B C D E
16.80 16.71 49.93 204.77 435.97
8.72 8.71 27.97 113.02 247.93
N h interaction surface of 100 360 grid (s)
A B C D E
3.76 3.86 9.78 43.30 112.64
2.13 2.22 5.32 22.99 59.50
N a interaction surface of 100 360 grid (s)
A B C D E
17.22 20.82 32.98 146.76 670.68
9.04 10.97 17.69 75.99 338.76
The results of Table 5 are depicted as ‘speed index vs. maximum error’ plots in Fig. 30. It is observed that the present method is
53
V.K. Papanikolaou/ Computers and Structures 98-99 (2012) 33–54 Table 5
Performance comparison between present method and MTLI method by Bonet et al. [1].
Normal strength f c = 20 MPa
Rectangular
Hollow
Double-T
10.21 20.21 (+98%) 0.087%
19.77 34.20 (+73%) 0.086%
8.81 11.70 (+33%) 0.016%
Present method – total 21 100,000 = 2.1 million integrations
Simulated Fafitis method (3 Gauss points) (s) Present method (s) Error of present method (%) MTLI method – total 21 10 50 = 10,500 integrations
Simulated Fafitis method (3 Gauss points) (s) MTLI with 2 Gauss points (s) Error of MTLI with 2 Gauss points (%) MTLI with 3 Gauss points (s) Error of MTLI with 3 Gauss points (%) MTLI with 4 Gauss points (s) Error of MTLI with 4 Gauss points (%) High strength f c = 80 MPa
1.65 4.39 (+166%) 0.54% 5.94 (+260%) 0.0068% 6.97 (+322%) 0.000096%
3.13 5.66 (+81%) 0.54% 7.25 (+132%) 0.0068% 8.95 (+186%) 0.000096%
5.21 12.19 (+134%) 0.12% 15.32 (+194%) 0.0016% 18.62 (+257%) 0.000024%
Rectangular
Hollow
Double-T
10.80 37.12 (+244%) 0.71%
20.70 55.17 (+167%) 1.04%
8.74 18.44 (+111%) 0.70%
1.65 6.21 (+276%) 0.25% 8.02 (+386%) 0.1223% 10.00 (+506%) 0.0448%
3.07 9.12 (+197%) 0.34% 11.76 (+283%) 0.244% 14.28 (+365%) 0.10%
5.21 13.34 (+156%) 0.80% 16.75 (+221%) 0.39% 20.38 (+291%) 0.145%
Present method – total 21 100,000 = 2.1 million integrations
Simulated Fafitis method (3 Gauss points) (s) Present method (s) Error of present method (%) MTLI method – total 21 10 50 = 10,500 integrations
Simulated Fafitis method (3 Gauss points) (s) MTLI with 2 Gauss points (s) Error of MTLI with 2 Gauss points (%) MTLI with 3 Gauss points (s) Error of MTLI with 3 Gauss points (%) MTLI with 4 Gauss points (s) Error of MTLI with 4 Gauss points (%)
Fig. 30. Performance comparison between present and MTLI method by Bonet et al. [1]for normal and high strength concrete.
clearly faster than its MTLI counterpart for normal strength concrete, exhibiting smaller error levels with faster execution times for all three sections under consideration. As far as the high strength concrete is concerned, the difference between the two methods is less pronounced probably due to the brittle material behavior (very steep softening part in the concrete stress–strain law), albeit the present method generally performs faster under acceptable error levels, which can be further reduced if tighter tolerance criteria are used (Eq. (16)). It is therefore concluded that the suggested methodology exhibits a very satisfactory performance in both terms of accuracy and execution speed, fulfilling all prerequisite conditions imposed beforehand.
load, featuring various new algorithmic procedures regarding geometry manipulation, stress integration and solution strategies. With the aid of validation tests, comparisons with the existing literature and benchmarks, it was demonstrated that the suggested method shows very satisfactory performance in terms of stability and speed. It is therefore believed that it can provide a solid basis for further developments in the fields of structural design, assessment and optimization. Various extensions of the present method are currently investigated, specifically in the directions of section assessment against external loading, design for reinforcement, prestressed tendons, jacketed sections with preloading and design optimization.
6. Closure
References
In this study, a new methodology was presented for the analysis of fully arbitrary composite sections in biaxial bending and axial
[1] Bonet JL, Romero ML, Miguel PF, Fernandez MA. A fast stress integration algorithm for reinforced concrete sections with axial loads and biaxial bending. Comput Struct 2004;82:213–25.
54
V.K. Papanikolaou / Computers and Structures 98-99 (2012) 33–54
[2] Furlong RW. Ultimate strength of square columns under biaxially eccentric loads. ACI J. 1961;57:1129–40. [3] Penelis GG. Analytical investigation of the biaxial bending problem of R/C sections using plasticity theory (in Greek). Aristotle University of Thessaloniki; 1969. [4] Werner H. Schiefe Biegung polygonal umrandeter Stahlbeton-Querschnitte. Beton Und Stahlbetonbau. 1974;69:92–7. [5] Rosati L, Marmo F, Serpieri R. Enhanced solution strategies for the ultimate strength analysis of composite steel–concrete sections subject to axial force and biaxial bending. Comput Methods Appl Mech Eng 2008;197:1033–55. [6] Pallarés L, Miguel PF, Fernandez MA. A numerical method to design reinforced concrete sections subjected to axial forces and biaxial bending based on ultimate strain limits. Eng Struct 2009;31:3065–71. [7] Sousa JJBM, Muniz CFDG. Analytical integration of cross section properties for numerical analysis of reinforced concrete, steel and composite frames. Eng Struct 2007;29:618–25. [8] Charalampakis AE, Koumousis VK. Ultimate strength analysis of composite sections under biaxial bending and axial load. Adv Eng Softw 2008;39:923–36. [9] EN1992-1-1. Eurocode 2: design of concrete structures – Part 1-1: general rules and rules for buildings. European Committee for Standardization; 2004. [10] EN1998-2. Eurocode 8: design of structures for earthquake resistance – Part 2: Bridges. European Committee for Standardization; 2005. [11] CEB. CEB/FIP model code 1990. Lausanne, Bulletin d’Information CEB 213/214; 1993. [12] Popovics S. A numerical approach to the complete stress–strain curve of concrete. Cem Concr Res 1973;3:583–99. [13] Tsao WH, Hsu CTT. A nonlinear computer analysis of biaxially loaded L-shaped slender reinforced concrete columns. Comput Struct 1993;49:579–88. [14] Spacone E, Filippou FC, Taucer FF. Fibre beam-column model for non-linear analysis of R/C frames: part I. Formulation. Earthquake Eng Struct Dynam 1996;25:711–25. [15] Sfakianakis MG. Biaxial bending with axial force of reinforced, composite and repaired concrete sections of arbitrary shape by fiber model and computer graphics. Adv Eng Softw 2002;33:227–42. [16] Chen SF, Teng JG, Chan SL. Design of biaxially loaded short composite columns of arbitrary section. J Struct Eng, ASCE 2001;127:678–85. [17] Fafitis A. Interaction surfaces of reinforced-concrete sections in biaxial bending. J Struct Eng, ASCE 2001;127:840–6. [18] Chiorean CG. Computerised interaction diagrams and moment capacity contours for composite steel-concrete cross-sections. Eng Struct 2010;32: 3734–57.
[19] Zupan D, Saje M. Analytical integration of stress field and tangent material moduli over concrete cross-sections. Comput Struct 2005;83:2368–80. [20] Rodriguez-Gutierrez JA, Aristizabal-Ochoa JD. Biaxial interaction diagrams for short RC columns of any cross section. J Struct Eng, ASCE. 1999;125: 672–83. [21] Yen JYR. Quasi-Newton method for reinforced concrete column analysis and design. J Struct Eng, ASCE 1991;117:657–66. [22] Dundar C, Sahin B. Arbitrarily shaped reinforced concrete members subject to biaxial bending and axial load. Comput Struct 1993;49:643–62. [23] Yau CY, Chan SL, So AKW. Biaxial bending design of arbitrarily shaped reinforced concrete column. ACI Struct J 1993;90:269–78. [24] Martínez FJ, González-Vidosa F, Hospitaler A, Yepes V. Heuristic optimization of RC bridge piers with rectangular hollow sections. Comput Struct 2010;88:375–86. [25] Brøndum-Nielsen T. Ultimate flexural capacity of cracked polygonal concrete sections under biaxial bending. ACI J 1985;82:863–9. [26] De Vivo L, Rosati L. Ultimate strength analysis of reinforced concrete sections subject to axial force and biaxial bending. Comput Methods Appl Mech Eng 1998;166:261–87. [27] Alfano G, Marmo F, Rosati L. An unconditionally convergent algorithm for the evaluation of the ultimate limit state of RC sections subject to axial force and biaxial bending. Int J Numer Meth Eng 2007;72:924–63. [28] Kwan KH, Liauw TC. Computerized ultimate strength analysis of reinforced concrete sections subjected to axial compression and biaxial bending. Comput Struct 1985;21:1119–27. [29] Roik K, Bergmann R. Design method for composite columns with unsymmetrical cross-sections. J Constr Steel Res 1990;15:153–68. [30] ACI. Building code requirements for structural concrete ACI 318-95. Detroit: American Concrete Institute; 1995. [31] An W, Saadatmanesh H, Ehsani MR. RC beams strengthened with FRP plates. II: analysis and parametric study. J Struct Eng, ASCE 1991;117:3434. [32] Brent RP. Algorithms for minimization without derivatives. Englewood Cliffs, NJ: Prentice-Hall; 1973. [33] Charalampakis AE. Full plastic capacity of equal angle sections under biaxial bending and normal force. Eng Struct 2011;33:2085–90. [34] Charalampakis AE. myBiaxial 2.0: analysis of arbitrary composite sections in biaxial bending and axial load. National Technical University of Athens; 2005. [35] Stefan L, Léger P. Multicriteria capacity envelopes for biaxial bending of concrete hydraulic structures. J Struct Eng, ASCE 2010;136:1035–43.
