# Comparison of the Von Kármán and Kirchhoff models for the post-buckling and vibrations of elastic beams

Sébastien Neukirch, Morteza Yavari, Noël Challamel, Olivier Thomas

► To cite this version:

Sébastien Neukirch, Morteza Yavari, Noël Challamel, Olivier Thomas. Comparison of the Von Kármán and Kirchhoff models for the post-buckling and vibrations of elastic beams. Journal of Theoretical, Computational and Applied Mechanics, 2021, Issue 2 | 2022, pp.1-18. 10.46298/jtcam.6828. hal-02957425v3

HAL Id: hal-02957425

https://hal.science/hal-02957425v3

Submitted on 19 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.6828

OAI hal-02957425v3

**History**

Received Oct 7, 2020

Accepted Mar 11, 2021

Published May 17, 2021

**Associate Editor**

Alexander Popp

**Reviewers**

Yury Vetyukov
Oliver Weeger
Christoph Meier

**Open Review**

DOI 10.5281/zenodo.4766972

**Supplementary Material**

Data
OAI hal-02957425v2

**Licence**

CC BY 4.0
©The Authors

# Comparison of the von Kármán and Kirchhoff models for the post-buckling and vibrations of elastic beams

Sébastien Neukirch<fnref n="1" />, Morteza Yavari<fnref n="2" />, Noël Challamel<fnref n="3" />, and Olivier Thomas<fnref n="4" />

We compare different models describing the buckling, post-buckling and vibrations of elastic beams in the plane. Focus is put on the first buckled equilibrium solution and the first two vibration modes around it. In the incipient post-buckling regime, the classic Woinowsky-Krieger non-linear model is known to grasp the behavior of the system. It is based on the von Kármán approximation, a 2nd order expansion in the rotation and vertical displacement of the buckled beam. But as the rotation and the vertical displacement in the beam become larger, the Woinowsky-Krieger model starts to show limitations and we introduce a 3rd order model, derived from the geometrically exact Kirchhoff model. We discuss and quantify the shortcomings of the Woinowsky-Krieger model and the contributions of the 3rd order terms in the new model, and we compare them both to the Kirchhoff model. Furthermore, we show that the limit in the validity range of the Woinowsky-Krieger model is only marginally affected by the slenderness ratio of the beam. Different ways to nondimensionalize the models are compared and we believe that, although this study is performed for specific boundary conditions, the present results have a general scope and can be used as abacuses to estimate the validity range of the simplified models.

**Keywords:** nonlinearities, postbuckling, natural frequencies

## 1 Introduction

Every model is wrong (George 1976), but a good model is both accurate and easy to handle. In mechanical engineering, a trade-off is usually made between accuracy and computability. When looking at the deformation of elastic structures, simplification in the kinematics or constitutive relations are for example performed to ease calculations. Here, we investigate the post-buckling and vibration behavior of elastic beams using both geometrically exact and approximate models. In particular, we question the validity of semi-linearized models and their efficiency to capture the nonlinear response of elastic beams. The equations of motion for extensible, geometrically-exact beams have been established by Kirchhoff (1876) and generalized by Reissner (1972) to include shear effects. There are several recent textbooks devoted to nonlinear structural models, see e.g. (Nayfeh and Pai 2004; Lacarbonara 2013; Luongo and Zulli 2013; Audoly and Pomeau 2010; Dill 1992) for a nice historical analysis, and (Meier et al. 2019; Thomas et al. 2016; Cottanceau et al. 2017) for recent reviews on finite elements approaches and continuation methods for elastic beams. Here, we deal with a nonlinear beam problem and account for both bending and extensional deformations (while neglecting shear and rotational inertia). Due to the difficulties to find exact solutions in the nonlinear case, approximate engineering models have been formulated. These models rely on simplified kinematics, either linearized or weakly nonlinear, and include a coupling between axial and bending motions. The so-called Woinowsky-Krieger model (Woinowsky-Krieger 1950; Mettler 1951; Eringen 1952; Burgreen 1951) assumes a linearized curvature calculation and a von Kármán-type axial strain measure, first introduced for the buckling of elastic plates (von Kármán 1907; Eisley 1964). This model has been widely used in the literature and has shown its efficiency for computing approximate

<footnote n="1">¹ Sorbonne Université, CNRS, Institut Jean Le Rond d'Alembert UMR 7190, Paris, France</footnote>

<footnote n="2">² Department of Physics, Islamic Azad University, Kashan Branch, Kashan, Iran</footnote>

<footnote n="3">³ Université de Bretagne-Sud, CNRS, Institut de Recherche Dupuy de Lôme UMR 6027, F-56321 Lorient, France</footnote>

<footnote n="4">⁴ Arts et Métiers Institute of Technology, LISPEN, HESAM Université, F-59046 Lille, France</footnote>
amplitude-frequency dependence of extensible elastic beams in the weakly nonlinear regime, see for example (Ray and Bert 1969; Lou and Sikarskie 1975; McDonald 1991; Nayfeh and Mook 1995; Nayfeh and Emam 2008; Thomas et al. 2016) and references therein. Furthermore, exact solutions have been derived for the Woinowsky-Krieger model (Nayfeh et al. 1995). However, as mentioned in (Gao 2000), it is based on a linearization of the transverse displacement equation, so that computations of the nonlinear behavior of beams are only valid under certain conditions and for small deflections. For the geometrically exact case, Kirchhoff’s equations for the extensible beam have been reformulated by Pflüger (1964) who gave the exact buckling load of the extensible column. Moreover, exact solutions for the equilibrium of extensible columns in term of elliptic integrals have been derived, see for instance (Goto et al. 1990).

In this paper, we investigate the range of validity of the Woinowsky-Krieger beam model by comparing it to the geometrically exact Kirchhoff model. We consider clamped-clamped boundary conditions and analytically and numerically compute the planar equilibrium and vibrations of a beam in a displacement-controlled loading. We use the Kirchhoff extensible model in Section 2, and the Woinowsky-Krieger model in Section 3, to compute equilibrium and vibrations in the post-buckling regime, and we then compare results from the two models in Section 4. As the Woinowsky-Krieger model is only 2nd order, we derive in Section 5 a new model, comprising third order terms, and we compare it to the two previous models. We discuss our findings and conclude in Section 6.

## 2 The Kirchhoff model

We look at the planar equilibrium and vibrations of a beam in the post-buckling domain, as illustrated in Figure 1. The beam has length L, with a homogeneous cross-section of area A and second moment of area I, cast from a homogeneous and isotropic material of Young’s modulus E and density ρ. An important parameter is the slenderness ratio

$$ \eta = I/(AL^2) \qquad (1) $$

that becomes $η = (h/L)^2/12$ for a rectangular cross-section of width w and thickness h. We adopt the Euler-Bernoulli assumptions, that is we neglect shear deformations and rotational inertia, which only become important in the high frequency domain (Timoshenko 1921). However, we keep extensional deformations because in displacement-controlled setups, extension plays a crucial role to describe the vibrations of the beam (Neukirch et al. 2012; Neukirch et al. 2014). The beam is naturally flat, the reference configuration being along the horizontal axis. We use the arc-length S of the beam in its reference configuration as a Lagrangian variable, that is $S \in [0, L]$ always. The beam is clamped horizontally at its left end $S = 0$, which lies at the origin. The right end at $S = L$ is constrained to lie on the horizontal axis, with a horizontal tangent, see Figure 1. An axial displacement D is imposed and we compute the equilibrium shape and the vibrations

<figure><img src="image_2.png" /><figcaption>**Figure 1** Clamped-clamped beam with imposed axial displacement D. Note that $V = Y$ and $U = X - S$. The external force at the right end has horizontal $P(T) = -N_x(L, T)$ and vertical $Q(T) = -N_y(L, T)$ components, with physical time T.</figcaption></figure>

around this shape. We then study how equilibrium and vibrations vary as D is changed.

In this section, we present the extensible Kirchhoff model, which we tend to regard as the reference with respect to which we compare the models of the subsequent sections. Kirchhoff’s
framework uses the current position (X, Y) and inclination angle θ as kinematic variables, and the internal bending moment M and force vector (Nₓ, Nᵧ) as stress variables. Linear bending M = EI dθ/dS and stretching Nθ = EAe constitutive relations are used, where e(S, T) is the extension of the beam and Nθ = Nₓ cos θ + Nᵧ sin θ is the tension in the beam. The motion of the beam is given by

$$X' = (1 + e) \cos \theta \qquad N'_{x} = \rho A \ddot{X} \qquad (2a,b)$$

$$Y' = (1 + e) \sin \theta \qquad N'_{y} = \rho A \ddot{Y} \qquad (2c,d)$$

$$EI\theta' = M \qquad M' = N_{x}Y' - N_{y}X' \qquad (2e,f)$$

$$EAe = N_{x} \cos \theta + N_{y} \sin \theta. \qquad (2g)$$

Every variable depends on both the arc-length S and the time T with the notations ()' = d()/dS and (·) = d()/dT. Unless otherwise stated, from now on, we work with non-dimensionalized variables, that is we use L as unit length, EI/L² as unit force, and τ = L²ρA/(EI) as unit time. Non-dimensionalized variables are written in lower case, e.g. x = X/L, or nₓ = NₓL²/(EI). Also note that in the remainder Ω = ω/τ is the physical angular frequency in radians per second. The non-dimensionalized version of System (2) is simply obtained by setting EI = 1, L = 1, ρA = 1, and EA = 1/η:

$$x' = (1 + e) \cos \theta \qquad n'_{x} = \ddot{x} \qquad (3a,b)$$

$$y' = (1 + e) \sin \theta \qquad n'_{y} = \ddot{y} \qquad (3c,d)$$

$$\theta' = m \qquad m' = n_{x}y' - n_{y}x' \qquad (3e,f)$$

$$e = \eta(n_{x} \cos \theta + n_{y} \sin \theta) \qquad (3g)$$

with ()' = L d()/dS and (·) = τ d()/dT, see (Neukirch et al. 2012; Neukirch et al. 2014) for more details. We stress that this beam model (3), and thus its solution, solely depends on η, defined in Equation (1).

The equilibrium solution (x_E, y_E, θ_E, m_E, n_{xE}, n_{yE}, e_E) is found by solving Equation (3) with x̄ = 0 and ŷ = 0. Once the equilibrium is known, we compute vibrations by using the ansatz functions

$$x(s, t) = x_E(s) + \delta\bar{x}(s) \cos \omega t, \\ y(s, t) = y_E(s) + \delta\bar{y}(s) \cos \omega t, \qquad (4)$$

...

with δ ≪ 1 and ω, the non-dimensionalized angular frequency. Injecting Equation (4) into System (3) and keeping only first order terms in δ yields the following linear differential system for the vibration modes x̄, ŷ, θ̄, m̄, nₓ, nᵧ:

$$\begin{align}
\bar{n}'_y &= -\omega^2 \bar{y} & \bar{m}' &= \bar{n}_x y'_E - \bar{n}_y x'_E + n_{xE} \bar{y}' - n_{yE} \bar{x}' \nonumber \\
\bar{n}'_x &= -\omega^2 \bar{x} & \bar{y}' &= (1 + e_E) \cos \theta_E \bar{\theta} + \bar{e} \sin \theta_E \nonumber \\
\bar{\theta}' &= \bar{m} & \bar{x}' &= -(1 + e_E) \sin \theta_E \bar{\theta} + \bar{e} \cos \theta_E \tag{5}
\end{align}$$

with ē = η[ n̄_y sin θ_E + n̄_x cos θ_E + (n_{yE} cos θ_E - n_{xE} sin θ_E)θ̄ ]. The clamped-clamped boundary conditions read

$$y_E(0) = 0 = y_E(1), \qquad \theta_E(0) = 0 = \theta_E(1), \qquad x_E(0) = 0 = x_E(1) - 1 + d \qquad (6a)$$

$$\bar{y}(0) = 0 = \bar{y}(1), \qquad \bar{\theta}(0) = 0 = \bar{\theta}(1), \qquad \bar{x}(0) = 0 = \bar{x}(1) \qquad (6b)$$

where d = D/L is the non-dimensionalized axial displacement. We note that in this displacement-controlled setup, the position x_E(1) is fixed, but the applied axial p(t) and shear q(t) forces vary with time, and we have nₓ(1, t) = -p(t) and nᵧ(1, t) = -q(t), see Figure 1.

We are eventually left with a nonlinear boundary value problem (3), (5) and (6) that we solve numerically for a beam with η = 1/4800 (that is L = 20h in the case of a rectangular cross-section).
We focus on the first buckling mode, which has $n_{yE} = 0$ (Domokos 1994; Neukirch et al. 2014), and on the first two vibration modes around it. We used the ordinary differential equations integrator of the Mathematica software, coupled with a shooting technique, to solve the boundary value problem, and a pseudo arc-length continuation procedure to compute the evolution of the solution as the parameter $D/L$ is varied. Curves in Figures 2, 3, 4 and 5 were generated in few minutes.

<figure><img src="image_3.png" /><figcaption>**Figure 2** Post-buckled equilibrium of the clamped-clamped beam for $\eta = 1/4800$: Kirchhoff model (Ki.), Woinowsky-Krieger model (WK), order 3 model of Section 5 (Order 3), and fourth order development of $p_E$ and $d$ (Ki.dev 4). Equilibrium shape for $d = 0.2$. Both the order 3 and Woinowsky-Krieger models tend to the classical Euler value $P_E L^2 / (EI) = 4\pi^2 \approx 39.5$ at $D/L = 4\pi^2\eta$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3** Post-buckled equilibrium of the clamped-clamped beam: transverse displacement at midspan $y_E(1/2)$ vs axial displacement $d$ for $\eta = 1/4800$. Kirchhoff model (Ki.), Woinowsky-Krieger model (WK), order 3 model of Section 5 (Order 3), and 4th-order development of $y_E(1/2)$ and $d$ (Ki.dev 4). Inset: Equilibrium shape for $d = 0.8$.</figcaption></figure>

In Figure 2, we plot the external axial force $p_E$ as a function of the axial displacement $d$. In Figure 3, we plot the transverse displacement at midspan $y_E(1/2)$ as a function of the axial displacement $d$. In Figures 4 and 5, we plot the angular frequency $\omega$ of the first two vibration modes, as a function of the axial displacement $d$. These plots will be analyzed in Section 4. Please see the supplementary material (Neukirch et al. 2021) for plots with different values of $\eta$.

## 3 The Woinowsky-Krieger model

We now turn to a simplified model to describe the same equilibrium and vibrations experiment. This model, which was introduced in (Woinowsky-Krieger 1950) to correct the fully linear approach, includes the axial/bending coupling that arises when the transverse displacement of the beam becomes finite. This model was also introduced earlier in the German book by Kirchhoff (1876, Eq. (16), p. 441).
<figure><img src="image_5.png" /><figcaption>**Figure 4** Vibration curve for the first vibration mode around the post-buckled equilibrium solution and $\eta = 1/4800$. Kirchhoff model (Ki.), Woinowsky-Krieger model (WK), order 3 model of Section 5 (Order 3). Inset: equilibrium shape and extremal modal shapes of the first vibration mode for $d = 0.2$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5** Vibration curve for the second vibration mode and $\eta = 1/4800$. Kirchhoff model (Ki.), Woinowsky-Krieger model (WK), and order 3 model of Section 5 (Order 3). Inset: equilibrium shape and extremal modal shapes of the second vibration mode and $d = 0.2$.</figcaption></figure>

It is based on the same assumption as the one used by von Kármán for the statics of plates (von Kármán 1907), which consists in keeping only the first nonlinear term in the expansion of the axial strain $e$ as a function of the cross-section rotation $\theta$. Namely, the term $\cos \theta$ in Equation (3a) is treated up to the second order, leading to

$$e = U' + \frac{1}{2}\theta^2, \qquad (7)$$

with $U(S, T) = X(S, T) - S$, the axial displacement of the cross-section. This assumption is energetically consistent with the approximation which replaces Equation (3f) with $N_y = -M' + N_x Y'$ (see Section C.2 for a variational approach to this model). As a second assumption, the axial inertia is neglected in this model. From Equation (3b), this omission leads to a uniform axial force, $N_x(S, T) = N_x(T) = -P(T)$. Consequently, Equation (3c) is treated linearly in $\theta$ and $e$ is neglected with respect to 1, yielding $Y' = \theta$. Finally, Equation (3g) is truncated to the zero-th order in $\theta$, that is $EAe = N_x$. Combining all these equations and keeping only $U(S, T)$ and $Y(S, T)$ as unknowns leads to

$$EIY''''(S, T) + \rho A \ddot{Y}(S, T) + P(T)Y''(S, T) = 0 \qquad (8a)$$

$$EA[U'(S, T) + \frac{1}{2}Y'^2(S, T)] = -P(T). \qquad (8b)$$

From here, we have two ways to write these equations in a dimensionless form. The first way consists in using the same dimensionless variables as for the Kirchhoff model of Section 2, and yields

$$\ddot{y} + y''' + py'' = 0 \qquad (9a)$$

$$u' + \eta p + \frac{1}{2}y'^2 = 0 \qquad (9b)$$
with $u = U/L$. In this case, the behaviour of the beam depends solely on the slenderness ratio $\eta$.

The second way to introduce dimensionless variables is to scale the transverse displacement $Y$ with the radius of gyration $r = \sqrt{I/A} = L\sqrt{\eta}$ and the axial displacement $U$ with $r^2/L$, one order of magnitude smaller. The physical meaning of $r$ is a characteristic thickness of the cross-section. In particular, for a rectangular cross section, $r = h/\sqrt{12}$. Writing $\hat{y} = Y/r = y/\sqrt{\eta}$ and $\hat{u} = UL/r^2 = u/\eta$ recasts Equation (9) as

$$ \ddot{y} + \hat{y}'''' + p\hat{y}'' = 0 \qquad (10a) $$

$$ \hat{u}' + p + \frac{1}{2}\hat{y}'^2 = 0 \qquad (10b) $$

which does not depend on any geometrical or material parameter. This shows that any beam modelled by the Woinowsky-Krieger model exhibits the same mechanical behaviour. However, this scaling cannot be applied to the Kirchhoff model for which the dependence on the slenderness ratio $\eta$ cannot be avoided. Consequently, in the following, when comparing models we will use the set of dimensionless variables of the Kirchhoff model and thus Equation (9).

The equilibrium version of Equation (9) is

$$ y_E''' + p_E y_E'' = 0 \qquad (11a) $$

$$ x_E' - 1 + \eta p_E + \frac{1}{2} y_E'^2 = 0 \qquad (11b) $$

which is completed by

$$ \bar{y}'''' + p_E \bar{y}'' + \bar{p} y_E'' = \omega^2 \bar{y} \qquad (12a) $$

$$ \bar{x}' + \eta \bar{p} + y_E' \bar{y}' = 0 \qquad (12b) $$

for the vibrations. The great advantage of this model is that, although nonlinear, it can be solved analytically (Nayfeh et al. 1995; Mamou-Mani et al. 2009). The equilibrium solution is

$$ y_E(s) = \frac{1}{2}\epsilon(1 - \cos 2\pi s) \qquad \Rightarrow \qquad y_E(1/2) = \epsilon \qquad (13a) $$

$$ x_E(s) = s(1 - \eta 4\pi^2) + \frac{1}{16}\epsilon^2\pi(\sin 4\pi s - 4\pi s) \qquad \Rightarrow \qquad d = \eta 4\pi^2 + \frac{1}{4}\pi^2\epsilon^2 \qquad (13b) $$

$$ p_E = 4\pi^2 \qquad (13c) $$

where the amplitude $\epsilon$ of the linear solution $y_E(s)$ is defined in Equation (16). The solution for the vibrations is

$$ \bar{y}(s) = c_1 \sin ns + c_2 \cos ns + c_3 \sinh ms + c_4 \cosh ms + \frac{2\bar{p}\pi^2\epsilon}{\omega^2} \cos 2\pi s \qquad (14) $$

with $n = (\sqrt{\omega^2 + 4\pi^4} + 2\pi^2)^{1/2}$ and $m = (\sqrt{\omega^2 + 4\pi^4} - 2\pi^2)^{1/2}$. Boundary conditions (6) yield the solvability condition

$$ 0 = 8(m^2 + n^2)\pi^4\epsilon^2 R_1(n, m) - 2mn(2\pi^4\epsilon^2 - \eta\omega^2)R_2(n, m) \qquad (15) $$

with $R_1(n, m) = n(\cosh m - 1)\sin n + m(\cos n - 1)\sinh m$ and $R_2(n, m) = mn(\cos n \cosh m - 1) + 2\pi^2 \sin n \sinh m$, which is in fact an equation for the frequency $\omega$. We plot $p_E$, $y_E(1/2)$, and $\omega$ in Figures 2, 3, and 4 to compare with the results from Kirchhoff's model.

As reported in (Nayfeh et al. 1995; Mamou-Mani et al. 2009), we remark in Figure 5 that the second vibration mode in this model has an angular frequency that does not depend on $d$: it is constant throughout the whole post-buckling regime. All the even vibration modes share this property. For these modes, Equation (15) is in fact fulfilled through a common zero of the functions $R_1(\omega)$ and $R_2(\omega)$. These common zeros do not depend on $\epsilon$, hence do not depend on $d$. See Section B for a study of the common zeros of $R_1(\omega)$ and $R_2(\omega)$.

Finally, in Section E, we re-plot Figures 2, 3, 4, and 5 with the load $p_E$ (instead of $D/L$) on the horizontal axis.
# 4 Models comparison and validity range of the Woinowsky-Krieger model

## 4.1 Limit in term of axial displacement

We first remind that, as explained in Section 2 and Section 3, with suitable choices of dimensionless variables, the Woinowsky-Krieger model does not depend on any parameter, but the Kirchhoff (reference) model depends on the slenderness ratio $\eta$. This is illustrated in Figure 6 which shows the plots of Figures 2, 3, 4, and 5 with rescaled axes. As expected, on each of the four plots the

<figure><img src="image_7.png" /><figcaption>**Figure 6** Comparison of the Woinowsky-Krieger and Kirchhoff models for various values of the parameter $\eta \in \{1/1200, 1/4800, 1/19200, 1/76800\}$, respectively corresponding to $L/h \in \{10, 20, 40, 80\}$ for a rectangular cross-section.</figcaption></figure>

curves for the Woinowsky-Krieger model fall into a single master curve, while the curves for the Kirchhoff model are seen to depend on $\eta$. For the transverse displacement at midspan $Y_E(L/2)$ and the natural frequency $\omega$ of the first mode, the Woinowsky-Krieger curves agree with the Kirchhoff curves until limiting values, that depend on $\eta$. In Figure 7, we plot the curves for $Y_E(L/2)$ and $\omega$ as functions of $d = D/L$, and the relative error of Woinowsky-Krieger model (as compared to the Kirchhoff model) for several values of the slenderness ratio $\eta$. It is observed that all relative error curves are almost superimposed, meaning that the error as a function of $d = D/L$ is almost independent of $\eta$, see also Section F. This result is interesting since most of the literature about geometrical nonlinearities traditionally gives the validity limit of the Woinowsky-Krieger model in term of $\hat{y}_E(1/2) = Y_E(L/2)/r = Y_E(L/2)/(L\sqrt{\eta})$, usually around $\hat{y}_E(1/2) \approx 2$ or 3. Here, we prove that a correct validity limit should be given in terms of $y_E(L/2) = Y_E(L/2)/L$ (or in terms of $d = D/L$) and thus that the correct scaling of $Y$ for the validity limit is not the radius of gyration $r$ but the length $L$ of the beam. In particular, an axial displacement of $D/L = 0.1$ (or equivalently a transverse displacement $Y_E(L/2)/L \approx 0.2$), gives errors of less than 4% on $Y$ and less than 3% on $\omega$. We then conclude that the Woinowsky-Krieger model tends to depart from the reference Kirchhoff model as soon as 10% of axial displacement or 20% of transverse displacement, and this for any slenderness ratio. Equivalently, in terms of $\hat{y}_E(1/2)$, for $\eta \in \{1/1200, 1/4800, 1/19200, 1/76800\}$ (that is $L/h \in \{10, 20, 40, 80\}$ in the case of a rectangular cross-section), the limit of the Woinowsky-Krieger model is then $\hat{y}_E(1/2) \approx \{7, 14, 28, 55\}$, which is larger than $\hat{y}_E(1/2) \approx 2$ or 3, as often claimed.

Another interesting result is that these errors are roughly linear functions of $D/L$, as long as $D/L < 0.5$. Finally, we note that, for each of the plots in Figure 7, one of the curves lies slightly apart from the others. This curve is associated with the largest $\eta = 1/1200$, corresponding to a beam with a thickness only ten times smaller than the length ($L/h = 10$). Such a case has to be considered with caution since the validity of the Euler-Bernoulli kinematics is then questionable.
<figure><img src="image_8.png" /><figcaption>**Figure 7** Comparison of the Woinowsky-Krieger and Kirchhoff models for various values of the parameter $\eta \in \{1/1200, 1/4800, 1/19200, 1/76800\}$ corresponding to $L/h \in \{10, 20, 40, 80\}$, respectively, for a rectangular cross section. Top: transverse displacement and frequency of the first mode. Bottom: relative error of the Woinowsky-Krieger model (compared to the Kirchhoff model), for the transverse displacement and the frequency of the first mode. Relative errors defined as $Err_y = (y_{WK}(1/2) - y_{Ki}(1/2))/y_{Ki}(1/2)$ and $Err_\omega = (\omega_{Ki}(1/2) - \omega_{WK}(1/2))/\omega_{Ki}(1/2)$.</figcaption></figure>

## 4.2 Offset in the critical load

Moreover, looking at the results for small $d$, we detect an offset between Woinowsky-Krieger and Kirchhoff curves, see insets in Figures 2, 3, 4 and 5. This offset exists right from buckling and to disclose it analytically we proceed to construct a series expansion in powers of $\epsilon$, a small parameter measuring the mean transverse displacement of the buckled solution (Arbocz et al. 1985). We find it convenient to use

$$ \epsilon = -4 \int_{0}^{1} y_{E}(s) \cos 2\pi s \, ds. \qquad (16) $$

For the Woinowsky-Krieger model, $\epsilon$ is exactly equal to $y_E(L/2) = Y_E(L/2)/L$, see Equation (13a). We inject the following expansions

$$ x_E(s) = x_{E0}(s) + \epsilon x_{E1}(s) + \epsilon^2 x_{E2}(s) + \epsilon^3 x_{E3}(s) + O(\epsilon^4) \qquad (17a) $$

$$ y_E(s) = \epsilon y_{E1}(s) + \epsilon^2 y_{E2}(s) + \epsilon^3 y_{E3}(s) + O(\epsilon^4) \qquad (17b) $$

$$ \theta_E(s) = \epsilon \theta_{E1}(s) + \epsilon^2 \theta_{E2}(s) + \epsilon^3 \theta_{E3}(s) + O(\epsilon^4) \qquad (17c) $$

$$ -n_{xE} = p_E = p_{E0} + \epsilon p_{E1} + \epsilon^2 p_{E2} + \epsilon^3 p_{E3} + O(\epsilon^4) \qquad (17d) $$

$$ d = d_0 + \epsilon d_1 + \epsilon^2 d_2 + \epsilon^3 d_3 + O(\epsilon^4) \qquad (17e) $$

into the Kirchhoff equilibrium system, and solve the equations at each order of $\epsilon$. At order $\epsilon^0$, we find $x_{E0}(s) = s(1 - \eta p_{E0})$, with $p_{E0}$ still unknown. At order $\epsilon^1$, we solve

$$ x'_{E1} = -\eta p_{E1} \quad \text{with} \quad x_{E1}(0) = 0 \qquad (18a) $$

$$ y'_{E1} = (1 - \eta p_{E0}) \theta_{E1} \quad \text{with} \quad y_1(0) = 0 = y_1(1) \qquad (18b) $$

$$ \theta''_{E1} = -p_{E0}(1 - \eta p_{E0}) \theta_{E1} \quad \text{with} \quad \theta_1(0) = 0 = \theta_1(1) \qquad (18c) $$

and find that the solution is

$$ x_{E1} = -\eta p_{E1}s \qquad (19a) $$

$$ y_{E1} = \frac{1}{2}(1 - \cos 2\pi s) \Rightarrow \theta_{E1} = \frac{\pi \sin 2\pi s}{1 - \eta p_{E0}} \qquad (19b) $$

$$ p_{E0} = \frac{1 - \sqrt{1 - 16\pi^2\eta}}{2\eta} \qquad (19c) $$
where we see that $p_{E1}$ is not defined at order $\epsilon^1$ and will only be selected when solving order $\epsilon^2$. In the same manner, finding $p_{E2}$ requires to solve order $\epsilon^3$. We shall come back to this remark when discussing the Woinowsky-Krieger model. We note that the same shift occurs for $x_E(s)$ (hence $d$) and $\omega$. Conducting the expansion up to order $\epsilon^3$ (included) yields

$$p_E = p_{E0} + \epsilon^2 \frac{p_{E0}^3 (16\pi^2 - 3p_{E0})}{128\pi^2 (8\pi^2 - p_{E0})} + p_{E3}\epsilon^3 + O(\epsilon^4) \quad (20a)$$

$$= 4\pi^2 + 16\pi^4\eta + O(\eta^2) + [\frac{1}{2}\pi^4 + 2\pi^6\eta + O(\eta^2)]\epsilon^2 + p_{E3}\epsilon^3 + O(\epsilon^4) \quad (20b)$$

$$y_E(1/2) = \epsilon + p_{E0}^2 \frac{16\pi^2 - 3p_{E0}}{4096\pi^4} \epsilon^3 + O(\epsilon^4) \quad (20c)$$

$$= \epsilon + \frac{1}{64}\pi^2 [1 - 4\pi^2\eta + O(\eta^2)]\epsilon^3 + O(\epsilon^4) \quad (20d)$$

$$d = \eta p_{E0} + \epsilon^2 p_{E0} \frac{2\pi(p_{E0}^2 + 4p_{E0}\pi^2 - 64\pi^4)}{256\pi^3(8\pi^2 - p_{E0})} + d_3\epsilon^3 + O(\epsilon^4) \quad (20e)$$

$$= 4\pi^2\eta + O(\eta^2) + \frac{1}{4}\pi^2[1 + 2\pi^2\eta + O(\eta^2)]\epsilon^2 + d_3\epsilon^3 + O(\epsilon^4). \quad (20f)$$

The solutions $x_E(s)$, $y_E(s)$, and $\theta_E(s)$ are listed in Section A. Please note that, as explained in Section A, further calculations lead to $p_{E3} = 0$ and $d_3 = 0$. A similar expansion for the first mode of vibration (with $\omega_0 = 0$) leads to

$$\omega = \epsilon \frac{p_0}{8\pi^2} \sqrt{\frac{64\pi^4 - p_0^2 - 4\pi^2 p_0}{12\eta}} + \omega_3\epsilon^3 + O(\epsilon^4) \quad (21a)$$

$$= \epsilon\pi^2 \sqrt{\frac{2}{3\eta}} [1 + O(\eta)] + \omega_3\epsilon^3 + O(\epsilon^4). \quad (21b)$$

Note that the results noted 'Ki.dev 4' in the figures involve expansions of the solutions up to order $\epsilon^4$ included. These expansions are given in the supplementary materials (Neukirch et al. 2021).

Comparing Woinowsky-Krieger and Kirchhoff solutions for $p_E$, $y_E(1/2)$, and $d$, listed in Equation (13) and Equation (20), reveals the nature of the offset mentioned earlier: as soon as order $\epsilon^0$, $p_E$ is not computed exactly in the Woinowsky-Krieger model which, strictly speaking, is then not a rigorous expansion of Kirchhoff's model. In this matter, the Woinowsky-Krieger model is equivalent to the traditional linear beam model — see e.g. (Bažant and Cedolin 1991) — and predicts the critical buckling load to be the classical Euler value $p_E = 4\pi^2$, see Equation (13c). On the contrary, Kirchhoff's model predicts this buckling load to be $p_E = 4\pi^2(1 + 4\pi^2\eta + O(\eta^2))$, see Equation (20b) with $\epsilon = 0$. The offset, though non-zero, is in most practical cases negligible as $\eta \ll 1$ for slender beams. The physical explanation of this offset is that the Woinowsky-Krieger model neglects the small shortening of the beam before the critical load. Following Equation (3a), this shortening is $(1 + e)$ for an infinitesimal axial element when $\theta \ll 1$ (that is $dx = (1 + e) ds$). This axial shortening is neglected in the Woinowsky-Krieger model which writes Equation (3c) as $y' = \theta$. To (artificially) correct the Woinowsky-Krieger model, one could use $y' = (1 + e)\theta$ and keep all other approximations. In doing this, Equation (3e) would become $m = y''/(1 + e)$ and since $e = -\eta p$, this would lead to replace $p$ by $p(1 - \eta p)$ in Equation (11a). Solving this equation would yield $p_E(1 - \eta p_E) = 4\pi^2$, which is the exact buckling load, Equation (19c). Please note nevertheless that it would not cure all Woinowsky-Krieger shortcomings.

Finally, we can also compare the two models in their prediction of the curvature of the curve $y_E(L/2) = f(d)$ just after the buckling point. For the Kirchhoff model, combining Equation (20d) and Equation (20f), we have

$$d \approx 4\pi^2\eta + y_E^2(1/2)\frac{1}{4}\pi^2[1 + 2\pi^2\eta] \quad \text{for small } \eta \text{ and } y_E \quad (22)$$

while for the Woinowsky-Krieger model, combining Equation (13a) and Equation (13b), we have

$$d = 4\pi^2\eta + \frac{1}{4}\pi^2 y_E^2(1/2). \quad (23)$$

Here also, the Woinowsky-Krieger model is wrong by a small term, proportional to $\eta$.
## 4.3 Second order in the axial load

There is yet another, more important, flaw in Woinowsky-Krieger approach: there is no order $\epsilon^2$ in the solution for $p_E$. Indeed, computing $p_{E2}$ would require an order 3 in its equation for transverse displacement (11a). In this sense, the Woinowsky-Krieger model does not yield a proper order 2 expansion of Kirchhoff's solutions. This implies that the Woinowsky-Krieger model predicts a constant load $p_E$ along the post-buckling path, see Figure 2, and is therefore unable to deal with load-controlled experiments. Nevertheless, in Equation (11b) the load $p_E$ is multiplied by the small parameter $\eta$, which tends to weaken the absence of the $p_{E2}$ term.

Another way to test for the order 2 conformity of the Woinowsky-Krieger equation is to take Kirchhoff's solutions $x_E(s)$, $y_E(s)$, and $p_E$ (see Section A) and inject them in Equation (11). We obtain

$$y_E''' + p_E y_E'' = 2\pi^2 p_{E0}^2 \eta \epsilon \cos 2\pi s + 0 \epsilon^2 + O(\epsilon^3) \quad (24a)$$

$$x_E' - 1 + \eta p_E + \frac{1}{2} y_E'' = \frac{\eta^2 p_{E0}^4}{32\pi^2} \epsilon^2 \sin^2 2\pi s + 0 \epsilon^3 + O(\epsilon^4) \quad (24b)$$

where we see that, strictly speaking, the transverse displacement equation is not fulfilled at order $\epsilon^1$ and the axial displacement equation is not fulfilled at order $\epsilon^2$.

For the vibrations, the picture is very much the same: if we make an expansion of the solution of Equation (15) in powers of $\epsilon$, we find

$$\omega = \epsilon\pi^2 \sqrt{\frac{2}{3\eta}} + 0 \epsilon^2 + O(\epsilon^3) \quad (25)$$

which misses an $O(\eta)$ term, but is the leading $\eta$ order of Kirchhoff's result, see Equation (21b).

# 5 An order 3 model

The motivation for a third order model is theoretical: following up on the remark of the previous section about the necessity of having an order 3 in the transverse displacement equation (8a), we set on investigating what precisely third order terms are bringing to the solution. Doing this, we keep in mind the elastic plates application: the (second order) Von Kármán assumption (7) is widely used in plate models and one is tempted to know what would third order terms bring to such plate models, and how they would compare to geometrically-exact approaches, as for example Naghdi's model (Naghdi 1963; Naghdi 1964).

We start from Kirchhoff's system (3) and proceed to derive an order 3 model. Owing to the remark that the parameter $\eta$ is small and that the Woinowsky-Krieger model is only exact in the limit $\eta \to 0$, we settle on removing $\eta$ as much as possible since this makes derivation much easier. Nevertheless, it has been shown (Neukirch et al. 2012; Neukirch et al. 2014) that in displacement-controlled experiments, $\eta$ should at least remain in the axial displacement of Equation (3a). This equation is then simplified to

$$x' = 1 + \eta n_x - \frac{1}{2} \theta^2. \quad (26)$$

This resembles Woinowsky-Krieger in Equation (9b), with the difference that we keep the s dependence in the load $n_x(s, t)$. For the transverse displacement, we know extension is only playing a minor role and we readily set $\eta = 0$. We then start with $y' = \sin \theta$ and proceed to develop the sinus up to order 3, $y' = \theta - \frac{1}{6}\theta^3$. Inverting this relation yields

$$\theta = y' + \frac{1}{6}y'^3 \quad \Rightarrow \quad \theta'' = y''' + y'y''^2 + \frac{1}{2}y'^2y''' \quad (27)$$

which is injected into Equation (3f). We end up with a system for the axial and transverse displacement $(x, y)$ and forces $(n_x, n_y)$:

$$x'(s, t) = 1 + \eta n_x(s, t) - \frac{1}{2} y'^2(s, t) \quad (28a)$$

$$n_x'(s, t) = \ddot{x}(s, t) \quad (28b)$$

$$y(s, t)''' + f_3 = n_x(s, t)y'(s, t) - n_y(s, t)x'(s, t) \quad (28c)$$

$$n_y'(s, t) = \ddot{y}(s, t) \quad (28d)$$
with $f_3 = y'(s, t)y''(s, t)^2 + \frac{1}{2}y'(s, t)^2y'''(s, t)$. To readily compare this new system with the Woinowsky-Krieger model, one has to differentiate Equation (28c) and find

$$y'''' - n_x y'' + \ddot{y} = -y''' - 3y'y''y''' - \frac{1}{2}y'^2y'''' + \ddot{x}y' - n_x x'' - \ddot{y}(\eta n_x - \frac{1}{2}y'^2) \quad (29)$$

which, in contrast to Equation (8a), has order 3 terms in the right-hand side.

Equilibrium equations for $x_E$, $y_E$, $n_{xE}$, $n_{yE}$ are obtained by setting $\ddot{x} = 0 = \ddot{y}$ in Equation (28) and vibration equations are then derived by injecting Equation (4) into Equation (28). We obtain

$$\bar{x}' = \eta \bar{n}_x - y_E \bar{y}' \quad (30a)$$

$$\bar{n}_x' = -\omega^2 \bar{x} \quad (30b)$$

$$\bar{y}''' + \bar{f}_3 = \bar{n}_x y' + n_{xE} \bar{y}' - \bar{n}_y x' - n_{yE} \bar{x}' \quad (30c)$$

$$\bar{n}_y' = -\omega^2 \bar{y} \quad (30d)$$

with $\bar{f}_3 = \bar{y}'\bar{y}''^2 + 2\bar{y}_E'\bar{y}''\bar{y}_E'' + \bar{y}'\bar{y}_E'\bar{y}''' + \frac{1}{2}\bar{y}_E'^2\bar{y}''''$. We solve this system numerically and compare the results to the two previous models in Figures 2, 3, 4, and 5. A first remark is that this new model suffers from the same offset as the Woinowsky-Krieger model: right from buckling a small shift exists in the curves. As explained earlier it arises from the setting of $\eta = 0$ in the transverse displacement equation. Next we see that the load curve is no longer flat (Figure 2), nor is the frequency of the 2nd vibration mode (Figure 4).

In order to make sure we indeed came up with a model exhibiting the correct terms up to order 3, we perform the expansion (17) and compute

$$y_E(s) = \frac{1}{2}\epsilon(1 - \cos 2\pi s) + \epsilon^3 \frac{1}{64}\pi^2 \sin^2 3\pi s + O(\epsilon^4) \quad (31a)$$

$$x_E(s) = s(1 - 4\pi^2\eta) - \frac{1}{16}\pi(4\pi s + 8\pi^3\eta s - \sin 4\pi s)\epsilon^2 + O(\epsilon^4) \quad (31b)$$

$$p_E = 4\pi^2 + \frac{1}{2}\pi^4\epsilon^2 + O(\epsilon^4) \quad (31c)$$

$$d = 1 - x_E(1) = 4\pi^2\eta + \frac{1}{4}\pi^2(1 + 2\pi^2\eta)\epsilon^2 + O(\epsilon^4) \quad (31d)$$

$$y_E(1/2) = \epsilon + \frac{1}{64}\pi^2\epsilon^3 + O(\epsilon^4) \quad (31e)$$

which, indeed, is correct up to order 3 when compared to Kirchhoff's results, see Equation (20) and Section A. For the first mode of vibration (with $\omega_0 = 0$) we find

$$\omega = \epsilon\pi^2 \sqrt{\frac{2(1 + 2\pi^2\eta)}{3\eta(1 - 4\pi^2\eta)}} + O(\epsilon^3) \quad (32)$$

which corresponds to the first $\eta$ order of Equation (21). See the supplementary material (Neukirch et al. 2021) for detailed calculations.

The interest of the present order 3 model lies in the fact that it efficiently corrects the Woinowsky-Krieger model for the axial load and the second natural frequency, in the small $\eta$ limit. Moreover, it is the extension of a well known order 3 model, commonly used for nonlinear vibrations of inextensible cantilever beams and first introduced by Crespo da Silva and Glynn (1978)—see (Thomas et al. 2016) for a list of other references: as shown in Section D, the present order 3 model reduces to the Crespo da Silva model in the case of clamped-free boundary conditions and inextensible beams. However, contrary to the Crespo da Silva model, which elegantly involves a single equation for the unique variable $y(s, t)$, the present order 3 model consists of a system of four equations in four variables, system (28), which might be complex to use in practice.

# 6 Conclusion

We have studied the range of validity of the Woinowsky-Krieger equations for the planar equilibrium and vibrations of post-buckled beams. The Woinowsky-Krieger equations are useful and widely used, especially when dealing with nonlinear vibrations, but are only valid in the weakly nonlinear regime and under displacement-controlled setups. We have shown that these equations are not a rigorous second order development of Kirchhoff's equations, but that they
nevertheless capture faithfully the post-buckling behavior of the beam up to 10% ($D = 0.1L$) of axial displacement and/or 20% of transverse displacement ($Y(L/2) = 0.2L$), and that these limits only weakly depend on the slenderness ratio of the beam. If the transverse displacement $Y(L/2)$ is written in units of the beam thickness $h$, we have shown that the validity limit of the Woinowsky-Krieger model then depends on the slenderness ratio of the beam and that it can be large: $Y(L/2)/h < 4$ for a thickness to length ratio of $h/L = 1/20$ and $Y(L/2)/h < 16$ for $h/L = 1/80$. Incidentally, we have also rigorously proved that every other vibration frequency in the Woinowsky-Krieger model is load-independent in the entire post-buckling regime. Finally, we have introduced a third order model capable of coping with load-controlled setups and more accurately predicting vibration modes in the moderate post-buckling regime.

## A Expansions for the solution of the Kirchhoff model

In the clamped-clamped case, buckling happens in a symmetrical pitchfork bifurcation. Consequently, with the chosen definition of $\epsilon$ in Equation (16), the developments of the axial variables $x_E$, $p_E = -n_{xE}$, and $e_E$ only comprise even terms in $\epsilon$, while the developments of the transverse variables $y_E$, $\theta_E$, $m_E$, and $n_{y_E}$ only comprise odd terms in $\epsilon$. Please see the supplementary material (Neukirch et al. 2021) for detailed calculations leading to

$$ 
\begin{aligned}
\theta_E(s) &= \epsilon \frac{p_0}{4\pi} \sin 2\pi s + \epsilon^3 \frac{p_0^3(16\pi^2 - 3p_0)(96\pi^2 \sin 2\pi s - (8\pi^2 - p_0) \sin 6\pi s)}{48(4\pi)^5(8\pi^2 - p_0)} + O(\epsilon^5) \\
&= \epsilon\pi \sin 2\pi s [1 + 4\pi^2\eta + O(\eta^2)] + \epsilon^3\pi^3 \frac{25 + 2\cos 4\pi s + 96\pi^2\eta + O(\eta^2)}{192} \sin 2\pi s + O(\epsilon^5)
\end{aligned}
$$

$$ 
\begin{aligned}
y_E(s) &= \frac{\epsilon}{2}(1 - \cos 2\pi s) + \epsilon^3 p_0^2 \frac{16\pi^2 - 3p_0}{4096\pi^4} \sin^2(3\pi s) + O(\epsilon^5) \\
&= \frac{\epsilon}{2}(1 - \cos 2\pi s) + \epsilon^3 \frac{\pi^2}{64} [1 - 4\pi^2\eta + O(\eta^2)] \sin^2(3\pi s) + O(\epsilon^5)
\end{aligned}
$$

$$ 
\begin{aligned}
x_E(s) &= s(1 - \eta p_0) - \epsilon^2 p_0 \frac{2\pi(p_0^2 + 4p_0\pi^2 - 64\pi^4)s + (8\pi^2 - p_0)^2 \sin 4\pi s}{256\pi^3(8\pi^2 - p_0)} + O(\epsilon^4) \\
&= [s - 4\pi^2\eta s + O(\eta^2)] + [\frac{\pi}{16}(\sin 4\pi s - 4\pi s) - \frac{\pi^4}{2}s\eta + O(\eta^2)]\epsilon^2 + O(\epsilon^4).
\end{aligned}
$$

## B Common zeros of $R_1$ and $R_2$

We replace $n = 2\pi\beta$, $m = 2\pi\alpha$, $R_1 = 4\pi\hat{R}_1$, and $R_2 = 4\pi^2\hat{R}_2$ to obtain

$$ \hat{R}_1(\alpha, \beta) = \beta \sin 2\pi\beta \sinh^2 \pi\alpha - \alpha \sin^2 \pi\beta \sinh 2\pi\alpha \quad (B.1a) $$

$$ \hat{R}_2(\alpha, \beta) = 2\alpha\beta(\cos 2\pi\beta \cosh 2\pi\alpha - 1) + \sin 2\pi\beta \sinh 2\pi\alpha \quad (B.1b) $$

$$ \beta^2 = 1 + \alpha^2 \quad (B.1c) $$

with $\beta = \frac{1}{2\pi}(\sqrt{\omega^2 + 4\pi^4} + 2\pi^2)^{1/2}$ and $\alpha = \frac{1}{2\pi}(\sqrt{\omega^2 + 4\pi^4} - 2\pi^2)^{1/2}$. We work with $n > 2\pi$ and $m > 0$, that is

$$ \beta > 1 \text{ and } \alpha > 0. \quad (B.2) $$

In this section, we show that

1. $\hat{R}_1(\omega)$ and $\hat{R}_2(\omega)$ have common zeros,

2. but also have separate zeros, with

2(a) $\hat{R}_1(\omega) = 0$ and $\hat{R}_2(\omega) \neq 0$ when $\beta = 2, 3, 4, \dots$,

2(b) $\hat{R}_2(\omega) = 0$ and $\hat{R}_1(\omega) \neq 0$ when $\hat{A} - 1/\hat{A} = \hat{B} - 1/\hat{B}$ with $\hat{A} \neq \hat{B}$, where $\hat{A}$ and $\hat{B}$ are defined in Section B.3.

### B.1 Individual zeros of $\hat{R}_1$

To prove 2(a), we factorize $\hat{R}_1$ as $\hat{R}_1 = 2 \sinh \pi\alpha \sin \pi\beta(y \cos \pi\beta \sinh \pi\alpha - \alpha \sin \pi\beta \cosh \pi\alpha)$ and see, from Equation (B.2), that $\hat{R}_1(\omega) = 0$ for $\beta = 2, 3, 4, \dots$. In such cases, $\hat{R}_2 = 2\alpha\beta(\cosh 2\pi\alpha - 1) > 0$, hence we have zeros of $\hat{R}_1$ which are not zeros of $\hat{R}_2$.
## B.2 Common zeros of $\hat{R}_1$ and $\hat{R}_2$

We first remark that the zeros of $\hat{R}_1$ or $\hat{R}_2$ are such that $\cos \pi\beta \neq 0$: If $\beta = 3/2, 5/2, 7/2, \dots$ then $\hat{R}_1 = -\alpha \sinh 2\pi\alpha \neq 0$ and $\hat{R}_2 = -2\alpha\beta(\cosh 2\pi\alpha + 1) \neq 0$. We may then divide by $\alpha, \beta, \cos \pi\beta$, and $\cosh \pi\alpha$ without any trouble and rewrite

$$ \hat{R}_1(\alpha, \beta) = 2 \sin \pi\beta \cos \pi\beta \sinh \pi\alpha \cosh \pi\alpha (\beta \tanh \pi\alpha - \alpha \tan \pi\beta) \quad (B.3a) $$

$$ \hat{R}_2(n, m) = 4 \cos^2 \pi\beta \cosh^2 \pi\alpha (\alpha\beta \tanh^2 \pi\alpha - \alpha\beta \tan^2 \pi\beta + \tanh \pi\alpha \tan \pi\beta). \quad (B.3b) $$

We then see that if $\beta \neq 2, 3, 4, \dots$ and $\hat{R}_1 = 0$, we have $\beta \tanh \pi\alpha = \alpha \tan \pi\beta$ and $\alpha\beta \tanh^2 \pi\alpha - \alpha\beta \tan^2 \pi\beta + \tanh \pi\alpha \tan \pi\beta = 0$, i.e. $\hat{R}_2 = 0$. This proves 1.

## B.3 Individual zeros of $\hat{R}_2$

If we have $\hat{R}_2 = 0$ and $\beta \neq 2, 3, 4, \dots$, then, using Equation (B.1c)

$$ \frac{\tan \pi\beta}{\tanh \pi\alpha} - \frac{\tanh \pi\alpha}{\tan \pi\beta} = \frac{\beta}{\alpha} - \frac{\alpha}{\beta} \quad (B.4) $$

which is $\hat{A} - 1/\hat{A} = \hat{B} - 1/\hat{B}$ with $\hat{A} = \tan \pi\beta/\tanh \pi\alpha$ and $\hat{B} = \beta/\alpha$. Solutions are either $\hat{A} = \hat{B}$ is which case we have a common zero of $\hat{R}_1$ and $\hat{R}_2$, or solutions with $\hat{A} < 0$ and $\hat{B} > 0$ in which case we have a zero of $\hat{R}_2$ such that $\hat{R}_1 \neq 0$. This proves 2(b).

## B.4 Summary

In Table B.1, we see that each common zero is followed by an individual zero of $\hat{R}_2$, then by an individual zero of $\hat{R}_1$.

<table><caption>**Table B.1** Lowest six zeros of $\hat{R}_1$ and $\hat{R}_2$. Please note that $\omega = 4\pi^2\alpha\beta$, but that only common roots correspond to the actual vibration frequencies of the beam.</caption><thead><tr><th>ω</th><th>α</th><th>β</th><th>R̂<sub>1</sub></th><th>R̂<sub>2</sub></th><th>Â</th><th>Б̂</th><th>Â - 1/Â</th><th>Б̂ - 1/Б̂</th></tr></thead><tbody><tr><td>44.4</td><td>0.85</td><td>1.31</td><td>0</td><td>0</td><td>1.54</td><td>1.54</td><td>0.89</td><td>0.89</td></tr><tr><td>103.5</td><td>1.47</td><td>1.78</td><td>-7.6 · 10<sup>3</sup></td><td>0</td><td>-0.83</td><td>1.21</td><td>0.38</td><td>0.38</td></tr><tr><td>136.8</td><td>1.73</td><td>2</td><td>0</td><td>9.2 · 10<sup>4</sup></td><td>0</td><td>1.15</td><td>∞</td><td>0.29</td></tr><tr><td>182.1</td><td>2.03</td><td>2.27</td><td>0</td><td>0</td><td>1.11</td><td>1.11</td><td>0.22</td><td>0.22</td></tr><tr><td>280.6</td><td>2.57</td><td>2.76</td><td>-1.4 · 10<sup>7</sup></td><td>0</td><td>-0.93</td><td>1.07</td><td>0.14</td><td>0.14</td></tr><tr><td>334.9</td><td>2.83</td><td>3</td><td>0</td><td>2.2 · 10<sup>8</sup></td><td>0</td><td>1.06</td><td>∞</td><td>0.12</td></tr></tbody></table>

# C Energies for the various models

We present a variational approach for the three different models used in this paper. We list the kinetic $\mathcal{T}$ and potential $\mathcal{V}$ energies and compute the first variation of the Action $S = \int_{t_1}^{t_2} \mathcal{L} dT$ where the Lagrangian $\mathcal{L} = \mathcal{T} - \mathcal{V}$.

## C.1 The Kirchhoff model

In this model the kinetic energy is computed as if the mass of each section were concentrated on the centerline, that is no rotational inertia is involved. We have

$$ \mathcal{T} = \frac{1}{2} \rho A \int_{0}^{L} (\dot{X}^2 + \dot{Y}^2) \, dS. \quad (C.1) $$

The potential energy comprises bending and extension deformations

$$ \mathcal{V} = \frac{1}{2} \int_{0}^{L} (EI\theta'^2 + EAe^2) \, dS \quad (C.2) $$

with boundary conditions (6) valid at all time $T$,

$$ X(0, T) = 0 = X(L, T) - L + D \quad (C.3a) $$

$$ Y(0, T) = 0 = Y(L, T) \quad (C.3b) $$

$$ \theta(0, T) = 0 = \theta(L, T). \quad (C.3c) $$
The Action $\mathcal{S}$ is then a functional of $q = (X, Y, \theta, e)$ and the principle of least Action then selects the dynamical evolution $q(S, T)$ of the system. This principle reads

$$ \mathcal{S}(q + \epsilon \bar{q}) \geq \mathcal{S}(q) \quad \text{for all small } \epsilon \text{ and for all admissible } \bar{q} \qquad (C.4) $$

under the pointwise kinematic constraints

$$ \phi_x = X' - (1 + e) \cos \theta = 0 \quad \text{and} \quad \phi_y = Y' - (1 + e) \sin \theta = 0. \qquad (C.5) $$

We introduce Lagrange multipliers to deal with constraints (C.5). We anticipate the multipliers to be the components $N_x$ and $N_y$ of the force vector and use $\mathcal{L} = \mathcal{T} - \mathcal{V} - N_x \phi_x - N_y \phi_y$ as the Lagrangian. The first order necessary condition for Equation (C.4) to hold is then

$$ \bar{\mathcal{S}}(q, \bar{q}) = \lim_{\epsilon \to 0} \frac{\mathcal{S}(q + \epsilon \bar{q}) - \mathcal{S}(q)}{\epsilon} = 0 \quad \text{for all admissible } \bar{q} \qquad (C.6) $$

with

$$ -\bar{\mathcal{S}} = \int_{t_1}^{t_2} \int_{0}^{L} \left\{ -\rho A (\dot{X}\bar{X} + \dot{Y}\bar{Y}) + EI\theta'\bar{\theta}' + EAe\bar{e} + N_x(\bar{X}' - \bar{e}\cos\theta + (1+e)\bar{\theta}\sin\theta) + N_y(\bar{Y}' - \bar{e}\sin\theta - (1+e)\bar{\theta}\cos\theta) \right\} dS \, dT. \qquad (C.7) $$

Using Equation (C.3), boundary conditions for the test functions $\bar{q}$ read

$$ \bar{X}(0, T) = 0 = \bar{X}(L, T) \qquad (C.8a) $$

$$ \bar{Y}(0, T) = 0 = \bar{Y}(L, T) \qquad (C.8b) $$

$$ \bar{\theta}(0, T) = 0 = \bar{\theta}(L, T). \qquad (C.8c) $$

Using Equation (C.8) and $\bar{X}(S, t_1) = 0 = \bar{X}(S, t_2)$, $\bar{Y}(S, t_1) = 0 = \bar{Y}(S, t_2)$ for all $S$, we perform integration by parts on $S$ and $T$ to find

$$ -\bar{\mathcal{S}} = \int_{t_1}^{t_2} \int_{0}^{L} \left\{ \rho A (\ddot{X}\bar{X} + \ddot{Y}\bar{Y}) - EI\theta''\bar{\theta} + EAe\bar{e} - N_x'\bar{X} - N_x\bar{e}\cos\theta + N_x(1+e)\bar{\theta}\sin\theta - N_y'\bar{Y} - N_y\bar{e}\sin\theta - N_y(1+e)\bar{\theta}\cos\theta \right\} dS \, dT. \qquad (C.9) $$

Finally, imposing that $\bar{\mathcal{S}} = 0$ for all test functions $\bar{X}(S, T)$, $\bar{Y}(S, T)$, $\bar{\theta}(S, T)$, and $\bar{e}(S, T)$ yields System (3).

## C.2 The Woinowsky-Krieger model

See also (Younis and Nayfeh 2003) for an energy approach of the Woinowsky-Krieger equations. In this model, the kinetic energy only comprises the transverse displacement term

$$ \mathcal{T} = \frac{1}{2} \rho A \int_{0}^{L} \dot{Y}^2 dS. \qquad (C.10) $$

Moreover, a linear formula is used for the curvature in the potential energy

$$ \mathcal{V} = \frac{1}{2} \int_{0}^{L} (EIY''^2 + EAe^2) dS \qquad (C.11) $$

and the kinematic constraint (C.5) reads

$$ \phi_x = U' + \frac{1}{2}Y'^2 - e = 0 \qquad (C.12) $$

to which is associated a continuous Lagrange multiplier which we call $N_x(S, T)$. The Action $\mathcal{S}$ is then a functional of $q = (Y, U, e)$ and the first variation reads, after several integrations by parts,

$$ -\bar{\mathcal{S}} = \int_{t_1}^{t_2} \int_{0}^{L} \{ \rho A \ddot{Y} \bar{Y} + EIY''' \bar{Y} + EAe\bar{e} - N_x' \bar{U} - N_x \bar{e} - N_x' Y' \bar{Y} - N_x Y'' \bar{Y} \} dS \, dT. \qquad (C.13) $$

Imposing $\bar{\mathcal{S}} = 0$ for all test functions $\bar{U}(S, T)$, $\bar{Y}(S, T)$, and $\bar{e}(S, T)$ yields System (8) with $P(T) = -N_x(T)$.
### C.3 The order 3 model

The kinetic and potential energies take the same form as in the Kirchhoff model, see Equation (C.1) and Equation (C.2). As we aim at an order 3 model, we need to use expansions at order 4 for the kinematic constraints (C.5)

$$X' = 1 + e^{-\frac{1}{2}\theta^2} + \frac{1}{24}\theta^4 \qquad (C.14a)$$

$$Y' = \theta - \frac{1}{6}\theta^3 \qquad (C.14b)$$

where, as explained in Section 5, we only kept the lowest order in the small extension $e$ limit in Equation (C.14a) and completely removed it from Equation (C.14b). The Action $S$ is then a functional of $q = (X, Y, \theta, e)$ and the first variation reads, after integration by parts,

$$-\bar{S} = \int_{t_1}^{t_2} \int_{0}^{L} \{ \rho A (\ddot{X}\bar{X} + \ddot{Y}\bar{Y}) - EI\theta''\bar{\theta} + EAe\bar{e} - N_x'\bar{X} - N_x\bar{e} + N_x\theta\bar{\theta} - \frac{1}{6}N_x\theta^3\bar{\theta} \\ - N_y'\bar{Y} - N_y\bar{\theta} + \frac{1}{2}N_y\theta^2\bar{\theta} \} dS \, dT. \qquad (C.15)$$

Imposing that $\bar{S} = 0$ for all test functions $\bar{X}, \bar{Y}, \bar{e}$, and $\bar{\theta}$ respectively yields

$$N_x' = \rho A \ddot{X} \qquad (C.16a)$$

$$N_y' = \rho A \ddot{Y} \qquad (C.16b)$$

$$N_x = EAe \qquad (C.16c)$$

$$EI\theta'' = N_x(\theta - \frac{1}{6}\theta^3) - N_y(1 - \frac{1}{2}\theta^2). \qquad (C.16d)$$

Injecting Equation (C.16c) into Equation (C.14a), we obtain Equation (26). Using Equation (C.14) and Equation (27), we rewrite Equation (C.16d) as

$$EI(y''' + y'y''^2 + \frac{1}{2}y'^2y''') = N_xY' - N_y(X' - e - \frac{1}{24}\theta^4) \qquad (C.17)$$

which is Equation (28c) in the small $e$ limit and up to order 3.

### D The order 3 model for a cantilever beam

System (30) is an order 3 approximation of the planar dynamics of an extensible beam with general boundary conditions. In the special case where one end is free and where the beam is considered inextensible ($\eta = 0$), we show here that System (30) reduces to the Crespo da Silva equation (Crespo da Silva and Glynn 1978). We start from Equation (28a) and express $x(s, t)$ with $y(s, t)$

$$x(s, t) = \int_{0}^{s} \left(1 - \frac{1}{2}y'^2\right) ds \Rightarrow \ddot{x}(s, t) = -\int_{0}^{s} \frac{1}{2}(\ddot{y}^2) ds. \qquad (D.1)$$

Next we use the free-end condition at $s = 1$ to integrate Equation (28b) and write

$$n_x(s, t) = -\int_{s}^{1} \ddot{x} \, ds. \qquad (D.2)$$

We then inject Equation (D.1) and Equation (D.2) into Equation (28c) and isolate $n_y(s, t)$:

$$n_y(s, t) = -\frac{y''' + f_3}{1 - \frac{1}{2}y'^2} - \frac{y'}{1 - \frac{1}{2}y'^2} \int_{s}^{1} \ddot{x} \, ds. \qquad (D.3)$$

Differentiating once more with respect to $s$, using Equation (28d), and restricting to order 3 finally yields

$$\ddot{y}(s, t) + (y''' + y'y''^2 + y'^2y''')' + \frac{1}{2}\{y' \int_{1}^{s} \int_{0}^{s} (\ddot{y}^2) ds \, ds\}' = 0 \qquad (D.4)$$

which is Equation (61) in (Thomas et al. 2016).
## E Bifurcation curves plotted with axial load

In Figure E.1, we plot the graphs of Figures 2, 3, 4, and 5 with the axial load $p_E$ on the horizontal axis. This re-plotting sheds light on the shortcomings of the Woinowsky-Krieger model, which fails to predict how the axial load $p_E$ is evolving is the post-buckling regime. One must keep in mind that the frequencies of the first and second modes shown here are not the ones obtained in a load-controlled experiment, since frequencies depend on the axial boundary conditions. For the case considered in this text, the axial distance between the ends of the beam is constant during vibrations, whereas in a load-controlled experiment, this distance varies (i.e. vibrates) since the axial load is prescribed constant. This case is well documented in (Neukirch et al. 2014).

<figure><img src="image_9.png" /><figcaption>**Figure E.1** Upper graphs show post-buckled equilibrium curves. Lower graphs show the frequency curves of the first and second vibration modes around the buckled equilibrium. Kirchhoff model (Ki), Woinowsky-Krieger model (WK), order 3 model of Section 5 (order 3), and 4th-order development of Kirchhoff model (Ki.dev 4). Curve 'Ki.dev 4' not plotted for the second vibration mode. $\eta = 1/4800$.</figcaption></figure>

## F Error of Woinowsky-Krieger model as function of $\eta$

In Figure F.2, we plot curves of constant relative error between the Woinowsky-Krieger and Kirchhoff models, in the plane $(1/\sqrt{\eta}, D/L)$. It is observed that the relative errors $Err_y$ and $Err_\omega$ only weakly depend on the slenderness ratio $\eta$ and grow steadily with the axial displacement $D/L$.

<figure><img src="image_10.png" /><figcaption>**Figure F.2** Comparison between the Woinowsky-Krieger and Kirchhoff models. Relative errors $Err_y = (y_{WK}(1/2) - y_{Ki}(1/2))/y_{Ki}(1/2)$ for the transverse displacement (left) and $Err_\omega = (\omega_{Ki}(1/2) - \omega_{WK}(1/2))/\omega_{Ki}(1/2)$ for the frequency of the first vibration mode (right). Curves of constant error for $Err_y = 2\%, 5\%, 10\%,$ and $20\%$ (left), and $Err_\omega = 1\%, 5\%,$ and $10\%$ (right).</figcaption></figure>
References

Arbocz, J., M. Potier-Ferry, J. Singer, and V. Tvergaard (1985). *Buckling and Post-Buckling*. Vol. 288. Lectures Notes in Physics. Springer-Verlag. [DOI].

Audoly, B. and Y. Pomeau (2010). *Elasticity and Geometry: From hair curls to the non-linear response of shells*. Oxford University Press. ISBN: 9780198506256.

Bažant, Z. J. and L. Cedolin (1991). *Stability of Structures*. Oxford University Press. ISBN: 9789814317030.

Burgreen, D. (1951). Free vibrations of a pin ended column with constant distance between ends. *Journal of Applied Mechanics* 18:135-139.

Cottanceau, E., O. Thomas, P. Véron, M. Alochet, and R. Deligny (2017). A finite element/quaternion/ asymptotic numerical method for the 3D simulation of flexible cables. *Finite Elements in Analysis and Design* 139:14-34. [DOI], [OAI].

Crespo da Silva, M. R. M. and C. C. Glynn (1978). Nonlinear flexural-flexural-torsional dynamics of inextensional beams. I. Equations of motion. *Journal of Structural Mechanics* 6(4):437-448. [DOI].

Dill, E. H. (1992). Kirchhoff's theory of rods. *Archive for History of Exact Sciences* 44(1):1-23. [DOI].

Domokos, G. (1994). Global description of elastic bars. *ZAMM - Journal of Applied Mathematics and Mechanics* 74(4):T289-T291. [DOI].

Eisley, J. G. (1964). Nonlinear vibration of beams and rectangular plates. *Zeitschrift für Angewandte Mathematik und Physik (ZAMP)* 15(2):167-175. [DOI], [HDL].

Eringen, A. C. (1952). On the non-linear vibration of elastic bars. *Quarterly of Applied Mathematics* 9(4):361-369. [DOI].

Gao, D. Y. (2000). Finite deformation beam models and triality theory in dynamical post-buckling analysis. *International Journal of Non-Linear Mechanics* 35:103-131. [DOI].

George, E. P. B. (1976). Science and Statistics. *Journal of the American Statistical Association* 71(356):791-799. [DOI].

Goto, Y., T. Yoshimisu, and M. Obata (1990). Elliptic integral solutions of plane elastica with axial and shear deformations. *International Journal of Solids and Structures* 26(4):375-390. [DOI].

von Kármán, T. (1907). Festigkeitsprobleme im Maschinenbau. *Encyklopädie der Mathematischen Wissenschaften*. Ed. by F. Klein and C. Müller. Mechanik. Teubner Verlag, pp 311-385. [DOI].

Kirchhoff, G. (1876). *Vorlesungen über Mathematische Physik: Mechanik*. Teubner. [ARK].

Lacarbonara, W. (2013). *Nonlinear Structural Mechanics: Theory, Dynamical Phenomena and Modeling*. Springer. [DOI].

Lou, C. and D. Sikarskie (1975). Nonlinear vibration of beams using a form-function approximation. *Journal of Applied Mechanics* 42(1):209-214. [DOI].

Luongo, A. and D. Zulli (2013). *Mathematical Models of Beams and Cables*. Wiley. [DOI].

Mamou-Mani, A., J. Frelat, and C. Besnainou (2009). Prestressed soundboards: analytical approach using simple systems including geometric nonlinearity. *Acta Acustica united with Acustica* 95:915-928. [DOI].

McDonald, P. H. (1991). Nonlinear motion of a beam. *Nonlinear Dynamics* 2:187-198. [DOI].

Meier, C., A. Popp, and W. A. Wall (2019). Geometrically exact finite element formulations for slender beams: Kirchhoff-Love theory versus Simo-Reissner theory. *Archives of Computational Methods in Engineering* 26(1):163-243. [DOI], [ARXIV].

Mettler, E. (1951). Zum Problem der Stabilität erzwungener Schwingungen elastischer Körper. *ZAMM - Journal of Applied Mathematics and Mechanics* 31(8-9):263-264. [DOI].

Naghdi, P. M. (1963). A new derivation of the general equations of elastic shells. *International Journal of Engineering Science* 1(4):509-522. [DOI].

Naghdi, P. M. (1964). Further results in the derivation of the general equations of elastic shells. *International Journal of Engineering Science* 2(3):269-273. [DOI].

Nayfeh, A. H. and S. Emam (2008). Exact solution and stability of postbuckling configurations of beams. *Nonlinear Dynamics* 54(4):395-408. [DOI].

Nayfeh, A. H. and D. T. Mook (1995). *Nonlinear Oscillations*. Wiley. ISBN: 9783527617593.

Nayfeh, A., W. Kreider, and T. Anderson (1995). Investigation of natural frequencies and mode shapes of buckled beams. *AIAA Journal* 33(6):1121-1126. [DOI].
Nayfeh, A. H. and P. F. Pai (2004). *Linear and Nonlinear Structural Mechanics*. Wiley. [DOI].

Neukirch, S., J. Frelat, A. Goriely, and C. Maurini (2012). Vibrations of post-buckled rods: The singular inextensible limit. *Journal of Sound and Vibration* 331(3):704–720. [DOI], [OAI].

Neukirch, S., A. Goriely, and O. Thomas (2014). Singular inextensible limit in the vibrations of post-buckled rods: Analytical derivation and role of boundary conditions. *Journal of Sound and Vibration* 333(3):962–970. [DOI], [OAI].

Neukirch, S., M. Yavari, N. Challamel, and O. Thomas (2021). Comparison of the Von Kármán and Kirchhoff models for the post-buckling and vibrations of elastic beams. Preprint and supplementary files. [OAI].

Pflüger, A. (1964). *Stabilitätsprobleme der Elastostatik*. 2nd ed. Springer-Verlag. [DOI].

Ray, J. D. and C. W. Bert (1969). Nonlinear vibrations of a beam with pinned ends. *Journal of Engineering for Industry* 91(4):997–1004. [DOI].

Reissner, E. (1972). On one-dimensional finite-strain beam theory: The plane problem. *Zeitschrift für Angewandte Mathematik und Physik* (ZAMP) 23:795–804. [DOI].

Thomas, O., A. Sénéchal, and J.-F. Deü (2016). Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams. *Nonlinear Dynamics* 86(2):1293–1318. [DOI], [OAI].

Timoshenko, S. P. (1921). On the correction for shear of the differential equation for transverse vibrations of prismatic bars. *Philosophical Magazine* 41:744–746. [DOI].

Woinowsky-Krieger, S. (1950). The effect of an axial force on the vibration of hinged bars. *Journal of Applied Mechanics* 17:35–36. [OAI].

Younis, M. I. and A. H. Nayfeh (2003). A study of the nonlinear response of a resonant microbeam to an electric actuation. *Nonlinear Dynamics* 31:91–117. [DOI].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

<figure><img src="image_11.png" /></figure>

**Authors' contributions** SN: Analysis, Computations, Writing. MY: Computations. NC: Analysis, Writing. OT: Analysis, Writing. All authors discussed the results, read and approved the final manuscript.

**Supplementary Material** The data used to plot the figures of the article can be found at the permalink (*See related files* tab): hal-02957425v2.

**Acknowledgements** SN would like to thank Corrado Maurini and Adrien Mamou-Mani for igniting his interest in Equation (8).

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.