Modal P-delta – simplified geometric nonlinear method by using modal and buckling analysis

Abstract Recent works show that in structures, there is a relationship between the natural period of vibration and global second-order effects. This relationship occurs because both depend on the stiffness matrix and the mass matrix of the structure. In this work, a non-linear geometric method - modal P-delta – is proposed that takes advantage of the relationship between the dynamic parameters and the non-linear effects. The methodology establishes an association between the buckling instability modes and the structural vibration modes. An interpolation of the vibration modes without axial loading and vibration modes with critical axial loading is proposed to approximate the vibration modes and frequencies of the loaded structure. In this way, through a simple formulation, the vibration modes and the natural frequencies of the loaded structure can be used to evaluate the displacements of the structures including the non-linear effects. Several numerical examples were simulated in regular structures in the plane, such as a free-fixed column and a plane frame with two different loading configurations. The results generated with the modal P-delta method provide information about the nonlinear behavior of the pre-buckling equilibrium path. These results are different from the findings in the literature, where the relationship between dynamic parameters and non-linear effects is used as a simple indicator or amplification factors to determine non-linear effects. Furthermore, our results indicate that the modal P-delta reduces the computational time spent considerably compared to the traditional P-delta iterative method.


INTRODUCTION
The movement of the structure to a deformed position in the analysis of structures subjected to vertical and lateral loads causes additional internal forces called second-order effects. According to Rutenberg [1], this type of additional effect is not considered in normal static analysis and can be important for structural analysis of tall and slender structures. The second-order effects are also called P-delta effects since additional internal moments of the structure can be evaluated by multiplying of structural weight "P" on each floor and the displacement "delta" caused by lateral loads.
Structures with lateral loads and high axial internal forces can require a non-linear geometric analysis -NLG for the correct evaluation of the internal forces in the structure. Traditionally, the NLG effects are evaluated in two ways: with iterative methods [2] or with direct methods [1]. Commonly, iterative methods have a high computational cost in the analysis of large structures and cannot be appropriate for dynamic analysis where there may be structural instability due to resonance vibrations and P-delta effects simultaneously [3]. These characteristics have made its use difficult in common engineering projects. Alternatively, direct methods try simplifying or approximate the geometric stiffness matrix to avoid iterations and reduce computational time with an acceptable degree of approximation [1]. Another alternative, that can be classified as direct methods, is the use of simplified amplification factors to the first-order values [4], [5]. Several works can be found in the literature exploring these two alternatives, detailed information of NL effects and methods can be found in Gaiotti and Smith [6]. Also, the work by Rezaiee-Pajand et al. [7], [8] presented a comparative study related to the formulation, characteristic, and efficiency of several well-known structural geometrical non-linear -NL -solutions in iterative methods.
Recently, new methodologies have emerged to consider the NLG effects in structures, based on changes in dynamic parameters by the axial loads in structures [9], [10]. These methodologies can be considered as direct methods and are mainly proposed with two aims 1) as a method to evaluate the amplification factors of first-order analysis to approximate second-order analysis and 2) as an indicator of the structure's susceptibility to second-order effects. Specifically, the work by Statler et al. [11] proposes using a single degree of freedom model of the structure to define second-order amplification factors from the first natural frequency of the system. In that work, it was considered that the stiffness of a single degree of freedom system can be modified by a moment amplification factor defined in a highly simplified linear deformed shape of a rigid structure. In this way, it is possible to evaluate the natural frequency of the structure as a function of the moment amplification factor. Reis et al. [12] uses a similar procedure applied in Statler [11], however, it uses a dynamic model of a flexible structure based on equivalent cantilever beam-column to estimate the first natural frequency. In another work by Leitão et al. [13], the proposed method by Reis et al. [12] was used to evaluate the applicability of the amplification factor based on natural frequencies in an asymmetric structure in which torsional vibration modes can become fundamental vibration modes.
The previous methodologies based on the natural frequency of the structure show facility to estimate the amplification factors due to second-order effects since that amplification factor depends only on the natural frequency and thus requires only a modal analysis. However, these methodologies are based on approximations of the amplification factors that are used to perform rapid verification of the susceptibility for second-order effects, commonly used in standards or codes for the design of reinforced concrete during preliminary stages of design. Therefore, this amplification factor is only valid for a particular type of structure (e.g. sway frame) and limited for small values of amplification, smaller than 1.3. Additionally, these methods cannot determine the equilibrium path curve before buckling load.
To overcome the various limitations of based methodologies on dynamic parameters to estimate second-order effects and the problems related to the computational cost in iterative methods, this paper proposed a Modal P-delta method to estimate second-order effects. The proposed methodology provides the NL pre-buckling equilibrium path, and not only an indication of the susceptibility to second-order effects or definition of an amplification factor. In this way, the proposed method takes advantage of the relationship between dynamical parameters and NL effects to create an alternative methodology, which can be computationally cheaper than iterative traditional methods. The Modal Pdelta establishes structural displacements from natural frequencies and vibration modes of the loaded structure, including P-delta effects. These dynamic parameters of the structure can be approximated from an interpolation of vibration modes obtained with the unloaded structure and with the structures under the critical load. Different structural examples were addressed to demonstrate the capability of the proposed method to evaluate the pre-buckling equilibrium path. The results and computational time-consuming of the proposed method were compared with the iterative traditional P-delta method.

THEORETICAL BACKGROUND
This section contains a description of the theoretical bases to define the relationship between NLG, modal analysis, and the buckling phenomenon. First, the effects of axial load in the vibration modes and frequencies will be shown. With this, a relationship between the buckling mode with the vibration mode of the structure under the critical load can be established. Initially, these concepts in a simple example are presented -a column element -and later they are applied in a frame structure. In summary, it was shown that any vibration mode and natural frequency can be defined using the orthogonality properties of modes, and an interpolation of the vibration modes obtained using the unloaded structure and the vibration modes obtained using the structure under the critical load. Finally, with these concepts, the modal P-delta method will be proposed.

Effect of axial load on vibration modes and frequencies
Here will be shown as the vibration modes and natural frequencies are altered when the internal loads are increased in a Euler-Bernoulli beam element. A solution for this problem was presented by Shaker [14], however, some modifications were made to this work to obtain vibration modes and natural frequencies in a simple manner. The governing equation of the beam shown in Figure 1 is given as where �( , ) is the transverse deflection; ̈= ∂ 2 u( , )/ ∂t 2 is the acceleration; ( , ) ′′ represents the second derivative related to ; � = is the mass per length unit; is the Young modulus; is the moment of inertia of the cross-section; represents the axial load (in this case the internal axial load and the axial external load have the same value = ); and ( , ) is the transverse external load. In Figure 1b represents the shear load and the bending moment. The dynamical characteristics of free vibration system are obtained with ( , ) = 0. Applying in Equation 2.1, we obtained the following expression: 2) The solution of Equation 2.2 according to Shaker [14] and Paz [17], is found by substituting ( , ) = ( ) ( ), which correspond to a variable separation principle, then it can be obtained: where, Then, considering Equations 2.6 and 2.7, the solution ( ) is rewrite as, 10) or, 11) where , , and are the constant of integration, that are determined using the boundary conditions of the beam. Therefore, knowing Equation 2.11, it can be found a specific solution for the beam. The boundary conditions are those given for a column element as shown in Figure 2, The characteristic equation is obtained, which is the non-trivial solution, applying the boundary conditions, Equation 2.12, in Equation 2.11, The roots of Equation 2.13 allow the definition of the analytical natural frequencies of the column considering the effect of axial load on the element. It is now intended to show that the characteristic equation, Equation 2.13, can be written as a function of a single variable, 1 . Thus, operating in Equations 2.6 and 2.7 is obtained the following expressions.
The roots of Equation 2.18 represent the natural frequencies of a free-fixed column under axial load, which are expressed as Vibration modes are obtained considering the boundary conditions shown in Equations 2.12 and 2.11. Assuming the constant = 1, rewrites Equation 2.11 as, where constant is given by the expression, Thus, the expression (2.20) represents the vibration modes of the structure. Equations 2.19 and 2.20 are used to analyze the effect of axial stress on the natural frequencies and vibration mode of a generic beam. Thus, in Figure 3 the variation of the first natural frequency is shown, , as the load increases. The values shown in the figure, in the vertical and horizontal directions, are normalized, for 1 and by the critical load , respectively. It can be noticed as the axial load increases, the frequency of the structure decreases, and when the equilibrium path approaches the bifurcation point ( = ), the frequency becomes zero.

Relationship between buckling mode and vibration mode in a free-fixed column
In this section a column will be analyzed, obtaining the first vibration mode for the critical buckling load and the first buckling mode, which will be compared to each other, showing the relationship between them. Initially, the first buckling mode will be analytically represented. Considering the column in Figure 2, the differential equation of deflection is given as, The solution of Equation 2.22 is, .20 the first vibration mode will be found when the column is under the action of the first critical buckling load. When the axial load is equal to the critical load ( ), then 4 = 0, which makes the first natural frequency equal to zero. With this result and considering Equation 2.6 it can be concluded that 1 is also equal to zero.
Substituting in Equation 2.8 the first critical buckling load, and considering Equation 2.17, it is shown that, (2.28) With the previous result, it is possible to find the first vibration mode by rewriting Equation 2.20, thus, Comparing the first buckling mode, Equation 2.27, with the first vibration mode at critical load, Equation 2.29, it is noted that the two equations are linearly dependent. This allows us to conclude that, in the case shown here, the first vibration mode in the critical load and the buckling mode are equivalent.

Modal analysis and linear-buckling analysis of a structure
In this section, the concepts described in sections 2.1 and 2.2 will be expanded to any reticular space structure.

Finite element model
The balance of the force of an element in the deformed condition is given by the following relationship [15], where [ ] is elastic stiffness matrix, � � is the geometric matrix, { } is the nodal load vector of the element, { } is the nodal displacement vector { }. In detail, the elastic stiffness matrix of the element is expressed as, where [ ] represents a matrix containing the derivatives of the shape functions [ ] and [ ] describes the matrix material constitutive; � � the geometric matrix is given by, Where is the axial internal force, and finally, the vector { } is the nodal force vector of the element, it is written as, where { } is the nodal force vector applied directly to the nodes of the element � � the equivalent nodal force vector due to distributed forces acting along of element.

Linear buckling analysis
In a buckling problem, internal axial loads are sought that cause the sum of the stiffness matrix and the geometric matrix to be a singular matrix [16], thus the problem is given by, .34 represents a homogeneous system of equations, which can be formulated as an eigenvalue problem. Where are the eigenvalues and are the buckling modes. A load { } must be defined to estimate the geometric matrix that depends on the internal axial stresses. Following the definition of critical load, the magnitude of { } is not important since the critical load will be dimensioned by the load multipliers .

Modal analysis
The modal analysis aims to obtain the response of the structure through vibration modes and frequencies [17]. In this way, it will be shown how these dynamic parameters are obtained when the geometrical matrix is included. Equation of motion can be written considering the geometrical matrix, To study the natural vibration of the structure, the load vector applied to the structure must be null { ( )} = 0 and it is assumed that the nodal displacements vector can be expressed as, The buckling problem shows that the total stiffness matrix is singular when the first critical buckling load is reached, therefore the buckling mode obtained for the first critical load � � belongs to the vibration mode when the first natural frequency is equal to zero. Modal analysis can also be performed when the structure is completely unloaded, in which case the internal axial stresses in the elements do not exist and the geometric matrix is null, thus, where and are the vibration modes and natural frequencies for unloaded structure, respectively. It can be noted that the vibration modes of the structure vary between the vibration modes of the unloaded structure { } and the vibration modes of the structure when subjected to critical loading � � . By considering this aspect, in this work, it is proposed to approximate the vibration modes of the loaded structure � � � from a simple linear interpolation between { } and � � in the following way, Where ∆ is a parameter that indicates the internal load magnitude and varies between 0, when the structure is unloaded, and 1 when the structure is under the critical load. If the internal load is equal to the critical load, the vibration mode will be the buckling mode � �, while if the internal load is 0 then the vibration mode will be the vibration mode of unloaded structure { }.
If both modes, and , are normalized to the mass matrix, { } [ ]{ } = 1 and using the orthogonality property of vibration modes, the natural frequency of the loaded structure � can be approximated, such as, (2.41) If the dynamic equation is pre-multiplied by [ ] and vibration modes are normalized to the mass matrix, uncoupled motion equations for each vibration mode are obtained. Besides, for static loads, the acceleration vector is {̈} = 0. In this way, the following equation is given by,

Static displacement from vibration modes and frequencies
Or uncoupled equations, (2.48) Considering the contribution of " " vibration modes, the following formulation for the calculation of the static displacements is obtained, . (2.49)

Modal p-delta analysis
It was shown through Equations 2.48 and 2.49 that it is possible to approximate the vibration modes and the natural frequencies of the loaded structure by interpolating the vibration modes between the state without axial load and the state with critical axial load. Therefore, it is possible to calculate the displacements of the structure for a determined load level using Equation 2.48 that depends on the vibration modes and natural frequencies. These displacements consider the effects of NLG. Based on these concepts, the method Modal P-delta is proposed, which uses modal and buckling analysis to perform a NLG analysis. This method allows us to obtain the behavior of the equilibrium path until the buckling of the structure occurs.
In this way, the basic routine for performing NLG analysis using the Modal P-delta method is described in Figure 5. The methodology is divided into two phases, pre-calculation, and calculation phase. In the first phase, some input parameters are determined and in the second, the displacements of the structure for each loading step are calculated.
In the first phase, the number of vibration modes "nmod" and the number of load steps "nstep" must be initially defined. After, the global stiffness matrix of the structure, the nodal force vector, and the global mass matrix of the structure must be calculated; the boundary conditions of the structure in all these parameters must be considered. With this, modal analysis can be performed to calculate the vibration modes of the unloaded structure [ ] and then calculate the geometric matrix � � by using the displacement due to external load { }. With the geometric matrix can be performed a buckling analysis to calculate the buckling modes� � and eigenvalues { }. The second phase begins with the definition of the total stiffness matrix of the structure that considers the contribution of the elastic matrix and the geometric matrix. The geometric matrix is considered to vary linearly between the unloaded case and the case of the fully-loaded structure, using the factor (∆= / ), where k is the load step. For each load, the step must be calculated the following parameters: the vibration mode of the loaded structure using a linear interpolation between { } and � � and the natural frequency for loaded structure. Thus, using Equation 2.48, the displacement in the current load step can be calculated. The previous procedure can be performed by adding various vibration and buckling modes, as shown in Figure 5.

Analysis of a column
The purpose of this item is to use the proposed method, Modal P-delta, and apply it to a column as shown in Figure 2. The equilibrium path obtained with the Modal P-delta method will be compared with that obtained using the traditional iterative P-delta method [3]. Thus, the characteristics of the column are defined as: a length of 3 , a cross-section of 0.04 ², a specific mass of 250 / 3, Young's modulus of 2.72 × 10 10 / ² and a moment of inertia equal to 0.000133 4 . The number of load steps has been set as equal to 200 and there is a lateral load at the top of the element of 10 . The total axial load applied was 994.27 which corresponds to the critical buckling load of the element. For numerical analysis, the discretized model was defined with 20 elements. Figure 6 shows the displacements at the top of the column, using the iterative P-delta methods and the proposed method -varying the number of vibration and buckling modes. A simple analysis of this figure shows a total agreement between the results obtained with the delta-P method and the results obtained with the modal P-delta even when there is a variation in the amount of vibration and buckling modes. Additionally, it is important to report that the P-delta analysis demonstrated a computational processing time is approximately the tenth part of time spent by the traditional method. All analysis were developed on a computer with AMD Ryzen 7 1700 Eight-Core processor 3.00GHz, physical memory (RAM) 16 GB.   Figure 7 shows the error of the Modal P-delta method, taking as reference the analytic P-delta analysis, which formulation is in the Appendix A. The error consists of a difference between the displacements of the two methods. In the results of the modal P-delta method, the number of vibration and buckling modes was varied between 1 and 6 modes. The analysis of the error shows that, while is increasing the number of modes used, the error tends to decrease for the first load cases and converge to a constant value in the last ones, to this test the maximum average error is around 0.0002143m in the last case of load.
In general, it is noted that to the first load steps, the resulting error is less as the number of used modes increases, but to the last load steps, the increasing the number of modes makes the increasing of amount of error, which tends to be a constant value. However, despite the existence of error in the calculation of displacements in all situations, all of these are considered minimal, in the order of 1×10 -4 m.

Analysis of a plane frame
In this item, it is shown the applicability of the proposed modal P-delta method in a plane frame structure. The results were compared to the traditional iterative P-delta method.
Two analyses were performed on the same model, performing a modification of the external loads. The plane frame uses the properties shown in Table 1. The number of load steps used was equal to 200. The first type of applied loading, called test 01, were applied axial loads at the top in each of the columns P01, P02, and P03 with a load equal to 0.28, 0.5, and 0.28 of the calculated critical load of 8.215 × 10 7 , respectively. In test 01 there is also a lateral load of 100 , as shown in Figure 8. Test 02 was generated with vertical loads distributed in the floor beams and lateral loads on each of the floors, the loads are shown in Figure 8. The loads are amplified with a multiplier factor of 10.5, to bring the structure to a condition of critical buckling loading. In both tests, the equilibrium path was constructed extracting the displacements at the top of the P03 column, with the P-delta and Modal P-delta methods. In the case of the Modal P-delta method, the number of vibration and buckling modes used was varied. The results of the analyses in both load tests are presented in Figures 9 and 10. In these figures, it can be noted that comparatively there are minimal differences between the use of P-delta and Modal P-delta.
In Figures 9 and 10, the calculation times of both load tests with p-delta and modal p-delta methodologies are also reported, where is possible to notice the considerable time reduction of the computational calculations using the proposed method opposite to the traditional P-delta method.
It is important to say that the calculation times for the iterative P-delta method have been verified in commercial computational software such as SAP-2000©, having a consideration that, the number of elements of discretization in each model tested with the proposed methodology was the same used in a verification with a software SAP-2000, which were 2 elements for the columns and beams. This verification shows calculation times like or greater than those shown with our computational code. These notable computational calculation time differences between P-delta and the proposed method are due to which the traditional P-delta method calculates an inverse of the global stiffness matrix to calculate the displacements at each load step or uses some optimization method. In contrast, the proposed method performs a single-time modal and buckling analysis for calculating displacements in each load step, the rest of the calculation procedures are simple multiplications and sums of vectors with no iterations. Figures 11 and 12 shown the error of the P-delta modal analysis, considering as reference the iterative P-delta method. It can be noted, as in the case of the column, the error between the two methodologies is minor to the first steps loads also as the number of vibration modes are increased in the calculation, the error tends to converge to a constant value, to the test 01 the maximum relative error is around 0.02447m and to the test 02 the maximum relative error is around 0.071864m, and these values are almost the same from using 4 modes onwards, for both tests. However, for all situations of the number of vibration modes the maximum relative error was 1 × 10 −2 .   Step  As seen in the Figures 11 and 12, the first two or three vibration modes have a direct and rapid influence on the error reduction, however it is not yet clear the ideal number of vibration modes to use. Because the results show that with considering the first modes, the method proves to be efficient, with low error values, but when higher modes are added, the proposed method shows a relative stabilization of the error, without improvement, which allows conclude that the methodology cannot converge the error to zero with an increase in the number of vibration modes. This characteristic of the error is due to the approximation performed by the interpolation of the vibration modes (structure without load and at critical load) to establish the displacements of the structure. It is assumed that this interpolation has an intrinsic error in the interpolation of each vibration mode. It is necessary in future research with different structural problems to check the error behavior for different modes and try to understand what error is produced in the interpolation of the higher modes after the first load step. Another characteristic seen in the error graphs is the increase in error as the number of steps increases. This may be since in the upper steps there is an accumulation of error from the previous steps.
Finally, it is important to highlight that to make the interpolation of vibration modes and buckling modes, it is necessary to identify longitudinal modes which will be omitted in the calculation due these do not contribute to the lateral displacements. The identification of the longitudinal vibration modes was performed using the intensities of the amplitudes of the vibration modes, considering as a condition that the absolute value of the amplitudes of the transverse vibration modes are close to zero.

CONCLUSIONS
This paper describes a new methodology -Modal P-delta -for NLG analysis of structures. This method considers the relationship between dynamic parameters (vibration modes and natural frequencies) and geometrical NL effects. It has been shown that the vibration mode in the critical buckling load and the buckling mode are linearly dependent. In this way, the vibration modes of the loaded structure were approximated by performing a linear interpolation of vibration modes between these two states of load, without and critical buckling load. Furthermore, the natural frequencies were approximated using the orthogonal properties of the approximate vibration modes of the loaded structure. With the approximate vibration modes and frequencies of the loaded structure, NL displacements were calculated.
Several examples were simulated on regular structures in the plane, showing the effectiveness of the proposed method. A column and a plane frame with two load cases were used to build the pre-buckling equilibrium path. Compared with the traditional iterative P-delta method, the results showed maximum errors in the order of 1 × 10 −2 for the frames and 1 × 10 −4 for the column analyzed. The Modal P-delta method allowed adding a certain number of vibration modes in the analyses and it was found that in the cases analyzed, while major will be the number of modes used, the error will tend to be a constant value, in this way using a maximum of six vibration modes is enough to obtain satisfactory results. Therefore, the use of higher vibration modes in the analyzed cases does not result in an error reduction.
The proposed method can be classified as a direct method since it approximates the displacements without an iterative process. Therefore, a relevant feature of the proposed methodology, compared to the traditional iterative P-delta method, is the reduction of computational time. In the different simulations it is showed a great reduction in the computational time. This is due to a drastic reduction in complex calculations that uses the iterative P-delta method at each loading step to obtain the displacement, such as the inverse of matrices or optimization methods (e.g., Newton-Raphson).
The efficiency of the modal P-delta method in irregular three-dimensional structures will be tested shortly. These future analyses will aim to understand the behavior of the method when the superior vibration and buckling modes have a greater influence on the displacements. Another interesting idea would be to analyze the computational time testing the proposed methodology in larger dynamical NL problems.

ACKNOWLEDGE
This work was conducted with the economic support of the Social Demand Program -UNILA, Strict Sense Postgraduate Scholarships.