# Fourier analysis of membrane locking and unlocking

René R. Hiemstra<EMAIL_ADDRESS>Federico Fuentes
<EMAIL_ADDRESS>Dominik Schillinger  dominik.schillinger@tu-
darmstadt.de Institute for Mechanics, Computational Mechanics Group, Technical
University of Darmstadt, Germany  Institute for Mathematical and Computational
Engineering (IMC), School of Engineering and Faculty of Mathematics,
Pontificia Universidad Católica de Chile, Chile
Dedicated to Thomas J.R. Hughes on the occasion of his 80th birthday.

###### Abstract

Membrane locking in finite element approximations of thin beams and shells has
remained an unresolved topic despite four decades of research. In this
article, we utilize Fourier analysis of the complete spectrum of natural
vibrations and propose a criterion to identify and evaluate the severity of
membrane locking. To demonstrate our approach, we utilize standard and mixed
Galerkin formulations applied to a circular Euler-Bernoulli ring discretized
using uniform, periodic B-splines. By analytically computing the discrete
Fourier operators, we obtain an exact representation of the normalized error
across the entire spectrum of eigenvalues. Our investigation addresses key
questions related to membrane locking, including mode susceptibility, the
influence of polynomial order, and the impact of shell/beam thickness and
radius of curvature. Furthermore, we compare the effectiveness of mixed and
standard Galerkin methods in mitigating locking. By providing insights into
the parameters affecting locking and introducing a criterion to evaluate its
severity, this research contributes to the development of improved numerical
methods for thin beams and shells.

###### keywords:

membrane locking , Fourier analysis , natural vibrations , Euler-Bernoulli
beam , Isogeometric analysis

††journal: Computer Methods in Applied Mechanics and Engineering (CMAME)

##  1 Introduction

It was late 2009 when I first heard about isogeometric analysis. As a fresh
M.Sc. student with a keen interest in research, particularly in finite element
analysis, I delved into the first few papers by Tom and his former students.
It became evident that this technology was a game-changer. I promptly acquired
a copy of the IGA book released earlier that year and absorbed its contents in
a matter of weeks. The prospect of working with such celebrated researchers
ignited my imagination.

Our paths converged two years later at the World Congress on Computational
Mechanics in Vienna. Tom, a prominent figure in the field, delivered the
opening talk to a crowd of roughly 4000 people. After his presentation, I
seized the opportunity to approach him, and to my delight, he generously
engaged in a conversation. It was during this interaction that I discovered
Tom’s genuine interest in my research topic. Two days later, as I presented my
work, I was pleasantly surprised to find Tom among the audience. Despite my
initial nerves, I managed to deliver a solid presentation. This pivotal moment
led to a remarkable turn of events—Tom extended an invitation to visit his
research group in Austin, an invitation that effectively offered me a Ph.D.
opportunity.

I had the incredible fortune of working alongside Tom, and it was a truly
remarkable experience. Our initial encounter turned out to be just the
beginning of a series of memorable interactions. Tom possesses a unique
quality—he always makes time for his students and others, engaging with them
not only on scientific matters but also on personal and amusing topics. He is
quite the collector of anecdotes, regaling us with tales from the “good old
days” of finite elements in the late sixties and seventies.

Allow me to share two of Tom’s captivating stories, a testament to his sense
of humor and ability to find joy in unexpected situations. Picture a young Tom
in an elevator with the late Oleg Zienkiewicz and a youthful Robert Taylor. In
this particular episode, Robert earnestly tries to explain to Oleg why his
penalty formulation is fundamentally flawed, emphasizing that “the solution is
zero to all problems”. However, Zienkiewicz interrupts him with a playful
remark, urging Robert not to bother him with details and to just implement the
idea. In another story from a conference in western Germany, Tom recounts how
a group of young scientists, including Ted Belytschko, Jurgen Bathe, and Wing
Kam Liu, found themselves scribbling equations in a café. Unwittingly, their
activity raised suspicion, and they were mistaken for terrorists, leading to
an encounter with the police. It’s these glimpses into Tom’s past that reveal
not only his expertise but also his wit and capacity to navigate scientific
collaborations with a light-hearted approach.

Dear Tom, on the momentous occasion of your 80th birthday, I extend my warmest
congratulations to you. As you gracefully embrace this respectable age, it’s
worth noting that you possess a vibrant spirit that defies the passing years.
I dedicate this contribution to you, considering your unwavering interest in
shell finite elements, spectral analysis, and the intriguing subject of
membrane locking. I believe you will thoroughly enjoy reading it and find it
enriching. With warm regards, Rene

Membrane locking in finite element formulations of curved beam and shell
models refers to a manifestation of artificial bending stiffness resulting
from the coupling between the bending and membrane response induced by the
curvature [1, 2]. In an extreme case, a curved beam, which would be expected
to exhibit a bending-dominated response, instead displays an in-plane
membrane-dominated discrete solution. Membrane locking is only one of several
sources of locking, the most well-known being transverse shear locking in
beam, plate and shell elements [1, 3] and volumetric locking due to
incompressibility in solid elements [4, 5]. All sources of locking may
severely affect accuracy and robustness of finite element discretizations. The
development of locking-preventing discretization technology has a history of
more than 40 years, first within classical finite elements and more recently
in isogeometric analysis. Common techniques include, but are not limited to,
reduced and selective integration techniques, strain projection methods, and
certain mixed formulations based on alternative variational principles. An
extensive discussion on several forms of locking in relation to finite element
design can be found in [1, Chapter 5.4] and [6, Chapter 9.4, and Chapter 10].
An analysis of locking effects from the mathematical point of view can be
found in [7, Chapter 6]. We also refer to [8] for a concise overview of the
relevant literature until 2013 and to [9, Table 1] for an overview of more
recent approaches, in particular in the context of isogeometric analysis.

From a numerical analysis viewpoint, locking has been described as a loss of
robustness (or uniform approximation) as a problem dependent parameter $t$
approaches a limiting value $t_{0}$ [4, 10]. For example, Poisson’s ratio
approaching the incompressibility limit of one half and plate or shell
thickness approaching zero. Most a priori error estimates for elliptic
boundary value problems yield optimal asymptotic convergence rates when the
parameter $t\neq t_{0}$ is fixed. This fact does not necessarily aid the study
of locking phenomena, because these often manifest at coarse discretizations
and vanish with mesh refinement. A series of early works [3, 4, 10, 11] are
devoted to the mathematical analysis of finite element methods in the context
of parameter dependency and locking. It has been shown how solid elements lose
uniform convergence in the incompressibility limit [4] and how uniform
approximation is lost due to shear locking as the beam or plate thickness
approaches zero [3, 10]. However, membrane locking as a standalone phenomenon
has not received the same level of mathematical scrutiny as the two other main
sources of locking. The limited attention to membrane locking can be
attributed, in part, to its physical rather than mathematical definition. It
refers to the propensity for a higher presence of membrane stresses or
displacements, as opposed to bending stresses and displacements, in comparison
to the exact solution. Consequently, its distinct characteristics and
mathematical analysis have not been thoroughly explored in the literature.

Benchmark problems designed to evaluate the performance and correctness of
beam and shell elements [12, 13] have proven indispensable in the development
of robust shell elements, but alone are not adequate to assess membrane
locking and unlocking. These benchmarks cannot possibly cover all conditions
and can therefore not positively rule out membrane locking. In addition, due
to the complexity in loading and boundary conditions, these tests cannot
isolate and separate the effect of membrane locking from other parasitic
behavior due to, e.g., boundary layers.

Recently, we presented a novel approach to enhance the understanding of
membrane locking and unlocking in thin beams and shells through spectral
analysis [9]. The eigenvalue decomposition is a valuable tool because it
offers a comprehensive understanding of a linear operator111Consider a well
posed linear elliptic boundary value problem: find $u\in V$ that satisfies
$a(u,v)=l(v)$ for any $v\in V$. Suppose $u_{n}\in V$ are given and satisfy
$a(u_{m},u_{n})=\lambda_{n}\delta_{mn}$, then the solution allows the
decomposition: $u=\sum\nolimits_{n}\alpha_{n}\cdot u_{n}$, where
$\alpha_{n}=l(u_{n})/\lambda_{n}$ whenever $\lambda_{n}\neq 0$ and free
otherwise. The functions $u_{n}\in V$ are called the eigenfunctions or
eigenmodes of $a$ and $\lambda_{n}=a(u_{n},u_{n})$ its eigenvalues.. By
analyzing the eigenvalues and eigenfunctions, we can gain valuable insights
into the properties and behavior of the operator, without being confined by
the specific details of a boundary value problem, such as boundary conditions
and the right-hand side. Eigenvalue and modal analysis provides a holistic
view of the discretization error associated with a method. Specifically, the
discretization error, such as the one occurring in the energy norm of an
elliptic boundary value problem, can be expressed solely in terms of the
errors made in eigenvalues and modes [14, 15].

In [9] we examined five different techniques to prevent or mitigate membrane
locking in the context of a curved Euler-Bernoulli circular ring. The
simplicity and periodic nature of the model problem made it possible to
isolate the effect of membrane locking from shear locking and spurious
boundary phenomena. Locking was measured, at fixed normalized mode number, as
the discrepancy between the spectrum error and the spectrum error in the limit
of refinement. The severity of locking and the performance of locking
mitigation or removal could thereby be measured in a quantitative manner,
which enabled an objective comparison. We note that the use of spectral or
Fourier analysis in relation to locking is not a completely new idea.
Inspection of the low eigenvalues to investigate the presence and effect of
locking has been considered in the design of efficient shell elements [16, 17]
or multigrid methods [18]. Uncoupled discrete Fourier operators and
corresponding characteristic equations have been used to identify spurious
effects and locking in early work [19].

In this article, we provide a comprehensive analysis of the discrete error
through Fourier (spectral) analysis, revealing important insights into locking
phenomena. We analytically show that the numerical solution of the isolated
membrane or bending response exhibits normalized spectra that remain invariant
regardless of mesh refinement. However, our model displays a coupled response
for which the normalized spectra no longer possess this invariance, except in
the limit of mesh refinement. This observation motivates the definition of a
concrete criterion that characterizes the presence and severity of membrane
locking as a function of the normalized spectrum. To illustrate the issue, we
investigate a standard Galerkin formulation employing quadratic, cubic, and
quartic spline function spaces, and demonstrate its susceptibility to membrane
locking. In contrast, a mixed formulation based on the Hellinger-Reissner
variational principle proves to be free from locking for all polynomial
degrees and offers significantly improved absolute accuracy across all levels
of mesh refinement.

Our solution approach revolves around an analytical technique to solve the
discrete eigenvalue problem, leveraging properties of circulant matrices that
arise during discretization with periodic uniform spline basis functions. We
establish the existence of a discrete analogue to the Fourier basis,
effectively decoupling the $2N\times 2N$ system of equations into $N$ systems
of $2\times 2$ equations. The solution approach, depicted in Figure 1, closely
resembles the exact solution obtained using a Fourier basis. This analytical
approach enables us to compute exact expressions of the normalized error
across the complete eigenvalue spectrum, incorporating parameters such as
element size, thickness, and radius of curvature. Our analysis reveals how
these relevant physical and discretization parameters contribute to membrane
locking within the standard formulation. Notably, while our focus lies on
quadratic discretizations of the standard formulation, our techniques apply to
general polynomial degrees and can be adapted to the mixed formulation as
well.

Although the dynamics of a circular Euler-Bernoulli ring are well-documented,
as detailed in [20, Chapter 5, Section 5.3, page 82], a comparison between the
analytical and discrete spectra, for both membrane-dominated and bending-
dominated eigenmodes, offers fresh insights into locking phenomena.
Significantly, our findings establish that membrane locking solely affects
eigenvalues corresponding to modes with a bending-dominated response,
effectively overdampening the contributions of these eigenmodes from the
discrete solution, while leaving the eigenvalues associated to a membrane-
dominated response mostly unaffected. We also find that the phenomenon of
membrane locking depends nontrivially on an interplay between a normalized
thickness and the discretization size, and in a particular regime, is
mathematically analogous to the phenomenon of shear locking (despite the fact
that there are no shear strains in our model).

Figure 1: The solution approach of the discrete eigenvalue problem mirrors the
solution approach to the exact solution.

The paper is organized as follows: in Section 2, we present the standard and
mixed formulation of the eigenvalue model problem and investigate the behavior
of the analytical solution for our model problem of a circular ring. In
Section 3, we present the standard and mixed Galerkin discretization using
periodic uniform splines, and discuss an analytical solution technique based
on the Discrete Fourier Transform (DFT). Section 4 introduces a criterion to
identify and measure the severity of membrane locking in thin beams and
shells. It is based on the discrepancy of the normalized spectrum of
eigenvalues in relation to the spectrum in the limit of mesh refinement. These
notions are applied in Section 5 to the numerical results stemming from the
computation of the spectra of the circular Euler-Bernoulli ring. Here we
juxtapose the standard formulation with the mixed formulation and discuss
performance with respect to accuracy and robustness. In Section 6, we perform
a mathematical analysis of the spectrum in the standard formulation and derive
exact expressions for the relative error and its dependence on different
parameters. In Section 7 we further discuss some observations on the
phenomenon of membrane locking in the context of our study. Finally, in
Section 8 we draw conclusions, and propose potential future work.

##  2 The continuous eigenvalue problem

In the scope of this paper, we study the phenomenon of membrane locking in the
special case of an Euler-Bernoulli ring. We highlight that this problem rules
out the presence of shear strains, so is ideal to isolate membrane-related
phenomena. In this section, we introduce the continuous mixed and standard
eigenvalue problem. We then discuss properties of the analytical solution.

Figure 2: Closed circular ring modeled as a curved Euler-Bernoulli beam.

###  2.1 The mixed formulation

Let  denote the circular ring of radius $R$ displayed in Figure 2,
parameterized by arc-length parameter $s\in[0,2\pi R)$. The field variables
are the circumferential displacement component $u$, the transverse
displacement component $w$, the membrane strain resultant $\varepsilon$, and
the bending strain resultant $\chi$. Displacements are considered as elements
from the space $V=H^{2}(\leavevmode\hbox to9.34pt{\vbox
to9.34pt{\pgfpicture\makeatletter\hbox{\hskip 4.66791pt\lower-2.66791pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}})\times
H^{2}(\leavevmode\hbox to9.34pt{\vbox
to9.34pt{\pgfpicture\makeatletter\hbox{\hskip 4.66791pt\lower-2.66791pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}})$ and strains as
elements in $S=L^{2}(\leavevmode\hbox to9.34pt{\vbox
to9.34pt{\pgfpicture\makeatletter\hbox{\hskip 4.66791pt\lower-2.66791pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}})\times
L^{2}(\leavevmode\hbox to9.34pt{\vbox
to9.34pt{\pgfpicture\makeatletter\hbox{\hskip 4.66791pt\lower-2.66791pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}})$.

We consider the following eigenvalue problem: find $\boldsymbol{u}=(u,w)\in
V$, $\boldsymbol{\varepsilon}=(\varepsilon,\chi)\in S$ and
$\lambda\in\mathbb{R}^{+}_{0}$ such that

$\displaystyle-c(\boldsymbol{\varepsilon},\,\boldsymbol{\delta\varepsilon})+c(B(\boldsymbol{u}),\,\boldsymbol{\delta\varepsilon})$
$\displaystyle=0$ $\displaystyle\forall\boldsymbol{\delta\varepsilon}$
$\displaystyle=(\delta\varepsilon,\delta\chi)\in S$ (1a) $\displaystyle
c(\boldsymbol{\varepsilon},\,B(\boldsymbol{\delta u}))$ $\displaystyle=\lambda
b(\boldsymbol{u},\,\boldsymbol{\delta u})$
$\displaystyle\forall\boldsymbol{\delta u}$ $\displaystyle=(\delta u,\delta
w)\in V.$ (1b)

where

$\displaystyle c(\boldsymbol{\varepsilon},\,\boldsymbol{\delta\varepsilon})$
$\displaystyle:=\int_{0}^{2\pi
R}\left(EA\,\varepsilon\,\delta\varepsilon+EI\,\chi\,\delta\chi\right)\,ds$
$\displaystyle\boldsymbol{\varepsilon},\boldsymbol{\delta\varepsilon}\in S$
(1c) $\displaystyle b(\boldsymbol{u},\,\boldsymbol{\delta u})$
$\displaystyle:=\int_{0}^{2\pi R}\rho A\left(u\,\delta u+w\,\delta
w\right)\,ds$ $\displaystyle\boldsymbol{u},\boldsymbol{\delta u}\in V.$ (1d)

Here $E,A,I$ and $\rho$ denote, respectively, the Young’s modulus, the cross
sectional area and moment of inertia, and the density of the ring.
Furthermore, $B\;:\;V\mapsto S$ is defined via the kinematic relations

$\displaystyle
B(\boldsymbol{u}):=\begin{cases}\varepsilon(\boldsymbol{u})&=\frac{\partial
u}{\partial s}+\frac{1}{R}w\\\
\chi(\boldsymbol{u})&=-\frac{\partial^{2}w}{\partial
s^{2}}+\frac{1}{R}\frac{\partial u}{\partial s}\end{cases}$ (1e)

Membrane locking occurs due to the disparate nature of work due to membrane
stresses (first term in (1c)) and work due to bending (second term in (1c)).
To study the phenomenon we are interested only in the relative magnitude of
these terms. In the following we transform the problem to polar coordinates
$(r,\theta)$, factor out the density and Young’s modulus, and assume a
rectangular cross section with constant thickness, that is,

$\displaystyle E/\rho=1,$ $\displaystyle t=\text{const},$ $\displaystyle
I=\frac{At^{2}}{12},$ $\displaystyle
B(\boldsymbol{u}):=\begin{cases}\varepsilon(\boldsymbol{u})&=\frac{1}{R}\frac{\partial
u}{\partial\theta}+\frac{1}{R}w\\\
\chi(\boldsymbol{u})&=-\frac{1}{R^{2}}\frac{\partial^{2}w}{\partial\theta^{2}}+\frac{1}{R^{2}}\frac{\partial
u}{\partial\theta}\end{cases}$ (1f)

The bilinear forms, $c\,:\,S\times S\mapsto\mathbb{R}$ and $b\,:\,V\times
V\mapsto\mathbb{R}$, then simplify to

$\displaystyle c(\boldsymbol{\varepsilon},\,\boldsymbol{\delta\varepsilon})$
$\displaystyle:=\int_{0}^{2\pi}\left(\varepsilon\,\delta\varepsilon+\tfrac{t^{2}}{12}\,\chi\,\delta\chi\right)\,R\,d\theta$
$\displaystyle\boldsymbol{\varepsilon},\boldsymbol{\delta\varepsilon}\in S$
(1g) $\displaystyle b(\boldsymbol{u},\,\boldsymbol{\delta u})$
$\displaystyle:=\int_{0}^{2\pi}\left(u\,\delta u+w\,\delta
w\right)\,R\,d\theta$ $\displaystyle\boldsymbol{u},\boldsymbol{\delta u}\in
V.$ (1h)

###  2.2 The standard formulation

Equivalently, we may consider the standard formulation obtained from the mixed
formulation by substituting $\boldsymbol{\varepsilon}=B(\boldsymbol{u})$ into
equation (1b). We seek $\boldsymbol{u}=(u,w)\in V$ and
$\lambda\in\mathbb{R}^{+}_{0}$ such that

$\displaystyle a(\boldsymbol{u},\,\boldsymbol{\delta u})=\lambda
b(\boldsymbol{u},\,\boldsymbol{\delta u})\qquad\forall\boldsymbol{\delta
u}=(\delta u,\delta w)\in V.$ (2a) Here $a\,:\,V\times V\mapsto\mathbb{R}$
denotes the semi-definite symmetric bilinear form defined as $\displaystyle
a(\boldsymbol{u},\,\boldsymbol{\delta
u}):=c(B(\boldsymbol{u}),\,B(\boldsymbol{\delta
u})),\qquad\boldsymbol{u},\boldsymbol{\delta u}\in V.$ (2b)

The standard and mixed formulations are equivalent, that is, they have the
same eigenvalues and eigenfunctions.

###  2.3 Analytical computation of exact eigenvalues and eigenfunctions

Let $\phi\in[0,2\pi)$ denote an arbitrary phase shift. An analytical solution
to (1) may be obtained by employing the following informed guess for
$\boldsymbol{u}_{n}=(u_{n},\,w_{n})$,

$\displaystyle u_{n}(\theta)$
$\displaystyle=U_{n}\sin{\left(n(\theta-\phi)\right)},$ (3a) $\displaystyle
w_{n}(\theta)$ $\displaystyle=W_{n}\cos{\left(n(\theta-\phi)\right)}.$ (3b)

The eigenvalues and amplitudes of the eigenfunctions satisfy, for each $n$, a
two by two matrix eigenvalue problem (see C or [20, Chapter 5, Section 5.3,
page 82])

$\displaystyle\begin{pmatrix}\tilde{A}_{n}-\lambda_{n}&\tilde{B}_{n}\\\
\tilde{B}_{n}^{*}&\tilde{D}_{n}-\lambda_{n}\end{pmatrix}\begin{pmatrix}U_{n}\\\
W_{n}\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix}.$ (4)

The values $\tilde{A}_{n}$, $\tilde{B}_{n}$ and $\tilde{D}_{n}$ are listed in
C in terms of the beam thickness $t$ and the radius $R$. For each $n$, there
are two eigenfunctions and corresponding eigenvalues. The eigenvalues are
computed as the roots of the characteristic equation

$\displaystyle\lambda_{n+}$ $\displaystyle=K_{n}\left(1+L_{n}\right)$ (5a)
$\displaystyle\lambda_{n-}$ $\displaystyle=K_{n}\left(1-L_{n}\right)$ (5b)

where $K_{n}$ and $L_{n}$ are constants that depend on $n$ , the normalized
thickness $\bar{t}=t/R$ and radius $R$, and are given by

$\displaystyle
K_{n}=\frac{\left(n^{2}+1\right)\,\left(n^{2}\,\bar{t}^{2}+12\right)}{24R^{2}},\qquad
M_{n}=\frac{n^{2}\,\bar{t}^{2}\,{\left(n^{2}-1\right)}^{2}}{12R^{4}},\qquad
L_{n}=\sqrt{1-\frac{M_{n}}{K_{n}^{2}}}.$ (6)

The eigenfunctions are defined up to a constant. However, the ratio of the
amplitudes of the eigenfunctions, $U_{n}/W_{n}$, satisfy

$\displaystyle\frac{U_{n+}}{W_{n+}}=\frac{\tilde{B}_{n}}{\lambda_{n+}-\tilde{A}_{n}}$
(7a)
$\displaystyle\frac{U_{n-}}{W_{n-}}=\frac{\tilde{B}_{n}}{\lambda_{n-}-\tilde{A}_{n}}$
(7b)

When $U_{n}/W_{n}<1$ the transverse displacement is larger than the
circumferential displacement, thus the bending response dominates. When
$U_{n}/W_{n}>1$, the circumferential displacements are larger than the
transverse displacements and, consequently, the membrane response dominates.
In what follows we will refer to these two types of modes as bending-dominated
and membrane-dominated. Figure 4 and 4 depict several examples of the two
different types of modes.

(a) $n=1$ (b) $n=2$

(c) $n=3$ (d) $n=4$

(e) $n=5$ (f) $n=6$

Figure 3: The first eight analytical bending-dominated eigenfunctions.

(a) $n=0$ (b) $n=1$

(c) $n=2$ (d) $n=3$

(e) $n=4$ (f) $n=5$

Figure 4: The first eight analytical membrane-dominated eigenfunctions.

Figure 5 provides a more detailed view of the behavior of the ratio of the
amplitudes, $U_{n}/W_{n}$, at two different choices of the normalized
thickness. Modes $\boldsymbol{u}_{n-}$, with a ratio of amplitudes
$U_{n-}/W_{n-}$ and corresponding eigenvalue $\lambda_{n-}$, have dominating
transverse displacements for $n<\hat{n}$, where $\hat{n}$ is a value that
depends on the normalized thickness $\bar{t}$. In contrast, modes
$\boldsymbol{u}_{n+}$, with a ratio of amplitudes $U_{n+}/W_{n+}$ and
corresponding eigenvalue $\lambda_{n+}$, have dominating circumferential
displacements for $n<\hat{n}$. When $n>\hat{n}$ this behavior reverses.

(a) $\bar{t}=t/R=0.1$ (b) $\bar{t}=t/R=0.01$

Figure 5: For $n<\hat{n}$ the modes corresponding to $\lambda_{n+}$ are
membrane-dominated since $U_{n+}/W_{n+}>1$ and the modes corresponding to
$\lambda_{n-}$ are bending-dominated since $U_{n-}/W_{n-}<1$. For $n>\hat{n}$
this behavior is reversed. This is shown for two different normalized
thicknesses $\bar{t}$, since $\hat{n}$ depends on $\bar{t}$.

(a) $\bar{t}=t/R=0.1$ (b) $\bar{t}=t/R=0.01$

Figure 6: The analytical spectrum plotted on a log-log scale. Typically
$\lambda_{n+}>>\lambda_{n-}$ except near $n=\hat{n}$, which is the location of
the local minimum of the function $K_{n}\cdot L_{n}$, denoted in blue.

The corresponding branches of eigenvalues, $\lambda_{n+}$ and $\lambda_{n-}$,
also behave differently, as may be observed in Figure 6, where they are
plotted as a continuous (rather than discrete) function of $n$. For reasonable
values of the material and cross sectional parameters it follows that
$\lambda_{n+}>>\lambda_{n-}$, except near $n=\hat{n}$, which is the location
of the local minimum of the function $K_{n}\cdot L_{n}$, shown in blue in
Figure 6. At this point the argument inside the square root,
$L^{2}_{n}=1-M_{n}/K_{n}^{2}$, attains a global minimum, which itself
approaches zero when viewed as a function of $\bar{t}$ with
$\bar{t}\rightarrow 0$, see Figure 7.

Figure 7: The function $L_{n}^{2}(n,\bar{t})=1-M_{n}/K^{2}_{n}$ is a smooth
positive function that is less than or equal to $1$. For $n\geq 1$, the
function attains a single minimum at mode number $\hat{n}$, which depends on
the value of $\bar{t}=t/R$.

More precisely, $L_{n}^{2}$ depends on $n$ and on the normalized thickness
$\bar{t}$, but not on the radius,

$\displaystyle
L^{2}_{n}=1-\frac{48n^{2}\,\bar{t}^{2}\,(n^{2}-1)^{2}}{(n^{2}+1)^{2}(n^{2}\,\bar{t}^{2}+12)^{2}}$
(8)

Thus, the mode number at which the minimum of $L_{n}^{2}$ is attained, denoted
by $\hat{n}$, can be determined as a function of the normalized thickness
$\bar{t}=t/R$. The minimum coincides with the maximum of the fraction in (8).
Let $f$ denote its numerator and $g$ the denominator. At an optimum we require
that $f^{\prime}(n)g(n)-f(n)g^{\prime}(n)=0$, or, since $g(n)>0$ for any
$n\geq 0$, $f^{\prime}(n)-f(n)g^{\prime}(n)/g(n)=0$, meaning that we must seek
the zeros of the function

$\displaystyle z(n)$ $\displaystyle=f^{\prime}(n)-f(n)g^{\prime}(n)/g(n)$
$\displaystyle=\frac{96\,n\,\bar{t}^{2}\,\left(n^{2}-1\right)\,\left(-\bar{t}^{2}\,n^{6}+(4\bar{t}^{2}+12)n^{4}+(\bar{t}^{2}+48)n^{2}-12\right)}{(n^{2}+1)(n^{2}\,\bar{t}^{2}+12)}\,.$

By inspection, the roots of $z(n)$ are $n=0$, $n=\pm 1$, and the roots of the
sixth order polynomial in the numerator. Focusing on the latter we substitute
$x=n^{2}$, and obtain a cubic polynomial,
$h(x)=-\bar{t}^{2}\,x^{3}+(4\,\bar{t}^{2}+12)x^{2}+(\bar{t}^{2}+48)x-12$,
whose only root above unity, $n^{2}=x>1$, is given by Cardano’s formula, and
corresponds precisely to $\hat{n}$ (i.e., the minimum of $L_{n}$ for $n\geq
1$). We do not show the explicit expression for $\hat{n}$, as it is very
lengthy, but the main point is that it is an algebraic expression depending
solely on the normalized thickness $\bar{t}$. Figure 8 plots the value of
$\hat{n}$ as a function of $\bar{t}$ in a log-scale, and we note the function
is very well approximated (with negligible mean squared error in a broad range
of $\bar{t}$) by a linear function in log-space given by

$\displaystyle\log_{10}\big{(}\hat{n}(\bar{t})\big{)}\approx-0.9985\log_{10}(\bar{t})+0.5421$
(9)

Figure 8: The position of $\hat{n}\geq 1$, where $L^{2}_{n}(n,\bar{t})$
achieves a minimum, as a function of the normalized thickness $\bar{t}=t/R$.

##  3 Discretization by periodic uniform splines

In this section we introduce the discrete function spaces used to discretize
the standard and mixed eigenvalue problem presented in the previous section.
We consider as trial space a maximally smooth periodic space of splines
defined on a uniform partition of the circle. In this special case there
exists a discrete Fourier basis that inherits many of the properties of
analytical Fourier modes (see Figure 1), allowing us analytically compute
expressions for discrete eigenvalues and the ratio of amplitudes of the modes.

###  3.1 Periodic uniform splines

This section introduces Cardinal and uniform splines and discusses how they
are used to develop a discrete trial and test space of splines.

#### Cardinal B-spline function

A convenient basis in which to represent polynomial splines is given by
B-splines. An important special case is the so called Cardinal B-spline.

######  Definition 3.1 (Cardinal B-spline).

Let $\boldsymbol{1}_{[0,1)}$ denote the characteristic function on $[0,1)$.
The Cardinal B-spline of polynomial degree $p$ is defined via the convolution

$\displaystyle\phi_{p}:=\phi_{p-1}\ast\boldsymbol{1}_{[0,1)}$ (10)

A graphical interpretation of Definition 3.1 is depicted in Figure 9.
Convolution increases the polynomial degree as well as the smoothness. Several
useful properties follow directly. Here, $\phi_{p}(x)$ is a degree $p$
polynomial on every element $(k,k+1)$, is $C^{p-1}$ smooth, positive inside
the interval $(0,p+1)$ and zero elsewhere.

$\phi_{0}$$\phi_{1}$$\phi_{2}$$\phi_{3}$ Figure 9: Convolution by the
characteristic function increases smoothness and polynomial degree.

#### Uniform splines

Consider a partitioning of the real line with mesh-size $h$

$\displaystyle h\mathbb{Z}:=\ldots-2h,\;-h,\;0,\;h,\;2h,\;\ldots$ (11)

A space of uniform B-splines on $h\mathbb{Z}$ is obtained by translated and
scaled copies of $\phi_{p}$.

######  Definition 3.2 (Uniform B-splines).

For $h>0$

$\displaystyle B_{i,p}(x)=\phi_{p}(x/h-i),\qquad i\in\mathbb{Z}$ (12)

is the $i$-th B-spline defined on the grid $h\mathbb{Z}$. Linear combinations
$\sum_{i\in\mathbb{Z}}\alpha_{i}\,B_{i,p}(x)$ are called uniform splines of
polynomial degree $p$ and mesh-size $h$.

#### Periodic uniform splines on the unit ring

Periodic spaces of splines of dimension $N$ can be constructed by taking $N+p$
sequential B-splines and applying suitable end-conditions to the last $p$
B-splines. In Definition 3.3 this is done for a smooth uniform space of
splines defined on the unit circle. Figure 10 shows a graphical
interpretation.

######  Definition 3.3 (Periodic uniform spline space ).

Let $h=2\pi/N$ denote the mesh size and $\theta\in[0,2\pi)$ denote the angular
coordinate. The periodic $N$-dimensional space of splines on the unit-circle
is defined as

$\displaystyle\mathbb{S}^{p}_{N}(0,2\pi)=\left\\{s\;:\;[0,2\pi)\mapsto\mathbb{R}\;:\;s(\theta)=\sum_{i=-p}^{N-1}\alpha_{i}\,B_{i,p}(\theta),\;\text{with
}\alpha_{-p}=\alpha_{N-p},\ldots,\alpha_{-1}=\alpha_{N-1}\right\\}$ (13)

Uniform, periodic B-splines inherit local support, positivity, and are
globally $C^{p-1}$ smooth. The periodic and translation invariant properties
of the basis lead to circulant system matrices when used in weak forms. Such
matrices are easily diagonalized via the Discrete Fourier Transform (DFT), see
A. Since we are dealing with piecewise polynomials, differentiation and
integration can be performed analytically using the standard spline formulae,
see [21].

Figure 10: The periodic space $\mathbb{S}^{2}_{8}(0,2\pi)$, consisting of $8$
quadratic uniform B-splines on the circle

###  3.2 The mixed formulation

We consider discrete function spaces of the form

$\displaystyle V^{h}=\mathbb{S}^{p}_{N}(\leavevmode\hbox to9.34pt{\vbox
to9.34pt{\pgfpicture\makeatletter\hbox{\hskip 4.66791pt\lower-2.66791pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}})\times\mathbb{S}^{p}_{N}(\leavevmode\hbox
to9.34pt{\vbox to9.34pt{\pgfpicture\makeatletter\hbox{\hskip
4.66791pt\lower-2.66791pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}),$ $\displaystyle
S^{h}=\mathbb{S}^{p-1}_{N}(\leavevmode\hbox to9.34pt{\vbox
to9.34pt{\pgfpicture\makeatletter\hbox{\hskip 4.66791pt\lower-2.66791pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}})\times\mathbb{S}^{p-1}_{N}(\leavevmode\hbox
to9.34pt{\vbox to9.34pt{\pgfpicture\makeatletter\hbox{\hskip
4.66791pt\lower-2.66791pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}).$ (14a)

For $p>1$ it holds that $V^{h}$ and $S^{h}$ are conforming subspaces of $V$
and $S$, respectively. The Galerkin mixed formulation reads: find
$\boldsymbol{u}^{h}=(u^{h},w^{h})\in V^{h}$,
$\boldsymbol{\varepsilon}^{h}=(\varepsilon^{h},\chi^{h})\in S^{h}$ and
$\lambda^{h}\in\mathbb{R}^{+}_{0}$ such that

$\displaystyle-c(\boldsymbol{\varepsilon}^{h},\,\boldsymbol{\delta\varepsilon}^{h})+c(B(\boldsymbol{u}^{h}),\,\boldsymbol{\delta\varepsilon}^{h})$
$\displaystyle=0$ $\displaystyle\forall\boldsymbol{\delta\varepsilon}^{h}$
$\displaystyle=(\delta\varepsilon^{h},\delta\chi^{h})\in S^{h}$ (15a)
$\displaystyle c(\boldsymbol{\varepsilon}^{h},\,B(\boldsymbol{\delta u}^{h}))$
$\displaystyle=\lambda^{h}\,b(\boldsymbol{u}^{h},\,\boldsymbol{\delta u}^{h})$
$\displaystyle\forall\boldsymbol{\delta u}^{h}$ $\displaystyle=(\delta
u^{h},\delta w^{h})\in V^{h}.$ (15b)

The mixed Galerkin formulation (15) results in a $4N\times 4N$ system of
equations. As shown in B, this system can be reduced to a $2N\times 2N$ system
in standard form using static condensation.

###  3.3 The standard formulation

As usual, the Galerkin method restricts the formulation to the finite
dimensional setting as follows: find $\boldsymbol{u}^{h}=(u^{h},w^{h})\in
V^{h}$ and $\lambda^{h}\in\mathbb{R}^{+}_{0}$ such that

$\displaystyle a(\boldsymbol{u}^{h},\,\boldsymbol{\delta
u}^{h})=\lambda^{h}\,b(\boldsymbol{u}^{h},\,\boldsymbol{\delta
u}^{h})\qquad\forall\boldsymbol{\delta u}^{h}=(\delta u^{h},\delta w^{h})\in
V^{h}.$ (16)

The standard formulation leads directly to a $2N\times 2N$ matrix eigenvalue
problem, see B. Although the standard and mixed formulation are equivalent at
the continuous level, they are not in the discrete setting. The mixed Galerkin
method offers additional flexibility in the discretization of $S^{h}$.

###  3.4 Analytical computation of discrete eigenvalues and eigenfunctions

The matrix system corresponding to the standard and mixed Galerkin
discretization are discussed in B.1 and B.2, respectively. The Fourier modes
decouple the $2N\times 2N$ system of equations to $N$ systems of $2\times 2$
equations, as discussed in A.2 and A.3. The eigenvalues and amplitudes of the
eigenfunctions satisfy, for each $n$, a two by two matrix eigenvalue problem

$\displaystyle\begin{pmatrix}\tilde{A}^{h}_{n}-\lambda^{h}_{n}&\tilde{B}^{h}_{n}\\\
(\tilde{B}_{n}^{h})^{\ast}&\tilde{D}^{h}_{n}-\lambda^{h}_{n}\end{pmatrix}\begin{pmatrix}U^{h}_{n}\\\
W^{h}_{n}\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix}.$ (17)

The values $\tilde{A}^{h}_{n}$, $\tilde{B}^{h}_{n}$ and $\tilde{D}^{h}_{n}$
have been analytically computed for both the standard and mixed formulation
(for $p=2,3,4$) using the approach outlined in B. They are constants that
depend on $n$ and $N$, the thickness $t$, and radius $R$. For each $n$ there
are two discrete eigenvalues

$\displaystyle\lambda^{h}_{n+}$
$\displaystyle=K^{h}_{n}\left(1+L^{h}_{n}\right)$ (18a)
$\displaystyle\lambda^{h}_{n-}$
$\displaystyle=K^{h}_{n}\left(1-L^{h}_{n}\right)$ (18b)

where $K^{h}_{n}$ and $L^{h}_{n}$ are derived in A.3 in terms of
$\tilde{A}^{h}_{n},\tilde{B}^{h}_{n}$ and $\tilde{D}^{h}_{n}$. Again, the
ratio of the amplitudes of the eigenfunctions, $U^{h}_{n}/W^{h}_{n}$,
satisfies

$\displaystyle\frac{U^{h}_{n+}}{W^{h}_{n+}}=\frac{\tilde{B}^{h}_{n}}{\lambda^{h}_{n+}-\tilde{A}^{h}_{n}}$
(19a)
$\displaystyle\frac{U^{h}_{n-}}{W^{h}_{n-}}=\frac{\tilde{B}^{h}_{n}}{\lambda^{h}_{n-}-\tilde{A}^{h}_{n}}$
(19b)

The discrete eigenvalues and the discrete ratio of the amplitudes may be
directly compared with the analytical solutions developed in the previous
section. Figure 11 illustrates the global behavior of the approximation of the
eigenvalues $\lambda_{n-}$ and relative amplitudes $U_{n-}/W_{n-}$ for a range
of meshes at a normalized thickness of $\bar{t}=0.1$. We can make the
following general observations

  * 1.

The eigenvalues, particularly $\lambda_{-}$, cover many orders of magnitude.
It should be challenging to capture this with a numerical method,
particularly, maintaining accuracy in the very small values. As we show later,
the standard formulation struggles to accurately capture these small
eigenvalues, which is a manifestation of membrane locking. Coarse
discretizations suffer more than fine ones.

  * 2.

At a discretization with $N>\hat{n}$ degrees of freedom, it is challenging to
capture the two different types of physical behavior exhibited by the exact
solution. We show later that the asymptotic accuracy in the amplitude ratio
reduces, in relative terms, compared with results on coarse meshes. Hence,
this phenomenon affects fine discretizations more than coarser ones.

In Section 5 we take a deeper look into the normalized eigenvalue and relative
amplitude errors. In Section 6 we perform a mathematical analysis that sheds
light on some of these observations.

(a) $\lambda^{h}_{-}$ \- standard formulation (b) $\lambda^{h}_{-}$ \- mixed
formulation

(c) $\rho^{h}_{n-}$ \- standard formulation (d) $\rho^{h}_{n-}$ \- mixed
formulation

$N=16$$N=32$$N=64$$N=128$$N=256$$N=512$exact

Figure 11: Approximation of eigenvalues $\lambda_{-}$ and relative amplitudes
$\rho_{n-}=U_{n-}/W_{n-}$ using $p=2$ at a normalized thickness of
$\bar{t}=0.1$.

##  4 A criterion for assessing membrane locking

In this section we present a criterion to assess the presence of membrane
locking in thin beam and shell formulations. It is motivated by the lack of
invariance of normalized eigenvalue spectra as a function of the normalized
mode number. First, we show that the isolated membrane response or isolated
bending response, modeled via a second order or fourth order eigenvalue
problem, respectively, has a discrete normalized eigenvalue error that is
invariant, that is, it is independent of the number of degrees of freedom of
the discrete space. Then, we show an example where the combined membrane and
bending response of a curved Euler-Bernoulli beam lead to a discrete
normalized eigenvalue spectrum that lacks this invariance. This motivates the
choice of the presented criterion to assess membrance locking rigorously in
the end of this section.

###  4.1 Invariance of the isolated membrane response

Let  denote the unit circle and consider the space of periodic functions
$V=H^{1}(\leavevmode\hbox to9.34pt{\vbox
to9.34pt{\pgfpicture\makeatletter\hbox{\hskip 4.66791pt\lower-2.66791pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}})$. We seek $u\in V$
and $\lambda\in\mathbb{R}$ such that

$\displaystyle(u_{,\theta},\,v_{,\theta})=\lambda(u,v)\quad\forall v\in V$

Assuming analytical modes of the form $u_{n}(\theta)=A_{n}\cos(n\theta)$,
$\theta\in[0,2\pi)$, we may determine the analytical eigenvalues as

$\displaystyle\lambda_{n}=n^{2},\quad n=0,1,\ldots$

Let $V^{h}\subset V$ denote a space of periodic, quadratic splines of
dimension $2N$, defined on a uniform partition of the unit circle with mesh
size $h=\pi/N$. It can be shown that the matrix eigenvalue problem is
represented by the $2N$ equations

$\displaystyle\frac{1}{6h}\left(-u_{A-2}-2u_{A-1}+6u_{A}-2u_{A+1}-u_{A+2}\right)-\frac{\lambda^{h}h}{120}\left(u_{A-2}+26u_{A-1}+66u_{A}+26u_{A+1}+u_{A+2}\right)=0$

for $A=1,2,\ldots,2N$. The discrete eigenvalues can then be computed using
elementary properties of circulant matrices. The result is

$\displaystyle\lambda^{h}_{n}=\frac{20}{h^{2}}\cdot\frac{6-4\cos{(nh)}-2\cos{(2nh)}}{66+52\cos{(nh)}+2\cos{(2nh)}},\qquad
n=0,1,\ldots,N-1.$

The error in the non-zero eigenvalues may then be computed analytically as

$\displaystyle\frac{\lambda^{h}_{n}-\lambda_{n}}{\lambda_{n}}=\frac{20}{(nh)^{2}}\cdot\frac{6-4\cos{(nh)}-2\cos{(2nh)}}{66+52\cos{(nh)}+2\cos{(2nh)}}-1.$

In classical eigenvalue error analysis [15, Chapter 6], one proceeds then by
fixing the mode number $n$ and deriving the order of approximation in terms of
the mesh size. Here it can be shown that the discrete eigenvalues converge at
a rate of $h^{4}$ for fixed $n$. We are, however, more interested in the
behavior of the eigenvalue error at fixed normalized mode number,
$\xi=n/N\in[0,1]$. Using $nh=\pi\xi$, the eigenvalue error may be expressed in
terms of the normalized mode number as

$\displaystyle\frac{\lambda^{h}(\xi)-\lambda(\xi)}{\lambda(\xi)}=\frac{20}{(\pi\xi)^{2}}\cdot\frac{6-4\cos{(\pi\xi)}-2\cos{(2\pi\xi)}}{66+52\cos{(\pi\xi)}+2\cos{(2\pi\xi)}}-1.$
(20)

_Note there is no dependence on mesh size!_ Hence, the normalized discrete
frequency spectrum is invariant, that is independent of the number of degrees
of freedom $2N$. Figure 12a illustrates this observation.

###  4.2 Invariance of the isolated bending response

To show that this is not merely a coincidence, we apply the same analysis to a
fourth order eigenvalue problem. Consider periodic functions in
$V=H^{2}(\leavevmode\hbox to9.34pt{\vbox
to9.34pt{\pgfpicture\makeatletter\hbox{\hskip 4.66791pt\lower-2.66791pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{}{{{}} {}{}{}{}{}{}{}{}
}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@moveto{4.26791pt}{2.0pt}\pgfsys@curveto{4.26791pt}{4.35713pt}{2.35713pt}{6.26791pt}{0.0pt}{6.26791pt}\pgfsys@curveto{-2.35713pt}{6.26791pt}{-4.26791pt}{4.35713pt}{-4.26791pt}{2.0pt}\pgfsys@curveto{-4.26791pt}{-0.35713pt}{-2.35713pt}{-2.26791pt}{0.0pt}{-2.26791pt}\pgfsys@curveto{2.35713pt}{-2.26791pt}{4.26791pt}{-0.35713pt}{4.26791pt}{2.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}})$. We seek $u\in V$
and $\lambda\in\mathbb{R}$ such that

$\displaystyle(u_{,\theta\theta},\,v_{,\theta\theta})=\lambda(u,v)\quad\forall
v\in V$

Assuming analytical modes of the form $u_{n}(\theta)=A_{n}\cos(n\theta)$,
$\theta\in[0,2\pi)$, we may determine the analytical eigenvalues as

$\displaystyle\lambda_{n}=n^{4},\quad n=0,1,\ldots$

Let $V^{h}\subset V$ denote the same space of periodic, quadratic splines of
dimension $2N$, defined on a uniform partition of the unit circle with mesh
size $h=\pi/N$. It can be shown that the matrix eigenvalue problem is
represented by the $2N$ equations

$\displaystyle\frac{1}{h^{3}}\left(u_{A-2}-4u_{A-1}+6u_{A}-4u_{A+1}+u_{A+2}\right)-\frac{\lambda^{h}h}{120}\left(u_{A-2}+26u_{A-1}+66u_{A}+26u_{A+1}+u_{A+2}\right)=0$

for $A=1,2,\ldots,2N$. The discrete eigenvalues are computed as

$\displaystyle\lambda^{h}_{n}=\frac{120}{h^{4}}\cdot\frac{6-8\cos{(nh)}+2\cos{(2nh)}}{66+52\cos{(nh)}+2\cos{(2nh)}},\qquad
n=0,1,\ldots,N-1.$

Analoguous to the second order eigenvalue problem, using $nh=\pi\xi$, the
eigenvalue error may be expressed as

$\displaystyle\frac{\lambda^{h}(\xi)-\lambda(\xi)}{\lambda(\xi)}=\frac{120}{(\pi\xi)^{4}}\cdot\frac{6-8\cos{(\pi\xi)}+2\cos{(2\pi\xi)}}{66+52\cos{(\pi\xi)}+2\cos{(2\pi\xi)}}-1.$
(21)

Importantly, once again, the expression is invariant with respect to mesh
refinement! This is illustrated in Figure 12b.

(a) Second order eigenvalue problem (b) Fourth order eigenvalue problem

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 12: Invariance of spectra as a function of the normalized mode number.

###  4.3 Definition of a rigorous criterion based on spectrum invariance

The isolated membrane and bending responses are invariant when viewed as
functions of the normalized mode number $\xi$. This may no longer be the case
when they are coupled in a curved Euler-Bernoulli beam. Figure 13 depicts a
preview of the normalized eigenvalue error obtained with quadratic splines.
The eigenvalue error is clearly not invariant with respect to the number of
degrees of freedom of the spline space. Particularly, the lower eigenvalues
suffer from a visible discrepancy compared to the rest, which can manifest
partly by a larger preasymptotic regime before converging. We identify this as
membrane locking. The previous discussion motivates the following
characterization of membrane locking

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 13: Membrane locking is manifested by a lack of invariance of the
discrete normalized spectrum. The discrepancy between the spectrum error and
the asymptotic spectrum error, denoted by $d$, measures the extent of locking.

> _Membrane locking is measured, at fixed normalized mode number, as the
> discrepancy between the spectrum error and the spectrum error in the limit
> of asymptotic refinement._

To make this idea more precise, use $n=\xi N$, and let

$\displaystyle e_{N}(\xi):=\frac{|\lambda^{h}_{\xi N}-\lambda_{\xi
N}|}{\lambda_{\xi N}}.$ (22)

In mathematical notation, the above statement may then be translated as
follows.

######  Definition 4.1 (Spectral characterization of membrane locking ).

Let $0<\varepsilon<<1$ denote a user-prescribed tolerance. A discretization is
characterized as locking-free if and only if

$\displaystyle\log_{10}{\left(e_{N}(\xi)\right)}-\lim_{\mathbb{Z}_{N}\ni
M\rightarrow\infty}\log_{10}{\left(e_{M}(\xi)\right)}<\varepsilon\qquad\forall\xi=n/N,\qquad
n=1,2,\ldots,N.$ (23)

The condition that $M\in\mathbb{Z}_{N}$, the set of integers modulo $N$, is
needed to ensure there exists an integer $m$ such that $m/M=n/N=\xi$. The
tolerance $\varepsilon$ is a distance in log-space, say $\varepsilon=0.01$,
beyond which we cannot observe (in log-space) a noticeable difference between
the spectrum error and its asymptotic result. Figure 14 depicts the distance
between the spectrum error and the asymptotic spectrum error,
$d=\log_{10}{\left(e_{N}(\xi)\right)}-\log_{10}{\left(e_{\infty}(\xi)\right)}$,
corresponding to the results in Figure 13 obtained using the standard
formulation with $p=2$. According to our characterization of locking, locking
is exhibited by all eigenmodes associated to the negative eigenvalue branch
for discretizations with $2N=32$ and $2N=64$ degrees of freedom. In
discretizations with $2N=128$ and $2N=256$, still about $20\%$ and $10\%$ of
the modes lock, respectively.

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 14: Criterion to assess locking according to definition 4.1 with
$\varepsilon=1e-2$. Values in the gray area denote mode numbers for which
there is measurable membrane locking.

##  5 Spectral analysis results

In this section, we assess the extent of membrane locking in the normalized
eigenvalue spectra of both the standard and mixed formulation, employing the
locking criterion outlined in Section 4. While our primary focus lies on
quadratic discretizations, we also present findings that validate membrane
locking in cubic and quartic discretizations. Furthermore, we explore
intriguing behaviors that manifest in fine discretizations and mode
amplitudes.

###  5.1 Presentation of the data

We present the results pertaining to the two types of eigenvalues, namely
$\lambda_{-}$ and $\lambda_{+}$, for both the standard and mixed formulation.
In certain instances, we focus on results associated with the eigenvalues
$\lambda_{-}$ of the standard formulation, as they exhibit pronounced locking
effects. The presented figures share the following attributes:

  * 1.

Normalized eigenvalue errors and mode amplitudes are shown for meshes with
$2N=32,\,64,\,128$, and $256$ degrees of freedom. The breathing mode ($n=0$)
and the rigid body mode ($n=1$), see Figure 4, are exact or well approximated
on any mesh and are not plotted. The remaining $N-2$ eigenvalues have
multiplicity two. Only the unique ones are plotted.

  * 2.

Eigenvalue error and mode amplitudes are normalized with respect to the
normalized mode number $\xi=n/N$.

  * 3.

Four asymptotic curves are drawn in each graph that allow a succinct
comparison: two for each method and two corresponding to each type of
eigenvalue. The asymptotic curve under consideration is plotted as a solid
curve, while the other three are dotted curves.

  * 4.

All results are obtained for a normalized thickness $\bar{t}=0.015$.

###  5.2 Normalized error in eigenvalues

Figure 15 showcases the normalized eigenvalue error on a logarithmic scale. As
discussed in the previous section, the normalized eigenvalue error,
$\frac{\lambda^{h}_{-}-\lambda_{-}}{\lambda_{-}}$, of the standard formulation
reveals significant membrane locking, as the errors deviate considerably from
the asymptotic limit curves. Conversely, in the mixed method, these
eigenvalues are accurately approximated, exhibiting minimal signs of locking.
Additionally, it is worth noting that the asymptotic accuracy of the mixed
method surpasses that of the standard formulation, evident from a comparison
of the limit curves between the two methods. For the eigenvalues
$\lambda_{+}$, both methods yield equally satisfactory approximations, with no
indications of locking as the results align with or slightly fall below the
limit curve.

(a) $\lambda^{h}_{-}$ \- Standard formulation (b) $\lambda^{h}_{-}$ \- Mixed
formulation

(c) $\lambda^{h}_{+}$ \- Standard formulation (d) $\lambda^{h}_{+}$ \- Mixed
formulation

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 15: Normalized eigenvalue error obtained with (a) standard and (b)
mixed quadratic discretizations.

By calculating the logarithmic distance between the eigenvalue errors and the
asymptotic curves, we can rigorously assess the extent of membrane locking
using the criterion outlined in the previous section. Figure 16 illustrates
the results. As previously discussed, the standard formulation exhibits
significant locking for coarse to medium refined discretizations. In contrast,
the mixed method demonstrates only minor locking in the first three
eigenvalues, regardless of the level of refinement.

(a) Standard formulation (b) Mixed formulation

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 16: Identification of eigenvalues $\lambda^{h}_{-}$ that display
membrane locking ($\varepsilon=1\cdot 10^{-2}$) in quadratic discretizations.

Figure 17 and 18 illustrate the severity of locking of the standard
formulation for cubics and quartics. Although the absolute accuracy improves
dramatically with polynomial degree, the severity of membrane locking does not
improve concomitantly. The degree of membrane locking is the same as for
quadratics. Note also that the absolute accuracy of the mixed method again far
exceeds that of the standard formulation, as may be observed from the
asymptotic curves in Figure 17.

(a) Standard formulation $p=3$ (b) Standard formulation $p=4$

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 17: Normalized error in the eigenvalues $\lambda^{h}_{-}$ for the
standard discretization obtained with cubics and quartics.

(a) Standard formulation $p=3$ (b) Standard formulation $p=4$

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 18: Identification of eigenvalues $\lambda^{h}_{-}$ that display
membrane locking ($\varepsilon=1\cdot 10^{-2}$) in cubic and quartic
discretizations.

###  5.3 Normalized error in mode amplitudes

Figure 19 illustrates the normalized error in the amplitudes of the eigenmodes
corresponding to the negative branch. Both methods exhibit the same level of
accuracy when applied to meshes with $2N=32,64,128$, and $256$ degrees of
freedom, consistent with earlier findings reported in [9]. However, the
asymptotic limit curves demonstrate distinct behavior. In relative terms, the
approximations deteriorate with increasing refinement, especially noticeable
in the standard formulation. This interesting phenomenon might be another
manifestation of membrane locking, occurring specifically in the higher modes.
It is important to note that the positive branch displays an identical pattern
of behavior.

(a) Standard formulation (b) Mixed formulation

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 19: Normalized error in the mode amplitudes of the negative branch,
$\rho^{h}_{n-}=U_{n-}^{h}/W_{n-}^{h}$, obtained with (a) standard and (b)
mixed quadratic discretizations.

###  5.4 Behaviour of normalized eigenvalue errors for fine discretizations

In addition to our investigation into membrane locking occurring in coarse to
medium refined discretizations, we present interesting and previously
unobserved findings that occur in fine discretizations of the circular Euler-
Bernoulli ring, see Figure 20. Notably, for discretizations where $N>\hat{n}$,
the curves of normalized eigenvalue errors deviate from their asymptotic
pattern. Although initially surprising, this behaviour can be anticipated upon
examining the relative amplitudes presented in Figures 5 and 11. We should
recall that for $n<\hat{n}$, the modes corresponding to $\lambda_{n+}$ are
membrane-dominated ($U_{n+}/W_{n+}>1$), whereas the modes corresponding to
$\lambda_{n-}$ are dominated by bending effects ($U_{n-}/W_{n-}<1$). However,
this behaviour is reversed for $n>\hat{n}$. While we specifically demonstrate
this behavior using quadratic discretizations of the standard formulation, it
holds true for polynomial degrees of any order and is also applicable to the
mixed formulation.

(a) Standard formulation - $\lambda^{h}_{-}$ (b) Standard formulation -
$\lambda^{h}_{+}$

$N=16$$N=32$$N=64$$N=128$$N\rightarrow\infty$

Figure 20: Normalized eigenvalue error obtained with quadratic, fine
discretizations of the standard formulation.

##  6 Mathematical analysis

In this section, we present a detailed mathematical analysis of the error
between the discrete eigenvalues and their exact counterparts as a function of
the relative mode number, the mesh size, and the relevant physical parameters.
This relative error is measured by the difference between $\lambda^{h}_{n\pm}$
and $\lambda_{n\pm}$ normalized by $\lambda_{n\pm}$ with $n=\xi N=\xi\pi/h$
(see (22) in Section 4.3). Our primary focus is on the standard formulation,
where membrane locking is more pronounced. However, our analysis is
comprehensive and applies to general polynomial degrees and could be adapted
to the mixed method as well. In light of our characterization of membrane
locking, we are particularly interested in studying the behavior of the
relative error as $h$ approaches zero while keeping $\xi$ fixed. Additionally,
we extend our investigation to explore the relationship between eigenvalue
error and normalized thickness, thereby providing a broader understanding of
the phenomenon.

###  6.1 Relative eigenvalue error as a function of the normalized mode number

We will focus on the expression

$\frac{\lambda^{h}_{n\pm}-\lambda_{n\pm}}{\lambda_{n\pm}}=\frac{\lambda^{h}_{n\pm}}{\lambda_{n\pm}}-1=\frac{K_{n}^{h}}{K_{n}}\left[\frac{1\pm\sqrt{1-M_{n}^{h}/(K_{n}^{h})^{2}}}{1\pm\sqrt{1-M_{n}/(K_{n})^{2}}}\right]-1\,,$
(24)

where we take the values of $\lambda_{n\pm}$ and $\lambda_{n\pm}^{h}$ from (5)
and (18).

First, note that rewriting $K_{n}$ in (6) with $n=\tfrac{\xi\pi}{h}$ yields
(25a). On the other hand, a careful analysis of B.1 yields an analogous
expression for the discrete case (25b),

$\displaystyle K_{n}$
$\displaystyle=\frac{1}{2}\frac{1}{R^{2}}\frac{1}{h^{4}}\left(\alpha_{0}+\alpha_{2}h^{2}+h^{4}\right),$
(25a) $\displaystyle K_{n}^{h}$
$\displaystyle=\frac{1}{2}\frac{1}{R^{2}}\frac{1}{h^{4}}\left(\alpha_{0}^{h}+\alpha_{2}^{h}h^{2}+h^{4}\right),$
(25b)

with

$\displaystyle\alpha_{0}$ $\displaystyle=\frac{\bar{t}^{2}}{12}(\xi\pi)^{4},$
$\displaystyle\alpha_{0}^{h}$
$\displaystyle=\frac{\bar{t}^{2}}{12}\frac{\tilde{\mathsf{K}}_{b}^{22}}{\tilde{\mathsf{M}}},$
(26a) $\displaystyle\alpha_{2}$
$\displaystyle=\left(1+\frac{\bar{t}^{2}}{12}\right)(\xi\pi)^{2},$
$\displaystyle\alpha_{2}^{h}$
$\displaystyle=\left(1+\frac{\bar{t}^{2}}{12}\right)\frac{\tilde{\mathsf{K}}_{b}^{11}}{\tilde{\mathsf{M}}}.$
(26b)

Here, the $\tilde{\mathsf{K}}_{b}^{ij}$ and $\tilde{\mathsf{K}}_{m}^{ij}$
represent the eigenvalues of the “normalized” circulant matrices, as presented
in (62), by which we mean that matrices have been stripped of any dependence
on $\bar{t}$, $R$, and $h$. Similarly, $\tilde{\mathsf{M}}$ is the eigenvalue
of the normalized circulant symmetric matrix described by (61). More
precisely, these are the eigenvalues corresponding to the matrices reported in
Table 1.

As highlighted in Section A, these eigenvalues can be analytically computed
and expressed as weighted sums of sines and cosines of integer multiples of
$\xi\pi=nh=n\tfrac{2\pi}{2N}$ (see equations (48) and (49)). This key
observation implies that the eigenvalues are solely determined by the
normalized mode number $\xi$ and are independent of other parameters. This
point holds significant importance in understanding the behavior of the
eigenvalues as functions of the normalized mode number.

A similar approach to $M_{n}$ and $M_{n}^{h}$ yield the following expressions,

$\displaystyle M_{n}$
$\displaystyle=\frac{1}{4}\frac{1}{R^{4}}\frac{1}{h^{6}}\left(a_{0}+a_{2}h^{2}+a_{4}h^{4}\right),$
(27a) $\displaystyle M_{n}^{h}$
$\displaystyle=\frac{1}{4}\frac{1}{R^{4}}\frac{1}{h^{6}}\left(a_{0}^{h}+a_{2}^{h}h^{2}+a_{4}^{h}h^{4}\right),$
(27b)

with

$\displaystyle a_{0}$ $\displaystyle=\frac{4\bar{t}^{2}}{12}\,(\xi\pi)^{6},$
$\displaystyle a^{h}_{0}$
$\displaystyle=\frac{4\bar{t}^{2}}{12}\frac{1}{\tilde{\mathsf{M}}^{2}}\,\Bigg{[}\bigg{(}1+\frac{\bar{t}^{2}}{12}\bigg{)}\tilde{\mathsf{K}}_{b}^{11}\tilde{\mathsf{K}}_{b}^{22}-\frac{\bar{t}^{2}}{12}\big{|}\tilde{\mathsf{K}}_{b}^{12}\big{|}^{2}\Bigg{]},$
(28a) $\displaystyle a_{2}$
$\displaystyle=-\frac{4\bar{t}^{2}}{12}\,2(\xi\pi)^{4},$ $\displaystyle
a^{h}_{2}$
$\displaystyle=-\frac{4\bar{t}^{2}}{12}\frac{1}{\tilde{\mathsf{M}}^{2}}\,2\tilde{\mathsf{K}}_{m}^{12}(\tilde{\mathsf{K}}_{b}^{12})^{\ast},$
(28b) $\displaystyle a_{4}$
$\displaystyle=\frac{4\bar{t}^{2}}{12}\,(\xi\pi)^{2},$ $\displaystyle
a^{h}_{4}$
$\displaystyle=\frac{4\bar{t}^{2}}{12}\frac{1}{\tilde{\mathsf{M}}^{2}}\,\frac{12}{\bar{t}^{2}}\Bigg{[}\Bigg{(}1+\frac{\bar{t}^{2}}{12}\Bigg{)}\tilde{\mathsf{K}}_{b}^{11}\tilde{\mathsf{K}}_{m}^{22}-\big{|}\tilde{\mathsf{K}}_{m}^{12}\big{|}^{2}\Bigg{]}.$
(28c)

Within this context we note that
$\tilde{\mathsf{K}}_{m}^{12}(\tilde{\mathsf{K}}_{b}^{12})^{\ast}$ (with
$(\cdot)^{\ast}$ denoting complex conjugation) is a real number in view of the
fact that they are purely imaginary eigenvalues coming from skew-symmetric
matrices.

Meanwhile, it follows directly from (25a) and (25b) that

$\displaystyle\begin{aligned}
K_{n}^{2}&=\frac{1}{4}\frac{1}{R^{4}}\frac{1}{h^{8}}\Big{(}\alpha_{0}^{2}+2\alpha_{0}\alpha_{2}h^{2}+\big{(}2\alpha_{0}+\alpha_{2}^{2}\big{)}h^{4}+2\alpha_{2}h^{6}+h^{8}\Big{)}\\\
&=\frac{1}{4}\frac{1}{R^{4}}\frac{1}{h^{8}}\Big{(}b_{0}+b_{2}h^{2}+b_{4}h^{4}+b_{6}h^{6}+h^{8}\Big{)}\end{aligned}$
(29a) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\begin{aligned}
(K_{n}^{h})^{2}&=\frac{1}{4}\frac{1}{R^{4}}\frac{1}{h^{8}}\Big{(}(\alpha_{0}^{h})^{2}+2\alpha_{0}^{h}\alpha_{2}^{h}h^{2}+\big{(}2\alpha_{0}^{h}+(\alpha_{2}^{h})^{2}\big{)}h^{4}+2\alpha_{2}^{h}h^{6}+h^{8}\Big{)}\\\
&=\frac{1}{4}\frac{1}{R^{4}}\frac{1}{h^{8}}\Big{(}b_{0}^{h}+b_{2}^{h}h^{2}+b_{4}^{h}h^{4}+b_{6}^{h}h^{6}+h^{8}\Big{)}\,.\end{aligned}$
(29b)

Putting everything together gives,

$\frac{\lambda^{h}_{n\pm}}{\lambda_{n\pm}}=\frac{K_{n}^{h}}{K_{n}}\left[\frac{1\pm\sqrt{1-\frac{M_{n}^{h}}{(K_{n}^{h})^{2}}}}{1\pm\sqrt{1-\frac{M_{n}}{K_{n}^{2}}}}\right]=\frac{\alpha_{0}^{h}+\alpha_{2}^{h}h^{2}+h^{4}}{\alpha_{0}+\alpha_{2}h^{2}+h^{4}}\left[\frac{1\pm\sqrt{1-h^{2}\frac{a_{0}^{h}+a_{2}^{h}h^{2}+a_{4}^{h}h^{4}}{b_{0}^{h}+b_{2}^{h}h^{2}+b_{4}^{h}h^{4}+b_{6}^{h}h^{6}+h^{8}}}}{1\pm\sqrt{1-h^{2}\frac{a_{0}+a_{2}h^{2}+a_{4}h^{4}}{b_{0}+b_{2}h^{2}+b_{4}h^{4}+b_{6}h^{6}+h^{8}}}}\right]\,.$
(30)

The main observation at this point is that there is no dependence on $R$
outside of the normalized thickness $\bar{t}=\tfrac{t}{R}$ (more generally
outside of $\bar{I}=\tfrac{I}{AR^{2}}=\frac{t^{2}}{12R^{2}}$).

###  6.2 Relative eigenvalue error in the asymptotic limit of refinement

We are particularly concerned with studying the relative error with respect to
the normalized mode number $\xi\in(0,1]$ in the asymptotic refinement limit,
where $h$ approaches zero. As previously mentioned, we hypothesize that
membrane locking is responsible for the disparity between the spectrum
computed at a fixed $h$ and the spectrum in this asymptotic limit.

The expression in (30) lends itself very well to taking this limit in both the
positive and negative branches of eigenvalues. For the positive branch, taking
the limit $h\to 0$ gives

$\lim_{h\to
0}\frac{\lambda^{h}_{n+}}{\lambda_{n+}}-1=\frac{\alpha_{0}^{h}}{\alpha_{0}}\left[\frac{1+\sqrt{1}}{1+\sqrt{1}}\right]-1=\frac{1}{(\xi\pi)^{4}}\frac{\tilde{\mathsf{K}}_{b}^{22}}{\tilde{\mathsf{M}}}-1\,,$
(31)

where the values of $\alpha_{0}$ and $\alpha_{0}^{h}$ were taken from (25a)
and (25b). It is important to observe that the normalized thickness does not
affect this expression as it cancels out. In simpler terms, this expression
solely relies on $\xi$ and the polynomial order $p$, which determine
$\tilde{\mathsf{K}}_{b}^{22}$ and $\tilde{\mathsf{M}}$.

With that being said, it is crucial to further examine this expression. For
the case of $p=2$, we can compute the circulant eigenvalues of
$\tilde{\mathsf{K}}_{b}^{22}$ and $\tilde{\mathsf{M}}$ using (48) along with
the values provided in Table 1. The resulting expression is given by:

$\lim_{h\to
0}\frac{\lambda^{h}_{n+}}{\lambda_{n+}}-1=\frac{60}{(\xi\pi)^{4}}\cdot\frac{6-8\cos{(\xi\pi)}+2\cos{(2\xi\pi)}}{33+26\cos{(\xi\pi)}+\cos{(2\xi\pi)}}-1.$
(32)

Remarkably, this expression exactly matches (21), which represents the
relative error resulting from solving the isolated bending response. Notably,
this correspondence holds generally true for any value of $p$. In essence, it
implies that in the limit of mesh refinement, the relative error of the
positive branch coincides with that observed in the isolated bending response!

To analyze the negative branch of eigenvalues, we employ the Taylor series to
avoid encountering a limit of the form $\frac{0}{0}$. However, in the end, the
following expression holds in that limit:

$\lim_{h\to
0}\frac{\lambda^{h}_{n-}}{\lambda_{n-}}-1=\frac{\alpha_{0}^{h}}{\alpha_{0}}\frac{a_{0}^{h}}{a_{0}}\frac{b_{0}}{b_{0}^{h}}-1=\frac{1}{(\xi\pi)^{2}}\left[\left(1+\frac{\bar{t}^{2}}{12}\right)\frac{\tilde{\mathsf{K}}_{m}^{11}}{\tilde{\mathsf{M}}}-\frac{\bar{t}^{2}}{12}\frac{\big{|}\tilde{\mathsf{K}}_{b}^{12}\big{|}^{2}}{\tilde{\mathsf{M}}\tilde{\mathsf{K}}_{b}^{22}}\right]-1\,.$
(33)

Here, we utilized the fact that
$\tilde{\mathsf{K}}_{m}^{11}=\tilde{\mathsf{K}}_{b}^{11}$. It is evident that
this expression is dependent on $\bar{t}$, and as $\bar{t}\to 0$, the relative
error becomes
$\frac{1}{(\xi\pi)^{2}}\frac{\tilde{\mathsf{K}}_{m}^{11}}{\tilde{\mathsf{M}}}-1$,
which precisely corresponds to the isolated membrane response, regardless of
the value of $p$. Indeed, for $p=2$, it is not difficult to notice that the
above expression coincides with (20) in the limit of $\bar{t}\to 0$.

###  6.3 Relation to classical eigenvalue error analysis

Typically, convergence in eigenvalue problems is studied by fixing a mode
number $n$ and refining the mesh. Interestingly, our analysis involving the
normalized mode number $\xi=\frac{n}{N}=\frac{nh}{\pi}$ incorporates this
aspect naturally. This is due to the fact that for a fixed mode number $n$, as
we decrease the mesh size $h$ (equivalently, increase the number of elements
$N$), the value $nh=\xi\pi$ becomes small. In other words, when we focus on
small values of $\xi$, we essentially recover the behavior of the relative
error at small $h$ for any fixed $n$. This implies that by examining the
relative error curves in the asymptotic limit of refinement, as derived in
(31) and (33), we can gain insights into the behavior of the relative error
for any fixed $n$.

To illustrate this point, let us consider the case of $p=2$. First, examine
the bending-dominated asymptotic limit of refinement described by (32). When
$\xi\pi$ is small (i.e., $\xi\pi\ll 1$), the expression should approach zero.
However, upon closer inspection, we notice that it exhibits an _inverse_ power
of $\xi\pi$, which may initially seem counterintuitive as it suggests a large
value for small $\xi$. This apparent contradiction can be resolved by
analyzing the expression more carefully using its Taylor series expansion,
which allows us to account for the cancellation of these inverse powers. This
technique is analogous to the methods employed in finite differences [22,
§3.4]. By performing the necessary calculations, we obtain the following
expression:

$\lim_{h\to
0}\frac{\lambda^{h}_{n+}}{\lambda_{n+}}-1=\frac{60}{(\xi\pi)^{4}}\frac{6-8\cos(\xi\pi)+2\cos(2\xi\pi)}{33+26\cos(\xi\pi)+\cos(2\xi\pi)}-1=\frac{1}{12}(\xi\pi)^{2}+\frac{1}{240}(\xi\pi)^{4}+\mathcal{O}\big{(}(\xi\pi)^{6}\big{)}\,.$
(34)

Therefore, when considering a fixed $n$ and small enough $h$, the relative
error exhibits a behavior of $\mathcal{O}\big{(}h^{2}\big{)}$. This result
aligns with expectations, as the anticipated convergence rate for eigenvalues
in a fourth-order problem is $\mathcal{O}\big{(}h^{2(p-1)}\big{)}$. It can be
easily verified that this behavior holds true for other values of $p$ as well.

Moving forward, let’s examine the membrane-dominated asymptotic limit of
refinement for the case of $p=2$, as presented in (33). To analyze this
expression, we need to consider the Taylor series expansions of its two
components:

$\displaystyle\frac{1}{(\xi\pi)^{2}}\frac{\tilde{\mathsf{K}}_{m}^{11}}{\tilde{\mathsf{M}}}$
$\displaystyle=\frac{60}{3}\frac{3-2\cos(\xi\pi)-\cos(2\xi\pi)}{33+26\cos(\xi\pi)+\cos(2\xi\pi)}$
(35)
$\displaystyle=1+\frac{1}{720}(\xi\pi)^{4}+\frac{1}{3360}(\xi\pi)^{6}+\mathcal{O}\big{(}(\xi\pi)^{8}\big{)}$
$\displaystyle\frac{1}{(\xi\pi)^{2}}\frac{\big{|}\tilde{\mathsf{K}}_{b}^{12}\big{|}^{2}}{\tilde{\mathsf{M}}\tilde{\mathsf{K}}_{b}^{22}}$
$\displaystyle=\frac{60\big{(}2\sin(\xi\pi)-\sin(2\xi\pi)\big{)}^{2}}{\big{(}6-8\cos(\xi\pi)+2\cos(2\xi\pi)\big{)}\big{(}33+26\cos(\xi\pi)+\cos(2\xi\pi)\big{)}}$
$\displaystyle=1-\frac{1}{12}(\xi\pi)^{2}-\frac{1}{180}(\xi\pi)^{4}+\mathcal{O}\big{(}(\xi\pi)^{6}\big{)}\,.$

It is important to note that the isolated membrane response, represented by
$\frac{1}{(\xi\pi)^{2}}\frac{\tilde{\mathsf{K}}_{m}^{11}}{\tilde{\mathsf{M}}}$,
exhibits eigenvalue convergence of $\mathcal{O}\big{(}h^{4}\big{)}$. This
convergence rate coincides with the expected rate for second-order problems,
which is $\mathcal{O}\big{(}h^{2p}\big{)}$. On the other hand, the bending
component, given by
$\frac{1}{(\xi\pi)^{2}}\frac{|\tilde{\mathsf{K}}_{b}^{12}|^{2}}{\tilde{\mathsf{M}}\tilde{\mathsf{K}}_{b}^{22}}$,
converges at a slower rate of $\mathcal{O}\big{(}h^{2}\big{)}$. Consequently,
the overall convergence rate of the relative error is determined by the
bending component, resulting in a convergence rate of
$\mathcal{O}\big{(}h^{2}\big{)}$ for small $h$. For arbitrary $p$, the
membrane-dominated eigenvalues will converge at a rate of
$\mathcal{O}\big{(}h^{2(p-1)}\big{)}$.

The inclusion of the relative error curves as functions of $\xi$ provides
additional insights. We can determine the value of $\xi^{*}$ below which the
relative error behaves consistently with
$\mathcal{O}\big{(}h^{2(p-1)}\big{)}$. This value serves as a measure of the
preasymptotic domain’s size, indicating the threshold at which the eigenvalues
begin to converge at the expected rate for a given $n$. These $\xi^{*}$ values
(or equivalently, $h^{*}=\tfrac{\xi^{*}\pi}{n}$) may vary depending on the
value of $p$. From the curves of the asymptotic limit of refinement shown in
Figure 15, we can estimate that the convergence for both branches of
eigenvalues begins at $\xi^{*}<0.1$. It is also noteworthy that the mixed
formulation exhibits a distinct asymptotic curve for the bending-dominated
eigenvalues, suggesting an earlier onset of convergence (i.e., at a larger
value of $\xi$).

###  6.4 Considerations for comparing relative errors

To accurately assess the relative errors, it is crucial to consider a subtle
aspect. When $h$ takes on moderate values such that $N<\hat{n}$ (i.e., the
number of elements in the discretization is less than the transition mode
$\hat{n}$), the eigenvalues in the positive branch correspond to a membrane-
dominated response, as depicted in Figure 5. However, as $h$ approaches zero
and the number of elements $N$ becomes much larger than the transition mode
$\hat{n}$ (i.e., $N\gg\hat{n}$), the positive branch of eigenvalues
predominantly consists of bending-dominated eigenvalues (see Figure 20). This
aligns with the observation made earlier regarding the curve derived in (31),
which is associated with the isolated bending response.

To effectively analyze the discrepancy in the relative error of an eigenvalue
at a finite value of $h$, it is essential to compare it against the asymptotic
limit of refinement associated with the same _class of physical response_ ,
rather than focusing solely on the same eigenvalue branch. In particular, at
$\bar{t}=0.01$, the first 100 eigenvalues corresponding to the positive branch
exhibit membrane-dominated behavior (see Figure 5). When comparing their
relative errors, we should consider the membrane-dominated asymptotic limit,
which corresponds to the asymptotic limit of the negative branch of
eigenvalues. Similar considerations apply to the eigenvalues in the negative
branch.

Lastly, it is worth recalling that the transition mode $\hat{n}$ depends
solely on the normalized thickness $\bar{t}$, as discussed in Section 2.3.
This implies that the behavior of the eigenvalues and their dominance by
membrane or bending responses are determined by a combination of mesh
refinement (controlled by $h$, see Figure 20) and the normalized thickness
$\bar{t}$. By considering the interplay between $h$ and $\bar{t}$, we can
accurately interpret and compare eigenvalues with the correct limit curves.
Therefore, in the locking criterion discussed in Section 4.3, which we propose
as a means of determining whether a particular method is affected by membrane
locking, the comparison being made might span across eigenvalue branches. This
detail should be taken into account when examining the definition.

###  6.5 Analysis of the physical parameters

Most physical parameters, such as Young’s modulus, density, and area, have a
multiplicative effect and therefore do not influence the behavior of either
the analytical or discrete solutions. However, the thickness $t$ and radius of
the ring $R$ do have a relative effect on certain terms. In fact, they can be
consolidated into a single parameter, the normalized thickness
$\bar{t}=\frac{t}{R}$ as defined previously, or alternatively, the normalized
moment of inertia $\bar{I}=\frac{I}{AR^{2}}$, where $A$ represents the cross-
sectional area of the ring. To see the equivalence, consider that for a
rectangular cross-section, $\bar{I}=\frac{\bar{t}^{2}}{12}$, and $\bar{I}\to
0$ if and only if $\bar{t}\to 0$. To simplify the resulting formulas as much
as possible, we will consistently use the parameter $\bar{I}$ throughout this
section. However, it should be noted that $\bar{I}$ is merely proportional to
$\bar{t}^{2}$.

The idea is to proceed similarly to Section 6.1, and obtain an expression
analogous to (30). Reorganizing the terms in (25) in terms of $\bar{I}$ yields

$\displaystyle
K_{n}=\frac{1}{2}\frac{1}{R^{2}}\frac{1}{h^{4}}\left(\beta_{0}+\beta_{1}\bar{I}\right),\qquad\qquad\quad
K_{n}^{h}=\frac{1}{2}\frac{1}{R^{2}}\frac{1}{h^{4}}\left(\beta_{0}^{h}+\beta_{1}^{h}\bar{I}\right),$
(36)

with

$\displaystyle\beta_{0}$ $\displaystyle=h^{2}(\xi\pi)^{2}+h^{4},$
$\displaystyle\beta_{0}^{h}$
$\displaystyle=h^{2}\frac{\tilde{\mathsf{K}}_{b}^{11}}{\tilde{\mathsf{M}}}+h^{4},$
(37a) $\displaystyle\beta_{1}$ $\displaystyle=(\xi\pi)^{4}+h^{2}(\xi\pi)^{2},$
$\displaystyle\beta_{1}^{h}$
$\displaystyle=\frac{1}{\tilde{\mathsf{M}}}\left(\tilde{\mathsf{K}}_{b}^{22}+h^{2}\tilde{\mathsf{K}}_{b}^{11}\right).$
(37b)

Proceeding similarly with (27) gives

$\displaystyle
M_{n}=\frac{1}{4}\frac{1}{R^{4}}\frac{1}{h^{8}}\left(\tilde{a}_{0}+\tilde{a}_{1}\bar{I}+\tilde{a}_{2}\bar{I}^{2}\right),\qquad\qquad\quad
M_{n}^{h}=\frac{1}{4}\frac{1}{R^{4}}\frac{1}{h^{8}}\left(\tilde{a}_{0}^{h}+\tilde{a}_{1}^{h}\bar{I}+\tilde{a}_{2}^{h}\bar{I}^{2}\right),$
(38)

with

$\displaystyle\tilde{a}_{0}$ $\displaystyle=0,$
$\displaystyle\tilde{a}_{0}^{h}$
$\displaystyle=\frac{4h^{2}}{\tilde{\mathsf{M}}^{2}}h^{4}\left[\tilde{\mathsf{K}}_{m}^{11}\tilde{\mathsf{K}}_{m}^{22}-\big{|}\tilde{\mathsf{K}}_{m}^{12}\big{|}^{2}\right],$
(39a) $\displaystyle\tilde{a}_{1}$
$\displaystyle=4h^{2}\left((\xi\pi)^{6}-2h^{2}(\xi\pi)^{4}+h^{4}(\xi\pi)^{2}\right),$
$\displaystyle\tilde{a}_{1}^{h}$
$\displaystyle=\frac{4h^{2}}{\tilde{\mathsf{M}}^{2}}\left[\tilde{\mathsf{K}}_{b}^{11}\tilde{\mathsf{K}}_{b}^{22}-2h^{2}\tilde{\mathsf{K}}_{m}^{12}(\tilde{\mathsf{K}}_{b}^{12})^{*}+h^{4}\tilde{\mathsf{K}}_{m}^{11}\tilde{\mathsf{K}}_{m}^{22}\right],$
(39b) $\displaystyle\tilde{a}_{2}$ $\displaystyle=0,$
$\displaystyle\tilde{a}_{2}^{h}$
$\displaystyle=\frac{4h^{2}}{\tilde{\mathsf{M}}^{2}}\left[\tilde{\mathsf{K}}_{b}^{11}\tilde{\mathsf{K}}_{b}^{22}-\big{|}\tilde{\mathsf{K}}_{b}^{12}\big{|}^{2}\right].$
(39c)

Putting everything together yields the following exact expression for the
relative magnitude of the eigenvalues,

$\frac{\lambda^{h}_{n\pm}}{\lambda_{n\pm}}=\frac{K_{n}^{h}}{K_{n}}\left[\frac{1\pm\sqrt{1-\frac{M_{n}^{h}}{(K_{n}^{h})^{2}}}}{1\pm\sqrt{1-\frac{M_{n}}{K_{n}^{2}}}}\right]=\frac{\beta_{0}^{h}+\beta_{1}^{h}\bar{I}}{\beta_{0}+\beta_{1}\bar{I}}\left[\frac{1\pm\sqrt{1-\frac{\tilde{a}_{0}^{h}+\tilde{a}_{1}^{h}\bar{I}+\tilde{a}_{2}^{h}\bar{I}^{2}}{(\beta_{0}^{h})^{2}+2\beta_{0}^{h}\beta_{1}^{h}\bar{I}+(\beta_{1}^{h})^{2}\bar{I}^{2}}}}{1\pm\sqrt{1-\frac{\tilde{a}_{1}\bar{I}}{(\beta_{0})^{2}+2\beta_{0}\beta_{1}\bar{I}+(\beta_{1})^{2}\bar{I}^{2}}}}\right]\,.$
(40)

To gain insights into the influence of $\bar{I}$ in (40), one can employ a
Taylor series expansion centered around $\bar{I}=0$. However, more interesting
conclusions can be drawn by directly examining (40) as $\bar{I}$ approaches
zero. Notably, the behavior of the eigenvalues in both branches exhibits
distinct characteristics.

Regarding the positive eigenvalue branch, as $\bar{I}\to 0$ (or equivalently,
$\bar{t}\to 0$), all terms in the numerators and denominators of (40) tend
towards a finite value. Consequently, the limit yields the relative error.

$\lim_{h\to
0}\frac{\lambda^{h}_{n+}}{\lambda_{n+}}-1=\frac{\beta_{0}^{h}}{\beta_{0}}\frac{1}{2}\left[1+\sqrt{1-\frac{\tilde{a}_{0}^{h}}{(\beta_{0}^{h})^{2}}}\,\right]-1\,.$
(41)

Next, in the limit of refinement ($h\to 0$) we obtain
$\tilde{a}_{0}^{h}/(\beta_{0}^{h})^{2}\to 0$, and
$\beta_{0}^{h}/\beta_{0}\to\frac{1}{(\xi\pi)^{2}}\frac{\tilde{\mathsf{K}}_{m}^{11}}{\tilde{\mathsf{M}}}$.
Interestingly, this means we recover (33) in the limit of $\bar{I}\to 0$,
which was an expression associated to the negative eigenvalue branch. Thus,
when taking limits as $h\to 0$ and $\bar{I}\to 0$, the order of the limits
matter.

The more interesting case arises when considering the negative eigenvalue
branch and examining the limit as $\bar{I}\to 0$. In this case, the
denominator involving the square root approaches zero, while the numerator
converges to the finite value of
$\frac{\tilde{a}_{0}^{h}}{2(\beta_{0}^{h})^{2}}$. Consequently, the relative
error becomes unbounded! Indeed, a proper Laurent series expansion reveals a
leading term of
$(\frac{\beta_{0}}{\beta_{0}^{h}})(\frac{\tilde{a}_{0}^{h}}{\tilde{a}_{1}})(\frac{1}{\bar{I}})$.
This indicates that for small values of $\bar{I}$, the error can become
exceedingly large. However, it is important to note that this effect
diminishes as the mesh is refined since it can be observed that
$\tilde{a}_{0}^{h}/\tilde{a}_{1}$ rapidly approaches zero as $h\to 0$. In
other words, the “blowup” resulting from small $\bar{I}$ or $\bar{t}$
disappears as the mesh is refined. Nevertheless, it is highly likely that this
phenomenon is responsible for the significant relative error observed in our
results for moderate values of $h$ (as seen in the negative eigenvalue branch
of the standard method in Figure 15). We will discuss this further in the next
section.

##  7 Discussion

In this section we discuss the qualitative and quantitative mathematical
results obtained in the previous sections. Specifically, we explore the
dependence of locking behavior on the physical parameters and provide
insightful discussions on this topic.

###  7.1 Characterization of membrane locking

To initiate our discussion, it is worth noting that membrane locking has
traditionally been described in physical rather than _mathematical_ terms.
However, this does preclude its mathematical description, as we aim to
demonstrate in the subsequent analysis, even if it deviates from other
conventional forms of locking phenomena. We reiterate the key advantages of
our one-dimensional Euler-Bernoulli circular beam model: the deliberate
absence of shear strains, a constant nonzero curvature, and periodic boundary
conditions that eliminate other spurious behaviors due to boundaries. These
unique features positions our model as highly suitable for isolating and
characterizing membrane locking.

To dispense with rigid-body transformations, consider $V$ to be the span of
eigenfunctions not associated with rigid-body motions (i.e., with nonzero
eigenvalues). In technical terms, considering a linear functional
$l:V\rightarrow\mathbb{R}$ with a given forcing within the range of the
operator, the system governed by the standard variational formulation can be
represented as:

$\displaystyle a(\boldsymbol{u},\,\boldsymbol{\delta u})=l(\boldsymbol{\delta
u})\qquad\forall\boldsymbol{\delta u}=(\delta u,\delta w)\in V,$ (42)

where the bilinear operator $a$ is defined as in equation (2b). This system
has a solution (up to rigid-body translations) expressed by:

$\displaystyle\boldsymbol{u}=\sum_{n=1}^{\infty}\alpha_{n}\,\boldsymbol{u}_{n}\qquad\text{with
}\alpha_{n}=l(\boldsymbol{u}_{n})/\lambda_{n}.$ (43)

This solution involves an expansion in terms of the eigenmodes
$\boldsymbol{u}_{n}$, with corresponding eigenvalues $\lambda_{n}\neq 0$. The
coefficients $\alpha_{n}$ are determined based on the eigenvalues and the
linear functional $l$. This expression also holds true in the discrete case,
but the sum is finite up to the total number of modes being discretized.

The total error in approximating $\boldsymbol{u}$ in (43) is a combination of
the error accrued in approximating the eigenmodes $\boldsymbol{u}_{n}$ and its
corresponding eigenvalues $\lambda_{n}$. The approximation error in the modes
is influenced by the accuracy of each component, namely the discrete and
analytical Fourier modes, as well as the approximation error in their relative
amplitudes. It is possible to demonstrate that the discrete Fourier modes
provide optimal approximations of the analytical modes. Additionally, in
Figure 19 we show that the relative amplitudes associated with the membrane-
and bending-dominated responses are well approximated, even when using coarse
meshes. Based on these observations, we argue that the primary source of error
stems from the misapproximation of the eigenvalues.

In Section 5.2 (Figures 15 and 17), we have observed that the discrete
eigenvalues corresponding to the negative branch significantly overestimate
the true eigenvalues, particularly for coarse discretizations. On the other
hand, all other eigenvalues, especially those in the positive branch, are well
approximated. This discrepancy in the negative eigenvalue branch, particularly
for low mode numbers $n$, results in a significant underestimation of the
associated coefficients $\alpha_{n}$ in (43) since they are inversely
proportional to the eigenvalues
($\alpha_{n}=l(\boldsymbol{u}_{n})/\lambda_{n}$). Figure 5 further confirms
that these eigenvalues correspond to a bending-dominated response, indicated
by $U_{n-}/W_{n-}<1$. Consequently, the bending-dominated response is unduly
suppressed in the discrete solution, leading to a relative over-representation
of the membrane-dominated contributions. In essence, the discrete solution
exhibits membrane-locking. In light of the fact that the bending part is
underestimated while the membrane response is accurately represented, one
could argue if “membrane locking” is “the right” term for the phenomenon being
observed.

This discussion serves to reinforce and validate the arguments presented
throughout this article, further justifying our focus on the spectrum and the
proposed criterion developed in Section 4.3. Furthermore, with an explicit
expression for the relative error in the eigenvalues associated with membrane
locking, we are well-equipped to further study this phenomenon. In the
following subsections, we will explore the implications of the observed
eigenvalue errors and investigate their causes.

###  7.2 Membrane locking as a function of mesh size and physical parameters

Upon examining expressions (38)–(40), we can gain insights into the phenomenon
of membrane locking as it relates to the mesh size $h$ and the physical
parameters, which, as we previously discussed, are captured by the normalized
thickness $\bar{t}=\frac{t}{R}$ or equivalently the normalized moment of
inertia $\bar{I}=\frac{\bar{t}^{2}}{12}$.

The occurrence of membrane locking can be attributed to a fundamental
difference among the expressions for the terms $M_{n}$ and $M_{n}^{h}$ in (38)
and (39): $M_{n}$ is a multiple of $\bar{I}=\frac{\bar{t}^{2}}{12}$ so, as a
polynomial in $\bar{I}$, it lacks constant and quadratic terms; however, in
the discrete case, $M_{n}^{h}$ is a polynomial that does include constant and
quadratic terms, specifically denoted by $\tilde{a}_{0}^{h}$ and
$\tilde{a}_{2}^{h}$ in (39). The nonzero values of $\tilde{a}_{0}^{h}$ and
$\tilde{a}_{2}^{h}$ have a profound impact on the behavior of the discrete
solution.

First, the presence of $\tilde{a}_{0}^{h}$ leads to a singularity in the
relative error as $\bar{I}$ approaches zero with the leading relevant
coefficient being $\tilde{a}_{0}^{h}/\tilde{a}_{1}$. As the mesh is refined,
this effect quickly subsides because this coefficient converges rapidly to
zero, specifically at a rate of $\mathcal{O}(h^{4})$. This is shown in Figure
21 where the $n=2$ eigenvalue of the negative branch is seen to be very large
at large values of $h$ and small values of $\bar{I}$. At a fixed value of
$\xi$, Figure 22 also shows how the negative branch of eigenvalues is severely
over-approximated at coarse $h$, but quickly recovers accuracy for small $h$.

On the other hand, the non-zero coefficient $\tilde{a}_{2}^{h}$ introduces a
term that is linear in $\bar{I}$ (quadratic in $\bar{t}$). The relevant
coefficient for this term is $\frac{\tilde{a}_{2}^{h}}{\tilde{a}_{1}}$. Unlike
$\tilde{a}_{0}^{h}$, this coefficient does not rapidly diminish as the mesh is
refined for a fixed value of $\xi$. However, it may decrease with
$p$-refinement, as discussed further below. In fact, according to (33), this
term eventually becomes part of the _membrane_ -dominated response in the
limit of mesh refinement (and this is indeed what happens to the left of the
diagonal dip observed in Figure 22).

Figure 21: Normalized eigenvalue error of the second eigenvalue, $n=2$, of the
negative branch as a function of mesh size $h=\pi/N$ and the normalized moment
of inertia $\bar{I}=\bar{t}^{2}/12$. Note the blowup that occurs on coarse
meshes when $\bar{I}\to 0$. This is a manifestation of membrane locking.
Figure 22: Landscape of the normalized eigenvalues at fixed $\xi=0.1$ as a
function of (continuous) mesh size $h$ and the normalized moment of inertia
$\bar{I}=\bar{t}^{2}/12$. Note the blowup that occurs on coarse meshes, which
gets more severe when $\bar{I}\to 0$. This is a manifestation of membrane
locking. Note the peculiar behavior that runs diagonally across the plot,
which occurs approximately for $h\approx\xi\pi\sqrt{\bar{I}}$, and corresponds
to when this negative eigenvalue branch becomes membrane-dominated (recall (9)
which translates to $\hat{n}(\bar{I})\approx 1/\sqrt{\bar{I}}$ and Section
5.4).

Thus, the non-trivial interplay between $\tilde{a}_{0}^{h}$ and
$\tilde{a}_{2}^{h}$ is crucial in understanding the occurrence of membrane
locking. When dealing with coarse meshes, the magnitude of $\tilde{a}_{0}^{h}$
is significant, resulting in the $\frac{1}{\bar{I}}$ term to dominate the
expression of relative error. As a consequence, the relative error becomes
large as $\bar{I}$ decreases, leading to membrane locking. Conversely, in the
case of fine meshes with a fixed $\xi$ value, $\tilde{a}_{0}^{h}$ approaches
zero, allowing the linear term in $\bar{I}$ to dominate. Consequently, the
relative error becomes large as $\bar{I}$ increases.

Finally, it is worth noting that there are mathematical similarities between
membrane locking and shear locking. Specifically, membrane locking occurs when
$\bar{I}$ is small, corresponding to a small thickness $t$. This resemblance
is not unexpected, as shear locking also introduces an artificial bending
stiffness, particularly when the thickness $t$ approaches zero ([6, Chapter
9.4, Page 368].).

###  7.3 Membrane locking as a function of polynomial degree

As pointed out in Section 6 the dependence on polynomial degree $p$ lies
hidden in the circulant eigenvalues $\tilde{\mathsf{K}}_{b}^{ij}$,
$\tilde{\mathsf{K}}_{m}^{ij}$ and $\tilde{\mathsf{M}}$. This dependence
implies that the asymptotic limits of mesh refinement are influenced by the
choice of $p$. Understanding these limits is crucial, as they not only
determine the convergence rate but also affect the absolute value of the
relative error.

While increasing the value of $p$ significantly improves the convergence rate
and the overall accuracy of the method, it does not eliminate the parasitic
behavior we identify as membrane locking as defined in Section 4.3. This is
because the relative error of the eigenvalues still deviates significantly
from the asymptotic curve, particularly at low values of $\xi$. However, it is
worth mentioning that the terms
$\tilde{\mathsf{K}}_{m}^{11}\tilde{\mathsf{K}}_{m}^{22}-|\tilde{\mathsf{K}}_{m}^{12}|^{2}$
and
$\tilde{\mathsf{K}}_{b}^{11}\tilde{\mathsf{K}}_{b}^{22}-|\tilde{\mathsf{K}}_{b}^{12}|^{2}$
appearing in $\tilde{a}_{0}^{h}$ and $\tilde{a}_{2}^{h}$ from (39), tend to
decrease as the value of $p$ increases. This reduction may potentially help in
mitigating the effects of membrane locking. However, it is important to note
that in the case of cubic and quartic discretizations, we have not observed a
significant reduction in membrane locking despite the presence of these terms.
Membrane locking remains prominent in these cases, as may be observed in
Figures 17 and 18.

Lastly, as discussed in Section 6.3, the preasymptotic domain can be
effectively described by the relative error in the asymptotic limit of
refinement, which is $p$-dependent. For instance, Figure 17 illustrates that
for $p=4$, the onset of “good” convergence initiates at an earlier stage (for
coarser values of $h$ at a fixed $n$). This occurs because the accelerated
decay of the error occurs at a larger value of $\xi$ compared to $p=2$ and
$p=3$.

###  7.4 Comparison of the standard and mixed formulation

The mixed Galerkin discretization based on the Hellinger-Reissner principle
was shown to be largely free of membrane locking (Figure 16). Furthermore, the
asymptotic accuracy of eigenvalues for bending-dominated modes was
significantly better in the mixed formulation compared to the standard
formulation, see Figure 15 and 17. A more detailed study of the relative error
in the case of this variational formulation is viable, but is left for a
future analysis.

##  8 Conclusion

This paper addresses the issue of membrane locking which affects finite
element approximations of thin beams and shells. Membrane locking, along with
transverse shear locking and volumetric locking, negatively impacts the
accuracy and reliability of finite element discretizations. However, membrane
locking is still not well understood, and the absence of a rigorous assessment
methodology presents challenges in the development of effective beam and shell
elements. Our research aims to achieve two main goals: firstly, to gain
insights into the mechanisms underlying membrane locking, and secondly, to
establish a rigorous criterion for accurately identifying and evaluating the
extent of membrane locking.

Our study involved a comprehensive Fourier analysis of the standard and mixed
Galerkin formulations of the Euler-Bernoulli beam model for a circular ring,
employing smooth splines for the discretization. Importantly, this model
selection enabled the effective isolation and investigation of membrane
locking due to several key factors: the model inherently lacks shear strains,
maintains a constant nonzero curvature, and incorporates periodic boundary
conditions to eliminate other undesired phenomena, including e.g. boundary
layers. Furthermore, our choice of model facilitated exact formulas for the
_discrete_ eigenvalues and corresponding eigenfunctions by exploiting
properties of circulant matrices. The exact expressions were presented in
terms of the relevant discretization and physical parameters, and compared to
the exact analytical solutions.

Our analysis revealed two distinct types of eigenmodes: bending-dominated and
membrane-dominated. Remarkably, we observed that all eigenmodes were well
approximated even with coarse to medium discretizations. However, the
eigenvalues associated with bending modes in the standard formulation
exhibited significant deviations, especially in coarse discretizations. This
motivated the development of a criterion to accurately identify and quantify
the occurrence of membrane locking based on accuracy of eigenvalues. This
criterion effectively illuminated the presence and severity of membrane
locking in the standard formulation for quadratics, cubics and quartics.
Despite the substantial improvement in absolute accuracy with increasing
polynomial order, we found that it did not mitigate the occurrence of locking
behavior. In contrast, when utilizing a mixed Galerkin discretization to
numerically solve the same problem, we observed minimal effects of membrane
locking on the solutions across all polynomial degrees and meshes. Moreover,
we discovered that the absolute accuracy of the mixed method was significantly
superior, particularly in relation to the eigenvalues associated with bending-
dominated modes.

Furthermore, we mathematically derived and analyzed exact expressions for the
relative error. We established the expected convergence rates for a fixed mode
number and determined that the normalized spectra in the asymptotic refinement
limit exhibited a relationship with the spectra associated with isolated
bending and membrane responses. Significantly, our investigation revealed that
in coarse discretizations of the standard formulation, the eigenvalues
associated with the bending-dominated response were grossly overestimated,
particularly for a small thickness and large radius of curvature. The result
is a suppression of the bending-dominated component of the discrete solution.
Hence, our careful eigenvalue analysis demonstrates that the phenomenon
observed is not simply an overestimation of the membrane response but rather a
suppression of bending, prompting a question regarding the suitability of the
term “membrane locking” for this phenomenon. We hypothesize that this effect
shares mathematical similarities with shear locking, despite our model not
incorporating shear strains.

Looking ahead, we plan to investigate membrane locking in doubly curved shells
and explore the extension of our analytical techniques to periodic, uniform
spline discretizations of a torus. Additionally, we would like to explore
other types of shell elements, such as non-conforming techniques based on
$C^{0}$ finite elements or DG. Despite the resulting submatrices not being
circulant, they exhibit a circulant block structure with block size
proportional to the polynomial order, which may enable analytical solution
techniques. Furthermore, we note that our techniques may provide insights into
other forms of locking, specifically shear locking in thin beams.

## Acknowledgments

R.R. Hiemstra and D. Schillinger gratefully acknowledge financial support from
the German Research Foundation (Deutsche Forschungsgemeinschaft) through the
DFG Emmy Noether Grant SCH 1249/2-1 and the standard DFG grant SCH 1249/5-1.
F. Fuentes’ work was partially supported by the National Center for Artificial
Intelligence CENIA FB210017, Basal ANID. The authors thank T.J.R Hughes and
T.-H. Nguyen for helpful discussions.

##  Appendix A Eigenvalue decomposition of matrices with circulant blocks

The discrete eigenvalue problems that arise in this paper involve two-by-two
block matrices with circulant blocks, sometimes referred to as circulant block
matrices (as opposed to block-circulant matrices). First, we describe the
special structure of circulant matrices and their eigenvalue decomposition
based on the Discrete Fourier Transform (DFT). Then we use these properties to
diagonalize two-by-two block matrices, which are the focus of this paper.

###  A.1 Eigenvalue decomposition of a circulant matrix

A circulant $N\times N$ real matrix is uniquely defined by $N$ real numbers,
$c_{0},c_{1},\ldots,c_{N-1}$, and has the following special structure

$\displaystyle\mathsf{C}=\begin{bmatrix}c_{0}&c_{1}&\cdots&c_{N-2}&c_{N-1}\\\
c_{N-1}&c_{0}&c_{1}&\cdots&c_{N-2}\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\
c_{2}&\cdots&c_{N-1}&c_{0}&c_{1}\\\
c_{1}&c_{2}&\cdots&c_{N-1}&c_{0}\end{bmatrix}$ (44)

Every circulant matrix has an eigenvalue decomposition

$\displaystyle\mathsf{C}=\mathsf{F}\mathsf{\Lambda}\mathsf{F}^{-1}$ (45)

where $\mathsf{F}\in\mathbb{C}^{N\times N}$, called the Fourier matrix,
denotes the collection of eigenvectors and $\mathsf{\Lambda}$ the
corresponding eigenvalues. The eigenvectors, the columns of the Fourier
matrix, are always the same, namely the Fourier modes,

$\displaystyle\mathsf{v}_{j}=\frac{1}{\sqrt{N}}\begin{pmatrix}1&\omega^{j}&\omega^{2j}&\ldots&\omega^{(N-1)j}\end{pmatrix}^{T},\;\;j=0,1,\ldots,N-1,$
with $\displaystyle\omega=\exp{(2\pi\cdot i/N)}$ (46)

The corresponding eigenvalues are defined via the following simple
relationship

$\displaystyle\lambda_{j}=c_{0}+c_{1}\cdot\omega^{j}+c_{2}\cdot\omega^{2j}+\cdots+c_{N-1}\cdot\omega^{(N-1)j},\qquad
j=0,1,\ldots,N-1.$ (47)

If $\mathsf{C}$ is symmetric then the eigenvalues are all real and are given
by

$\displaystyle\lambda_{j}=c_{0}+2c_{1}\cdot\cos(2\pi\cdot
j/N)+2c_{2}\cdot\cos(4\pi\cdot j/N)+\ldots+2c_{\lfloor
N/2\rfloor}\cdot\cos(\lfloor N/2\rfloor\pi\cdot j/N).$ (48)

On the other hand, if $\mathsf{C}$ is skew-symmetric then the eigenvalues are
complex and are given by

$\displaystyle\lambda_{j}=2ic_{1}\cdot\sin(2\pi\cdot
j/N)+2ic_{2}\cdot\sin(4\pi\cdot j/N)+\ldots+2ic_{\lfloor
N/2\rfloor}\cdot\sin(\lfloor N/2\rfloor\pi\cdot j/N).$ (49)

###  A.2 Similarity transformations of matrices with circulant blocks

Let $\mathsf{A},\mathsf{B},\mathsf{C},\mathsf{D}\in\mathbb{R}^{N\times N}$ be
circulant matrices. We are interested in the eigenvalue decomposition of the
$2\times 2$ block matrix

$\displaystyle\mathsf{K}=\begin{bmatrix}\mathsf{A}&\mathsf{B}\\\
\mathsf{C}&\mathsf{D}\end{bmatrix}$ (50)

Using the Fourier matrix $\mathsf{F}$ we can bring matrix $\mathsf{K}$ to the
following simpler form

$\displaystyle\tilde{\mathsf{K}}=\begin{bmatrix}\tilde{\mathsf{A}}&\tilde{\mathsf{B}}\\\
\tilde{\mathsf{C}}&\tilde{\mathsf{D}}\end{bmatrix}=\begin{bmatrix}\mathsf{F}^{-1}&\\\
&\mathsf{F}^{-1}\end{bmatrix}\begin{bmatrix}\mathsf{A}&\mathsf{B}\\\
\mathsf{C}&\mathsf{D}\end{bmatrix}\begin{bmatrix}\mathsf{F}&\\\
&\mathsf{F}\end{bmatrix}$ (51)

where $\tilde{\mathsf{A}},\tilde{\mathsf{B}},\tilde{\mathsf{C}}$ and
$\tilde{\mathsf{D}}$ are diagonal with entries the eigenvalues of the
respective matrices. This matrix is easier to diagonalize than the original
matrix $\mathsf{K}$.

###  A.3 Eigenvalue decomposition of a $2\times 2$ block matrix with diagonal
blocks

Let $\tilde{\mathsf{K}}$ denote a $2\times 2$ block matrix with diagonal
blocks. Since diagonal matrices commute,
$\tilde{\mathsf{C}}\tilde{\mathsf{D}}=\tilde{\mathsf{D}}\tilde{\mathsf{C}}$,
we have that the characteristic equation is similar to that of a $2\times 2$
matrix (see [23]),

$\displaystyle\det{\left(\tilde{\mathsf{K}}-\lambda\mathsf{I}\right)}$
$\displaystyle=\det\left(\begin{bmatrix}\tilde{\mathsf{A}}-\lambda_{n}\mathsf{I}&\tilde{\mathsf{B}}\\\
\tilde{\mathsf{C}}&\tilde{\mathsf{D}}-\lambda_{n}\mathsf{I}\end{bmatrix}\right)$
(52)
$\displaystyle=\det{\left(\left(\tilde{\mathsf{A}}-\lambda_{n}\mathsf{I}\right)\left(\tilde{\mathsf{D}}-\lambda_{n}\mathsf{I}\right)-\tilde{\mathsf{B}}\tilde{\mathsf{C}}\right)}$
(53)
$\displaystyle=\det{\left(\tilde{\mathsf{A}}\tilde{\mathsf{D}}-\tilde{\mathsf{B}}\tilde{\mathsf{C}}-\lambda_{n}\left(\tilde{\mathsf{A}}+\tilde{\mathsf{D}}\right)+\lambda_{n}^{2}\mathsf{I}\right)}$
(54)
$\displaystyle=\prod_{n=1}^{N}\left(\tilde{A}_{n}\tilde{D}_{n}-\tilde{B}_{n}\tilde{C}_{n}-\lambda_{n}\left(\tilde{A}_{n}+\tilde{D}_{n}\right)+\lambda_{n}^{2}\right)$
(55)

The characteristic equation of the $2N\times 2N$ system is simply the product
of $N$ characteristic equations of $2\times 2$ matrices. Hence, instead of the
$2N\times 2N$ matrix eigenvalue problem, we may consider, for each $n$, the
$2\times 2$ matrix eigenvalue problem

$\displaystyle\begin{bmatrix}\tilde{A}_{n}-\lambda_{n}&\tilde{B}_{n}\\\
\tilde{C}_{n}&\tilde{D}_{n}-\lambda_{n}\end{bmatrix}\begin{bmatrix}U_{n}\\\
W_{n}\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix}.$ (56)

The eigenvalues follow as the roots of the characteristic equation

$\displaystyle\lambda_{n\pm}$ $\displaystyle=K_{n}\left(1\pm L_{n}\right)$
(57a) where $K_{n}$ and $L_{n}$ are constants that depend on $n$
$\displaystyle K_{n}=\frac{\tilde{A}_{n}+\tilde{D}_{n}}{2},\qquad
M_{n}=\tilde{A}_{n}\cdot\tilde{D}_{n}-\tilde{C}_{n}\cdot\tilde{B}_{n},\qquad
L_{n}=\sqrt{1-M_{n}/K^{2}_{n}}.$ (57b)

The eigenfunctions are defined up to a constant. However, the ratio of the
amplitudes the eigenfunctions, $U_{n}/W_{n}$, satisfy

$\displaystyle\rho_{n\pm}=\frac{U_{n\pm}}{W_{n\pm}}=\frac{\tilde{B}_{n}}{\lambda_{n\pm}-\tilde{A}_{n}}$
(58)

##  Appendix B Matrix eigenvalue problem

We compute exact expressions for the entries of the block-circulant stiffness
and mass matrices that arise in the standard and mixed formulation of the
discrete eigenvalue problem based on approximation with uniform spline basis
functions of degree $p=2,3,4$. The expressions are given in terms of rational
numbers and symbolic variables, such as the mesh size $h$, the radius $R$ and
the normalized beam thickness $\bar{t}=\tfrac{t}{R}$.

###  B.1 Standard formulation

The matrix eigenvalue problem is

$\displaystyle\mathbf{K}\,\mathbf{u}=\lambda^{h}\mathbf{M}\,\mathbf{u}.$ (59)

The mass and stiffness matrices are $2\times 2$ block-circulant matrices that
have the structure

$\displaystyle\mathbf{K}=\mathbf{K}_{m}+\mathbf{K}_{b}=\begin{bmatrix}\mathsf{K}^{11}_{m}&\mathsf{K}^{12}_{m}\\\
\mathsf{K}^{21}_{m}&\mathsf{K}^{22}_{m}\end{bmatrix}+\begin{bmatrix}\mathsf{K}^{11}_{b}&\mathsf{K}^{12}_{b}\\\
\mathsf{K}^{21}_{b}&\mathsf{K}^{22}_{b}\end{bmatrix},$
$\displaystyle\mathbf{M}=\begin{bmatrix}\mathsf{M}&\\\
&\mathsf{M}\end{bmatrix}.$ (60)

Here $\mathbf{K}_{m}$ denote the contribution to the membrane stiffness and
$\mathbf{K}_{b}$ the portion related to the bending. The blocks $\mathsf{M}$
are circulant, symmetric positive definite matrices, and are computed as

$\displaystyle[\mathsf{M}]_{ij}=R\int_{0}^{2\pi}B_{i,p}(\theta)B_{j,p}(\theta)\,\mathrm{d}\theta.$
(61)

In the stiffness matrix we make a distinction between the action of membrane
and bending terms, respectively. The submatrices on the diagonal are circulant
and symmetric. The off-diagonal submatrices are skew-symmetric circulant
matrices. These are computed according to the relations

ij
$\displaystyle=\frac{1}{R}\int_{0}^{2\pi}B^{\prime}_{i,p}(\theta)B^{\prime}_{j,p}(\theta)\,\mathrm{d}\theta,\qquad\quad$
$\displaystyle[\mathsf{K}^{11}_{b}]_{ij}$
$\displaystyle=\frac{1}{R}\frac{\bar{t}^{2}}{12}\int_{0}^{2\pi}B^{\prime}_{i,p}(\theta)B^{\prime}_{j,p}(\theta)\,\mathrm{d}\theta,$
(62) $\displaystyle[\mathsf{K}^{12}_{m}]_{ij}$
$\displaystyle=\frac{1}{R}\int_{0}^{2\pi}B^{\prime}_{i,p}(\theta)B_{j,p}(\theta)\,\mathrm{d}\theta,\qquad\quad$
$\displaystyle[\mathsf{K}^{12}_{b}]_{ij}$
$\displaystyle=\frac{1}{R}\frac{\bar{t}^{2}}{12}\int_{0}^{2\pi}-B^{\prime}_{i,p}(\theta)B^{\prime\prime}_{j,p}(\theta)\,\mathrm{d}\theta$
$\displaystyle[\mathsf{K}^{22}_{m}]_{ij}$
$\displaystyle=\frac{1}{R}\int_{0}^{2\pi}B_{i,p}(\theta)B_{j,p}(\theta)\,\mathrm{d}\theta,\qquad\quad$
$\displaystyle[\mathsf{K}^{22}_{b}]_{ij}$
$\displaystyle=\frac{1}{R}\frac{\bar{t}^{2}}{12}\int_{0}^{2\pi}B^{\prime\prime}_{i,p}(\theta)B^{\prime\prime}_{j,p}(\theta)\,\mathrm{d}\theta$

Table 1 lists the non-zero values of the circulant blocks of the mass and
stiffness matrix in the case of $p=2,3$ and $4$. Inverses, transposes and
products of circulant matrices are circulant. Hence, the blocks in the matrix
$\tilde{\mathbf{K}}=\mathbf{M}^{-1}\mathbf{K}$ remain circulant and the
technique in A can be used to compute the discrete eigenvalues and
eigenvectors of the system.

Table 1: Circulant submatrices of the mass and stiffness matrix of the standard formulation $\frac{1}{R}\frac{1}{h}\mathsf{M}=R\frac{1}{h}\mathsf{K}^{22}_{m}$ |  $c_{0}$ |  $c_{1}=c_{N-1}$ |  $c_{2}=c_{N-2}$ |  $c_{3}=c_{N-3}$ |  $c_{4}=c_{N-4}$
---|---|---|---|---|---
$p=2$ |  $\frac{11}{20}$ |  $\frac{13}{60}$ |  $\frac{1}{120}$ |  |
$p=3$ |  $\frac{151}{315}$ |  $\frac{397}{1680}$ |  $\frac{1}{42}$ |  $\frac{1}{5040}$ |
$p=4$ |  $\frac{15619}{36288}$ |  $\frac{44117}{181440}$ |  $\frac{913}{22680}$ |  $\frac{251}{181440}$ |  $\frac{1}{362880}$
$Rh\mathsf{K}^{11}_{m}=Rh\frac{12}{\bar{t}^{2}}\mathsf{K}^{11}_{b}$ |  $c_{0}$ |  $c_{1}=c_{N-1}$ |  $c_{2}=c_{N-2}$ |  $c_{3}=c_{N-3}$ |  $c_{4}=c_{N-4}$
$p=2$ |  $1$ |  $-\frac{1}{3}$ |  $-\frac{1}{6}$ |  |
$p=3$ |  $\frac{2}{3}$ |  $-\frac{1}{8}$ |  $-\frac{1}{5}$ |  $-\frac{1}{120}$ |
$p=4$ |  $\frac{35}{72}$ |  $-\frac{11}{360}$ |  $-\frac{17}{90}$ |  $-\frac{59}{2520}$ |  $-\frac{1}{5040}$
$R\mathsf{K}^{12}_{m}$ |  $c_{0}$ |  $c_{1}=-c_{N-1}$ |  $c_{2}=-c_{N-2}$ |  $c_{3}=-c_{N-3}$ |  $c_{4}=-c_{N-4}$
$p=2$ |  $0$ |  $\frac{5}{12}$ |  $\frac{1}{24}$ |  |
$p=3$ |  $0$ |  $\frac{49}{144}$ |  $\frac{7}{90}$ |  $\frac{1}{720}$ |
$p=4$ |  $0$ |  $\frac{809}{2880}$ |  $\frac{289}{2880}$ |  $\frac{41}{6720}$ |  $\frac{1}{40320}$
$Rh^{2}\frac{12}{\bar{t}^{2}}\mathsf{K}^{12}_{b}$ |  $c_{0}$ |  $c_{1}=-c_{N-1}$ |  $c_{2}=-c_{N-2}$ |  $c_{3}=-c_{N-3}$ |  $c_{4}=-c_{N-4}$
$p=2$ |  $0$ |  $1$ |  $-\frac{1}{2}$ |  |
$p=3$ |  $0$ |  $\frac{19}{24}$ |  $-\frac{1}{3}$ |  $-\frac{1}{24}$ |
$p=4$ |  $0$ |  $\frac{217}{360}$ |  $-\frac{67}{360}$ |  $-\frac{3}{40}$ |  $-\frac{1}{720}$
$Rh^{3}\frac{12}{\bar{t}^{2}}\mathsf{K}^{22}_{b}$ |  $c_{0}$ |  $c_{1}=c_{N-1}$ |  $c_{2}=c_{N-2}$ |  $c_{3}=c_{N-3}$ |  $c_{4}=c_{N-4}$
$p=2$ |  $6$ |  $-4$ |  $1$ |  |
$p=3$ |  $\frac{8}{3}$ |  $-\frac{3}{2}$ |  $0$ |  $\frac{1}{6}$ |
$p=4$ |  $\frac{19}{12}$ |  $-\frac{43}{60}$ |  $-\frac{4}{15}$ |  $\frac{11}{60}$ |  $\frac{1}{120}$

###  B.2 Mixed formulation

The generalized matrix eigenvalue problem is

$\displaystyle\begin{bmatrix}-\mathbf{K}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}&\mathbf{K}_{\mathbf{\varepsilon}\mathbf{u}}\\\
\mathbf{K}_{\mathbf{\varepsilon}\mathbf{u}}^{T}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{\varepsilon}\\\
\mathbf{u}\end{bmatrix}=\begin{bmatrix}\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\lambda^{h}\mathbf{M}\end{bmatrix}\begin{bmatrix}\mathbf{\varepsilon}\\\
\mathbf{u}\end{bmatrix}.$ (63)

The matrices have the following structure

$\displaystyle\mathbf{K}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}$
$\displaystyle=\begin{bmatrix}\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}&0\\\
0&\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}\end{bmatrix},$
$\displaystyle\mathbf{K}_{\mathbf{\varepsilon}\mathbf{u}}$
$\displaystyle=\begin{bmatrix}\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{u}}&\mathsf{K}^{12}_{\mathbf{\varepsilon}\mathbf{u}}\\\
\mathsf{K}^{21}_{\mathbf{\varepsilon}\mathbf{u}}&\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{u}}\end{bmatrix},$
$\displaystyle\mathbf{M}=\begin{bmatrix}\mathsf{M}&0\\\
0&\mathsf{M}\end{bmatrix}.$ (64)

The non-zero blocks of the mass and stiffness matrix are all circulant. The
entries of $\mathsf{M}$ are computed as

$\displaystyle[\mathsf{M}]_{ij}=R\int_{0}^{2\pi}B_{i,p}(\theta)B_{j,p}(\theta)\,\mathrm{d}\theta,$
(65)

and the entries of the submatrices of
$\mathbf{K}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}$ and
$\mathbf{K}_{\mathbf{\varepsilon}\mathbf{u}}$ are given by

$\displaystyle[\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}]_{ij}$
$\displaystyle=R\int_{0}^{2\pi}B_{i,p-1}(\theta)B_{j,p-1}(\theta)\,\mathrm{d}\theta,$
$\displaystyle[\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}]_{ij}$
$\displaystyle=R^{3}\frac{\bar{t}^{2}}{12}\int_{0}^{2\pi}B_{i,p-1}(\theta)B_{j,p-1}(\theta)\,\mathrm{d}\theta$
(66) $\displaystyle[\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{u}}]_{ij}$
$\displaystyle=\int_{0}^{2\pi}B_{i,p-1}(\theta)B^{\prime}_{j,p}(\theta)\,\mathrm{d}\theta,$
$\displaystyle[\mathsf{K}^{12}_{\mathbf{\varepsilon}\mathbf{u}}]_{ij}$
$\displaystyle=\int_{0}^{2\pi}B_{i,p-1}(\theta)B_{j,p}(\theta)\,\mathrm{d}\theta$
(67) $\displaystyle[\mathsf{K}^{21}_{\mathbf{\varepsilon}\mathbf{u}}]_{ij}$
$\displaystyle=R\frac{\bar{t}^{2}}{12}\int_{0}^{2\pi}B_{i,p-1}(\theta)B^{\prime}_{j,p}(\theta)\,\mathrm{d}\theta,$
$\displaystyle[\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{u}}]_{ij}$
$\displaystyle=R\frac{\bar{t}^{2}}{12}\int_{0}^{2\pi}-B_{i,p-1}(\theta)B^{\prime\prime}_{j,p}(\theta)\,\mathrm{d}\theta$
(68)

The above generalized eigenvalue problem can be transformed to standard form
using static condensation, that is, by solving for $\mathbf{\varepsilon}$ in
terms of $\mathbf{u}$, leading to

$\displaystyle\mathbf{K}\,\mathbf{u}=\lambda^{h}\mathbf{M}\mathbf{u}$ with
$\displaystyle\mathbf{K}=\mathbf{K}_{\mathbf{\varepsilon}\mathbf{u}}^{T}\mathbf{K}^{-1}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}\,\mathbf{K}_{\mathbf{\varepsilon}\mathbf{u}}.$
(69)

As before, the contribution to the stiffness matrix can be factored into a
contribution related to membrane strains and a contribution related to bending

$\displaystyle\mathbf{K}=\mathbf{K}_{m}+\mathbf{K}_{b}=\begin{bmatrix}\mathsf{K}^{11}_{m}&\mathsf{K}^{12}_{m}\\\
\mathsf{K}^{21}_{m}&\mathsf{K}^{22}_{m}\end{bmatrix}+\begin{bmatrix}\mathsf{K}^{11}_{b}&\mathsf{K}^{12}_{b}\\\
\mathsf{K}^{21}_{b}&\mathsf{K}^{22}_{b}\end{bmatrix}.$ (70)

where,

$\displaystyle\mathsf{K}^{11}_{m}$
$\displaystyle=(\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{u}})^{T}\,(\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{\varepsilon}})^{-1}\,\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{u}},$
$\displaystyle\mathsf{K}^{11}_{b}$
$\displaystyle=(\mathsf{K}^{21}_{\mathbf{\varepsilon}\mathbf{u}})^{T}\,(\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{\varepsilon}})^{-1}\,\mathsf{K}^{21}_{\mathbf{\varepsilon}\mathbf{u}}$
(71) $\displaystyle\mathsf{K}^{12}_{m}$
$\displaystyle=(\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{u}})^{T}\,(\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{\varepsilon}})^{-1}\,\mathsf{K}^{12}_{\mathbf{\varepsilon}\mathbf{u}},$
$\displaystyle\mathsf{K}^{12}_{b}$
$\displaystyle=(\mathsf{K}^{21}_{\mathbf{\varepsilon}\mathbf{u}})^{T}\,(\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{\varepsilon}})^{-1}\,\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{u}}$
(72) $\displaystyle\mathsf{K}^{22}_{m}$
$\displaystyle=(\mathsf{K}^{12}_{\mathbf{\varepsilon}\mathbf{u}})^{T}\,(\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{\varepsilon}})^{-1}\,\mathsf{K}^{12}_{\mathbf{\varepsilon}\mathbf{u}},$
$\displaystyle\mathsf{K}^{22}_{b}$
$\displaystyle=(\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{u}})^{T}\,(\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{\varepsilon}})^{-1}\,\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{u}}.$
(73)

Table 2 lists the non-zero values of the circulant blocks of the mass and
stiffness matrix in the case of $p=2,3$ and $4$. Inverses, transposes and
products of circulant matrices are circulant. Hence, the blocks in the matrix
$\tilde{\mathbf{K}}=\mathbf{M}^{-1}\mathbf{K}$ remain circulant and the
technique in A can be used to compute the discrete eigenvalues and
eigenvectors of the system.

Table 2: Circulant submatrices of the mass and stiffness matrix of the mixed formulation $\frac{1}{R}\frac{1}{h}\mathsf{M}$ |  $c_{0}$ |  $c_{1}=c_{N-1}$ |  $c_{2}=c_{N-2}$ |  $c_{3}=c_{N-3}$ |  $c_{4}=c_{N-4}$
---|---|---|---|---|---
$p=2$ |  $\frac{11}{20}$ |  $\frac{13}{60}$ |  $\frac{1}{120}$ |  |
$p=3$ |  $\frac{151}{315}$ |  $\frac{397}{1680}$ |  $\frac{1}{42}$ |  $\frac{1}{5040}$ |
$p=4$ |  $\frac{15619}{36288}$ |  $\frac{44117}{181440}$ |  $\frac{913}{22680}$ |  $\frac{251}{181440}$ |  $\frac{1}{362880}$
$\frac{1}{R}\frac{1}{h}\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}=\frac{1}{R^{3}}\frac{1}{h}\frac{12}{\bar{t}^{2}}\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{\varepsilon}}$ |  $c_{0}$ |  $c_{1}=c_{N-1}$ |  $c_{2}=c_{N-2}$ |  $c_{3}=c_{N-3}$ |  $c_{4}=c_{N-4}$
$p=2$ |  $\frac{2}{3}$ |  $\frac{1}{6}$ |  |  |
$p=3$ |  $\frac{11}{20}$ |  $\frac{13}{60}$ |  $\frac{1}{120}$ |  |
$p=4$ |  $\frac{151}{315}$ |  $\frac{397}{1680}$ |  $\frac{1}{42}$ |  $\frac{1}{5040}$ |
$\mathsf{K}^{11}_{\mathbf{\varepsilon}\mathbf{u}}=\frac{1}{R}\frac{12}{\bar{t}^{2}}\,\mathsf{K}^{21}_{\mathbf{\varepsilon}\mathbf{u}}$ |  $c_{0}=-c_{N-1}$ |  $c_{1}=-c_{N-2}$ |  $c_{2}=-c_{N-3}$ |  $c_{3}=-c_{N-4}$ |  $c_{4}=-c_{N-5}$
$p=2$ |  $-\frac{1}{2}$ |  $-\frac{1}{6}$ |  |  |
$p=3$ |  $-\frac{1}{3}$ |  $-\frac{5}{24}$ |  $-\frac{1}{120}$ |  |
$p=4$ |  $-\frac{35}{144}$ |  $-\frac{17}{80}$ |  $-\frac{17}{720}$ |  $-\frac{1}{5040}$ |
$\frac{1}{h}\mathsf{K}^{12}_{\mathbf{\varepsilon}\mathbf{u}}$ |  $c_{0}=c_{N-1}$ |  $c_{1}=c_{N-2}$ |  $c_{2}=c_{N-3}$ |  $c_{3}=c_{N-4}$ |  $c_{4}=c_{N-5}$
$p=2$ |  $\frac{11}{24}$ |  $\frac{1}{24}$ |  |  |
$p=3$ |  $\frac{151}{360}$ |  $\frac{19}{240}$ |  $\frac{1}{720}$ |  |
$p=4$ |  $\frac{15619}{40320}$ |  $\frac{477}{4480}$ |  $\frac{247}{40320}$ |  $\frac{1}{40320}$ |
$\frac{1}{R}h\frac{12}{\bar{t}^{2}}\mathsf{K}^{22}_{\mathbf{\varepsilon}\mathbf{u}}$ |  $c_{0}=c_{N-1}$ |  $c_{1}=c_{N-2}$ |  $c_{2}=c_{N-3}$ |  $c_{3}=c_{N-4}$ |  $c_{4}=c_{N-5}$
$p=2$ |  $\frac{1}{2}$ |  $-\frac{1}{2}$ |  |  |
$p=3$ |  $\frac{5}{12}$ |  $-\frac{3}{8}$ |  $-\frac{1}{24}$ |  |
$p=4$ |  $\frac{49}{144}$ |  $-\frac{21}{80}$ |  $-\frac{11}{144}$ |  $-\frac{1}{720}$ |

##  Appendix C Analytical eigenvalues and eigenfunctions

We represent the circumferential and radial displacement, $u$ and $w$
respectively, by the following Fourier series [20, Chapter 5, Section 5.3,
page 82]

$\displaystyle u(\theta)$ $\displaystyle=\sum_{n=1}^{N}U_{n}\sin(n\theta),$
(74) $\displaystyle w(\theta)$
$\displaystyle=\sum_{n=1}^{N}W_{n}\cos(n\theta).$ (75)

The resulting stiffness and mass matrix have the structure

$\displaystyle\mathbf{K}=\begin{bmatrix}\mathsf{A}&\mathsf{B}\\\
\mathsf{C}&\mathsf{D}\end{bmatrix},$
$\displaystyle\mathbf{M}=\begin{bmatrix}\mathsf{M}&\\\
&\mathsf{M}\end{bmatrix}.$ (76)

Due to the orthogonality of the Fourier basis the submatrices of $\mathbf{K}$
and $\mathbf{M}$ are diagonal, with the diagonal entries of the mass matrix
being $\pi R$. The eigenvalues and eigenmodes of the system
$\mathbf{K}\mathbf{u}=\lambda\mathbf{M}\mathbf{u}$ are computed by performing
an eigenvalue decomposition of the $2\times 2$ block matrix

$\displaystyle\tilde{\mathbf{K}}=\mathbf{M}^{-1}\mathbf{K}=\begin{bmatrix}\tilde{\mathsf{A}}&\tilde{\mathsf{B}}\\\
\tilde{\mathsf{C}}&\tilde{\mathsf{D}}\end{bmatrix}.$ (77)

where each of the blocks remains diagonal by virtue of $\mathbf{M}$ being
diagonal. The values at the diagonals are

$\displaystyle\tilde{A}_{n}=\frac{n^{2}}{R^{2}}\left(1+\frac{t^{2}}{12R^{2}}\right),$
(78a)
$\displaystyle\tilde{B}_{n}=\tilde{C}_{n}=\frac{n}{R^{2}}\left(1+\frac{t^{2}n^{2}}{12R^{2}}\right),$
(78b)
$\displaystyle\tilde{D}_{n}=\frac{1}{R^{2}}\left(1+\frac{t^{2}n^{4}}{12R^{2}}\right).$
(78c)

The eigenvalues and amplitudes of the Fourier modes may now be determined
using the technique in A.3.

## References

  * [1] M. Bischoff, W. A. Wall, K.-U. Bletzinger, E. Ramm, Models and Finite Elements for Thin-walled Structures, in: E. Stein, R. de Borst, T. J. R. Hughes (Eds.), Encyclopedia of Computational Mechanics, Vol. 2, John Wiley & Sons, 2004, Ch. 3, pp. 59–137.
  * [2] H. Stolarski, T. Belytschko, Membrane locking and reduced integration for curved elements, Journal of Applied Mechanics 49 (1982) 172.
  * [3] M. Suri, I. Babuška, C. Schwab, Locking effects in the finite element approximation of plate models, Mathematics of Computation 64 (210) (1995) 461–482.
  * [4] I. Babuška, M. Suri, Locking effects in the finite element approximation of elasticity problems, Numerische Mathematik 62 (1) (1992) 439–463.
  * [5] T. J. R. Hughes, The finite element method: linear static and dynamic finite element analysis, Courier Corporation, 2012.
  * [6] P. Wriggers, Nonlinear finite element methods, Springer, 2008.
  * [7] D. Braess, Finite elements: Theory, fast solvers, and applications in solid mechanics, Cambridge University Press, 2007.
  * [8] A. Ambroziak, Locking effects in the finite element method, Shell structures: Theory and application (2013) 369.
  * [9] T.-H. Nguyen, R. R. Hiemstra, D. Schillinger, Leveraging spectral analysis to elucidate membrane locking and unlocking in isogeometric finite element formulations of the curved Euler–Bernoulli beam, Computer Methods in Applied Mechanics and Engineering 388 (2022) 114240.
  * [10] D. N. Arnold, Discretization by finite elements of a model parameter dependent problem, Numerische Mathematik 37 (1981) 405–421.
  * [11] I. Babuška, M. Suri, On locking and robustness in the finite element method, SIAM Journal on Numerical Analysis 29 (5) (1992) 1261–1293.
  * [12] T. Belytschko, H. Stolarski, W. K. Liu, N. Carpenter, J. S. Ong, Stress projection for membrane and shear locking in shell finite elements, Computer Methods in Applied Mechanics and Engineering 51 (1-3) (1985) 221–258.
  * [13] R. H. Macneal, R. L. Harder, A proposed standard set of problems to test finite element accuracy, Finite Elements in Analysis and Design 1 (1) (1985) 3–20.
  * [14] T. J. R. Hughes, J. A. Evans, A. Reali, Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems, Computer Methods in Applied Mechanics and Engineering 272 (2014) 290–320.
  * [15] G. Strang, G. Fix, An Analysis of the Finite Element Method, 2nd Edition, Wellesley-Cambridge Press, 2008.
  * [16] Y. Lee, H.-M. Jeon, P.-S. Lee, K.-J. Bathe, The modal behavior of the mitc3+ triangular shell element, Computers & Structures 153 (2015) 148–164.
  * [17] Z. Zou, T. J. R. Hughes, M. A. Scott, R. A. Sauer, E. J. Savitha, Galerkin formulations of isogeometric shell analysis: Alleviating locking with greville quadratures and higher-order elements, Computer Methods in Applied Mechanics and Engineering 380 (2021) 113757.
  * [18] D. Braess, A multigrid method for the membrane problem, Computational Mechanics 3 (5) (1988) 321–329.
  * [19] K. C. Park, D. L. Flaggs, A Fourier analysis of spurious mechanisms and locking in the finite element method, Computer Methods in Applied Mechanics and Engineering 46 (1) (1984) 65–81.
  * [20] W. Soedel, Vibrations of Shells and Plates, Marcel Dekker, Inc., New York, 2004\.
  * [21] C. de Boor, A practical guide to splines, Springer, 1978.
  * [22] R. J. LeVeque, Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems, SIAM, 2007.
  * [23] J. R. Silvester, Determinants of block matrices, The Mathematical Gazette 84 (501) (2000) 460–467.