Mathematical Problems in Engineering
Volume 2015 (2015), Article ID 325969, 17 pages
http://dx.doi.org/10.1155/2015/325969

Research Article

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

Zhaohui Qi,1,2 Xianchao Kong,1,2 and Gang Wang1,2

¹State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian 116024, China
²Department of Engineering Mechanics, Dalian University of Technology, Dalian 116024, China

Received 26 April 2015; Accepted 8 July 2015

Academic Editor: Roman Lewandowski

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.

Abstract.

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.

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 load-carrying 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 load-control 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 arc-length 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 arc-length 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–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 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 large-scale 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 large-scale 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 yields

where 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 load-control 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 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 becomes singular.

(a)

(b)

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

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) yields

It indicates that if matrix becomes singular, the local buckling will occur.

Figure 2: Local buckling of a frame.

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 element-independent 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 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.

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 two-node beam element, a point on the left end-section can serve as the third point III, as shown in Figure 3(b)

(a)

(b)

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

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 , its skew symmetric matrix is expressed as

From (8), it can obtain the conclusion that is parallel to and has no component along , 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 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.

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 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).

(a)

(b)

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

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 . The current orientation matrix of the cross section of node in the global coordinate system can be described as

where 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 by

where 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 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]

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 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 end-section of the beam, so is not independent but related to and 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.

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 as

where 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 as

where , , , 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 as

where is the density of surface forces, 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, is the position of the mass center in the local coordinate system, and and 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 and are the submatrices of gravity influence coefficient matrix of the substructure. The remaining parts of can be rewritten as

where 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 by

where and are two coefficient matrices defined by

It 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 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).

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.

Figure 5: Three groups of nodes in the substructure .

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 block-partitioned form as follows:

After doing this, (44) can be rewritten in the form

In the above equation, is independent because the virtual powers of the remaining substructures have nothing to do with it, which yields

When more than six components in are constrained, 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

where , , , 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 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 can be determined by the virtual velocity of the local coordinate system origin and the virtual angular velocity of the left end-section represented by (33), and the force in (56) is carved up between and resulting in increments to them by

Meanwhile, 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 system

In 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 yields

where the derivatives of the local generalized internal nodal forces and torques with respect to the global parameters can be obtained by

According 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 by

Representing 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 , , 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 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.

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)

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

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 elements are shown, respectively, in Figure 7, where the meaning of legends is defined by Figure 6(b).

(a)

(b)

Figure 7: Load-displacement curves of the nodes on the beam when the beam is modeled by different numbers of elements.

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.

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

Euler buckling load
The number of elements along its axial direction

2
4
8
16

Critical load (kN)
139.953
147.316
141.791
140.434
140.096

Relative error (%)
—
5.26
1.31
0.34
0.10

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 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.

5.2. 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 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)

Figure 8: Geometry of 24-member star-shaped shallow dome.

By using elements for each member, the equilibrium paths of node obtained by the presented method are shown, respectively, in Figure 9.

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

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.

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

2
4
8
16

Critical load (N)
561.4
626.1
582.9
571.9
569.3

Relative error (%)
—
11.52
3.83
1.87
1.41

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.

Figure 10: Frame with a subsidiary beam.

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 will buckle first while the remaining parts remain stable, as illustrated in Figure 11.

Figure 11: Local buckling of the frame.

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)

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

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 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 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.

Figure 13: Load-displacement curves of and when all beams are modeled by twenty elements.

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³ kg/m³. The gravity acceleration is 9.807 m/s², 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)

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

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)

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

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).

Figure 16: Contour plots displacement vector sum of nodal solution of the slender frame structure.

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.

Figure 17: Load-displacement curves of the point .

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 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.

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.

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.

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. 16-18, pp. 2169–2185, 2003.

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. 108-109, pp. 110–117, 2012.

M. A. Crisfield, “A fast incremental/iterative solution procedure that handles ‘snap-through’,” Computers & Structures, vol. 13, no. 1–3, pp. 55–62, 1981.

M. A. Crisfield, “An arc-length method including line searches and accelerations,” International Journal for Numerical Methods in Engineering, vol. 19, no. 9, pp. 1269–1289, 1983.

M. Rezaiee-Pajand and J. Alamatian, “Automatic DR structural analysis of snap-through and snap-back using optimized load increments,” Journal of Structural Engineering, vol. 137, no. 1, pp. 109–116, 2011.

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.

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.

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.

M. A. Crisfield and G. F. Moita, “A unified co-rotational framework for solids, shells and beams,” International Journal of Solids and Structures, vol. 33, no. 20–22, pp. 2969–2992, 1996.

C. A. Felippa and B. Haugen, “A unified formulation of small-strain corotational finite elements: I. Theory,” Computer Methods in Applied Mechanics and Engineering, vol. 194, no. 21–24, pp. 2285–2335, 2005.

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 steady-state deformation and free vibration analysis of a rotating-inclined beam,” Mathematical Problems in Engineering, vol. 2011, Article ID 146505, 29 pages, 2011.

J. Rungamornrat, S. Sihanartkatakul, and P. Kanchanakitcharoen, “Adaptive, small-rotation-based, corotational technique for analysis of 2D nonlinear elastic frames,” Mathematical Problems in Engineering, vol. 2014, Article ID 640983, 15 pages, 2014.

J.-M. Battini and C. Pacoste, “Co-rotational beam elements with warping effects in instability problems,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 17-18, pp. 1755–1789, 2002.

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.

Y. He, X. Zhou, and P. Hou, “Combined method of super element and substructure for analysis of ILTDBS reticulated mega-structure with single-layer latticed shell substructures,” Finite Elements in Analysis and Design, vol. 46, no. 7, pp. 563–570, 2010.

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. 3-4, pp. 235–254, 1998.

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.

J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, 2008.

L. F. Shampine and M. W. Reichelt, “The MATLAB ODE suite,” SIAM Journal on Scientific Computing, vol. 18, no. 1, pp. 1–22, 1997.

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.

J. H. Argyris, “An excursion into large rotations,” Computer Methods in Applied Mechanics and Engineering, vol. 32, no. 1–3, pp. 85–155, 1982.