Direct plastic analysis of steel structures by flexibility-based element with initial imperfection

Second-order direct analysis has been used in some regions for reliable analysis and design of steel structures. Currently, the stiffness-based element is widely used with accuracy improved by enforcing equilibrium along mid-span or “stations” along the member length in order to achieve equilibrium which is not guaranteed along an element. In this paper, a flexibility-based beam-column element considering member imperfection based on Hellinger-Reissner functional is developed and used for practical second-order direct analysis. This new element is a flexibility-based element with member initial bowing at the element level for direct analysis of three-dimensional frame analysis whereas previous flexibility-based elements assumed perfectly straight geometry for the element. The fiber plastic hinge approach is adopted to account for the distributed plasticity of a section. The new flexibility-based element performs excellently for modeling of members under high stress with material yielded as the conventional stiffness-based element has less accuracy when few elements are used in modeling a plastic member. This will significantly enhance accuracy and computational efficiency for direct plastic analysis which can then be more widely used in practical design. Several examples are employed to validate the accuracy and efficiency of the proposed element along this line of thought.


Introduction
Steel structure is mainly made of steel, and is one of the main kinds of building structure.It has many advantages, such as, high strength and good ductility.Stability problems are important during the design of steel structure.However, the traditional linear analysis cannot satisfy the requirements of current steel structure design.Nowadays, design codes of steel structure in many countries have incorporated the theory of direct analysis, and recommend this new method which can replace the conventional linear analysis and the design method based on effective length.The design method of calculating effective length is rarely used or eliminated in Eurocode-3 2005 [1], LRFD 2010 [2] and HKSC 2011 [3], and is displaced by second-order direct analysis.
The traditional design method needs to classify types of steel structure.These structure cannot be designed by linear analysis when their critical elastic coefficient is less than a particular value, for example, 3 in Eurocode-3 2005 [1], 4 in HKSC 2011 and AS4100(1995).However, elastic critical load has limited applications, such as, regular buildings mainly under gravity load.It is not able to measure many complex structures, for example, large power transmission tower, scaffolding and spatial latticed structures, which lead to the invalidation of linear analysis and design.Second-order direct analysis has been a preferred method in current design of steel structure.LRFD 2010 [2] includes the secondorder direct analysis method in the core chapter, while the linear analysis method in appendix.Eurocode-3 2005 [1] puts it in front of the linear analysis method.Theoretically, the Liu, Y.P., Shu, G.P. and Chan, S.L. 2018, Universitat Politècnica de València variation of structure and members under load action is not considered and internal force results are inaccurate by the linear analysis.Its assumption is that all stiffness comes from material and geometrical features of structures, neglecting the effect of load applied on members which makes no difference of member's capacity under tension and compression, while the long practical experience shows that compression members' bearing capacity significantly lower than tensile ones.Therefore, LFRD 2010 [2] requires using reduction factor to cover the shortage of linear analysis.More importantly, the value of effective length is uncertain, maybe overestimated or too low, and the effect of eccentricity is hard to be measured.That is why current design codes try to avoid using effective length method to structures' design.
Chan and Gu [4] derived a stability equation of bending according to the curvature shown in the table 5.1 of Eurocode-3 2005.Chan [5] applied this equation in second-order direct analysis and design of semi-rigid steel frame.Fig. 1 shows load-deflection curves obtained by different analysis and design methods.These methods' ultimate bearing capacities are also compared with u got in experiment.The failure criterion of structure is appearance of the fist plastic hinge when designed by secondorder elastic analysis.The first plastic hinge method can capture failure of the first member in structure, which is loading coefficient y when the first member begin to yield.This compromise method is relatively conservative because traditional plastic hinge method can only consider plastic hinge happened at member's two ends and disturbed plastic zone method needs huge computation.
In the present paper, a beam-column element based on flexibility method, for using secondorder direct analysis of structure, has been introduced, where disturbed plastic zone has been considered.This element has the capacity of directly modelling geometrical initial imperfection of members, and has a high accuracy so that it is able to model one member by one element.These features can significantly reduce engineers' workload of modelling and calculation time of analysis.Especially, this method eliminates inconvenience of modelling one member's geometrical initial imperfection by multi-elements.At the same time, considering nonlinearity of material by disturbed plastic hinge method, it can reflect development of yielding zone so as to capture the response of members relative to traditional plastic hinge element.Compared with beamcolumn elements with disturbed zone method, it has advantages on computational speed because there is no section integration during every cycle.

Description of beam-column element based on flexibility method
The shape function of beam-column element based on flexibility method is the equation of force and error of this kind of element only comes from the integration along the member.For this reason, they have higher accuracy than these based on stiffness method.The proposed element incorporates effect of P-and material nonlinearity, as well as geometrical initial imperfection in the elemental stiffness matrix for practical use.Derivation of this element is as follows.

Hellinger-Reissner vibrational method
The displacement-based elements are generally derived by the principle of minimum potential energy while the flexibility-based elements are commonly based on the Hellinger-Reissner (HR) variational principle which is expressed in equation ( 1

2018, Universitat Politècnica de València
Taking the first variation for equation (1) with regard to the displacement and the stress resultant and setting it to zero, the stationary of the Hellinger-Reissner potential is expressed as Thus, the weak form of equilibrium and compatibility equations can be obtained in equations ( 3) and ( 4) respectively.

Equilibrium equation
The equilibrium equations defined in equation ( 3) can be expanded as with  0 and  0 are geometrical initial imperfection.The section forces () in related with end forces  can be presented as

Compatibility equation
The compatibility equations given in equation ( 4) can be further expressed as From the virtual work principle, there exists a relation between the increment of virtual internal forces and virtual end forces given by ∫ ()  () =     (8) The relation between the end displacements  and the section deformations corresponding to the generalized strains  along the member can be obtained as

Elemental flexibility matrix
Taking derivative of the end nodal displacements  in equation ( 9) with respect to end nodal forces  , the element flexibility matrix can be established as From above, once the displacements () and () as well as initial imperfections  0 () and  0 () are known the element flexibility matrix can be determined.

Transformation from basic to global system
To incorporate this new flexibility-based element into the conventional stiffness-based Liu, Y.P., Shu, G.P. and Chan, S.L. 2018, Universitat Politècnica de València software package like NIDA [6], the element stiffness matrix   in basic coordinate system can be obtained as in which   is the flexibility matrix defined in equation (10).
In this paper, the co-rotational method used in [7] is adopted to carry out the transformation between the basic and global systems.Under the co-rotational framework, the tangent stiffness matrix   of beam-column element can be calculated as in which  is transformation matrix from basic to local system,  is a matrix considering the work due to the initial force and the translational displacements,  is transformation matrix from local to global system.The details of the matrixes ,  and  can be found in [8].

Yield surface of frame section
The concentrated plasticity method assumes that plastic hinge only happens at the end of members, and imports internal degree of freedom to describe material nonlinearity approximately.This method The procedures of this method can be described as follows.
Step 1: Defining plastic hinge model before analysis.
Step 2: Determination of elemental state such as elastic or plastic.
Step 3: Insert a plastic hinge in the element.
Step 4: Condensation of internal degrees of freedom, and determination of whole element state.The predefined plastic hinge model can save computational time and should be more efficient.However, the location of plastic hinge is limited at the end of a member and as a result it is hard to capture actual material nonlinearity behaviour along a member.Some researchers have proposed arbitrarily-located-plastic-hinge beam-column element.This kind of element condenses internal degrees of freedom essentially, which has a complex deduced process and a poor expression.
In order to get   () in equation ( 10), traditional disturbed plastic zone method discretises frame section into many fibers.Generally, the precision of section internal force and stiffness depends on the number of fibers.Each fiber needs individual parameters to record material status.This method has two disadvantages: (1) the computer time increases as the number of fibers increases, especially in nonlinear dynamic analysis; (2) generous fibers also increase the requirement of memory.Because of these, a section constitutive model is proposed in this paper to simplify the procedures of second-order inelastic analysis with acceptable accuracy.

Elemental flexibility matrix
Following the concept of constitutive relation of metal, Krenk [9] and Powell [10] proposed constitutive relation at the section level.Lu [11] deduced refined plastic hinge formulation, which is able to capture plasticity development.Although the yield criterion by P-My-Mz function is used to detect the sectional state, the sectional forces such as P, My and Mz do not consider the coupling effect between the axis force P and bending moments My and Mz.For the steel member with compact wide-flange section, Orbison [12] proposed a yield surface reproduced in equation ( 17) to trace the material nonlinearity.Fig. 2 shows the typical full yield surface proposed by Orbison [12].
in which z-axis and y-axis are the major and minor axes respectively;   ,   and   are the axial force capacity of the cross-section, full plastic moments about the major and minor principal axes respectively.
Van Long and Hung [13] proposed a strain hardening rule for frame section, which is applied in a beam-column element with traditional plastic hinge.Their strain hardening rule is able to describe three ranges, i.e. elastic range, strain hardening range and flowed range.The rule can be expressed as where  is the parameter for hardening and given in equation ( 21 in which ̅  is the effective strain; ℎ and  are, respectively, the depth and the wide of the section;  is the member length; ̅   is the limit effective strain;  is the strain hardening modulus.The differential of equation ( 19) is given below

Derivation of section stiffness matrix
The incremental sectional force can be calculated as in which [  ] is the stiffness matrix of section.When the section is in elastic state, the elastic stiffness matrix can be written (torsion-stiffness is not included for simplicity) For the elastic case, the matrix [] is same as that generated by fiber section method.It should be pointed out that the fiber section method requires much computer resources since many members remain in elastic range in practical projects at low load level.
The sectional deformation composes of elastic part and plastic part.
The incremental section internal force is only determined by the incremental elastic deformation.Thus, one has When the section is in plastic range, the force-deformation relation is expressed as in which [  ] is the plastic matrix and calculated by Finally, the section stiffness matrix can be obtained as Liu, Y.P., Shu, G.P. and Chan, S.L.

Determination of section elastic-plastic status
Incremental-iterative method is the general approach to solve nonlinear problems.The elasto-plastic behaviour is reflected by nonlinear constitutive law on the basis of stress resultant.The nodal displacements by solving system equations and residual forces will be extracted for each element and further converted to basic system to calculate the section deformation.The detailed process of incremental-iterative method to determine sectional state is introduced as follows.
Step 1: calculation of trial internal force at section level.
Assuming that the monitoring section is elastic, the initial elastic stiffness matrix is used to calculate trial incremental sectional forces as d[S()] = [  ]d[d()](37)Then, the total trial sectional forces of the section is Step 2: calculation of section yield function.
The value ϕ(F, K) of section yield function can be calculated using the trial sectional forces obtained in previous step.If section state is in elastic range in last step, go to step 3, otherwise jump to step 4.
Step 3: section state is elastic in last step.
If ϕ(F, K) ≤ 0 , the section is still in elastic range, the trial sectional forces are the final forces, and exits the loop.
If ϕ(F, K) > 0 , the coefficient λ is obtained by iterative solving, and then it goes to step 5.
Step 4: section state is in plastic range in last step.
If ϕ(F, K) > 0, the coefficient λ is equal to zero, and then it goes to step 5.
If ϕ(F, K) ≤ 0, the section is in unloading phase and assumed to be elastic.The trial sectional forces are determined; exits the loop.
Step 5: calculation of section internal forces and tangent stiffness.

Verification examples 4.1 A single column under compression
Fig. 3 shows two columns with different boundary conditions.The section is W12x4, with A = 4.13x10 -3 m 2 , I = 1.94x10 -6 m 4 and r = 0.02166 m.The yield stress and elastic modulus are taken as 275 kN/m 2 and 2.05x10 8 kN/m 2 .The hardening ratio is 0.01.The column is subjected to a concentrated point load at the top point.The member is modelled by one proposed element, and the section is represented with the proposed model.From Fig. 4, it can be seen that the initial imperfection will significantly affect the member behavior.

Six-story frame
The six-story frame as shown in Fig. 5 subjected to distributed gravity loads and concentrated lateral loads is studied here.

Conclusions
In this paper, a beam-column element, with disturbed plastic hinges, based on flexibility method is proposed to account for geometrical and material non-linearity behaviour of frame structures.This element incorporates geometric initial imperfection with high accuracy and therefore it can use one element to model one member.This element can significantly reduce engineers' modelling workload as well as computational time of nonlinear analysis.The section constitutive model derived from section yield function can replace the section integration by fiber section approach to represent the relationship between the internal forces and deformations.This technique has advantages on computational efficiency and simulation on plasticity development and is ready for static and dynamic nonlinear analysis.

Fig. 3 .
Fig. 3. Layout of columns.The rotation of the top point is shown in Fig. 4 (a) and (b) against the applied force.

Fig. 5 .
Fig. 5. Layout of six-story frame.The displacement of top right joint in horizontal direction is shown in Fig.6.The ultimate load factor is around 1.1 which is close to previous research.The equivalent geometrical imperfection of "L/500" and the separate consideration of geometrical imperfection "L/1000" with residual stress are considered in this example.It can be seen that the two methods produce well-agreed results.