Consistent co-rotational framework for Euler-Bernoulli and Timoshenko beam-column elements under distributed member loads

Abstract

Conventional co-rotational formulations for geometrically nonlinear analysis are based on the assumption that the finite element is only subjected to nodal loads and as a result, they are not accurate for the elements under distributed member loads. The magnitude and direction of member loads are treated as constant in the global coordinate system, but they are essentially varying in the local coordinate system for the element undergoing a large rigid body rotation, leading to the change of nodal moments at element ends. Thus, there is a need to improve the co-rotational formulations to allow for the effect. This paper proposes a new consistent co-rotational formulation for both Euler-Bernoulli and Timoshenko two-dimensional beam-column elements subjected to distributed member loads. It is found that the equivalent nodal moments are affected by the element geometric change and consequently contribute to a part of geometric stiffness matrix. From this study, the results of both eigenvalue buckling and second-order direct analyses will be significantly improved. Several examples are used to verify the proposed formulation with comparison of the traditional method, which demonstrate the accuracy and reliability of the proposed method in buckling analysis of frame structures under distributed member loads using a single element per member.

Keywords

beam-column element buckling analysis co-rotational formulation distributed member load

Introduction

The total Lagrangian (TL) and updated Lagrangian (UL) formulations based on the Green-Lagrangian strains have been widely used in finite element analysis of geometrically nonlinear structural and solid mechanics problems, which can be found in many textbooks such as Zienkiewicz (2005) and Bathe (2006). The co-rotational (CR) formulation is another well-accepted method (Felippa and Haugen 2005) in nonlinear analysis, especially in frame analysis due to significant simplification of element formulations.

The early geometrically nonlinear analysis procedures were mainly based on TL and UL formulations and directly used the Green-Lagrangian strains to consider large displacements and large rotations, since the nonlinear strains can remove rigid body motions. For example, Bathe and his coworkers (Bathe and Bolourchi 1980; Bathe et al. 1983; Dvorkin and Bathe 1984) proposed geometrically degenerated shell elements using the Green-Lagrangian strains. However, these works are limited to nonlinear small rotations between two load increments, as reported by Surana (1983). He proposed an expression of three-dimensional (3D) large rotations to improve the previous works, but his element formulation should determine the rotation sequence firstly, causing inconvenient in practical applications.

How to deal with 3D large rotations in geometrically nonlinear finite element method attracted considerable research interest such that numerous different expressions have been proposed. One of the most accepted methods was introduced by Argyris (1982). His work is based on Euler theorem and still often used in the studies of geometrically nonlinear elements, for example, the work by Jeon (2015).

All these mentioned geometrically nonlinear elements are mainly derived using total displacement functions and the Green-Lagrangian strains, while the pertaining analysis algorithms are mainly conducted in the fixed global coordinate system. Thus, their displacement functions and element formulations are very complicated, especially when an expression of 3D large rotation is used.

The co-rotational concept can be integrated with the TL and UL formulations, in which the total motions of an element are divided into the rigid body part and the pure deformational part by attaching a local coordinate system to the element. Only the pure deformational motions, which can be described by the local element formulation, produce strain energy. In this analysis procedure, formulations are established in the local element system firstly and then transformed into the global coordinate system. Thus, the co-rotational TL and UL formulations are simpler than the general formulations mentioned before. For example, Hsiao (1987) and Hsiao and Hung (1989) proposed geometrically nonlinear shell elements based on co-rotational TL and UL formulations, in which a co-rotational procedure is used to eliminate the rigid body motions and get the element internal nodal forces. Yang and Chiou (1987) proposed a rigid body motion test to verify whether a finite element can be used in a geometrically nonlinear analysis, and this test is consistent with the co-rotational concept. The other applications of the co-rotational TL and UL formulations can be found in the works by Tang et al. (2016), Jiang et al. (1994), and Khosravi et al. (2007, 2008).

The co-rotational (CR) formulation discussed in this paper is also called the element independent co-rotational (EICR) formulation which was firstly introduced by Rankin and Brogan (1986). This kind of formulations is derived through establishing the relationship between local and global coordinate systems. Thus, this CR formulation has an explicit geometric stiffness matrix expressed with end nodal forces rather than stresses at the integration points for TL and UL formulations. This feature leads to a unified geometric stiffness matrix for the elements having the same nodes and degrees of freedom. For this reason, this co-rotational formulation is element-independent. Moreover, the material nonlinearity can be simply considered in the local element formulation with significant enhancement of numerical efficiency.

A strict EICR formulation for 3D finite elements is also very complicated, in which the changes in element dimensions, the rigid body rotations of the element and the non-commutative nature of 3D large rotations should be considered in the derivation (Battini 2002; Crisfield 1990; Nour-Omid and Rankin 1991; Skallerud and Haugen 1999). In addition, the tangent stiffness and the internal forces are strictly consistent since the tangent stiffness is derived through the partial differentiation of the internal forces. To simplify EICR formulations, some simplified methods were proposed and still have good performance with high accuracy and efficiency. For example, Tang et al. (2017) neglected the changes in element dimensions and significantly simplified the EICR formulation for shell elements; Izzuddin (2005) proposed a co-rotational method for shell elements without removing rigid body motion completely. This method was further improved by Izzuddin and Liang (2016) and Li et al. (2017, 2018).

Before the introduction of EICR formulation by Rankin and Brogan (1986), a simplified EICR formulation had been used in the geometrically nonlinear analysis of frames, such as in the work by Meek and Tan (1983). Their CR formulation has been extensively used in the derivation of nonlinear beam-column elements (Chan and Zhou 1994, 1995; Du et al. 2017; Tang et al. 2015, 2019). Compared with the conventional TL and UL formulations, the most prominent advantage of the EICR formulation is that some complicated beam-column elements allowing for local second-order effect can be readily extended into the analyses of large displacements and large rotations.

It should be pointed out that the previous EICR formulations assumed that the applied loads are nodal. This assumption is valid for the geometrically nonlinear analysis of structures subjected to constant distributed loads without significant geometric change. However, it is not reasonable for the case when the element is undergoing large rotations. To solve this problem, the traditional methods should use several elements per member. It is no doubt that the computational time will be significantly increased. It will also bring problem in modeling of member geometrical imperfection using one element per member. Thus, the changes of the equivalent nodal moments should be considered in the exact CR formulation. Note that the distributed member load is not only due to the self-weight of structural element but also those from non-structural elements such as claddings, liquid in pipe, ice loading, and so on. It is worth to propose a method to improve the accuracy of EICR formulations for nonlinear analysis.

In this paper, a new co-rotational framework is proposed, which is suitable for both Euler-Bernoulli and Timoshenko beam-column elements. The 2D problem is studied here for clear illustration of the proposed method. First, the CR formulation for 2D beam-column elements only subjected to nodal forces are introduced and its physical meanings are discussed. After that, the exact CR formulation for beam-column elements under distributed member loads is derived. The local beam-column elements under distributed load according to both Euler-Bernoulli and Timoshenko beam theories will be introduced. Finally, several numerical examples are presented to illustrate the accuracy and efficiency of the present study.

Co-rotational beam-column element under member loads

Unlike the traditional co-rotational formulation for beam-column elements which only consider the nodal loads at element ends, the distributed member loads are taken into account in this paper.

Figure 1 shows the deformations of the beam-column element in the global coordinate system (o^gx^gy^g), in which the member load is assumed as the uniformly distributed load. The local coordinate system (oxy) moving with the beam-column element is defined by the two end nodes, dividing the element total motions into two parts, that is, rigid body movement and pure deformation. The magnitude (p) and the direction (−y^g) of the distributed member load defined in the global coordinate system remain unchanged after the element deforms, but they are varying in the local coordinate system as the element deflects.

Figure 1.

Co-rotational beam-column element under member load.

As shown in Figure 1, the beam-column element moves from the initial position to the current position with the global displacement vector $d^{g} = {\begin{matrix} u_{1}^{g} v_{1}^{g} θ_{1}^{g} u_{2}^{g} v_{2}^{g} θ_{2}^{g} \end{matrix}}^{T}$ under the point loads $f^{g} = {\begin{matrix} f_{x 1}^{g} f_{y 1}^{g} m_{1}^{g} f_{x 2}^{g} f_{y 2}^{g} m_{2}^{g} \end{matrix}}^{T}$ and the uniformly distributed load −p due to gravity, in which $u_{j}^{g}$ and $v_{j}^{g}$ are the global displacements at node j (j = 1,2) along axis x^g and axis y^g, respectively; $θ_{j}^{g}$ is the global rotation at node j; $f_{xj}^{g}$ , $f_{yj}^{g}$ , and $m_{j}^{g}$ are the nodal forces corresponding to $u_{j}^{g}$ , $v_{j}^{g}$ , and $θ_{j}^{g}$ , respectively. The initial and current nodal coordinates in of the beam-column element in the global coordinate system are ${{}^{i}x}_{j}^{g} = {\begin{matrix} {{}^{i}x}_{j}^{g} & {{}^{i}y}_{j}^{g} \end{matrix}}^{T}$ and ${{}^{i + 1}x}_{j}^{g} = {\begin{matrix} {{}^{i + 1}x}_{j}^{g} & {{}^{i + 1}y}_{j}^{g} \end{matrix}}^{T}$ respectively, in which the subscript j is the node number; the superscripts i and i+1 means that the pertaining parameters belong to the initial and current times, respectively.

Then, the pure deformations produced in the local frame (oxy) can be given by

e = \sqrt{{({{}^{i + 1}x}_{2}^{g} - {{}^{i + 1}x}_{1}^{g})}^{2} + {({{}^{i + 1}y}_{2}^{g} - {{}^{i + 1}y}_{1}^{g})}^{2}} - \sqrt{{({{}^{i}x}_{2}^{g} - {{}^{i}x}_{1}^{g})}^{2} + {({{}^{i}y}_{2}^{g} - {{}^{i}y}_{1}^{g})}^{2}}

(1a)

θ_{1} = θ_{1}^{g} - α' + α

(1b)

θ_{2} = θ_{2}^{g} - α' + α

(1c)

α' = \arctan (\frac{{{}^{i + 1}y}_{2}^{g} - {{}^{i + 1}y}_{1}^{g}}{{{}^{i + 1}x}_{2}^{g} - {{}^{i + 1}x}_{1}^{g}})

(1d)

α = \arctan (\frac{{{}^{i}y}_{2}^{g} - {{}^{i}y}_{1}^{g}}{{{}^{i}x}_{2}^{g} - {{}^{i}x}_{1}^{g}})

(1e)

where ${{}^{i + 1}x}_{j}^{g} = {{}^{i}x}_{j}^{g} + u_{j}^{g}$ , ${{}^{i + 1}y}_{j}^{g} = {{}^{i}y}_{j}^{g} + v_{j}^{g}$ , $α$ , and $α'$ are the angles between axis-x^g and the initial and current configurations of the beam-column element respectively.

By taking a variation of equation (1), the relation between the virtual pure deformations and the virtual global displacements is

δ d = {[T]}^{T} {[L]}^{T} δ d^{g}

(2a)

{[T]}^{T} = [\begin{matrix} - 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 / L' & 1 & 0 & 1 / L' & 0 \\ 0 & 1 / L' & 0 & 0 & 1 / L' & 1 \end{matrix}], L' = L + e

(2b, c)

[L] = [\begin{matrix} [L_{e}] & 0 \\ 0 & [L_{e}] \end{matrix}]

(2d)

[L_{e}] = [\begin{matrix} \cos α' & - \sin α' & 0 \\ \sin α' & \cos α' & 0 \\ 0 & 0 & 1 \end{matrix}]

(2e)

in which $d = {\begin{matrix} e & θ_{1} & θ_{2} \end{matrix}}^{T}$ is the pure deformations in the basic system with only three degrees of freedom, while the corresponding force vector is $f = {\begin{matrix} P & M_{1} & M_{2} \end{matrix}}^{T}$ .

Nodal applied loads

When the external loads only contain nodal loads at the element ends, the relation between the basic and global force vectors can be derived through the principle of virtual work as follows

δ d^{g^{T}} f^{g} = δ d^{T} f

(3a)

f^{g} = [L] [T] f

(3b)

The previous derivation is the traditional element-independent co-rotational formulation for 2D beam-column element with only considering nodal loads. The basic to local transformation matrix [ T ] extends three basic forces to six forces in the local coordinate system. Note that the basic and local systems are coincident, but they have different numbers of degrees of freedom. The local to global transformation matrix [ L ] is completely determined by the angle between the local and global coordinate systems.

In fact, the force relation between basic and global systems in equation (3b) can be given directly through force equilibrium equations, in which the physical interpretation of the transformation from basic to local plays an important role. Specifically, the basic forces f can be regarded as applied loads at the free degrees of freedom of the beam-column element with a simple support as shown in Figure 2, and these forces do not change after deformation. The transformation from basic to local system is a procedure to derive reaction forces of the simple support and get rid of the simple support. Further, it should be noticed that the reaction forces are related to the deformed length.

Figure 2.

Basic beam-column element with a simple support.

Distributed applied member loads

To derive the exact co-rotational formulation for beam-column elements under distributed member loads, we can resort to the force equilibrium equations and the physical interpretation introduced in the last subsection. Figure 3 shows the beam-column element under nodal and distributed member loads in the basic system.

Figure 3.

Beam-column element under distributed load.

In this system, the beam-column element is derived based on the total Lagrangian description. Then, the total potential energy functional Π, which is the sum of strain energy U, work done by nodal forces W_f, and work done by gravity load W_p, can be given by

Π = U - W_{f} - W_{p}

(4a)

W_{f} = f^{T} d = Pe + M_{1} θ_{1} + M_{2} θ_{2}

(4b)

W_{p} = \int_{0}^{L} p_{x} u (x) dx + \int_{0}^{L} p_{y} v (x) dx

(4c)

in which $p_{x} = - p \cos α'$ and $p_{y} = - p \sin α'$ are the distributed member load along local axis-x and axis-y respectively, $u (x)$ and $v (x)$ are the corresponding displacement functions along local axes.

By taking the first variation of the total potential energy functional in equation (4a), it is easy to get the nodal forces in the basic system as follows

f = f_{U} + f_{p}

(5)

where $f_{U}$ is the nodal force vector due to strain energy U and $f_{p}$ is the equivalent nodal force vector of distributed member load, derived from W_p.

It can be seen from equation (4c) that $f_{p}$ is related to displacement functions. In the next section, the beam-column element formulations based on Euler-Bernoulli and Timoshenko theories in the basic system is derived, respectively.

Local beam-column element formulations

Several local beam-column elements considering distributed loads and containing the loads in the displacement or force interpolation functions have been proposed by some researchers. For example, Zhou and Chan (1997) proposed an effective element for second-order analysis of steel frame allowing for axial and lateral distributed loads, while the local second-order effect was considered at the same time. Also, the family of flexibility-based beam-column elements (Du et al., 2017; Neuenhofer and Filippou, 1997, 1998) are able to consider the distributed loads in the force fields or stress resultants.

This paper mainly focuses on the co-rotational formulation rather than the local elements. To illustrate the present study clearly, the proposed element is assumed as prismatic and elastic and neglects the warping and the local second-order effect. Then, based on the classical beam-column element with linear axial and cubic transverse displacement interpolations, the uniformly distributed load is exactly considered in the displacement functions herein.

To consider the distributed loads, the axial and transverse displacement interpolations in the basic system are assumed as quadric and quartic functions respectively in the following

u = \sum_{i = 0}^{2} a_{i} x^{i}

(6a)

v = \sum_{i = 0}^{4} b_{i} x^{i}

(6b)

As the local second-order effect is neglected, the axial and bending parts of the element are not coupled and then can be treated separately in the derivation. With respect to the axial part of the element, the pertaining equilibrium differential equation and the boundary conditions can be given by

EA \frac{d^{2} u}{d x^{2}} = p_{x}

(7a)

u = 0 at x = 0, u = e at x = L

(7b, c)

Through equations (7) and (6a), the coefficients a_i in equation (6a) can be solved and then the expression of the axial displacement function is given by

u = (\frac{e}{L} + \frac{L p_{x}}{2 EA}) x - \frac{p_{x}}{2 EA} x^{2}

(8)

In terms of the bending part of the element, two different beam elements based on Euler-Bernoulli and Timoshenko beam theories are studied respectively, by using the quartic displacement function in equation (6b).

Euler-Bernoulli beam theory

According to Euler-Bernoulli beam theory, the relation between rotation and transverse displacement at arbitrary point x of the beam is given by

θ = \frac{dv}{dx}

(9)

Then, the equilibrium differential equation of the beam under uniformly distributed load can be derived as follows

p_{y} = - \frac{dQ}{dx} = - \frac{d^{2} M}{d x^{2}} = EI \frac{d^{3} θ}{d x^{3}} = EI \frac{d^{4} v}{d x^{4}}

(10)

where Q and M are the shear force and moment of the beam element, respectively.

Also, the boundary conditions of the beam element are

v = 0 at x = 0, L

(11a, b)

\frac{dv}{dx} = θ_{1} at x = 0, \frac{dv}{dx} = θ_{2} at x = L

(11c, d)

in which E is the elastic Young’s modulus, A is the cross-sectional area, and I is the second moment of area.

Substituting equations (10) and (11) into equation (6b), the coefficients b_i can be solved and then the expression of the transverse displacement function is derived out in the following

\begin{matrix} v = θ_{1} x - (2 \frac{θ_{1}}{L} + \frac{θ_{2}}{L} - \frac{L^{2} p_{y}}{24 EI}) x^{2} \\ + (\frac{θ_{1}}{L^{2}} + \frac{θ_{2}}{L^{2}} - \frac{L p_{y}}{12 EI}) x^{3} + \frac{p_{y}}{24 EI} x^{4} \end{matrix}

(12)

Also, based on Euler-Bernoulli beam theory, the strain energy of the beam-column element can be given by

U = \frac{1}{2} \int_{0}^{L} {EA (\frac{du}{dx})}^{2} dx + \frac{1}{2} \int_{0}^{L} {EI (\frac{d^{2} v}{d x^{2}})}^{2} dx

(13)

Through equations (4), (5), (8), (12), and (13), the nodal forces due to strain energy and the equivalent nodal forces of the uniformly distributed load in the basic system can be derived out and given by

f_{U} = {\begin{matrix} \frac{EA}{L} e \\ \frac{EI}{L} (4 θ_{1} + 2 θ_{2}) \\ \frac{EI}{L} (2 θ_{1} + 4 θ_{2}) \end{matrix}}, f_{p} = {\begin{matrix} - \frac{L}{2} p_{x} \\ - \frac{L^{2}}{12} p_{y} \\ \frac{L^{2}}{12} p_{y} \end{matrix}}

(14a, b)

Then, the material stiffness matrix of the beam-column element in the basic system is

δ f_{U} = [k_{e}] δ d

(15a)

[k_{e}] = [\begin{matrix} \frac{EA}{L} & 0 & 0 \\ 0 & 4 \frac{EI}{L} & 2 \frac{EI}{L} \\ 0 & 2 \frac{EI}{L} & 4 \frac{EI}{L} \end{matrix}]

(15b)

Timoshenko beam theory

Based on Timoshenko beam theory, the transverse shear strain γ and the curvature κ can be given by

γ = \frac{dv}{dx} - θ, κ = - \frac{d θ}{dx}

(16a, b)

Besides, the equilibrium differential equations of Timoshenko beam are

p_{y} = - \frac{dQ}{dx} = - kGA \frac{d γ}{dx}

(17a)

Q = \frac{dM}{dx} = EI \frac{d κ}{dx}

(17b)

in which k is the sectional shear correction factor and G is the shear modulus of elasticity.

Through these equations and the boundary conditions given in equation (11), the expression of transverse displacement and rotation can be obtained as follows

\begin{matrix} v = \frac{(L - x) x (((2 + s) L - 2 x) θ_{1} - (sL + 2 x) θ_{2})}{2 (1 + s) L^{2}} \\ - \frac{p_{y} (L - x) x ((- 1 + s) s L^{2} + (- 1 + 7 s) Lx + (1 - 7 s) x^{2})}{24 (1 + s) EI} \end{matrix}

(18a)

\begin{matrix} θ = \frac{((- 2 + s) L + 3 x) x θ_{2} + (L + sL - 3 x) (L - x) θ_{1}}{(1 + s) L^{2}} \\ - \frac{2 p_{y} (- 1 + 5 s) (L - 2 x) (L - x) x}{12 (1 + s) EI} \end{matrix}

(18b)

in which $s = \frac{12 EI}{kGA L^{2}}$ .

The strain energy of the beam-column element based on Timoshenko beam theory is

\begin{matrix} U = \frac{1}{2} \int_{0}^{L} {EA (\frac{du}{dx})}^{2} dx + \frac{1}{2} \int_{0}^{L} {EI (\frac{d θ}{dx})}^{2} dx \\ + \frac{1}{2} \int_{0}^{L} kGA γ^{2} dx \end{matrix}

(19)

Following the derivation presented in the last subsection, the nodal forces due to strain energy and equivalent nodal forces of the uniformly distributed load in the basic system become

f_{U} = {\begin{matrix} \frac{EA}{L} e \\ \frac{EI ((4 + s) θ_{1} + (2 - s) θ_{2})}{(1 + s) L} \\ \frac{EI ((2 - s) θ_{1} + (4 + s) θ_{2})}{(1 + s) L} \end{matrix}}, f_{p} = {\begin{matrix} - \frac{L}{2} p_{x} \\ - \frac{L^{2}}{12} p_{y} \\ \frac{L^{2}}{12} p_{y} \end{matrix}}

(20a, b)

The effect of the transverse shear deformation is considered in f_U , while f _p are identical for both Euler-Bernoulli and Timoshenko beam theories. Therefore, the corresponding material stiffness matrix in the basic system is

δ f_{U} = [k_{e}] δ d

(21a)

[k_{e}] = [\begin{matrix} \frac{EA}{L} & 0 & 0 \\ 0 & \frac{4 + s}{1 + s} \frac{EI}{L} & \frac{2 - s}{1 + s} \frac{EI}{L} \\ 0 & \frac{2 - s}{1 + s} \frac{EI}{L} & \frac{4 + s}{1 + s} \frac{EI}{L} \end{matrix}]

(21b)

Transformation from basic to global systems

In the last section, the nodal force vectors of the beam-column element in the basic system based on Euler-Bernoulli and Timoshenko beam theories have been presented, respectively. To derive the exact co-rotational formulation of the beam-column element, the relation between basic and global systems, which is involved with two transformations, should be established.

The first transformation of the nodal forces is from basic to local systems. As pointed out in Section 2, it is a procedure to extend three element qualities into six qualities and remove the element constraints. In addition, the change in element length should be considered in this transformation. Thus, the nodal forces in the local system can be given by

f^{l} = f_{U}^{l} + f_{p}^{l}

(22a)

f_{U}^{l} = [T] f_{U}

(22b)

in which the transformation matrix [ T ] is given in equation (2b), making the local nodal forces due to strain energy, $f_{U}^{l}$ , affected by the change in element length.

Meanwhile, the equivalent nodal forces of distributed member load in the basic system after element deformation are

f_{p} = {\begin{matrix} - \frac{L'}{2} p'_{x} & - \frac{L'^{2}}{12} p'_{y} & \frac{L'^{2}}{12} p'_{y} \end{matrix}}^{T}

(23a)

p'_{x} = - p' \cos α', p'_{y} = - p' \sin α'

(23b, c)

in which the variables with a symbol “′” on the right side means that they belong to the current element configuration.

Assuming that the distributed member load is gravity, we have $p' L' = p L$ since the gravity of the beam-column element is constant. In addition, taking advantage of the equilibrium conditions, the equivalent nodal forces of gravity load in the local system can be derived in the following

f_{p}^{l} = {\begin{matrix} - \frac{L p_{x}}{2} & - \frac{L p_{y}}{2} & - \frac{L p_{y}}{12} L' & - \frac{L p_{x}}{2} & - \frac{L p_{y}}{2} & \frac{L p_{y}}{12} L' \end{matrix}}^{T}

(24)

The second transformation of the nodal forces is from local to global systems. The global nodal forces of the element can be expressed as

f^{g} = [L] f^{l} = [L] (f_{U}^{l} + f_{p}^{l}) = [L] [T] f_{U} + f_{p}^{g}

(25a)

\begin{matrix} f_{p}^{g} = [L] f_{p}^{l} = \\ {0 - \frac{L p}{2} - \frac{L p}{12} L' \cos α' 0 - \frac{L p}{2} \frac{L p}{12} L' \cos α'}^{T} \end{matrix}

(25b)

If the equivalent nodal forces of the gravity load, $f_{p}^{g}$ , are neglected, equation (25a) is identical to equation (3b) for the case when the applied loads are only nodal. Also, it should be noted that $L' \cos α' = {{}^{i}x}_{2}^{g} + u_{2}^{g} - {{}^{i}x}_{1}^{g} - u_{1}^{g}$ .

Then, by taking a variation of equation (25a), we have

\begin{matrix} δ f^{g} = δ [L] [T] f_{U} + [L] δ [T] f_{U} + [L] [T] δ f_{U} + δ f_{p}^{g} \\ = ([k_{U}^{g}] + [k_{p}^{g}]) δ d^{g} \end{matrix}

(26a)

[k_{U}^{g}] = [L] ([T] [k_{e}] {[T]}^{T} + [k_{g}]) {[L]}^{T}

(26b)

[k_{g}] = [\begin{matrix} 0 & \frac{M_{1} + M_{2}}{L'^{2}} & 0 & 0 & - \frac{M_{1} + M_{2}}{L'^{2}} & 0 \\ \frac{P}{L'} & 0 & 0 & - \frac{P}{L'} & 0 \\ 0 & 0 & 0 & 0 \\ 0 & \frac{M_{1} + M_{2}}{L'^{2}} & 0 \\ sym . & \frac{P}{L'} & 0 \\ 0 \end{matrix}]

(26c)

[k_{p}^{g}] = [\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{L p}{12} & 0 & 0 & - \frac{L p}{12} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ - \frac{L p}{12} & 0 & 0 & \frac{L p}{12} & 0 & 0 \end{matrix}]

(26d)

It can be seen that the global tangent stiffness matrix has two parts. The matrix $[k_{U}^{g}]$ is identical to the global tangent stiffness for the co-rotational method only considering nodal loads, while the matrix $[k_{p}^{g}]$ is produced by the variations of the equivalent nodal forces of gravity load due to element motions. Thus, if the variations of the equivalent moments due to member load are not properly considered, the eigenvalue buckling analysis using one element per member cannot accurately predict the buckling loads of the structures.

Numerical examples

To verify the accuracy and efficiency of the proposed exact co-rotational formulation for beam-column elements under distributed member loads, several numerical examples are presented in the section. From the last section, the expressions of $f_{p}^{g}$ and $[k_{p}^{g}]$ related to the member load are identical for both Euler-Bernoulli and Timoshenko beam theories. Thus, the following examples only adopt Euler-Bernoulli beam theory for demonstration.

To validate the proposed co-rotational formulation, the traditional method to deal with member loads is also investigated in these examples, which adopts the equivalent nodal loads of the member load computed according to the initial configuration of structure and omits their variations due to the large rigid body rotations. The global nodal forces due to member load in equation (25b) are constant in the whole analysis procedure and given by

f_{p}^{g} = {\begin{matrix} 0 & - \frac{L p}{2} & - \frac{L^{2} p}{12} \cos α & 0 & - \frac{L p}{2} & \frac{L^{2} p}{12} \cos α \end{matrix}}^{T}

(27)

There is no doubt that the pertaining tangent stiffness matrix due to the gravity load, $[k_{p}^{g}]$ is null.

Besides the geometrically nonlinear analyses, the eigenvalue buckling analyses for Examples 1 and 3 are also presented. The eigenvalue problem is generally expressed as

(K_{e} + λ K_{g}) ϕ = 0

(28)

where $K_{e}$ is the assembled elastic material stiffness matrix of the whole structure and $K_{g}$ is the assembled geomatic stiffness matrix of the whole structure under preloads; $ϕ$ is the buckling mode shape (eigenvectors); and $λ$ is the elastic buckling load factor (eigenvalues). It should be noted that the proposed geometric stiffness matrix $K_{g}$ is different from those of the traditional methods, leading to more accurate prediction of buckling load factors.

Example 1: Column under gravity load

As shown in Figure 4, a simple column under gravity load and horizontal concentrated load at the top is used to test the proposed co-rotational formulation in the large deflection analysis. In this example, the column is 100 cm in length (L), 1 cm² in area (A), 1 cm⁴ in second moment of area (I), and 2 × 10⁷ N/cm² in elastic Young’s modulus. The maximum values of the horizontal concentrated load and the pressure due to gravity load are 100 N and 300 N/cm, respectively, while this study does not consider the load applied sequence and therefore these loads are applied simultaneously to the column. Also, it should be noted that the concentrated force at the top is small and produces small deflection when it is applied on the column independently, with the purpose to initiate imperfection and sideway of the column.

Figure 4.

Column under gravity load.

The column is modeled by 1, 2, and 4 beam-column elements and then solved by the presented co-rotational method and the traditional load method, respectively. The results obtained by these different methods are plotted in Figure 5. As shown in Figure 5, the results of the present study are in good agreement with the traditional method when the column is modeled with four elements and can be regarded as reference solutions. For the column model with two elements, the results of both two methods are close to the reference solutions. However, significant differences between these two methods arise for the column using a single element, and it can be seen that the present study is evidently better than the traditional load method.

Figure 5.

Load versus displacement curves for cantilever column: (a) Load vs. displacement u, and (b) Load vs. displacement v.

The results of the eigenvalue buckling analyses by different models are provided in Table 1 which shows that the proposed formulation is more accurate than the traditional method when only one element per member is used.

Table 1.

Results of eigenvalue buckling analyses for Example 1.

	Buckling load factor (first mode)
	Oneelement	Twoelements	Fourelements
Presentstudy	0.533	0.557	0.532
Traditionalmethod	0.400	0.521	0.525

It can be seen that the one-element model using the present method can accurately predict the buckling load, though its results of nonlinear analysis gradually deviate from the reference solutions in post-buckling stage as the local linear beam-column element is limited to small local pure deformations. This problem can be improved by using higher performance beam-column elements allowing for large deformations. Regarding the one-element column model using the traditional load method, the results are far away from the reference solutions. This is because the equivalent moments are neglected according to the initial vertical configuration, which in fact arise after deformation and help the column against buckling. In this example, the advantage of the proposed co-rotational formulation for beam-column elements under gravity load is remarkable.

Example 2: Cantilever beam under gravity load

A cantilever beam under gravity load, which is a classic example for checking the capability of an element undergoing large displacements and large rotations, is shown in Figure 6. This model adopts the same geometric and material properties as those in Example 5.1. Different from the last problem of a column, there is no buckling behavior in this example and the gravity load is mainly resisted by the bending stiffness of the cantilever beam.

Figure 6.

Cantilever beam under gravity load.

Following the last example, different numbers of beam-column elements and geometrically nonlinear analysis formulations are used to simulate the cantilever beam under gravity load. The results obtained from these methods are plotted in Figure 7. In addition, the horizontal and vertical displacements corresponding to the gravity load p of 400 N/cm are listed in Table. 2.

Figure 7.

Load versus displacement curves for cantilever beam: (a) Load vs. displacement u, and (b) Load vs. displacement v.

Table 2.

Numerical results for cantilever beam with p = 400 N/cm.

	One element		Two elements		Four elements
	u (cm)	v (cm)	u (cm)	v (cm)	u (cm)	v (cm)
Present study	55.57	89.69	58.00	87.17	55.84	84.15
Traditional method	45.24	83.78	55.08	86.52	55.03	83.99

The performance of the present study and the traditional method are similar to the last example. The differences between the results obtained from these methods decrease with the increasing number of elements, while the present study is still better than the traditional load method when only one beam-column element is used to simulate this example. For the one-element model, the discrepancy of the present study comes from the local beam-column element, which could be reduced by considering the local second-order effect in the element formulation, such as the element proposed by Tang et al. (2015). In general, compared with the traditional method, the improvement of the presented co-rotational formulation under gravity load is also significant.

Example 3: Three-story frame under gravity load

The three-story frame shown in Figure 8 is investigated in this example. All members of this frame adopt the same square cross-section with 1 cm² in area (A) and 0.0833 cm⁴ in second moment of area (I), while their material properties are the same as those in the last two examples. The notional force F representing the frame imperfection is applied at the top left corner point of the frame, and its magnitude is 1000 N. Both beam and column members of the frame are subjected to gravity load p of 10 N/cm. All connections between members are rigid while the column bases are fixed.

Figure 8.

Three-story frame under gravity load.

Similar to the first example, the notional force F and the gravity load p are applied to the structure simultaneously during analysis process, whereas the load factor λ is defined by the current loads divided by input loads. Here, the frame models with one element and two elements per member are, respectively, solved by the proposed co-rotational formulation and the traditional load method. The horizontal displacement u and vertical displacement −v at Point A are plotted in Figure 9, while the results of eigenvalue buckling analyses are shown in Table 3.

Figure 9.

Load versus displacement curves for three-story frame: (a) Load vs. horizontal displacement u, and (b) Load vs. vertical displacement -v.

Table 3.

Results of eigenvalue buckling analyses for Example 3.

	Buckling load factor (first mode)
	One elementper member	Two elementsper member
Proposed method	0.305	0.295
Traditionalmethod	0.359	0.295

It can be seen from Figure 9 and Table 3 that the results of the present study and the traditional method agree well with each other when each frame member is divided into two elements. In terms of the models of one element per member, the results given by the present study are close to those by the two-elements per member models. It demonstrates that the proposed method is better than the traditional one.

Example 4: Plastic analysis for three-story frame under gravity load

The three-story frame presented in Example 3 is also used here to illustrate the performance of the proposed method considering material yielding by plastic hinge method. The yield stress of the material is assumed as 23,500 N/cm² and the elastic-perfectly plastic model is adopted in the plastic hinge method. The influence of axial force to the full yield surface is also neglected herein since it is small compared to the end moments. The plastic moment capacity of all members is 5875 N cm. Once the bending moment is larger than the plastic moment, a plastic hinge will form at the member end. The well-known arc-length method for nonlinear solution is used to trace the full load-displacement curve.

The horizontal displacements u at Point A of the three-story frame are plotted in Figure 10. It shows that the results obtained from the proposed method using one element per member is more accurate than the traditional method. The deviation compared to the models using two elements per member is due to that the local beam-column element formulation is based on the assumption of small strains without consideration of local second-order effect. This example clearly shows the advantage of the proposed co-rotational framework. To demonstrate the merit of the proposed method, only the conventional beam-column elements are used for easy comparison. To increase the analysis accuracy, more complicated second-order beam-column elements can be introduced.

Figure 10.

Load versus displacement curves for three-story frame.

Conclusion

The conventional co-rotational formulation adopted in eigenvalue buckling and second-order direct analyses cannot accurately consider the distributed member loads, especially for member undergoing large rotations. In this paper, an exact and consistent co-rotational formulation is first proposed for both Euler-Bernoulli and Timoshenko beam-column elements under distributed member loads. This work will significantly improve the accuracy and reliability of the commonly used structural buckling methods in practical applications. Several examples have been used to verify the accuracy of the proposed method. By contrast, the traditional method cannot consider the variations of the equivalent nodal moments due to large rotations. From this study, the following conclusions can be made:

Several elements per member are required for the structures subjected to distributed member loads undergoing large rotations when using conventional co-rotational formulations.

The element independent co-rotational (EICR) method provides a unified approach for buckling analysis using both Euler-Bernoulli and Timoshenko beam-column elements.

The proposed new co-rotational formulation significantly improve the accuracy of both eigenvalue buckling and second-order direct analyses.

It should be pointed out that this paper mainly focuses on the co-rotational framework. To predict the nonlinear behaviors more accurate, more advanced local element including initial member imperfections should be used. For practical application of the proposed method, the co-rotational formulations for 3D beam-column elements will be presented in our future work.

Footnotes

Authors’ Note

Yao-Peng Liu is also affiliated with NIDA Technology Company Limited, Hong Kong Science Park, Shatin, N.T., Hong Kong, China.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors are grateful for financial support from the Research Grant Council of the Hong Kong SAR Government on the project “Second-order analysis and design of scaffolds with scissor braces and allowing for kink imperfections (PolyU).” The first author is also grateful for the financially supported by the National Natural Science Foundation of China (52008094) and Natural Science Foundation of Jiangsu Province (SBK20200402).

ORCID iDs

Yao-Peng Liu

Siu-Lai Chan

References

Argyris

(1982) An excursion into large rotations. Computer Methods in Applied Mechanics and Engineering 32(1): 85–155.

Bathe

(2006) Finite Element Procedures. USA: Prentice Hall, Pearson Education, Inc.

Bathe

Bolourchi

(1980) A geometric and material nonlinear plate and shell element. Computers and Structures 11(1–2): 23–48.

Bathe

Dvorkin

(1983) Our discrete-Kirchhoff and isoparametric shell elements for nonlinear analysis-an assessment. Computers and Structures 16(1–4): 89–98.

Battini

(2002) Co-rotational beam elements in instability problems. PhD Thesis, Department of Mechanics, Royal Institute of Technology, Sweden.

Chan

Zhou

(1994) Pointwise equilibrating polynomial element for nonlinear analysis of frames. Journal of Structure Engineering 120(6): 1703–1717.

Chan

Zhou

(1995) Second-order elastic analysis of frames using single imperfect element per member. Journal of Structure Engineering 121(6): 939–945.

Crisfield

(1990) A consistent co-rotational formulation for non-linear, three-dimensional, beam-elements. Computer Methods in Applied Mechanics and Engineering 81(2): 131–150.

Liu

Chan

(2017) A second-order flexibility-based beam-column element with member imperfection. Engineering Structures 143: 410–426.

10.

Dvorkin

Bathe

(1984) A continuum mechanics based four-node shell element for general non-linear analysis. Engineering computations 1(1): 77–88.

11.

Felippa

Haugen

(2005) A unified formulation of small-strain corotational finite elements: I. Theory. Computer Methods in Applied Mechanics and Engineering 194(21): 2285–2335.

12.

Hsiao

(1987) Nonlinear analysis of general shell structures by flat triangular shell element. Computers and Structures 25(5): 665–675.

13.

Hsiao

Hung

(1989) Large-deflection analysis of shell structure by using corotational total lagrangian formulation. Computer Methods in Applied Mechanics and Engineering 73(2): 209–225.

14.

Izzuddin

(2005) An enhanced co-rotational approach for large displacement analysis of plates. International Journal for Numerical Methods in Engineering 64(10): 1350–1374.

15.

Izzuddin

Liang

(2016), Bisector and zero-macrospin co-rotational systems for shell elements. International Journal for Numerical Methods in Engineering 105(4): 286–320.

16.

Jeon

Lee

, et al. (2015) The MITC3+shell element in geometric nonlinear analysis. Computers and Structures 146: 91–104.

17.

Jiang

Chernuka

Pegg

(1994) A co-rotational, updated Lagrangian formulation for geometrically nonlinear finite element analysis of shell structures. Finite Elements in Analysis and Design 18(1–3): 129–140.

18.

Khosravi

Ganesan

Sedaghati

(2007) Corotational non-linear analysis of thin plates and shells using a new shell element. International Journal for Numerical Methods in Engineering 69(4): 859–885.

19.

Khosravi

Ganesan

Sedaghati

(2008) An efficient facet shell element for corotational nonlinear analysis of thin and moderately thick laminated composite structures. Computers and Structures 86(9): 850–858.

20.

Izzuddin

Vu-Quoc

, et al. (2017) A 3-nodeco-rotational triangular elasto-plastic shell element using vectorial rotational variables. Advanced Steel Construction 13(3): 206–240.

21.

Vu-Quoc

, et al. (2018) A nine-node corotational curved quadrilateral shell element for smooth, folded, and multishell structures. International Journal for Numerical Methods in Engineering 116(8): 570–600.

22.

Meek

Tan

(1983) Large deflection and post-buckling analysis of two and three dimensional elastic spatial frames. Report, Department of Civil Engineering, University of Queensland, Australia.

23.

Neuenhofer

Filippou

(1997) Evaluation of nonlinear frame finite-element models. Journal of Structure Engineering 123(7): 958–966.

24.

Neuenhofer

Filippou

(1998) Geometrically nonlinear flexibility-based frame finite element. Journal of Structure Engineering 124(6): 704–711.

25.

Nour-Omid

Rankin

(1991) Finite rotation analysis and consistent linearization using projectors. Computer Methods in Applied Mechanics and Engineering 93(3): 353–384.

26.

Rankin

Brogan

(1986) An element independent corotational procedure for the treatment of large rotations. Journal of Pressure Vessel Technology 108(2): 152–165.

27.

Skallerud

Haugen

(1999) Collapse of thin shell structures-stress resultant plasticity modelling within a co-rotated ANDES finite element formulation. International Journal for Numerical Methods in Engineering 46(12): 1961–1986.

28.

Surana

(1983) Geometrically nonlinear formulation for the curved shell elements. International Journal for Numerical Methods in Engineering 19(4): 581–615.

29.

Tang

Liu

Chan

, et al. (2019) An innovativeco-rotational pointwise equilibrating polynomial element based on Timoshenko beam theory for second-order analysis. Thin-Walled Structures 141: 15–27.

30.

Tang

Zhou

Chan

(2015) Nonlinear beam-column element under consistent deformation. International Journal of Structural Stability and Dynamics 15(05): 1450068.

31.

Tang

Zhou

Chan

(2016) Geometrically nonlinear analysis of shells by quadrilateral flat shell element with drill, shear and warping. International Journal for Numerical Methods in Engineering 108(10): 1248–1272.

32.

Tang

Zhou

Chan

(2017) A simplified co-rotational method for quadrilateral shell elements in geometrically nonlinear analysis. International Journal for Numerical Methods in Engineering 112(11): 1519–1538.

33.

Yang

Chiou

(1987) Rigid body motion test for nonlinear analysis with beam elements. Journal of Engineering Mechanics 113(9): 1404–1419.

34.

Zhou

Chan

(1997) Second-order analysis of slender steel frames under distributed axial and member loads. Journal of Structure Engineering 123(9): 1187–1193.

35.

Zienkiewicz

Taylor

(2005) The Finite Element Method for Solid and Structural Mechanics. UK: Elsevier.