A computer method for nonlinear inelastic analysis of 3D composite steel–concrete frame structures

Engineering Structures 57 (2013) 125–152

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

A computer method for nonlinear inelastic analysis of 3D composite steel–concrete frame structures C.G. Chiorean ⇑ Faculty of Civil Engineering, Technical University of Cluj-Napoca, 15 C. Daicoviciu Str., RO-400020 Cluj-Napoca, Romania

a r t i c l e

i n f o

Article history: Received 21 March 2013 Revised 12 September 2013 Accepted 16 September 2013 Available online 18 October 2013 Keywords: Distributed plasticity Large deflection Tangent flexural rigidity Residual stresses 3D frameworks Advanced inelastic analysis

a b s t r a c t This paper presents an efficient computer method for nonlinear inelastic analysis of three-dimensional composite steel–concrete frameworks. The proposed formulation is intended to model the geometrically nonlinear inelastic behaviour of composite frame elements using only one element per physical member. The behaviour model accounts for material inelasticity due to combined bi-axial bending and axial force, gradual yielding is described through basic equilibrium, compatibility and material nonlinear constitutive equations. In this way, the states of strain, stress and yield stress are monitored explicitly during each step of the analysis, the arbitrary cross-sectional shape, various stress–strain relationships for concrete and steel and the effect of material imperfections such as residual stresses are accurately included in the analysis. Tangent flexural rigidity of cross-section is derived and then using the flexibility approach the elasto-plastic tangent stiffness matrix and equivalent nodal loads vector of 3-D beam-column element is developed. The method ensures also that the ultimate strength capacity of the cross-section is nowhere exceeded once a full plastified section develops. The proposed nonlinear analysis formulation has been implemented in a general nonlinear static purpose computer program. Several computational examples are given to validate the effectiveness of the proposed method and the reliability of the code to approach large-scale spatial frame structures. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, have witnessed significant advances in nonlinear inelastic analysis methods for steel and composite steel–concrete framed structures and integrate them into the new and more rational advanced analysis and design procedures [1,2]. Reliable nonlinear analysis tools are, for instance, essential in performance-based earthquake engineering, and advanced analysis methodologies, that involves accurate predictions of inelastic limit states up or beyond to structural collapse. A number of approaches have been proposed in the last years to model the nonlinear response of composite steel–concrete elements [3–24]. A detailed discussion about this issue can be found in [3,4]. There currently exist several methods and computer programs concerning the nonlinear inelastic analysis that calculate strength limit states of steel and composite steel–concrete frame structures. At one extreme, two-and three dimensional finite elements enhanced with advanced material constitutive laws [14,15,17] were used to investigate the nonlinear response of steel and composite steel–concrete frame members. Currently the available tools for such analysis are general purpose FE programs that require very

⇑ Tel./fax: +40 264 594967. E-mail address: [email protected] 0141-0296/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engstruct.2013.09.025

fine-grained modelling that is often impractical to the structural engineer. At the other extreme, the line elements approach in conjunction with either distributed or concentrated plasticity models, have been devoted to the development of nonlinear analysis tools for frames that provide a desirable balance between accuracy and computational efficiency [5–13,16–31]. In the concentrated plasticity approach [16,18–20] which is usually based on the plastic hinge concept, the effect of material yielding is ‘‘lumped’’ into a dimensionless plastic hinge. Regions in the beam-column elements other than at the plastic hinges are assumed to behave elastically. In the plastic hinge locations if the cross-section forces are less than cross-section plastic capacity, either elastic behaviour or gradual transition (refined plastic hinge) from elastic to plastic behaviour is assumed. The plastic hinge approach could eliminate the integration process on the cross section and permits the use of fewer elements for each member, and hence greatly reduces the computing effort. Unfortunately, as plastification in the member is assumed to be concentrated at the member ends, the plastic hinge model is usually less accurate in formulating the member stiffness, requires calibration procedures, but make possible to use only one element per physical member to simulate geometric and material nonlinearities in composite building frameworks [16,20]. In the distributed plasticity models gradual yielding and spread of plasticity is allowed throughout

126

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

cross-section and along the member length [5–14,17,22–32]. There are two main approaches that have been used to model the gradual plastification of members in a second-order inelastic analysis, one based on the displacement method or finite element approach [17,25] and the other based on the force or flexibility method [5– 14,22–24,26–32]. Because displacement based elements implicitly assumed linear curvatures along the element length, accuracy in this approach when material nonlinearity is taken into account can be obtained only using several elements in a single structural member, thus the computational effort is greatly enhanced and the method becomes prohibited computational in the case of large scale frame structures. On the other hand in the flexibility based approach only one element per physical member can be used to simulate the gradual spread of yielding throughout the volume of the members but the complexity of these methods derives from their implementation in a finite element analysis program and the inclusion of the element geometrical effects [26]. Mixed finite element approaches have been developed also to model composite beams with bond slip [3,5,6,8,22]. However in order to allow the concrete and steel to have independent displacements all these methods include additional degrees of freedom at the element ends. When modelling the semi-rigid composite frameworks some difficulties may arise enforcing the compatibility conditions at the semi-rigid composite connections [23]. In the efforts to develop an intermediate solution that has the computational efficiency of plastic hinge methods and the accuracy of distributed plasticity methods several researchers developed quasi-plastic hinge [23,30–32] or stress-resultant constitutive models [21]. Although subject to some limitations of required calibration these methods have been shown to make distributed plasticity analyses practical for large scale 3D steel [30,31] and composite steel–concrete frameworks [21,23], usually only one element per member is necessary to analyze. In spite of the availability of such nonlinear inelastic algorithms and powerful computer programs, the advanced nonlinear inelastic analysis of real large-scale composite steel–concrete frame structures still posses huge demands on the most powerful of available computers and still represents unpractical tasks to most designers. The present work attempts to develop accurate yet computational efficient tools for the nonlinear inelastic analysis of real large-scale 3D composite steel–concrete frameworks fulfilling the practical and advanced analysis requirements. Essentially, nonlinear inelastic analysis employed herein uses the accuracy of the fibre elements approach for inelastic frame analysis and address its efficiency and modelling shortcomings both to element level, through the use of only one element to model each physical member of the frame, and to cross-sectional level through the use of path integral approach to numerical integration of the cross-sectional nonlinear characteristics. This is an essential requirement to approach real large spatial frame structures, combining modelling benefits, computational efficiency and reasonable accuracy. Recent studies show that the slip effect between the steel and concrete interface has negligible influence on the global behaviour of multi-storey and high-rise fully connected steel–concrete composite frames [14]. In the proposed approach perfect bond between steel beam and concrete slab is assumed. Within the framework of flexibility based formulation a 3D frame element with 12 DOF able to take into account the distributed plasticity and element second order effects is developed. Comparing the proposed method with the related methods developed in [26–29] the present approach has several features that make the proposed element more practical in the context of implementation in finite element analysis program and posses accuracy comparable to that of fibre-flexibility or fibre-displacement finite elements.

As will be briefly described in the following sections, the element incremental stiffness matrix and the equivalent nodal loads are derived directly from energetic principles. In this way the elements of the stiffness matrix and equivalent nodal loads can be obtained analytically and readily evaluated by computing the correction coefficients that affect the elastic flexibility coefficients and equivalent loads. In this way numerical integrations are required only to evaluate these correction coefficients and not the entirely flexibility or stiffness matrix elements as in [26–29]. Besides, the effect of the transverse shear deformation can be readily included in the element formulation, both in stiffness matrix and equivalent nodal loads. The resulting flexibility matrix of the element may have both elastic and plastic contributions. During the loading process unsymmetrical distribution of plastic zone throughout the cross-section may occur and consequently there are coupling between axial force and bending moments in elastoplastic domain. The present formulation does not consider the plastic interaction terms relating the axial and bending terms in the flexibility/stiffness matrix of the element. However the neglected terms in flexibility matrix have only plastic contributions and may be ignored. This is obviously a simplification of the proposed approach whose acceptance must be justified by verification studies, but the resulting stiffness matrix does not incurring the expense of a detailed flexibility based methods [26–29]. In its final computerized implementation the proposed method is very similar to quasi-plastic hinge approaches [32]. Moreover, in the proposed approach, the effects of the discontinuity and/or discrete loading along element can be efficiently taken into account by writing a single moment equation in such a way that it becomes continuous for entire length of the element in spite of the discontinuity of loading. Thus the separate moment equation for each change of loading point is not required. The element stiffness matrix are evaluated in [26–29] by an iterative procedure carried out at the element level, nested in the iterative procedure adopted to solve the nonlinear global structural response [33]. Thus approximations in the strain distribution along the element length, in the control sections, are required in the force-based frame elements. This fact makes these methods to be more complicated in implementation in finite element analysis framework. On the other hand, in the proposed approach the element force fields are described by the second order transfer matrix as function of the nodal and applied element forces and the inelastic response of the cross-sections (control points) is rigorously evaluated by enforcing the equilibrium between external and internal forces for each cross-section by a global convergence iterative procedure. In this way gradual yielding throughout the crosssection subjected to combined action of axial force and bi-axial bending moments is described through basic equilibrium, compatibility and material nonlinear constitutive equations, the states of strain, stress and yield stress are monitored explicitly during each step of the analysis, the arbitrary cross-sectional shape and the effect of material imperfections such as residual stresses are accurately included in the analysis. Tangent flexural and axial rigidity of the cross-sections are explicitly derived and the inelastic response at the element level is determined by integrating the variable section flexural EIy and EIz and axial EA rigidity along the member length, depending on the bending moments and axial force level, cross-sectional shape and nonlinear constitutive relationships. Comparing the proposed formulation with those described in [23,25–29] another difference of the proposed approach refers to how the element geometrical effects are taken into account. In displacement-based formulations [23,25], the deformed shape of the element is obtained directly based on the nodal displacement values and the adopted shape functions. Thus the implementation of the element second-order effects are straightforward, but the accu-

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

racy is dependent by the number of the finite elements involved. On the other hand in flexibility based element there are no deformation shape functions to relate the deformation field inside the element to the nodal displacements hence more elaborate and more time consuming procedures are required. For instance, in order to capture the element geometrical effects in [27] is developed a curvature-based displacement interpolation (CBDI) function whereas in [29] Simpson integration scheme along with piecewise interpolation of the curvature is applied. In the proposed approach because the method requires tangent flexural rigidities to be computed along the element length, the stability stiffness functions can be easily applied as corrections of the element elasto-plastic tangent stiffness matrix, and updating at each load increment the length, axial force and the flexural rigidity about of each principal axes of the element. This method minimizes the modelling and solution time and generally only one element are needed per member in order to simulate the second order effects. Although subjected to some limitations of required calibration this approach has been shown to make inelastic second-order analyses practical for large building frames [30,31]. Recent studies [34–37] show that one of the most critical element that affects both the accuracy and speed of the inelastic section response and consequently the inelastic element analysis procedure is the stress integration scheme. Most of the existing methods implemented in flexibility based or displacement based fibre elements for the inelastic analysis of cross-sections rely on the numerical integration of stress resultants using the well known ‘‘fibre decomposition method’’ where the cross section is decomposed in filaments and the section response is computed by composing the uniaxial behaviour of each filament. These techniques are not numerically efficient due to large amount of information needed to characterize the section and the high number of operations required by stress integration with an allowable error level. A further complication is represented by the nonlinear constitutive laws that are usually assumed for concrete in compression and tension. These laws are defined by piecewise functions and there is no continuity in the derivative and the classical integration methods can produce important integration errors [35,36]. As will be briefly described in this paper in the proposed approach an improved adaptive Gauss–Lobatto numerical integration scheme on a Green path integral is applied demonstrating fast execution and accuracy for stress integration throughout the cross-sections. This approach is extremely rapid because stress integrals need only be evaluated at a small number of points on the section boundary and rapid convergence is assured by the inclusion of exactly determined tangent stiffnesses. A closely related issue to the aforementioned stress integration represents the incorporation and evaluation of the residual stress effects for encased steel section on the carrying capacity and inelastic behaviour during the loading process of the composite steel–concrete sections. Only a few studies available in the literature [38,39] have addressed this effect and detailed finite element studies are rather scarce in the open literature. During the loading process the equilibrium condition may be violated as the applied bending moments at the section are greater than the moment capacities, which is not acceptable for maximum strength analysis of members. An important feature of the present approach is represented by the capacity to determine directly the ultimate bending moments in order to check that they fulfil the ultimate limit state condition. Such a procedure is absolutely necessary in order to ensures that the ultimate strength capacity of the cross-section is nowhere exceeded once a full plastified section develops and to enforce the member equilibrium forces. Using an updated Lagrangian formulation, the global geometrical effects are considered updating the geometry of the structure at each load increment. The combined effects of material and

127

geometrical nonlinearity sources are simulated into a computer program automatically. In order to verify the efficiency and the accuracy of the proposed procedure and the developed code, this computer program was used to study the nonlinear behaviour of several composite steel–concrete frame structures that has been studied previously by other researchers. Several examples including simply column, simply-supported beams, continuous beams, plane frames and large scale 3D frames were analysed and compared with available experimental and numerical results. The examples run and the comparisons made prove the effectiveness of the proposed numerical procedure. The proposed software is presented as an efficient, reliable tool ready to be implemented into design practice for advanced analysis and seismic performance evaluation of spatial composite frame structures. 2. Formulation of the proposed analysis method The following assumptions are adopted in the formulation of analytical model: (1) plane section remain plane after flexural deformation; no slip occurs at the steel–concrete interface, full composite action is considered (2) warping and cross-section distortion are not considered; (3) torsional buckling do not occur; and (4) small strain but arbitrarily large displacements and rotations are considered. The proposed approach is based on the most refined type of second order inelastic analysis, distributed plasticity model, where elasto-plastic behaviour is modelled accounting for spread-of plasticity effect in sections and along the element and employs modelling of structures with only one line element per member, which reduces the number of degree of freedom involved and the computational time. The above assumptions allow the formulation details to be considered on two distinct levels, namely, the cross-sectional level and the member longitudinal axis level. Thus the nonlinear inelastic response of a beam-column element can be computed as a weighted sum of the inelastic response of a discrete number of cross-sections (i.e. stations) that are located at the numerical integration scheme points. Although the present study concerns mainly on frames with rigid joints and the effect of the floor slab action is ignored, in the proposed approach can be easily implemented the effects of the nonlinear behaviour of semi-rigid connections, with proper nonlinear moment-relative rotation models for composite connections, and the floor slab action. For instance the mathematical model described in [31] can be useful to include the effects of the semirigidity and the penalty element method can be used, as described in [31], to include the effect of the rigid floor diaphragm action. Alternatively, the floor slab may be modelled as flat shell elements (with 6 degrees of freedom per node) and coupled in this way with the beam-column element developed in this paper. Usually it is assumed that the floor slab remain elastic until the collapse of the frame, consequently only one flat shell element is required to model the stiffness of the slab. However, under extreme loadings which may cause severe damage in buildings, tensile membrane action of composite slab, which may significantly affect the capacity and response of the frame, is not taken into account in the proposed approach. 2.1. Elasto-plastic cross section analysis The elasto-plastic cross-section stiffness may be modelled by explicit integration of stresses and strains over the cross-section area (e.g., as micro model formulation) or through calibrated parametric equations that represent force-generalized strain curvature response (e.g. macro model formulation). In the macro model approach the gradual plastification through the cross-section sub-

128

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

jected to combined action of axial force and biaxial bending moments may be described by moment–curvature-thrust (M–U–N), and moment-axial deformation-thrust (M–e–N) analytical type curves that are calibrated either by numerical or experimental tests. Due to the fact that the inelastic behaviour of composite cross-sections is influenced by a various effects such as the shape of stress–strain relationship for concrete and steel, the geometrical shape of the cross-section, and material imperfections, rigorous analytical relationships are difficult to develop. In the micro model formulation, as proposed in the present paper, gradual plastification through the cross-section subjected to combined action of axial force and bi-axial bending moments is described through basic equilibrium, compatibility and material nonlinear constitutive equations. In this way, the states of strain, stress and yield stress are monitored explicitly during each step of the analysis, the arbitrary cross-sectional shape and the effect of material imperfections such as residual stresses are accurately included in the analysis. 2.1.1. Basic assumptions and constitutive material models Consider the cross-section subjected to the action of the external bending moments about both global axes and axial force as shown in Fig. 1. The cross-section may assume any shape with multiple polygonal or circular openings. It is assumed that plane section remains plane after deformation. This implies a perfect bond between the steel and concrete components of a composite concrete-steel cross section. Shear and torsional interaction effects are not accounted for in the steel and concrete constitutive models. 2.1.1.1. Behaviour of concrete in compression. The constitutive relation for concrete under compression is represented by a combination of a second-degree parabola (for ascending part) and a straight line (for descending part), Eq. (1), as depicted in Fig. 2:

8   2 > < fc 2 ee  ee2 ; e 6 e0 c0 fco ¼   c0  > : fc 1  c eec0 ; e0 < e ecu ec0

ð1Þ

Fig. 2. Stress–strain relationship for concrete in compression.

where fc represents the prism compressive strength and c represents the degree of strain-softening in the concrete and allows for the modelling of strain-softening effect and creep in the concrete by simply varying the crushing strain ec0, ultimate compressive strain ecu and c respectively. The prism compressive strength fc is taken as 0.76fcu, where fcu represents the cubic compressive strength and can be approximately evaluated as 1:25fc0 where fc0 represents the cylinder compressive strength. 2.1.1.2. Behaviour of concrete in tension. Neglecting tension strength of concrete could lead to a loss in the smoothness of moment–curvature curves due to the sudden drop in stress from the cracking strength to zero at the onset cracking. In addition, tension strength of concrete results in a small change in peak strength, but this is usually negligible. The model to account for tension strength, developed by Vecchio and Collins [40] is taken into account in the present investigation. The model of concrete in tension can be given in the following analytical form (Fig. 3):

Fig. 1. Arbitrary composite steel–concrete cross-section.

129

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

Residual stresses may be incorporated in the analysis. The magnitude and distribution of residual stresses in hot-rolled members depend on the type of cross-section and manufacturing processes and different patterns are proposed. In the US, the residual stress is considered constant in the web although when the depth of a wide flange section is large, it varies more or less parabolically (Fig. 5b). Another possible residual stress pattern in the web is the one simplified by a linear variation as used in European calibration frames (Fig. 5a).

Fig. 3. Stress–strain relationships for concrete in tension.

( ft ¼

Et e;

e 6 ecr pffiffiffiffiffiffiffi f ; e > ecr 1þ 500e cr

ð2Þ

a1 a22

where Et denotes the modulus of elasticity of concrete in tension; fcr represents the tensile strength of concrete; ecr is the concrete cracking strain; a1 is a factor that takes into account the bounding characteristics of the reinforcement and a2 represents a factor that takes into account the effects of load duration and cyclic loads. The tensile strength of concrete fcr is obtained as 1.4 (fc/10)2/3 and the tensile elastic modulus before cracking Et is assumed the same as compressive elastic modulus. As illustrated in Fig. 3 a slow rate of tension softening is assumed for the concrete in tension. 2.1.1.3. Behaviour of structural steel. A multi-linear elastic–plastic stress–strain relationship, both in tension and in compression, is assumed for the structural steel. In this way the strain-hardening effect may be included in analysis. The analytical model can be given in the following form (Fig. 4):

8 > < Es e; jej 6 jesy1 j   fs ¼ sgnðeÞfsy1 þ Esh1 e  sgnðeÞ  esy1 ; esy1 < jej 6 esy2   > : sgnðeÞfsy2 þ Esh2 e  sgnðeÞ  esy2 ; esy2 < jej 6 esu

ð3Þ

where Es is the Young modulus, fsyi denotes the yield stresses, esyi represents the yield strains, esu the ultimate steel strain and Eshi represents the slopes of the yielding branch, and sgn() represents the signum function, returns 1 for negative values and +1 for positive values.

2.1.1.4. Behaviour of reinforcing bars. When the steel bar is subjected to tension, the crack in concrete will lead to the inhomogeneous distribution of stress of the steel bar along the longitudinal direction. Based on experimental results and theoretical analysis in [41] has been proposed a method for considering the inhomogeneous distribution of stress and smeared crack model. The stress– strain relationship of the steel bar in tension concrete region could be calculated as (Fig. 6):

( fr ¼

Es e;   fyr ð0:91  2BÞ þ ð0:02 þ 0:25Be=eyr Þ ;

e 6 enr e > enr

ð4Þ

1:5

ðfcr =fyr Þ where the parameter B ¼ , q = longitudinal reinforcement q steel ratio (limited to a minimum of 0.25%), the modified yield strain of the steel bar is enr = eyr(0.93 - 2B)/(1 - 0.25B), the modified yield stress is fnr = Esenr, and the ultimate average strain is eur = eyr(0.07 + 2B)/(0.25B). The hardening effect is not considered when the steel bar bearings compression load, therefore, a perfect elasto-plastic model for compression is assumed. 2.1.2. Elasto-plastic flexural rigidity of cross-section Considering the cross-section subjected to the action of the external bending moments (My, Mz) about each global axes and axial force (N) as shown in Fig. 1. The origin of the reference axis is usually considered in the geometric centroid of the cross-section. Under the above assumptions the resultant strain distribution corresponding to the curvatures about global axes U ¼ bUy Uz c and the axial strain e0 can be expressed in point r ¼ ½y z in a linear form as:

e ¼ e0 þ Uy z þ Uz y þ er ¼ e0 þ UrT þ er

ð5Þ

where er represents the initial deformation produced by residual stresses, and is taken into account in above equation only for the structural steel. The equilibrium is satisfied when the external

Fig. 4. Stress–strain relationship for steel.

130

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

computed. According to the Newton iterative method, the iterative changes of unknowns vector X can be written as: 1

Xkþ1 ¼ Xk  F0 ðXk Þ FðXk Þ;

k0

ð8Þ

where F’ represents the Jacobian (or tangent cross-section stiffness matrix) of the nonlinear system (6) and can be expressed as:

2 F0 ¼



@N int

6 @ e0 6 @Mint @F y ¼6 6 @ e0 @X 4 int

@M z @ e0

@Nint @/y

@Nint @/z

@Mint y @/y

@Mint y @/z

@Mint z

@Mz @/z

@/y

3

7 7 7 7 5 int

ð9Þ

Explicitly the expressions of the Jacobian’s coefficients are given in Eq. (10). k11 ¼ Fig. 5. Residual stress patterns.

FðXÞ ¼ f

int

f

ext

¼0

ð6Þ T

where X ¼ ½ e0 /y /z  , the external and internal loading vectors can then be represented by:

2 f

ext

N

6 7 ¼ 4 M y 5; Mz

2

3 f

int

R

3 rðeðe0 ; /y ; /z ÞÞdAcs Acs 6 int R 7 7 ¼6 4 My ¼ Acs rðeðe0 ; /y ; /z ÞÞzdAcs 5 R M int z ¼ Acs rðeðe0 ; /y ; /z ÞÞydAcs Nint ¼

ð7Þ

in which e0, /y, /z represent the unknowns and the surface integral is extended over concrete and structural steel areas (Acs). The above system can be solved numerically using, for instance, the load-controlled Newton method and taking into account the fact that the stresses are implicit functions of the axial strain and curvatures through the resultant strain distribution given by Eq. (5). Eq. (6) are solved numerically using the Newton–Raphson method, and results in three recurrence relationships to obtain the unknowns e0, /y, /z and then flexural EI and axial EA rigidity modulus can be

Z



rðeðe0 ; /y ; /z ÞÞdAcs ¼

Acs

Z



Z Acs

@N int @ ¼ @/y @/y

rðeðe0 ; /y ; /z ÞÞdAcs ¼

Z

@r @e dAs ¼ @ e @ e0

Z

ET dAcs

Acs

Z @r @e dAcs ¼ ET zdAcs @ e @/ Acs Acs Acs y  Z Z Z @N int @ @r @e k13 ¼ ¼ rðeðe0 ; /y ;/z ÞÞdAcs ¼ dAs ¼ ET ydAcs @/z @/z Acs Acs @ e @/z Acs Z Z Z int @M y @ @r @e k21 ¼ ¼ rðeðe0 ;/y ; /z ÞÞzdAcs ¼ zdAs ¼ ET zdAcs @ e0 Acs @ e0 Acs @ e @ e0 Acs Z Z Z int @M y @ @r @e k22 ¼ ¼ rðeðe0 ;/y ; /z ÞÞzdAcs ¼ zdAs ¼ ET z2 dAcs @/y Acs @/y Acs @ e @/y Acs Z Z Z @M int @ @r @e y k23 ¼ ¼ rðeðe0 ; /y ; /z ÞÞzdAcs ¼ zdAcs ¼ ET yzdAcs @/z Acs @/z Acs @ e @/z Acs  Z Z Z @M int @ @r @e z k31 ¼ ¼ rðeðe0 ;/y ; /z ÞÞydAcs ¼ ydAcs ¼ ET ydAcs @ e0 Acs @ e0 @ e @ e 0 Acs Acs Z Z Z int @M z @ @r @e k32 ¼ ¼ rðeðe0 ;/y ; /z ÞÞydAcs ¼ ydAcs ¼ ET yzdAcs @/y Acs @/y Acs @ e @/y Acs Z Z Z int @M z @ @r @e k33 ¼ ¼ rðeðe0 ; /y ; /z ÞÞydAcs ¼ ydAcs ¼ ET y2 dAcs @/z Acs @/z Acs @ e @/z Acs k12 ¼

forces (N, My, Mz) are equal to the internal ones. These conditions can be represented mathematically in terms of the following nonlinear system of equations as:

@N int @ ¼ @ e0 @ e0

ð10Þ

These coefficients are expressed in terms of the tangent modulus of elasticity Et = dr/de. The iterative procedures starts with the all unknowns e0, /y, /z set to zero and the solutions will be computed in just a few iterations by applying the rapid locally convergent Newton iterative procedure given by Eq. (8) enhanced with a line search algorithm such that global convergence can be achieved. By a globally convergence algorithm we mean that for

Fig. 6. Stress–strain relationship for reinforcing bars.

131

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

any initial iterate the iteration either converges to a solution or fails to do so in one way. When iteration fails to converge, within a prescribed number of iterations, this means that the external loads exceeded ultimate strength capacity of the cross-section. In this case ultimate strength capacity procedure is applied as will be described in the next sections. The convergence criterion is expressed as a ratio of the norm of the out-of-balance force vector to the norm of the total applied load. So the solution is assumed to have converged if:

pffiffiffiffiffiffiffiffi FT F pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 TOL ext T ext f f

ð11Þ

where TOL is the specified computational tolerance, usually taken as 1E  5. In this way, for given bending moments and axial force we can obtain the strain and stress distribution throughout the cross-section and then the axial and flexural rigidity of the cross-section can be computed, as will be briefly described below. The tangent stiffness matrix for the cross section relates small changes in the deformations to small changes in the actions sustained. The incremental relationships between incremental efforts and incremental deformations can be expressed as:

2

k11 6 4 k21 k31

k12 k22 k32

3 2 3 De0 DN 7 6 7 6 7 k23 5  4 D/y 5 ¼ 4 DM y 5 D/z DM z k33 k13

3 2

ð12Þ

ð13Þ

where kt represents the tangent stiffness matrix of cross-section, DF represents the incremental forces (axial force and bending moments about each principal axes of cross-section) and Du represents the incremental deformations (axial deformation and curvatures). The coefficients of the tangent stiffness matrix kij are given by Eq. (10) and are evaluated considering the strains and stresses at equilibrium for given external actions. Inverting the relationship (13) gives:

f t DF ¼ Du

DN ¼ 0;

DM z ¼ 0

DN ¼ 0;

DM y ¼ 0

ð15Þ

Considering the system (14) with the unknown DF the flexural rigidity about y axis is determined solving the system with the following imposed restrictions:

DN ¼ 0 DM z ¼ 0

ð16Þ

The incremental forces DF can be expressed, in this case, as:

2 3 0   ^ ¼ 4 1 5 DM y DF ¼ TDF 0

ð17Þ

Substituting Eq. (17) for DF into Eq. (14) gives:

^ ¼ Du f t TDF

ð18Þ

^ we Multiply both members of Eq. (18) with T and solving for DF obtain: T

1

^ ¼ ðTT f t TÞ TT Du DF

ð19Þ

Eqs. (17) and (19) can now be combined by realizing that the basic force–deformation relationship is given by: 1

DF ¼ TðTT f t TÞ TT Du

ð20Þ

or in explicitly matrix form:

or in a condensed matrix form:

kt  Du ¼ DF

DM y ; D/y DM z EItz ¼ ; D/z

EIty ¼

ð14Þ

2

3 2 DN 0 0 6 7 6 4 DMy 5 ¼ 4 0 EIty 0 0 DM z

3 De0 7 6 7 0 5  4 D/y 5 0 D/z 0

3 2

where EIty represents the tangent flexural rigidity about x axis of cross-section and is given by the following relation: 2

EIty ¼

2

2

k11 k22 k33  k11 k23  k33 k12 þ 2k12 k13 k23  k22 k13 k11 k33 

2 k13

ð22Þ

Applying a similar approach the tangent flexural rigidity about z axis and tangent axial rigidity can be obtained as: 2

1

where f t ¼ kt represents the tangent flexibility matrix of crosssection. We define the tangent flexural rigidity about one axis of crosssection as a ratio between incremental bending moment and incremental curvature about that axis while keeping constant the axial force and bending moment about the other axis (Fig. 7):

ð21Þ

EItz ¼

2

2

k11 k22 k33  k11 k23  k33 k12 þ 2k12 k13 k23  k22 k13 k11 k22 

2 k12

k12 k21 k12 k23 k31 k22  k32 k21  k22 k22 k33 k22  k32 k23 k31 k22  k32 k21  k13 k33 k22  k32 k23

ð23Þ

EAt ¼ k11 

ð24Þ

It is important to note that, despite the fact that the position of the reference axis is kept fixed in evaluation of the elasto-plastic cross-section rigidities, the effect of unsymmetrical development of the plastic zones throughout the cross-section is efficiently taken into account by Eqs. (22)–(24). For instance, in the case of uniaxial bending for symmetric cross-sections (considering the bending about z axis) it can be readily derived from Eq. (23) that the tangent flexural rigidity of the cross section can be computed using the following relationship:

EItz ¼

Z Acs

R

Et y2 dAcs 

Acs

R

Et ydAcs

Acs

Et dAcs

2

ð25Þ

Now it can be easily observed that the tangent flexural rigidity can be further expressed as:

EItz ¼

Z

Acs

Fig. 7. Tangent flexural rigidity definition.

where

Et ðy  yp Þ2 dAcs

ð26Þ

132

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

Fig. 8. Composite beam with full composite action.

R

yp ¼ RAcs

Et ydAcs

ð27Þ

E dAcs Acs t

and represents the coordinate of the plastic centroid of the crosssection computed about the fixed reference axis. The above relationships can be applied even for evaluation of the elastic bending stiffness (about z axis) of the composite cross-section, as shown in Fig. 8. Indeed, considering the origin of the reference axis in the geometric centre of the steel beam, when apply Eq. (25) the expression of the bending stiffness can be reduced at the following relationship:

ðEIÞ0 ¼ Ec Ic þ Es Is þ

E c Ac  E s As 2 r Ec Ac þ Es As

ð28Þ

which represents the elastic flexural rigidity of the composite beam with full composite action and where EsAs, EcAc represents the axial stiffness of the steel and concrete components respectively and EsIs, EcIc are flexural rigidity of steel and concrete components respectively computed about their own reference axis. 2.1.3. Ultimate strength capacity evaluation A particularly important feature of the present approach is represented by the capacity to determine directly the ultimate bending moments for a given value of axial force and bending moment ratio in order to check that they fulfil the ultimate limit state condition (Fig. 9). Such a procedure is absolutely necessary in order to ensures that the ultimate strength capacity of the cross-section is nowhere exceeded once a full plastified section develops and to

enforce the member forces to move on the plastic surface during the loading steps. The cross-section subjected to bi-axial bending moments and axial force, reaches its failure limit state when the strain in the extreme fibre, attains the ultimate value. Consequently, at ultimate strength capacity the equilibrium is satisfied when the external forces are equal to the internal ones and in the most compressed concrete or tensioned steel fibre the ultimate strain is attained. These conditions can be represented mathematically in terms of the following nonlinear system of equations as:

  8R r eðe0 ; /y ; /z Þ dAcs  N ¼ 0 > Acs > >   >R < r eðe0 ; /y ; /z Þ ydAcs  Mz ¼ 0 Acs   R > r eðe0 ; /y ; /z Þ zdAcs u  My ¼ 0 > > Acs > : e0 þ /y zc ð/y ; /z Þ þ /z yc ð/y ; /z Þ  ecu ¼ 0

ð29Þ

in which N, My, Mz, e0, /y, /z represent the unknowns, the surface integral is extended over concrete and steel area (Acs). In Eqs. (29) the first three relations represent the basic equations of equilibrium for the axial load N and the biaxial bending moments My, Mz respectively, given in terms of the stress resultants. The last equation represents the ultimate strength capacity condition in which yc(/y, /z) and zc(/y, /z) represent the coordinates of the point in which this condition is imposed (Fig. 1). The coordinates of the ‘‘constrained’’ point can be always determined for each inclination of the neutral axis defined by the parameters /y and /z, and ecu represents the ultimate strain. The numerical procedure consists on the solution of the nonlinear system (29) for the following linear constraint:

Fig. 9. Interaction diagram for given bending moment’s ratio. The plastic surface requirements.

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

L1 ðN; My ; M z Þ  N  N0 ¼ 0 L2 ðN; My ; M z Þ  M z  tanðaÞM y ¼ 0

ð30Þ

where N0, represents the given axial force and bending moment ratio is given by tan(a). Corresponding to the linear constraint (30) we can define a point on the failure surface (Fig. 9) for a (given) fixed axial load (N0) and a given bending moment’s ratio. For each inclination of the neutral axis defined by the parameters /y and /z the farthest point on the compressed or tensioned side is determined (i.e. the point with co-ordinates yc, zc). Assuming that the failure condition is achieved in this point, the resulting strain distribution corresponding to the curvatures /y and /z can be expressed in linear form as:

eð/y ; /z Þ ¼ ecu þ /z ðy  yc Þ þ /y ðz  zc Þ

ð31Þ

Then, substituting the strain distribution given by Eq. (31) in the basic equations of equilibrium, the unknown e0 together with the failure constraint equation can be eliminated from the nonlinear system (29). Thus, the nonlinear system of Eq. (29) is reduced to an only three basic equations of equilibrium and together with the linear constraints (Eq. (30)), forms a determined nonlinear system of equations (with only two nonlinear equations), and the solutions can be obtained iteratively following the procedure described in [34]. This procedure will be used further in the present formulation (see Section 2.2.3) to determine the plastic bending moment capacities associated to a given value of axial force. 2.1.4. Evaluation of tangent stiffness and stress resultant Based on Green’s integration formula according to which the domain integrals appearing in the evaluation of internal resultant efforts and tangent stiffness matrix coefficients of the section can be evaluated in terms of boundary integral. This approach is extremely rapid because stress integrals need only be evaluated at a small number of points on the section boundary and rapid convergence is assured by the inclusion of exactly determined tangent stiffnesses and, of great importance, it is assure convergence for any load case. For this purpose, it is necessary to transform the variables first, so that the stress field is uniform in a particular direction, given by the current position of the neutral axis [34]. This is achieved by rotating the reference axes x, y to n, g oriented parallel to and perpendicular to the neutral axis, respectively. As the integration area contour is approximated by a polygon, the integral over the perimeter L, can be obtained by decomposing this integral side by side along the perimeter:

I

hðgÞnp dg ¼

L

nL Z X

giþ1

hðgÞnp dg

ð32Þ

gi

i¼1

where nL is the number of sides that forms the integration area. The sides are defined by the ng co-ordinates of the end-points as shown in Fig. 1. When the integration area is a circle with radius R, the integral over the perimeter L can be obtained by decomposing this integral as:

I L

hðgÞnp dg ¼

Z

R

hðgÞðR2  g2 Þ

p=2

dg þ ð1Þp

R



Z

R

p=2

hðgÞðR2  g2 Þ

dg

ð33Þ

R

This leads to a significant saving in imputing the data to describe the circular shapes, without the need to decompose the circular shapes as a series of straight lines and approximate the correct solution when circular boundaries are involved and allows efficiently to handle various circular shapes such as fillet regions which define the exact geometry of the structural steel profiles

133

(Fig. 8). In order to perform the integral on a determined side of the contour (Li), polygonal or circular, of the integration area, an adaptive interpolatory Gauss–Lobatto method is used [34]. The effect of residual stresses may be included in the analysis providing that the residual stress can be linearized for individual zones in the steel section associated to variations of residual stresses throughout the height of cross-section. For instance assuming the EC3 distribution of residual stresses the cross-section has to be divided in six regions as depicted in Fig. 10. In this way for each region the total strain in can be expressed as:

e ¼ e0 þ /y  z þ /z  y þ er

ð34Þ

where er represents a linear residual strain field which can be expressed for each particularly region as:

er ¼ a1 þ a2 z þ a3 y

ð35Þ

Next, the integration of the stress resultant and stiffness coefficients over the steel cross-section will be transformed into line integrals along the perimeter of the cross-section as already described, but in this case the reference axes are rotated for each region using the following value for angle h:

tan h ¼

/y þ a2 /z þ a3

ð36Þ

2.2. Second-order inelastic member analysis 2.2.1. Elasto-plastic tangent stiffness matrix and equivalent nodal loads Flexibility-based method is used to formulate the distributed plasticity model of a 3D frame element (12 DOF) under the above assumptions. An element is represented by several cross sections (i.e. stations) that are located at the numerical integration scheme points (Fig. 11). The spread of inelastic zones within an element is captured considering the variable section flexural EIy and EIz and axial EA rigidity along the member length, depending on the bending moments and axial force level, cross-sectional shape and nonlinear constitutive relationships. The elasto-plastic sectional rigidities are evaluated based on the iterative procedure already described at Section 2.1.2. Fig. 12 shows the 3D beam-column element in local system attached to the initially straight centre line, with the rigid body modes removed. Nonlinear analysis by the stiffness method requires incremental loading, i.e. the inelastic behaviour is approximated by a series of elastic analysis. The element incremental flexibility matrix fr which relates the end displacements to the actions Dsr and the elastoplastic equivalent nodal forces transferred to the nodes, can be derived directly from energetic principles. The elasto-plastic equivalent nodal forces transferred to the nodes, from the member loads, will not be constant during the analysis, and will be dependent on the variable flexural rigidity along the member according with the process of gradual formation of plastic zones. The equivalent nodal forces will be computed in order to accommodate member lateral loads and the plastic strength surface requirements. Let us consider the element in Fig. 12 subjected to nodal bending moments (DMi, DMj), uniform distributed loads (Dq), the concentrated forces at the ‘‘a’’ location (DP) and the concentrated bending moments at the ‘‘b’’ location and at the ‘‘j’’ end of the member (DMpb, DMpj). The concentrated bending moments are considered here as plastic correction moments, taken into account only when in the respective cross-section forces get to the full plastic capacity and constraining in this way that the section forces to move on the plastic surface during the loading process. It is obvious that in this case the loading conditions change along the span

134

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

Fig. 10. Discretization of the structural steel. Incorporation of residual stresses.

Fig. 11. Beam column element with 12 DOF.

of beam, and consequently there is corresponding change in moment and shear forces equations. This requires that a separate mo-

ment equation be written between each change of load point. These complications can be avoided by writing single moment

135

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

DW ¼

Z

L 2

0

L þ 2

Z 2 DM 2z ðnÞ L 1 DM y ðnÞ dn þ dn EIz ðnÞ 2 0 EIy ðnÞ 0 Z Z 2 DM2x L 1 DT y ðnÞ L 1 DT 2z ðnÞ dn þ þ 2 0 GAy ðnÞ 2 0 GAz ðnÞ GIt ðnÞ

DN 2 L dn þ 2 EAðnÞ

1

Z

1

0

Z

1

ð42Þ

Using the second theorem of Castigliano the relationship between incremental deformations and efforts can be readily calculated and partitioned as follows:

2

u

2

3

@ DW @ DN

3

3 2 3 DN 0 6 DM 7 6 d 7 iy 7 iy 7 6 6 # 6 7 6 7 6 DMjy 7 6 djy 7 0ð33Þ 7 6 7 6 6 DM 7 þ 6 d 7 f 2ð33Þ iz 7 6 6 iz 7 7 6 7 6 4 DM jz 5 4 djz 5 2

6 @ DW 7 6h 7 6 @ DM iy 7 7 6 iy 7 6 6 " 6 7 6 @ DW 7 6 hjy 7 6 @DMjy 7 f 1ð33Þ 7 6 7¼6 ¼ 6 h 7 6 @ DW 7 7 0ð33Þ 6 iz 7 6 @DMiz 7 6 7 6 4 hjz 5 6 @ DW 7 7 4 @ DMjz 5 hx @ DW

Fig. 12. Beam column element with rigid body modes removed.

DM x

ð43Þ

0

@ DM x

or in a condensed form: equation in such a way that it becomes continuous for entire length of the beam in spite of the discontinuity of loading. In this paper, Macaulay’s method [42] is applied for cases of discontinuous and/or discrete loading. Macaulay functions represent quantities that begin at a point a. Before point a the function has zero value, after point a the function has a defined value:

F n ðxÞ ¼ hx  ain ¼

0;

when x 6 a

ðx  aÞn

when x > a

ð37Þ

where n = 0, 1, 2, . . . When the exponent n = 0, we have:

F 0 ðxÞ ¼ hx  ai0 ¼

0; when x 6 a 1

when x > a

ð38Þ

Dur ¼ f r  Dsr þ dr

ð44Þ

where fr represents the incremental flexibility matrix of the beamcolumn element without rigid body modes, and in which the matrices fi (i = 1, 2) have the following expressions: 2 R1 1 3 L 0 EAðnÞ dn 0 0 6 7 R R R R 2 1 1 1 1 6 7 f1 ¼ 6 0 L 0 ðn1Þ dn þ 1L 0 GAdn L 0 nðn1Þ dn þ 1L 0 GAdn 7 EIy ðnÞ EIy ðnÞ y ðnÞ y ðnÞ 5 4 R 1 nðn1Þ R R R 2 1 n 1 1 dn 1 1 dn 0 L 0 EIy ðnÞ dn þ L 0 GAy ðnÞ L 0 EIy ðnÞ dn þ L 0 GAy ðnÞ 2 R1 L 6 R0 6 f 2 ¼ 6 L 01 4

ðn1Þ2 EIz ðnÞ

dn þ 1L

ðn1Þ2 EIz ðnÞ

dn þ 1L

R1

dn 0 GAz ðnÞ

L

R1

R1

dn 0 GAz ðnÞ

L

R1

0

0

x

F nþ1 ðxÞ hx  ainþ1 F n ðxÞdx ¼ ¼ nþ1 nþ1

8 n1 > nhx  ai ; n P 1 @F n ðxÞ @hx  ain < 0 ¼ ¼ hx  ai ; n¼1 > @x @x : 0; n¼0

ð39Þ

ð40Þ

0

ðn1Þ2 EIz ðnÞ

dn þ 1L

R1

dn 0 GAz ðnÞ

0

R1

dn 0 GAz ðnÞ

0 L

R1

dn 0 GIt ðnÞ

3 7 7 7 ð45:bÞ 5

shear diðjÞyðzÞ ¼ dbending iðjÞyðzÞ þ diðjÞyðzÞ

ð46Þ

where dbending ¼ iyðzÞ

DqyðzÞ L2 nðn  1Þ DM yðzÞ ðnÞ ¼ DM iyðzÞ ð1  nÞ  DMjyðzÞ n þ 2 h i 0 þ DM pbyðzÞ n  hLn  bi þ DMpjyðzÞ n h i þ DPyðzÞ hLn  ai1  nðL  aÞ DT yðzÞ ðnÞ

dshear iðjÞy ¼

where n = x/L. Assuming elastic behaviour within a load increment, and no coupling of axial and flexural responses at the section level, the increment of the strain energy DW can be written as follows, including the additional shear and torsional deformations, Fig. 12:

dn þ 1L

and dr is a term resulting from loading actions. The equivalent nodal displacements dr of the beam-column element without rigid body modes are computed according to Eq. (43) and their elements can be written as:

With these functions we can determine the element force fields (bending moments and shear forces) as functions of the nodal element forces (DMi, DMj), lateral loads and concentrated bending moments applied along the beam column element. With positive convention as shown in Fig. 12 the internal forces, about each principal axis (y or z), are given by the following expressions:

DMiyðzÞ þ DMjyðzÞ DqyðzÞ Lð2n  1Þ dM yðzÞ ðnÞ ¼ þ ¼ 2 dn L  1 1 La þ DM pbyðzÞ þ DM pjyðzÞ þ DPyðzÞ hLn  ai0  ð41Þ L L L

ðn1Þ2 EIz ðnÞ

0

These functions can be integrated and differentiated as follows:

Z

0

ð45:aÞ

Z DqyðzÞ L3 1 nðn  1Þ2 dn þ DMpbyðzÞ 2 EIyðzÞ ðnÞ 0 Z 1 ð1  nÞbn  hLn  bi0 c  dn þ DM pjyðzÞ EIyðzÞ ðnÞ 0 Z 1 Z 1 nð1  nÞ ð1  nÞ½hLn  ai1  nðL  aÞ  dn þ DPyðzÞ EIyðzÞ ðnÞ EIyðzÞ ðnÞ 0 0

Z DqyðzÞ L 1 ð1  2nÞ DMpbyðzÞ dn  GA ðnÞ L 2 yðzÞ 0 Z 1 D M 1 pjyðzÞ dn   L 0 GAyðzÞ ðnÞ Z 1 Z 1 ½hLn  ai0  La  1 L dn  DP yðzÞ dn  GAyðzÞ ðnÞ 0 GAyðzÞ ðnÞ 0

dbending ¼ jyðzÞ

Z 1 n2 ðn  1Þ dn þ DMpbyðzÞ EIyðzÞ ðnÞ 0 0 Z 1 nbn  hLn  bi0 c n2  dn þ DMpjyðzÞ dn EIyðzÞ ðnÞ EI yðzÞ ðnÞ 0 Z 1 n½hLn  ai1  nðL  aÞ þ DPyðzÞ EIyðzÞ ðnÞ 0

DqyðzÞ L3 2

Z

ð47Þ

ð48Þ

1

ð49Þ

136

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

To produce the deformational-stiffness relation, Eq. (44) is inverted, resulting the following deformational-stiffness equation:

Dsr ¼ kr  Dur  qr

ð50Þ

where the vector qr is the equivalent load vector, whereas kr represents the instantaneous element stiffness matrix of the beam-column element without rigid body modes, determined by matrix inversion of the flexural matrix fr:

" 1

krð66Þ ¼ f r ¼

1

f1

0

0

f2

1

#

" ¼

k1ð33Þ

0ð33Þ

0ð33Þ

k2ð33Þ

# ð51Þ

qr ¼ kr dr The integrals of the flexural coefficients and equivalent load vector are calculated numerically using Eqs. (22)–(24) to express flexural and axial rigidity EIy(z), EA. The cross-sections are located at control points whose number and location depends on the numerical integration scheme. In this work, the Gauss–Lobatto rule for element quadrature is adopted because it has integration points at each ends of the element, where the plastic deformations is important, and hence performs better in detecting yielding. However for members subjected to transversal loads within the member length important plastic deformations could be developed inside the element length where the bending moment has a maximum value. In these situations the integration rule is applied by dividing the interval into two subintervals, with nodes at the element ends and at the section with maximum bending moment (at the section where the shear force has zero value) and use the quadrature rule in each subinterval. In this context the Lobatto integration scheme has another advantage over the Legendre integration scheme, given by the fact that the point corresponding to the left end in one interval is the same as the point corresponding to right end in the next, the cost of evaluating a Lobatto rule is reduced by about one integrand evaluation compared with Legendre rule. The resulting element stiffness matrix is a 6  6 matrix. To include rigid body modes, the stiffness matrix is pre-and post multiplied by a transformation matrix to result in the required 12  12 matrix [31]. 2.2.2. The second-order effects on element tangent stiffness matrix The geometrical nonlinear effects for each element are taken into account in the present analysis, in a beam column approach, by the use of the stability stiffness functions and updating at each load increment the length, axial force and the flexural rigidity about of each principal axes of the element. This way minimizes modelling and solution time, generally only one element is needed per member [16,30]. The effect of axial force on torsional stiffness is ignored in the present formulation. The element force fields are described by the second order transfer matrix as function of the nodal element forces [30]. The equivalent nodal forces are calculated taking into account the second-order effects in a similar way. 2.2.3. The plastic surface requirements If the state of forces at any cross-section along the beam column element equals or exceeds the plastic section capacity (i.e. when the strain in the extreme fibre, attains the ultimate value), the flexural stiffness at the respective location approaches zero or becomes negative when the strain-softening effect for the concrete in compression is taken into account. In order to avoid numerical instabilities, for this cases, the sectional hardening is activated, thus a residual value of the tangent flexural rigidity of the cross EIt section is considered to be EI ¼ 0:001. This effect is reflected in 0 the element tangent stiffness matrix coefficients given by Eqs. (45). It is important to highlight that negative values of tangent flexural rigidities are correlated with the strain softening degree of the concrete in compression and with the magnitude of the

compressive axial force. For steel sections (elastic–perfect plastic) or for composite steel–concrete sections when the strain softening effect is not taken into account the tangent flexural rigidities are always positive at ultimate strength capacity of the cross-section. For the cross-sections when the ultimate strength capacity is attained it is considered that a plastic hinge is developed. Once a plastic hinge is formed, the internal forces must satisfy the condition of plastic strength but the combination of internal forces may change. Once the member forces get to the full plastic surface, the equilibrium condition may be violated as the applied bending moments at the section are greater than the moment capacities, which is not acceptable for maximum strength analysis of members. This reflects the condition of the force point lying outside the failure surface indicated in Fig. 9 and a procedure has to be applied to bring back this force point onto the failure surface. Therefore, when the axial force of a member increases at the following loading step, the incremental force–displacement relationships at the element level has to be modified such that the loading result in motion along the interaction surfaces and the plastic strength surface requirement of the full plastified sections is always satisfied. This approach derives from the plastic theory of work hardening (or stable) material as one in which the work done during incremental loading is positive and associated flow rule can be applied. An important characteristic of concrete that cannot be adequately treated by the classical work-or strain hardening theory of plasticity is the full stress– strain behaviour beyond the peak stress called stress-softening effects. Various types of non-associated flow rules have also been proposed, but the associated flow rule concept is applied here for practical reasons, since the question of strain softening is highly controversial [43] and since there is very little experimental evidence reported in the literature on the direction of flow of the plastic-strain-increment vector under two and three dimensional states of stress. Assuming that at the end of a loading step, at a cross-section, the forces violates the plastic surface, the bending moment is adjusted to be reduced, considering the member loaded by the correction bending moment DMp, where DMp is the change in the plastic moment capacity at the respective cross section as axial force changes (Fig. 9). The correction plastic bending moments are determined as follows. Let us consider that at the current loading step in a specified location the forces ðN  ; M y ; M z Þ exceed the plastic surface, this means that for that loads the maximum strain in most compressed or tensioned fibre of the cross-section exceeds the assumed ultimate strain. For given value of axial force (N⁄) the iterative procedure described at Section (2.1.3) is applied in order to determine the ultimate bending moments ðM yp ; Mzp Þ for a given bending moMz ment ratio tan a ¼ M  . Next the correction plastic bending moment y is determined as:

DM yp ¼ M y  M yp DM zp ¼ M z  M zp

ð52Þ

The additional incremental moments are applied as loadings on the element and are considered through equivalent nodal forces as already described in Section (2.2.1). Since the force-point movement remains on the plastic strength surface of a member, the plastic strength surface requirement of a section is not violated by the change of member forces after the full plastic strength of cross-section is reached. 2.2.4. Geometry updating and analysis algorithm In order to trace the equilibrium path, for proportionally and non-proportionally applied loads, the proposed model has been implemented in a simple incremental matrix structural-analysis program. The simple Euler stepping algorithm is used in conjunction with constant work-load increments and the natural deforma-

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

137

tion approach. This analysis is simple, reliable and is not sensitive to convergence failures that can occur in incremental-iterative schemes, and can also give the full nonlinear load-deformation response including the ultimate load and post-critical response. The incremental change in the displacements can be written as the solution of:

tions at each load increment. The natural deformation approach (NDA) in conjunction with the geometrical ‘‘rigid body qualified’’ stiffness matrix [44] is adopted for the element force recovery and the web plane vector approach is effectively used to update the frame element coordinates [30].

KiT  DUi ¼ Dki  DFil þ DFip

3. Computational examples

ð53Þ

where within a particular load cycle, i, KiT the tangent stiffness matrix; DUi the incremental displacement vector, DFil the incremental nodal force vector including the member loads; DFip the additional self-equilibrating nodal force vector. Only those elements that have moments in excess of the strength capacities will contribute to DFip . The incremental load factor Dk is computed so as to keep constant the incremental work DW, performed by the applied external loads, at each load cycle:

Dki ¼

DW  DFip  DUi DFil  DUi

ð54Þ

The accuracy of the solution is controlled via multiple analyses based on convergence of system response.Using an updated Lagrangian formulation (UL) the nonlinear geometrical effects are considered updating the element forces and geometry configura-

Based on the analysis algorithm just described, a computer program, NEFCAD, has been developed to study the combined effects of material and geometric nonlinear behaviour on the load-versus-deflection response for spatial composite steel–concrete framed structures. It combines the structural analysis routine with a graphic routine to display the final results. The computational engine was written in Compaq Visual Fortran. The graphic interface was created using Microsoft Visual Basic 6. Dynamic Link Libraries (DLL) are used to communicate between the interface and engine. The accuracy of the analytic procedure and the computer program developed here, has been evaluated using several benchmark problems analyzed previously by other researchers using independent finite element solutions. The selected problems consist of simply column, proposed here, simply-supported beams [45], continuous beams [46], plane frames [18,23,47] whose nonlinear response is dominated by spreading plasticity effects in individual members,

Fig. 13. Composite steel–concrete column.

Fig. 14. Moment–curvature analysis for different values of compressive axial loads with and without residual stress effects: (a) weak axis bending; and (b) strong axis bending.

138

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

and two large-scale space frames [18,23] in which both geometric and material nonlinear effects contributed to the failure. In the present approach, one element has been used to model each column and beam in all computational examples. 3.1. Example 1: Composite steel–concrete column Two sets of nonlinear analyses have been conducted for axially compressed element, proposed in this work, as shown in Fig. 13, to verify the performance and the efficiency of the proposed method in capturing spread of plasticity, second-order effects and also the effects of residual stresses over ultimate strength capacity, deformability and flexural rigidity of composite steel–concrete elements. In the first set of tests inelastic cross-section analyses have been conducted. In this respect parametric studies concerning moment– curvature response, moment-tangent flexural rigidity evaluation under condition of constant axial load and ultimate strength capacity evaluation for different loading orientation have been performed. The automatic drawing of the moment -curvature diagram and subsequently the bending moment-tangent flexural rigidity

diagram, is accomplished here by applying the iterative procedure described at Section 2.1.2, scaling the bending moments through the load factor under constant axial force and then an automatic step length adaptation scheme, such as arc length strategy, for the loading factor is applied [34]. The cross-section consists of a concrete core and a symmetrically placed USA wide flange section W12  120 (Fig. 13). Characteristic strength for concrete in compression is fc = 20 MPa and the stress–strain curve which consists of a parabolic and linear part was used in the calculation (Eq. (1)), with crushing strain e0 = 0.002, ultimate strain ecu = 0.0035. The strain softening of concrete is taken into account with c = 0.15. The Young modulus for encased structural steel and reinforcements is 200 GPa. A bi-linear elasto-perfect plastic stress–strain relationship for the reinforcement bars and structural steel, both in tension and in compression, is assumed with the yield strength of fyr = 400 MPa and fys = 300 MPa respectively. Fig. 14 shows the comparative bending moment–curvature diagrams of cross-section, considering uniaxial bending about X-axis (strong) and Y-axis (weak) respectively under several compressive axial loads with and without residual stress effect. The effect of residual stresses is taken into account considering two types of

Fig. 15. Variation of tangent flexural rigidity: (a) weak axis; and (b) strong axis.

Fig. 16. Variation of tangent flexural rigidity: biaxial bending-(a) strong axis; and (b) weak axis.

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

Fig. 17. Interaction diagrams for different values of bending moment’ s ratio.

Fig. 18. Lateral load–displacement curves (a) Minor axis bending; and (b) Major axis bending.

139

140

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

residual stress patterns for structural steel element, EC3 and AISCLRFD, respectively. The comparison reveals that the influences of residual stresses play an important role on both carrying capacity and inelastic behaviour, during the loading process, on composite cross section and this effect becomes more effective when EC3 distribution is taken into account and for bending about weak axis. The moment-tangent flexural rigidity curves, for weak and strong axis bending, under condition of constant axial load, shown in Fig. 15 demonstrate that the stiffness degradation of composite section, in presence of residual stresses, is evident and more pronounced especially for weak axis bending, in the case when EC3 distribution for residual stresses is assumed. This section was also analysed for different compressive axial loads and with moments applied about an axis at a = 30° to the strong axis. Tangent flexural rigidities for weak and strong axis bending are shown in Fig. 16. The results indicate similar behaviour as already described in the cases of uniaxial bending. Fig. 17

presents the interaction diagrams for different values of bending moment ratio (a = 0-strong axis bending, a = 30°; a = 60° and a = 90°-weak axis bending) and with and without the presence of residual stresses. As it can be seen the presence of residual stresses indicates lower capacity of cross section, even for low compressive axial load, and this effect is more pronounced when bending take place about weak axis of the cross-section and in the case of EC3 distribution for residual stresses. Fig. 18 shows the lateral applied load–displacement curves for the column considering uniaxial bending about either weak or strong axis respectively and for different values of compressive axial load N. The column is subjected to non-proportional loads. The compressive axial force N is first applied and kept constant whereas the lateral loads, either Hx (bending on weak axis) or Hy (bending on strong axis) are then applied and progressively increased. For the encased structural steel EC3 distribution for residual stresses has been considered. As it can be seen when the residual stresses effect is taken into account

Fig. 19. Simply supported composite beam subjected to two concentrated loads.

Fig. 20. Simply supported composite beam subjected to concentrated load at midspan length.

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

141

Fig. 21. Two span continuous composite beam.

tion of the column is more pronounced in the presence of the residual stresses especially in the case of bending about weak axis. One important conclusion from this case study can be drawn that the influence of residual stresses, for encased steel section, on the carrying capacity and inelastic behaviour during the loading process is important, especially for bending about weak axis, and neglecting this effect may overestimate the structural stiffness and ultimate capacity of composite steel–concrete elements.

3.2. Example 2: Simply supported composite beams Fig. 22. Composite beam frame.

the ultimate load factor is reduced and also the stiffness degrada-

Two simply supported composite beams subject to saging moments, tested by Nie and Cai [45], are presented here, for numerical verification of the proposed approach. The composite beams have been also analysed numerically by Nie et al. [14] using an advanced mixed finite-element approach combining the fibred beam and

Fig. 23. Load–displacement curves for composite beam frame.

142

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

Fig. 24. Composite portal frame: (a) the geometry and loading arrangement; and (b) cross-sections details.

Fig. 25. Load–displacement curves for composite portal frame.

Fig. 26. Composite steel–concrete portal frame.

layered shell elements. The geometry, material and section properties of the simple composite beams are depicted in Figs. 19 and 20. The cubic compressive strength of the concrete in compression

is fcu = 27.7 N/mm2, the yield stress of the steel beams is fsy1 = 310 N/mm2 , Young’s modulus is Es = 20 000 MPa, Esh1 = 0, esy1 = 2.5% (Fig. 4) and the strain hardening modulus is

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

143

Fig. 27. Load–displacement curves of composite portal frame.

Fig. 28. Tangent flexural rigidity (1  EIt/EI0) distribution at ultimate load factor: (a) about weak axis; and (b) about strong axis.

Esh2 = 100 MPa. The behaviour of the reinforced bars, both in tension and compression, are modelled according with the model described in Section 2.1.1 with the yield strength fyr = 290 N/mm2. The tensile strength of the concrete is taken into

account according with the model described in Fig. 3, with a1 = 1 and a2 = 0.75. According to the present analysis, the midspan deflection is plotted on Figs. 19 and 20 where the experimental results from

144

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

[45] and fibre finite element analyses [14] are also compared. As it can be seen the behaviour of composite beams predicted by the present analysis is consistent with that of experimental tests and is in close agreement with that of mixed-finite element analysis. It is important to note that when the contribution of the rebars is neglected the proposed numerical approach gives conservative predictions compared to the experimental results. This fact has been also observed by Nie et al. [14] based on their advanced finite element simulations. However, because in the present model, the

shear slip effect between concrete slab and steel beam is not taken into account, the load deflection curves predicted by the proposed model indicates a slightly stiffened behaviour of the beams as compared with the experimental tests. 3.3. Example 3. Continous composite beam The two span continous composite beam with a loading arrangement as shown in Fig. 21 has been tested by Slutter and

Fig. 29. Orbison’s six story rigid space frame: (a) plan view; (b) perspective view; and (c) member cross-sections.

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

Driscoll [46] and numerically analysed by the Liew et al. [23] and Nie et al. [14]. The geometry and section properties of the continous composite beam are depicted in Fig. 21. The cylinder compressive strength of the concrete in compression is fc0 ¼ 16 MPa, the yield stress of the structural steel is fsy1 = 252.4 MPa, Young’s modulus is Es = 20 000 MPa, Esh1 = 0, esy1 = 2.5% (Fig. 4) and the strain hardening modulus is Esh2 = 100 MPa. As it can be seen in Fig. 21 the behaviour of continuous composite beam predicted by the present analysis is in close agreement with that of experimental test [46] and mixed-finite element analysis [14].

145

3.4. Example 4. Composite beam frame The fully connected composite frame tested by Bursi and Gramola [47] is depicted in Fig. 22. This frame has been analysed numerically by Nie et al. [14] using an advanced mixed fibre finite-element approach. The geometry, section properties and loading arrangement of the composite frame are depicted in Fig. 22. The reinforcement in the concrete slab is doubled near the column. The cylinder compressive strength of the concrete in compression is fc0 ¼ 39 MPa, the yield stress of the steel beams is fsy1 = 300 MPa, Young’s modulus is Es = 20 000 MPa, Esh1 = 0, esy1 = 2.5% (Fig. 4) and

Fig. 30. Load–deflection curves at node A in X and Y directions.

Fig. 31. Tangent flexural rigidity (1  EIt/EI0) distribution along the member’s length at ultimate load factor: (a) composite frame; and (b) pure steel frame.

146

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

Fig. 32. Details of tangent flexural rigidity (1  EIt/EI0) distribution for composite frame: (a) third floor beams (strong axis); and (b) first level columns (strong and weak axis).

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

147

Fig. 33. Details of tangent flexural rigidity (1  EIt/EI0) distribution for pure steel frame: (a) third floor beams (strong axis); and (b) first level columns (strong and weak axis).

148

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

the strain hardening modulus is Esh2 = 100 MPa. The behaviour of the reinforced bars, both in tension and compression, are modelled according with the model described in Section 2.1.1 with the yield strength fyr = 480 MPa. The tensile strength of the concrete is taken into account according with the model described in Fig. 3, with a1 = 1 and a2 = 0.75. As it can be seen in Fig. 23 the results obtained by the present analysis correlate well with that of experimental test and mixed-finite element analysis. It is important to note that in absence of the strain hardening for the structural steel components (elastic–perfect plastic, Esh1 = 0, Esh2 = 0) the proposed numerical approach gives conservative predictions compared to the experimental results (Fig. 23). 3.5. Example 5. Composite portal frame The steel–concrete portal frame shown in Fig. 24 consists of a composite steel–concrete beam with full composite action and two wide flange W12  50 steel profiles completely encased in concrete. This frame has been analysed by Liew et al. [23] considering the steel columns and composite beam and further modified by Iu et al. [18] in which the columns have been completely encased in concrete as shown in Fig. 24. A refined plastic hinge model for composite elements was employed in [18]. In this study the behaviour of a steel frame with or without composite action for beam and columns is investigated. The geometry, section properties and loading arrangement of the fully composite frame are depicted in Fig. 24. The frame is subjected to both vertical and lateral loading, proportionally applied. The yield strength of all steel members is 248.2 MPa and Young’s modulus is E = 20 000 MPa. The stran hardening effect of the steel is ignored. Characteristic strength for concrete in compression is fc = 16 MPa and the stress–strain curve which consists of a parabolic and linear part was used in the calculation, with crushing strain e0 = 0.002 and ultimate strain ecu = 0.0035. The tensile strength of the concrete is ignored in the proposed analysis. The inelastic behaviour represented by the load–deflections curves calculated by the proposed approach and those given in [18,23] is compared in Fig. 25. As it can be seen the results agree closely.

(elastic–perfect plastic, Esh1 = 0) the ultimate load factor is drastically reduced and a clear plastic mechanism is revealed. For this case the Ref. [16] does not presents comparative results. EIt Fig. 28 shows the variation of the flexural rigidities (1  EI ) 0 along the member lengths, at the ultimate load factor. As it can be seen plastic deformations are concentrated at the midspan and right-side of the beam and at the column’s base. The failure of the frame is due to progressive yielding of the composite beam and columns leading to significant stiffness degradation and sideway deflection. Running the present computer program on a laptop computer with 2 GHz processor, the present analysis was performed in only 4 s, with over 100 load cycles, whereas the load–displacement curve obtained by Abaqus and reported in [16] requires almost 48 min running on a similar computer. This result demonstrates the computational efficiency and time saving of the proposed approach.

3.6. Example 6: Steel portal frame with composite beam The steel–concrete portal frame shown in Fig. 26 consists of a composite steel–concrete beam with full composite action and two wide flange W12  50 steel columns. The geometric configuration and loading arrangement is shown in Fig. 26. This frame has been analyzed by the Ngo-Huu and Kim [16] using a fibre plastic hinge approach and an advanced fine-grained finite element numerical model developed in the Abaqus software. The frame is subjected to the combined action of concentrated load (P = 150 kN) at midspan of the beam and a lateral load as shown in Fig. 26. Both loads are considered to be applied proportionally until the collapse of the structure. The yield strength of all steel members is 252.4 MPa, Young’s modulus is E = 20 000 MPa and the strain hardening modulus is Esh1 = 6000 MPa. The beam section consists of a W12  27 steel section and a concrete slab 102 mm  1219 mm. Characteristic strength for concrete in compression is fc = 16 MPa and the stress–strain curve which consists of a parabolic and linear part was used in the calculation, with crushing strain e0 = 0.002 and ultimate strain ecu = 0.00 806. The strain softening of concrete is not taken into account (c = 0). Fig. 27 shows the comparative load factor-lateral displacement curves obtained using the proposed approach, with and without strain hardening effect, and the results retrieved from [16]. As it can be seen the proposed model agrees fairly well with references, when strain hardening effect is taken into account. However in absence of the strain hardening for the structural steel components

Fig. 34. Twenty story composite space frame: (a) plan view; and (b) perspective view.

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

3.7. Example 7: Six-story composite space frame The Orbison’s six story rigid space frame, studied previously by other researchers, has been modified and included in this verification study (Fig. 29a–c). In order to demonstrate the effectiveness of the numerical procedure developed here, this example converts the pure steel frame to a composite frame [18]. The steel beams in the frame rigidly support the concrete floor slab, and the steel columns are partially encased in the concrete (Fig. 29c) [18]. The yield strength of all steel members is 250 MPa, Young’s modulus is E = 20 700 MPa and shear modulus G = 79 293 MPa. The compressive yield stress of concrete is fc = 16 MPa and the stress–strain curve which consists of a parabolic and linear part was used in the calculation, with crushing strain e0 = 0.002 and ultimate strain ecu = 0.0035. The strain softening of concrete is not taken into account (c = 0). The frame is subjected to the combined action of

149

gravity and lateral loads acting in the Y-direction. Uniform floor pressure is 9.6 kN/m2 (simulated here as uniform distributed loads on each beam) and the wind loads are simulated by point loads of 53.376 kN in the Y direction at every beam-column joints. The inelastic behaviour represented by the load–deflections curves at node A in X and Y directions calculated by the proposed approach and those given in the literature is compared in Fig. 30. A refined plastic hinge model for composite elements was employed in [18], whereas a fibre finite element approach for pure steel elements was included in [17]. Good correlation of the numerical results can be observed. However in the case of full composite frame (beams and columns) the proposed approach indicates the ultimate load factor 1.23 which is slightly higher than ultimate load factor obtained by Iu et al. (1.17) using a refined plastic hinge approach [18]. Some differences can be observed also over the lateral stiffness of the frame. In the analysis of frame while considering

Fig. 35. Load–deflection curves at node A.

Fig. 36. Load–deflection curves at node B.

150

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

the composite beams and columns the ultimate load factor is slightly higher than in the case of steel frame with composite beams only, in the last case the ultimate load factor resulting 1.19. As it can be seen the lateral stiffness of the frame slightly increases for full composite case (beams and columns) as compared with the case in which the composite action is taken into account only for beams. Fig. 31 shows the variation of flexural rigidities along the member lengths, at the ultimate load factor, for the pure steel frame and for the composite steel–concrete frame, respectively. Details of the flexural rigidity distribution (1  EIt/EI0) for the columns (both major and minor axes) at the first level and for the beams (about major axis) at the third level are depicted in Fig. 32 (composite frame) and Fig. 33 (pure steel frame), respectively. It can be clearly observed the flexural stiffness degradation for beams and columns in these cases. Running the present computer program on a laptop computer with 2 GHz processor, the present analysis was performed in almost 20 s, considering 5 integration points along the member’s length.

the columns (both major and minor axes) at the first level and for the beams (about major axis) at the third level are depicted in Figs. 38 and 39 respectively. It can be clearly observed the flexural stiffness degradation for beams and columns in these cases. As it can be seen in Fig. 39 the stiffness degradation is more accentuated for the beams in the y direction and plastic deformations are concentrated significantly toward the end region of the beams. The same observation has been reported by Liew et al. [23] when a distributed plasticity approach was involved to study the behaviour of the beams of this frame. Running the present computer program on a laptop computer with 2 GHz processor, the present analysis, for the case when composite columns and composite beams is considered, was performed in almost 2 min, despite the fact that the analysis has been launched within the graphical user interface. This result proves the high computational efficiency of the proposed approach.

4. Conclusions 3.8. Example 8: Twenty-story composite space frame Twenty-story steel space frame with dimensions and properties shown in Fig. 34 has been studied previously by Liew et al. [23], Jiang et al. [17] and Ngo-Huu et al. [48]. Liew et al. [23] analyzed also this frame considering composite beams with full shear connection. A50 steel is used for all steel sections, the yield strength of steel is assumed to be fy = 344.8 N/mm2, Young modulus E = 2  105 N/mm2 and elastic–perfect plastic constitutive relationship is considered. Overall slab depth is assumed to be 127 mm, and the compressive yield stress of concrete is fc = 27.6 N/mm2. The concrete slab is considered within the range of effective flange width (1/4 span length for interior beams and 1/8 span length for exterior beams). The frame is analyzed for the combination of gravity loads = 4.8 kN/m2 (simulated here as uniform distributed loads on each beam) and wind loads = 0.96 kN/m2, acting in the y-direction (simulated here as nodal loads). In Liew’s et al. [23] study, rigid floor diaphragm action is assumed in the global analysis. Liew employed one plastic-hinge beam-column element to model each steel column and four elements for each beam. The inelastic behaviour of beams is taken into account considering M–U relationship. The limit load of the frame is reached at the load ratio of 1.031, for bare steel frame whereas a load factor of 1.338 is obtained when composite steel– concrete beams are considered. In Ngo-Huu’s et al. [48] study, the fibre plastic hinge concept has been used to predict the second-order inelastic behaviour of the pure steel frame, and the inelastic limit point reported is 1.003. In present analysis one element with five integration points has been used to model each column and beam and the stress–strain curve for concrete in compression which consists of a parabolic and linear part was used in the calculation, with crushing strain e0 = 0.002, ultimate strain ecu = 0.0035. The strain softening of concrete is not taken into account (c = 0). The comparative load–deflection curves of nodes A and B at the top of the frame calculated by the previous researchers and the results obtained from the new method are shown in Figs. 35 and 36. As it can be seen the results obtained by the proposed model agree fairly well with the references. In the analysis of frame while considering the composite beams and columns the ultimate load factor is practically the same as in the case of steel frame with composite beams only, but the lateral stiffness of the frame is slightly increased for full composite case (beams and columns) as compared with the case in which the composite action is taken into account only for beams. Fig. 37 shows the variation of flexural rigidities (1  EIt/EI0) along the member lengths, at the ultimate load factor, for the steel–concrete frame with composite beams. Details of the flexural rigidity distribution (1  EIt/EI0) for

A reliable and robust nonlinear inelastic analysis method for composite steel–concrete space frames has been developed. The proposed model is based on the most refined type of second order inelastic analysis, the plastic zone analysis. Gradual yielding is described through basic equilibrium, compatibility and material non-

Fig. 37. Tangent flexural rigidity (1  EIt/EI0) distribution along the member’s length at ultimate load factor.

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

151

Fig. 38. Details of tangent flexural rigidity (1  EIt/EI0) distribution for first level columns (strong and weak axis).

Fig. 39. Details of tangent flexural rigidity (1  EIt/EI0) distribution for third floor beams (strong axis).

linear constitutive equations and path integral approach is applied to numerical integration of the cross-sectional nonlinear characteristics. The proposed analysis can practically account for major factors influencing nonlinear composite steel–concrete space frame behaviour: gradual and distributed yielding associated with biaxial

bending and axial force, stress–strain relationships for concrete in compression and tension, residual stress distribution in encased steel elements, shear deformations, local and global second order effects with computational efficiency, and the necessary degree of accuracy, usually only one element per member is necessary to analyze. Lateral loads acting along the member length can be

152

C.G. Chiorean / Engineering Structures 57 (2013) 125–152

directly input into the analysis without the need to divide a member into several elements, leading to a consistency in the linear and nonlinear structural models. The model has been implemented in a simple incremental matrix structural analysis program and allows proportionally and non-proportionally loading, and has been verified by comparing the predicted results with the established results available from the literature. The studies show that the proposed analysis compares very well to experimental and finite fibre element solution with much less computational effort. The proposed software is presented as an efficient, reliable tool ready to be implemented into design practice for advanced analysis and pushover analysis of spatial composite steel–concrete frame structures. Future work is envisaged, considering the extension of the proposed method for nonlinear inelastic analysis of 3D composite steel–concrete frameworks with partial shear connection of composite beams and semi-rigid behaviour of composite connections. Acknowledgements The author gratefully acknowledges the support from Romanian National Authority for Scientific Research (ANCS and CNCSIS Grant PNII-IDEI No. 193/2008) for this study. References [1] Li GQ, Li JJ. Advanced analysis and design of steel frames. Wiley; 2007. [2] Ziemian RD. Guide to stability design criteria for metal structures. 6th ed. Wiley; 2010. [3] Ayoub A. Analysis of composite frame structures with mixed elements-state of the art. Struct Eng Mech 2012;41:157–81. [4] Spacone E, El-Tawil S. Nonlinear analysis of steel-concrete composite structures: state of the art. J Struct Eng ASCE 2004;130:159–68. [5] Nguyen QH, Hjiaj M, Uy B, Guezouli S. Analysis of composite beams in the hogging moment regions using a mixed finite element formulation. J Constr Steel Res 2009;65:737–48. [6] Valipour HR, Bradford MA. A steel-concrete composite beam element with material nonlinearities and partial shear interaction. Finite Elem Anal Des 2009;45:966–72. [7] Hjiaj M, Battini JM, Nguyen QH. Large displacement analysis of shear deformable composite beams with interlayer slips. Int J Non-Linear Mech 2012;47:895–904. [8] Salari MR, Spacone E. Analysis of steel-concrete composite frames with bondslip. J Struct Eng ASCE 2001;127:1243–50. [9] Zona A, Barbato M, Conte JP. Nonlinear seismic response analysis of steelconcrete composite frames. J Struct Eng ASCE 2008;134:986–97. [10] Elghazouli AY, Treadway J. Inelastic behaviour of composite members under combined bending and axial loading. J Constr Steel Res 2008;64:1008–19. [11] Pi YL, Bradford B, Uy B. Second order nonlinear analysis of composite steelconcrete members. I: theory. J Struct Eng ASCE 2006;132:751–61. [12] Fang LX, Chan SL, Wong YL. Strength analysis of semi-rigid steel-concrete composite frames. J Constr Steel Res 1999;52:269–91. [13] Bursi OS, Sun FF, Postal S. Non-linear analysis of steel concrete composite frames with full and partial shear connection subjected to seismic loads. J Constr Steel Res 2005;61:67–92. [14] Nie J, Tao M, Cai CS, Chen GE. Modeling and investigation of elasto-plastic behavior of steel-concrete composite frame systems. J Construct Steel Res 2011;67:1973–84. [15] Ellobody E, Young B. Numerical simulation of concrete encased steel composite columns. J Construct Steel Res 2011;67:211–22. [16] Ngo-Huu C, Kim SE. Practical nonlinear analysis of steel-concrete composite frames using fiber-hinge method. J Construct Steel Res 2012;74:90–7. [17] Jiang XM, Chen H, Liew JYR. Spread of plasticity analysis of three-dimensional steel frames. J Construct Steel Res 2002;58:193–212. [18] Iu CK, Bradford MA, Chen WF. Second-order inelastic analysis of composite framed structures based on the refined plastic hinge method. Eng Struct 2009;31:799–813.

[19] Iu CK. Inelastic finite element analysis of composite beams on the basis of the plastic hinge approach. Eng Struct 2008;30:2912–22. [20] Liu SW, Lui YP, Chan SL. Advanced analysis of hybrid steel and concrete frames. Part I: cross section analysis technique and second-order analysis. J Construct Steel Res 2012;70:326–36. [21] El-Tawil S, Deierlein GG. Nonllinear analysis of mixed steel-concrete frames. I: element formulation. J Struct Eng ASCE 2001;127:647–55. [22] Ashraf A, Filippou FC. Mixed formulation of nonlinear steel-concrete composite beam element. J Struct Eng ASCE 2000;126:371–81. [23] Liew JYR, Chen H, Shanmugam NE. Inelastic analysis of steel frames with composite beams. J Struct Eng ASCE 2001;127:194–202. [24] Nukala PKV, White DW. A mixed finite element for three dimensional nonlinear analysis of steel frames. Comp Methods Appl Mech Eng 2004;193:2507–45. [25] Izzudin BA, Siyam AAFM, Lloyd-Smith D. An efficient beam-column formulation for 3D RC frames. Comput Struct 2002;80. [26] Neuenhofer A, Filippou FC. Evaluation of nonlinear frame finite-element models. J Struct Eng ASCE 1997;123:958–66. [27] Neuenhofer A, Filippou FC. A geometrically nonlinear flexibility-based frame finite element. J Struct Eng ASCE 1998;124:704–11. [28] Sivaselvan MV, Reinhorn AM. Collapse analysis: large inelastic deformations analysis of planar frames. J Struct Eng ASCE 2002;128:1575–83. [29] Valipour HR, Foster SJ. A total secant flexibility-based formulation for frame elements with physical and geometrical nonlinearities. Finite Elem Anal Des 2010;46:288–97. [30] Chiorean CG, Barsan GM. Large deflection distributed plasticity analysis of 3D steel frameworks. Comput Struct 2005;83:1555–71. [31] Chiorean CG. A computer program for nonlinear inelastic analysis of 3D semirigid steel frameworks. Eng Struct 2009;31:3016–33. [32] Attala MR, Deierlein GG, McGuire W. Spread of plasticity: quasi-plastic hinge approach. J Struct Eng ASCE 1994;120:2451–73. [33] Marmo F, Rosati L. An improved flexibility-based nonlinear frame element endowed with the fiber-free formulation. In: European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012). Viena, Austria; 2012. p. 1–17. [34] Chiorean CG. Computerised interaction diagrams and moment capacity contours for composite steel-concrete cross-sections. Eng Struct 2010;32:3734–57. [35] Papanikolau VK. Analysis of arbitrary composite sections in biaxial bending and axial load. Comput Struct 2012;98-99:33–54. [36] 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. [37] Marmo F, Serpieri R, Rosati L. Ultimate strength analysis of prestressed reinforced concrete sections under axial force and biaxial bending moments. Comput Struct 2011;89:91–108. [38] Skrabek BW, Mirza SA. Strength reliability of short and slender composite steel-concrete columns. In: Civil Engineering Report Senes. No. CE-90-1. Lakehead University, Thunder Bay, Ontario; 1990. 323 p. [39] Virdi KS, Dowling PJ. The ultimate strength of composite columns in biaxial bending. In: Procedures institution of civil engineers (London), vol. 56(May). 1973. p. 251–72. [40] Vecchio FJ, Collins MP. The modified compression field theory for reinforced concrete elements subjected to shear. ACI J 1986;83:219–31. [41] Mansour M, Lee JY, Hsu TTC. Cyclic stress-strain curves of concrete and steel bars in membrane elements. J Struct Eng ASCE 2001;127:1402–11. [42] Yavari A, Sarkani S, Moyer ET. On applications of generalized functions to beam bending problems. Int J Solids Struct 2000;37:5675–705. [43] Zhao XM, Wu YF, Leung AYT. Analyses of plastic hinge regions in reinforced concrete beams under monotonic loading. Eng Struct 2012;34:466–82. [44] Yang YB, Yau JD, Leu LJ. Recent developments in geometrically nonlinear and postbuckling analysis of framed structures. Appl. Mech Rev 2003;56:431–49. [45] Nie JG, Cai CS. Steel concrete composite beams considering shear slip effects. J Struct Eng ASCE 2003;129(4):495–506. [46] Slutter RG, Driscoll GC. Flexural strength of steel-concrete composite beams. J Struct Div ASCE 1965;91:71–99. [47] Bursi OS, Gramola G. Behavior of composite substructures with full and partial shear connection under quasi-static cyclic and pseudodynamic displacements. Mater Struct Paris 2000;33:154–63. [48] Ngo-Huu C, Kim SE, Oh JR. Nonlinear analysis of steel frames using fiber plastic hinge concept. Eng Struct 2007;29:649–57.