Second-order flexibility-based model for nonlinear inelastic analysis of composite steel-concrete frameworks with partial composite action and semi-rigid connections

This paper presents an efficient computer method for large deflection distributed plasticity analysis of 3D semi-rigid composite steel-concrete frameworks. A novel second-order inelastic flexibility-based element has been developed by combining the Maxwell-Mohr rule and the second-order force based functions for computation of the generalized displacements. The proposed model allows explicit and efficient modeling of the combined effects of nonlinear geometrical effects, gradual spread-of-plasticity, partial shear connection of composite beams, finite-size joints and joint flexibility by using only one 2-noded beam-column element per physical member. For composite beams, based on elasto-plastic cross-sectional analyses the model is able to take into account the effects of partial composite action between the concrete slab and the steel beam. At the crosssectional level the proposed method addresses computational efficiency through the use of path integral approach to numerical integration of the cross-sectional nonlinear characteristics and residual stresses, enabling in this way the accurate geometrical specifications and precise modeling of cross-sections. The proposed nonlinear analysis formulation has been implemented in a general nonlinear static purpose computer program, NEFCAD. Several computational examples are given to validate the accuracy and efficiency of the proposed method.


Introduction
In recent years, have witnessed significant advances in nonlinear inelastic analysis methods for composite steel-concrete beams and framed structures and integrate them into the new and more rational advanced analysis and design procedures [1,2].Currently the available tools for such analysis are general purpose FE programs that require extensive calibration and mesh generation studies and still possess huge demands on the most powerful of available computers and represents unpractical tasks for structural engineer.The present work attempts to develop accurate yet computational efficient tools for the nonlinear inelastic analysis of partially connected composite steel-concrete frameworks fulfilling the practical and advanced analysis requirements.

Mathematical formulation
In this paper, the following general assumptions are adopted in the formulation of analytical model for inelastic analysis of composite beams: (1) Plane sections remain plane for entire cross-section after flexural deformation; throughout the depth of the crosssection, the strain distribution is linear, but a discontinuity exists at the concrete slab-steel beam interface due to slip, frictional effects and uplift are neglected, the interface slip is small; (2) The vertical displacement and the curvature of the different subcomponents (concrete slab and structural steel) are assumed to be the same; (3) Discretely located interlayer connectors with uniform spacing are regarded as continuous and ductile with an nonlinear elastic-plastic behavior.In the formulation of the inelastic behavior of composite columns full composite action between concrete matrix and steel profile is assumed.Small strains but large displacements and rotations can be considered.Transverse shear deformations, associated to the transverse shear forces are neglected in the plastic constitutive relationships.The model suggested by the CEB-FIB Model Code 90, is adopted in the present paper to model the concrete under compression and tension.Multilinear elastic-plastic stress-strain relationships, both in tension and in compression, are assumed for the structural steel and the conventional reinforcing bars.Gradual yielding of cross-sections and elasto-plastic tangent flexural and axial rigidities, for composite steelconcrete columns with arbitrary shape, when full shear connection is assumed, is described through basic equilibrium, compatibility and material nonlinear constitutive equations following the procedure described in [2] where the discussions concerning the numerical integration of stresses and stiffnesses over cross sections and inclusion of residual stresses for structural steel are also addressed.Therefore, this paper is focused on partial composite action effect over inelastic response of beam crosssections.As will be briefly described in the following sections, the incremental forcedisplacement relationships at the element level are derived by applying the Maxwell-Mohr rule for computation of generalized displacements in the second-order geometrically nonlinear analysis and using an updated Lagrangian formulation the nonlinear global geometrical effects are considered updating the element forces and geometry configurations at each load increment [3].

Elasto-plastic analysis of beam crosssections
The elastic and inelastic behavior of steelconcrete composite beams is quite complex because the shear connectors generally permit the development of only partial composite action between the individual components of the member, and their analysis requires the consideration of the interlayer slip between the subcomponents.Usually, for a given composite beam, the full shear connection is defined as the least number of shear connectors (n f ) such that the bending resistance of the beam would not be affected if more shear connectors are provided, whereas partial shear connection occurs when the number of connectors (n) used in a beam is lower than (n f ).In order to analyze and design the composite beams, simplified methods are very useful and such methods are proposed in international literature and in some codes.For instance in the Eurocode 4 the concept of the degree of shear connection η = n/n f is used and the ultimate bending strength capacity of cross-section is evaluated by simple equilibrium of stresses with a prescribed compressive axial force in the concrete slab N c = nP sc , where n<n f represents the number of shear connectors and P sc is the connector strength.The composite beam cross-section considered here consists of a concrete solid slab connected to a steel beam as presented in Fig 1.Let us consider the cross-section subjected to the action of the external bending moment (M), and axial force (N) as shown in Fig. 1a.The resultant strain distribution corresponding to the curvature  and the axial strain u can be expressed, in full composite action, throughout the depth of cross-section in a linear form as: The equilibrium is satisfied when the external forces (N, M) are equal to the internal ones.These conditions can be represented mathematically in terms of the nonlinear system of equations as in Eq. ( 2).The surface integrals in Eq. ( 2) are extended over steel (A s ) and concrete areas (A c ), the contribution of slab a.
b. (2) The above system can be solved numerically using the Newton iterative method and taking into account the fact that the stresses are implicit functions of the axial strain and curvature through the resultant strain distribution given by Eq. ( 1).In this way, for given bending moment 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.In this way the internal axial force in concrete ( cf N ) in which the contribution of the conventional steel reinforcements is included, can be evaluated [4].Let us consider now the cross-section, in partial composite action, subjected to the action of the external bending moment (M), and axial force (N) as shown in Fig. 1b.Under the above assumptions, the resultant strain distribution, corresponding to the curvature  and the axial strains u c and u s evaluated at the centroid of concrete slab and structural steel respectively, can be expressed in a linear form as: where ε c and ε s represents the strains in concrete slab and steel beam respectively and r represents the distance from the central axis of the concrete slab to that of the steel beam.The equilibrium is satisfied when the external forces are equal to the internal ones: (4) in which u c , u s and  represent the unknowns.In order to solve the above nonlinear system, the internal axial force in the concrete slab (N c int ), under partial composite action, is assumed to be a fraction of the axial force in the concrete slab under full composite action (N cf ) and the amount is defined by a function of the degree of composite action γ eff as in Fig. 1b: ( 5) where N cf represents the internal axial force in the concrete slab of the cross-section subjected to the same external bending moment M and axial force N but under the assumption of full composite action between the steel beam and concrete slab, f(γ eff ) represents a function of the effective degree of composite action [4].Thus, the nonlinear system (5) becomes: (6) The above system can be solved numerically using, the Newton iterative method.The incremental relationships between incremental efforts and incremental deformations can be expressed as: where the coefficients of the tangent stiffness matrix k ij can be evaluated as in [4].We define the tangent flexural rigidity of cross-section as a ratio between incremental bending moment and incremental curvature while keeping constant the axial force (∆N=0) as: (8) Solving the system given by Eq. ( 7) and taking into account Eq. ( 8), the tangent flexural rigidity of the cross section with partial composite action can be developed as: where ∆N cf represents the incremental axial force in the concrete slab under the assumption of full composite action computed as a difference between the axial force in the concrete slab associated at given value of the bending moment M and the axial force associated at an incremented bending moment M+∆M [4].The value of effective degree of composite action defined in [4] is assumed to be constant over the length of the member and is updated at each loading step according with the existing state of stress [4].Need to be mentioned that for shear forces (P) less than 50% of P sc a constant value for the shear connection stiffness (k 50% ) is considered (Fig. 2), assuming a secant connector stiffness corresponding to 50% of P sc while for shear forces greater than50% of P sc , a secant value for the shear connection stiffness (k sec ) is considered (Fig. 2) [4].As already mentioned above the approach discussed in this paper assume for the entire member length a constant value for the function used to introduce the partial composite action f(γ eff )=constant.Hence the effects of non-uniform distribution of the shear connectors can not be taken into account and also the accuracy in detecting stress distribution along the member length can be only determined in a approximately manner.A more efficient approach able to overcome the above mentioned drawbacks is currently under development.Such an approach implies an explicit solution of the second-order differential equilibrium equation of the composite beam with partial composite action in which the axial force in the concrete slab represents the main unknown.This axial force, the solution of the differential equilibrium equation, can be expressed in function of the axial force under the assumption of the full composite action multiplied with a function of the degree of composite action which include also the exact distribution of the bending moment along the member length.Moreover, by simply dividing the beam according with the variable distribution of shear connectors and solving the second-order differential equilibrium equations for each segment considered as beams with uniform distribution of shear connectors, and then imposing the boundary and continuity conditions could represents a direct and simple way to extend the proposed approach to the cases of non-uniform distribution of shear connection along the beam length.In this way both uniform and non-uniform distribution of the shear connectors can be efficiently considered in the proposed formulation but this issues requires further investigations and calibrations and will be treated in a future work.

Second-order flexibility-based element
Flexibility-based method is used to formulate the distributed plasticity model of a 3D frame element (12 DOF) (Fig. 3).The spread of inelastic zones within an element is captured considering the variable section flexural EI y and EI z and axial EA rigidity along the member length, depending on the bending moments and axial force level, cross-sectional shape and nonlinear constitutive relationships as already described.Assuming elastic behavior within a load increment, and no coupling of axial and flexural responses at the section level, the generalized displacement in point i of the member produced by the force ∆P (∆M iy(z) ,∆M jy(z) ,∆T iy(z) ,∆T jy(z) ) could be expressed as in Eq. ( 12) where with superior indices (I) and (II) represent the first order and second order efforts respectively.The second term in Eq. ( 12) introduces the additional effect of shear deformations.In this respect, for composite beams with partial composite action, equivalent transverse shear stiffness has been derived by using the energy relations [4].For column cross-sections equivalent transverse shear stiffness is assumed to be computed taken into account only the contribution of the steel component.The second order bending moments and shear forces can be evaluated in function of the nodal bending moments and uniform distributed loads as described in [3] and then the relationship between nodal displacements (∆u r ) and the nodal efforts (∆s r ) could be further expressed by defining the element flexibility matrix (f r ) and nodal displacements given by the uniform distributed loads (δ r ) as: The detailed expressions for f r, δ r can be found in [3].To produce the deformational-stiffness relation, the Eq. ( 13) is inverted, obtaining the following deformational-stiffness equation: where the vector ∆q r is the incremental equivalent load vector, whereas k r represents the instantaneous element stiffness matrix of the beam-column element without rigid body modes, determined by matrix inversion of the flexural matrix f r .If the state of forces at any cross-section along the beam column element equals or exceeds the plastic section capacity, the flexural stiffness at the respective location approaches zero.Once the member forces get to the full plastic surface they are assumed to move on the plastic surface.Considering the member in Fig. 3 and that at the end of the member at node i the forces (N, M y , M z ) get to the full plastic surface the incremental bending moments and the incremental axial force N  can be linearly related as (Fig. 4):  For given value of axial force N the iterative procedure described in [2] is applied in order to determine the ultimate bending moments (M yp * , M zp * ) for a given bending moment ratio (tan = M z / M y ) and then Similarly, for given values of the bending moments M y , M z the axial resistance N * associated with a failure criterion is evaluated.Consequently the basic nodal element forces for the beam column element can be expressed in matrix form as in Eq. ( 15) or symbolically in condensed matrix form: where the transformation matrix T c introduces the correlation between nodal forces such that the plastic strength surface requirement at section "i" is not violated by the change of member forces after the full plastic strength of cross-section is reached (Fig. 4).Denoting r s  , u  and r q  as finite changes in the force vector, displacement vector and fixed-end force vector, respectively, the incremental forcedisplacement relationship for the element including the equivalent nodal loads can be expressed as: (17) or symbolically in condensed matrix form: represents the stiffness matrix and equivalent nodal loads vector for the element when a full plastified cross section forms at the i-th end of the element: Following a similar approach we can obtain the elasto-plastic stiffness matrix and equivalent nodal loads for the cases when a full plastified sections forms at j-th end of the member or at both ends.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 momentrelative rotation models for composite connections, and the floor slab action.For instance the mathematical model described in [3] can be useful to include the effects of the semi-rigidity and the penalty element method can be used to include the effect of the rigid floor diaphragm action.Alternatively, the floor slab may be modelled as flat shell elements (with 6 DOF's per node) and coupled in this way with the beam-column element developed in this paper.The member is assumed to consist of three zones as shown in Fig. 5.The first zone is the rigid block zone (finite-size joint) located at each end of the member.The second zone is the semi-rigid connection assumed to have zero length, and the remaining part of the member represents the beam-column element.The incremental change in the displacements and nodal forces for the semi-rigid beam-column element, considered without finite joint size (rigid block), may be symbolically expressed as: in which k r and r q  represents the stiffness matrix and equivalent nodal loads vector for the semi-rigid beam-column element, including the effects of material and geometric nonlinearity sources as described in the earlier sections, and in which the rows and columns corresponding to axial and torsional deformations have been removed, r u  represents the end element rotation vector (Fig. 5): and r s  represents the incremental nodal bending moments about the cross-sectional axes y and z acting at ends i and j of the element: Let us assume that the element nodal rotations are collected in the vector u  whereas the incremental nodal bending moments about the cross-sectional axes y and z acting at ends 1 and 2 of the element are collected in the vector s  : From geometry and assuming small rotations during each load increment the relationships between rotations at the ends (i and j) of the semi-rigid beam-column element, and the same quantities at the element nodes (1 and 2) can be explicitly expressed in matrix form as:  gives: Hence, the basic equation which define the incremental change in the rotations and nodal bending moments at the nodes 1 and 2 of the element can be expressed as: in which k and eq q  represents the stiffness matrix and equivalent nodal loads vector for the semi-rigid beam-column element including the effects of finite-sizes of the nodes.The resulting stiffness matrix is a 4x4 matrix, and does not include torsional and axial degrees of freedom.Torsional and axial stiffness components are then added to result in the required 6x6 stiffness matrix.To include rigid body modes, the stiffness matrix is pre-and post multiplied by a transformation matrix to result in the required 12 x12 matrix [3].It should be noted that panel zones are not necessarily rigid components as assumed in the proposed approach and high shear forces and deformations in panel zones may have a pronounced effect on the stiffness and ultimate strength of frame structures.

Computational examples
The accuracy of the analytic procedure and the computer program (NEFCAD) developed here has been evaluated using several benchmark problems.In the present approach, one element has been used to model each column and beam in all computational examples and the advanced numerical simulation is conducted by using the specialized software for nonlinear analysis of structures, ABAQUS as described in [4].

Simply supported composite beam with partial shear connection
The proposed numerical model is validated by comparisons against Chapman and Balakrishnan [5] experimental tests on E1 simply supported composite beam, as well as against the results predicted by ABAQUS software and those published by other authors.Details about geometry and material properties can be found in [4,5].The mid-span deflection vs. applied load comparative curves is plotted in Fig. 6a.As it can be seen the behavior of the composite beam predicted by the proposed analysis procedure (k 50% =2028N/mm 2 , γ=11.34) and the advanced finite element model developed in Abaqus is in close agreement with that of experimental test.The effectiveness of the proposed procedure is further assessed by varying the degree of shear connection (by means of varying the number of shear connectors-partial shear connection) and comparing the predicted curves with those obtained with more complex finite element analysis, as shown in Fig. 6b.By comparing the curves depicted in Fig. 6b, it can be observed, that the proposed method (Nefcad) predicts fairly well the nonlinear behavior and ultimate load capacity of the system when different levels (ranging from 47% to 136%) of shear connection are considered.

Six-story composite frame
The geometry loading conditions of Vogel's six-story two-bay frame are reported in [4].The yield strength of all steel members is 235 MPa while Young and shear Modulus are E=20500 MPa and G=7885 MPa.The compressive yield strength of concrete for columns is f c =26MPa and for concrete slab f c =16MPa.In order to evaluate the effects of the finite size of the joints and to make possible comparisons with more advanced nonlinear FEM solutions (Abaqus), in the proposed approach (Nefcad) the frame has been modelled considering member end offsets and assuming effective rigid joint size of one-half the true size of the joint.In this study a value of P sc =64.34 kN was considered for the shear connector capacity.The parameters that describe the shape of the studs constitutive law have been selected as proposed by Ollgaard: β=1 and α=0.558, obtaining in this way a value of 94375 N/mm for the connector stiffness corresponding to 50% of P sc , namely K 50% (see Fig. 2).The values of shear connection stiffness, k 50% , are evaluated based on the number of shear connectors corresponding to the desired degree of shear connection [4].Fig. 7a presents the comparative load-deflection curves obtained by the proposed approach when full composite action is assumed and those obtained with the advanced FEM model considering different levels of shear connection, ranging from 100% to 300%.It can be observed that, as expected, increasing the level of shear connection the system became more rigid with increased strength and stiffness and the inelastic behavior becomes similar with the response predicted by the proposed approach (full composite action).The effectiveness of the proposed procedure is further assessed by varying the degree of shear connection and comparing the predicted curves with those obtained with more complex finite element analysis, as shown in Fig. 7b.

Conclusions
A reliable and robust nonlinear inelastic analysis method for composite steel-concrete frames with partial composite action has been developed.The proposed formulation has been found to be effective in predicting the global behavior of composite beams and frames with partial shear connection, both in elastic and post-elastic field and, the numerical results agrees fairly well with the experimental results and those obtained by advanced nonlinear FEM approaches but with much less computational effort.The proposed formulation takes advantage of using only one 2-noded beamcolumn element and features, in this way, the ability to be used for practical applications by combining modeling benefits, computational efficiency and reasonable accuracy and may circumvent some of the difficulties that may arise enforcing the compatibility conditions at the semi-rigid composite connections.

Fig. 3 .
Fig. 3. Beam-column element with 12 DOF (a) beam element; (b) column element.The basic incremental force-displacement relationships are determined considering the element represented in natural coordinate system with rigid body modes removed[3].The total bending moment associated to given value of axial force N and N * represents the ultimate axial force for given value of total bending moment  and N * ) are determined as follows.Let us consider that at the current loading step in a specified location the forces exceed the plastic surface.

Fig. 4 .
Fig. 4. Interaction diagram for given bending moment's ratio.The plastic surface requirements.
2018, Universitat Politècnica de València Chiorean, C.G. and Buru, S.M. 2018, Universitat Politècnica de València which l represents updated length of the element computed during each load increment, and a, b represents the sizes of the rigid blocks located at each end of the member.Symbolically matrix L introduces the effect of finite joint sizes.Substituting Eq. (30) for r u  into Eq.(24)