On a remarkable geometric-mechanical synergism based on a novel linear eigenvalue problem

The vertices of two specific eigenvectors, obtained from a novel linear eigenvalue problem, describe two curves on the surface of an N-dimensional unit hypersphere. N denotes the number of degrees of freedom in the framework of structural analysis by the Finite Element Method. The radii of curvature of these two curves are 0 and 1. They correlate with pure stretching and pure bending, respectively, of structures. The two coefficient matrices of the eigenvalue problem are the tangent stiffness matrix at the load level considered and the one at the onset of loading. The goals of this paper are to report on the numerical verification of the aforesaid geometric-mechanical synergism and to summarize current attempts of its extension to combinations of stretching and bending of structures.


Introduction
Generally, it is preferable that the loads are mainly carried by membrane and axial forces instead of bending moments. A measure to which extent this goal is reached is the percentage of the "non-membrane" energy of the total strain energy. It is defined as (U − U M )/U , where U M denotes the membrane (stretching) energy and U stands for the total strain energy. The lower bound of this ratio is zero, and it refers to pure stretching. The upper bound is one, and it refers to pure bending.
The basic idea of this work is to find a geometric quantity such that its lower and upper bound agrees with the bounds of (U − U M )/U . This is done with the help of a novel linear eigenvalue problem in the framework of the Finite Element Method (FEM). Such a quantity is the radius of the first Frenet-curvature of a curve on the surface of an N -dimensional unit hypersphere, denoted as ρ. The reason why ρ is associated with such a hypersphere is that the vertex of the unit vector, describing this curve, is an eigenvector of a novel linear eigenvalue problem referring to a structure with N degrees of freedom, discretized by the FEM. One of the two coefficient matrices of this eigenvalue problem is the tangent stiffness matrix. The other one needs to be determined such that the two aforementioned limiting values of ρ can be realized for pure stretching and pure bending, respectively.
The significance of these two limiting cases of (U − U M )/U stems from a priori knowledge of the values of this ratio. Hence, there is no need to verify the hypothetically assumed correlation between (U − U M )/U and ρ by computing U and U M separately by the FEM. If, however, an expanded form of the geometricmechanical synergism described in this work could be shown to exist, which so far has not been the case, direct computation of (U − U M )/U would be advantageous. It is mentioned, in passing, that the deformations do not directly enter into the calculation of ρ.
The first forerunner of the present work was a paper by Mang et al. [8]. The linear eigenvalue problem used for the determination of ρ was the so-called "Consistently Linearized Eigenvalue Problem" (CLE). It was originally proposed by Helnwein [3]. The two coefficient matrices of this eigenvalue problem are the tangent stiffness matrix and its derivative with respect to a dimensionless load parameter in the framework of proportional loading. A disadvantage of this model was the explicit dependence of ρ on the relevant eigenvalue, defined as the zero eigenvalue at the stability limit. This resulted in the vanishing of ρ at the stability limit. Thus, the hypothesis ρ = (U − U M )/U [8] incorrectly signaled pure stretching at this point instead of pure bending. Another weak spot of the CLE was the ab initio indefiniteness of the derivative of the tangent stiffness matrix, resulting in positive as well as negative eigenvalues in the prebuckling domain and in conjugate complex eigenvalues when the tangent stiffness matrix became indefinite after the stability limit [16], notwithstanding the insignificance of the primary load-displacement path after this has happened. A numerical challenge of the CLE was a sufficiently accurate finite-difference approximation of the derivative of the tangent stiffness matrix [4,9].
Another forerunner of the present work was a paper by Mang [6], characterized by the replacement of ρ by k a , where a denotes the acceleration of a fictitious particle moving along a curve on the surface of an N -dimensional unit hypersphere, obtained by the CLE, and where k stands for a proportionality factor. Although good numerical results were obtained for the two aforementioned limiting cases, in retrospect, lack of the invariance of k a with respect to the chosen parameter is viewed as a shortcoming of that approach.
The present work is organized as follows: In Section 2, the FEM-based linear eigenvalue problem is introduced. One of its two coefficient matrices is the tangent stiffness matrix. The other one is yet to be determined. Herein, it is considered to be a constant, symmetric, positive definite matrix. Section 3 is devoted to the numerical implementation of the theoretical concept. This includes the numerical solution of the underlying eigenvalue problem and the numerical determination of ρ with the help of finite-difference approximations of the first and the second derivative of the relevant eigenvector with respect to the parameter, to be defined in Section 2. Section 4 deals with the numerical verification of the asserted geometric-mechanical synergism for the limiting cases of pure stretching and pure bending. In Section 5, current attempts to numerically prove the existence of a wider range of validity of the geometric-mechanical synergism, delineated in this work, are described.

Linear eigenvalue problem for the determination of ρ
The mathematical formulation of the chosen eigenvalue problem reads as where K T is the tangent stiffness matrix in the framework of the FEM and B is a constant, symmetric, positive definite matrix. B must enable determination of the first eigenpair (χ 1 , r 1 ), with χ 1 denoting the smallest eigenvalue and r 1 standing for the corresponding eigenvector, normalized to 1, such that The parameter ξ represents an arc length in the context of the FEM. It depends on the vector of differential node displacements, following from the equilibrium relation with P standing for the vector of work-equivalent node forces. If snap-through can be ruled out, ξ is replaced by the dimensionless load parameter λ.
Once r 1 is known, ρ 1 can be computed from [2,5,14] In order to check whether a specific matrix B enables verification of the hypothesized correlation of U = U M with r 1 (ξ(λ)) = const., it must be shown that this matrix allows for where the subscript 0 indicates the onset of loading, i. e. λ = 0. Since the eigenvectors represent a complete basis,ṙ 1 can be expressed in terms of r j , with j = 2, 3, . . . , N : where Appendix 1 contains the derivation of (7). It would be unfeasible to use (6) for the numerical computation oḟ r 1 . Instead of doing this,ṙ 1 is approximated by a central finite-difference expression. In the following, it will be shown that setting B as where (K T ) 0 ≡ K T (ξ (λ = 0)) enables (ṙ 1 ) 0 = 0. Substitution of (8) into (1) gives Specialization of (9) for λ = 0 yields Since (K T ) 0 is a positive definite matrix, Consequently, Thus, the initial eigenvalues are an N -fold eigenvalue, equal to 1. Specialization of (7) for λ = 0 gives Differentiation of with respect to the chosen parameter yields Specialization of (15) for λ = 0 and consideration of results in Substitution of the orthogonality relations which follow from (17), into (13), and consideration of give where L'Hôpital's rule has been used. In general, c 1 j 0 = 0 ⇔ (ṙ 1 ) 0 = 0. As follows from (20), for the special case of c 1 j 0 = 0 ⇔ (ṙ 1 ) 0 = 0, the following relation is fulfilled: Differentiation of (15) with respect to the chosen parameter and specialization of the obtained equation for resulting in the orthogonality relations in addition to the orthogonality conditions following from (1), and to the orthogonality conditions (18). Computation of the second derivative of (15) with respect to the chosen parameter and specialization of the result for λ = 0, (ṙ 1 ) 0 = 0 and (r 1 ) 0 = 0 yields resulting in the orthogonality relations in addition to the three aforementioned orthogonality relations. Hence, If B had alternatively been chosen e. g. as the unit matrix I, (17) would have had to be replaced by where the symbols marked with an asterisk have replaced the corresponding symbols without an asterisk, reserved for B = (K T ) 0 . In this case, (ṙ * In Section 4, it will be shown numerically that ρ 1 = 0 correlates with pure stretching and that ρ 1 = 1 correlates with pure bending.

Convergence studies concerning spatial and "temporal" discretizations
The linear eigenvalue problem with is solved with the help of six different finite beam elements. Information about these elements is given in Appendix 2. To underline the significance of the coefficient matrix (K T ) 0 in (30), the analysis results are compared with the ones obtained with The numerical analysis involves spatial discretizations in the framework of the FEM and "temporal" discretizations resulting from finite-difference approximations ofṙ 1 andr 1 appearing in the expression for ρ 1 , see (5). In order to separate convergence studies concerning these discretizations from the numerical verification of the asserted geometric mechanical synergism, they have been moved forward to this Section. Table 1 refers to a convergence study concerning the dependence of ρ 1 on the number of finite elements used for numerical analysis of a thrust-line arch and a beam subjected to pure bending, treated in Section 4.
The results of the thrust-line arch were obtained by the Abaqus finite element B32. The deviation of ρ 1 from the hypothesized value 0 is less than 0.005709 in analyses with more than 20 elements. The beam subjected to pure bending was analyzed by the Abaqus finite elements B32OS and B32OSH. The results obtained by the former did converge, however, not to the hypothesized value 1. The results obtained by the latter are close to 1 up to 50 elements, with the absolute difference smaller than 0.0034. One reason why fewer degrees of freedom result in a value closer to 1 might be that the projection of a curve onto a unit-sphere in an infinite-dimensional frame onto a finite-dimensional frame leads to an increase of ρ 1 . In the extreme case of only two degrees of freedom, ρ 1 becomes equal to the value of 1.
As follows from (4), computation of ρ 1 requires knowledge ofṙ 1 andr 1 . These vectors are approximated by the following central finite-difference expressions: Tables 2 and 3 refer to a convergence study concerning the dependence of the first Frenet-radius ρ 1 and the second Frenet-curvature κ 2 , respectively, on the size of the load step λ. κ 2 (κ 2 is part of the expression foṙ ρ 1 , see (39)) is considerably more sensitive to the size of λ than ρ 1 because of depending also on the third derivative of r 1 , see (41). For the thrust-line arch, analyzed with the Abaqus element B32, the median value of κ 2 is numerically stable for λ ≥ 0.005. For pure bending, the median value of κ 2 for the Abaqus element B32OS is numerically stable for λ ≥ 0.002.   Table 1, has shown that the discretization with 20 finite elements provides sufficiently accurate results. Therefore, further analyses were performed with 20 finite elements. Preliminary analyses, see Table 1, have shown that this discretization provides sufficiently accurate results. Loss of stability of the arch by flexural bifurcation buckling, characterized by det (K T ) = 0 andλ > 0, occurred at p = p S ≈ 2.77 · 10 6 N/m. The existence of a stability limit is relevant Table 4 Verification (in green and bold-face) and falsification (in red and italics), respectively, of the hypothetically asserted geometric-mechanical synergism for B = (K T ) 0 and B = I to this work only insofar as the choice of the coefficient matrix B should not have a significantly larger influence on the solution for the radius of the first Frenet-curvature at, and in the vicinity of, the point on the surface curve on the N -dimensional unit hypersphere that corresponds to the stability limit than on the remaining part of this curve. Following the formulation given in Section 2, the eigenvalue problem (1) was solved for B = (K T ) 0 as well as for B = I. According to the hypothetically asserted geometric-mechanical synergism, pure stretching should correlate with ρ 1 = 0.
It is seen that for B = (K T ) 0 the hypothesis is verified for the first four out of the six finite elements considered. Conversely, for B = I, the hypothesis is verified just for the last two of these elements. The conclusion from the given numerical results that B = (K T ) 0 is superior to B = I would be premature. Nevertheless, the orthogonality of the eigenvectors with respect to B = (K T ) 0 in addition to the one with respect to K T is viewed as an advantage because of leading to an a priori known N -fold initial eigenvalue (χ i ) 0 = 1, representing a constraint on χ 1 (ξ(λ)), which is the basis for computation of ρ 1 (ξ(λ)). Figure 2 shows ρ 1 − p diagrams, obtained with B = (K T ) 0 and B = I and with the Abaqus-element B32. Apart from small deviations of ρ 1 from 0, ranging from 0.0011 to 0.17, B = (K T ) 0 provides the correct result. The upper bound, ρ 1 = 0.17, refers to a load level well above the stability limit, and out of the region of p in Fig. 2a. This confirms the main theoretical finding, reported in Section 2, that B = (K T ) 0 enables r 1 = const., which is the basis for the hypothesized geometric-mechanical synergism of ρ 1 = 0 and U = U M . Figure 3 shows a simply supported beam of length L = 5 m. The IPE 400 beam is subjected to bending moments λM y , withM y = 500 kNm, at both ends. The values of E, ν, and of the second moment of area, I y , are 210 · 10 9 Pa, 0.3, and 220 · 10 −6 m 4 , respectively.

Pure bending of a simply supported beam subjected to equal bending moments at both ends
The numerical investigation of the beam was performed with the same finite elements that were used for the analysis of the thrust-line arch. The convergence study, documented in Table 1, has shown that the  discretization with 20 finite elements provides sufficiently accurate results. Loss of stability of the beam by flexural-torsional buckling, characterized by det (K T ) = 0, occurred at M y = M y,S ≈ 286 · 10 3 N m.
The linear eigenvalue problem (1) was solved with B = (K T ) 0 as well as with B = I. According to the hypothetically asserted geometric-mechanical synergism, U M = 0 should correlate with ρ 1 = 1. Table 5 refers to verification and falsification, respectively, of this assertion for the two coefficient matrices B and the six finite beam elements described in Appendix 2.
It is seen that for B = (K T ) 0 the hypothesis is only verified for the accurate Pi&Bradford element and the Abaqus elements B32OSH and B31OSH. Since the approximate Pi&Bradford element and the Abaqus elements B32OS and B33H provided the correct results for the displacements of the beam and the von Mises stress, membrane locking was ruled out as a possible reason for the falsification of the hypothesis with these three elements. The fact that the Abaqus element B32OSH is an expansion of the Abaqus element B32OS, characterized by additional degrees of freedom to avoid volumetric locking for materials with values of Poisson's ratio larger than 0.4999999 [1, p. 357], see Appendix 2, might be the reason for its success in the given case, although the value of Poisson's ratio was chosen as 0.3. A counterargument to the assumed significance of avoiding volumetric locking is the failure of the hypothesis for the Abaqus element B32OSH in case of B = I. Interestingly, for that case, both Pi&Bradford elements and the Abaqus element B32OS provide the correct result. Figure 4 shows ρ 1 − M diagrams obtained with B = (K T ) 0 and B = I and with the Abaqus element B32OSH.
In order to demonstrate that the existence of a stability limit may affect the solution for ρ 1 in case of an inappropriate choice of the coefficient matrix B, the ρ 1 − λ diagram obtained with the CLE, i. e. with the variable, symmetric, indefinite matrix B = −K T , is shown in Fig. 5.
The reason why the hypothesis fails at the stability limit is shown hereafter. Analogous to the derivation of (7), see Appendix 1, the coefficient is obtained as [8] Because of χ 1 = 0 at λ = λ S ,ṙ instead of 1.

Numerical falsification of the hypothetically extended range of validity of the geometric-mechanical synergism
The range of the geometric-mechanical synergism was hypothetically extended to a variable ratio (U − U M )/U . It has been hypothesized that the ratio (U − U M )/U , depending on the load level, is equal to a variable radius of the first Frenet-curvature, ρ 1 . This calls for a mathematical expression for the rate of change of ρ 1 , which is given asρ Appendix 3 contains the derivation of (39). In this relation, denotes the speed of a fictitious particle, moving along the curve on the surface of the N -dimensional unit hypersphere, described by the vector r 1 (ξ(λ)), κ 2 e 3 is a vector of length κ 2 in the direction of the binormal vector e 3 , which is a unit vector; κ 2 stands for the second Frenet-curvature [2,5] of the surface curve; κ 2 e 3 is given as [2,5,14] κ 2 e 3 = κ 1 r 1 − κ 1 r 1 In (41), κ 1 = 1/ρ 1 , with κ 1 denoting the first Frenet-curvature of the surface curve, and withṡ 1 according to (40) and In the subsequent numerical investigation,ṙ 1 andr 1 are approximated by the central finite-difference expressions (33) and (34);r 1 is approximated by the central finite-difference expression ...
The two examples presented in Section 4 are characterized by The condition for an extreme value of ρ 1 , for which κ 2 = 0, is given as (46) Figure 6 shows a bar of length 5 m. The IPE 400 bar is subjected to an eccentric force λP, withP = 1 kN as the reference compressive force. The eccentricity e, E, and ν were chosen as 40.447 · 10 −3 m, 210 · 10 9 Pa and 0.3, respectively.
With the exception of (U − U M )/U at λ = 0, for which 2 nd -order theory provides the correct result, the given problem must be solved by a full nonlinear analysis, occasionally called 3 rd -order theory [7], for which  an analytic solution does not exist. The initial value of (U − U M )/U is obtained by specializing the general solution for this ratio by means of 2 nd -order theory for λ = 0. This yields where A, given as 8.0678 · 10 −3 m 2 , stands for the area of the cross-section, and I ζ , given as 13.1985 · 10 −6 m 4 , denotes the second moment of area with respect to the cross-sectional axis ζ . The value of the eccentricity of the compressive force was chosen such that (U − U M )/U at λ = 0 is equal to 0.5. This resulted in e = 4.0447 cm. Table 6 shows that the initial values of ρ 1 , obtained with B = (K T ) 0 and B = I, with four different finite elements each, do not agree with the initial value of (U − U M )/U . Irrespective of the incorrect initial values of ρ 1 , the course of the functions ρ 1 (P), where P = λP, will be discussed in the following. Figure 7 shows that ρ(P) is a non-monotonic function. This correlates with the expected non-monotony of the ratio (U − U M )/U , indicating that initially the bending energy is increasing more strongly than the membrane energy. After reaching its maximum value, ρ(P) is decreasing. This correlates with the decrease of (U − U M )/U , indicating that after reaching its maximum value the membrane energy is increasing more strongly than the bending energy. The reason for termination of the ρ − P diagram relatively soon after reaching the maximum value of ρ is the irrelevance of the assumption of a linear elastic material for large deformations. Since the value of ρ(P = 0) is significantly smaller than the initial value of (U − U M )/U , i. e. 0.5, the hypothetically extended range of validity of the geometric-mechanical synergism could not be verified. This result is corroborated by the supposition that the load level of ρ max does not agree with the one of ((U − U M )/U ) max .
The expression forρ 1 , see (39), refers to a curve on the surface of an N -dimensional unit hypersphere. The (r 1 ·e 3 )−P diagrams shown in Fig. 8 correspond to the ρ 1 − P diagrams illustrated in Fig. 7. This proves that (46), i. e. r 1 ·e 3 = 0, is the condition for an extreme value of ρ 1 .
One attempt to improve the original hypothesis was to replace ρ 1 by 1 − (r 1 ·e 3 ) 2 . The rationale for this modification was that (48) Another attempt was to explicitly consider κ 2 in the hypothesis, albeit without abandoning the restriction of the constancy of the coefficient matrix B in the linear eigenvalue problem (1). Both attempts were not successful. The second one did not rule out the possibility of values of ρ 1 larger than 1, which obviously do not correlate with corresponding values of (U − U M )/U .

Conclusions
-The proposed hypothesis was numerically verified, for some elements, for the investigated limiting cases of ρ 1 = 0 and 1. For a problem with a variable percentage of the non-membrane energy, the result for (U − U M )/U was qualitatively reasonable but quantitatively incorrect. -ρ 1 is invariant with respect to the parameter chosen for its computation as well as to the selected non-rotating Cartesian system of reference to which the element stiffness matrices are transformed before assemblage to the global stiffness matrix. -A reason why the hypothesis does not hold for all elements considered, may be the norm of the eigenvector r 1 , noting that its elements, in general, have different dimensions. This may endanger the mechanical objectivity of ρ 1 . -The rigid-body rotations of the basis of element tangent stiffness matrix correlate with the Frenet-radii of curvature of the eigenvectors. -Evidence of this situation was the numerical verification of the hypothesis for pure stretching, for four out of the six different finite elements considered, but only for the two remaining elements if (K T ) 0 was replaced by the unit matrix I. -For the limiting case of pure bending, the hypothesis was verified for three finite elements each, for (K T ) 0 and I. The only element for which the hypothesis was verified for both (K T ) 0 and I was the accurate Pi&Bradford element. One of the reasons for the falsification of the hypothesis seems to be the consequence of the inability of some finite elements to satisfy a subsidiary condition for r 1 for ρ 1 = 1 in case of (K T ) 0 and/or I. -The validity of the hypothetically asserted geometric-mechanical synergism for a variable ratio of (U − U M )/U could not be verified for the example of a bar subjected to eccentric compression. From a qualitative viewpoint, the obtained non-monotonic ρ 1 − P diagram was reasonable. It reflected the expected non-monotony of (U − U M )/U . However, the value of ρ 1 (P = 0) was significantly smaller than the one of ((U − U M )/U )(P = 0), for which an analytic solution exists. Current modifications of the original hypothesis were not successful. -The thrust of future research on the topic of this work will be the search for a mechanically objective solution of the underlying linear eigenvalue problem.
All curvatures except the top curvature, i. e. the torsion τ , are non-negative reals [2].
The radius of the first curvature, ρ 1 (s), is the inverse of the first curvature κ 1 (s):