Accurate Timoshenko Beam Elements For Linear Elastostatics and LPB Stability

Abstract

Several methods to derive accurate Timoshenko beam finite elements are presented and compared. Two application problems are examined: linear elastostatics and linearized prebuckling (LPB) stability analysis. Accurate elements can be derived for both problems using a well known technique that long preceeds the Finite Element Method: using homogeneous solutions of the governing equations as shape functions. An interesting question is: can accurate elements be derived with simpler assumptions? In particular, can linear-linear interpolation of displacements and rotations with one-point integration reproduce those elements? The answers are: no if standard variational tools based on classical functionals are used, but yes if modified functionals are introduced. The connection of modified functionals to newer methods, in particular templates, modified differential equations and Finite Increment Calculus (FIC) are examined. The results brings closure to a 50-year conumdrum centered on this particular finite element model. In addition, the discovery of modified functionals provides motivation for extending these methods to full geometrically nonlinear analysis while still using inexpensive numerical integration.

Appendices

Appendix 1: Abbreviations and Nomenclature

This Appendix contains tables pertaining to abbreviations and symbols used in paper (Tables 3, 4, 5, 6, 7).

Table 3 Abbreviations used in paper

Table 4 Upper case roman and italics symbols

Table 5 Notation: lower case roman, italics and caligraphic symbols

Table 6 Notation: bold symbols

Table 7 Notation: greek, caligraphic and special symbols

Appendix 2: Legendre Polynomials

In this Appendix, which focuses on individual elements, superscript e is omitted to reduce clutter.

1.1 2.1 Legendre Polynomial Basis

If the nonlinear effect of the axial force is ignored, the homogeneous solutions of the equilibrium equations listed in Sect. 2.4 are standard cubic polynomials. These have been extensively used in the HODES method since the early 1920s; cf. Appendix Sect. 8.2. These functions and their derivatives possess limited element-integrated orthogonality properties. As pointed out in [30, 31] those properties can be improved by passing to the Legendre polynomial basis (Legendre, 1785).

The first four Legendre polynomials, expressed in the natural isoparametric coordinate $\xi =2x/\ell$, are used in this and the next Appendix, except that the cubic one is modified by introducing a free parameter $\chi$ in its linear term. (This modification helps to derive optimal geometric stiffness matrices.) In this basis the transverse displacement assumption over the element is written

$$\begin{aligned} v(\xi )=c_{0} P_{0} +c_{1} P_{1}+c_{2} P_{2}+c_{3} P_{(\chi )}= \mathbf {P}\mathbf {c}, \end{aligned}$$

(67)

in which $\xi$ is the isoparametric coordinate, $\mathbf {c}=\left[ \begin{array}{cccc}c_{0}&c_{1}&c_{2}&c_{3}\end{array}\right] ^T$ and $\mathbf {P}=\left[ \begin{array}{cccc}P_{0}&P_{1}&P_{2}&P_{(\chi )}\end{array}\right]$, with

$$\begin{aligned} P_{0}=1, \quad P_{1}=\xi , \quad P_{2}={\textstyle {1\over 2}}(3\xi ^2-1), \quad P_{(\chi )}={\textstyle {1\over 2}}(5 \xi ^3-\chi \xi ). \end{aligned}$$

(68)

This will be called the parametrized Legendre basis. $P_{(\chi )}$ becomes the conventional Legendre cubic polynomial $P_3={\textstyle {1\over 2}}(5\xi ^3-3\xi )$ if $\chi =3$. Figure 8a illustrates the classical $P_{0}$ through $P_{3}$; note that $P_{n}(1)=1$, $P_{n}(-1)=(-1)^n$ for $n=0,1\ldots$. Figure 8b shows $P_{(\chi )}$ for various values of $\chi$.

For beam elements, each function in (68) has clear physical meaning in the spirit of natural modes (also called deformation patterns) advocated by Argyris since the early 1960s [4, 5]:

• $P_{0}$ and $P_{1}$ represent translational and rotational rigid motions, respectively.

• $P_{2}$ is a constant-curvature bending mode, which is symmetric about $\xi =0$. For the BE model it represents pure bending with zero shear.

• $P_{(\chi )}$ is a linear-curvature bending mode, which is antisymmetric about $\xi =0$. For the Timoshenko model it represents a linearly varying moment with constant mean shear strain.

For a deeper study of mathematical properties of Legendre polynomials, as well as their generation by recurrent relations, see [1, Ch. 22]. An important application in numerical analysis is the construction of 1D Gauss integration rules: the abscissas of the n-point rule are the zeros of $P_{n}(\xi )$.

1.2 2.2 Orthogonality Properties and Covariance Matrices

Orthogonality properties of (67) can be exhibited by introducing the row matrix derivatives

$$\begin{aligned} \mathbf {P}'= \left[ \begin{array}{cccc} P_{0}'&P_{1}'&P_{2}'&P_{(\chi )}' \end{array}\right] , \quad \mathbf {P}''= \left[ \begin{array}{cccc} P_{0}''&P_{1}''&P_{2}''&P_{(\chi )}'' \end{array}\right] , \quad \mathbf {P}'''=\left[ \begin{array}{cccc} P_{0}'''&P_{1}'''&P_{2}'''&P_{(\chi )}''' \end{array}\right] , \end{aligned}$$

(69)

where primes denote derivatives with respect to x. Since the parameter $\chi$ affects only the linear $\xi$ term in $P_{(\chi )}$, $\mathbf {P}''$ and $\mathbf {P}'''$ are independent of $\chi$. Four covariance matrices, denoted by $\mathbf {V}$ here, are obtained by taking dyadic product integrals over the element length $\ell$:

$$\begin{aligned} \mathbf {V}_0=\int _{-\ell /2}^{\ell /2} \mathbf {P}^T \mathbf {P}dx, \mathbf {V}_1=\int _{-\ell /2}^{\ell /2} (\mathbf {P}')^T \mathbf {P}' dx, \mathbf {V}_2=\int _{-\ell /2}^{\ell /2} (\mathbf {P}'')^T \mathbf {P}'' dx, \mathbf {V}_3=\int _{-\ell /2}^{\ell /2} (\mathbf {P}''')^T \mathbf {P}''' dx. \end{aligned}$$

(70)

Fig. 8

The Legendre basis over $\xi \in [-1,1]$. a classical Legendre polynomials $P_0$ through $P_3$; b cubic $P_{(\chi )}$ drawn for $\chi =0,3,5,8$—note that the vertical scale is different

Recalling that $dx={\textstyle {1\over 2}}\ell d\xi$, evaluation of (70) gives

$$\begin{aligned} \mathbf {V}_0&= \ell \left[ \begin{array}{cccc} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1/3 &{} 0 &{} (3-\chi )/6 \\ 0 &{} 0 &{} 1/5 &{} 0 \\ 0 &{} (3-\chi )/6 &{} 0 &{} (75-42\chi +7\chi ^2)/84 \end{array}\right] , \nonumber \\ \mathbf {V}_1&= {4\over \ell }\left[ \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} (5-\chi )/2 \\ 0 &{} 0 &{} 3 &{} 0 \\ 0 &{} (5-\chi )/2 &{} 0 &{} (45-10\chi +\chi ^2)/4\end{array}\right] , \nonumber \\ \mathbf {V}_2&= {48\over \ell ^3}\left[ \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 3 &{} 0 \\ 0 &{} 0 &{} 0 &{} 25 \end{array}\right] =\mathbf{diag }[0, 0, 144, 1200]/\ell ^3, \nonumber \\ \mathbf {V}_3&= {14400\over \ell ^5}\left[ \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 \end{array}\right] = \mathbf{diag }[0, 0, 0, 14400]/\ell ^5, \end{aligned}$$

(71)

It can be observed that

• $\mathbf {V}_2$ and $\mathbf {V}_3$ are diagonal and independent of $\chi$, hence displaying full orthogonality of entries of $\mathbf {P}''$ and $\mathbf {P}'''$ (the last one being trivial)

• Setting $\chi =3$ (whence $P_{(\chi )}$ becomes the classical Legendre cubic polynomial $P_3$) diagonalizes $\mathbf {V}_0=\ell \mathbf{diag }[ 1, 1/3, 1/5, 1/7]$.

• Setting $\chi =5$ diagonalizes $\mathbf {V}_1= (4/\ell ) \mathbf{diag }[ 0, 1, 3, 5]$.

• No value of $\chi$ makes all four covariant matrices diagonal.

For future use, $\mathbf {V}_2$ is split as $\mathbf {V}_2=\mathbf {V}_{2b}+\mathbf {V}_{2h}$, where

$$\begin{aligned} \mathbf {V}_{2b}=\mathbf{diag }[0, 0, 144, 0]/\ell ^3, \quad \mathbf {V}_{2h}=\mathbf{diag }[0, 0, 0, 1200]/\ell ^3. \end{aligned}$$

(72)

The subscripts b and h in (72) denote “basic” and “higher order,” respectively, in accordance with the terminology of Bergan’s Free Formulation [10, 11]. Note that $\mathbf {V}_{2b}$ and $\mathbf {V}_{2h}$ are associated with the symmetric and antisymmetric Legendre bending modes, respectively. All of the foregoing results are independent of $\Phi$. Consequently, they are applicable to both Timoshenko and BE models.

1.3 2.3 Material Stiffness Matrix

The material stiffness matrix will be derived in two steps: first as a natural stiffness in terms of natural coordinates (the Legendre polynomial coefficients), and then transformed to physical coordinates (the node displacements). This staged process allows extra flexibility in several derivations, notably those in Appendix Sect. 2.6. It is extended to templates in Appendix 3.

1.3.1 2.3.1 Natural Material Stiffness

A beam material stiffness matrix expressed in terms of the Legendre coefficients $c_i$ defined in (67) will be called a natural material stiffness, and denoted by $\mathbf {S}_M$. [The term “natural” is explained in Appendix Sect. 2.1.] Since $\mathbf {S}_M$ depends only on $\mathbf {V}''$ and $\mathbf {V}'''$, it is independent of $\chi$.

Matrix $\mathbf {S}_M$ is the Hessian of the natural strain energy denoted by $\Upsilon =\Upsilon _f+\Upsilon _s$, where $\Upsilon _f$ comes from flexure (bending) and $\Upsilon _s$ from shear. The former is further split as $\Upsilon _f=\Upsilon _b+\Upsilon _h$, due to symmetric bending $v''_b=\mathbf {V}''_b \mathbf {c}$, (basic bending mode) and antisymmetric bending $v''_h=\mathbf {V}''_h \mathbf {c}$ (higher order bending mode). Denote by $R_{b}$, $R_{h}$ and $R_{s}$ the rigidities pertaining to $\Upsilon _b$, $\Upsilon _h$, and $\Upsilon _s$, respectively, which are assumed constant along the element. Then

$$\begin{aligned} \Upsilon _b&={\textstyle {1\over 2}}\int _{-\ell /2}^{\ell /2} v''_b R_{b} v''_b dx= {\textstyle {1\over 2}}\mathbf {c}^T R_b\int _{-\ell /2}^{\ell /2} (\mathbf {P}'')^T \mathbf {P}'' dx \mathbf {c}= {\textstyle {1\over 2}}\mathbf {c}^T R_b \mathbf {V}_{2b} \mathbf {c}, \nonumber \\ \Upsilon _h&={\textstyle {1\over 2}}\int _{-\ell /2}^{\ell /2} v''_h R_{h} v''_h dx= {\textstyle {1\over 2}}\mathbf {c}^T R_h\int _{-\ell /2}^{\ell /2} (\mathbf {P}'')^T \mathbf {P}'' dx \mathbf {c}= {\textstyle {1\over 2}}\mathbf {c}^T R_h \mathbf {V}_{2h} \mathbf {c}, \nonumber \\ \Upsilon _s&={\textstyle {1\over 2}}\int _{-\ell /2}^{\ell /2} v''' R_s v''' dx= {\textstyle {1\over 2}}\mathbf {c}^T R_s\int _{-\ell /2}^{\ell /2} (\mathbf {P}''')^T \mathbf {P}''' dx \mathbf {c}= {\textstyle {1\over 2}}\mathbf {c}^T R_s \mathbf {V}_3 \mathbf {c}. \end{aligned}$$

(73)

The natural material stiffness matrix accordingly splits as

$$\begin{aligned} \mathbf {S}_{M}=\mathbf {S}_{Mb}+\mathbf {S}_{Mh}+\mathbf {S}_{Ms}, \end{aligned}$$

(74)

The components of $\mathbf {S}_M$ are obtained by taking the Hessians

$$\begin{aligned} \mathbf {S}_{Mb}= {\partial ^2 \Upsilon _b\over \partial \mathbf {c} \,\, \partial \mathbf {c}}=R_b \mathbf {V}''_b, \quad \mathbf {S}_{Mh}= {\partial ^2 \Upsilon _h\over \partial \mathbf {c}\,\, \partial \mathbf {c}}=R_h \mathbf {V}''_h, \quad \mathbf {S}_{Ms}= {\partial ^2 \Upsilon _s\over \partial \mathbf {c} \,\, \partial \mathbf {c}}=R_s \mathbf {V}''_s, \end{aligned}$$

(75)

All these matrices are diagonal, independent of $\chi$, and have rank at most one. Normally $R_b=R_h=EI$ but sometimes it is useful to adjust these values, as done to construct hinge-equipped elements later. $R_s$ is a shear rigidity coefficient defined as

$$\begin{aligned} R_s={EI\over GA_s} R_h= {\textstyle {1\over 12}}R_h \Phi \ell ^2, \quad \hbox { where } \Phi =12EI/(GA_s\ell ^2). \end{aligned}$$

(76)

$R_s$ vanishes for the BE model.

1.3.2 2.3.2 Physical Material Stiffness

The physical and natural DOF are linked by $\mathbf {u}=\mathbf {G}\mathbf {c}$ and $\mathbf {c}=\mathbf {G}^{-1} \mathbf {u}=\mathbf {H}\mathbf {u}$, where $\mathbf {u}=\left[ \begin{array}{cccc}v_1&\theta _1&v_2&\theta _2\end{array}\right] ^T$ and $\mathbf {c}=\left[ \begin{array}{cccc}c_0&c_1&c_2&c_3\end{array}\right] ^T$. Both $\mathbf {G}$ and $\mathbf {H}$ are $4\times 4$ matrices. For the Timoshenko model the total section rotation is $\theta =v'-\gamma =v'+EI/(GA_s) v'''=v'+\Phi \ell ^2 v'''/12$. Evaluating at nodes and inverting gives

$$\begin{aligned} \mathbf {G}=\left[ \begin{array}{cccc} 1 &{} -\,1 &{} 1 &{} G_{14} \\ 0 &{} 2/\ell &{} -\,6/\ell &{} G_{24} \\ 1 &{} 1 &{} 1 &{} G_{34} \\ 0 &{} 2/\ell &{} 6/\ell &{} G_{44} \end{array}\right] , \quad \mathbf {H}=\mathbf {G}^{-1}={1\over 60} \left[ \begin{array}{cccc} 30 &{} 5\ell &{} 30 &{} -\,5\ell \\ H_{21} &{} H_{22} &{} H_{23} &{} H_{24} \\ 0 &{} -\,5\ell &{} 0 &{} 5\ell \\ H_{41} &{} H_{42} &{} H_{43} &{} H_{44} \end{array}\right] , \end{aligned}$$

(77)

where $G_{14}=( \chi -5 )/2$, $G_{34}=-G_{14}$, $G_{24}=G_{44}=(15-\chi +10\Phi )/\ell$, $H_{21}= -30+5(\chi -5)/(1+\Phi )$, $H_{23}=-H_{21}$, $H_{22}=H_{24}=3\ell (\chi -5)/(2(1+\Phi ))$, $H_{41}=6/(1+\Phi )$, $H_{43}=-H_{41}$, and $H_{42}=H_{44}=3\ell /(1+\Phi )$. For the BE model set $\Phi =0$.

The physical material matrix is obtained by applying the $\mathbf {H}$ transformation matrix (77) to (74). We assume that all rigidities are constant along the element. Thus

$$\begin{aligned} \mathbf {K}_M=\mathbf {H}_c^T (\mathbf {S}_{Mb}+\mathbf {S}_{Mh}+\mathbf {S}^g_{Ms}) \mathbf {H}_c. =\mathbf {K}_{Mb}+\mathbf {K}_{Mh}+\mathbf {K}_{Ms}, \end{aligned}$$

(78)

all of which are independent of $\chi$. It is instructive to express each component as a rank-one matrix in dyadic (spectral) form. Define the basic (symmetric) and higher-order (antisymmetric) eigenvectors

$$\begin{aligned} \mathbf {v}_b=\left[ \begin{array}{cccc}-1&0&1&0\end{array}\right] ^T, \quad \mathbf {v}_h=\left[ \begin{array}{cccc} 2/\ell&1&-\,2/\ell&1 \end{array}\right] ^T. \end{aligned}$$

(79)

(These can be normalized to unit length by dividing through by $\sqrt{2}$ and $\sqrt{2}\sqrt{4+\ell ^2}/\ell$, respectively.) Then

$$\begin{aligned} \mathbf {K}_{Mb}&=\mathbf {H}_c^T \mathbf {S}_{Mb} \mathbf {H}_c=C_b \mathbf {v}_b \mathbf {v}_b^T =C_b \left[ \begin{array}{cccc}0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} -\,1 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} -\,1 &{} 0 &{} 1 \end{array}\right] , \nonumber \\ \mathbf {K}_{Mh}&=\mathbf {H}_c^T \mathbf {S}_{Mb} \mathbf {H}_c=C_h \mathbf {v}_h \mathbf {v}_h^T =C_h \left[ \begin{array}{cccc} 4/\ell ^2 &{} 2/\ell &{} -\,4/\ell ^2 &{} 2/\ell \\ 2/\ell &{} 1 &{} -\,2/\ell &{} 1 \\ -4/\ell ^2 &{} -\,2\ell &{} 4/\ell ^2/ &{} 2\ell \\ -2/\ell &{} 1 &{} 2/\ell &{} 1 \end{array}\right] , \nonumber \\ \mathbf {K}_{Ms}&=\mathbf {H}_c^T \mathbf {S}_{Ms} \mathbf {H}_c=C_s \mathbf {v}_h \mathbf {v}_h^T= \mathbf {K}_{Mh} \Phi . \hbox { (a new result)}. \end{aligned}$$

(80)

For a prismatic element $R_b=R_h=EI$. If so $C_b=2EI/\ell$, $C_h= 3 EI (4+\ell ^2)/(\ell (1+\Phi )^2)$, and $C_s=C_h \Phi$. Adding these matrices: $\mathbf {K}_M=C_b \mathbf {v}_b \mathbf {v}_b^T+C_h (1+\Phi ) \mathbf {v}_h\mathbf {v}_h^T$ yields the optimal, nodally-exact, material stiffness (29). Since the eigenvectors (79) are orthogonal, they survive in $\mathbf {K}_M$, which has rank 2, with two nonzero eigenvalues: $C_b$ and $C_h (1+\Phi )$.

The template interpretation is: $\mathbf {K}_{Mb}$ is the basic stiffness, whereas $\mathbf {K}_{Mh}+\mathbf {K}_{Ms}=\mathbf {K}_{Mh} (1+\Phi )$ is the higher order stiffness. This was first noted in [27] for the BE model, in which $\mathbf {K}_{Ms}$ vanishes. Further development of templates is done in Appendix 3.

1.4 2.4 Geometric Stiffness Matrix

The geometric stiffness matrix can be derived from the $W^e_v$ work functional given in (14)– (15), where it is also called the geometric energy functional. By analogy with the natural strain energy notation in (73), evaluation of this term in the parametrized Legendre basis (68) will be denoted by $\Upsilon _F$. We have

$$\begin{aligned} \Upsilon _F= {\textstyle {1\over 2}}\int _{-\ell /2}^{\ell /2} F (v')^2 dx = {\textstyle {1\over 2}}F \int _{-\ell /2}^{\ell /2} \mathbf {c}^T (\mathbf {P}')^T \mathbf {P}' \mathbf {c}dx= {\textstyle {1\over 2}}F \mathbf {c}^T \mathbf {V}_1 \mathbf {c}. \end{aligned}$$

(81)

The natural geometric stiffness is the Hessian of $\Upsilon _F$ with respect to $\mathbf {c}$:

$$\begin{aligned} \mathbf {S}_{G}= {\partial ^2 \Upsilon _f\over \partial \mathbf {c}\partial \mathbf {c}}=F \mathbf {V}_1= {4F\over \ell }\left[ \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} (5-\chi )/2 \\ 0 &{} 0 &{} 3 &{} 0 \\ 0 &{} (5-\chi )/2 &{} 0 &{} (45-10\chi +\chi ^2)/4\end{array}\right] , \end{aligned}$$

(82)

Transformation to physical coordinates gives $\mathbf {K}_{Gv}= \mathbf {H}^T \mathbf {S}_{G} \mathbf {H}$, which is independent of $\chi$. This matrix agrees with the consistent geometric stiffness displayed in (36), derived using the shape functions of Sect. 2.7. Setting $\Phi =0$ yields the BE version (38). The latter can be traced back to [58]. It was noted in Sect. 4.1 that those matrices are NOT optimal. Reaching that goal requires more advanced mathematical machinery, which is deployed in Appendix 3.

1.5 2.5 Hinged Timoshenko Beam Elements

This section derives stiffness matrices for Timoshenko beam elements with an internal hinge. Some of the following results have been used in graduate courses since 1989 but not published.

1.6 2.6 Midhinged Beam Element

The Timoshenko beam element shown in Fig. 9a has a mechanical hinge at midspan ($\xi =0$). Cross sections on both hinge sides can freely rotate respect to each other. Figure 9b sketches in 3D view a fabrication method sometimes used in short-span pedestrian bridges: gaps on either side of the hinged section cuts are filled with a bituminous material that permits small relative rotations.

Fig. 9

Beam element with a midspan hinge. a 2D element idealization; b possible fabrication method in 3D view

This configuration can be treated with two ordinary elements of length $\ell /2$ with a special 3-DOF node at the hinge, where lateral displacements are identical but rotations are not. This can introduce complications in programming. Those can be avoided by developing a single two-noded hinged element, which may be useful for articulated members that appear in robotics, as well as modeling midspan damage in reinforced concrete beams.

Both the curvature $\kappa =v''$ and the bending moment $M=EI \kappa$ must vanish at midspan. But in a element built via cubic interpolation of v(x), $\kappa$ must vary linearly in both $\xi$ and x. Consequently the mean curvature, which is controlled by the Legendre polynomial $P_{2}$ (that shown in blue in Fig. 8) must be zero. That constraint can be enforced by setting the symmetric (basic) bending rigidity $R_{b}$ to zero whereas the antisymmetric bending rigidity is the normal one: $R_{h}=EI$. Plugging into (78) yields the nodally-exact material stiffness matrix

$$\begin{aligned} \mathbf {K}_{MH}=\mathbf {K}_{Mh} (1+\Phi )=C_h (1+\Phi ) \mathbf {v}_h^{} \mathbf {v}_h^T= {EI\over \ell ^3(1+\Phi )} \left[ \begin{array}{cccc} 12 &{} 6 \ell &{} -\,12 &{} 6\ell \\ 6\ell &{} 3\ell ^2 &{} -\,6\ell &{} \ell ^2 \\ -12 &{} -\,6 \ell &{} 12 &{} -\,6 \ell \\ 6\ell &{} \ell ^2 &{} -\,6\ell &{} 3\ell ^2 \end{array}\right] . \end{aligned}$$

(83)

Plainly $\mathbf {K}_{MH}^e$ has rank one; its nonzero eigenvalue is $6 EI (4+\ell ^2) (1+\Phi )/\ell ^3$. The result may be gotten by more advanced methods, such as mixed variational principles (assuming transverse displacements and bending moments), but the above derivation uses only undergraduate mathematics.

Fig. 10

Two generalizations of the element of Fig. 9. a spring-midhinged element; b hinge moved to an arbitrary location $\xi =\xi _H$ within element

1.7 2.7 Spring-MidHinged Beam Element

A generalization of the previous element is shown in Fig. 10a: the midpoint hinge is restrained by a torsional (a.k.a. rotational) spring and an extensional spring with spring constants $k_T^{}=\mu EI/\ell$ and $k_E^{}=\delta EI/\ell ^3$, respectively. The rigidity parameters $\mu$ and $\delta$ are dimensionless.

The stiffness derivation is more involved. The transverse displacement is piecewise split into $v^-(\xi )$ going from $\xi =-1$ to $\xi =0$ (left halfspan), and $v^+(\xi )$ going from $\xi =0$ to $\xi =1$ (right halfspan). See Fig. 11 for visualization of this kinematic split. Both portions are expressed in terms of Legendre polynomials, giving 8 natural coefficients. Eight symbolic equations are then established for kinematic and static matching at node 1 ($\xi =-1$), node 2 ($\xi =1$) and the hinge node ($\xi =0$). Solving these equations via Mathematica links the natural and physical coordinates.

The resulting material stiffness matrix can be presented as a sum of two contributions:

$$\begin{aligned} \mathbf {K}_{MHk}= \mathbf {K}_{MH}+\mathbf {K}_{Mk}, \end{aligned}$$

(84)

where $\mathbf {K}_{MH}$ is (83) for the free hinge while $\mathbf {K}_{Mk}$ brings in the effect of the springs:

$$\begin{aligned} \mathbf {K}_{Mk}= \left[ \begin{array}{cccc}K_{Mk11} &{} K_{Mk12} &{} -\,K_{Mk11} &{} K_{Mk12} \\ K_{Mk12} &{} K_{Mk22} &{} -\,K_{Mk12} &{} K_{Mk22} \\ -K_{Mk11} &{} -\,K_{Mk12} &{} -\,K_{Mk11} &{} -\,K_{Mk12} \\ K_{Mk12} &{} K_{Mk12} &{} -\,K_{Mk12} &{} K_{Mk12} \end{array}\right] , \end{aligned}$$

(85)

in which

$$\begin{aligned} K_{Mk11}&= 12 C_k \delta (1+\mu ) (1+\Phi ), \quad K_{Mk12}= 3 C_k \ell \delta (2+\mu ) (1+\Phi ), \nonumber \\ K_{Mk22}&= C_k \ell ^2 (1+\Phi ) \big (48 \mu +\delta (3+\mu +\mu \Phi )\big ), \nonumber \\ C_k&={4 EI\over \ell ^3 (1+\Phi ) \big (192 (1+\mu ) +\delta (4+\mu +4 (1+\mu ) \Phi )\big )}. \end{aligned}$$

(86)

This matrix simplifies if one spring is missing. For example, if $\delta =0$ so $k_E=0$, it reduces to

$$\begin{aligned} \mathbf {K}_{Mk}= {EI \mu \over (1+\mu ) \ell } \left[ \begin{array}{cccc}0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} -\,1 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} -\,1 &{} 0 &{} 1 \end{array}\right] ,\end{aligned}$$

(87)

which is independent of $\Phi$. $\mathbf {K}_{Mk}$ has rank equal to the number of springs present. The rank of $\mathbf {K}_{MHk}$ is then 3 or 2. One rigid body mode: rigid rotation about the hinge, always remains.

Some limits can be checked to verify the above results: $\delta =0, \mu \rightarrow \infty$: the hinge disappears, the factor in (87) $\rightarrow EI/\ell$, and $\mathbf {K}_{MHk}$ becomes (29).

$\mu =0, \delta \rightarrow \infty$: $\mathbf {K}_{MHk}$ splits into two uncoupled $2\times 2$ blocks, accounting for the rigid support.

Fig. 11

Visualization of the piecewise-split deflection $v(\xi )$ and rotation $\theta (\xi )$ for the spring-midhinged element of Fig. 10a. Data: $EI=1$, $\ell =1$, $k_T=k_E=0$; nodal data and function expressions shown in figure. Dashed portions are extensions of $v^-(\xi )$ and $\theta ^-(\xi )$ into the right halfspan, and of $v^+(\xi )$ and $\theta ^+(\xi )$ into the left halfspan

1.8 2.8 Geometric Stiffness Matrix Reduction

The derivation of the geometric stiffness for the element of Fig. 10a triggers a challenge: how to eliminate (reduce) its three hinge freedoms. Since this operation is not covered in the FEM literature, a brief description of two model reduction methods is given.

Static Condensation. Let $\mathbf {K}_{Mu}$ and $\lambda \mathbf {K}_{Gu}$ denote the unreduced $7\times 7$ material and geometric stiffness matrices, respectively, that include the 3 hinge freedoms. Note that the coefficient $\lambda$ is explicitly taken out as dimensionless load factor. This can be done by writing $F=\lambda F_{ref}$, where $F_{ref}$ is a reference force, and then putting F back in the reduced equations.

The tangent stiffness combination $\mathbf {K}_{Tu}=\mathbf {K}_{Mu}+\lambda \mathbf {K}_{Gu}$ is treated by the well known static condensation procedure [80, 105] to eliminate the hinge freedoms. This yields a reduced tangent matrix $\mathbf {K}_{Tr}$ with rational entries in $\lambda$, which makes subsequent buckling analysis through a linear eigenproblem impossible. To recover that capability, that matrix is expanded in Taylor series in $\lambda$ up to linear terms:

$$\begin{aligned} \mathbf {K}_{Tr}=\mathbf {K}_{M}+\lambda \mathbf {K}_{G}+H.O.T. \end{aligned}$$

(88)

Higher order terms in (88) are discarded. Then $\mathbf {K}_{M}$ and $\mathbf {K}_{G}$ are taken as the reduced material and geometric stiffness matrices, respectively, and renamed $\mathbf {K}_{MHk}$ and $\mathbf {K}_{GHk}$ for the present element. This $\mathbf {K}_{MHk}$ does reproduce (84) but $\mathbf {K}_{GHk}$ is new. This procedure uses explicit linearization.

Kinematic Reduction. The optimal material stiffness (29) is assembled over the two beam portions and the spring stiffnesses directly added to the hinged freedoms. Partial condensation of this $7\times 7$ unreduced material matrix provides a $7\times 4$ transformation $\mathbf {T}_r$ matrix that connects the 3 hinge freedoms to the 4 retained ones. Let $\mathbf {u}=\left[ \begin{array}{cccc}v_1&\theta _1&v_2&\theta _2\end{array}\right] ^T$ and $\mathbf {u}_H=\left[ \begin{array}{cccc}v_0&\theta ^{-}_0&\theta ^{+}_0\end{array}\right] ^T$. The transformation can be presented in partitioned form:

$$\begin{aligned} \left[ \begin{array}{c}\mathbf {u}\\ \mathbf {u}_H\end{array}\right] = \left[ \begin{array}{c}\mathbf {I}_4 \\ \mathbf {T}_{H} \end{array}\right] \mathbf {u}\,\buildrel {\mathrm{def}}\over =\,\mathbf {T}_r \mathbf {u}. \end{aligned}$$

(89)

Here $\mathbf {I}_4$ is the $4\times 4$ identity matrix whereas

$$\begin{aligned} \mathbf {T}_H&={1\over C_T^{}} \left[ \begin{array}{cccc} T_{H11} &{} T_{H12} &{} T_{H11} &{} -\,T_{H12} \\ T_{H21} &{} T_{H22} &{} T_{H23} &{} T_{H24} \\ -T_{H23} &{} T_{H24} &{} -\,T_{H21} &{} T_{H22} \end{array}\right] , \end{aligned}$$

(90)

in which

$$\begin{aligned} T_{H11}&=384 \ell (1+\mu ) (1+\Phi ), \quad T_{H12}=96 \ell ^2 (2+\mu ) (1+\Phi ), \nonumber \\ T_{H21}&=-6 (192 (1+\mu )+\delta (8+\mu +8 \Phi +4 \mu \Phi )), \nonumber \\ T_{H22}&= \ell (192 (1-\mu +2 (2+\mu ) \Phi )+\delta (2 \Phi -1) (8+\mu +4 (2+\mu ) \Phi )), \nonumber \\ T_{H23}&= 6 (192{+}\mu (192{+}\delta {+}4 z \Phi )), \quad T_{H24}= \ell (\mu (2 \Phi {-}1) (192+\delta +4 \delta \Phi )-576), \nonumber \\ C_T&=4 \ell (1+\Phi ) (192 (1+\mu )+\delta (4+\mu +4 (1+\mu ) \Phi )). \end{aligned}$$

(91)

The reduced geometric stiffness is obtained by a congruential transformation:

$$\begin{aligned} \mathbf {K}_{GHk} = \mathbf {T}_r^T \mathbf {K}_{Gu} \mathbf {T}_r. \end{aligned}$$

(92)

The similarity to the well known Craig-Bampton mass reduction scheme [18] is plain. It can be shown that kinematic reduction produces the same geometric stiffness as static condensation followed by $\lambda$-linearization. Therefore (92) may be viewed as implicit linearization. The proof involves straightforward matrix algebra using the Schur complement [106] and is omitted.

1.9 2.9 Geometric Stiffness For Spring-MidHinged Beam

The kinematic reduction method is simpler to implement, and is the one used here. It was applied using the assembled quasi-optimal geometric stiffness (126) to built the $7 \times 7$ matrix $\mathbf {K}_{GH}$ . The result of the reduction (92) can be presented as a two-matrix combination:

$$\begin{aligned} \mathbf {K}_{GHk}= \mathbf {K}_{GH}+\mathbf {K}_{Gk}. \end{aligned}$$

(93)

Here $\mathbf {K}_{GH}$ is the geometric stiffness for a free hinge:

$$\begin{aligned} \mathbf {K}_{GH}= {F\over 320 \ell (1+\Phi )^2} \left[ \begin{array}{cccc} K_{GH11} &{} K_{GH12} &{} -\,K_{GH11} &{} K_{GH12} \\ K_{GH12} &{} K_{GH22} &{} -\,K_{GH12} &{} K_{GH24} \\ -K_{GH11} &{} -\,K_{GH12} &{} K_{GH11} &{} -\,K_{GH12} \\ K_{GH12} &{} K_{GH24} &{} -\,K_{GH12} &{} K_{GH22} \end{array}\right] . \end{aligned}$$

(94)

in which $K_{GH11}= 4 (97+180 \Phi +80 \Phi ^2)$, $K_{GH12}=2 \ell (17+20 \Phi )\ell$, $K_{GH22}= (97+ 180 \Phi +80 \Phi ^2) \ell ^2$, and $K_{GH24}= -(63+140 \Phi +80 \Phi ^2) \ell ^2$. The second matrix brings the effect of the springs:

$$\begin{aligned} \mathbf {K}_{Gk}= C_{Gk} \left[ \begin{array}{cccc} K_{Gk11} &{} K_{Gk12} &{} K_{Gk11} &{} -\,K_{Gk12} \\ K_{Gk12} &{} K_{Gk22} &{} K_{Gk12} &{} -\,K_{Gk22} \\ K_{Gk11} &{} K_{Gk12} &{} K_{Gk11} &{} -\,K_{Gk12} \\ -K_{Gk12} &{} -\,K_{Gk22} &{} -\,K_{Gk12} &{} K_{Gk22} \end{array}\right] , \end{aligned}$$

(95)

in which

$$\begin{aligned} K_{Gk11}&=48 \delta ^2 (97{+}44 \mu {+}7 \mu ^2{+}60 (1{+}\mu ) (3{+}\mu ) \Phi {+} 80 (1{+}\mu )^2 \Phi ^2). \nonumber \\ K_{Gk12}&=12 \ell \delta (-960 (8{+}5 \mu {+}\mu ^2{+}4 (1{+}\mu ) (2{+}\mu ) \Phi ){+} \delta (34{+}40 \Phi {+}\mu (11{+}2 \mu {+}20 (4{+}\mu ) \Phi ))),\nonumber \\ K_{Gk22}&=-\ell ^2 ((737280 \mu (3{+}2 \mu ){+}1920 \delta (48 (1{+}\Phi ){+} \mu (54{+}96 \Phi {+}11 \mu (1{+}4 \Phi ))) \nonumber \\&\quad +\delta ^2 (12 \mu (33{+}160 \Phi (1{+}\Phi )){+} 12 (63{+}20 \Phi (7{+}4 \Phi )){+}\mu ^2 (49{+}20 \Phi (19{+}44 \Phi ))))) \nonumber \\ C_{Gk}&={F\over 240 \ell (192 (1+\mu )+\delta (4+\mu +4 \Phi +4 \mu \Phi ))^2} \nonumber \\ \end{aligned}$$

(96)

This matrix simplifies considerably if the extensional spring is missing. Setting $\delta =0$ gives

$$\begin{aligned} \mathbf {K}_{Gk}= {F \mu (3+2 \mu ) \ell \over 12 \ell (1+\mu )^2} \left[ \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} -\,1 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} -\,1 &{} 0 &{} 1 \end{array}\right] . \end{aligned}$$

(97)

Fig. 12

Buckling problems solved with a single spring-midhinged element in Example 3. BE model used for both cases

Example 3

Two buckling problems involving a spring-midhinged Euler column are solved with the element developed in Appendix Sect. 2.7– Appendix Sect. 2.9. The BE model ($\Phi =0$) is assumed for both. Only one hinged element is used to discretize the member. Two discretizations are considered:

• Unreduced Model. The hinge DOF are retained in the LPB eigenproblem. This involves $7 \times 7$ stiffness matrices, which become to $5\times 5$ upon applying BC. Five critical loads are obtained, ordered $\lambda _{cr1}\le \ldots \lambda _{cr5}$.

• Reduced Model. The hinge DOF are eliminated from the geometric stiffness by kinematic reduction—as stated in Appendix Sect. 2.9, this gives identical results to static condensation followed by linearization. The $4 \times 4$ matrices are cut to $2\times 2$ upon applying BC, giving two critical loads: $\lambda _{cr1}\le \lambda _{cr2}$.

The first problem, shown in Fig. 12a, involves a pinned-pinned beam member of length L and constant bending rigidity EI, which is subjected to a compressive axial load $P=\lambda ^2 EI/L^2$. The member has a midspan hinge stabilized with a torsional spring of rigidity $k_T=\mu EI/L$, $\mu \ge 0$. The quasi-optimal geometric stiffness (126) is paired with the optimal BE material stiffness (30).

Results for the lowest critical load $P_{cr}=\lambda _{cr1}^2 EI/L^2$ are shown in Fig. 13a, where $\lambda _{cr1}^2$ is plotted as function of $\mu$ over [0, 1000]. They are compared with the exact solution with the buckling-exact element derived in Appendix 5. The exact characteristic equation is

$$\begin{aligned} (EI)^2\lambda ^4 \big (\lambda (\cos \lambda -1)+2\mu \sin \lambda \big )/ (2\mu L^4)=0, \end{aligned}$$

(98)

This agrees with [101, Sec. 2.5.2]. It gives symmetric buckling modes, starting with $\lambda _{cr1}$, which depend on $\mu$, as well as antisymmetric mode solutions $\lambda _{cr}=2 n \pi$, $n=1,2,\ldots$ for which the hinge is a node and the spring has no effect. As $\mu \rightarrow \infty$, (98) approaches $\sin \lambda =0$, which gives the classical Euler solutions. The unreduced model gives at least 3 correct places for all $\mu$ so no difference can be seen in the plots. The reduced model, however, loses accuracy as $\mu$ increases. The results for the second mode are way off and not shown.

Fig. 13

Results from problems of Example 3. left: pinned-pinned beam-column with rotationally constraint midhinge; right: pinned-pinned beam-column with extensionally restrained midhinge. s-mode and a-mode are appreviations for symmetric and antisymmetric buckling mode, respectively. Values of $\delta$ that maximize critical load in left figure are marked by crossed circles

The variation pictured in Fig. 12b has an extensional spring of rigidity $k_E=\delta EI/L^3$, $\delta \ge 0$. This configuration features eigenvalue crossings, with the attendant possibility of choosing $\delta$ to maximize the lowest critical load. The exact characteristic equation is

$$\begin{aligned} (EI)^2\lambda ^2 (\delta -4\lambda ^2) \sin (\lambda /2)/ (2 L^5)=0. \end{aligned}$$

(99)

There is a single symmetric mode with $\lambda _{cr}^2=\delta /4$, for which the submembers remain straight, taking no bending energy. It depends only on $\delta$ and gives the lowest critical load for $\delta \le 16\pi ^2$. The $\sin (\lambda /2)$ factor provides the antisymmetric mode solutions $\lambda _{cr}=2n\pi$, $n=1,2,...$ with a node at the hinge, which do not depend on $\delta$. The symmetric and first antisymmetric mode cross at $\delta _{\times }=16\pi ^2\approx 157.91$.

The FEM results for the first two modes are shown in Fig. 13b, along with the exact solution. The unreduced model captures the symmetric mode exactly, whereas for the second mode gives $\lambda _{cr}^2=48$ instead of $4\pi ^2$. Crossing occurs at $\delta _{\times }=192$. The reduced model overshoots the symmetric critical load over most of the range, gives a second critical load at $\lambda _{cr}^2=960/17\approx 56.47$, and predicts a crossing at $\delta _{\times }\approx 226$.

These simple but thoughtworthy examples underscore the quirky nature of model reduction. Condensation of the material stiffness has no side effects: equilibrium and nodal exactness are preserved. But reduction at the geometric stiffness level may cause significant accuracy loss. The underlying reason, explained in Appendix Sect. 2.8, is linearization—whether explicit or implicit—that discards higher order terms in the load factor. Thus errors can be expected to grow as the factor increases.

1.10 2.10 Off-Center Hinged Beam Element

Another generalization of the element in Fig. 9a is shown in Fig. 10b. The free hinge is now located at arbitrary isoparametric coordinate $-1\le \xi _H\le 1$. The derivation requires adjusting two Legendre polynomial bases over variable lengths and is omitted. The resulting material stiffness is the rank-one matrix

$$\begin{aligned} \mathbf {K}_{MH}= {3 EI \over \ell (1+3\xi _H^2+\Phi )} \left[ \begin{array}{cccc} {\displaystyle 4\over \displaystyle \ell ^2} &{} {\displaystyle 2 (1{+} \xi _H)\over \displaystyle \ell } &{} -{\displaystyle 4\over \displaystyle \ell ^2} &{} {\displaystyle 2 (1{-}\xi _H)\over \displaystyle \ell } \\ {\displaystyle 2 (1{+}\xi _H)\over \displaystyle \ell } &{} (1{+}\xi _H)^2 &{} -{\displaystyle 2 (1{+}\xi _H)\over \displaystyle \ell } &{} 1{-}\xi _H^2 \\ -{\displaystyle 4\over \displaystyle \ell ^2} &{} -{\displaystyle 2 (1{+}\xi _H)\over \displaystyle \ell } &{} {\displaystyle 4\over \displaystyle \ell ^2} &{} -{\displaystyle 2 (1{-}\xi _H)\over \displaystyle \ell } \\ {\displaystyle 2 (1{-}\xi _H)\over \displaystyle \ell } &{} 1{-}\xi _H^2 &{} -{\displaystyle 2 (1{-}\xi _H)\over \displaystyle \ell } &{} (1{-}\xi _H)^2\end{array}\right] . \end{aligned}$$

(100)

This reduces to (83) for $\xi _H=0$. The limit cases $\xi _H=\pm 1$ are not excluded. If springs are attached to the hinge, the matrices derived in previous subsections must be adjusted for a variable $\xi _H$. The entries become exceedingly complicated and are not given here.

Fig. 14

Hinge position optmization example. a centrally point-loaded fixed-fixed BE beam with 2 symmetrically placed hinges; b one element discretization by taking advantage of symmetry

Example 4

The material stiffness (100) for an arbitrarily placed hinge is applied to a simple design optimization problem that can be entirely done by hand in closed form. This is pictured in Fig. 14a. A Bernoulli–Euler (BE) fixed-fixed prismatic beam AB of total span L is loaded by a point force P (positive upward) at its center C. Two hinges are symmetrically placed at distances ${\textstyle {1\over 2}}\alpha L$ from the end supports. Because of the symmetry about C, it is possible to model the problem with just one hinged beam element, say the left half-span, as illustrated in Fig. 14b. The displacement support conditions are: $v_1=\theta _1=\theta _2=0$. This leaves only one unknown DOF: the midspan deflection $v_2=v_C$ under ${\textstyle {1\over 2}}P$. Question: for which hinge configuration, defined by $\alpha$ as design variable, is $v_C$ minimized in magnitude?

Mathematically this is a constrained optimization problem:

$$\begin{aligned} \hbox {find}\quad \min \limits _\alpha \left| v_C(\alpha )\right| \quad \hbox { for }\alpha \in [0,1]. \end{aligned}$$

(101)

Does such a minimum exist? For both $\alpha =0$ (a simply supported beam) and $\alpha =1$ (a doubly hinged cantilever), textbooks give $v_C=PL^3/(48EI)$, which provides checks on the computations by setting $\xi _H=\pm 1$. For $\alpha$ in between, $|v_C|$ gets smaller, so there is in fact an optimal value, which is readily computable by hand.

The element stiffness matrix (100), upon setting $\xi _H\rightarrow 2\alpha -1$, $\ell \rightarrow L/2$ and $\Phi =0$, is used as master stiffness. Applying the BCs $v_1=\theta _1=\theta _2=0$ leaves only one equation: $12 EI/(\ell ^3 (1+3\xi _H^2)) v_2= P/2$, which solved for $v_2$ gives

$$\begin{aligned} v_2=v_C= {P/2\over 12 EI /(\ell ^3 (1+3\xi _H^2)) }= { P \ell ^3 (1+3\xi _H^2) \over 24 EI }= {P L^3 (1+3\xi _H^2) \over 192 EI }. \end{aligned}$$

(102)

By inspection $|v_C|$ is minimized for $\xi _H=0$, or $\alpha =1/2$. That is, the hinges should be placed at the quarterspans. If this is done, the center deflection is $PL^3/(192EI)$ as compared to $PL^3/(48EI)$ when $\xi _H=\pm 1$. That is a gain of 4 in deflection stiffness under the center point load. If shear flexibility is considered: $\Phi >0$, the optimal $\xi _H$ moves closer to the supports.

Note: $PL^3/(192EI)$ is the same as the center defection of a centrally-point-loaded, fixed-fixed beam without hinges. The advantages of placing hinges are: (1) the absolute maximum moment is minimized, as can be easily verified, and (2) thermal expansion effects can be abated with appropriate sliding devices at the hinges.

Appendix 3: Templates and Modified Equations

This Appendix presents two mathematical techniques for deriving optimal elements: templates and modified differential equations. While the combination can be used on a variety of problems, they are better exhibited in the context of specific applications. In accordance with the main subject of the paper, we construct optimal elements for the Timoshenko plane beam model, with BE as special case. The rationale is plain: if the combination does not work for beams there is little hope for plates and shells, which carry more complicated forms of rotational DOF. References for more general applications are provided in Appendix 8 along with a brief historical summary.

The element superscript e is again omitted for brevity.

1.1 3.1 Templates: Concept and Examples

A finite element template, or template for short, is a set of FEM equations that contains symbolic free parameters. Their numerical values are picked to improve performance, and in some cases to achieve optimality in a defined sense. As an example, take a two-node linear elastostatics plane beam element with the standard 4 degrees of freedom. The generic material stiffness equations are

$$\begin{aligned} \left[ \begin{array}{cccc} K_{M11} &{} K_{M12} &{} K_{M13} &{} K_{M14} \\ K_{M12} &{} K_{M22} &{} K_{M23} &{} K_{M24} \\ K_{M13} &{} K_{M23} &{} K_{M33} &{} K_{M34} \\ K_{M14} &{} K_{M24} &{} K_{M34} &{} K_{M44} \end{array}\right] \left[ \begin{array}{cccc}v_1 \\ \theta _1 \\ v_2 \\ \theta _2\end{array}\right] = \left[ \begin{array}{cccc}f_1 \\ m_1 \\ f_2 \\ m_2\end{array}\right] \quad \Rightarrow \quad \mathbf {K}^e\mathbf {u}^e=\mathbf {f}^e. \end{aligned}$$

(103)

where the only pre-assumption is matrix symmetry. The values of the 10 independent stiffness entries: $K_{M11}$ through $K_{M44}$ may be viewed as parameters. This is called a universal template, since it include all possible symmetric stiffness forms for a linear model. But for optimality studies this is far too general. The number of free parameters may be cut down by applying behavioral and fabrication constraints. Two examples for the above element:

• Applying the rigid body motion $\mathbf {v}^e=\left[ \begin{array}{cccc}1&0&1&0\end{array}\right]$ must give a zero force vector. Consequence: the third column (row) must be the negated first column (row): $K_{M13}=K_{M31}=-K_{M11}$, etc.

• If the element is prismatic (constant contitutive and fabrication properties), symmetry requires $K_{M33}=K_{M11}$, etc.

Instead of going over individual constraints a more systematic procedure is to use natural modes. For beam elements, this is done in Appendix 2 using Legendre polynomials with a minor tweak in the cubic term. For both the Timoshenko and BE models of prismatic (uniform) elements it gives the material stiffness template

$$\begin{aligned} \mathbf {K}_M = \mathbf {K}_b+\mathbf {K}_h = {EI\over \ell }\left[ \begin{array}{cccc}0&{}0&{}0&{}0\\ 0&{}1&{}0&{}-1\\ 0&{}0&{}0&{}0\\ 0&{}-1&{}0&{}1\end{array}\right] +\beta {3EI\over \ell ^3}\left[ \begin{array}{cccc} 4 &{}-2\ell &{}-4 &{}-2\ell \\ -2\ell &{} \ell ^2 &{} 2\ell &{} \ell ^2 \\ -4 &{} 2\ell &{}4 &{} 2\ell \\ -2\ell &{} \ell ^2 &{} 2\ell &{} \ell ^2 \end{array}\right] . \end{aligned}$$

(104)

This form, presented in [27], has one free parameter: the scalar $\beta$. Plainly it is easier to deal with one parameter than with ten.

Some notation: $\mathbf {K}_b$ and $\mathbf {K}_h$ are the basic and higher order stiffness matrices, respectively, a terminology originated as described in Appendix Sect. 8.5 Their sum $\mathbf {K}_M$ has the correct rank of 2 and is positive semidefinite if $\beta >0$. Parameter $\beta$ scales only the higher order stiffness $\mathbf {K}_h^e$.

1.2 3.2 Geometric Stiffness Templates

Templates for the geometric element stiffness necessarily contain more parameters. Their generic configuration is

$$\begin{aligned} \mathbf {K}_{G}^e = \mathbf {H}^T \mathbf {S}_{G} \mathbf {H}, \end{aligned}$$

(105)

where $\mathbf {S}_{G}$ is a natural stiffness matrix endowed with $n_\beta$ free parameters, which are denoted $\beta _i$, $i=1,2,\ldots n_\beta$, and $\mathbf {H}$ is the transformation matrix given in (77) that maps natural coordinates to node displacements.

Matrix $\mathbf {S}_{G}$ is set up by including free parameters in the $\mathbf {V}_1$ covariant matrices derived in Appendix Sect. 2.2. For arbitrary values of $\chi$, only four independent entries of $\mathbf {V}_1$ are nonzero. It is convenient to parametrize those as

$$\begin{aligned} \mathbf {S}_{G}={4F \over \ell } \left[ \begin{array}{cccc}0 &{} 0 &{} 0 &{} 0 \\ 0 &{} \beta _1 &{} 0 &{} \beta _4\\ 0 &{} 0 &{} \beta _2 &{} 0 \\ 0 &{} \beta _4 &{} 0 &{} \beta _3\end{array}\right] , \end{aligned}$$

(106)

which will be called the natural geometric stiffness template. A simple computation shows that the factor $4F/\ell$ forces $\beta _1=1$, leaving four free parameters: $\beta _2,\beta _3,\beta _4,\chi$ for the optimal $\mathbf {K}_G$ search. Results are reported in Appendix Sect. 3.5 and following sections.

Fig. 15

Flowchart of the modified differential equation (MoDE) method applied to FEM. Adapted from [36], which addressed the convection-diffusion problem. Step 4 may be a difficult one requiring human decisions; for instance Taylor series in an patch dimension parameter and truncation. Steps 1–3 are automatic but symbolic processing normally requires a computer algebra program

1.3 3.3 Optimal BE Beam Material Stiffness

How do we pick the best $\beta$ in (104)? As announced, the modified differential equation method will be used. This is based on knowledge at the Mechanics of Materials level: the governing ODE for computing deflections (called strong form in the literature) for a prismatic BE beam is (9) with $\Phi =0$:

$$\begin{aligned} EI \, v''''(x) = q(x). \end{aligned}$$

(107)

where v(x) is the lateral deflection and q(x) the applied distributed load—see Fig. 2—and primes denote derivatives respect to x. Task: find the optimal $\beta$, with only (104) and (107) as given. By “optimal” is meant achieving a “best fit” in the sense explained below. No variational or weak forms are considered. Although both (104) and (107) are undergraduate level material, the fitting task is far from trivial.

The key idea is to construct an ODE that is satisfied exactly by the discrete system at the nodes. This is called the modified differential equation or MoDE. Typically it has infinite order, in which case it is labeled IOMoDE. If an instance of $\beta$ in (104) is found that

• Reduces the IOMoDE to a finite-order modified differential equation, called FOMoDE, and

• The FOMoDE matches (107),

the FEM model is nodally exact and it is happy hour. But often an exact match is impossible. If so the goal is to get the best fit possible in low powers of the element length $\ell$. This is done by expanding the IOMoDE in Taylor series in $\ell$ about $\ell =0$, and truncating the expansion after the first nonzero term.

The procedure is flowcharted in Fig. 15. The process begins from a two-element patch extracted from an infinite regular lattice pictured in Fig. 16.

Since computations are necessarily symbolic and error prone if done by hand, they are better carried out with the assistance of a computer algebra system. An algorithmic description is given in Appendix 4 along with a supporting Mathematica implementation. Running the steps of Fig. 15 as detailed in that Appendix gives the infinite order modified differential equation (IOMoDE)

$$\begin{aligned} EI \, v''''(x)=q(x) +{1-\beta \over \beta }\left[ -{\ell ^2 q''(x)\over 12}+ {\ell ^4 q''''(x)\over 720}-{\ell ^6 q''''''(x)\over 30240}+\ldots \right] . \end{aligned}$$

(108)

The optimal choice is clear: setting $\beta =1$ makes (108) a FOMoDE identical to (107). Replacing into (104) gives the well known result (30), Since we have exactly matched (107), the element endowed with this matrix is nodally exact (NE). There are two additional points:

• Nodal exactness is preserved if EI and $\ell$ vary from element to element. This requires more elaborate symbolic computations.

• The high order derivatives in (108) must be interpreted in a formal sense; see [46, Ch. X].

Fig. 16

Extracting a two-element patch from an infinite beam lattice

1.4 3.4 Timoshenko Beam Material Stiffness

Next, try the template (104) on a Timoshenko beam element. The governing ODE for the no-axial-load case is (9), reproduced here for convenience:

$$\begin{aligned} EI \, v''''(x) = q(x)-{\Phi h^2\over 12} q''(x). \end{aligned}$$

(109)

Going through the steps of Fig. 15, we get the IOMoDE

$$\begin{aligned} EI v''''= q(x)-{1-\beta \over 12\beta }\ell ^2 q''(x)+ {1-\beta -\beta \Phi \over \beta (1+\Phi )} \left[ {\ell ^4 q''''(x)\over 720}+{\ell ^6 q''''''(x)\over 30240}+\ldots \right] . \end{aligned}$$

(110)

Setting $\beta =1+\Phi$ matches (109). Inserting in (104) gives (29), which is nodally exact.

1.5 3.5 Optimal BE Geometric Stiffness

We consider first the BE model under LPB (Linearized Prebuckling) assumptions. Assume that the beam element is under a uniform axial force $F=-P$ (P is positive compression), which produces a state of initial axial stress. The governing ODE in terms of deflections becomes $EI v''''(x)+F v''(x)=q(x)$. For the following analysis the lateral load q(x) is irrelevant and is set to zero so the target ODE is

$$\begin{aligned} EI v''''(x)+F v''(x)=0, \end{aligned}$$

(111)

The search for optimal $\mathbf {K}_G$ starts with the natural geometric stiffness 4-parameter template $\mathbf {S}_G$ in (106). The physical geometric stiffness is $\mathbf {K}_G=\mathbf {H}^T \mathbf {S}_G \mathbf {H}$, where the matrix $\mathbf {H}$ given in (77) brings one more parameter: $\chi$, from the parametrized Legendre basis; see (68). The five parameters can be reduced to four by setting $\beta _1=1$, which follows automatically from the $\mathbf {S}_G$ scaling factor. The modified equation analysis produces the IOMode

$$\begin{aligned} EI v''''+Fv'' = v'' \big (C_2 (F/EI)^2 \ell ^2+C_4 (F/EI)^4 \ell ^4+C_6 (F/EI)^6 \ell ^6+H.O.T.\big ) \end{aligned}$$

(112)

where the “higher order terms” are $O(\ell ^8)$ and higher. The $C_i$’s are dimensionless functions of the template parameters. Two coefficients that can be zeroed by proper choice of parameters are

$$\begin{aligned} C_2&= (-100 \beta _2-3 (25+4 \beta _4^2+4 \beta _4 (-15+\chi )- 30 \chi +\chi ^2))/3600, \nonumber \\ C_4&= (10000 \beta _2^2+300 \beta _2 (-175+4 \beta _4^2+4 \beta _4 (-15+\chi )- 30 \chi +\chi ^2) \nonumber \\&+9 (6625+16 \beta _4^3 (-15+\chi )-1500 \chi +950 \chi ^2-60 \chi ^3+ \chi ^4+4 \beta _3 (-15+2 \beta _4+\chi )^2\nonumber \\&+\, 20 \beta _4^2 (145-30 \chi +\chi ^2)+8 \beta _4 (-375+475 \chi - 45 \chi ^2+\chi ^3)))/12960000. \end{aligned}$$

(113)

Solving $C_2=C_4=0$ for $\beta _2$ and $\beta _3$, both of which appear linearly there, gives

$$\begin{aligned} \beta _2&= (-3 (25+4 \beta _4^2+4 \beta _4 (-15+\chi )-30 \chi +\chi ^2))/100, \nonumber \\ \beta _3&= -(11625+16 \beta _4^3 (-15+\chi )-7500 \chi +1150 \chi ^2-60 \chi ^3+\chi ^4\nonumber \\&\quad +20 \beta _4^2 (185-30 \chi +\chi ^2)+8 \beta _4 (-1875+575 \chi \nonumber \\&\quad -45 \chi ^2+\chi ^3))/(4 (-15+2 \beta _4+\chi )^2). \end{aligned}$$

(114)

Upon choosing $\beta _2$ and $\beta _3$ as per (114), the $\ell ^6$ coefficient in (112) becomes

$$\begin{aligned} C_6&= -(469125+336 \beta _4^4+672 \beta _4^3 (-15+\chi )-187500 \chi +25150 \chi ^2-1260 \chi ^3\nonumber \\&+ 21 \chi ^4+8 \beta _4^2 (12575-1890 \chi +63 \chi ^2)+8 \beta _4 (-46875+12575 \chi \nonumber \\&- 945 \chi ^2+21 \chi ^3))/(6048000 (-15+2 \beta _4+\chi )^2). \end{aligned}$$

(115)

Examination of $C_6$ shows that it does not vanish for any real $\beta _4$ and $\chi$. Thus optimality reduces to finding minima of $C_6$ as function of $\beta _4$ and $\chi$. Solving $\partial C_6/\partial \chi =0$ for $\beta _4$ yields four roots, two of which are real. For convenience define

$$\begin{aligned} \upsilon _1=15/2-\root 4 \of {375}\approx 3.099441316033033, \quad \upsilon _2=15/2+\root 4 \of {375}\approx 11.900558683966967. \end{aligned}$$

(116)

The two real roots are

$$\begin{aligned} \beta _{4|1} =\upsilon _1-{\textstyle {1\over 2}}\chi , \quad \beta _{4|2}=\upsilon _2-{\textstyle {1\over 2}}\chi . \end{aligned}$$

(117)

Replacing either root into (115) gives

$$\begin{aligned} C_{6min}= {1\over 1890}-{1\over 480}\sqrt{15}\approx -8.8138245\times 10^{-6}\approx {1/113458.} \end{aligned}$$

(118)

This minimum is independent of $\chi$ so one cannot do better than this. The optimum—within the template assumption (112)—has been reached but we are free to pick $\chi$, which produces two $\beta _4$ as per (117). Inserting these into (114) gives

$$\begin{aligned} \beta _{2|1}&= \beta _{2|2} = 6-3 \sqrt{3/5}, \nonumber \\ \beta _{3|1}&=-5 \sqrt{15}+\root 4 \of {375} (\chi -15)+(425+(\chi -30) \chi )/4, \nonumber \\ \beta _{3|2}&=-5 \sqrt{15}-\root 4 \of {375} (\chi -15)+(425+(\chi -30) \chi )/4, \nonumber \\ \end{aligned}$$

(119)

The two parameter sets: $[\beta _2,\beta _{3|1},\beta _{4|1},\chi ]$, and $[\beta _2,\beta _{3|2},\beta _{4|2},\chi ]$, can be used to form two optimal geometric stiffness matrices called $\mathbf {K}_{Gop1}$ and $\mathbf {K}_{Gop2}$, respectively. Their expressions are:

$$\begin{aligned} \mathbf {K}_{Gop1}=c_{op} \left[ \begin{array}{cccc} a_1 &{} a_2 \ell &{} -\,a_1 &{} a_2 \ell \\ a_2 \ell &{} a_3 \ell ^2 &{} -\,a_2 \ell &{} a_4 \ell ^2 \\ -a_1 &{} -\,a_2 \ell &{} a_1 &{}-a_2 \ell \\ a_2 \ell &{} a_4 \ell ^2 &{} -\,a_2 \ell &{} a_3 \ell ^2\end{array}\right] , \mathbf {K}_{Gop2}=c_{op} \left[ \begin{array}{cccc} a_1 &{} a_5 \ell &{} -\,a_1 &{} a_5 \ell \\ a_5 \ell &{} a_6 \ell ^2 &{} -\,a_5 \ell &{} a_7 \ell ^2 \\ -a_1 &{}-a_5 \ell &{} a_1 &{}-a_5 \ell \\ a_5 \ell &{} a_7 \ell ^2 &{} -\,a_5 \ell &{} a_6 \ell ^2\end{array}\right] . \end{aligned}$$

(120)

where $c_{op}=F/(12 \sqrt{15} \ell )$, $a_1=24 \sqrt{15}-36$, $a_2=12 \sqrt{15}-18-6 \root 4 \of {135}$, $a_3=11 \sqrt{15}-12-6 \root 4 \of {135}$, $a_4=7 \sqrt{15}-6-6 \root 4 \of {135}$, $a_5=6 (2 \sqrt{15}+\root 4 \of {135}-3)$, $a_6=11 \sqrt{15}-12+6 \root 4 \of {135}$, and $a_7=7 \sqrt{15}-6+6 \root 4 \of {135}$. Both $\mathbf {K}_{Gop1}$ and $\mathbf {K}_{Gop2}$ are independent of $\chi$ (which cancels out in the $\mathbf {H}^T \mathbf {S}_G \mathbf {H}$ transformation), so for convenience one may preset $\chi =0$ in (119). In the example below both versions gave exactly the same results; it is not known if this holds for other problems. In any case $\mathbf {K}_{Gop1}$ has the better condition number and is the recommended one.

These optimal matrices do not seem to appear in the literature.

1.6 3.6 Quasi-Optimal BE Geometric Stiffness

The optimal matrices (120) have complicated coefficients, which leads to unwieldy expressions when extending to the Timoshenko model. It is possible to find “quasi-optimal” forms that keep $C_2=C_4=0$ with a slightly larger $C_6$, but have parameter values that are simple integers. Two instances are

$$\begin{aligned}&\beta _1=1, \quad \beta _2=3, \quad \beta _3=10, \quad \beta _4=0, \quad \chi =5. \nonumber \\&\beta _1=1, \quad \beta _2=3, \quad \beta _3=11, \quad \beta _4=1, \quad \chi =3. \end{aligned}$$

(121)

Both choices produce the same physical geometric stiffness:

$$\begin{aligned} \mathbf {K}_{Gqo}={F\over 60 \ell }\left[ \begin{array}{cccc} 84 &{} 12 \ell &{} -\,84 &{} 12 \ell \\ 12 \ell &{} 11 \ell ^2 &{} -\,12 \ell &{} \ell ^2 \\ -84 &{} -\,12 \ell &{} 84 &{} -\,12 \ell \\ 12 \ell &{} \ell ^2 &{} -\,12 \ell &{} 11 \ell ^2\end{array}\right] \end{aligned}$$

(122)

The modified equation analysis gives $C_6=-1/37800$, which is about 20% bigger than for the optimal. This matrix will be labeled quasi-optimal. It was derived for an advanced FEM course in 1998 but never published.

Fig. 17

Buckling of cantilever beam-column used as benchmark in Examples 5 and 6

Table 8 First two critical loads of cantilevered beam-column, BE model

Example 5

Table 8 collects results for the buckling analysis of a prismatic cantilever beam-column using the BE model. See Fig. 17. The number of elements $N_e$ varies from 1 to 16. All elements have the same length $\ell =L/N_e$. The material matrix is the nodally-exact one (29). Several geometric stiffness matrices are compared: the two optimal forms (120) — (both versions give identical results), quasi-optimal (122); consistent (38); and “bar” (123). The latter is obtained by setting $\beta _1=1$, $\beta _2=\beta _3=\beta _4=0$ and $\chi =5$ in (112), which produces

$$\begin{aligned} \mathbf {K}_{Gbar}={2F\over \ell }\left[ \begin{array}{cccc} 1 &{} 0&{} -\,1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ -1 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array}\right] \end{aligned}$$

(123)

This instance has historic importance as being the first geometric stiffness published in the open literature by Turner et. al. [94] and Argyris [4]. It was derived for the bar element—using linear shape functions–but applied to beams. [The term “stringer” was used instead of “bar” in the older literature.] Note that both “consistent” and “bar” instances converge from above; the former as $O(\ell ^4)$ and the latter as $O(\ell ^2)$. Both “optimal” and “quasi-optimal” converge as $O(\ell ^6)$ but not monotonically.

The computed accuracy for the first two critical loads is graphically summarized in the log-log plots of Fig. 18 as correct digits in $\lambda _{cr}$ versus number of elements. The impressive accuracy delivered by the optimal and quasi-optimal instances should be noted. For example both optimal versions give the first critical load to nearly 4 places with one element! Differences between optimal and quasi-optimal results are practically meaningless.

Fig. 18

Log-log plots of BE-modeled buckling prediction accuracy of cantilever beam-column member, for modes 1 and 2 in a and b, respectively. This summarizes results collected in Table 8. Horizontal axis: number of elements along member (1 through 16) in 2-log spacing. Vertical axis: number of correct digits of predicted $\lambda _{cr}$, obtained by taking the negated decimal log of the computed value absolute error

1.7 3.7 Quasi-Optimal Timoshenko Geometric Stiffness

For the Timoshenko model, the governing ODE to be matched (again ignoring the lateral load) is $(EI(1+F/GA_s) v''''(x)+ Fv''=0$, which may be rearranged to the canonical form

$$\begin{aligned} v''''(x)+{F\over EI+F \Phi \ell ^2/12} v''(x), \quad \hbox {or}\quad v''''(x)+k^2 v''(x)=0. \end{aligned}$$

(124)

As noted in Appendix Sect. 3.6, when passing to the Timoshenko beam model it is convenient to generalize the quasi-optimal BE geometric stiffness (122), which has simpler entries. The MoDE analysis shows that only the $\beta _3$ coefficient needs to be adjusted:

$$\begin{aligned}&\beta _1=1, \quad \beta _2=3, \quad \beta _3=10+25 \Phi , \quad \beta _4=0, \quad \chi =5. \nonumber \\&\beta _1=1, \quad \beta _2=3, \quad \beta _3=11+25 \Phi , \quad \beta _4=1, \quad \chi =3. \end{aligned}$$

(125)

Both choices produce the same matrix:

$$\begin{aligned} \mathbf {K}_{Gqo}={P\over 60 \ell (1+\Phi )^2} \left[ \begin{array}{cccc} b_1 &{} b_2\ell &{} -\,b_1 &{} b_2\ell \\ b_2\ell &{} b_3\ell ^2 &{} -\,b_2\ell &{} b_4\ell ^2 \\ -b_1 &{} -\,b_2\ell &{} b_1 &{} -\,b_2\ell \\ b_2\ell &{} b_4\ell ^2 &{} -\,b_2\ell &{} b_3\ell ^2 \end{array}\right] \end{aligned}$$

(126)

where $b_1=12(7+5 \Phi (3+\Phi ))$, $b_2=6(2+5\Phi )$, $b_3=11+5\Phi (5+\Phi )$, and $b_4=1+5(1-\Phi )\Phi$. This reduces to (122) if $\Phi =0$.

Table 9 First two critical loads of cantilevered beam-column, Timoshenko model

Fig. 19

Log-log plots of Timo-modeled buckling prediction accuracy of cantilever beam-column member, for modes 1 and 2 in a and b, respectively. See legend of Fig. 18 for details

Fig. 20

Arch frame optimization problem for Example 7. a structure; b finite element idealization

Example 6

Table 8 show results for the benchmark of Fig. 17, but now for a Timoshenko model with $\Phi _0=12EI/(GA_s L^2)=1/2$, which represents a high shear flexibility. Again 1 through 16 elements are used along the span. The element-level shear flexibility $\Phi$ need to be adjusted as per element length; thus

$$\begin{aligned} \Phi ={12 EI\over GA_s \ell ^2}= \Phi _0 N_e^2, \end{aligned}$$

(127)

where $N_e$ is the number of elements. This coefficient grows rapidly as $N_e$ increases; e.g., $\Phi =128$ for $N_e=16$ and $\Phi _0=1/2$. This has a significant degradation effect on the convergence rate, as shown below.

Results are shown in Table 9 for two geometric stiffnesses: quasi-optimal (126) and consistent (36). Accuracy results are graphically summarized as log-log plots in Fig. 19, which is displayed in the same vertical scale as Fig. 18, A glance plainly shows that both accuracy and convergence rate are significantly degraded. This can be explained by observing that the first nonzero $C_i$ coefficient in the IOMoDE are

$$\begin{aligned} \hbox {consistent:}\quad C_4={P^3 \Phi \over 144 (EI)^3},\quad \hbox {quasi-optimal:}\quad C_6={P^4 (8+50\Phi +105\Phi ^2)\over 302400 (EI)^4 (1+\Phi )}, \end{aligned}$$

(128)

so both grow as $O(\Phi )$ as $\Phi \rightarrow \infty$. From (127) note that a nonzero $\Phi$ grows as $1/\ell ^2$ as $\ell \rightarrow 0$. That cuts the convergence of the consistent version from $O(\ell ^4)$ to $O(\ell ^2)$ and that of the quasi-optimal one from $O(\ell ^6)$ to $O(\ell ^4)$. For the bar version (results not shown) it cuts $O(\ell ^2)$ to O(1), whence it does not converge for the Timoshenko model.

Fig. 21

Results for arch frame buckling optimization problem for selected values of $\rho$ and $\Phi$. Here $\lambda _1$ is associated with a symmetric buckling mode in which node 2 displaces vertically without rotation. On the other hand, $\lambda _2$ and $\lambda _3$ are lateral modes in which node 2 rotates and displaces horizontally. Note that the vertical scale of the rightmost two plots is different

Example 7

This example illustrates the application of an accurate Timoshenko geometric stiffness to a buckling optimization problem. The arch frame shown in Fig. 20 consists of two identical beam-column members clamped at their supports, and loaded by a vertical load $f_y=\lambda P_ref$ at the crown, where $\lambda >0$ if the load goes downward. Material and fabrication properties are indicated in the figure. The FEM idealization is shown in Fig. 20b. Each member is modeled by one Timoshenko beam-column element. The optimal material stiffness (29) and the quasi-optimal geometric stiffness (126) are used for the FEM model. The $6\times 6$ beam-column element stiffness matrices incorporate the effect of axial deformations and are appropriately rotated to the global positions through appropriate kinematic transformations.

The optimization question is: which rise angle $\varphi$ maximizes the critical buckling load? Note that since L is taken as unity—see (129) below—the member fabrication cost is the same for any rise angle.

It is convenient to reduce all equations to dimensionless form since the optimum angle $\varphi$ is also a dimensionless variable. This is done as follows:

$$\begin{aligned} E=1,\quad A=1,\quad L=1,\quad S=2 L \cos \varphi ,\quad H=L \sin \varphi , \quad I=\rho ^2 A=\rho ^2. \end{aligned}$$

(129)

Here $\rho$ plays the role of radius of gyration of the cross section. The shear coefficient $\Phi =12 EI/(GA_s L^2)=12\gamma ^2/(GA_s)$ is kept as a variable that indirectly defines $GA_s=12\gamma ^2/\Phi$. Three dimensionless variables are left: $\varphi$, $\gamma$, and $\Phi$. The load $f_Y$ is made dimensionless by taking $P_{ref}=-1$. The master stiffness matrices are $9\times 9$, but reduce to $3\times 3$ upon deleting all equations associated with nodes 1 and 3. Denoting $c=\cos \varphi$, $s=\sin \varphi$ and $c_2=\cos ^2\varphi$, the reduced material and geometric stiffness matrices are

$$\begin{aligned} \mathbf {K}_M=\left[ \begin{array}{cccc}2c^2+{\displaystyle 24 s^2 \rho ^2\over \displaystyle 1+\Phi } &{} 0 &{} {\displaystyle 12s \rho ^2 \over \displaystyle 1+\Phi } \\ 0 &{} 2s^2+{\displaystyle 24 c^2 \rho ^2\over \displaystyle 1+\Phi } &{} 0 \\ {\displaystyle 12s \rho ^2\over \displaystyle 1+\Phi } &{} 0 &{} {\displaystyle 2\rho ^2 (4+\Phi )\over \displaystyle 1+\Phi }\end{array}\right] , \quad \mathbf {K}_G=\left[ \begin{array}{cccc}K_{G11} &{} 0 &{} K_{G13} \\ 0 &{} K_{G22} &{} 0 \\ K_{G13} &{} 0 &{} K_{G33}\end{array}\right] , \end{aligned}$$

(130)

where for the geometric stiffness (126):

$$\begin{aligned} K_{G11}&={s(-12-25\Phi -10\Phi ^2+c_2 (2+5\Phi ))\over 10(1+\Phi ) (12 c^2 \rho ^2+s^2 (1+\Phi )) }, \ \ \ K_{G13}={-s^2 (2+5 \Phi )\over 10 (1+\Phi ) (12 c^2 \rho ^2+s^2 (1+\Phi ))},\nonumber \\ K_{G22}&={-s (12+25 \Phi +10 \Phi ^2+c_2 (2+5 \Phi ))\over 10 (1+\Phi ) (12 c^2 \rho ^2+s^2 (1+\Phi )) }, \ \ K_{G33}= { -(s (11+25 \Phi +5 \Phi ^2)\over 60 (1+\Phi ) (12 c^2 \rho ^2+s^2 (1+\Phi )) }. \end{aligned}$$

(131)

The determinant d of $K_M+\lambda K_G$ is taken. d is a cubic polynomial in $\lambda$, with coefficients that depend on $\varphi$, $\rho$ and $\Phi$. On giving numeric values to $\rho$ and $\Phi$, $d=0$ provides three critical load factors $\lambda _1$, $\lambda _2$ and $\lambda _3$ as functions of $\varphi$. These are plotted in Fig. 21 over $\varphi \in [30^\circ ,90^\circ ]$ for $\rho =1/5, 1/10,1/20$ combined with $\Phi =0,1/3$. (Shallow arches with $\varphi <30^\circ$ are omitted since they are subject to another type of instability: snap-through, which is not covered by the LPB model.) The optimal $\varphi$ is marked on the plots; it is the max distance from the $\lambda =0$ axis to the lowest plot. It may occur at a crossover, or at a local maximum. As can be expected, the influence of $\rho$ is important. The reduction in buckling strength from $\Phi =0$ to $\Phi =1/3$ is significant because that value implies a high shear flexibility; however, the optimum locations angles do not change significantly.

Rerunning the $\Phi =0$ case with the optimal BE geometric stiffness hardly changes the plots.

Appendix 4: MoDE Implementation Example

This Appendix provides computational details on the modified differential equation method (MoDE) flowcharted in Fig. 15. The first section goes over the analysis for a material stiffness template, in equation-by-equation detail. The second section presents the actual Mathematica implementation.

1.1 4.1 MoDE Analysis Example

The MoDE analyis details are presented here for the material stiffness of the BE beam model. Consider the lattice shown in Fig. 16a, from which a generic two-element patch is extracted as indicated in Fig. 16b. The identical elements are labeled as (1) and (2), respectively. The lateral load q(x) is not necessarily periodic (this is NOT Fourier analysis). The FEM patch equations at node j are

$$\begin{aligned} {EI\over \ell ^3}\left[ \begin{array}{ccccccc}-12&{} -\,6 \ell &{} 24&{} 0&{} -12&{} 6 \ell \\ 6 \ell &{} 2 \ell ^2&{} 0&{} 8 \ell ^2&{} -\,6 \ell &{} 2 \ell ^2\end{array}\right] \mathbf {v}_p= \left[ \begin{array}{cccc}f_j\\ m_j \end{array}\right] , \end{aligned}$$

(132)

Here $\mathbf {v}_p= \left[ \begin{array}{cccccc}v_i&\theta _i&v_j&\theta _j&v_k&\theta _k\end{array}\right] ^T$ is the 6-vector collecting the patch DOF, while $f_j$ and $m_j$ are consistent nodal loads to be computed as function of the lateral load q(x) and its derivatives at $x_j$ as shown below. Expand v(x), $\theta (x)$ and q(x) in Taylor series about $x=x_j$, truncating at $n+1$ terms for v and $\theta$, and $m+1$ terms for q. (Typically $m=n-4$ when doing a material stiffness but $m=0$ for a geometric stiffness.) Using $\zeta =(x-x_j)/\ell$ as dimensionless coordinate going from i to k, the series are

$$\begin{aligned} v(x)&=v_j+ \zeta \ell v_j'+{(\zeta ^2 \ell ^2/2!)} v_j''+ \ldots +{(\zeta ^n \ell ^n/n!)} v_j^{[n]}, \nonumber \\ \theta (x)&=\theta _j+ \zeta \ell \theta _j'+{(\zeta ^2\ell ^2/2!)} \theta _j''+ \ldots +{(\zeta ^n \ell ^n/n!)} \theta _j^{[n]}, \nonumber \\ q(x)&=q_j+\zeta \ell q_j'+{(\zeta ^2\ell ^2/2!)} q_j''+ \ldots +{(\zeta ^m\ell ^m/m!)} q_j^{[m]}. \end{aligned}$$

(133)

Here $v_j^{[n]}$, etc., is an abbreviation for $d^n v(x_j)/dx^n$. [Note: Brackets are used instead of parentheses to avoid confusion with element superscripts. If derivatives are indexed by primes or roman numerals, brackets are omitted.] Evaluate the v(x) and $\theta (x)$ series at i and k by setting $\zeta =\pm 1$, and insert in (132). To compute $f_j$ and $m_j$ use the q(x) series evaluated over elements (1) and (2). Denote by $\xi ^{(1)}$ and $\xi ^{(2)}$ the isoparametric coordinates over elements (1) and (2), respectively, linked to $\zeta$ by $\zeta ^{(1)}=-{\textstyle {1\over 2}}(1-\xi ^{(1)})$ and $\zeta ^{(2)}={\textstyle {1\over 2}}(1+\xi ^{(2)})$. The lateral loads over (1) and (2) are $q^{(1)}=q(\zeta ^{(1)})$, and $q^{(2)}=q(\zeta ^{(2)})$, respectively. Then

$$\begin{aligned} f_j&={\ell \over 2}\left( \int _{-1}^1 N_3^{(1)} q^{(1)} d{\xi ^{(1)}}+\int _{-1}^1 N_1^{(2)} q^{(2)} d{\xi ^{(2)}}\right) , \nonumber \\ m_j&={\ell \over 2}\left( \int _{-1}^1 N_4^{(1)} q^{(1)} d{\xi ^{(1)}}+\int _{-1}^1 N_2^{(2)} q^{(2)} d{\xi ^{(2)}}\right) . \end{aligned}$$

(134)

Here

$$\begin{aligned} N_3^{(1)}&={\textstyle {1\over 4}}(1+{\xi ^{(1)}})^2 (2-{\xi ^{(1)}}), \quad N_4^{(1)}=-{\textstyle {1\over 8}}\ell (1+{\xi ^{(1)}})^2 (1-{\xi ^{(1)}}), \nonumber \\ N_1^{(2)}&= {\textstyle {1\over 4}}(1-{\xi ^{(2)}})^2 (2+{\xi ^{(2)}}), \quad N_2^{(2)}= {\textstyle {1\over 8}}\ell (1-{\xi ^{(2)}})^2 (1+{\xi ^{(2)}}), \end{aligned}$$

(135)

are the BE Hermitian shape functions associated with the j node trial function. (For the Timoshenko beam model those functions must be adjusted; see Sect. 2.7). To show the resulting system in matrix form it is convenient to collect values and derivatives at node j into column vectors:

$$\begin{aligned} \mathbf {v}_j=\left[ \begin{array}{cccc}v_j&v_j'&v_j''&\ldots v_j^{[n]} \end{array}\right] ^T, \quad \varvec{\theta }_j=\left[ \begin{array}{cccc}\theta _j&\theta _j'&\theta _j''&\ldots \theta _j^{[n]} \end{array}\right] ^T, \quad \mathbf {q}_j=\left[ \begin{array}{cccc}q_j&q_j'&q_j''&\ldots q_j^{[m]} \end{array}\right] ^T. \end{aligned}$$

(136)

Both $\mathbf {v}_j$ and $\varvec{\theta }_j$ have length $n+1$ whereas $\mathbf {q}_j$ has length $m+1$. The resulting differential system can be compactly written in matrix form:

$$\begin{aligned} \left[ \begin{array}{cccc}\mathbf {R}_{vv} &{} \mathbf {R}_{v\theta }\\ \mathbf {R}_{\theta v}&{} \mathbf {R}_{\theta \theta }\end{array}\right] \left[ \begin{array}{cccc}\mathbf {v}_j\\ \varvec{\theta }_j\end{array}\right] = \left[ \begin{array}{cccc}\mathbf {L}_v\\ \mathbf {L}_\theta \end{array}\right] \mathbf {q}_j. \end{aligned}$$

(137)

Here $\mathbf {R}_{vv}$, $\mathbf {R}_{v\theta }$, $\mathbf {R}_{\theta v}$ and $\mathbf {R}_{\theta \theta }$ are upper triangular Toeplitz matrices of order $(n{+}1)\times (n{+}1)$ whereas $\mathbf {L}_v$ and $\mathbf {L}_\theta$ are generally rectangular matrices of order $(n{+}1)\times (m{+}1)$. Here are these matrices for $n=6$, $m=2$:

$$\begin{aligned} \mathbf {R}_{vv}&={EI \beta \over \ell } \left[ \begin{array}{ccccccc} 0 &{} 0 &{} -\,12 &{} 0 &{} -\,\ell ^2 &{} 0 &{} -\,\ell ^4/30 \\ 0 &{} 0 &{} 0 &{} -\,12 &{} 0 &{} -\,\ell ^2 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} -12 &{} 0 &{} -\ell ^2 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -12 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -12 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array}\right] , \nonumber \\ \mathbf {R}_{\theta v}&={EI \beta \over \ell }\left[ \begin{array}{ccccccc} 0 &{} -\,12 &{} 0 &{} -\,2\ell ^2 &{} 0 &{} -\,\ell ^4/10 &{} 0 \\ 0 &{} 0 &{} -\,12 &{} 0 &{} -\,2\ell ^2 &{} 0 &{} -\,\ell ^4/10 \\ 0 &{} 0 &{} 0 &{} -12 &{} 0 &{} -2\ell ^2 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} -12 &{} 0 &{} -2\ell ^2 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -12 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -12 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array}\right] , \nonumber \\ \mathbf {R}_{\theta \theta }&={EI\over \ell }\left[ \begin{array}{ccccccc} 12\beta &{} 0 &{} \ell ^2\phi &{} 0 &{} \ell ^4\phi /12 &{} 0 &{} \ell ^6\phi /360 \\ 0 &{} 12\beta &{} 0 &{} \ell ^2\phi &{} 0 &{} \ell ^4\phi /12 &{} 0 \\ 0 &{} 0 &{} 12\beta &{} 0 &{} \ell ^2\phi &{} 0 &{} \ell ^4\phi /12 \\ 0 &{} 0 &{} 0 &{} 12\beta &{} 0 &{} \ell ^2\phi &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 12\beta &{} 0 &{} \ell ^2\phi \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 12\beta &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 12\beta \end{array}\right] ,\nonumber \\ \mathbf {R}_{v\theta }&=-\mathbf {R}_{\theta v}, \quad \mathbf {L}_f = {\ell \over 15} \left[ \begin{array}{ccc} 15 &{} 0 &{} \ell ^2 \\ 0 &{} 15 &{} 0 \\ 0 &{} 0 &{} 15 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \end{array}\right] , \quad \mathbf {L}_m= {\ell ^3\over 15} \left[ \begin{array}{ccc} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \end{array}\right] , \quad \end{aligned}$$

(138)

In the $\mathbf {R}_{\theta \theta }$ expression, $\phi =3\beta -1$. Elimination by static condensation of $\varvec{\theta }_j$ gives the system

$$\begin{aligned} \mathbf {R}\mathbf {v}_j = \mathbf {L}\mathbf {q}_j, \quad \hbox {with}\quad \mathbf {R}= \mathbf {R}_{vv}^{}-\mathbf {R}_{v\theta }^{} \mathbf {R}_{\theta \theta }^{-1} \mathbf {R}_{\theta v}^{}, \quad \hbox {and}\quad \mathbf {L}= \mathbf {L}_{v}^{}-\mathbf {R}_{v\theta }^{} \mathbf {R}_{\theta \theta }^{-1} \mathbf {L}_{\theta }^{}. \end{aligned}$$

(139)

For the matrices (138) we get

$$\begin{aligned} \mathbf {R}= EI \ell \left[ \begin{array}{ccccccc} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} -\,\ell ^2 (\beta -1)/(12\beta )\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array}\right] ,\quad \mathbf {L}= \ell \left[ \begin{array}{cccc} 1 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \end{array}\right] \end{aligned}$$

(140)

Postmultiplication by $\mathbf {v}_j$ and $\mathbf {q}_j$ yields

$$\begin{aligned} \mathbf {R}\mathbf {v}_j&=\left[ \begin{array}{cccccc}12\beta v''''(x)-\ell ^2(\beta -1) v^{[vi]}(x)&v^{[v]}(x)&v^{[vi]}&0&0&0\end{array}\right] ^T, \nonumber \\ \mathbf {L}\mathbf {q}_j&=\left[ \begin{array}{ccccccc} q(x)&q'(x)&q''(x)&0&0&0&0 \end{array}\right] ^T \end{aligned}$$

(141)

in which x-derivatives of order 5 or higher are denoted with bracketed Roman numeral superscripts. The nontrivial equations in $\mathbf {R}\mathbf {v}_j = \mathbf {L}\mathbf {q}_j$ are $EI \,\, v''''(x)=q(x)-(1-\beta ) \ell ^4 q''''(x)/(720 \beta )$ and $EI v^v=q'(x)$. Solving for $v''''(x)$, which is uncoupled, gives the finite order ODE

$$\begin{aligned} EI \, v''''(x)=q(x)- {1-\beta \over \beta } { \ell ^2 q''(x)\over 12 \,EI}. \end{aligned}$$

(142)

If n and m are increased in steps of 2 so that $m=n-4$, additional terms appear. For example:

$$\begin{aligned} n=8, m=4: \quad EI v''''(x)&=q(x)- {1-\beta \over \beta } \left( {\ell ^2 q''(x)\over 12} -{\ell ^4 q''''(x)\over 720} \right) , \nonumber \\ n=10, m=6: \quad EI v''''(x)&=q(x)- {1-\beta \over \beta } \left( {\ell ^2 q''(x)\over 12} -{\ell ^4 q''''(x)\over 720}+{\ell ^6 q^{(vi)}(x)\over 30240}\right) . \end{aligned}$$

(143)

The infinite order ODE (IOMoDE) that emerges as $n\rightarrow \infty$ is (108). As noted in Appendix Sect. 3.3, this reduces to the exact governing ODE (107) if $\beta =1$, thus validating (30) as the nodally exact stiffness.

1.2 4.2 MoDE Implementation in Mathematica

The Mathematica coding of the MoDE analysis flowcharted in Fig. 15 is presented here. It can do material and geometric stiffness matrices that cover both Timoshenko and BE models.

Implementation assumption: the two elements in the patch of Fig. 16b must be identical, with same EI, $\ell$ and $\Phi$. Relaxing this constraint, which is sufficient for our purposes, would need changes to the logic of modules PatchNodeForces and BeamFOMoDE.

Fig. 22

Group listing of primitive modules

Five primitive (= self-contained) modules are group-listed in Fig. 22. (Actually one: TaylorSeries, is technically an inline function.) They are invoked as follows:

$$\begin{aligned} {\texttt {fs = TaylorSeries[f,x,h,n];}} \end{aligned}$$

(144)

This function returns the Taylor series of f[x] about x=0 with increment h, up to and including the $n^{th}$ derivative, where $n\ge 0$. For example, fseries = TaylorSeries[f,x,h,3] returns f[x]+h * f[x]+h² * f[x]/2+h³ * f[x])/6 in fseries.

$$\begin{aligned}&{\texttt {G = BeamGMatrix[Le},\Phi ,\chi \texttt {];}} \nonumber \\&{\texttt {H = BeamHMatrix[Le},\Phi ,\chi \texttt {];}} \end{aligned}$$

(145)

Module BeamGMatrix returns the $\mathbf {G}$ matrix defined in (77) that maps physical to natural DOF whereas module BeamHMatrix returns its inverse $\mathbf {H}=\mathbf {G}^{-1}$ thap maps natural to physical DOF. Both are $4\times 4$ matrices. The arguments are the element length Le (called $\ell$ in the text), the shear flexibility coefficient $\Phi$ and the adjusted Legendre cubic polynomial parameter $\chi$.

$$\begin{aligned} {\texttt {SG = BeamSGMatrix[tag,Le},\Phi ,F,\beta \texttt {pars];}} \end{aligned}$$

(146)

This module forms the natural geometric stiffness template $\mathbf {S}_G$. The first argument, tag, is a character string. Legal tags are “KG4” or “KG6”. If “KG4”, it builds the template matrix (106). using the element length Le, the axial force F and the four $\beta$ parameters supplied in the list argument $\beta$pars. If tag is “KG6” it fills the bottom $3\times 3$ principal minor of that matrix with two additional parameters: $\beta _5$ and $\beta _6$—this particular template is not used in the present work. Argument $\Phi$ is not used in the listed implementation; it is provided for possible extensions. The module returns $\mathbf {S}_G$ as a $4\times 4$ matrix. If tag is illegal it returns Null.

$$\begin{aligned} {\texttt {Nf = BeamShapeFunctions[model,Le},\Phi ,\xi \texttt {];}} \end{aligned}$$

(147)

This module returns polynomial shape functions for beam lateral displacements, as listed in Sect. 2.7. It is used to compute consistent loads as per (134) when doing a material stiffness template. Argument model is a character string: either “BE” or “Timo”. Arguments Le and $\Phi$ are as described above. Argument $\xi$ is the isoparametric natural coordinate; this is normally a symbol.

If model is “BE” it returns the four classical Hermitian shape functions, in which case $\Phi$ is not used. If model is “Timo” it returns the 4 cubic shape functions listed in (20). The module returns a 4-item list. If model is illegal it returns a list of four Null values.

Fig. 23

Modules to generate element stiffness matrix and two-element patch stiffness

Two modules that construct a beam element stiffness matrix and assemble the stiffness of a two-element patch are listed in Fig. 23. They are invoked as follows.

$$\begin{aligned} {\texttt {Ke = BeamStiffMatrix[tag,EI,Le},\Phi ,\texttt {F},\beta \texttt {pars},\chi \texttt {];}} \end{aligned}$$

(148)

This module forms the stiffness of a beam element. The first argument, tag, is a character string that can be either “KM” to form the material stiffness matrix, or “KG4” to form the geometric stiffness matrix with the 4-parameter template $\mathbf {S}_G$ of (106) and the $\chi$-parametrized matrix $\mathbf {H}$ of (77).

The other arguments are: the bending rigidity EI (constant over element), the element length Le, the shear flexibility coefficient $\Phi$, the axial force F, positive if tension (used only for a geometric stiffness matrix), the template parameter list $\beta$pars, and the Legendre basis parameter $\chi$.

If tag is “KM”, only one parameter is needed, so $\beta$pars is simply $\beta$ , and $\chi$ is not used since $\mathbf {K}_M$ is formed directly. If tag is “KG4” the module calls BeamSGMatrix and BeamHMatrix to get $\mathbf {K}_G$ as $\mathbf {H}^T\mathbf {S}_G\mathbf {H}$. The second to last arguments may be symbolic or numeric. The module returns a $4\times 4$ stiffness matrix. If tag is illegal an error message is printed, and the module returns Null.

$$\begin{aligned} {\texttt {Kp = TwoElemPatchStiff[Ke1,Ke2];}} \end{aligned}$$

(149)

This module assembles the patch stiffness equations of a patch of two adjacent elements, as pictured in Fig. 16. The inputs are their $4\times 4$ stiffness matrices computed via BeamStiffMatrix. In the present work the two elements are assumed identical, so Ke1=Ke2.

The assembled patch stiffness matrix is $6\times 6$ as it involves three nodes with two DOF each, cf. Fig. 16. But rows pertaining to i and k are not needed in the ensuing analysis because their DOF are eliminated via the series (133). Thus the first two and the last two rows can be removed. The module actually returns a $2 \times 6$ matrix Kp, such as the one shown in the left hand side of (132).

Fig. 24

Patch nodal forces computation module

Module PatchNodeForces, listed in Figs. 24, 25, constructs the matrices $\mathbf {L}_f$, $\mathbf {L}_m$ and vector $\mathbf {q}_f$ defined in (155) and (156), which are used for the matrix computation of the patch consistent node forces $f_j$ and $m_j$ of (151). These are used only for material stiffness analysis. The module is invoked as

$$\{\text{Lf},\text{Lm},\text{qj}\} = PatchNodeForces[model,q,Le,\Phi ,\xi ,x,n,m];$$

(150)

Here model is either “BE” or “Timo”, q is a symbol for the lateral load q(x), Le, $\Phi$ and $\xi$ are as defined before, and x is the symbol for the axial coordinate used in the Taylor series expansion of q(x). Finally, integers n and m specify the highest derivative order used in the Taylor expansions of v(x), $\theta (x)$, and q(x) as indicated in (136).

The module returns the list Lf,Lm,qj. Here the matrices Lf and Lm are defined in (137) whereas vector qj is defined in (136). In geometric stiffness template analysis it is usual to set q to 0, in which case null matrices are returned and the consistent force computation skipped.

Fig. 25

Finite order modified differential equation (FOMoDE) module

Module BeamFOMoDE, listed in Fig. 28, carries out the computation of the Finite Order Modified ODE. The module is invoked as

$$\text{FOMoDE, RHS, err = BeamFOMoDE}[ \text{model, MorG, EI, Le}, \Phi, {\text{F, v}}, \uptheta,\upxi, {\text{x, n, nh, m, Kp, q, simp, inv}}]; $$

(151)

Several arguments have been defined before: model, EI, $\Phi$, Le, F, $\xi$, x, n, m, and q. Additional ones are

MorG:

A letter string: “M” for material stiffness analysis, “G” for geometric stiffness analysis.

Kp:

The $2\times 6$ patch stiffness matrix returned by TwoElemPatchStiff.

nh:

In geometric stiffness analysis of Timo model, highest power of h=Le in a series expansion of the $k^2$ term of (124). Normally set to 6. See RHS below.

simp:

A logical flag. If True, additional simplifications are done internally.

inv:

A logical flag. If True, the Inverse function is used in the condensation process (139). If False, the LinearSolve function is used.

The module returns a 3-item list: FOMoDE.RHS,err:

FOMoDE :

The computed FOMoDE, represented as the value of $v''''(x)$. More precisely, its series expansion in even powers of h=Le up to order O(nh).

RHS :

The right hand side of the canonical governing equation. For material stiffness analysis it is the RHS of (109). For geometric stiffness analysis it is the series expansion of the coefficient $k^2$ in (124) in even powers of h=Le, up to order O(nh). Note that if $\Phi =0$ (the BE model), $k^2$ is simply F/EI and no expansion is needed.

err :

The difference FOMoDE-RHS. If zero, the correct matching is achieved. If nonzero, further analysis is done by examining its Taylor series in the element length.

Fig. 26

Driver script to analyze the one-parameter material stiffness matrix for the BE model

Fig. 27

Driver script to analyze the oparameter geometric stiffness matrix for the BE model using the “KG4” template. Four free parameters: $\beta _2$, $\beta _3$, $\beta _4$ and $\chi$ ($\beta _1$ is preset to 1)

Fig. 28

Module containing the inlined results for the IOMoDE errror coefficients $C_2$, $C_4$ and $C_6$, as output by the script of Fig. 27

1.3 4.3 MoDEAnalysis Examples

The script of Fig. 26 does the MoDE analysis of the one-parameter material stiffness matrix (104). It uses $n=10$ and $m=6$. The printed results (not shown) corroborate the result (108).

The script of Fig. 27 is the first step in finding an optimal geometric stiffness for the BE model. It uses the “KG4” template with $\beta _1$ preset to 1, three free $\beta$ parameters: $\beta _2$, $\beta _3$ and $\beta _4$, plus the Legendre basis parameter $\chi$.

The primary result of the MoDE analysis are the error coefficients $C_2$, $C_4$, and $C_6$, which are complicated functions of the four parameters. These are printed in InputForm format and inlined to flesh out the TruncBECoeffs module listed in Fig. ??. The subsequent parameter-optimization analysis is done with ad-hoc scripts (not shown) that call this module to retrieve the coefficients and perform the algebraic manipulations described in Appendix Sect. 3.5.

Appendix 5: Buckling Exact Stiffness

This Appendix develops exact stiffness equations for buckling analysis with plane Timoshenko beam elements. This is possible when exact homogeneous solutions of the governing stability equations are known. The case treated below is that of a prismatic element subjected to constant compressive axial force, under linearized prebuckling (LPB) assumptions. The resulting stiffness matrix will be labeled buckling exact (BuE) and tagged with an asterisk: $()^*$.

The homogeneous solutions of the equilibrium equation (11) are trigonometric functions (sine and cosine) if $k^2$ is real positive and constant over the element, and so is the dimensionless parameter $\lambda =k\ell$. The solutions are nonlinear in $\lambda$, and so are the resulting stiffness matrix entries. These complicated expressions have limited practical use as explained in Appendix Sect. 5.3, but serve to corroborate previous results when expanded in Taylor or Padé series.

The beam-column element studied here is pictured in Fig. 2. Of the six degrees of freedom (DOF) shown there only four are used in the BuE stiffness equations derived below:

$$\begin{aligned} \mathbf {u}=\left[ \begin{array}{cccc}v_1&\theta _1&v_2&\theta _2\end{array}\right] ^T, \end{aligned}$$

(152)

The axial displacements $u_1$ and $u_2$, which uncouple at the individual element level, may be included later to complete a beam-column element, able to transform to global coordinates for assembly.

As in previous Appendices, the element supercript: $()^e$ is omitted for brevity.

1.1 5.1 Buckling Exact Physical Stiffness

The axial force is denoted $P=-F$, which is positive if compressive. Under the LPB assumptions, the stiffness entries depend nonlinearly on P whereas the generic $4\times 4$ matrix form (103) still holds. The template method is used to construct $\mathbf {K}^*$. The first step is to introduce dimensionless free parameters $C_{ij}$:

$$\begin{aligned} \mathbf {K}^*={EI\over \ell ^3} \left[ \begin{array}{cccc}C_{11} &{} C_{12} \ell &{} C_{13} &{} C_{14} \ell \\ C_{12} \ell &{} C_{22} \ell ^2 &{} C_{23} \ell &{} C_{24} \ell ^2 \\ C_{13} &{} C_{23} \ell &{} C_{33} &{} C_{34} \ell \\ C_{14} \ell &{} C_{24} \ell ^2 &{} C_{34} \ell &{} C_{44} \ell ^2 \end{array}\right] \end{aligned}$$

(153)

This is called a universal template, since it include all possible elements with a symmetric stiffness. For a prismatic beam element, additional symmetries require $C_{11}=C_{33}$, $C_{22}=C_{44}$, etc. Let $\mathbf {v}_{T}=\left[ \begin{array}{cccc} 1&0&1&0\end{array}\right]$ be the translational rigid body mode. Then $\mathbf {K}^* \mathbf {v}_T=\mathbf {0}$, the zero force 4-vector. Imposing these constraints the number of free parameters is reduced to four, renamed $c_1$ through $c_4$:

$$\begin{aligned} \mathbf {K}^*={EI\over \ell ^3} \left[ \begin{array}{cccc}c_{1} &{} c_{3} \ell &{} -\,c_{1} &{} c_{3} \ell \\ c_{3} \ell &{} c_{2} \ell ^2 &{} -\,c_{3} \ell &{} c_{4} \ell ^2 \\ -c_{1} &{} -\,c_{3} \ell &{} c_{1} &{} -\,c_{3} \ell \\ c_{3} \ell &{} c_{4} \ell ^2 &{} -\,c_{3} \ell &{} c_{2} \ell ^2 \end{array}\right] \end{aligned}$$

(154)

Next consider the rigid “element tilt” $\mathbf {v}_{R}=\left[ \begin{array}{cccc} -\ell /2&1&\ell /2&1\end{array}\right]$. Now $\mathbf {K}^* \mathbf {v}_R=\left[ \begin{array}{cccc}-P&0&P&0\end{array}\right]$, which is consistent with the LPB assumptions. This supplies two additional constraints:

$$\begin{aligned} c_1+2c_3=P L^2/EI=\lambda ^2(1-\Phi \lambda ^2/12), \quad c_2-c_3+c_4=0, \end{aligned}$$

(155)

where the dimensionless axial force parameter $\lambda$ is defined as in (124), except that L becomes $\ell$:

$$\begin{aligned} \lambda =k \ell , \hbox { with } k^2={P\over EI \left( 1-{\displaystyle P\over \displaystyle GA_s}\right) }= {P\over EI-{\textstyle {1\over 12}}\Phi P \ell ^2}. \end{aligned}$$

(156)

As noted above, we assume that $k^2\ge 0$, so the $+$ square root is taken. Solving (155) for $c_3$ and $c_4$ gives

$$\begin{aligned} c_3 = {\textstyle {1\over 2}}c_1+{\lambda ^2/2 \over 1-\lambda ^2\Phi /12}, \quad c_4 = {\textstyle {1\over 2}}c_1-c_2+{\lambda ^2/2 \over 1-\lambda ^2\Phi /12}, \end{aligned}$$

(157)

whence the number of independent template parameters is cut to two: $c_1$ and $c_2$. Matching the homogeneous solutions given above to the nodal DOF leads to two trascendental nonlinear equations, which are solved with Mathematica, The solution can be presented as

$$\begin{aligned} \psi&=1/(1+\lambda ^2\Phi /12), \quad \eta =2-2\cos \lambda -\psi \lambda \sin \lambda , \nonumber \\ c_1&=\psi ^2\lambda ^3 \sin \lambda /\eta , \quad c_2=\lambda (\sin \lambda -\psi \lambda \cos \lambda )/\eta , \nonumber \\ c_3&=\psi \lambda ^2 (1-\cos \lambda )/\eta , \quad c_4=\lambda (\psi \lambda -\sin \lambda )/\eta , \end{aligned}$$

(158)

in which all Greek symbols identify dimensionless variables. If the node displacements are precluded: $v_1=v_2=0$, removal of rows and columns 1 and 3 reduces (154) to the moment-rotation stiffness matrix:

$$\begin{aligned} \mathbf {K}^*_R={EI\over \ell } \left[ \begin{array}{cccc} c_{2} &{} c_{4} \\ c_{4} &{} c_{2} \end{array}\right] = {EI \lambda \over \eta \ell }\left[ \begin{array}{cccc} \sin \lambda -\psi \lambda \cos \lambda &{} \psi \lambda -\sin \lambda \\ \psi \lambda -\sin \lambda &{} \sin \lambda -\psi \lambda \cos \lambda ) \end{array}\right] . \end{aligned}$$

(159)

As noted by Bazant and Cedolin [8, Sec. 2.1]. for the BE model ($\Phi =0$) these equations appeared in a 1935 NACA note [50], where it was used for the moment-distribution, hand-computation method of Hardy Cross [20]. Inversion provides the moment-rotation flexibility matrix.

$$\begin{aligned} \mathbf {F}^*_R=(\mathbf {K}^*_R)^{-1}= \lambda ^2\left[ \begin{array}{cccc} 1-\lambda \cot \lambda +\lambda ^2 \Phi /12 &{} 1-\lambda \csc \lambda +\lambda ^2 \Phi /12\\ 1-\lambda \csc \lambda +\lambda ^2 \Phi /12 &{} 1-\lambda \cot \lambda +\lambda ^2 \Phi /12\end{array}\right] . \end{aligned}$$

(160)

The entries of this matrix are known as stability functions in the older literature. For the BE model these can be traced back to von Mises [97]; see also Timoshenko and Gere [91, Ch. 1]. Additional references are [43, 54]. Stability functions for the Timoshenko model were first derived by Absi [2], and are further developed in [7, 81, 82, 100, 101].

Note that geometric and material properties are strongly coupled in $\mathbf {K}^*$. No segregation into material and geometric stiffness is possible. That separation appears, however, when expanding in Taylor series, as done in Appendix Sect. 5.4 below.

1.2 5.2 Buckling Exact Natural Stiffness

In Appendix 2 we started with the natural stiffness matrix: $\mathbf {S}_M$ for material and $\mathbf {S}_G$ for geometric, and passed to the physical stiffness through congruential transformations with the $\mathbf {H}$ of (77). Here we follow the reverse process: starting with the BuE physical matrix $\mathbf {K}^*$, we pass to the natural stiffness via $\mathbf {S}^*=\mathbf {G}^T \mathbf {K}^* \mathbf {G}$, where $\mathbf {G}=\mathbf {H}^{-1}$ is also given in (77). It is convenient to take $\chi =5$ when forming $\mathbf {G}$ to get a diagonal natural stiffness:

$$\begin{aligned} \mathbf {S}^*=(\left. \mathbf {G}\right| _{\chi =5})^T \mathbf {K}^* \left. \mathbf {G}\right| _{\chi =5}= \mathbf{diag }[0, S_{22}, S_{33}, S_{44}], \end{aligned}$$

(161)

where

$$\begin{aligned} S_{22}&={-48 EI \lambda ^2\over \ell ^3 (12+\lambda ^2 \Phi )}, \quad S_{33} ={72 EI \lambda \cot (\lambda /2)\over \ell ^3}, \nonumber \\ S_{44}&= {2400 EI \lambda ^2 (1{+}\Phi )^2 \sin (\lambda /2)^2\over \ell ^3 (12+\lambda ^2 \Phi - (12+\lambda ^2 \Phi ) \cos \lambda -6 \lambda \sin \lambda )} \end{aligned}$$

(162)

The Taylor $\lambda$-expansion of this matrix centered at $\lambda =0$ is

$$\begin{aligned} \mathbf {S}^*= \mathbf {S}^*_0+\lambda ^2\mathbf {S}^*_2+ \ldots = {EI\over \ell ^3}\big (\mathbf{diag }[0, 0, 144, 1200 (1{+}\Phi )] -4\lambda ^2 \mathbf{diag }[0, 1, 3, 5]+\ldots \big ) \end{aligned}$$

(163)

Applying a congruential transformation with $\left. \mathbf {H}\right| _{\chi =5}$ to this series reproduces the physical Taylor expansion given in (166). These results are new.

Fig. 29

Buckling of pinned-pinned beam column (a.k.a. Euler column): symmetric buckling: $P_{cr}=n^2 \pi ^2 EI/L^2$, $n=1,3,\ldots$; antisymmetric buckling: $P_{cr}=n^2 \pi ^2 EI/L^2$, $n=2,4,\ldots$

1.3 5.3 Fragility of the Buckling Exact Stiffness

Two advantages are quoted in the FEM literature for the buckling exact stiffness:

• One element per member is sufficient (assuming a prismatic member fabrication)

• All critical loads are captured by solving the resulting eigenproblem.

In summary, a multielement-per-member FEM discretization is avoided while keeping the DOF count small. The following simple counterexample shows that statement to be questionable. Consider the classical Euler beam-column shown in Fig. 29, with one BuE element dscretization. Introducing the node displacement BC: $v_1=v_2=0$ the stiffness matrix of the member reduces to $\mathbf {K}^*_R$ of (159). Its determinant is

$$\begin{aligned} d_R=\det (\mathbf {K}_R)=\left( {EI\over \ell }\right) ^2 (c_2^2-c_4^2)= \left( {EI\over \ell }\right) ^2 {\lambda ^3\over 2 \tan (\lambda /2)-\lambda } \end{aligned}$$

(164)

where the last expression is obtained by Mathematica’s FullSimplify. Now $d_R=0$ for $\lambda =\pi ,3\pi ,\ldots$, which provide the symmetric buckling modes, as expected. However, $d_R$ does not vanish for $\lambda =2\pi ,4\pi ,\ldots$. These are associated with axisymmetric buckling nodes, which are therefore lost. To recover these modes it is necessary to supply one more constraint: $\theta _1=\theta _2$. But in a complicated structure such post-facto corrections are unfeasible since they would depend on a priori knowledge of unknown buckling eigenmodes.

An additional complication in frame analysis is the possible occurrence of positive axial force F, that is, negative P, in some members. In that case the ordinary since and cosine functions must be replaced by their hyperbolic counterparts, which further complicates programming.

Singularities bring additional problems in purely numerical work. First, if $P\rightarrow 0$, $\lambda \rightarrow 0$ and all entries become 0/0. Second, the denominator $\eta$ vanishes for other nonzero values of $\lambda$ whence all entries become $\infty$. This can destroy the convergence of solution methods (e.g. Newton’s) for the transcendental eigenvalue problem.

1.4 5.4 Series Expansion

The aforementioned defects can be cured by taking Taylor series in $\lambda$, with the tradeoff of accepting approximations and the consequent need to put multiple elements per member. Expanding up to $\lambda ^4$ with Mathematica gives

$$\begin{aligned} c_1&={12\over 1+\Phi }- \lambda ^2 {6+10\Phi +5\Phi ^2 \over 5(1+\Phi )^2}+ \lambda ^4 {-3+165 \Phi +525 \Phi ^2+525 \Phi ^3+175 \Phi ^4\over 2100 (1+\Phi )^3}+ O(\lambda ^6), \nonumber \\ c_2&={4+\Phi \over 1+\Phi }-\lambda ^2 {8+10\Phi +5\Phi ^2 \over 60(1+\Phi )^2}+ \lambda ^4 {44+135 \Phi +105 \Phi ^2+35 \Phi ^3\over 25200 (1+\Phi )^3} +O(\lambda ^6), \nonumber \\ c_3&={6\over 1+\Phi }+\lambda ^2 {1 \over 10(1+\Phi )^2}+ \lambda ^4 {3+10 \Phi \over 4200 (1+\Phi )^3}+O(\lambda ^6), \nonumber \\ c_4&={2-\Phi \over 1+\Phi }+\lambda ^2 {2+10\Phi +5\Phi ^2 \over 60(1+\Phi )^2}+ \lambda ^4 {26+75 \Phi +105 \Phi ^2+35 \Phi ^3\over 25200 (1+\Phi )^3} +O(\lambda ^6) , \end{aligned}$$

(165)

Insertion into (154) gives the element stiffness matrix expansion in $\lambda$:

$$\begin{aligned} \mathbf {K}^*=\mathbf {K}_0+ \lambda ^2 \mathbf {K}_2+\lambda ^4\mathbf {K}_4+\ldots \end{aligned}$$

(166)

The first two terms have been encountered before. $\mathbf {K}_0$ is the material stiffness matrix given in (29), denoted by $\mathbf {K}_M$ there. $\mathbf {K}_2$ is the geometric stiffness given in (36), where it is called $\mathbf {K}_{Gv}$.

The names initial stress matrix and load-geometric stiffness were used for $\mathbf {K}_2$ before 1970. The geometric qualifier refers to an accidental feature of the BE model: that term is independent of material properties, since if $\Phi =0$, $\lambda ^2$ is simply P/EI, which cancels out EI. That is not the case, however, for the Timoshenko model, since those properties enter indirectly through $\Phi$.

$\mathbf {K}_4$ and following expansion terms are higher order stiffness matrices containing both material and geometric properties. These are rarely used for FEM buckling analysis since they lead to nonlinear eigenproblems.

Appendix 6: Fitting Element Fields

This Appendix looks at the problem of Polynomial Field Fitting over a finite element domain. This is a subset of the more general curve- and surface-fitting problem, as covered for example in [52]. The fitting methods described in Appendix Sect. 6.1 are applicable to one through three spatial dimensions, but only the 1D case is worked out in Appendix Sect. 6.2. The results listed there are used in Sect. 5.4 for fitting beam shape functions to the DOF configurations of Fig. 6. The problem is defined as follows.

• (I) A polynomial interpolation of a scalar function f, called the base function, is defined over a finite element domain using n nodal values of f and n shape functions expressed in natural coordinates.

• (II) A polynomial interpolation of a scalar function g, called the fitted function, is defined over the same element domain in terms of $m\le n$ nodal values and m shape functions. This function approximates the base function f in a least square sense.

The interpolations in (I) and (II) will be called the base and fitted interpolants of f and g, respectively.

The case $m=n$ is not excluded, as it covers situations in which the same polynomial order is defined from different nodal values. A well known 2D example is the linear triangle with $n=3$ corner values (the Turner triangle) versus that with $m=3$ midpoint values (the Veubeke triangle).

Passing from f to g involves least square fitting. This is an operation that has been extensively studied in the applied mathematics literature. However, most of the work available there does not consider the use of natural coordinates as well as integration over the element domain. These two features are essential for finite element applications.

Passing from g to f is a collocation operation, which simply involves evaluating g at the f nodes.

1.1 6.1 Least-Square Fitting

Let the node values for the base and reduced interpolant be arranged in column vectors denoted by $\mathbf {f}$ and $\mathbf {g}$ of lengths n and m, respectively. The associated shape functions are arranged in row vectors $\mathbf {N}_f$ and $\mathbf {N}_g$, respectively. The shape functions of both f and g will be assumed to be linearly independent. The interpolations may be expressed as the dot products

$$\begin{aligned} f =\mathbf {N}_f \mathbf {f}, \qquad g=\mathbf {N}_g \mathbf {g}. \end{aligned}$$

(167)

The distance between g and f is $d=g-f=\mathbf {N}_g \mathbf {g}-\mathbf {N}_f \mathbf {f}$. Its square is

$$\begin{aligned} d^2 =(\mathbf {N}_g \mathbf {g}-\mathbf {N}_f \mathbf {f})^T (\mathbf {N}_g \mathbf {g}-\mathbf {N}_f \mathbf {f})= \mathbf {g}^T \mathbf {N}_g^T \mathbf {N}_g\mathbf {g}-2 \mathbf {g}^T \mathbf {N}_g^T \mathbf {N}_f \mathbf {f}+ \mathbf {f}^T \mathbf {N}_f^T \mathbf {N}_f\mathbf {f}. \end{aligned}$$

(168)

As the distance norm we take

$$\begin{aligned} s={\textstyle {1\over 2}}\int _\Omega d^2 d\Omega = {\textstyle {1\over 2}}\left( \mathbf {f}^T \mathbf {A}\mathbf {f}-2 \mathbf {g}^T \mathbf {B}\mathbf {f}+ \mathbf {g}^T \mathbf {C}\mathbf {g}\right) = {\textstyle {1\over 2}}\left[ \begin{array}{cccc}\mathbf {f}^T&\mathbf {g}^T\end{array}\right] \left[ \begin{array}{cccc} \mathbf {A}&{}-\mathbf {B}^T \\ -\mathbf {B}&{}\mathbf {C}\end{array}\right] \left[ \begin{array}{cccc}\mathbf {f}\\ \mathbf {g}\end{array}\right] . \end{aligned}$$

(169)

in which $\Omega$ denotes the element domain, and

$$\begin{aligned} \mathbf {A}=\int _\Omega \mathbf {N}_f^T \mathbf {N}_f \,\, d\Omega , \quad \mathbf {B}=\int _\Omega \mathbf {N}_g^T \mathbf {N}_f \,\, d\Omega , \quad \mathbf {C}=\int _\Omega \mathbf {N}_g^T \mathbf {N}_g \,\, d\Omega . \end{aligned}$$

(170)

Matrices $\mathbf {A}$ and $\mathbf {C}$ are called the self kernels for f and g, respectively. They are symmetric and positive definite, because the shape functions are assumed to be linearly independent. Matrix $\mathbf {B}$ of order $m\times n$ is the coupled kernel, which is rectangular unless $n=m$. (Readers familiar with computational dynamics may recognize $\mathbf {A}$ and $\mathbf {C}$ as the numerical part of the consistent mass matrices pertaining to the base and fitted elements, respectively.) Minimizing s with respect to $\mathbf {g}$ yields

$$\begin{aligned} {\partial s\over \partial \mathbf {g}}=\mathbf {C}\mathbf {g}-\mathbf {B}\mathbf {f}=\mathbf {0}, \quad \hbox {whence}\quad \mathbf {g}=\mathbf {C}^{-1} \mathbf {B}\mathbf {f}=\mathbf {G}\mathbf {f}. \end{aligned}$$

(171)

$\mathbf {G}=\mathbf {C}^{-1} \mathbf {B}$ will be called the $f\rightarrow g$ fitting matrix. When identification of orders m and n are convenient, the matrix will be denoted $\mathbf {G}_m^n$.

Minimizing s with respect to $\mathbf {f}$ yields

$$\begin{aligned} {\partial s\over \partial \mathbf {f}}=\mathbf {A}\mathbf {f}-\mathbf {B}^T \mathbf {g}=\mathbf {0}, \quad \Leftarrow \quad \mathbf {f}= \mathbf {A}^{-1} \mathbf {B}^T \mathbf {g}=\mathbf {F}\mathbf {f}. \end{aligned}$$

(172)

$\mathbf {F}=\mathbf {A}^{-1} \mathbf {B}$ will be called the $g\rightarrow f$ collocation matrix. It is a generalized inverse of $\mathbf {G}$ since

$$\begin{aligned} \mathbf {G}\mathbf {F}=\mathbf {I}, \qquad \mathbf {F}\mathbf {G}=\mathbf {H}, \end{aligned}$$

(173)

in which $\mathbf {I}$ is the $m\times m$ identity matrix whereas $\mathbf {H}$ is an $n\times n$ projector (idempotent) matrix satisfying $\mathbf {H}=\mathbf {H}^2$. Matrix $\mathbf {H}$, which has rank m, is not generally symmetric.

Replacing the expressions of $\mathbf {G}$ and $\mathbf {F}$ into $\mathbf {G}\mathbf {F}=\mathbf {I}$ yields the property

$$\begin{aligned} \mathbf {C}^{-1}\mathbf {B}\mathbf {A}^{-1} \mathbf {B}^T=\mathbf {I}, \quad \hbox { or } \quad \mathbf {B}^T \mathbf {A}^{-1} \mathbf {B}=\mathbf {C}, \end{aligned}$$

(174)

whereas replacing the expressions of $\mathbf {G}$ and $\mathbf {F}$ into $\mathbf {F}\mathbf {G}=\mathbf {H}$ yields the property

$$\begin{aligned} \mathbf {A}^{-1}\mathbf {B}^T \mathbf {C}^{-1} \mathbf {B}=\mathbf {H}, \quad \hbox { or } \quad \mathbf {B}^T \mathbf {C}^{-1} \mathbf {B}=\mathbf {A}\mathbf {H}. \end{aligned}$$

(175)

Notice that $\mathbf {B}^T \mathbf {A}^{-1} \mathbf {B}$ and $\mathbf {B}^T \mathbf {C}^{-1} \mathbf {B}$ are the Schur complements of $\mathbf {C}$ and $\mathbf {A}$, respectively, in the supermatrix shown in (169). Since $\mathbf {B}^T \mathbf {C}^{-1} \mathbf {B}$ is symmetric, so is $\mathbf {A}\mathbf {H}$. It can be readily shown that the reverse product $\mathbf {H}\mathbf {A}$ is antisymmetric.

Inserting (171) into (169) gives the minimum distance as

$$\begin{aligned} s_{min}={\textstyle {1\over 2}}\mathbf {f}^T (\mathbf {A}-\mathbf {B}^T \mathbf {C}^{-1} \mathbf {B}) \mathbf {f}= {\textstyle {1\over 2}}\mathbf {f}^T (\mathbf {A}-\mathbf {A}\mathbf {H}) \mathbf {f}= {\textstyle {1\over 2}}\mathbf {f}^T \mathbf {A}\hat{\mathbf {H}} \mathbf {f}, \end{aligned}$$

(176)

in which $\hat{\mathbf {H}}=\mathbf {I}-\mathbf {H}$. Since $\mathbf {H}$ is a projection matrix, so is $\hat{\mathbf {H}}$. As a quick check, take $\mathbf {N}_f\equiv \mathbf {N}_g$, whence $m=n$ and $\mathbf {C}=\mathbf {B}=\mathbf {A}$; if so $\mathbf {H}=\mathbf {I}$ and $s_{min}=0$, as can be expected.

1.2 6.2 Line Elements

Line elements are used for one-dimensional FEM models. In structural mechanics: bars, beams, shafts, etc. This Section tabulates results for the seven line elements shown in Fig. 30. The line segment is assumed straight with constant metric and length $\ell$ along x. As a consequence, the integrals over the line can be evaluated exactly, without resource to numerical integration.

Elements labeled “LLn” are the standard Lagrangian ones, with function values at equally spaced n nodes, which include the end points if $n\ge 2$. Elements labeled “LLnG” have internal function values given at $\texttt {n>=2}$ Gauss points. (The one-point line element is denoted simply by “LL1” rather than “LL1G” since no confusion can arise. These “Gauss elements” have specialized uses.) That labeled “HL4” is the standard Hermite cubic-interpolated element with 4 DOF: two end-node values and two length-scaled (The scaling by $\ell$ is to reduce clutter in the fitting formulas. Scaling is supressed when node rotations are used.) x-derivatives at the end nodes: $[f_1,(df/dx)_1 \ell ,f_2,(df/dx)_2 \ell ]$.

Fig. 30

Line element notation. LL and HL denote Lagrangian and Hermitian interpolation, respectively. Nodes shown as circles have only one DOF: the function value there. Nodes pictured with a rotation (“harpoon”) symbol have two DOF there: the function value and its x-derivative scaled by the line length $\ell$

The shape functions for the line elements shown in Fig. 30 are returned by the Mathematica module LineShapeFunctions listed in Fig. 31. The module is invoked as

$$\begin{aligned} \{ \, {\texttt {Nsf,m},\xi \texttt {nc} \, \} =\texttt{LineShapeFunctions[typ},\xi ,\texttt {Le]}} \end{aligned}$$

(177)

Here typ is one of the identifiers shown in Fig. 30 supplied as a character string (e.g., “LL3”), $\xi$ is the isoparametric natural coordinate, and Le is the element length $\ell$. The latter is only used for type “HL4”. The module returns a list of shape functions in Nfs, the polynomial order m of the shape functions, and a list $\xi$nc of isoparametric coordinates of the nodes. If argument typ is not implemented, Nsf returns Null.

The polynomial fit operation is carried out by the Mathematica module LTFit listed in Fig. 32. The module is invoked as

$$\begin{aligned} \{\,{\texttt {A,B,C,G,F,FG,GF} \, \} =\texttt{LineFit[ftyp,gtyp},\xi \texttt {,Le]}} \end{aligned}$$

(178)

in which ftyp and gtyp are identifiers of the base and fitted line element, respectively, $\xi$ is the isoparametric natural coordinate, and Le is the element length. The latter is only relevant if ftyp or gtyp is “HL4”, but often is set to 1 to reduce clutter in the fitting matrices. The module returns a list of seven matrices: $\mathbf {A}$, $\mathbf {B}$, $\mathbf {C}$, $\mathbf {G}$, $\mathbf {F}$, $\mathbf {F}\mathbf {G}$ (which should be the identity), and $\mathbf {G}\mathbf {F}=\mathbf {H}$ (which should be a projector). These matrices were introduced in Appendix Sect. 6.1.

Fig. 31

Line element shape function module

Fig. 32

Line element fitting module

1.3 6.3 Line Self-Kernel Matrices

Self kernels are the matrices $\mathbf {A}$ and $\mathbf {C}$ introduced in Appendix Sect. 6.1. Since $\mathbf {C}=\mathbf {A}$ if $\mathbf {N}_f\equiv \mathbf {N}_g$, only $\mathbf {A}$ is listed. Its subscript identifies the line element.

$$\begin{aligned} \mathbf {A}_{LL1}&=\left[ \begin{array}{cccc}1\end{array}\right] , \quad \mathbf {A}_{LL2}={1\over 6}\left[ \begin{array}{cccc}2 &{} 1\\ 1 &{} 2\end{array}\right] , \quad \mathbf {A}_{LL3}={1\over 30}\left[ \begin{array}{cccc}4 &{} -\,1 &{} 2\\ -1 &{} 4 &{} 2\\ 2 &{} 2 &{} 16\end{array}\right] , \nonumber \\ \mathbf {A}_{LL4}&={1\over 1680}\left[ \begin{array}{cccc}128 &{} 19 &{} 99 &{} -\,36\\ 19 &{} 128 &{} -\,36 &{} 99\\ 99 &{} -\,36 &{} 648 &{} -\,81\\ -36 &{} 99 &{} -\,81 &{} 648\end{array}\right] , \quad \mathbf {A}_{LL2G}={1\over 2}\left[ \begin{array}{cccc}1 &{} 0\\ 0 &{} 1\end{array}\right] , \nonumber \\ \mathbf {A}_{LL3G}&={1\over 18}\left[ \begin{array}{cccc}5 &{} 0 &{} 0\\ 0 &{} 8 &{} 0\\ 0 &{} 0 &{} 5\end{array}\right] , \quad \mathbf {A}_{HL4}={\ell \over 420}\left[ \begin{array}{cccc}156 &{} 22 \ell &{} 54 &{} -\,13 \ell \\ 22 \ell &{} 4 \ell ^2 &{} 13 \ell &{} -\,3 \ell ^2\\ 54 &{} 13 \ell &{} 156 &{} -\,22 \ell \\ -13 \ell &{} -\,3 \ell ^2 &{} -\,22 \ell &{} 4 \ell ^2\end{array}\right] , \end{aligned}$$

(179)

1.4 6.4 Line Coupling Matrices

Line coupling matrices are the matrices $\mathbf {F}$, $\mathbf {G}$ and $\mathbf {F}$ defined in Appendix Sect. 6.1.

LL2 to LL1:

$$\begin{aligned} \mathbf {B}={1\over 2} \left[ \begin{array}{cccc}1&1 \end{array}\right] , \quad \mathbf {G}={1\over 2} \left[ \begin{array}{cccc}1&1 \end{array}\right] , \quad \mathbf {F}= \left[ \begin{array}{cccc}1 \\ 1 \end{array}\right] . \end{aligned}$$

(180)

LL2 to LL2G:

$$\begin{aligned} \mathbf {B}={1\over 12} \left[ \begin{array}{cccc}a_1 &{} a_2 \\ a_2 &{} a_1 \end{array}\right] , \quad \mathbf {G}=2 \mathbf {B}, \quad \mathbf {F}={1\over 2} \left[ \begin{array}{cccc}a_3 &{} a_4 \\ a_4 &{} a_3 \end{array}\right] , \end{aligned}$$

(181)

in which $a_1=3+\sqrt{3}$, $a_2=3-\sqrt{3}$, $a_3=1+\sqrt{3}$, and $a_4=1+\sqrt{3}$.

LL3 to LL2:

$$\begin{aligned} \mathbf {B}={1\over 6} \left[ \begin{array}{cccc}1 &{} 0 &{} 2 \\ 0 &{} 1 &{} 2 \end{array}\right] , \quad \mathbf {G}={1\over 3} \left[ \begin{array}{cccc}2 &{} -\,1 &{} 2 \\ -1 &{} 2 &{} 2 \end{array}\right] , \quad \mathbf {F}={1\over 2} \left[ \begin{array}{cccc}2 &{} 0 \\ 0 &{} 2 \\ 1 &{} 1 \end{array}\right] . \end{aligned}$$

(182)

LL3 to LL1:

$$\begin{aligned} \mathbf {B}=\mathbf {G}={1\over 6} \left[ \begin{array}{cccc}1&1&4 \end{array}\right] , \quad \mathbf {F}= \left[ \begin{array}{cccc}1 \\ 1 \\ 1 \end{array}\right] \end{aligned}$$

(183)

LL3 to LL3G:

$$\begin{aligned} \mathbf {B}={1\over 36} \left[ \begin{array}{cccc}b_1 &{} b_1 &{} 4 \\ 0 &{} 0 &{} 16 \\ b_1 &{} b_1 &{} 4 \end{array}\right] , \quad \mathbf {G}={1\over 10} \left[ \begin{array}{cccc}b_1 &{} b_1 &{} 4 \\ 0 &{} 0 &{} 10 \\ b_1 &{} b_1 &{} 4 \end{array}\right] , \quad \mathbf {F}={1\over 6} \left[ \begin{array}{cccc}b_3 &{} -\,4 &{} b_4 \\ b_4 &{} -\,4 &{} b_3 \\ 0 &{} 6 &{} 0 \end{array}\right] , \end{aligned}$$

(184)

in which $b_1=3+\sqrt{15}$, $b_2=3-\sqrt{15}$, $b_3=5+\sqrt{15}$, and $b_4=5-\sqrt{15}$.

LL3 to LL2G:

$$\begin{aligned} \mathbf {B}={1\over 12} \left[ \begin{array}{cccc}a_3 &{} a_4 &{} 4 \\ a_4 &{} a_3 &{} 4 \end{array}\right] , \quad \mathbf {G}=2 \mathbf {B}, \quad \mathbf {F}={1\over 2} \left[ \begin{array}{cccc}a_3 &{} a_4 \\ a_4 &{} a_3 \\ 1 &{} 1 \end{array}\right] , \end{aligned}$$

(185)

in which $a_1=3+\sqrt{3}$, $a_2=3-\sqrt{3}$, $a_3=1+\sqrt{3}$, and $a_4=1+\sqrt{3}$.

LL4 to LL3:

$$\begin{aligned} \mathbf {B}={1\over 120} \left[ \begin{array}{cccc}11 &{} 0 &{} 18 &{} -\,9 \\ 0 &{} 11 &{} -\,9 &{} 18 \\ 4 &{} 4 &{} 36 &{} 36 \end{array}\right] , \mathbf {G}={1\over 80} \left[ \begin{array}{cccc}62 &{} 18 &{} 54 &{} -\,54 \\ 18 &{} 62 &{} -\,54 &{} 54 \\ -5 &{} -\,5 &{} 45 &{} 45 \end{array}\right] , \mathbf {F}={1\over 9} \left[ \begin{array}{cccc}9 &{} 0 &{} 0 \\ 0 &{} 9 &{} 0 \\ 2 &{} -\,1 &{} 8 \\ -1 &{} 2 &{} 8 \end{array}\right] . \end{aligned}$$

(186)

LL4 to LL2:

$$\begin{aligned} \mathbf {B}={1\over 120} \left[ \begin{array}{cccc}13 &{} 2 &{} 36 &{} 9 \\ 2 &{} 13 &{} 9 &{} 36 \end{array}\right] , \mathbf {G}={1\over 20} \left[ \begin{array}{cccc}8 &{} -\,3 &{} 21 &{} -\,6 \\ -3 &{} 8 &{} -\,6 &{} 21 \end{array}\right] , \mathbf {F}={1\over 3} \left[ \begin{array}{cccc}3 &{} 0 \\ 0 &{} 3 \\ 2 &{} 1 \\ 1 &{} 2 \end{array}\right] . \end{aligned}$$

(187)

LL4 to LL1:

$$\begin{aligned} \mathbf {B}=\mathbf {G}={1\over 8} \left[ \begin{array}{cccc}1&1&3&3 \end{array}\right] , \quad \mathbf {F}= \left[ \begin{array}{cccc}1 \\ 1 \\ 1 \\ 1 \end{array}\right] . \end{aligned}$$

(188)

LL4 to LL3G:

$$\begin{aligned} \mathbf {B}&={1\over 720} \left[ \begin{array}{cccc} h_1 &{} h_2 &{} h_3 &{} h_4 \\ -20 &{} -20 &{} 180 &{} 180 \\ h_2 &{} h_1 &{} h_4 &{} h_3 \end{array}\right] , \mathbf {G}={1\over 400} \left[ \begin{array}{cccc} h_5 &{} h_6 &{} h_7 &{} h_8 \\ -25 &{} -25 &{} 225 &{} 225 \\ h_6 &{} h_5 &{} h_8 &{} h_7 \end{array}\right] , \nonumber \\ \mathbf {F}&={1\over 54} \left[ \begin{array}{cccc} h_9 &{} -36 &{} h_{10} \\ h_{10} &{} -36 &{} h_9 \\ h_{11} &{} 44 &{} h_{12} \\ h_{12} &{} 44 &{} h_{11} \end{array}\right] . \end{aligned}$$

(189)

in which $h_1=55+11 \sqrt{15}$, $h_2=55-11 \sqrt{15}$, $h_3=45+27 \sqrt{15}$, $h_4=45-27 \sqrt{15}$, $h_5=110+22 \sqrt{15}$, $h_6=110-22 \sqrt{15}$, $h_7=90+54 \sqrt{15}$, $h_8=90-54 \sqrt{15}$, $h_9=45+9 \sqrt{15}$, $h_{10}=45-9 \sqrt{15}$, $h_{11}=5+3 \sqrt{15}$, and $h_{12}=5-3 \sqrt{15}$.

LL4 to LL2G:

$$\begin{aligned} \mathbf {B}={1\over 240} \left[ \begin{array}{cccc} g_1 &{} g_2 &{} g_3 &{} g_4 \\ g_2 &{} g_1 &{} g_4 &{} g_3 \end{array}\right] , \quad \mathbf {G}=2 \mathbf {B}, \quad \mathbf {F}={1\over 2} \left[ \begin{array}{cccc} g_5 &{} g_6 \\ g_6 &{}g_5 \\ g_5 &{} g_6 \\ g_6 &{} g_5 \end{array}\right] \end{aligned}$$

(190)

in which $g_1=15+11 \sqrt{3}$, $g_2= 15-11 \sqrt{3}$, $g_3=45+27 \sqrt{3}$, $g_4=45-27 \sqrt{3}$, $g_5=1+\sqrt{3}$, and $g_6=1-\sqrt{3}$.

Appendix 7: Classical And Modified Functionals Analysis

This Appendix presents Mathematica tools used to test classical functionals and construct modified functionals, as covered in the main text.

1.1 7.1 Testing Classical Functionals

This section presents the symbolic code used to verify the conclusions arrived at in Sect. 5.4, namely that lower order assumptions and one-point integration cannot produce an optimal material stiffness for the Timoshenko beam element.

Fig. 33

Modules that operate at the individual element level

Figure 33 lists four modules that operate at the individual element level. They are invoked as follows.

$$\begin{aligned} {\texttt {Ke=OptTimoMateStiff[Le,EI},\Phi \texttt {]}} \end{aligned}$$

(191)

This module receives as arguments the bending rigidity EI, the element length Le, and the shear flexibility coefficient $\Phi$. It returns the $4\times 4$ optimal stiffness matrix (29) in Ke.

$$\begin{aligned} {\texttt {Ke=StiffMatrixFromEnergy[Ue,ue] } } \end{aligned}$$

(192)

Module StiffMatrixFromEnergy returns the element stiffness matrix Ke as the Hessian of the element generating internal energy functional Ue expressed in terms of the element node displacements collected in array ue. The entries of ue must be symbolic variables to allow differentiation.

$$\begin{aligned} {\texttt {match=StiffMatrixCompare[Ke,Kexact] } } \end{aligned}$$

(193)

Module StiffMatrixCompare compares two stiffness matrices: Ke and Kexact, supplied as arguments. Of these, Ke is normally an allegedly exact stiffness matrix whereas Kexact is the exact one. Both are typically symbolic. The module returns a logical variable match, which is True if an exact entry-by-entry match is found, else False.

$$\begin{aligned} {\texttt {Kc=CondLastFreedom[Ke] } } \end{aligned}$$

(194)

Module CondLastFreedom[Ke] receives an $n \times n$ element stiffness matrix Ke as argument, and statically condenses its last freedom. It returns the $(n-1)\times (n-1)$ condensed stiffness Kc as function value. Two abnormal conditions are checked for:

• If $n\le 1$, the module returns Null.

• If the last diagonal entry—the so called pivot–is exactly zero, normally an error exit should be taken. In the present studies, that can only happen if the last row and column are null, whence the pivot is set to 1 and the condensation effectively eliminates that row and column.

Fig. 34

Module that returns an interpolating function for a beam element as defined by nodal value configuration and polynomial order

Figure 34 lists module InterFun that implements the theory given in Sect. 5.2 by returning an interpolating function given its nodal value configuration and polynomial order. It is invoked as

$$\begin{aligned} {\texttt {fun=InterFun[nodval,order},\xi \texttt {] } } \end{aligned}$$

(195)

Here nodval is the nodal value array, which may contain 2, 3 or 4 values; see Fig. 6. Argument order is the polynomial interpolation order: 0 through 3. Finally $\xi$ is the symbol for the interpolation variable. The module gives the interpolation function as return value. If either the length of nodval or the value of order is out of range, it returns Null.

The number of combinations implemented is $3\times 4=12$. Of these 4 are isoparametric, whereas the others are either least-square fits or isoparametric collocations—see theory in Appendix 6. In the stiffness analysis of Timoshenko models, separate functions for transverse displacements and rotations are retrieved from this module, giving $12\times 12=144$ possibilities.

Examples: InterFun[{v1,v2,v0},0,$\xi$] returns (4v0+v1+v2)/6, which is Simpson’s rule. InterFun[{v1,v2,v0},1,$\xi$] returns the linear least-square fit (4v0+v1(1-$\xi$)+v2(1+$\xi$))/6. InterFun[{v1,v2,v0},2,$\xi$] returns v0(1-$\xi \hat{\,}$^2)/2-v1(1-$\xi$)$\xi$/2+v2(1+$\xi$)$\xi$/2, which is the isoparametric interpolation. Raising the order to 3 returns the same function.

Fig. 35

Modules that generate a Timoshenko beam stiffness matrix (material or geometric) given a functional generating kernel (or kernels), node value configuration, displacement and rotation polynomial interpolation orders, and integration rule

Figure 35 lists two modules that generate a Timoshenko beam stiffness matrix (material or geometric) given a functional generating kernel, node value configutation displacement and rotation polynomial interpolation orders, and integration rule. The calling sequences are similar:

$$\begin{aligned}&{\texttt {Ke=TimoBeamStiff[Qlis,de},\xi \texttt {,Le,ue,pw,p}\theta ,\texttt {intlis,cond] } } \end{aligned}$$

(196)

$$\begin{aligned}&{\texttt {Ue=TimoBeamFunctional[Q,de},\xi ,\texttt {Le,ue,pw,p}\theta \texttt {,int] }} \end{aligned}$$

(197)

The arguments of TimoBeamStiff are

Qlis:

A list that contains one or two functional generating kernels:

For a material stiffness, Qlis is {QMb,QMs}, where kernels QMb and QMs account for the bending and shear energies, respectively.

For a geometric stiffness, Qlis is simply QG, where QG is the kernel for the geometric energy.

de:

The kernel-compatible functional generating vector, such as {v[x], $\theta$[x],v’[x], $\theta '$[x]}.

$\xi$:

The isoparametric coordinate.

Le:

The element length.

ue:

Node displacement vector in symbolic form; e.g., {u1, $\theta$1, v2, $\theta$2}.

pv,p$\theta$:

Polynomial interpolation orders for displacements and rotations, respectively.

intlis:

A list that contains one or two integration rule specifications:

For a material stiffness, intlis is {intb,ints}, where intb and ints specify the integration rules for bending and shear energy contributions, respectively. Each specification is a character string: “A”, “G1”, “G2” or “G3”, for analytical, 1-point Gauss, 2-point Gauss and 3-point Gauss, respectively.

For a geometric stiffness, intlis is simply intg , where intg is the character string that specifies the integration rule for the geometric energy.

cond:

A logical flag. If True and the element stiffness order (the length of ue) is greater that 4, statically condense it to $4\times 4$. If False, return the uncondensed stiffness matrix.

The module returns the element stiffness matrix Ke as function value.

The arguments of TimoBeamFunctional are similar to those of TimoBeamStiff with three differences: Q is the functional generating matrix to be used (not a matrix list), int is the integration rule (not a specification list), and cond is missing, since the returned functional Ue is a scalar.

For a material stiffness, TimoBeamStiff calls TimoBeamFunctional twice to get the bending and shear contributions separately, which allows for testing different interpolation orders and integration rules. It adds the two contributions and takes the Hessian via StiffMatrixFromEnergy, described previously. For a geometric stiffness, it calls TimoBeamFunctional once.

Fig. 36

Modules that return classical functionals for the Timoshenko beam

Figure 36 lists three modules that return classical functionals for the Timoshenko beam specified as generating kernels. They are invoked as

$$\begin{aligned}&{\texttt {QMb=ClassicalQMb[EI,m] }} \end{aligned}$$

(198)

$$\begin{aligned}&{\texttt {QMs=ClassicalQMb[GAs,m] } } \end{aligned}$$

(199)

$$\begin{aligned}&{\texttt {QG=ClassicalQG[F,m] } } \end{aligned}$$

(200)

Modules ClassicalQMb, ClassicalQMs, and ClassicalQG return the functional generating kernels QMb, QMs, and QG for the bending, shear and geometric energy, respectively. These are sparse matrices. The first argument of the modules are the bending rigidity EI, shear rigidity GAs, and axial force F, respectively. The second argument m, specifies the kernel matrix dimension, which is normally 4 for the functional generating vector $\left[ \begin{array}{cccc}v(x),\theta (x),v'(x),\theta '(x)\end{array}\right]$. The shear rigidity may be expressed as a function of EI and $\Phi$ by replacing GAs by 12EI/(Le$^2$).

Fig. 37

Upper cell is driving script that runs 432 input combinations to form the material stiffness matrix of a Timoshenko beam element using the classical functionals. Lower cell reports which combinations do match the nodally exact stiffness (29)

Finally, the upper cell of Fig. 37 lists a driving script that tests the material stiffness classical functionals with 432 input combinations that involve five nested For loops. The combinations that match the nodally exact stiffness (29) are reported in the lower output cell. The negative conclusions of this study as regards low order assumptions are summarized at the end of Sect. 5.4.

1.2 7.2 Finding Modified Functionals

This section gives implementation details on the discovery of some of the modified functionals presented in Sect. 5.5ff. The scripts use several of the modules of the previous section but require new ones, which are presented next.

Fig. 38

Symmetric matrix constructor utilities

Figure 38 lists two utility modules that construct symbolic symmetric matrices and extract their entry information. They are culled from a more comprehensive set presented in [38]. The constructor is invoked as

$$\begin{aligned} {\texttt {S=SymbSymmMatrix[letter,n,zeros] } } \end{aligned}$$

(201)

Argument letter is the first letter (Roman or Greek) of the generated entry names, while $\texttt {n>0}$ is the matrix order.

Argument zeros is a list of entries that are to be set to zero; each entry is of the form i or {i,j}. If the latter, the (i,j) and (j,i) are set to zero. If the former, the i-th row and i-th column are cleared. If a blank list is supplied for this argument, no entries are zeroed. Example: S=SymbSymmMatrix[“Q”,4,zeros], where zeros contains {{1,2} , {4}, {3,3}}, returns {{Q11,0,Q13,0} , {0,Q22,Q23,0} , {Q13,Q23,0,0} , {0,0,0,0}.

The matrix entry extractor is invoked as

$$\begin{aligned} {\texttt {vars=SymmMatrixSymbols[S] } } \end{aligned}$$

(202)

This module returns a list of upper triangle entries of the symmetric matrix S. (Those are typically variables to be supplied to the built-in Solve function.) The list includes only unique symbolic entries (that is, entries with Symbol head), while eliminating any duplicates. Example: suppose S is the above example matrix. Then SymmMatrixSymbols[S] returns Q11,Q13,Q22,Q23 .

Figure 39 lists several utility modules that operate at the individual element level, complementing those in Fig. 33. They return beam element stiffness matrices for the Timoshenko model that become targets for deriving modified functionals. They are invoked as

$$\begin{aligned}&{\texttt {KM=ShearTimoBeamMateStiff[Le,EI},\Phi \texttt {] } } \end{aligned}$$

(203)

$$\begin{aligned}&{\texttt {KG=ConsTimoBeamGeomStiff[Le,F},\Phi \texttt {]}} \end{aligned}$$

(204)

$$\begin{aligned}&{\texttt {KG=QOptTimoBeamGeomStiff[Le,F},\Phi \texttt {]}} \end{aligned}$$

(205)

$$\begin{aligned}&{\texttt {KM=ArbTimoBeamMateStiff[Le,EI},\beta \texttt {]}} \end{aligned}$$

(206)

$$\begin{aligned}&{\texttt {KG=ArbTimoBeamGeomStiff[Le,F}, \{ \, \beta _1,\beta _2,\beta _3,\beta _4 \,\} ,\chi \texttt {]}} \end{aligned}$$

(207)

Fig. 39

Additional individual element utilities, extending those in Fig. 33

Arguments EI, Le, $\Phi$, and F have been defined before. Argument $\beta$ in ArbTimoBeamMateStiff is the material template parameter in (104). The list {$\beta _1$,$\beta _2$,$\beta _3$,$\beta _4$} in ArbTimoBeamGeomStiff fits the parameter list for the geometric natural stiffness template $\mathbf {S}_G$ of (106), whereas $\chi$ is the parameter in (68). The latter affects the matrix $\mathbf {H}$ that produces $\mathbf {K}_G=\mathbf {H}^T \mathbf {S}_G \mathbf {H}$.

The following matrices are provided by the top 3 modules. ShearTimoBeamMateStiff returns the pure shear material stiffness (31) (this module twice calls OptTimoMateStiff, which is listed in Fig. 33), ConsTimoBeamGeomStiff returns the consistent geometric stiffness (36), and QOptTimoBeamGeomStiff returns the quasi-optimal geometric stiffness (126). To specialize to the BE model, set $\Phi$ to zero.

The “arbitrary” modules: ArbTimoBeamMateStiff and ArbTimoBeamGeomStiff, are used to investigate parametrized templates for the material and geometric stiffness, respectively.

Note: it is common to replace either argument EI or GAs with the corresponding expression that introduces the shear flexibility coefficient $\Phi$. For example: GAs=12*EI/Le$\hat{\,}$2, makes the kernel a function of EI, Le and $\Phi$.

Fig. 40

Script to find the modified functional for the optimal Timoshenko material stiffness using linear-linear interpolation and one-point Gauss integration for both bending and shear

The script of Fig. 40 searches for the modified functional that will produce the optimal material stiffness matrix (29) using linear-linear interpolation and one-point Gauss quadrature for both bending and shear. The result is that reported in ().

Some comments on the code: QMb and QMs are the parametrized kernels for the bending and shear functional, respectively. They are formed as combining the classical functional kernels QMbcla and QMscla with the deviatoric kernels QMbdev and QMsdev. Each of the latter has 10 symbolic variables for a total of 20, which are collected in vars. The beam stiffness matrix Kem is a function of those variables, The built-in Solve is used to find those by zeroing the difference between Kem and the target exact stiffness. The result is replaced into the parametrized kernels to obtain the modified functional kernel QMmod and into the parametrized stiffness matrices for verification.

Fig. 41

Script to find the modified functional for the parametrized Timoshenko geometric stiffness using linear-linear interpolation and one-point Gauss integration

The script of Fig. 41 searches for the modified functional that will produce the template-parametried geometric stiffness matrix returned by ArbTimoGeomStiff using linear-linear interpolation and one-point Gauss quadrature. The result is that reported in Sect. 5.7.

The code is somehat simpler than that of the previous script because no energy splitting into bending and shear is necessary, whence only 10 symbolic variables have to be determined by the match.

Appendix 8: Historical Summary

This Appendix links the theme of the paper and the methods listed in Fig. 1 to their origins.

1.1 8.1 Beam Models

The Bernoulli–Euler (BE) beam model synthesizes pioneer work by Jacob and Daniel Bernoulli as well as that of Leonhard Euler in the XVIII Century [88, p. 25]. Although the model was first enunciated by 1750, it was not seriously applied in structural design and analysis before 1820. (Structural theory was largely developed during the period 1820–1870, as prompted by the rapid growth of metal construction.) While Galileo Galilei is credited with first attempts at a theory, recent studies collected in [6] argue that Leonardo da Vinci made crucial observations a century before Galileo. However, da Vinci lacked Hooke’s law and calculus to complete the theory.

The Timoshenko beam model includes a first order correction for transverse shear flexibility. It was introduced in a short 1921 note [86]. The material was expanded in the 1955 textbook [89].

The governing differential equations for a linearly elastic, shear-flexible beam had actually been presented decades earlier by Bresse in his 1859 book [14], and are in fact referred to as the “Navier-Bresse (N-B) equations” by French authors. Those suffer from an omission: the effective shear area that takes into account the non uniformity of the transverse shear stress distribution. In his history book [88] Timoshenko quotes several passages of [14] in the original French, showing familiarity with that work, as well as emphasizing Navier’s contributions in his Introduction. No references to the N-B equations, however, appear in [86, 89].

The Timoshenko model is thoroughly examined in [47], wherein it is compared to three simpler ones for vibration and dynamic response analyses. Those are important application areas, insofar as the simpler BE model suffers from defects such as an infinite wavespeed [44].

1.2 8.2 HODES Methods

Ordinary differential equations (ODE) were applied to problems in mechanics shortly after the invention of calculus by Newton and Leibniz by the end of the XVII Century. Celebrated early contributions are those by Euler on the stability of beam-columns and elastica theory [23, 24]. The “HODES method” entry in the table of Fig. 1, however, focuses on a narrower and more recent topic: use of ODE solutions in close ancestors of the current FEM. More specifically: structural analysis methods that select kinematic unknowns such as displacements and deformations.

Following the historical outline by Kurrer [51], priority is given to Ostenfeld as the developer of the Deformation Method (Die Deformationsmethode) in the 1920s [75]. This was put in matrix form by Duncan and Collar in the 1930s for aeronautical applications [21, 22, 42]; a detailed account is provided in [29]. It became the Matrix Displacement Method through linkage to XIX Century energy theorems by Levy [53] and Argyris [3]. (The latter formulated it as dual of the Force Method.) The evolution culminated in the Direct Stiffness Method (DSM) of Turner [95, 96], which is the FEM variant in universal use today.

A theorem by Tong [92] confirmed the HODES protocol for 1D problems: to get nodal exactness the homogeneous solutions of the Euler-Lagrange equations must be used as shape functions.

HODES is also the backbone of the Transfer Matrix Method (TMM), which was popular on the limited-memory computers of the late 1950s [80]. Although restricted to linear, one-dimensional, chain-like models, TMM is still in use for special applications, not necessarily structural.

1.3 8.3 Classical Functionals

This broad class of methods is by far the most popular today for structural FEM: discrete function bases (typically polynomials) are inserted in a governing functional, integrated over the element, and variational derivatives taken to produce the discrete stiffness equations. Surprisingly, linkage to the pre-FEM, older Rayleigh-Ritz (RR) methods was not explicitly proclaimed until Melosh’s 1962 thesis under Martin [59] and the subsequent AIAA paper [60]. This interpretation had important consequences:

• Regularized the Matrix Displacement Method of the 1950s, which was based on application of energy theorems such as those named after Maxwell and Castigliano. Such theorems, intended for truss and framework structures common in the XIX Century, become questionable when applied to 2D and 3D continuum models; see [48, Ch. 5].

• Provided a platform for integrating FEM into RR approximation theory. This attracted the interest of applied mathematicans and launched a wave of theoretical research. The first monograph on the mathematics of FEM appeared in 1973 [85].

• Facilitated the extension of FEM to non-structural problems. The first paper of this nature appeared in 1965 [107] and inaugurated a soaring applications period. For field problems lacking classical functionals, such as Navier-Stokes fluids, the Weighted Residual Method of Crandall [19], which unified older variants such as Galerkin, least-squares, subdomain and collocation, provided the necessary tools once appropriately tailored to fit FEM discretizations.

The main drawback in complying with the classical RR rules: conformity and exact integration, was slow convergence of low-order elements and the emergence of locking phenomena. These were addressed by a plethora of “variational crimes” that bypassed the RR rituals: nonconforming elements, reduced and selective integration, and the ever-present patch test.

To achieve higher accuracy while keeping the number of elements reasonably small, the popular approach has been to move onto higher order elements. Those are put to use in hierarchical schemes known in the literature as h, p and m refinement. That methodology is highly developed for linear static analysis. In nonlinear and dynamic analyses, however, simple elements remain popular for obvious reasons: computational efficiency, programming simplicity, and easier control of numerical stability. By far the simplest nonlinear elements are those evaluated with one point Gauss integration, complemented with ad-hoc stabilization devices as needed [9].

For elements involving rotational degrees of freedom (DOF), such as beams, plates and shells, use of one-point integration brings difficulties that have attracted much attention since the early 1970s [49]. This paper examines those frustrations for the simplest problem of that nature: the Timoshenko plane beam, which includes the BE model as a limit case.

Key question: can the material stiffness of that element be made rank sufficient and nodally exact in the linear static case under linear-linear assumptions and one point integration? It is shown in Sect. 5.5 that the answer is yes. However, nonstandard mathematical tools (some old, some newer) are required. Those tools are extended to the geometric stiffness matrix in Sect. 5.6 for the first time. These acomplishments bring hope for extension to plates and shells, as well as to geometrically nonlinear problems.

1.4 8.4 Modified Differential Equations

Modified differential equations, acronymed as MoDE, is a subset of a larger and highly active field: backward error analysis or BEA. Paraphrasing Griffiths and Sanz-Serna [45]: given a (mathematical) problem $\mathcal {P}$ with true solution ${\mathcal {S}}$, and given an approximate solution $\tilde{\mathcal {S}}$, the classical forward error analysis (FEA) compares the distance between $\tilde{\mathcal {S}}$ and ${\mathcal {S}}$. On the other hand, backward error analysis (BEA) shows that the approximation $\tilde{\mathcal {S}}$ is the true solution of a neighboring problem $\tilde{\mathcal {P}}$, and is concerned with the distance from $\tilde{\mathcal {P}}$ to ${\mathcal {P}}$.

Curiously, a BEA-like interpretation was in the minds of early pioneers of FEM. For example, in [17, p. 87] Clough observes that the FEM solution—ignoring computational errors—is exact for the FEM-idealized structure. Thus one may think about the distance between the two structures as an accuracy measure. In Clough’s words: “a modified structural system is substituted for the actual continuum.” Such interpretation clearly distinguishes FEM from the Finite Difference Method. In the latter, the exact equations of the actual physical system are solved approximately.

BEA has by now become a standard tool of the linear algebra community. It has impeccable roots: early use (in implicit form) by von Neumann and Goldstine in 1947 [98], and Turing in 1948 [93]. Following its popularization by Wilkinson [103], BEA has become standard material in linear algebra textbooks dealing with matrix computations.

MoDE is BEA applied to ODEs. This marriage has been less popular, because forward error has been the traditional accuracy measure for the numerical treatment of ODEs. The combination was called “modified equation” by Warming and Hyett in 1974 [102]. This paper analyzes dissipation properties of finite difference methods in computational fluid dynamics (CFD) and was the first to properly do higher derivative elimination. (References to previous work by Hirt, Lomax, Morton, Ritchmyer, Roache, and Tyler may be found there.) Its use in dynamics has been most productive in topics such as long time integration of planetary orbits, and, more generally, Hamiltonian systems treated with symplectic integrators. See, e.g., the textbook [46] and references therein.

The use of MoDE to analyze FEM discretizations can be traced back to [99], which actually preceeds [102] by 6 years. Unfortunately the higher order derivative elimination procedure used there was incorrect, leading to false conclusions on error order and convergence. More successful were the investigations of Park and Flaggs in the 1980s [40, 76,77,78]. These focused on the analysis of shear-flexible beam, plates and shell elements with the chief objective of detecting and tracing deficiencies such as locking. The study featured the use of two computer algebra systems: Macsyma and SMP; the latter being the ancestor of Mathematica.

The use of MoDE in this paper differs from the above ones in that the focus is on obtaining maximum accuracy of template configurations. In particular, it exhibits the first use of a powerful combination: templates-plus-MoDE, for improving the beam geometric stiffness matrix.

A parting question: can MoDE be effectively used for multiple space dimensions? Specifically, moving from ODEs to PDEs. There is no fundamental obstacle: a two-element patch such as that shown in Fig. 16 can be obviously generalized to a 2D or 3D one. But challenges emerge:

• Assessing how the modified PDE will depend on the patch configuration, and

• Significant higher expense in symbolic computations.

Whereas the second challenge may be eventually overcome with increasing computer power and more efficient computer algebra systems, the first one remains unanswered.

1.5 8.5 Templates

The template approach to finite element construction did not come suddenly out of the blue, but gradually emerged from a series of discoveries. The starting point was work of the first author (CAF) with Bergan in 1982 while he was on a sabbatical at Stanford. The work focused on the development of high performance 3-node triangle elements for plane stress [13] and Kirchhoff plate bending [25] using the Free Formulation (FF) of Bergan and Nygård [11, 12]. The elements were eventually combined to form a shell element and implemented into the corotational nonlinear program FENRIS used by Der Norske Veritas for analyzing marine structures under extreme conditions. The plane stress element was the first 3-node membrane triangle with drilling DOF that was rank sufficient and free of aspect ratio locking.

As the name implies, FF allows substantial freedom in the development of elements. In particular conformity can be ignored; all that matters is satifying the Individual Patch Test that can be tested on a single element. As previously discussed in Appendix Sect. 2.2, the material stiffness naturally splits into a basic component, which takes care of convergence, and a higher order component, which takes care of stability (full rank) and accuracy. Ensuing developments showed that free scalar parameters appeared naturally in both components. Which values should be given? Obviously that decision cannot be left to the discretion of program users, who might have no idea of the background theory. They were assigned using higher order patch tests. Free parameters kept appearing in FF succesors, such as the Assumed Natural Deviatoric Strain (ANDES) method [26, 61], which combined the matrix split of FF with the Assumed Natural Strain (ANS) method of Park and Stanley [79].

From these special cases a general approach eventually took shape. It was named the template method in a 1994 journal paper [27]. A fuller historical account is provided in a tutorial chapter [33]. The general concept of template as parametrized forms of FEM matrix equations is discussed in [28, 32, 35].

The main advantage of templates is flexibility. This brings the attendent ability to customize elements without worries about complying with restrictions such as interelement conformity. The approach can generate all possible convergent elements; the main question left is which one to pick up. This rosy picture, however, is clouded by practical difficulties:

• Element development computations must be done symbolically, carrying along material, fabrication and geometric information as variables. Such manipulations lie beyond human powers, and require a computer algebra system (CAS).

• Symbolic computational work effort grows exponentially in the number of free parameters as well as matrix order.

With current CAS the treatment of 1D elements is possible, as well as that of 2D elements of simple geometry. But 2D and 3D elements of arbitrary geometry (for instance, doubly curved shells) lie outside the scope of current computers and CAS power.

These daunting limitations explain why it is so important to reduce the number of free parameters by taking advantage of orthogonality properties. This technique is clearly illustrated in Appendices 2 and  3. The use of Legendre polynomials as a natural basis allowed the use of only 1 parameter for the material stiffness and 4 for the geometric stiffness studies carried out there.

1.6 8.6 Finite Increment Calculus

The Finite Increment Calculus (FIC) is a residual-based formulation of computational methods in science and engineering, proposed in 1998 [62]. The original name: Finite Calculus, was recently changed to avoid clashing with a particular mathematical subdiscipline. For recent application surveys see [64, 74].

The original motivation for FIC was stabilization of flow problems involving advection and diffusion. A unique feature is that stabilizion is done at the governing equation level. This is in contrast to other residual based stabilization methods, in which such terms are introduced a posteriori in the governing functional or weak form. The use of steplengths to improve accuracy has also received attention [36, 67, 73].

FIC has been particularly successful for stabilization of advective (transport) incompressible flow processes [66,67,68,69,70,71, 73]. The method has been also tried for quasi-incompressible solid mechanics set up in Galerkin form [63, 65] since models used in such problems were not encompassed by classical TPE functionals.

Only one particular feature of FIC is used in Sect. 5.7: introducing a steplength-scaled second derivative in a modified functional. A limitation is that the use of a linear-linear kinematic assumption is precluded because second derivatives vanish. Therefore a quadratic-linear assumption with one point integration was tried, which worked successfully.

1.7 8.7 Modified Functionals

The modified functional method debuts in this paper, so a historical track is lacking. But some closely related, largely empirical, schemes used since 1960 deserve mention. It was observed in Sect. 5.5 that passing from the classical TPE functional (17) to the modified one (53) involves a simple modification of the constitutive equations. This idea has in fact been used several times in FEM applications. Three instances are mentioned.

The first FEM analysis of a civil structure concerned the analysis of crack propagation of the Norfork gravity dam, carried out in 1961 by Wilson, then a doctoral student under Clough [16]. A “crack element” was constructed by a simple device: inserting a pseudo orthotropic material into a plane strain 3-node triangular element. The elastic modulus normal to the crack direction, which was found through an iterative scheme, was set to zero. No other change was made and the reported results were reasonable. Incidentally, the program used for this analysis was the first FEM code listed in a publicly available report [104]. A closely related idea is the empirical construction of “compressible fluid elements” by setting the shear modulus to zero.

A more recent event was the formulation of accurate shear-flexible beam, plate and shell elements by MacNeal [55], by adjusting the constitutive coefficients to eliminate shear locking. The technique has been labeled residual bending flexibility or RBF; see [56] or [72, Ch. 2]. These elements have been included in the NASTRAN program since 1980.

The last example is the construction of the nodally exact Timoshenko beam element through a bending-shear energy splitting, as narrated in Sect. 5.9. When constructing the equivalent web as a isoparametric 4-node plane stress element, the input $3\times 3$ stress-strain matrix has only one nonzero element: the scaled shear rigidity $C_w GA_s$. This trick works perfectly as long as the web subelement is rectangular.

The advantage of the method introduced here is that there is no need to investigate locking: the correct modified functional emerges directly by simply matching optimal forms. Furthermore the dependence of the functional on kinematic assumptions (linear-linear, quadratic-linear, etc.) is explicitly found without the need of recourse to empirical arguments. As noted in the Conclusions, however, an air of paradox emerges since finding the functional requires a priori knowledge of the optimal matrix form. Conclusion: modified functionals cannot be derived in vacuo.