Assessment of solutions from the consistently linearized eigenproblem by means of ﬁnite difference approximations

The consistently linearized eigenproblem (CLE) plays an important role in stability analysis of structures. Solution of the CLE requires computation of the tangent stiffness matrix e K T and of its ﬁrst derivative with respect to a dimensionless load parameter k , denoted as _ e K T . In this paper, three approaches of computation of _ e K T are discussed. They are based on (a) an analytical expression for the derivative of the element tangent stiffness matrix e K eT , (b) a load-based ﬁnite difference approximation (LBFDA), and (c) a displace-ment-based ﬁnite difference approximation (DBFDA). The convergence rate, the accuracy, and the computing time of the LBFDA and the DBFDA are compared, using the analytical solution as the benchmark result. The numerical investigation consists of the analysis of a circular arch subjected to a vertical point load at the vertex, and of a thrust-line arch under a uniformly distributed load. The main conclusion drawn from this work is that the DBFDA is superior to the LBFDA.


Introduction
Buckling is one of the most important causes of loss of the integrity of structures.Hence, the investigation of this phenomenon is very important in structural analysis and design.The socalled consistently linearized eigenproblem (CLE), originally proposed in [1,2], was initially used for ab initio estimates of stability limits by means of the Finite Element Method (FEM) [3].Herein, ab initio means ''without incremental analysis''.Compared to other modes of accompanying linear eigenvalue analysis, as proposed e.g. in [4,5], for specific problems, such as, for example the von Mises truss subjected to a point load at the vertex, the CLE provides a better estimation of the buckling load (see Fig. 2.19 in [2]).Later on, Mang and collaborators employed the CLE for assessing the initial postbuckling behavior of elastic structures [6][7][8].The CLE also plays a role in the energetical classification of limiting cases of loss of stability such as lateral torsional buckling [9].
The mathematical formulation of the CLE reads as where e K T denotes the tangent stiffness matrix of a structure in the frame of FEM, evaluated along the primary path; where k stands for a dimensionless load factor; ðk Ã The relevant eigenpair is the one that is associated with loss of stability, characterized by the semi-positive definiteness of e K T .It is given as To solve the CLE, _ e K T needs to be computed in addition to e K T .Proposing an efficient strategy for calculation of _ e K T is the main task of this work.
Three approaches for calculation of _ e K T will be discussed.The first one is based on an analytical expression for the first derivative of the tangent stiffness matrix of element e; e K e T , with respect to k for the special case of a co-rotational beam element.The second one is a finite difference approach, herein referred to as load-based finite difference approximation (LBFDA); the third one is also a finite difference approach, herein designated as displacementbased finite difference approximation (DBFDA).These two designations reflect the specific character of the two finite difference approximations.Solutions of the CLE, based on the first approach, are considered as the benchmark results.Results obtained by means of the LBFDA and the DBFDA will be compared with corresponding results obtained with the help of the analytical method.The comparison involves the convergence rate, the accuracy, and the computing time.
This paper is organized as follows: in Section 2, basic mathematical properties of the CLE will be presented.Section 3 is devoted to the derivation of an analytical expression for _ e K T for a co-rotational beam element.In Section 4, the two finite difference approaches will be delineated, including consideration of programming aspects.In Section 5, the theoretical findings will be verified in the frame of a numerical investigation.In particular, a circular arch subjected to a vertical point load at the vertex, and a thrustline arch subjected to a uniformly distributed load will be investigated with special emphasis on loss of stability.Conclusions from this investigation will be drawn in Section 6.

Basic mathematical properties of the consistently linearized eigenproblem
The The k Ã 1 -k curve refers to bifurcation buckling.S ¼ B denotes a stability limit of the form of a bifurcation point.At this point [6], Writing (1) for the first eigenpair and differentiating the obtained relation with respect to k gives (7) holds for buckling from a general stress state, and c 1j denotes the contribution of the eigenvector v Ã For buckling from a membrane stress state In Section 5.2, this special case will be verified numerically.

Analytical expression for _ e K T based on a co-rotational beam element
For a static, conservative system with N degrees of freedom (DOF), the infinitesimally incremental form of the equilibrium equations can be written as where P represents the vector of reference node forces, and with q denoting the vector of node displacements.In the frame of the FEM, e K T is obtained by assembling the element tangent stiffness matrices e K e T ; e ¼ 1; 2; . . .; M, i.e., where M denotes the number of elements, and A e is the connectivity matrix of element e.In ( 13), e K e T is a n Â n matrix, with n standing for the number of DOF of an element, and A e is a n Â N matrix.Since A e is constant, Differentiation of ( 13) with respect to k and consideration of ( 14) yields Eq. (15) shows that _ e K T can be obtained by assembling _ e K e T .Hence, differentiation of the N Â N matrix e K T is reduced to differentiation of the n Â n matrix e K e T .Because of N ) n, this results in a very significant reduction of the analytical work.
Because of its simplicity, the co-rotational approach is widely used in nonlinear finite element analysis [10][11][12].Herein, a twodimensional co-rotational beam element, based on the Euler-Bernoulli assumptions, is developed.The displacement of the element is decomposed into the rigid-body part and the deformational part, as shown in Fig. 2. The rigid-body displacement includes two parts, a translation (from 12 to 12 ) and a rotation (from 12 to 12 , described by the angle a).In the global coordinate system ðx; zÞ the displacement vector is given as and in the local coordinate system ð x; zÞ as q e ¼ u;  where u describes the change of length of the beam, which is related to the axial force N; h 1 ( h 2 ) denotes the rotation of the axis of the beam at point 1 ( 2), which results in the bending moment M 1 (M 2 ).Herein, the upper bar denotes quantities in the local coordinate system ð x; zÞ.

With the help of the principle of virtual work the matrix
In (18), X denotes the matrix for the transformation from local coordinates ð x; zÞ to global coordinates ðx; zÞ, which is given as Àsin b= l cos b= l 1 sin b= l Àcos b= l 0 Àsin b= l cos b= l 0 sin b= l Àcos b= l 1 is the length of the chord of the deformed beam and b is the angle enclosed by this chord and the x-axis, resulting in The vectors r and z were introduced for the sake of a more concise notation.They are given as T is the element tangent stiffness matrix in the local coordinate system.Its dimension is 3 Â 3. Differentiation of (18) with respect to k yields and _ q e ¼ A e Á _ q; ð27Þ with _ q denoting the vector of nodal displacement rates.
In ( 18) and (23), X; _ X; l; _ l; b; _ b; r; _ r; z, and _ z are purely geometric quantities.They are independent of the beam theory used for derivation of the expressions for N; M 1 ; M 2 , and K e T and of their derivatives with respect to k.
For the Euler-Bernoulli theory, the axial displacement u and the transverse deflection w are given as With the help of ( 28) and ( 29 where A denotes the area of the cross-section of the beam.Herein it is assumed to be constant along the axis of this beam.I y is the moment of inertia about the y-axis.The vector of internal forces Differentiation of (30) and (31) with respect to k, assuming the cross-section and the material properties to be constant, yields and Thus, all quantities appearing in (23) are known.All of them are functions of either q e or _ q e .

Finite difference approximations of
The approach described in Section 3 depends on the chosen element.Therefore, it is impractical in view of the great number of types of finite elements and the large variety of technical problems, requiring the choice of one or more problem-specific elements.
Alternatively, _ e K T can be approximated by means of the finite difference method [13,14] where h denotes a small change of k, which is positive for loading and negative for unloading.

Displacement-based finite-difference approximation of
which was already defined in (2), is redefined as a directional derivative [6], i.e.
where h is a small positive number.Following from (35), _ e K T can be approximated by a forward two-point finite-difference expression, i.e.
where K T ðq þ h _ qÞ, contrary to e K T ðqÞ, does not refer to points located on primary equilibrium paths.
Fig. 3 refers to a comparison of the two alternative finite-difference approximations of _ e K T for the special case of a system with one DOF.For this special case, (34) is replaced by where k T requires an iteration.Alternatively, (36) is replaced by In contrast to computation of k T , computation of k ð1Þ T does not require an iteration.This was one of the reasons for implementing the algorithm for determination of the DBFDA of _ e K T into the commercial FE program MSC.MARC 2012 [15].Moreover, as follows on closer inspection of Fig. 3, the right-hand side of (38) is a better T than the right-hand side of (37).Only if h is so small that the test of convergence of the iteration is satisfied already after the first iteration step, the two approaches are equivalent.
The analytical expression and the finite difference approximations of _ e K T were implemented in FEMv2, which is a Matlab-based nonlinear finite element program.The vertical displacement of the central node is chosen as the representative degree of freedom.The load-displacement paths obtained from the two FE codes are shown in Fig. 5.The result obtained from MSC.MARC 2012 agrees very well with the one obtained from FEMv2.S ¼ B denotes a bifurcation point which, in the given case, is the stability limit.Hence, the snap-through point D has no physical significance.

Circular arch
The quality of the DBFDA and LBFDA of _ e K T is assessed by comparing the dependence of a suitable error norm of _ e K T (see Fig. 6) and of the error of k Ã 1 (see Fig. 7) on h.In Fig. 6, where _ e K T EX indicates calculation of _ e K T from the analytical expression derived in Section 2, and _ e K T DBFDAðLBFDAÞ refers to its calculation as a DBFDA and LBFDA, respectively.The norm of the two matrices in (39) is defined as follows: where c ij is the coefficient in the ith row and the jth column of C .The red (black) curve illustrates the dependence of the error norm of _ e K T , based on the DBFDA (LBFDA) of this matrix, on h.For logðhÞ > $ À15, the error of the LBDFA is larger than the one of the DBFDA, which corroborates a statement in connection with the explanation of Fig. 3.
The convergence rate of a numerical approach for calculation of a quantity is directly related to the slope of the curve in a log-log plot [16].The red line in Fig. 6 is parallel to the black line, indicating that the convergence rate of the error norm based on the q 0 0 1 DBFDA of _ e K T is the same as the one of the error norm based on the LBFDA of this matrix.For logðhÞ < $ À14, the error is increasing, which is the consequence of the limitation of the number of digits in the representation of numbers in the computing machine.Hence, the lowest point of each curve characterizes the accuracy of the respective FDA.According to Fig. 6, the accuracy of the DBFDA is higher than the one of the LBFDA.
The second quantity used for assessing the accuracy of the two finite difference approximations is k Ã 1 .The results of this assessment are shown in Fig. 7. Herein, where k Ã 1 EX indicates calculation of k Ã 1 from the analytical solution, and k Ã 1 DBFDAðLBFDAÞ refers to calculation of k Ã 1 from the DBFDA and the LBFDA, respectively, of k Ã 1 .For logðhÞ > $ À13, the error of the DBDFA is larger than the one of the LBFDA.This can be explained by investigating the special case of a system with one DOF (see Fig. 3).Solving the CLE for this special case yields respectively, where _ k - T , is larger than that based on the LBFDA.
Nevertheless, the DBFDA of k Ã 1 has the same convergence rate as the LBFDA.For the same error tolerance within the accuracy range of each approximation, t DBFDA < t LBFDA , where t stands for the computer time.E.g., for s tol ¼ 10 À10 ; t DBFDA ¼ 16s < 72s ¼ t LBFDA .Hence, concerning computing time, the DBFDA of k Ã 1 is superior to the LBFDA of this quantity, which overcompensate the slightly larger error of the former for relevant values of h.In Fig. 9, the product of v Ã 1 ðk ¼ 0Þ Á v Ã 1 ðkÞ, which is related to the angle enclosed by the vectors v Ã 1 ðkÞ and v Ã 1 ðk ¼ 0Þ, is plotted.The curve has a minimum at B, which indicates that the prebuckling state involves both membrane and bending stresses.

Thrust-line arch
Fig. 5 shows a two-hinged arch subjected to a vertical uniform load.The geometric form of the axis of the arch is given as Geometric parameters and material data are shown in Fig. 10.Contrary to a three-hinged arch of the same geometric configuration, subjected to the same load as the two-hinged arch shown in Fig. 10, the latter, strictly speaking, is no thrust-line arch.However, since its bending moments are negligibly small [14], it is justified to speak of a thrust-line arch.For such an arch, buckling occurs from a membrane stress state.
Herein, the structure is analyzed by means of FEMv2 as well as MSC.MARC 2012.The arch is discretized by 100 two-node beam elements.The vertical displacement of the central node is chosen as the representative degree of freedom.The load-displacement paths obtained from FEMv2 and MSC.MARC 2012 are shown in Fig. 11.The results obtained from MSC.MARC 2012 agree very well with the ones obtained from FEMv2.Fig. 11(a) shows the entire computed load-displacement curve, containing a bifurcation point S ¼ B, followed by a snap-through point, D. Hence, the latter is physically insignificant.

Conclusions
Motivated by visualization of the concept of energy-based categorization of buckling problems by means of spherical geometry, methods for numerical solutions of the CLE, representing the mathematical tool for this categorization, were proposed and evaluated in this paper.
The two-dimensional co-rotational beam element, used in FEMv2, was found to be suitable for static structural stability analysis.The structural response agrees very well with the one obtained by means of the commercial finite element program MSC.MARC 2012.
The DBFDA of _ e K T is more accurate than the LBFDA of this matrix.
The convergence rate of the DBFDA of _ e K T and of the numerical solution for k Ã 1 based on this approximation is the same as that of the LBFDA of _ e K T and of the numerical solution for k Ã 1 based on this approximation.
The numerical solution for k Ã 1 based on the DBFDA of _ e K T is less accurate than the one based on the LBFDA of this matrix.However, for the same error tolerance, the former is superior to the latter as regards computing time.
The eigenvector v Ã 1 for the circular arch subjected to a vertical point load at the vertex changes its direction in the prebuckling regime, which reflects the existence of bending.The angle enclosed by v Ã 1 ðkÞ and v Ã 1 ðk ¼ 0Þ becomes a maximum at the stability limit.
The eigenvector v Ã 1 for the thrust-line arch subjected to a vertical uniform load is constant, which indicates that there is no bending before buckling.
The present work paves the way for the implementation of a stable and accurate numerical approach for the solution of the CLE in MSC.MARC 2012, which is a necessary prerequisite for the numerical realization of a new concept of categorization of buckling problems with the help of spherical geometry.

T and k ð2Þ T
denote the slopes of the tangents to the curve kðqÞ at point 0 and point 2, respectively, and h 6 Dk ð0Þ .Computation of k

Fig. 4
Fig. 4 shows a circular arch subjected to a vertical point load P at the vertex.The geometric properties and the material parameters are also shown in this figure.The arch is discretized by 100 two-node beam elements.It is analyzed by means of FEMv2 as well as by MSC.MARC 2012.The vertical displacement of the central node is chosen as the representative degree of freedom.The load-displacement paths obtained from the two FE codes are shown in Fig.5.The result obtained from MSC.MARC 2012 agrees very well with the one obtained from FEMv2.S ¼ B denotes a bifurcation point which, in the given case, is the stability limit.Hence, the snap-through point D has no physical significance.

Fig. 4 .
Fig. 4. Configuration and material properties of a circular arch subjected to a vertical point load.

Fig. 6 .
Fig. 6. logðhÞ-logðsÞ diagram for error assessment of the norm of _ e K T , based on the DBFDA and the LBFDA, respectively, of this matrix.

Fig. 8 (
a) shows the k Ã 1 -k and the k Ã 2 -k diagram, related to the first two eigenvalues of the CLE.The k Ã 1 -k curve has a horizontal tangent at the bifurcation point B (see Fig. 8(b)), whereas the k Ã 2 -k curve has a cusp of second kind, with a tangent normal to the diagonal, i.e. to the dashed line in Fig. 8(a) and (c).
Fig. 11(b)  illustrates the initial part of the load-displacement path.

Fig. 12
Fig.12shows a plot of the function v Ã 1 ðk ¼ 0Þ Á v Ã 1 ðkÞ for the thrust-line arch subjected to a uniformly distributed load.The straight line indicates that v Ã 1 is constant in the prebuckling regime, which confirms a remarkable mathematical property of the CLE, meaning that v Ã 1 is constant for a membrane stress state in the prebuckling regime.

Fig. 10 .Fig. 11 .Fig. 12
Fig. 10.Configuration and material properties of a two-hinged arch subjected to a vertical uniform load.
. Based on two different definitions of the first derivative of e K T with respect to k, two alternative finite- T is approximated by a forward two-point finite-difference expression, i.e._ e K T % e K T ðk þ hÞ À e K TðkÞh ;