Research Article  Open Access
Identifying Critical Loads of Frame Structures with Equilibrium Equations in Rate Form
Abstract
Frame structures are widely used in engineering practice. They are likely to lose their stability before damage. As an indicator of loadcarrying 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 welldeveloped 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 largescale 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.
1. 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 loadcontrol 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 loadcontrol 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 steplength 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–10], and the corotational formulations [11–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 largescale frame structures, all components in these structures should be treated as corotational elements, which will eventually lead to a set of highdimension 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 largescale 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 largescale 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 loadcontrol 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.
2. Equilibrium Equations in Rate Form
In general, the nonlinear finite element analysis of structures leads to a set of nonlinear algebraic equilibrium equations, which yieldswhere is the vector of internal forces depending on the displacement . Under the condition of proportional loading, the magnitudes of the external forces are represented by a unit load incremental vector and a loadcontrol parameter , respectively [2].
Governed by (1), displacement 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 asIt 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 loadcarrying 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 loadcontrol 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 welldeveloped 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 loadcontrol 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 becomes singular.
(a)
(b)
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 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 and the singularity of the tangent stiffness matrix .
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 can be related to the existing displacement by Differentiating (6) yieldsIt indicates that if matrix 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 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 elementindependent for more extensive applications.
3. 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 wellknown 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 socalled elementindependent corotational formulation, and Crisfield and Moita [11] described a unified approach for geometric nonlinear corotational formulation of various elements.
3.1. 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 twonode beam element, a point on the left endsection can serve as the third point III, as shown in Figure 3(b)
(a)
(b)
Its origin is placed at point I and its base vectors are defined bywhere a vector with a symbol “~” on the top is its skew symmetric matrix. For example, if a vector , its skew symmetric matrix is expressed as
From (8), it can obtain the conclusion that is parallel to and has no component along , which yieldsIn 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 byBy definition, the local coordinate system angular velocity can be expressed asIt can be verified that the above equation can be expanded aswhere
The local coordinate system defined in this manner is elementindependent in the sense that it applies to substructure elements, so long as the displacement relative to the local coordinate system is small enough.
3.2. Relationship between Global and Local Nodal Parameters
In the local coordinate system, a node located at will occupy a position in the global coordinate system, which yieldswhere a vector or a matrix with a superscript “0” represents its initial values. The global and local displacement of the node is defined, respectively, byThe 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).
(a)
(b)
In accordance with (15), the local displacement can be extracted from the global one asTheir rates of change are given byIf 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 . The current orientation matrix of the cross section of node in the global coordinate system can be described aswhere 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 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 is given bywhere 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 thatwhere 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 represents the rotation of the coordinate system relative to the coordinate system , defined as such is the angular velocity expressed in the coordinate system . In terms of the rotational vector and its rate of change , the angular velocity can be written as [23]whereIts inverse is given byWhen 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 asOn the basis of (18), (29), and (13), relationships between the local and global virtual velocities can be given bywhere
According to (30) and (32), we can prove that all components of and the last two components of as well as the last component of 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 endsection of the beam, so is not independent but related to and by the following equation:
Moreover, the relationships described by (30)–(32) are elementindependent 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.
3.3. 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 aswhere is the set of all nodes in this substructure. When the local displacement and rotations are small enough, the internal forces and torques can be obtained through some adequate linear elements aswhere , , , and are the submatrices of stiffness matrix 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 aswhere is the density of surface forces, is the gravity acceleration, and is the mass density. According to (30), they can be rewritten aswhere is the total mass of the substructure, is the position of the mass center in the local coordinate system, and and are defined byBy virtue of element shape functions, can be written in the following form:and the integration in the expression of can be written aswhere and are the submatrices of gravity influence coefficient matrix of the substructure. The remaining parts of can be rewritten aswhere is the set of three nodes I, II, and III for defining the local coordinate system, and the generalized forces resulted from gravities and are given bywhere and are two coefficient matrices defined byIt can be seen that only affects the generalized forces of the nodes relevant to the local coordinate system, whereas 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 thatThe principle of virtual power can be expressed aswhere the summation range includes all substructures in the system . According to a standard process, the equilibrium equations can be derived from (45).
3.4. Static Condensation Procedure for a Substructure Element
Some insights can be revealed by the characteristics of the nodes in a substructure. The nodes in a substructure can be divided into three groups. The first group comprises the three nodes I, II, and III for defining the local coordinate system. The second group consists of the nodes carrying surface loads and the nodes in connection with other substructures, except for the three nodes I, II, and III. The union of these two groups is defined as boundary nodes. The remaining nodes in the substructure are sorted in the third group and defined as interior nodes , as shown in Figure 5.
Local displacement and rotational vectors of the nodes in the substructure can be regrouped into or according to the criterion of whether the node belongs to the set of boundary nodes or the set of interior nodes . Similarly, the stiffness matrix and gravity influence coefficient matrix of the substructure element can be written in blockpartitioned form as follows:After doing this, (44) can be rewritten in the formIn the above equation, is independent because the virtual powers of the remaining substructures have nothing to do with it, which yieldsWhen more than six components in are constrained, will be an invertible matrix. In this case, (47) can be simplified aswhere the equivalent stiffness matrixand the equivalent gravity influence coefficient matrixIt 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 aswhere is an arbitrary node in boundary node group . The generalized forces and torques in the local coordinate system can be expressed aswhere , , , and and , are the submatrices of and related to the translational and rotational displacement vectors.
Substituting (30)–(32) into (41) and (52), we can obtain thatIn this equation, the generalized forces and torques in the global coordinate system are given byfor the nodes in the set where
In the case that the substructure is actually a beam element, the virtual velocity can be determined by the virtual velocity of the local coordinate system origin and the virtual angular velocity of the left endsection represented by (33), and the force in (56) is carved up between and resulting in increments to them byMeanwhile, and should be updated to be and , 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 systemIn this equation, the generalized internal forces and torques 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.
4. 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 yieldswhere the derivatives of the local generalized internal nodal forces and torques with respect to the global parameters can be obtained byAccording to (43), it can be verified that, for the vector expressed in the local coordinate system , the following equations are valid:where and .
By the definitions of , , and as shown by (10), their derivatives are given byRepresenting in the above equation by (13) and substituting the resulting equation into (62) yield thatwhere the coefficients matrices are defined by
In accordance with (56), the derivatives of the global generalized internal and gravity forces , , and with respect to the global parameters are obtained by According to (57), the derivatives of and with respect to the global parameters in (66) can be given byThe 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.
5. 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 and Poisson’s ratio .
5.1. 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 m, its area is , its polar moment of inertia is , and its two moments of inertia are .
(a)
(b)
According to Euler’s formula of fixedfree beam under axial compression, the critical load of this simple structure is given byA small perturbation moment of 1 N·m is applied at the beam topend 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 elements are shown, respectively, in Figure 7, where the meaning of legends is defined by Figure 6(b).
(a)
(b)
The critical load can be detected by the maximum norm of the displacement derivatives with respect to axial force . 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 m corresponding to the obtained critical load, so the more precise critical load is given byIf 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.
5.2. 24Member Hexagonal StarShaped Shallow Dome
The starshaped 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 and shear modulus . The geometric parameters are , , , , and mass density . 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.
(a)
(b)
By using 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 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.
5.3. A Frame with One Weaker Beam
As shown in Figure 10, a frame is composed of four main beams and a subsidiary beam , 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 , polar moment of inertia , and two moments of inertia , 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 fixedfixed beam.
It can be easily predicated that the subsidiary beam will buckle first while the remaining parts remain stable, as illustrated in Figure 11.
When the subsidiary beam is modeled by only one element and the remaining beams are divided into twenty elements, the equilibrium path of 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.
(a)
(b)
It is predicted by this model that the critical load is 3514.966 kN, and 1372.245 kN is carried by the subsidiary beam . However, the first buckling load of the beam is only 74.142 kN according to Euler’s formula of fixedfixed 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 and 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 is taken as the identification parameter. However, the local buckling of does not lead to collapse of the frame. When we take as the identification parameter, the critical load of the frame will reach 1223.360 kN.
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.
5.4. 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.
(a)
(b)
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).
(a)
(b)
For all boundary nodes displacement values, the absolute values of all components of the vector , 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. Loaddisplacement 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.
6. Conclusion
As an indicator of loadcarrying 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 welldeveloped numerical methods for ordinary differential equations can be used to facilitate the structural stability analysis. Considering the geometrical nonlinear effects of largescale 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 elementindependent 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.
Notations
Base vectors of the global coordinate system  
:  Global position vectors of three noncollinear points I, II, and III for defining the substructure local coordinate system 
:  Base vectors of the substructure local coordinate system 
:  3 × 3 identity and zero matrices 
:  Rotational matrix and angular velocity of the substructure local coordinate system with respect to the global coordinate system 
:  Local and global initial position vectors of a node of a substructure 
:  Local and global current position vectors of a node of a substructure 
:  Local and global translational displacement vectors of a node of a substructure 
:  Local and global rotational vectors of a node of a substructure 
:  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 
:  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 and 
:  Group of the interior nodes of substructure 
:  Group of all nodes of substructure , which is the union of the three groups , , and 
:  Gravity acceleration in the global coordinate system 
:  Stiffness and gravity influence coefficient matrices of substructure 
Condensation stiffness and gravity influence coefficient matrices of substructure  
:  Generalized internal forces and torques in the local coordinate system of substructure 
:  Generalized internal forces and torques in the global coordinate system of substructure 
:  Position of the mass center in the local coordinate system 
:  Mass density of the frame structures. 
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research was supported by the China National Science Foundation under Grant no. 11372057.
References
 J. Shi, “Computing critical points and secondary paths in nonlinear structural stability analysis by the finite element method,” Computers & Structures, vol. 58, no. 1, pp. 203–220, 1996. View at: Publisher Site  Google Scholar
 P. Wriggers and J. C. Simo, “A general procedure for the direct computation of turning and bifurcation points,” International Journal for Numerical Methods in Engineering, vol. 30, no. 1, pp. 155–176, 1990. View at: Publisher Site  Google Scholar
 J.M. Battini, C. Pacoste, and A. Eriksson, “Improved minimal augmentation procedure for the direct computation of critical points,” Computer Methods in Applied Mechanics and Engineering, vol. 192, no. 1618, pp. 2169–2185, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 R. Kouhia, M. Tůma, J. Mäkinen, A. Fedoroff, and H. Marjamäki, “Implementation of a direct procedure for critical point computations using preconditioned iterative solvers,” Computers & Structures, vol. 108109, pp. 110–117, 2012. View at: Publisher Site  Google Scholar
 M. A. Crisfield, “A fast incremental/iterative solution procedure that handles ‘snapthrough’,” Computers & Structures, vol. 13, no. 1–3, pp. 55–62, 1981. View at: Publisher Site  Google Scholar
 M. A. Crisfield, “An arclength method including line searches and accelerations,” International Journal for Numerical Methods in Engineering, vol. 19, no. 9, pp. 1269–1289, 1983. View at: Publisher Site  Google Scholar
 M. RezaieePajand and J. Alamatian, “Automatic DR structural analysis of snapthrough and snapback using optimized load increments,” Journal of Structural Engineering, vol. 137, no. 1, pp. 109–116, 2011. View at: Publisher Site  Google Scholar
 K. J. Bathe, E. Ramm, and E. L. Wilson, “Finite element formulations for large deformation dynamic analysis,” International Journal for Numerical Methods in Engineering, vol. 9, no. 2, pp. 353–386, 1975. View at: Publisher Site  Google Scholar
 Y. B. Yang, S. P. Lin, and L. J. Leu, “Solution strategy and rigid element for nonlinear analysis of elastically structures based on updated Lagrangian formulation,” Engineering Structures, vol. 29, no. 6, pp. 1189–1200, 2007. View at: Publisher Site  Google Scholar
 G. R. Joldes and B. Zwick, “Numerical simulation of anisotropic tissue growth using a total Lagrangian formulation,” Mathematical Problems in Engineering, vol. 2015, Article ID 674639, 5 pages, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 M. A. Crisfield and G. F. Moita, “A unified corotational framework for solids, shells and beams,” International Journal of Solids and Structures, vol. 33, no. 20–22, pp. 2969–2992, 1996. View at: Publisher Site  Google Scholar
 C. A. Felippa and B. Haugen, “A unified formulation of smallstrain corotational finite elements: I. Theory,” Computer Methods in Applied Mechanics and Engineering, vol. 194, no. 21–24, pp. 2285–2335, 2005. View at: Publisher Site  Google Scholar
 M. H. Tsai, W. Y. Lin, Y. C. Zhou, and K. M. Hsiao, “A corotational finite element method combined with floating frame method for large steadystate deformation and free vibration analysis of a rotatinginclined beam,” Mathematical Problems in Engineering, vol. 2011, Article ID 146505, 29 pages, 2011. View at: Publisher Site  Google Scholar
 J. Rungamornrat, S. Sihanartkatakul, and P. Kanchanakitcharoen, “Adaptive, smallrotationbased, corotational technique for analysis of 2D nonlinear elastic frames,” Mathematical Problems in Engineering, vol. 2014, Article ID 640983, 15 pages, 2014. View at: Publisher Site  Google Scholar
 J.M. Battini and C. Pacoste, “Corotational beam elements with warping effects in instability problems,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 1718, pp. 1755–1789, 2002. View at: Publisher Site  Google Scholar
 G. Wang, Z. Qi, and X. Kong, “Geometrical nonlinear and stability analysis for slender frame structures of crawler cranes,” Engineering Structures, vol. 83, pp. 209–222, 2015. View at: Publisher Site  Google Scholar
 Y. He, X. Zhou, and P. Hou, “Combined method of super element and substructure for analysis of ILTDBS reticulated megastructure with singlelayer latticed shell substructures,” Finite Elements in Analysis and Design, vol. 46, no. 7, pp. 563–570, 2010. View at: Publisher Site  Google Scholar
 J. L. Meek and Q. Xue, “A study on the instability problem for 3D frames,” Computer Methods in Applied Mechanics and Engineering, vol. 158, no. 34, pp. 235–254, 1998. View at: Publisher Site  Google Scholar
 J. Baumgarte, “Stabilization of constraints and integrals of motion in dynamical systems,” Computer Methods in Applied Mechanics and Engineering, vol. 1, pp. 1–16, 1972. View at: Publisher Site  Google Scholar  MathSciNet
 J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, 2008. View at: Publisher Site  MathSciNet
 L. F. Shampine and M. W. Reichelt, “The MATLAB ODE suite,” SIAM Journal on Scientific Computing, vol. 18, no. 1, pp. 1–22, 1997. View at: Publisher Site  Google Scholar  MathSciNet
 C. C. Rankin and F. A. Brogan, “An element independent corotational procedure for the treatment of large rotations,” Journal of Pressure Vessel Technology, vol. 108, no. 2, pp. 165–174, 1986. View at: Google Scholar
 J. H. Argyris, “An excursion into large rotations,” Computer Methods in Applied Mechanics and Engineering, vol. 32, no. 1–3, pp. 85–155, 1982. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2015 Zhaohui Qi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.