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

Authors: Zhaohui Qi; Xianchao Kong; Gang Wang
Journal: Mathematical Problems in Engineering
 (2015)

---

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


---

## Body

## 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 Section2, 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(1)δu,λ=f-λq0=0,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 q0 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):(2)δ˙=∂f∂ududλ-q0=0,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:(3)δ˙+αδ=0,where α is called a stabilization factor. In terms of the components of δ, its solution can be expressed as(4)δiδi0=δ˙iδ˙i0=e-αλ-λ0.It indicates that the violations in δ and δ˙ would be decreased exponentially when the displacement along the equilibrium path is solved from (3); that is,(5)∂f∂ududλ=q0-αf-λq0.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 Figure1(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.Figure 1
(a) Limit points on the equilibrium path; (b) bifurcation point on the equilibrium path.
(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 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 displacementdu/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 Figure2, without loss of generality, their displacement uin can be related to the existing displacement u by (6)ψuin,u=0.Differentiating (6) yields(7)∂ψ∂uinu˙in=-∂ψ∂uu˙.It indicates that if matrix ∂ψ/∂uin 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∂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.


## 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 Figure3(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)Figure 3
Local coordinate systems for (a) a substructure element and (b) a beam element.
(a)
(b)Its origin is placed at point I and its base vectors are defined by(8)e1=rII-rIrII-rI,e3=e~1rIII-rIe~1rIII-rI,e2=e~3e1,where a vector with a symbol “~” on the top is its skew symmetric matrix. For example, if a vector a=a1a2a3T, its skew symmetric matrix is expressed as(9)a~=0-a3a2a30-a1-a2a10.From (8), it can obtain the conclusion that rII-rI is parallel to e1 and rIII-rI has no component along e3, which yields(10)rII-rI=s0e1,rIII-rI=s1e1+s2e2.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(11)Acr=e1e2e3.By definition, the local coordinate system angular velocity ωcr can be expressed as(12)ωcr=e˙2Te3e1+e˙3Te1e2+e˙1Te2e3.It can be verified that the above equation can be expanded as(13)ωcr=Tcr2r˙II-r˙I+Tcr3r˙III-r˙I,where(14)Tcr2=1s0e3e2T-e2e3T-s1s0s2e1e3T;Tcr3=1s2e1e3T.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 nodei located at r-i will occupy a position in the global coordinate system, which yields(15)ri=rI+Acrr-i;ri0=rI0+Acr0r-i0,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(16)ui=ri-ri0;u-i=r-i-r-i0.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).Figure 4
Corotational kinematic descriptions for (a) a substructure element and (b) a beam element.
(a)
(b)In accordance with (15), the local displacement can be extracted from the global one as(17)u-i=AcrTui-uI+ri0-rI0-r-i0.Their rates of change are given by(18)u-˙i=AcrTu˙i-u˙I+r~i-r~Iωcr.If the node i 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 Ai. The current orientation matrix of the cross section of node i in the global coordinate system can be described as(19)Ai=RθiS0,where S0 is the global initial orientation matrix of the cross section of the node i, which is a constant; θi is the global rotational vector of node i in the sense that it is independent of the element. In terms of a rotational vector θ, we can express an orthogonal matrix as [23](20)Rθ=I+1-cos⁡θθ2θ~θ~+sin⁡θθθ~,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 θi is also large, and the expression of the corresponding rotational matrix cannot be simplified.Another way to describe the orientation matrixAi is given by(21)Ai=AcrRθ-iS-0,where S-0 is the local initial orientation matrix of the cross section of the node i, which is a constant; θ-i is the local rotational vector of node i, 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(22)Rθ-i=I+θ-~i.From (19) and (21), one can conclude that(23)Rθ-i=AcrTRθiAcr0,where θ-i and θi are the rotational vectors of node i 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:(24)ω~=R˙RT.If R represents the rotation of the coordinate system C1 relative to the coordinate system C0, ω defined as such is the angular velocity expressed in the coordinate system C0. In terms of the rotational vector θ and its rate of change θ˙, the angular velocity ω can be written as [23](25)ω=Tθθ˙,where(26)Tθ=I+1-sin⁡θθθ~θ~θ2+1-cos⁡θθ2θ~.Its inverse is given by(27)Tθ-1=I+1-θ/2tan⁡θ/2θ~θ~θ2-12θ~.When the norm of the rotational vector θ is small enough, the two matrices can be approximated as(28)Tθ≈I+12θ~;Tθ-1≈I-12θ~.According to (23)-(24), the relationship between the local angular velocity ω-i and global angular velocity ωi can be described as(29)ω-i=AcrTωi-ωcr.On the basis of (18), (29), and (13), relationships between the local and global virtual velocities can be given by(30)δu-˙i=AcrTδu˙i-δu˙I+r~i-r~Iδωcr,(31)δω-i=AcrTδωi-δωcr,where(32)δωcr=Tcr2δu˙II-δu˙I+Tcr3δu˙III-δu˙I.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:(33)δu˙III=δu˙I-r~III-r~IδωI.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 substructureSn can be expressed in the local coordinate system as(34)δpen=∑i∈Nnδu-˙iTf-i+δω-iTm-i,where Nn is the set of all nodes in this substructure. When the local displacement and rotations are small enough, the internal forces f-i and torques m-i can be obtained through some adequate linear elements as(35)f-i=∑j∈NnKijuuu-j+Kijuθθ-j,m-i=∑j∈NnKijθuu-j+Kijθθθ-j,where Kijuu, Kijuθ, Kijθu, and Kijθθ are the submatrices of stiffness matrix Ke 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δpsn and δpgn, respectively, are calculated in the global coordinate system and represented as(36)δpsn=∫sδu˙Ttsds;δpgn=∫vδu˙Tgρdv,where ts is the density of surface forces, g is the gravity acceleration, and ρ is the mass density. According to (30), they can be rewritten as(37)δpsn=δu˙ITAcr∫st-sds+δωcrTAcr∫sρ-~t-sds+∫sδu-˙Tt-sds,δpgn=mδu˙ITg+mδωcrTAcrr-~cg-+∫vδu-˙Tg-ρdv,where m is the total mass of the substructure, r-c is the position of the mass center in the local coordinate system, and t-s and g- are defined by(38)t-s=AcrTts;g-=AcrTg.By virtue of element shape functions, δpsn can be written in the following form:(39)δpsn=∑i∈Nnδu˙iTfis+δωiTmisand the integration in the expression of δpgn can be written as(40)δpginn=∫vδu-˙Tg-ρdv=∑i∈Nnδu-˙iTWiug-+δω-iTWiθg-,where Wiu and Wiθ are the submatrices of gravity influence coefficient matrix We of the substructure. The remaining parts of δpgn can be rewritten as(41)δpgcrn=mδu˙ITg+mδωcrTAcrr-~cg-=∑j∈Ncrnδu˙jTfjg,where Ncrn is the set of three nodes I, II, and III for defining the local coordinate system, and the generalized forces fjg resulted from gravities and are given by(42)fIg=mg-mAcrΓII+ΓIIIr-~cg-,fIIg=mAcrΓIIr-~cg-,fIIIg=mAcrΓIIIr-~cg-,where ΓII and ΓIII are two coefficient matrices defined by(43)ΓII=AcrTTcr2TAcr=1s0000001-s1s2-10;ΓIII=AcrTTcr3TAcr=1s2000000100.It can be seen that δpgcrn only affects the generalized forces of the nodes relevant to the local coordinate system, whereas δpginn 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(44)δpen-δpginn=∑i∈Nnδu-˙iTf-i-Wiug-+δω-iTm-i-Wiθg-.The principle of virtual power can be expressed as(45)∑Sn⊂Sδpen-δpginn=∑Sn⊂Sδpgcrn+δpsn,where the summation range includes all substructures in the system S. 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 substructureSn can be divided into three groups. The first group Ncrn comprises the three nodes I, II, and III for defining the local coordinate system. The second group Nsn 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 Nbn=Nsn∪Ncrn is defined as boundary nodes. The remaining nodes in the substructure are sorted in the third group and defined as interior nodes Ninn, as shown in Figure 5.Figure 5
Three groups of nodes in the substructureSn.Local displacement and rotational vectors of the nodes in the substructureSn can be regrouped into u-b or u-in according to the criterion of whether the node belongs to the set of boundary nodes Nbn or the set of interior nodes Ninn. Similarly, the stiffness matrix and gravity influence coefficient matrix of the substructure element can be written in block-partitioned form as follows:(46)Ke=KbbKbiKibKii;We=WbWin.After doing this, (44) can be rewritten in the form(47)δpen-δpginn=δu-˙bTKbbu-b+Kbiu-in-Wbg-+δu-˙inTKibu-b+Kiiu-in-Wing-.In the above equation, δu-˙in is independent because the virtual powers of the remaining substructures have nothing to do with it, which yields(48)Kibu-b+Kiiu-in-Wing-=0.When more than six components in u-b are constrained, Kii will be an invertible matrix. In this case, (47) can be simplified as(49)δpen-δpginn=δu-˙bTKsubu-b-Wsubg-,where the equivalent stiffness matrix(50)Ksub=Kbb-KbiKii-1Kiband the equivalent gravity influence coefficient matrix(51)Wsub=Wb-KbiKii-1Win.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(52)δpen-δpginn=∑k∈Nbnδu-˙kTf-k+δω-kTm-k,where k is an arbitrary node in boundary node group Nbn. The generalized forces and torques in the local coordinate system can be expressed as(53)f-k=∑j∈NbnK-ijuuu-j+K-ijuθθ-j-W-kug-,m-k=∑j∈NbnK-ijθuu-j+K-ijθθθ-j-W-kθg-,where K-kjuu, K-kjuθ, K-kjθu, and K-kjθθ and W-ku, W-kθ are the submatrices of Ksub and Wsub related to the translational and rotational displacement vectors.Substituting (30)–(32) into (41) and (52), we can obtain that(54)δpen-δpgn=∑k∈Nbnδu˙kTfk+δωkTmk.In this equation, the generalized forces and torques in the global coordinate system are given by(55)fk=Acrf-k,k∈Nsnmk=Acrm-k,k∈Nbnfor the nodes in the set Ncrn(56)fI=Acrf-I+ΓII+ΓIIIm-Σ-f-ΣfII=Acrf-II-ΓIIm-ΣfIII=Acrf-III-ΓIIIm-Σ,where(57)f-Σ=∑k∈Nbnf-k+mg-m-Σ=∑k∈Nbnr-~kf-k+m-k+mr-~cg-.In the case that the substructureSn 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 fIII in (56) is carved up between fI and mI resulting in increments to them by(58)ΔfI=fIIIΔmI=r~III-r~IfIII.Meanwhile, fI and mI should be updated to be fI+ΔfI and mI+ΔmI, 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(59)∑Sn⊂S∑k∈Nbnδu˙kTfk-fks+δωkTmk-mks=0.In this equation, the generalized internal forces fk and torques mk 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.


## 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 Figure3(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)Figure 3
Local coordinate systems for (a) a substructure element and (b) a beam element.
(a)
(b)Its origin is placed at point I and its base vectors are defined by(8)e1=rII-rIrII-rI,e3=e~1rIII-rIe~1rIII-rI,e2=e~3e1,where a vector with a symbol “~” on the top is its skew symmetric matrix. For example, if a vector a=a1a2a3T, its skew symmetric matrix is expressed as(9)a~=0-a3a2a30-a1-a2a10.From (8), it can obtain the conclusion that rII-rI is parallel to e1 and rIII-rI has no component along e3, which yields(10)rII-rI=s0e1,rIII-rI=s1e1+s2e2.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(11)Acr=e1e2e3.By definition, the local coordinate system angular velocity ωcr can be expressed as(12)ωcr=e˙2Te3e1+e˙3Te1e2+e˙1Te2e3.It can be verified that the above equation can be expanded as(13)ωcr=Tcr2r˙II-r˙I+Tcr3r˙III-r˙I,where(14)Tcr2=1s0e3e2T-e2e3T-s1s0s2e1e3T;Tcr3=1s2e1e3T.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 nodei located at r-i will occupy a position in the global coordinate system, which yields(15)ri=rI+Acrr-i;ri0=rI0+Acr0r-i0,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(16)ui=ri-ri0;u-i=r-i-r-i0.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).Figure 4
Corotational kinematic descriptions for (a) a substructure element and (b) a beam element.
(a)
(b)In accordance with (15), the local displacement can be extracted from the global one as(17)u-i=AcrTui-uI+ri0-rI0-r-i0.Their rates of change are given by(18)u-˙i=AcrTu˙i-u˙I+r~i-r~Iωcr.If the node i 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 Ai. The current orientation matrix of the cross section of node i in the global coordinate system can be described as(19)Ai=RθiS0,where S0 is the global initial orientation matrix of the cross section of the node i, which is a constant; θi is the global rotational vector of node i in the sense that it is independent of the element. In terms of a rotational vector θ, we can express an orthogonal matrix as [23](20)Rθ=I+1-cos⁡θθ2θ~θ~+sin⁡θθθ~,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 θi is also large, and the expression of the corresponding rotational matrix cannot be simplified.Another way to describe the orientation matrixAi is given by(21)Ai=AcrRθ-iS-0,where S-0 is the local initial orientation matrix of the cross section of the node i, which is a constant; θ-i is the local rotational vector of node i, 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(22)Rθ-i=I+θ-~i.From (19) and (21), one can conclude that(23)Rθ-i=AcrTRθiAcr0,where θ-i and θi are the rotational vectors of node i 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:(24)ω~=R˙RT.If R represents the rotation of the coordinate system C1 relative to the coordinate system C0, ω defined as such is the angular velocity expressed in the coordinate system C0. In terms of the rotational vector θ and its rate of change θ˙, the angular velocity ω can be written as [23](25)ω=Tθθ˙,where(26)Tθ=I+1-sin⁡θθθ~θ~θ2+1-cos⁡θθ2θ~.Its inverse is given by(27)Tθ-1=I+1-θ/2tan⁡θ/2θ~θ~θ2-12θ~.When the norm of the rotational vector θ is small enough, the two matrices can be approximated as(28)Tθ≈I+12θ~;Tθ-1≈I-12θ~.According to (23)-(24), the relationship between the local angular velocity ω-i and global angular velocity ωi can be described as(29)ω-i=AcrTωi-ωcr.On the basis of (18), (29), and (13), relationships between the local and global virtual velocities can be given by(30)δu-˙i=AcrTδu˙i-δu˙I+r~i-r~Iδωcr,(31)δω-i=AcrTδωi-δωcr,where(32)δωcr=Tcr2δu˙II-δu˙I+Tcr3δu˙III-δu˙I.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:(33)δu˙III=δu˙I-r~III-r~IδωI.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 substructureSn can be expressed in the local coordinate system as(34)δpen=∑i∈Nnδu-˙iTf-i+δω-iTm-i,where Nn is the set of all nodes in this substructure. When the local displacement and rotations are small enough, the internal forces f-i and torques m-i can be obtained through some adequate linear elements as(35)f-i=∑j∈NnKijuuu-j+Kijuθθ-j,m-i=∑j∈NnKijθuu-j+Kijθθθ-j,where Kijuu, Kijuθ, Kijθu, and Kijθθ are the submatrices of stiffness matrix Ke 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δpsn and δpgn, respectively, are calculated in the global coordinate system and represented as(36)δpsn=∫sδu˙Ttsds;δpgn=∫vδu˙Tgρdv,where ts is the density of surface forces, g is the gravity acceleration, and ρ is the mass density. According to (30), they can be rewritten as(37)δpsn=δu˙ITAcr∫st-sds+δωcrTAcr∫sρ-~t-sds+∫sδu-˙Tt-sds,δpgn=mδu˙ITg+mδωcrTAcrr-~cg-+∫vδu-˙Tg-ρdv,where m is the total mass of the substructure, r-c is the position of the mass center in the local coordinate system, and t-s and g- are defined by(38)t-s=AcrTts;g-=AcrTg.By virtue of element shape functions, δpsn can be written in the following form:(39)δpsn=∑i∈Nnδu˙iTfis+δωiTmisand the integration in the expression of δpgn can be written as(40)δpginn=∫vδu-˙Tg-ρdv=∑i∈Nnδu-˙iTWiug-+δω-iTWiθg-,where Wiu and Wiθ are the submatrices of gravity influence coefficient matrix We of the substructure. The remaining parts of δpgn can be rewritten as(41)δpgcrn=mδu˙ITg+mδωcrTAcrr-~cg-=∑j∈Ncrnδu˙jTfjg,where Ncrn is the set of three nodes I, II, and III for defining the local coordinate system, and the generalized forces fjg resulted from gravities and are given by(42)fIg=mg-mAcrΓII+ΓIIIr-~cg-,fIIg=mAcrΓIIr-~cg-,fIIIg=mAcrΓIIIr-~cg-,where ΓII and ΓIII are two coefficient matrices defined by(43)ΓII=AcrTTcr2TAcr=1s0000001-s1s2-10;ΓIII=AcrTTcr3TAcr=1s2000000100.It can be seen that δpgcrn only affects the generalized forces of the nodes relevant to the local coordinate system, whereas δpginn 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(44)δpen-δpginn=∑i∈Nnδu-˙iTf-i-Wiug-+δω-iTm-i-Wiθg-.The principle of virtual power can be expressed as(45)∑Sn⊂Sδpen-δpginn=∑Sn⊂Sδpgcrn+δpsn,where the summation range includes all substructures in the system S. 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 substructureSn can be divided into three groups. The first group Ncrn comprises the three nodes I, II, and III for defining the local coordinate system. The second group Nsn 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 Nbn=Nsn∪Ncrn is defined as boundary nodes. The remaining nodes in the substructure are sorted in the third group and defined as interior nodes Ninn, as shown in Figure 5.Figure 5
Three groups of nodes in the substructureSn.Local displacement and rotational vectors of the nodes in the substructureSn can be regrouped into u-b or u-in according to the criterion of whether the node belongs to the set of boundary nodes Nbn or the set of interior nodes Ninn. Similarly, the stiffness matrix and gravity influence coefficient matrix of the substructure element can be written in block-partitioned form as follows:(46)Ke=KbbKbiKibKii;We=WbWin.After doing this, (44) can be rewritten in the form(47)δpen-δpginn=δu-˙bTKbbu-b+Kbiu-in-Wbg-+δu-˙inTKibu-b+Kiiu-in-Wing-.In the above equation, δu-˙in is independent because the virtual powers of the remaining substructures have nothing to do with it, which yields(48)Kibu-b+Kiiu-in-Wing-=0.When more than six components in u-b are constrained, Kii will be an invertible matrix. In this case, (47) can be simplified as(49)δpen-δpginn=δu-˙bTKsubu-b-Wsubg-,where the equivalent stiffness matrix(50)Ksub=Kbb-KbiKii-1Kiband the equivalent gravity influence coefficient matrix(51)Wsub=Wb-KbiKii-1Win.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(52)δpen-δpginn=∑k∈Nbnδu-˙kTf-k+δω-kTm-k,where k is an arbitrary node in boundary node group Nbn. The generalized forces and torques in the local coordinate system can be expressed as(53)f-k=∑j∈NbnK-ijuuu-j+K-ijuθθ-j-W-kug-,m-k=∑j∈NbnK-ijθuu-j+K-ijθθθ-j-W-kθg-,where K-kjuu, K-kjuθ, K-kjθu, and K-kjθθ and W-ku, W-kθ are the submatrices of Ksub and Wsub related to the translational and rotational displacement vectors.Substituting (30)–(32) into (41) and (52), we can obtain that(54)δpen-δpgn=∑k∈Nbnδu˙kTfk+δωkTmk.In this equation, the generalized forces and torques in the global coordinate system are given by(55)fk=Acrf-k,k∈Nsnmk=Acrm-k,k∈Nbnfor the nodes in the set Ncrn(56)fI=Acrf-I+ΓII+ΓIIIm-Σ-f-ΣfII=Acrf-II-ΓIIm-ΣfIII=Acrf-III-ΓIIIm-Σ,where(57)f-Σ=∑k∈Nbnf-k+mg-m-Σ=∑k∈Nbnr-~kf-k+m-k+mr-~cg-.In the case that the substructureSn 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 fIII in (56) is carved up between fI and mI resulting in increments to them by(58)ΔfI=fIIIΔmI=r~III-r~IfIII.Meanwhile, fI and mI should be updated to be fI+ΔfI and mI+ΔmI, 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(59)∑Sn⊂S∑k∈Nbnδu˙kTfk-fks+δωkTmk-mks=0.In this equation, the generalized internal forces fk and torques mk 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(60)∂fk∂uI=Acr∂f-k∂uI+f~kTcr2+Tcr3∂fk∂uII=Acr∂f-k∂uII-f~kTcr2∂fk∂uIII=Acr∂f-k∂uIII-f~kTcr3∂fk∂uj=Acr∂f-k∂uj,j∈Nsn∂fk∂θj=Acr∂f-k∂θj,j∈Nbn,k∈Nsn∂mk∂uI=Acr∂m-k∂uI+m~kTcr2+Tcr3∂mk∂uII=Acr∂m-k∂uII-m~kTcr2∂mk∂uIII=Acr∂m-k∂uIII-m~kTcr3∂mk∂uj=Acr∂m-k∂uj,j∈Nsn∂mk∂θj=Acr∂m-k∂θj,j∈Nbn,k∈Nbn,where the derivatives of the local generalized internal nodal forces and torques with respect to the global parameters can be obtained by(61)∂f-k∂uI=K-k,IuuAcrT-∂f-k∂ωcrTcr2+Tcr3-∑j∈NbnK-kjuuAcrT∂f-k∂uII=K-k,IIuuAcrT+∂f-k∂ωcrTcr2∂f-k∂uIII=K-k,IIIuuAcrT+∂f-k∂ωcrTcr3∂f-k∂ωcr=∑j∈NbnK-kjuuAcrTr~j-r~I-K-kjuθTθ-j-1AcrT-W-kuAcrTg~∂f-k∂uj=K-kjuuAcrT,j∈Nsn∂f-k∂θj=K-kjuθTθ-j-1AcrTTθj,j∈Nbn,k∈Nsn∂m-k∂uI=K-k,IθuAcrT-∂m-k∂ωcrTcr2+Tcr3-∑j∈NbnK-kjθuAcrT∂m-k∂uII=K-k,IIθuAcrT+∂m-k∂ωcrTcr2∂m-k∂uIII=K-k,IIIθuAcrT+∂m-k∂ωcrTcr3∂m-k∂ωcr=∑j∈NbnK-kjθuAcrTr~j-r~I-K-kjθθTθ-j-1AcrT-W-kθAcrTg~∂m-k∂uj=K-kjθuAcrT,j∈Nsn∂m-k∂θj=K-kjθθTθ-j-1AcrTTθj,j∈Nbn,k∈Nbn.According to (43), it can be verified that, for the vector expressed in the local coordinate system m-Σ=m-Σ,1e1+m-Σ,2e2+m-Σ,3e3, the following equations are valid:(62)
Γ
˙
I
I
m
-
Σ
=
1
s
0
2
s
2
s
2
m
-
Σ
,
1
+
s
1
m
-
Σ
,
2
s
˙
0
-
s
0
m
-
Σ
,
1
s
˙
2
-
s
0
m
-
Σ
,
2
s
˙
I
Γ
˙
I
I
I
m
-
Σ
=
-
1
s
2
2
m
-
Σ
,
2
s
˙
2
,where m-Σ,1=0-m-Σ,3m-Σ,2T and m-Σ,2=00m-Σ,1T.By the definitions ofs0, s1, and s2 as shown by (10), their derivatives are given by(63)s˙0=e1Tu˙II-u˙Is˙1=e1Tu˙III-u˙I+s2e3Tωcrs˙2=e2Tu˙III-u˙I-s1e3Tωcr.Representing ωcr in the above equation by (13) and substituting the resulting equation into (62) yield that(64)Γ˙IIm-Σ=BΣ,IIu˙II-u˙I+BΣ,IIIu˙III-u˙IΓ˙IIIm-Σ=CΣ,IIu˙II-u˙I+CΣ,IIIu˙III-u˙I,where the coefficients matrices are defined by(65)BΣ,II=1s02s2s2m-Σ,1+s1m-Σ,2e1T+1s02s2s1m-Σ,1-s2m-Σ,2e2T;BΣ,III=-1s0s2m-Σ,1e2T+m-Σ,2e1T;CΣ,II=1s0s22s1m-Σ,2e2T;CΣ,III=-1s0s22s0m-Σ,2e2T.In accordance with (56), the derivatives of the global generalized internal and gravity forces fI, fII, and fIII with respect to the global parameters are obtained by (66)
∂
f
I
∂
u
1
=
A
c
r
∂
f
-
I
∂
u
I
+
Γ
I
I
+
Γ
I
I
I
∂
m
-
Σ
∂
u
I
-
B
Σ
,
I
I
+
B
Σ
,
I
I
I
+
C
Σ
,
I
I
+
C
Σ
,
I
I
I
-
∂
f
-
Σ
∂
u
I
+
f
~
I
T
c
r
2
+
T
c
r
3
∂
f
I
∂
u
I
I
=
A
c
r
∂
f
-
I
∂
u
I
I
+
Γ
I
I
+
Γ
I
I
I
∂
m
-
Σ
∂
u
I
I
+
B
Σ
,
I
I
+
C
Σ
,
I
I
-
∂
f
-
Σ
∂
u
I
I
-
f
~
I
T
c
r
2
∂
f
I
∂
u
I
I
I
=
A
c
r
∂
f
-
I
∂
u
I
I
I
+
Γ
I
I
+
Γ
I
I
I
∂
m
-
Σ
∂
u
I
I
I
+
B
Σ
,
I
I
I
+
C
Σ
,
I
I
I
-
∂
f
-
Σ
∂
u
I
I
I
-
f
~
I
T
c
r
3
∂
f
I
∂
u
j
=
A
c
r
∂
f
-
I
∂
u
j
+
Γ
I
I
+
Γ
I
I
I
∂
m
-
Σ
∂
u
j
-
∂
f
-
Σ
∂
u
j
,
j
∈
N
s
n
∂
f
I
∂
θ
j
=
A
c
r
∂
f
-
I
∂
θ
j
+
Γ
I
I
+
Γ
I
I
I
∂
m
-
Σ
∂
θ
j
-
∂
f
-
Σ
∂
θ
j
,
j
∈
N
b
n
∂
f
I
I
∂
u
I
=
A
c
r
∂
f
-
I
I
∂
u
I
-
Γ
I
I
∂
m
-
Σ
∂
u
I
+
B
Σ
,
I
I
+
B
Σ
,
I
I
I
+
f
~
I
I
T
c
r
2
+
T
c
r
3
∂
f
I
I
∂
u
I
I
=
A
c
r
∂
f
-
I
I
∂
u
I
I
-
Γ
I
I
∂
m
-
Σ
∂
u
I
I
+
B
Σ
,
I
I
-
f
~
I
I
T
c
r
2
∂
f
I
I
∂
u
I
I
I
=
A
c
r
∂
f
-
I
I
∂
u
I
I
I
-
Γ
I
I
∂
m
-
Σ
∂
u
I
I
I
+
B
Σ
,
I
I
I
-
f
~
I
I
T
c
r
3
∂
f
I
I
∂
u
j
=
A
c
r
∂
f
-
I
I
∂
u
j
-
Γ
I
I
∂
m
-
Σ
∂
u
j
,
j
∈
N
s
n
∂
f
I
I
∂
θ
j
=
A
c
r
∂
f
-
I
I
∂
θ
j
-
Γ
I
I
∂
m
-
Σ
∂
θ
j
,
j
∈
N
b
n
∂
f
I
I
I
∂
u
I
=
A
c
r
∂
f
-
I
I
I
∂
u
I
-
Γ
I
I
I
∂
m
-
Σ
∂
u
I
+
C
Σ
,
I
I
+
C
Σ
,
I
I
I
+
f
~
I
I
I
T
c
r
2
+
T
c
r
3
∂
f
I
I
I
∂
u
I
I
=
A
c
r
∂
f
-
I
I
I
∂
u
I
I
-
Γ
I
I
I
∂
m
-
Σ
∂
u
I
I
+
C
Σ
,
I
I
-
f
~
I
I
I
T
c
r
2
∂
f
I
I
I
∂
u
I
I
I
=
A
c
r
∂
f
-
I
I
I
∂
u
I
I
I
-
Γ
I
I
I
∂
m
-
Σ
∂
u
I
I
I
+
C
Σ
,
I
I
I
-
f
~
I
I
I
T
c
r
3
∂
f
I
I
I
∂
u
j
=
A
c
r
∂
f
-
I
I
I
∂
u
j
-
Γ
I
I
I
∂
m
-
Σ
∂
u
j
,
j
∈
N
s
n
∂
f
I
I
I
∂
θ
j
=
A
c
r
∂
f
-
I
I
I
∂
θ
j
-
Γ
I
I
I
∂
m
-
Σ
∂
θ
j
,
j
∈
N
b
n
.According to (57), the derivatives of f-Σ and m-Σ with respect to the global parameters in (66) can be given by(67)
∂
f
-
Σ
∂
u
I
=
-
m
A
c
r
T
g
~
T
c
r
2
+
T
c
r
3
+
∑
k
∈
N
b
n
∂
f
-
k
∂
u
I
∂
f
-
Σ
∂
u
I
I
=
m
A
c
r
T
g
~
T
c
r
2
+
∑
k
∈
N
b
n
∂
f
-
k
∂
u
I
I
∂
f
-
Σ
∂
u
I
I
I
=
m
A
c
r
T
g
~
T
c
r
3
+
∑
k
∈
N
b
n
∂
f
-
k
∂
u
I
I
I
∂
f
-
Σ
∂
u
j
=
∑
k
∈
N
b
n
∂
f
-
k
∂
u
j
,
j
∈
N
s
n
∂
f
-
Σ
∂
θ
j
=
∑
k
∈
N
b
n
∂
f
-
k
∂
θ
j
,
j
∈
N
b
n
∂
m
-
Σ
∂
u
I
=
-
m
r
-
~
c
A
c
r
T
g
~
T
c
r
2
+
T
c
r
3
-
f
-
~
I
A
c
r
T
+
∑
k
∈
N
b
n
r
-
~
k
∂
f
-
k
∂
u
I
+
∂
m
-
k
∂
u
I
+
f
-
~
k
A
c
r
T
+
f
-
~
k
A
c
r
T
r
~
k
-
r
~
I
T
c
r
2
+
T
c
r
3
∂
m
-
Σ
∂
u
I
I
=
m
r
-
~
c
A
c
r
T
g
~
T
c
r
2
-
f
-
~
I
I
A
c
r
T
+
∑
k
∈
N
b
n
r
-
~
k
∂
f
-
k
∂
u
I
I
+
∂
m
-
k
∂
u
I
I
-
f
-
~
k
A
c
r
T
r
~
k
-
r
~
I
T
c
r
2
∂
m
-
Σ
∂
u
I
I
I
=
m
r
-
~
c
A
c
r
T
g
~
T
c
r
3
-
f
-
~
I
I
I
A
c
r
T
+
∑
k
∈
N
b
n
r
-
~
k
∂
f
-
k
∂
u
I
I
I
+
∂
m
-
k
∂
u
I
I
I
-
f
-
~
k
A
c
r
T
r
~
k
-
r
~
I
T
c
r
3
∂
m
-
Σ
∂
u
j
=
-
f
-
~
j
A
c
r
T
+
∑
k
∈
N
b
n
r
-
~
k
∂
f
-
k
∂
u
j
+
∂
m
-
k
∂
u
j
,
j
∈
N
s
n
∂
m
-
Σ
∂
θ
j
=
∑
k
∈
N
b
n
r
-
~
k
∂
f
-
k
∂
θ
j
+
∂
m
-
k
∂
θ
j
,
j
∈
N
b
n
.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 Section5.2, all the members in the remaining three examples are made of the material with elasticity modulus E=2.1×1011Pa and Poisson’s ratio ν=0.3.


### 5.1. A Slender Axial Compression Beam

The first example is a slender beam clamped at its bottom and subjected to an axial forceF at its top, as shown in Figure 6(a). The length of the beam is L=10 m, its area is A=6×10-3m2, its polar moment of inertia is J=5.4020×10-5m4, and its two moments of inertia are Ix=Iy=2.7010×10-5m4.Figure 6
(a) Axial compression beam and (b) its model withn elements.
(a)
(b)According to Euler’s formula of fixed-free beam under axial compression, the critical load of this simple structure is given by(68)Fcr=π2EI4L2=139.953kN.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 n∈2,4,8,16 elements are shown, respectively, in Figure 7, where the meaning of legends is defined by Figure 6(b).Figure 7
Load-displacement curves of the nodes on the beam when the beam is modeled by different numbers of elements.
(a)
(b)The critical load can be detected by the maximum norm of the displacement derivatives with respect to axial forcedu/dF. 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.10It 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 erroris 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δ=0.0012 m corresponding to the obtained critical load, so the more precise critical load is given by(69)Fcr=π2EI4L-δ2=139.987kN.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 E=3.030×109Pa and shear modulus G=1.096×109Pa. The geometric parameters are A=3.170×10-4m2, Iy=2.377×10-8m4, Iz=2.95×10-9m4, J=9.18×10-9m4, and mass density ρ=5×104N/m3. The structure is pinned and restrained against translational motion. The load case is shown in Figure 8, where the node k is under a vertical load F along the negative y-axis direction.Figure 8
Geometry of 24-member star-shaped shallow dome.
(a)
(b)By usingn∈2,4,8,16 elements for each member, the equilibrium paths of node k obtained by the presented method are shown, respectively, in Figure 9.Figure 9
Load-displacement curves of the nodek by using n elements for each member.In this example, the derivativeduk/dF 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.41This 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 Figure10, a frame is composed of four main beams and a subsidiary beam b0, carrying a vertical load F 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 A=7.0686×10-4m2, polar moment of inertia J=3.6226×10-7m4, and two moments of inertia Ix=Iy=1.8113×10-7m4, 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 Figure10. 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 beamb0 will buckle first while the remaining parts remain stable, as illustrated in Figure 11.Figure 11
Local buckling of the frame.When the subsidiary beamb0 is modeled by only one element and the remaining beams are divided into twenty elements, the equilibrium path of u1 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.Figure 12
Load-displacement curve ofu1 when b0 is modeled by one element.
(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 beamb0. However, the first buckling load of the beam b0 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 ofu1 and u2 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 du2/dF is taken as the identification parameter. However, the local buckling of b0 does not lead to collapse of the frame. When we take du1/dF as the identification parameter, the critical load of the frame will reach 1223.360 kN.Figure 13
Load-displacement curves ofu1 and u2 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 loadF 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 × 103 kg/m3. The gravity acceleration is 9.807 m/s2, which is parallel to the negative y-axis. A small perturbation moment of 10 N·m is applied on the two nodes a and b (Figure 14(a)) of its top section to make the detection of the critical load easier.Figure 14
(a) Slender space frame structure carrying an axial load; (b) geometric parameters of the periodic insert section.
(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 Figure15(a). Consequently, the whole frame structure is divided into ten same corotational substructure elements, as shown in Figure 15(b).Figure 15
(a) Corotational substructure element of an insert section; (b) substructure model of the slender frame structure.
(a)
(b)For all boundary nodes displacement values, the absolute values of all components of the vectordu/dF, 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 Figure16 (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 pointc (Figure 16) obtained in this paper are compared with the ANSYS results, as shown in Figure 17.Figure 17
Load-displacement curves of the pointc.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.


## 5.1. A Slender Axial Compression Beam

The first example is a slender beam clamped at its bottom and subjected to an axial forceF at its top, as shown in Figure 6(a). The length of the beam is L=10 m, its area is A=6×10-3m2, its polar moment of inertia is J=5.4020×10-5m4, and its two moments of inertia are Ix=Iy=2.7010×10-5m4.Figure 6
(a) Axial compression beam and (b) its model withn elements.
(a)
(b)According to Euler’s formula of fixed-free beam under axial compression, the critical load of this simple structure is given by(68)Fcr=π2EI4L2=139.953kN.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 n∈2,4,8,16 elements are shown, respectively, in Figure 7, where the meaning of legends is defined by Figure 6(b).Figure 7
Load-displacement curves of the nodes on the beam when the beam is modeled by different numbers of elements.
(a)
(b)The critical load can be detected by the maximum norm of the displacement derivatives with respect to axial forcedu/dF. 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.10It 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 erroris 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δ=0.0012 m corresponding to the obtained critical load, so the more precise critical load is given by(69)Fcr=π2EI4L-δ2=139.987kN.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 E=3.030×109Pa and shear modulus G=1.096×109Pa. The geometric parameters are A=3.170×10-4m2, Iy=2.377×10-8m4, Iz=2.95×10-9m4, J=9.18×10-9m4, and mass density ρ=5×104N/m3. The structure is pinned and restrained against translational motion. The load case is shown in Figure 8, where the node k is under a vertical load F along the negative y-axis direction.Figure 8
Geometry of 24-member star-shaped shallow dome.
(a)
(b)By usingn∈2,4,8,16 elements for each member, the equilibrium paths of node k obtained by the presented method are shown, respectively, in Figure 9.Figure 9
Load-displacement curves of the nodek by using n elements for each member.In this example, the derivativeduk/dF 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.41This 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 Figure10, a frame is composed of four main beams and a subsidiary beam b0, carrying a vertical load F 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 A=7.0686×10-4m2, polar moment of inertia J=3.6226×10-7m4, and two moments of inertia Ix=Iy=1.8113×10-7m4, 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 Figure10. 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 beamb0 will buckle first while the remaining parts remain stable, as illustrated in Figure 11.Figure 11
Local buckling of the frame.When the subsidiary beamb0 is modeled by only one element and the remaining beams are divided into twenty elements, the equilibrium path of u1 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.Figure 12
Load-displacement curve ofu1 when b0 is modeled by one element.
(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 beamb0. However, the first buckling load of the beam b0 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 ofu1 and u2 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 du2/dF is taken as the identification parameter. However, the local buckling of b0 does not lead to collapse of the frame. When we take du1/dF as the identification parameter, the critical load of the frame will reach 1223.360 kN.Figure 13
Load-displacement curves ofu1 and u2 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 loadF 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 × 103 kg/m3. The gravity acceleration is 9.807 m/s2, which is parallel to the negative y-axis. A small perturbation moment of 10 N·m is applied on the two nodes a and b (Figure 14(a)) of its top section to make the detection of the critical load easier.Figure 14
(a) Slender space frame structure carrying an axial load; (b) geometric parameters of the periodic insert section.
(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 Figure15(a). Consequently, the whole frame structure is divided into ten same corotational substructure elements, as shown in Figure 15(b).Figure 15
(a) Corotational substructure element of an insert section; (b) substructure model of the slender frame structure.
(a)
(b)For all boundary nodes displacement values, the absolute values of all components of the vectordu/dF, 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 Figure16 (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 pointc (Figure 16) obtained in this paper are compared with the ANSYS results, as shown in Figure 17.Figure 17
Load-displacement curves of the pointc.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.



---

*Source: 325969-2015-08-23.xml*