A Co-Rotational Meshfree Method for the Geometrically Nonlinear Analysis of Structures

: This paper presents a co-rotational beam formulation, which is used for geometric nonlinear analysis with the differential reproducing kernel (DRK) approximation collocation method. The present formulation, based on the Timoshenko beam hypothesis, is capable of effectively solving geometrically nonlinear problems such as large deformation, postbuckling, lateral buckling, and snap-through problems. The kinematics have been constructed with the concept of co-rotational formulation adopted in the ﬁnite element method (FEM). A meshfree method based on the differential reproducing kernel (DRK) approximation collocation method, combined with the Newton–Raphson method, is employed to solve the strong forms of the geometrically nonlinear problems. The DRK method takes full advantage of the meshfree method. Moreover, only a scattered set of nodal points is necessary for the discretization. No elements or mesh connectivity data are required. Therefore, DRK will be able to completely circumvent the problems of mesh dependence and mesh distortion. The effectiveness of this study and its performance are shown through several numerical applications.


Introduction
The geometrically nonlinear analysis of a beam structure undergoing large displacement remains a strong topic of interest, with applications in various fields such as civil engineering over the past few decades. Various formulation strategies and numerical procedures using FEM have been developed to deal with the large deformation problems . For the large displacements and finite rotation problems of a geometrically nonlinear analysis, the description of the kinematics can be regarded as the dominant factor. The formulations used to describe the kinematics of beam structures in the literature might be divided into three categories: total Lagrangian (TL) formulation [2,5,9,11,12,15], updated Lagrangian (UL) formulation [3,7,18], and co-rotational (CR) formulation [4,8,10]. The main difference between these formulations lies in their definitions of the reference of the coordinates adopted in FEM. In the total Lagrangian formulation, the initial undeformed configuration is adopted as the reference configuration. Consequently, the system equations are defined with respect to the fixed global coordinate system through the analysis. In contrast, the updated Lagrangian formulation adopted the local coordinate system, updated with the current configuration, as the reference configuration-that is, the system equations were constructed at the current deformed beam configuration. Generally, the system equations expressed in the updated Lagrangian formulation will be simpler than the corresponding equations in the total Lagrangian formulation. However, if the displacement from the current configuration to the last equilibrium configuration is large, it will be in violation of the basic assumption and will lead to unreliable results. To circumvent this problem, the co-rotational formulation was introduced. In the co-rotational framework, there are global coordinates and co-rotational coordinates, kept fixed at the deformed beam segment during the analysis procedure. The co-rotational coordinate rotates and translates with each element but does not deform with it. Therefore, the total motions of the element can be divided into the rigid body part and the purely deformational part.
For the finite element formulation, another important issue is the choice of the primary interpolated variables. No matter which of the above formulations is adopted to describe the kinematics, most of the large deformation analyses performed with the FEM use the conventional displacement approach with the principle of virtual work or minimum potential energy [2][3][4][5]. Both displacements and rotations or rotations alone are taken as the interpolated variables in the finite element implementation. Hence, the choice of the primary interpolated variables is also an important issue in the finite element formulation. In standard finite element formulations, the stress resultants solved by the equilibrium equations, and those computed by the constitutive equations, are not equal since the differentiated quantities are, as a rule, one order less accurate than the quantity itself. When formulations based on displacement and rotation were employed, the strains had to be computed by differentiating toward the assumed kinematic field. The internal forces and moments will be evaluated by the constitutive law with the recovered strain field. However, the accuracy of the differentiated quantities may decline. In particular, employing low-order shape functions to solve the shear flexible models may yield locking problems. Therefore, the performance of the element may deteriorate rapidly while dealing with shear flexible models such as thinner beam problems. However, the numerical effect can be corrected with the aid of reduced integration. Still, the displacement-rotation-based formulations may acquire inaccurate stress fields due to recovery by the differentiation procedure, even if the displacement fields can be solved accurately. To circumvent the numerical problem, some researchers have taken strain vectors [13] or stress resultants [16] as the primary variables. For instance, Zupan and Saje [13] presented a formulation in terms of the strain vectors such that the additive interpolation has no restrictions. On the contrary, Cannarozzi and Molari [16] took the stress vectors as the only interpolated degree of freedom. In this way, the displacement can be obtained by integration of the kinematic equations. Magisano et al. [21] developed a large rotation finite element analysis of three-dimensional beams based on the Mixed Integration Point (MIP) Newton iterative scheme [22], in which the co-rotational nodal rotations defined by incremental nodal rotation vectors are interpolated for the evaluation of the nonlinear strains. In the above studies, whether the stress or strain vectors have been employed to be the interpolated variables in the finite element formulation, it still deals with the weak form based on the principle of virtual work or the principle of minimum potential energy.
As opposed to FEM, adopted to solve numerical analysis, in this study a meshfree method based on the Differential Reproducing Kernel (DRK) Approximation collocation method combined with the Newton-Raphson method was adopted to analyze geometrically nonlinear problems. The kinematics of beam structures undergoing large deformation are derived with the concepts of the co-rotational framework. It is well-known that the primary characteristic of meshfree methods is that the domain is completely discretized by a scattered set of nodal points instead of elements. However, this naturally circumvents mesh problems may result in numerical problems in FEM. Therefore, in the present scheme the solution of nonlinear equilibrium equations such as displacements, rotations, and force resultants of a beam that has undergone large deformation will be directly solved toward the strong form of the governing equations. The DRK implementation of the threedimensional co-rotational beam formulation, leading to differentiation of the assumed kinematic field, is not necessary while evaluating strains, internal forces, and moments. This also shows that the primary interpolated variables will not be the dominant factor in the present scheme.
Onate et al. [29] adopted weighted least squares interpolations to develop a finite point method to analyze elasticity problems. Yang et al. [30] proposed a radial basis collocation method with a strong-form generalized displacement control method to solve geometrically nonlinear problems. Most point collocation methods are based on the idea of reproducing kernel (RK)approximations. However, conventional RK methods compute the derivatives of the shape functions by differentiating toward the shape functions directly. In particular, this may lead to complex expressions and computations while computing high-order derivatives. Several methods were developed to solve the complicated problems. Wang et al. [31] employed a trivial numerical implementation of second-order smoothed gradients to develop a superconvergent gradient smoothing meshfree collocation method. Wang et al. [32], Yang et al. [33], and Yeh [34] developed the DRK method for one-dimensional and two-dimensional elasticity analysis.
In comparison with the conventional RK approximations, the DRK approximations involve the calculation of the derivatives of shape functions in the conventional RK. In the DRK approximations, the derivatives will be determined by the differential reproducing conditions instead of differentiating toward the shape functions directly in the conventional RK. In the present paper, the meshfree method based on the DRK approximations proposed by Wang et al. [32,33] and Yeh [34] is adopted to solve the geometrically nonlinear problems. In addition to the advantage of having no mesh, the DRK method also allows for easy computation of the derivative of the shape function. This characteristic enables the DRK approximations to construct system equations easily and directly from the governing equations. The aim of the present paper is to develop a three-dimensional co-rotational beam formulation with the aid of the DRK approximations combined with the Newton-Raphson method. The effectiveness of this approach and its performance are shown by numerical examples.

Basic Assumptions
The following assumptions are made in the derivation of the three-dimensional corotational beam formulation:

1.
The Timoshenko beam hypothesis is valid.

2.
The theory is limited to elasticity problems. 3.
The twist rate and unit extension of the centroid axis of the beam are uniform. 4.
The displacements and rotations of the beam can be large.

5.
The out-of-plane warping of the cross section is excluded.

Coordinate System
In order to describe the geometry and kinematics of the beam in the undeformed and deformed configurations, we defined two sets of right-handed rectangular Cartesian coordinate systems: 1. Fixed global coordinates: As shown in Figure 1, the global coordinate, which remains fixed during the entire analysis, describes the displacements, rotation, and curvature of the structure. The governing equations of the system are also constructed with respect to the global coordinates.
2. Co-rotational local rectangular Cartesian coordinates: The DRK collocation method adopts a scattered set of nodal points to fulfill the discretization instead of the element adopted in FEM. Hence, a set of co-rotational local coordinates has to be defined, associated with the nodal points collocated on the centroidal axis of each cross section. The origin of the co-rotational local coordinate system is tied to the centroid of the cross section. The x o and y o axes of the co-rotational local coordinates are taken as the principal directions of the cross section. In addition, the z o axis is selected to coincide with the normal direction of the cross section. In this way, this co-rotational local coordinate will be co-rotated with the beam undergoing large deformation. cretization instead of the element adopted in FEM. Hence, a set of co-rotational local coordinates has to be defined, associated with the nodal points collocated on the centroidal axis of each cross section. The origin of the co-rotational local coordinate system is tied to the centroid of the cross section. The xo and yo axes of the co-rotational local coordinates are taken as the principal directions of the cross section. In addition, the zo axis is selected to coincide with the normal direction of the cross section. In this way, this co-rotational local coordinate will be co-rotated with the beam undergoing large deformation.

Kinematics of Deformed Beam
From the kinematic assumptions made in the present paper, the deformations of the beam can be described by the displacements and rotation of the beam centroid axis. As shown in Figure 1, the position vector of an arbitrary point P on the beam centroidal axis at the undeformed and deformed configuration can be expressed as follows: where eX, eY, and eZ are the unit basis vector of the global coordinates; S and s are the arc length of beam centroidal axis in undeformed and deformed configuration, respectively. Moreover, Xo(S), Yo(S), Zo(S) and X(s), Y(s), Z(s) are the coordinate components of the position vector associated with undeformed and deformed configuration.
From the geometric relationship between the global and co-rotational local coordinates at deformed configuration, Equation (3) can be given as follows: The relationship between the global and co-rotational local coordinate associated with each nodal point by collocation point method on the deformed beam axis can be represented by the rotation matrix R(s): where R is the rotation matrix defined according to the Rodrigues formula [35], as follows:

Kinematics of Deformed Beam
From the kinematic assumptions made in the present paper, the deformations of the beam can be described by the displacements and rotation of the beam centroid axis. As shown in Figure 1, the position vector of an arbitrary point P on the beam centroidal axis at the undeformed and deformed configuration can be expressed as follows: And where e X , e Y , and e Z are the unit basis vector of the global coordinates; S and s are the arc length of beam centroidal axis in undeformed and deformed configuration, respectively.
The relationship between the global and co-rotational local coordinate associated with each nodal point by collocation point method on the deformed beam axis can be represented by the rotation matrix R(s): where R is the rotation matrix defined according to the Rodrigues formula [35], as follows: The rotation ϑ is a vector that lies on the axis of rotation, with its length equal to the angle of rotation. The rotation vector ϑ can be written as a column matrix composed of components ϕ, χ, and ψ in the fixed global coordinates. Furthermore, Θ is a skewsymmetrical matrix composed from the rotation vector. For the convenience of later derivation, the curvature vector is defined to represent the bending and twisting of the beam in terms of the co-rotational local coordinate as follows: where ω x , ω y , and ω z are curvatures with respect to the x, y, and z axes of the co-rotational coordinates, respectively. Hence, the change rates of unit vectors with respect to the chord length of beam axis at deformed configuration can be expressed as follows: where  (4) and (7), the curvature matrix can be represented in terms of the rotation matrix as follows: Then, substituting Equation (5) into Equation (8) will yield the curvature of the beam at the deformed configuration as follows: The change rate of rotation with respect to the arc length at deformed configuration can be obtained by solving Equations (9)-(11) as follows: Appl. Sci. 2021, 11, 6647 6 of 20 In addition, from Equations (3)-(5), the relationships between the fixed global coordinate and the rotation vector ϑ at the deformed configuration can be expressed in terms of the rotation vector as follows:

Equilibrium and Constitutive Equations
As shown in Figure 2, the fulfilment of equilibrium of the beam segment with distributed forces and moments has to satisfy the following conditions: In addition, from Equations (3)-(5), the relationships between the fixed global coordinate and the rotation vector J at the deformed configuration can be expressed in terms of the rotation vector as follows:

Equilibrium and Constitutive Equations
As shown in Figure 2, the fulfilment of equilibrium of the beam segment with distributed forces and moments has to satisfy the following conditions: The co-rotational scheme means the equilibrium condition can be constructed directly with respect to the local coordinate system. Hence, Equations (21) and (22) can be given by the following transformation: where G XYZ and G xyz represent any vector with reference to the global coordinate and local coordinate, respectively: In the present study, we established the constitutive equations with the Timoshenko beam hypothesis. Therefore, the constitutive laws between the external forces and deformation can be expressed as follows: where E is the Young's modulus, A is the cross-sectional area, J is the torsional constant, G is the shear modulus, ε is the axial strain of beam axis, γ x and γ y are the average shear Appl. Sci. 2021, 11, 6647 7 of 20 strain, k x and k y are the shear modificatory factors, and I xx and I yy are the moments of inertia about the x and y axes, respectively.

Governing Equations
The above equations have been constructed from the arc length of the beam axis, s, at the deformed configuration. For the convenience of later derivation, the equations can be represented with respect to the arc length of the undeformed beam axis, S, by the relationship of the chord length between the undeformed and deformed configuration as follows: Substituting Equation (29) into the kinematics and equilibrium equations leads to the governing equations: where The 12 nonlinear ordinary differential equations, Equations (30)-(41), are derived to solve the spatial beam with naturally curved and twisted beams undergoing large deformations.

Reproducing Kernel Approximation
Most numerical methods transform a differential equation into the integration form while solving it, such as the Galerkin method [23,27]. However, the DRK approximations adopted in this paper can circumvent the numerical mesh for numerical integration. The DRK collocation method is developed based on the DRK approximation [32][33][34]. The DRK approximation can be derived from the kernel function with a reproducing condition or the moving least square method. Two kinds of derivations represent different characteristics of the DRK approximation. In this section, a brief introduction is given with regard to the moving least square method. Consider a set of N nodal points scattered in a domain Ω and let x* be the reference point. The DRK approximation of a variable u(x) in the neighborhood of the reference point x*, denoted by u h (x), is is a vector of undetermined coefficients. Hence, the error in the weighted least-squares of N nodes at neighborhood with respect to the reference point x* can be expressed as follows: where w a (x * − x) is the local weight function and Ne is the number of neighboring nodes of x*. The index a is the radius of support and u(x j ) is the approximation of u(x) at node x j . By minimizing the error of the weighted least-square, the approximation u h (x) can be obtained as follows: where is the shape function that can be defined as follows: The approximation u h (x) in Equation (47) is a local approximation in the neighborhood with respect to the refence point x*. However, the global approximation can be fulfilled by the moving least-square method by imposing x = x* such that point x* will move with the point that has to be computed. In this way, the global approximation can be expressed as follows:

Differential Reproducing Kernel Approximation
Different from conventional RK approximants, the derivatives of the shape function in DRK are computed by differentiating toward the reproducing condition instead of directly with a shape function. Hence, no special treatment is needed in DRK while computing a different order derivative of the shape function. The derivative of the reproducing kernel function can be obtained by differentiating Equation (47) directly: Appl. Sci. 2021, 11, 6647 9 of 20 By moving the least squares method and imposing x = x*, the derivatives of the global approximation can be obtained as follows: Thus, the procedures for computing the different order derivative in DRK are the same. In addition, by imposing the point collocation method with the moving least squares method, a boundary valued problem can be expressed as in the following equations: where Ω and ∂Ω are the domain and boundary of the problem; L Ω and L ∂Ω are linear differential operators. If the number of nodes in the domain is N p , the number of nodes on the boundary is N b . Equations (52) and (53) lead to The variables in the above equation are unknown, so the least squares method is adopted to find the optimum solutions.
Similar to most meshfree methods, the influence of the nodal values of the neighboring nodes is estimated by the weight function. Therefore, domain Ω I should contain the nodal point x * itself and also make the shape function to be zero outside. In the following analysis, the cubic spline is adopted to be the weight function as follows: where the normal distance d = x j − x * /a and a is the radius of support. However, the radius of support also has to be small enough to exhibit the local character and cannot be too small to yield an ill-conditioned problem at the same time. Hence, it is necessary to determine the optimum range of radius of support. The present scheme is capable of tuning with the density of nodes automatically such that the radius of support can exhibit the local character and not lead to an ill-conditioned problem. For a point at x * , the number of neighboring nodes, Ne, is determined by the order of basis functions. Then we can sort the neighboring nodes by the distance to the point x * . The distance of the Ne-th neighboring node is adopted as the radius of support. In this way, the radius of support can be tuned automatically with the density of nodes instead of being fixed.

Strategy for Incremental-Iterative Procedure
An incremental-iterative method combining the DRK approximations with the Newton-Raphson method was adopted to solve the 12 coupled, nonlinear ordinary differential equations, Equations (30)-(41). In the present incremental-iterative procedure, force or displacement will be adopted as the prescribed increment. Whichever kind of increment is adopted, the solution from the preceding iteration will be taken as the new initial guess for the next one. Let X k , Y k , Z k , ϕ k , χ k , ψ k , M x,k , M y,k , M z,k , Q x,k , Q y,k , and Q z,k denote the solutions of the Kth incremental step. The incremental-iterative steps adopted in the present scheme are as shown in Figure 3. Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 21 Figure 3. The incremental-iterative steps adopted in the present scheme.

Numerical Examples
Five numerical examples, including large deformation, postbuckling of two-dimensional beam structures, stability, and lateral buckling problems of three-dimensional beam structures, are presented to show the effectiveness and performance of the proposed three-dimensional co-rotational beam formulation.

Cantilever Curved Beam Subjected to End Moment
Hence, a curved beam with an initial circular shape will bend into a straight beam when the end moment, , is applied. If the end moment reaches L , the beam will bend into a circle in the opposite direction.

Numerical Examples
Five numerical examples, including large deformation, postbuckling of two-dimensional beam structures, stability, and lateral buckling problems of three-dimensional beam structures, are presented to show the effectiveness and performance of the proposed threedimensional co-rotational beam formulation.

Cantilever Curved Beam Subjected to End Moment
Hence, a curved beam with an initial circular shape will bend into a straight beam when the end moment, M = 2πEI L , is applied. If the end moment reaches M = 4πEI L , the beam will bend into a circle in the opposite direction. and more nodal points can effectively improve the solution accuracy. Moreover, the higher degree basis functions will achieve better solution accuracy even if the same number of nodal points is adopted. A combination of quadratic basis functions and five neighboring nodes was adopted in the following examples.
The moment-displacement curves of the curved cantilever beam are plotted in Figure 7, solved by 31 nodal points with 10 loading increments. The number of iterations of each increment is less than three. As can be seen, the present theory can provide accurate solutions. To find the optimal number of neighboring nodes and order of basis function, the L2 error for a combination of neighboring nodes 3, 5, and 7 with respect to the basis functions of quadratic, cubic, and quartic was plotted in Figure 5. It shows that the optimal relationship between the number of the neighboring nodes and basis functions can be given as follows: where Ne is the number of the neighboring nodes and m is the number of basis functions. The computational time of the above combinations shown in Figure 5 was less than 10 s. Additionally, the above optimal combination of neighboring nodes with order of basis function was adopted to plot Figure 6, which shows the convergence corresponding to the number of nodal points of 21, 31, 41, and 51 with quadratic, cubic, and quartic basis functions, respectively. It should be noted that h, shown in Figure 6, means the spacing of nodal points. The results show that the DRK approximations with a higher degree basis and more nodal points can effectively improve the solution accuracy. Moreover, the higher degree basis functions will achieve better solution accuracy even if the same number of nodal points is adopted. A combination of quadratic basis functions and five neighboring nodes was adopted in the following examples. To find the optimal number of neighboring nodes and order of basis function, the L2 error for a combination of neighboring nodes 3, 5, and 7 with respect to the basis functions of quadratic, cubic, and quartic was plotted in Figure 5. It shows that the optimal relationship between the number of the neighboring nodes and basis functions can be given as follows: where Ne is the number of the neighboring nodes and m is the number of basis functions. The computational time of the above combinations shown in Figure 5 was less than 10 s. Additionally, the above optimal combination of neighboring nodes with order of basis function was adopted to plot Figure 6, which shows the convergence corresponding to the number of nodal points of 21, 31, 41, and 51 with quadratic, cubic, and quartic basis functions, respectively. It should be noted that h, shown in Figure 6, means the spacing of nodal points. The results show that the DRK approximations with a higher degree basis and more nodal points can effectively improve the solution accuracy. Moreover, the higher degree basis functions will achieve better solution accuracy even if the same number of nodal points is adopted. A combination of quadratic basis functions and five neighboring nodes was adopted in the following examples. The moment-displacement curves of the curved cantilever beam are plotted in Figure 7, solved by 31 nodal points with 10 loading increments. The number of iterations of each increment is less than three. As can be seen, the present theory can provide accurate solutions.

A cantilever Beam Initially Curved to π/4 Subjected to a Concentrated Tip Load
As shown in Figure 8, a cantilever beam initially curved to an arch of radius 100 m with 45-degree bending was subjected to a concentrated end load. This example was originally proposed by Bathe [1] and has also been adopted by a number of other researchers [2,3,12,13,21]. The beam has a unit square cross section and lies in the X-Y plane with initial tip position 70.71 m, 29.29 m, 0 m of the fixed X-Y-Z global coordinates. The material was linear and elastic, with the following properties: Elastic modulus, E = 10 7 N/m 2 , and Shear modulus, G = 5 × 10 6 N/m 2 . In the above reference, the beam was modeled by straight and curved beam elements with the nodes on the centroidal axis of the beam. In the present study, DRK was fulfilled by nodal points to simulate a beam undergoing large deformation. Figure 9, solved by 31 nodal points with 10 loading increments, shows the evolution of tip displacement as a function of the applied load modulus. The number of iterations is fewer than five within each increment. Table 1 shows the free end position of the initial curved cantilever beam under out-of-plane force. A comparison of the free end position from these references is also given in Table 1. The results show that the present theory can provide accurate solutions. The moment-displacement curves of the curved cantilever beam are plotted in Figure 7, solved by 31 nodal points with 10 loading increments. The number of iterations of each increment is less than three. As can be seen, the present theory can provide accurate solutions.

A cantilever Beam Initially Curved to π/4 Subjected to a Concentrated Tip Load
As shown in Figure 8, a cantilever beam initially curved to an arch of radius 100 m with 45-degree bending was subjected to a concentrated end load. This example was originally proposed by Bathe [1] and has also been adopted by a number of other researchers [2,3,12,13,21]. The beam has a unit square cross section and lies in the X-Y plane with initial tip position 70.71 m, 29.29 m, 0 m of the fixed X-Y-Z global coordinates. The material was linear and elastic, with the following properties: Elastic modulus, E = 10 7 N/m 2 , and Shear modulus, G = 5 × 10 6 N/m 2 . In the above reference, the beam was modeled by straight and curved beam elements with the nodes on the centroidal axis of the beam. In the present study, DRK was fulfilled by nodal points to simulate a beam undergoing large deformation. Figure 9, solved by 31 nodal points with 10 loading increments, shows the evolution of tip displacement as a function of the applied load modulus. The number of iterations is fewer than five within each increment. Table 1 shows the free end position of the initial curved cantilever beam under out-of-plane force. A comparison of the free end position from these references is also given in Table 1. The results show that the present theory can provide accurate solutions.

A cantilever Beam Initially Curved to π/4 Subjected to a Concentrated Tip Load
As shown in Figure 8, a cantilever beam initially curved to an arch of radius 100 m with 45-degree bending was subjected to a concentrated end load. This example was originally proposed by Bathe [1] and has also been adopted by a number of other researchers [2,3,12,13,21]. The beam has a unit square cross section and lies in the X-Y plane with initial tip position 70.71 m, 29.29 m, 0 m of the fixed X-Y-Z global coordinates. The material was linear and elastic, with the following properties: Elastic modulus, E = 10 7 N/m 2 , and Shear modulus, G = 5 × 10 6 N/m 2 . In the above reference, the beam was modeled by straight and curved beam elements with the nodes on the centroidal axis of the beam. In the present study, DRK was fulfilled by nodal points to simulate a beam undergoing large deformation.

Lateral Buckling of a Cantilever Beam
This example presents the lateral buckling of a cantilever beam with a narrow rectangular cross section. As shown in Figure 10, a transverse force, Fz, is loaded on the direction of the strong axis of the cross section. The geometry and material properties of the   Table 1 shows the free end position of the initial curved cantilever beam under out-of-plane force. A comparison of the free end position from these references is also given in Table 1. The results show that the present theory can provide accurate solutions.    [21] 4 (CR-Incremental) --53.47

Lateral Buckling of a Cantilever Beam
This example presents the lateral buckling of a cantilever beam with a narrow rectangular cross section. As shown in Figure 10, a transverse force, Fz, is loaded on the direction of the strong axis of the cross section. The geometry and material properties of the

Lateral Buckling of a Cantilever Beam
This example presents the lateral buckling of a cantilever beam with a narrow rectangular cross section. As shown in Figure 10, a transverse force, Fz, is loaded on the direction of the strong axis of the cross section. The geometry and material properties of the cantilever beam are: E = 10 7 kN/m 2 , G = 5 × 10 6 kN/m 2 , b = 0.1 m, h = 1.0 m, and L = 10.0 m. The theoretical buckling load is given by Timoshenko [36] as follows: where EI x is the flexural rigidity of the weak axis of the cross section, a is the torsional rigidity, and L is the length of the beam. The deformation processes of the cantilever beam with a downward concentrated load up to instability are shown in Figure 11. Figure 12 shows the load deflection curves obtained by 31 nodal points, with 62 loading increments up to three times the buckling load. The number of iterations is fewer than five within each increment. The buckling load, Pcr = 149.557 kN, obtained from Timoshenko [36], was indicated. It can be seen that the present numerical solution agrees well with the analytical solution. The theoretical buckling load is given by Timoshenko [36] as follows: where EIx is the flexural rigidity of the weak axis of the cross section, a is the torsional rigidity, and L is the length of the beam. The deformation processes of the cantilever beam with a downward concentrated load up to instability are shown in Figure 11. Figure 12 shows the load deflection curves obtained by 31 nodal points, with 62 loading increments up to three times the buckling load. The number of iterations is fewer than five within each increment. The buckling load, Pcr = 149.557 kN, obtained from Timoshenko [36], was indicated. It can be seen that the present numerical solution agrees well with the analytical solution.  The theoretical buckling load is given by Timoshenko [36] as follows: where EIx is the flexural rigidity of the weak axis of the cross section, a is the torsional rigidity, and L is the length of the beam. The deformation processes of the cantilever beam with a downward concentrated load up to instability are shown in Figure 11. Figure 12 shows the load deflection curves obtained by 31 nodal points, with 62 loading increments up to three times the buckling load. The number of iterations is fewer than five within each increment. The buckling load, Pcr = 149.557 kN, obtained from Timoshenko [36], was indicated. It can be seen that the present numerical solution agrees well with the analytical solution.

Postbuckling of a Clamped-Hinged Circular Arch
This famous example presents the geometrically nonlinear analysis of the pre-and postbuckling deformation of the circular arch, clamped at one end and hinged at the other end, with a concentrated load applied at the apex as shown in Figure 13. The material and geometric data for the arch are R = 100, EI = 10 6 , EA = GA = 100EI, 145°= f . Different numbers of nodal points were adopted to evaluate the buckling load, and 72 nodal points were capable of acquiring good results. The deformation processes of the circular arch subject to a concentrated load up to instability are plotted in Figure 14. An analytical solution by DaDeppo and Schmidt is also available in [37]. Figure 15 shows that the load factor-deflection curves of the pre-and postbuckling analysis were solved by 72 nodal points and 136 displacement increments. The number of iterations is fewer than ten within each increment. Moreover, a comparison of numerical solutions of the buckling load with FEM by Simo [2], Ibrahimbegovic [7], and Kapania [12] is given in Table 2.

Postbuckling of a Clamped-Hinged Circular Arch
This famous example presents the geometrically nonlinear analysis of the pre-and postbuckling deformation of the circular arch, clamped at one end and hinged at the other end, with a concentrated load applied at the apex as shown in Figure 13. The material and geometric data for the arch are R = 100, EI = 10 6 , EA = GA = 100EI, f = 145 • .

Postbuckling of a Clamped-Hinged Circular Arch
This famous example presents the geometrically nonlinear analysis of the pre-and postbuckling deformation of the circular arch, clamped at one end and hinged at the other end, with a concentrated load applied at the apex as shown in Figure 13. The material and geometric data for the arch are R = 100, EI = 10 6 , EA = GA = 100EI, 145°= f . Different numbers of nodal points were adopted to evaluate the buckling load, and 72 nodal points were capable of acquiring good results. The deformation processes of the circular arch subject to a concentrated load up to instability are plotted in Figure 14. An analytical solution by DaDeppo and Schmidt is also available in [37]. Figure 15 shows that the load factor-deflection curves of the pre-and postbuckling analysis were solved by 72 nodal points and 136 displacement increments. The number of iterations is fewer than ten within each increment. Moreover, a comparison of numerical solutions of the buckling load with FEM by Simo [2], Ibrahimbegovic [7], and Kapania [12] is given in Table 2.  Different numbers of nodal points were adopted to evaluate the buckling load, and 72 nodal points were capable of acquiring good results. The deformation processes of the circular arch subject to a concentrated load up to instability are plotted in Figure 14. An analytical solution by DaDeppo and Schmidt is also available in [37]. Figure 15 shows that the load factor-deflection curves of the pre-and postbuckling analysis were solved by 72 nodal points and 136 displacement increments. The number of iterations is fewer than ten within each increment. Moreover, a comparison of numerical solutions of the buckling load with FEM by Simo [2], Ibrahimbegovic [7], and Kapania [12] is given in Table 2.

Lateral Buckling of a Curve Beam
This example presents the performance of the three-dimensional nonlinear analysis of the lateral buckling of a curve beam. As shown in Figure 16, there is a simply supported curve beam with a subtended angle of 60 • . Note that the rotation with respect to the principal axes of inertia of the cross section is allowable, but it is unable to rotate about the normal direction of the cross section. The geometric and material properties of the curved beam are: curve length L = 240 mm, b = 0.6 mm, h = 30 mm, E = 71,240 N/mm 2 , and G = 27,190 N/mm 2 .
beam with a subtended angle of 60°. Note that the rotation with respect to the principal a of inertia of the cross section is allowable, but it is unable to rotate about the normal direc of the cross section. The geometric and material properties of the curved beam are: cu length L = 240 mm, b = 0.6 mm, h = 30 mm, E = 71,240 N/mm 2 , and G = 27,190 N/mm 2 .
A comparison of the effect of imperfect loading is given by several imperfect loadi FYB = α MYA, applied at central point B. The coefficients of imperfection loads, α , adopte trigger the lateral deformation, are 10 −5 and 10 −6 , respectively. Figures 17 and 18 show the ment-displacement curves associated with different imperfection loads, FYB, under posi and negative moments, in which 82 nodal points and 36 loading increments were adop The maximum number of iterations was seven during the iteration computation. The cri moments, Mcr = 835.82 N-mm and −411.70 N-mm, obtained from the curve beam theor Yang [38], are also indicated. Moreover, the postbuckling deformation processes of the cur beam from buckling moment, Mcr, to 1.04 Mcr are also shown in Figures 19 and 20.   A comparison of the effect of imperfect loading is given by several imperfect loadings, F YB = αM YA , applied at central point B. The coefficients of imperfection loads, α, adopted to trigger the lateral deformation, are 10 −5 and 10 −6 , respectively. Figures 17 and 18 show the moment-displacement curves associated with different imperfection loads, F YB , under positive and negative moments, in which 82 nodal points and 36 loading increments were adopted. The maximum number of iterations was seven during the iteration computation. The critical moments, Mcr = 835.82 N-mm and −411.70 N-mm, obtained from the curve beam theory of Yang [38], are also indicated. Moreover, the postbuckling deformation processes of the curved beam from buckling moment, Mcr, to 1.04 Mcr are also shown in Figures 19 and 20.

Lateral Buckling of a Curve Beam
This example presents the performance of the three-dimensional nonlinear analysis of the lateral buckling of a curve beam. As shown in Figure 16, there is a simply supported curve beam with a subtended angle of 60°. Note that the rotation with respect to the principal axes of inertia of the cross section is allowable, but it is unable to rotate about the normal direction of the cross section. The geometric and material properties of the curved beam are: curve length L = 240 mm, b = 0.6 mm, h = 30 mm, E = 71,240 N/mm 2 , and G = 27,190 N/mm 2 .
A comparison of the effect of imperfect loading is given by several imperfect loadings, FYB = α MYA, applied at central point B. The coefficients of imperfection loads, α , adopted to trigger the lateral deformation, are 10 −5 and 10 −6 , respectively. Figures 17 and 18 show the moment-displacement curves associated with different imperfection loads, FYB, under positive and negative moments, in which 82 nodal points and 36 loading increments were adopted. The maximum number of iterations was seven during the iteration computation. The critical moments, Mcr = 835.82 N-mm and −411.70 N-mm, obtained from the curve beam theory of Yang [38], are also indicated. Moreover, the postbuckling deformation processes of the curved beam from buckling moment, Mcr, to 1.04 Mcr are also shown in Figures 19 and 20.

Conclusions
A three-dimensional co-rotational beam formulation for the geometrically nonlinear analysis of a spatial beam was constructed in the context of a co-rotational formulation and Timoshenko beam hypothesis. In this paper, the Rodrigues formula was adopted to represent the finite rotation in a co-rotational scheme; this also circumvents the famous Gimbal lock problem that may occur with rotation in terms of the Euler angle. The numerical implementation of the DRK approximation collocation method, combined with the Newton-Raphson method, directly solved the strong forms of the geometrically nonlinear problems-different from FEM, which solves the weak form by means of the variation method. The present scheme means that the stress field can be solved directly from the governing equations instead of recovering through the differentiation procedure toward the kinematic field. That is, the primary interpolated variables will not be the dominant factor. The scheme presented is also capable of effectively solving geometrically nonlinear analyses such as large deformation, snap-through, postbuckling, and lateral buckling problems. Numerical examples are given to demonstrate the validity and effectiveness.