Identifying Critical Loads of Frame Structures with Equilibrium Equations in Rate Form

Frame structures are widely used in engineering practice. They are likely to lose their stability before damage. As an indicator of load-carrying capability, the first critical load plays a crucial role in the design of such structures. In this paper, a new method of identifying this critical load is presented, based on the governing equations in rate form. With the presented method, a great deal of well-developed numerical methods for ordinary differential equations can be used. As accurate structural tangent stiffness matrices are essential to stability analysis, the method to obtain them systematically is discussed. To improve the computational efficiency of nonlinear stability analysis in large-scale frame structures, the corotational substructure elements are formulated as well to reduce the dimension of the governing equations. Four examples are studied to illustrate the validity and efficiency of the presented method.


Introduction
Frame structures are widely used in engineering practice and are progressively optimized towards sparsity and slenderness.These frame structures often show strong nonlinear behaviors after carrying heavy loads.Instability is one of the causal factors for damages and injuries, which results in permanent disability.In general, the more slender such structures become, the more likely they are to lose stability and even result in damage.To improve the stable loadcarrying capacity, a simplified and efficient computational method for determining the first critical loads is required for the structural design, application, and safety of these structures.
In the stability analysis of structures, the load on the structure is usually assumed to be varied in a proportional manner; that is, it can be represented by a single loadcontrol parameter.In this case, identifying critical loads requires the solution of the equilibrium equations under the variation of the load-control parameter and the detection of critical points where the tangent stiffness matrix becomes singular.Generally, the critical points are detected by two main processes called direct and indirect methods [1].
In the direct method, the equilibrium equations are augmented with a criticality condition [2,3].Since the augmented equations are highly nonlinear, a nice preconditioner [4] is essential for the successful application of the method.In the indirect method, the behavior of the structure is monitored by following the equilibrium path and detecting the singularity of the tangent stiffness matrix associated with each point on the path.When the research target is the first critical load, the load-control parameter is monotone increasing.If all critical loads are being sought, an arclength technique is presented [5].The major advantages of this technique are its better convergence property and the potential to reveal the whole picture of the structural stability.Although the path following methods combined with arclength techniques are the most popular methodology in the stability analysis of frame structures, they still suffer from the difficulties in adjusting the step-length adequately to ensure the desired accuracy [5,6].
Both of the direct and indirect methods need the tangent stiffness matrix which is not easy to obtain when the effect of geometric nonlinearity is considered.Such difficulties can be avoided by the Dynamic Relaxation (DR) method, in which the static system is treated as a fictitious vibration system through adding fictitious inertia and damping forces to the static equilibrium equations [7].However, the accuracy and efficiency of the DR method are decreased compared with the methods using the tangent stiffness matrix.
A proper and efficient finite element model is necessary for describing the geometrical nonlinear behaviors of the frame structures.Three kinematic descriptions are used for the finite element analysis of nonlinear structural mechanics: the total Lagrangian formulations, the updated Lagrangian formulations [8][9][10], and the corotational formulations [11][12][13][14].Specifically, the corotational formulation is well suitable for describing the geometrical nonlinear behaviors of the frame structures, since it leads to an artificial separation of the geometric nonlinearities [12].The corotational formulation has been adopted by many researchers for the geometric nonlinearity and stability analyses of frame structures [15,16].
However, when the corotational formulation is adopted for modeling the large-scale frame structures, all components in these structures should be treated as corotational elements, which will eventually lead to a set of high-dimension nonlinear equilibrium equations to be solved.To reduce the dimension of nonlinear equilibrium equations, the largescale structure can be divided into several substructures [17].Meanwhile, a corotational substructure element formulation is required to improve the geometrical nonlinear and stability analysis of large-scale frames, based on the corotational method and the substructure method.By applying substructure elements with the static condensation technique, a large number of internal nodal variables of the structures are controlled by their global variables.Consequently, the internal nodal displacement of each substructure is removed, resulting in the fact that the total number of degrees of freedom decreases consequentially.However, the static condensation technique only applies to linear systems and cannot be directly used for the geometrical nonlinear and stability analysis of large-scale frame structures.Therefore, the static condensation must be implemented within the local coordinate system, and structural gravity needs to be considered in this static condensation procedure.
The paper is organized as follows.In Section 2, a new method of identifying the first critical load of frame structures is presented by solving the rate form of the equilibrium equations and detecting the singularity of the tangent stiffness matrices.Such differential equations can be obtained by differentiating the equilibrium equations with respect to the load-control parameter.In Section 3, the corotational substructure elements are formulated to make the presented method available in the nonlinear stability analysis of largescale frame structures.In consideration of gravity forces, a static condensation technique for corotational substructure elements is implemented under the local coordinate system.As the accurate tangent stiffness matrix is an essential ingredient in this method, Section 4 is devoted to the systematic formulation of this matrix for a general frame structure, which is given by taking the derivative of the generalized nodal forces with respect to the global parameters.In Section 5, four numerical examples are given to illustrate the validity and efficiency of the presented method.

Equilibrium Equations in Rate Form
In general, the nonlinear finite element analysis of structures leads to a set of nonlinear algebraic equilibrium equations, which yields where f is the vector of internal forces depending on the displacement u.Under the condition of proportional loading, the magnitudes of the external forces are represented by a unit load incremental vector q 0 and a load-control parameter , respectively [2].
Governed by (1), displacement u varies with the parameter , and its rates are determined by the derivatives of (1): where the dot over a variable represents its derivative with respect to .
Theoretically, ( 1) and ( 2) should be satisfied simultaneously, but this is hard to achieve in reality.In analogy to the constraint stabilization method [19], the two individual equations can be integrated into the following equation: where  is called a stabilization factor.In terms of the components of , its solution can be expressed as It indicates that the violations in  and δ would be decreased exponentially when the displacement along the equilibrium path is solved from (3); that is, As an indicator of load-carrying capacity, the critical load associated with the first limit or bifurcation point on the equilibrium path plays a crucial role in the structural design.In this paper, we focus on the method for identifying such critical loads and suggest using (5) instead of (1) as the governing equation.In this case, the load-control parameter  grows continually until the first critical point is arrived at, and the parameter  in (5) should be a positive constant.
From a mathematical point of view, following the equilibrium paths is the same thing as describing the history of the solutions varying with a parameter, which can be achieved better by the differential equations with the same solutions [20].One of the benefits of transforming the algebraic governing equation to its rate form resides in that a great deal of well-developed numerical methods can be used for reference, especially the time step adjustment procedure.With a suitable solver of ordinary differential equations, such as Matlab ode15s [21], we can pass over the smooth path segments with few solution points, whereas the steep path segments can be described in detail by more solution points.Another benefit of such a transformation is the facilitation for identifying the critical loads.As is well known, at limit points on the equilibrium path, the derivatives of displacement with respect to the load-control parameter are infinite, as shown in Figure 1(a), and these derivatives become nonunique at bifurcation points, as shown in Figure 1(b).At limit or bifurcation points, the tangent stiffness matrix [f/u] becomes singular.
In this case, if the right side of (5) happens to be a linear combination of column vectors of the tangent stiffness matrix, the derivatives of displacement du/d might not be large.This is just the condition for the presence of bifurcation points.However, such a condition is so frail that it can be destroyed by a tiny random perturbation in the load or numerical errors in the solution process.
In engineering practice, when the point in question is to find the critical load, the types of critical points are less important, and we can identify the critical load by monitoring the rapid changes in the derivatives of displacement du/d and the singularity of the tangent stiffness matrix [f/u].
It is worth noting that there is always the danger of omitting the local buckling when finite element methods are applied in the critical load analysis.This is chiefly because any finite element model of structures includes the displacement of only part of the points in the structure.Assuming some new points between the element nodes are added in the model as shown in Figure 2, without loss of generality, their displacement u in can be related to the existing displacement u by Differentiating ( 6) yields It indicates that if matrix [/u in ] becomes singular, the local buckling will occur.
In a frame structure, the shorter a beam element is, the greater the required force for the local buckling would be.So the elements should be divided fine enough in order to prevent the local buckling force from being smaller than the critical load.
From the above discussion, it can be seen that the tangent stiffness matrix [f/u] plays a crucial role in identifying the critical load.It has to be displacement dependent so that the features of geometric nonlinearity can be represented adequately.Moreover, the method to obtain it should be element-independent for more extensive applications.

Corotational Formulation of a Substructure Element
In general, the maximum load sustained by a frame structure can be determined by two specific loads.The first one, called the failure load, is the load that can make some components broken.The second one, called the critical buckling load, is the load that can turn the structure from stable to unstable states while structural strains remain small.The latter is more important to structural design.If this load is less than the failure load, the structure would be not optimal.Consequently, we can assume that the critical buckling load is mainly caused by geometric nonlinearity.In this case, the corotational formulation is more desirable among many methods of representing the influence of geometric nonlinearity in virtue of its acceptance of abundant linear elements [15].The corotational formulation in geometrical nonlinear problems is a well-known technique.Its most important advantage is that it leads to artificial separation of the geometric nonlinearity, and a linear strain definition in the local coordinate system is used [14,15].Using this property, Rankin and Brogan [22] introduced the so-called element-independent corotational formulation, and Crisfield and Moita [11] described a unified approach for geometric nonlinear corotational formulation of various elements.

Local Coordinate System
. A local coordinate system of a substructure element can be established with the positions of three noncollinear points, denoted by I, II, and III, respectively, as shown in Figure 3(a).When the substructure is actually a beam, such as a two-node beam element, a point on the left end-section can serve as the third point III, as shown in Figure 3(b) Its origin is placed at point I and its base vectors are defined by where a vector with a symbol "∼" on the top is its skew symmetric matrix.For example, if a vector a = ( 1  2  3 ) T , its skew symmetric matrix is expressed as From ( 8), it can obtain the conclusion that r II −r I is parallel to e 1 and r III − r I has no component along e 3 , which yields In pace with the motions of the three points, the local coordinate system will rotate relative to the global coordinate system, and its orientation matrix is given by By definition, the local coordinate system angular velocity  cr can be expressed as It can be verified that the above equation can be expanded as where The local coordinate system defined in this manner is element-independent in the sense that it applies to substructure elements, so long as the displacement relative to the local coordinate system is small enough.

Relationship between Global and Local Nodal Parameters.
In the local coordinate system, a node  located at r  will occupy a position in the global coordinate system, which yields where a vector or a matrix with a superscript "0" represents its initial values.The global and local displacement of the node is defined, respectively, by The initial and current states of a substructure element are shown in Figure 4(a).When the substructure is actually a beam, such states are shown in Figure 4(b).
In accordance with (15), the local displacement can be extracted from the global one as Their rates of change are given by If the node  belongs to a substructure, it is associated with a cross section whose orientation in the global coordinate system can be represented by an orthogonal matrix A  .The current orientation matrix of the cross section of node  in the global coordinate system can be described as where S 0 is the global initial orientation matrix of the cross section of the node , which is a constant;   is the global rotational vector of node  in the sense that it is independent of the element.In terms of a rotational vector , we can express an orthogonal matrix as [23] where  is the norm of  and I is the identity matrix.In the case that the element undergoes a large global rotation, the norm of   is also large, and the expression of the corresponding rotational matrix cannot be simplified.Another way to describe the orientation matrix A  is given by where S 0 is the local initial orientation matrix of the cross section of the node , which is a constant;   is the local rotational vector of node , and its norm is small enough to make the linear elements survive in the geometrical nonlinear situation.So the rotational matrix can be simplified as From ( 19) and ( 21), one can conclude that where   and   are the rotational vectors of node  expressed in the local and global coordinate systems, respectively.When the relative displacement and rotations of element nodes are all small enough in a substructure, we can attach a local coordinate system to this substructure and treat the whole substructure as a single generalized element.In this case, ( 23) is still valid.
The angular velocity of a rotation can be defined by the following equation: If R represents the rotation of the coordinate system  1 relative to the coordinate system  0 ,  defined as such is the angular velocity expressed in the coordinate system  0 .In terms of the rotational vector  and its rate of change θ , the angular velocity  can be written as [23] where Its inverse is given by When the norm of the rotational vector  is small enough, the two matrices can be approximated as According to ( 23)-( 24), the relationship between the local angular velocity   and global angular velocity   can be described as On the basis of ( 18), (29), and (13), relationships between the local and global virtual velocities can be given by where According to (30) and (32), we can prove that all components of  u I and the last two components of  u II as well as the last component of  u III are zeroes, which is consistent with the definition of the local coordinate system.
In the case that the substructure is actually a beam element, the third point in the definition of the local coordinate system is located on the left end-section of the beam, so  u III is not independent but related to  u I and  I by the following equation: Moreover, the relationships described by (30)-(32) are element-independent and valid in a substructure.Since a beam element can be viewed as a specific substructure, in the following section, we will study how to obtain the nonlinear equilibrium equations for a substructure by means of these relationships.

Virtual Powers of Internal and External
Forces for a Substructure Element.As rigid motions play no role in internal forces, the virtual power of internal forces for a substructure   can be expressed in the local coordinate system as where   is the set of all nodes in this substructure.When the local displacement and rotations are small enough, the internal forces f  and torques m  can be obtained through some adequate linear elements as where K   , K   , K   , and K   are the submatrices of stiffness matrix K  of the substructure.
In general, the external forces of substructures are composed of surface and body forces.Concentrated forces are regarded as the specific surface forces, and in most cases, body forces are gravity forces.Their contributions to the virtual power of external forces, denoted by    and    , respectively, are calculated in the global coordinate system and represented as where t  is the density of surface forces, g is the gravity acceleration, and  is the mass density.According to (30), they can be rewritten as where  is the total mass of the substructure, r  is the position of the mass center in the local coordinate system, and t  and g are defined by By virtue of element shape functions,    can be written in the following form: and the integration in the expression of    can be written as where W   and W   are the submatrices of gravity influence coefficient matrix W  of the substructure.The remaining parts of    can be rewritten as where   cr is the set of three nodes I, II, and III for defining the local coordinate system, and the generalized forces f   resulted from gravities and are given by where Γ II and Γ III are two coefficient matrices defined by ] ; It can be seen that   cr only affects the generalized forces of the nodes relevant to the local coordinate system, whereas   in can be treated together with the virtual power of internal forces due to their same expression form.
By means of (34) and (40), one can conclude that The principle of virtual power can be expressed as where the summation range includes all substructures in the system .According to a standard process, the equilibrium equations can be derived from (45).sorted in the third group and defined as interior nodes   in , as shown in Figure 5.

Static Condensation
Local displacement and rotational vectors of the nodes in the substructure   can be regrouped into u  or u in according to the criterion of whether the node belongs to the set of boundary nodes    or the set of interior nodes   in .Similarly, the stiffness matrix and gravity influence coefficient matrix of the substructure element can be written in block-partitioned form as follows: After doing this, (44) can be rewritten in the form In the above equation,  u in is independent because the virtual powers of the remaining substructures have nothing to do with it, which yields When more than six components in u  are constrained, K  will be an invertible matrix.In this case, (47) can be simplified as where the equivalent stiffness matrix and the equivalent gravity influence coefficient matrix It is indicated by (49) that only the displacement and rotational vectors of the boundary nodes are required to form the equilibrium equations.Owing to such a feature, we can greatly reduce the unknowns to be solved.In expanded form, (49) can be written as where  is an arbitrary node in boundary node group    .The generalized forces and torques in the local coordinate system can be expressed as Substituting (30)-( 32) into (41) and (52), we can obtain that In this equation, the generalized forces and torques in the global coordinate system are given by for the nodes in the set where In the case that the substructure   is actually a beam element, the virtual velocity  u III can be determined by the virtual velocity of the local coordinate system origin  u I and the virtual angular velocity of the left end-section  I represented by (33), and the force f III in (56) is carved up between f I and m I resulting in increments to them by Meanwhile, f I and m I should be updated to be f I + Δf I and m I + Δm I , respectively.Substituting (54) and (39) into (45) yields the system virtual power equations in terms of the displacement and rotational vectors of the boundary nodes of each substructure in the system In this equation, the generalized internal forces f  and torques m  are nonlinear functions of the displacement and rotational vectors, even though they originate from linear elements.Their partial derivatives with respect to the global displacement or rotational vectors compose the tangent stiffness matrix that plays a crucial role in the identification of critical loads.

Tangent Stiffness Matrices
Formulations of the tangent stiffness matrices can be obtained by taking the derivatives of the generalized internal nodal forces and torques in (55) with respect to the global parameters, which yields According to (43), it can be verified that, for the vector expressed in the local coordinate system m Σ =  Σ,1 e 1 +  Σ,2 e 2 +  Σ,3 e 3 , the following equations are valid: where By the definitions of  0 ,  1 , and  2 as shown by (10), their derivatives are given by Representing  cr in the above equation by ( 13) and substituting the resulting equation into (62) yield that where the coefficients matrices are defined by ; In accordance with (56), the derivatives of the global generalized internal and gravity forces f I , f II , and f III with respect to the global parameters are obtained by According to (57), the derivatives of f Σ and m Σ with respect to the global parameters in (66) can be given by The Jacobian matrix of the whole system can be obtained by assembling the tangent stiffness matrix of each substructure.However, the Jacobian matrix so obtained may not be the final result because of the constraints among virtual velocities arising from boundary conditions of the structure.Some modifications are required in these cases.Boundary conditions have to be treated in accordance with their details.

Numerical Examples
In this section, four examples are studied to illustrate the validity and efficiency of the proposed method.Except for Section 5.2, all the members in the remaining three examples are made of the material with elasticity modulus  = 2.1 × 10 11 Pa and Poisson's ratio ] = 0.3.

A Slender Axial Compression Beam.
The first example is a slender beam clamped at its bottom and subjected to an axial force  at its top, as shown in Figure 6(a).The length of the beam is  = 10 m, its area is  = 6×10 −3 m 2 , its polar moment of inertia is  = 5.4020 × 10 −5 m 4 , and its two moments of inertia are   =   = 2.7010 × 10 −5 m 4 .According to Euler's formula of fixed-free beam under axial compression, the critical load of this simple structure is given by A small perturbation moment of 1 N⋅m is applied at the beam top-end to make the detection of the critical points easier, as shown in Figure 6(a).The equilibrium paths obtained with the presented method by using  ∈ [2,4,8,16] elements are shown, respectively, in Figure 7, where the meaning of legends is defined by Figure 6(b).
The critical load can be detected by the maximum norm of the displacement derivatives with respect to axial force d/d.When this norm reaches the given threshold, the associated load is taken as the critical load.In this example, the threshold is set to be 0.001, which implies that 1 N increment in the load would result in 1 mm increment in the horizontal displacement.The obtained critical loads are compared with the first Euler buckling load in Table 1.
It can be seen that, with the increasing number of divided elements, the relative error becomes smaller, and once the number of elements exceeds four, the relative error is less than 1%, which is acceptable in engineering practice.
However, Euler buckling load is not the analytical solution, because the effect of axial displacement is ignored.In this example, the axial displacement of the beam top is near Table 2: The comparison of the first critical loads from the presented method and Meek and Xue [18].
Meek and Xue [18] The number of elements per member  = 0.0012 m corresponding to the obtained critical load, so the more precise critical load is given by If the buckling is not taken into consideration, such an axial displacement would require 150.420 kN compression force, which implies a dramatic reduction in the tensile stiffness when buckling occurs.This example also indicates that the presented method can identify the critical load with satisfying accuracy and the relative errors are decreased with the number of elements increasing.

24-Member Hexagonal Star-Shaped Shallow Dome.
The star-shaped shallow dome has been studied by a number of researchers as a benchmark problem due to its obvious geometrical nonlinear characteristics before collapse [18].The dome is made of the material with elasticity modulus  = 3.030 × 10 9 Pa and shear modulus  = 1.096 × 10 9 Pa.The geometric parameters are  = 3.170 × 10 −4 m 2 ,   = 2.377 × 10 −8 m 4 ,   = 2.95 × 10 −9 m 4 ,  = 9.18 × 10 −9 m 4 , and mass density  = 5 × 10 4 N/m 3 .The structure is pinned and restrained against translational motion.The load case is shown in Figure 8, where the node  is under a vertical load  along the negative -axis direction.
By using  ∈ [2,4,8,16] elements for each member, the equilibrium paths of node  obtained by the presented method are shown, respectively, in Figure 9.
In this example, the derivative |d  /d| of the vertical displacement at the apex of the structure is taken as the identification parameter, and the threshold is set to 0.1, which implies that 1 N increment in the load results in 100 mm increment in the vertical displacement.The first critical loads are compared with the results from Meek and Xue [18] in Table 2.
This example indicates that the presented method can also adapt to the space frame with limit point, and the relative errors are decreased with the number of elements increasing.Figure 8: Geometry of 24-member star-shaped shallow dome.

A Frame with One Weaker Beam.
As shown in Figure 10, a frame is composed of four main beams and a subsidiary beam  0 , carrying a vertical load  at the middle point of its top line.The constitutive constants of the main beams are the same as those in the first example, whereas those of the subsidiary beam are  = 7.0686 × 10 −4 m 2 , polar moment of inertia  = 3.6226 × 10 −7 m 4 , and two moments of inertia   =   = 1.8113 × 10 −7 m 4 , which are much weaker than the others in the buckling resistance.
A small perturbation moment of −1 N⋅m is applied at the middle point of its top line, as shown in Figure 10.Because all beams in this model are clamped to each other, each component can be approximately equivalent to a fixed-fixed beam.
It can be easily predicated that the subsidiary beam  0 will buckle first while the remaining parts remain stable, as illustrated in Figure 11.
When the subsidiary beam  0 is modeled by only one element and the remaining beams are divided into twenty  elements, the equilibrium path of  1 is shown in Figure 12, where the meaning of legend is defined in Figure 11.In this example, the threshold is set to 0.001.It is predicted by this model that the critical load is 3514.966kN, and 1372.245kN is carried by the subsidiary beam  0 .However, the first buckling load of the beam  0 is only 74.142 kN according to Euler's formula of fixed-fixed beam under axial compression.The capacity of carrying axial loads for the subsidiary beam is overestimated in this model.
When all beams are divided into twenty elements, the corresponding equilibrium paths of  1 and  2 are shown in Figure 13, where the meaning of legends is defined in Figure 11.It is seen from the corresponding equilibrium paths that the critical load of the frame is only 75.660 kN when d 2 /d is taken as the identification parameter.However, the local buckling of  0 does not lead to collapse of the frame.When we take d 1 /d as the identification parameter, the critical load of the frame will reach 1223.360kN.The results of this example imply that each component in a structure has to be modeled by enough elements to capture the possible local buckling and identify the critical load precisely.

A Slender Space
Frame Structure.The forth example is a slender space frame structure clamped at its bottom and carrying axial load  equally distributed on the four nodes of its top section, as shown in Figure 14(a).The slender frame structure is successively combined with 10 periodic insert sections.The resulting outline of this frame structure is a cuboid, whose length, width, and height are 60 m, 2 m, and 2 m, respectively.Each insert section is composed of four chord members and a series of web members, whose geometric parameters are depicted in Figure 14(b).The cross section of its chord members is a hollow circular shape with inner diameter 0.174 m and outer diameter 0.194 m, and the cross sections of its web members are 0.079 m and 0.089 m.All the members are made of the material with mass density 7.85 × 10 3 kg/m 3 .The gravity acceleration is 9.807 m/s 2 , which is parallel to the negative -axis.A small perturbation moment of 10 N⋅m is applied on the two nodes  and  (Figure 14(a)) of its top section to make the detection of the critical load easier.
Along the longitudinal direction, each insert section is modeled as one corotational substructure element in the presented method to reduce the dimension of the governing equations, as shown in Figure 15(a).Consequently, the whole frame structure is divided into ten same corotational substructure elements, as shown in Figure 15(b).
For all boundary nodes displacement values, the absolute values of all components of the vector du/d, as well as the singularity of the tangent stiffness matrices, are taken as the identification parameters, from which the first critical load is identified as 3212 kN.
Under the first critical load, contour plot of the displacement vector sum of all nodal solutions is shown in Figure 16 (it is noted that all displacement values are under 20 times' magnification).
The structure is studied with ANSYS software, and Beam188 element is used for modeling each member of the frame structure in the ANSYS model.NLGEOM command in ANSYS software is selected to consider geometrical nonlinear behavior of the structure, and the number of loading substeps is chosen as 45.Load-displacement curves of the point  (Figure 16) obtained in this paper are compared with the ANSYS results, as shown in Figure 17.In the ANSYS model, 1680 nonlinear equations are solved to identify the first critical load.Only 264 nonlinear equations are required by using the corotational substructure element formulation in the presented method.The corotational substructure element formulation improves the numerical efficiency by means of reducing the degrees of freedom.Meanwhile, the accuracy of the presented procedure is acceptable while comparing with the ANSYS results.The unknowns in the governing equations are greatly reduced with the proposed method.

Conclusion
As an indicator of load-carrying capability, the first critical load plays a crucial role in the structural design.This critical load can be identified by the path following method presented in this paper on the basis of governing equations in rate form.In doing so, the well-developed numerical methods for ordinary differential equations can be used to facilitate the structural stability analysis.Considering the geometrical nonlinear effects of large-scale frame structures, a corotational formulation for substructure elements is presented, by which the number of the governing equations can be greatly reduced.The accurate structural tangent stiffness matrices required in the method can be obtained systematically through the proposed element-independent corotational formulation.The presented method can be extended to the fields of mechanical and civil engineering, and so forth, such as the slender lattice boom and tower and jib structures of cranes and derricks.g: Gravity acceleration in the global coordinate system K  , W  :

Notations
Stiffness and gravity influence coefficient matrices of substructure   K sub , W sub : Condensation stiffness and gravity influence coefficient matrices of substructure   f  , m  : Generalized internal forces and torques in the local coordinate system of substructure   f  , m  : Generalized internal forces and torques in the global coordinate system of substructure   r  : Position of the mass center in the local coordinate system : Mass density of the frame structures.

Figure 1 :
Figure 1: (a) Limit points on the equilibrium path; (b) bifurcation point on the equilibrium path.

Figure 3 :
Figure 3: Local coordinate systems for (a) a substructure element and (b) a beam element.

Figure 4 :
Figure 4: Corotational kinematic descriptions for (a) a substructure element and (b) a beam element.

Figure 5 :
Figure 5: Three groups of nodes in the substructure   .
submatrices of K sub and W sub related to the translational and rotational displacement vectors.

Figure 6 :
Figure 6: (a) Axial compression beam and (b) its model with  elements.

Figure 7 :
Figure 7: Load-displacement curves of the nodes on the beam when the beam is modeled by different numbers of elements.
Vertical loading F (N)Using n = 2 elements per member Using n = 4 elements per member Using n = 8 elements per member Using n = 16 elements per member Meek and Xue[18]

Figure 9 :FFigure 10 :
Figure 9: Load-displacement curves of the node  by using  elements for each member.

Figure 12 :Figure 13 :
Figure 12: Load-displacement curve of  1 when  0 is modeled by one element.

Figure 14 :
Figure 14: (a) Slender space frame structure carrying an axial load; (b) geometric parameters of the periodic insert section.

Figure 15 : 5 Figure 16 :
Figure 15: (a) Corotational substructure element of an insert section; (b) substructure model of the slender frame structure.

Figure 17 :
Figure 17: Load-displacement curves of the point .

Table 1 :
The comparison between the first Euler buckling load and the critical loads obtained from presented method.

g 1
, g 2 , g 3 : Base vectors of the global coordinate system r I , r II , r III : Global position vectors of three noncollinear points I, II, and III for defining the substructure local coordinate system e 1 , e 2 , e 3 : Base vectors of the substructure local coordinate system I, 0: 3× 3 identity and zero matrices A cr ,  cr : Rotational matrix and angular velocity of the substructure local coordinate system with respect to the global coordinate system r 0  , r 0  : Local and global initial position vectors of a node  of a substructure r  , r  : Local and global current position vectors of a node  of a substructure u  , u  : Local and global translational displacement vectors of a node  of a substructure   ,   : Local and global rotational vectors of a node  of a substructure R   , R   : Local and global rotational matrix of a node  of a substructure   ,   : Local and global angular velocities of a node  of a substructure   : A substructure in the frame system    cr : Group of the three nodes I, II, and III for defining the local coordinate system of substructure      : Group of the interface nodes of substructure   , except for the three nodes in the definition of the local coordinate system    : Group of the boundary nodes of substructure   , which is the union of the two groups   cr and      in : Group of the interior nodes of substructure     : Group of all nodes of substructure   , which is the union of the three groups   cr ,    , and in Mathematical Problems in Engineering