Second-order two-cycle analysis of frames based on interpolation functions from the solution of the beam-column differential equation

Abstract In geometrically nonlinear problems solved using the Finite Element Method (FEM), the structure response is directly influenced by the level of discretization and the nonlinear solution algorithm used. To reduce the discretization dependence, exact solutions are developed based on the deformed infinitesimal element equilibrium. To deal with the nonlinear solution problem, the two-cycle method can be used, since it is not dependent on load or displacement steps. The two-cycle method developed by Chen & Lui (1991) uses the classical geometric matrix and is not accurate for high axial loads. This happens because the geometric matrix is obtained using Hermitian polynomials which are approximate solutions. To circumvent this issue, the frame element’s tangent matrix is obtained using interpolation functions that match the homogeneous solution of the differential equation of the beam-column problem. The main objective of this study is to carry out a second order analysis of the frames and obtain equilibrium paths using the two-cycle method and the tangent stiffness matrix based on solutions of the differential equations obtained from the element’s deformed configuration. The results in terms of displacements and rotations for the examples studied are identical to the analytical solutions, showing that the combination of the two-cycle method with the exact element formulation is promising and can diminish the need for discretization.

The evolution of engineering demands the development of economic structures with less weight and material consumption without any trade-off in durability and safety.Hence, the use of high strength materials and slender elements becomes more common in new ventures.Under this view, the consideration of geometric nonlinearity is indispensable, especially for slender elements.
In context of the geometric nonlinearity of frames, a second-order stability analysis can be carried out in many ways, from simple checks based on first-order linear analyses to higher-order nonlinear numerical analyses.Despite the existence of several analytical solutions to solve stability problems, numerical solutions using the Finite Element Method (FEM) are often employed.From this numerical approach, the continuous structure response is directly influenced by the engineer's experience in choosing the most suitable number of elements to be used in structural analysis (discretization).All these influences occur because the FEM discrete solution approximates the analytical solution, i.e., the interpolation functions that define the structure deformed configuration do not always agree with the problem's exact solution (Burgos & Martha, 2013;Rodrigues et al., 2019).
In literature, to reduce element discretization dependence using the FEM in nonlinear geometric analysis, several studies have been carried out, such as Burgos et al., (2005), who used additional degrees of freedom within the elements in a classi-cal linearization of the stability problem.Also, Chen & Lui (1991) and Aristizabal-Ochoa (1997, 2008, 2012), used classical or modified stability functions.Other authors have sought solutions based on infinitesimal element equilibrium in their deformed or undeformed configurations considering a combination of effects, such as transversal and axial loads, shear deformation (Timoshenko beam theory) and elastic foundation.Davis et al., (1972) and Nukulchai et al., (1981) formulated the exact solution for beam elements considering shear deformation.Zhaohua & Cook (1983) and Ting & Mockry (1984) developed the finite beam element considering two and one-parameter elastic foundations respectively, but with no shear deformation.Shrima & Giger (1992), in turn, formulated Timoshenko's beam finite elements considering a twoparameter elastic foundation.Onu (2008) was the first author to formulate a finite element combining all effects (transversal and axial loads, shear deformation and two-parameter elastic foundation).Burgos & Martha (2013) presented exact shape functions including transversal and axial loading and shear deformation.Rodrigues et al. (2021) developed the three-dimensional tangent stiffness matrix and nodal equivalent loads, using exact interpolation functions, considering axial load, shear deformation and high-order terms of the Green-Lagrange strain tensor.
In addition to the exact finite elements previously addressed in literature, a nonlinear analysis methodology is also needed in numerical methods to describe the equilibrium path of a structural system.Thus, incremental, and iterative methods are commonly used, and some stand out among them: Newton-Raphson method, secant method and arc length method (Wempner, 1971;Ricks, 1972Ricks, , 1979;;Ramm, 1981;Crisfield, 1983Crisfield, , 1986)), displacement control method (Yang & Shieh, 1990;Yang & Kuo, 1994) and indirect displacement control method (de Borst, 1986(de Borst, , 1987)), among others.Along with these traditional general methods, there are simplified approaches capable of solving nonlinear problems with less computational effort, such as the γ z coefficient method (NBR 6118, 2014), moment amplification method (Chen & Lui, 1991), P-Delta method (Chen & Lui, 1991) and, finally, the two-cycle method (Chen & Lui, 1991), which is applied in this study.The main objective of this research is to carry out a second order analysis of frames and obtain equilibrium paths using the two-cycle method and the tangent stiffness matrix based on exact solutions of the differential equations obtained for the element's deformed configuration.
The article is organized in the following manner: the next section introduces the differential equations that govern the problem.Section 3 presents the exact tangent stiffness matrix.The twocycle method is presented in section 4 and some examples of second order analysis of frames are presented in section 5. Finally, in section 6, the authors conclude this article and suggest ideas for future research.

Introduction
The shape functions are used to interpolate the nodal values of the discrete element, obtaining displacements and rotations within the element of the continu-ous solution of the problem, according to Equation (9): Inverting matrix [H] in Equation ( 8) and substituting in Equation ( 7), shape functions can be easily written: According to Burgos & Martha (2013), P can be either a tensile axial load (P>0), in which case μ is a real number and the expression for v h (x) can be written in terms of hyperbolic functions.On the other hand, in the case of a compressive axial load (P<0), μ is a complex number and the expression for v h (x) can be written in terms of trigonometric functions.All those expressions are available in Rodrigues et al. (2019).In matrix form, the conditions given by Equations ( 5) and ( 6) are written as: The same process can be performed using the expressions in terms of hyperbolic or trigonometric functions and the stiffness matrices are then obtained accordingly.
With boundary conditions given by:

Tangent stiffness matrix based on interpolation functions from the beam-column differential equation
The differential equation solution results in lateral displacement v h (x) and the rotation can be obtained using θ Herein, the Navier (Euler-Bernoulli) theory is employed.Thus, the rotation and transversal load can be written as the derivative of the lateral displacement and shear force respectively, leading, after some algebraic manipulation, to Equation (3): The previous equation can be rewritten considering only its homogeneous part (q(x)=0), which leads to v h (x): (3) (2) In the case of μ being a complex number (P<0), coefficients in terms of trigonometric functions are: The axial components of the tangent stiffness matrix of the frame element are the same as those of the classical elastic stiffness matrix (EA/L), easily found in literature, in which E is the material Young's modu-lus, L is the element length and A is the cross-section area.
In the case of μ being a real number (P>0), coefficients in terms of hyperbolic functions are: (µL) 3 (e µL + 1) (µL) 2 (e µL -1) µL [µL(e 2µL + 1) -(e 2µL -1)] µL [(e 2µL + 1) -2µLe µL ] P EI D = µL (e µL + 1) -2(e µL -1), µ = Burgos & Martha (2013) presented the expressions and plots of these shape functions for all cases described above, as well as the development to obtain the tangent stiffness matrices described in Equations ( 11) to ( 13), based on the Principle of Minimum Potential Energy.In case the homogeneous solution is given by exponential functions, the coefficients of the tangent stiffness matrix are given by: Chen & Lui (1991), the two-cycle method uses only two iterations to obtain second-order ef-fects.In this method, the classical Equation ( 14) is utilized in two steps.In the first step, an elastic linear analy-sis is carried out using the classical elastic stiffness matrix -easily found in literature.
In this study, MATLAB ® ( 2018) was used to obtain the solutions for the models shown in Figure 2. The equilibrium paths for second-order analysis using the two-cycle method with exact shape functions were compared to the analytical solutions and results obtained using classical frame elements (Hermitian shape functions) through MASTAN2 software (McGuire et al., 2000), which adopts the predictorcorrector method for the nonlinear analysis.For the two-cycle solution, each load value results in a displacement vector which is independent from previous steps, since the method is used to obtain the exact displacement solution, and not in an iterative manner.
Figure 2a shows a laterally restrained (non-sway) frame subjected to a vertical load P and bending moments M = αPL = 0.001PL.Properties of all bars are: L = 6 m, E = 10 8 kN/m 2 , I = 10 -5 m 4 , and EA = 10 6 kN.The critical load is P cr = 2.5515 π 2 EI ⁄ L 2 = 699.51kN and the analytical solution in terms of the free rotation is given by: A look at Figure 3, shows that the results for the two-cycle method using the exact tangent stiffness matrix is equivalent to the equilibrium path of the analytical solution and better than the MASTAN2 solution discretized in 1, 2 and 4 elements (M1, M2 and M4).In this example, it is clear how the engineer's experience would be required in discretization, especially considering that the difference between the results produced using discretization M1 and M2 is 44.10%.One of the main objectives of this study is to reduce the need for discretization in simple problems, such as nonlinearities in planar frames.From the axial load components obtained in the previous step, the second step starts with the substitution of these loads in the tangent stiffness matrix coefficients.A second solution is obtained using Equation ( 14) again, but replacing the elastic stiffness matrix with the tangent stiffness matrix.Thus, in this second step, displacements and forces will be of second order, i.e., considering geometric nonlinearity.Herein, the tangent stiffness matrix coefficients are given by Equations ( 11) to (13).In the original two-cycle method, these coef-ficients are obtained by the sum of the elastic and geometric matrices and are not as accurate.One of the problems that arise from this approach is the possibility of numerical singularities (division by zero) that come from low values in axial loads.This can be easily solved by adding a tolerance that flags the need for a nonlinear analysis.
It is interesting to point out that in obtaining the differential equation, it is assumed that the axial load is known.Thus, in this method, the first step is necessary only to calculate and increment the axial loads for the second step.In short: the major contribution to achieve reliable results closer to the analytical solution is in the exact approach of the tangent stiffness matrix in the second cycle.Since the differential equation is obtained assuming that the axial load is constant within the element, the discretization using one element per bar is not only advantageous, but also mandatory.Discretization using more than one element per bar is not comparable with the analytical solution, for the same reasons exposed before.Figure 4 shows that the results for two-cycle method using the exact tangent stiffness matrix are identical to the analytical solution.In other words, the combination of the two-cycle method with exact element formulation provides excellent results in predicting the nonlinear geometric behavior of the sway frame.This article presented the combination of a simple method for secondorder analysis with a formulation of the element tangent stiffness that uses shape functions that match the homogeneous solution of the differential equation in deformed configuration to solve nonlinear geometric problems.
From the examples shown, the analyses using the two-cycle method with the suggested formulation of the tangent stiffness proved to be efficient when compared to traditional numerical solutions; furthermore, results are also independent of model discretization.It is interesting to notice that result accuracy in numerical solutions by finite elements strongly depends on the boundary conditions of the models.The sway frame and cantilever beam achieved better results than the nonsway frame in one element discretization (M1) for the numerical solution in MASTAN2.This issue alone shows the need for a tool that does not rely on user experience in discretizing the model and providing nonlinear parameters (load step, solution algorithm).This problem is no longer present when the combination of the two-cycle method with exact tangent stiffness is applied.It is important to emphasize that all the analyses are performed for small displacements and loads that are below critical values, i.e., postcritical behavior is not considered.
In future studies, the authors suggest the development of analogous analyses that consider shear deformation, elastic foundation and/or variable cross-sections.Another suggestion is the use of distributed loads.
Figure2bshows an unbraced (sway) frame subjected to a vertical load P and horizontal loads