A Thorough Look at the (In)Stability of Piston-theoretic Beams

We consider a beam model representing the transverse deflections of a one dimensional elastic structure immersed in an axial fluid flow. The model includes a nonlinear elastic restoring force, with damping and non-conservative terms provided through the flow effects. Three different configurations are considered: a clamped panel, a hinged panel, and a flag (a cantilever clamped at the leading edge, free at the trailing edge). After providing the functional framework for the dynamics, recent results on well-posedness and long-time behavior of the associated dynamical system for solutions are presented. Having provided this theoretical context, in-depth numerical stability analyses are provided, focusing both at the onset of flow-induced instability (flutter), and qualitative properties of the post-flutter dynamics across configurations. Modal approximations are utilized, as well as finite difference schemes.


Introduction
In this treatment we provide a multifaceted numerical study-from the dynamical systems point of view-for a simplified class of 1-D models representing the physical phenomenon known as aeroelastic flutter. Flutter is fluid-structure (feedback) instability that occurs between elastic displacements of a structure and responsive aerodynamic pressure changes at the flow-structure interface [7,21,11,12]. The onset of flutter, for particular flow parameters, represents a bifurcation in the linearization of the system [26], and the qualitative properties of the post-flutter dynamics can be analyzed from an infinite dimensional/control-theoretic point of view [27,12,16].
Beam (or plate) flutter is a topic of great interest in the engineering literature [7,18,1,38,44,49,45,46], as well as, more recently, the mathematical literature [12,28,34,35] (and many references therein). Though there is a vast (mostly engineering-oriented) literature on flow-structure instabilities, we provide a straight-forward analysis utilizing a 1-D structural model in multiple configurations of interest. In particular, we consider both linear and nonlinear extensible beam models with pistontheoretic flow terms that, although simplified, demonstrate good correspondence of the dynamics with

Focus of This Study and Relation to Previous Literature
This study considers both pre-and post-flutter extensible, piston-theoretic beams in three configurations, from two numerical points of view: (i) a "modal" analysis, and (ii) a traditional finite difference scheme in space, both described in detail below.
We consider the effect of the piston-theoretic flow terms on the linear and nonlinear beam. We offer a thorough study of each key parameter's effect on stability, in the sense of: convergence to an equilibrium, onset of dynamic instability, convergence to a limit cycle oscillation (or chaos), and qualitative properties of non-transient dynamics. The analysis spans the three most physically relevant configurations appearing in the literature, and, unlike most previous work (e.g., [26,46,30]), we provide comparisons of the dynamics across the configurations. The inclusion of structural nonlinearity distinguishes this work from [46,45,43,30], allowing us to discuss qualitative properties of post-flutter behaviors. Our detailed parameter studies are done while bearing in mind the large body of recent theoretical results, in particular those in [28,12], and those in [29] for cantilevers with piston-theoretic terms.
Our analysis is done in the spirit of the impressive (by now classical) references [26,27] that study a piston-theoretic 1-D, H panel, allowing for in-plane pre-stressing. These papers consider various unstable behaviors, including divergence/buckling (static instability) and flutter (dynamic instability): the earlier [26] performs a tremendously thorough ODE analysis for a truncated-low mode-system, while [27] provides an infinite dimensional framework for discussions instability. We also mention the classical reference [19], where the same pre-stressed panel model is reduced to a two degree of freedom dynamical system via spectral methods, and considered as a deterministic example exhibiting chaos (for sufficiently large flow and stressing parameters).
For more general studies of fluttering beams and plates (with extensive references), we point to [21]. In the discussion of piston theory, we mention the modern engineering references [45,46], classic engineering references [1,7], and the recent mathematical-minded surveys [11,12]. Other mathematical references addressing piston-theoretic models include [28,13,6,14]. The article [28]: (i) performs numerical simulations (which guide the simulations herein); (ii) theoretically addresses various classes of attractors arising in the C beam/plate configuration. We follow much of that analysis there, in particular tracking stability/instability of dynamics with respect to piston-theoretic terms, as well as types and size of damping effects. The recent [29] follows the general theory/numerics approach in [28], but focuses specifically on the CF configuration and the effect of rotational inertia in the beam.

Modeling and Discussion
We now address the modeling pertinent to the theoretical discussion and numerical studies below.
Though there are various ways to consider flow-beam coupling, the simplest is to eliminate the fluid dynamic variables altogether. Such a tack has the benefit of reducing stability questions to a single, non-conservative dynamics. Additionally, the reduction eliminates moving boundary issues for the flow (particularly for the cantilever). On the other hand, such a reduction is a dramatic simplification of complex, multi-physics phenomena; however, focusing on the simple model allows us use robust mathematical theory and perform a thorough numerical study that can be exposited straight-forwardly. (More sophisticated flow-structure models are certainly explored in the rigorous mathematical literature [15,47,14,16]).
Motivated by the bulk of engineering studies cited thus far, we consider beam dynamics interacting with a potential flow. For high Mach numbers, the dynamic pressure on the surface of the structure, given by p(x, t), can be approximated point-wise in x by the pressure on the head of a piston moving through a column of fluid [46,1]. The pressure can then be written in terms of the down-wash of the fluid W = (∂ t + U ∂ x )w, where w(x, t) is the transverse displacement of the beam, and U is the unperturbed axial flow velocity. This results in a nonlinear expression [1] which is then linearized to produce the so called linear piston theory used on the RHS of the beam: Above, p 0 (x) is a static pressure on the surface of the beam (for the numerical portion of this paper we will take p 0 (x) = 0). The parameter β > 0 is a fluid density parameter that typically depends on U , though numerically we consider β decoupled from U . 3 For the structural model, we assume a slightly extensible beam [33]. Taking the standard Kirchhofftype hypotheses [4,32] leads to an Euler-Bernoulli beam, or the so called Rayleigh beam when allowing for rotational inertia in beam filaments [4]. For nonlinear effects, we invoke a large deflection elasticity model, taking into account quadratic terms in the strain-displacement relation [32]. From this, we obtain the beam analog of the full von Karman plate equations [33], that take into account both in-plane u and out-of-plane w (Lagrangian) dynamics. The principal nonlinear effect accounts for the local effect of stretching on bending (thus requiring beam extensibility). Let (w, w t ) correspond to the out-of-axis (transverse) dynamics, and (u, u t ) the in-axis (longitudinal) dynamics. Rotational inertia in beam filaments is represented by α ≥ 0, and D 1 , D 2 ≥ 0 are elastic coefficients. 4 The boundary conditions BC(u, w) are determined by the configuration to be studied.
The beam model under consideration here is a simplification of the above system (1.2): take u tt ≈ 0 and solve the resulting equation for u in terms of w, reducing the system to transverse dynamics. In eliminating the u dynamics, we enforce the condition that u(L) is zero-an unproblematic assumption for C and H, but certainly an imposition for CF (see Remark 1.3). A non-zero choice corresponds (mathematically) to pre-stressing the beam at equilibrium. Completing the simplification, one obtains the beam dynamics considered here: Above, we renamed or introduced various parameters: b 1 ∈ R measures in-plane stretching (b 1 < 0) or compression (b 1 > 0) at equilibrium (pre-stressing); b 2 > 0 gives the strength of the effect of stretching on bending (i.e., the nonlinear restoring force). We forgo a detailed discussion of the inertial parameter α ≥ 0, as well as associated structural damping in the model, until the next section.
Remark 1.2. This model has been historically referred to as the Krieger-Woinowsky beam [48,37], but has also been referred to as Kirchhoff or Berger-the latter used owing to the plate equation of the same nonlinear structure [28,13]. Classic references [2,3,5,17,23] discuss well-posedness and long-time behavior for extensible beams and Berger plates, typically with C or H boundary conditions.
The notation BC(w) will provide the boundary conditions corresponding to the pertinent configuration (C, H, or CF) given mathematically in the next section. These are standard boundary conditions, except in the CF case when imposing u(L) ≈ 0 necessitates a new condition: It is certainly noteworthy that the inherited boundary condition is nonlinear and nonlocal (but collapses to the standard free boundary condition if b 1 = b 2 = 0). (1.4) is mathematically non-trivial, as well as numerically challenging- [29] addresses this condition, and we elaborate on it below in Section 6.
Remark 1.3. For the cantilever, the restriction u(L) = const. throughout deflection is not entirely physical, at least in the absence of external mechanical effects. Indeed, (1.4) implicitly confines the free end of the beam to a fixed vertical line throughout deflection, thereby inducing stretching that provides the nonlinear restoring force. For a post-flutter cantilever, free end displacements can be quite large, with non-trivial inertial effects [22,44]. However, in using piston theory we assume larger values of U , known to reduce the magnitude of u(L), thereby improving the nonlinear structural model's viability. With this in mind, a "good" cantilever flutter model should allow for in-plane displacements, and should not primarily consider extensible effects. However, the added complexity of the nonlinearly coupled u/w system in (1.2) or the recent inextensible beam models [22,44,41,49] dramatically increases computational and analytical difficulties. Rigorous results for such models are few.
In summary: Our studies here provide a thorough analysis of simple, nonlinear fluttering beams, across various configurations. We acknowledge that the model considered here is greatly simplified: (i) it is only 1-D, not accounting for in-axis dynamics; (ii) nonlinear effects studied here are strictly due to extensibility, and we do not consider more sophisticated nonlinear models; (iii) we do not fully model flow dynamics, instead relying on piston-theoretic simplifications. Yet, even this simplified model yields interesting behaviors worthy of rigorous study, complemented by a thorough numerical examination done from a dynamical systems point of view. Our conclusions verify empirical conclusions from the engineering literature, as well as provide validation for recent abstract (infinite dimensional) results for beams and plates [28,29]. Our studies also serve as a baseline for future studies of more complex models, testing conjectured similarities between these dynamics in certain regimes.

PDE Model and Configurations
Recalling from Section 1.2: D > 0 is the elastic stiffness parameter, while the b 1 ∈ R corresponds to inaxis forces and b 2 ≥ 0 to the nonlinear restoring force; the function p 0 (x) represents static pressure on the surface of the beam, with flow parameters β > 0 (density) and U ≥ 0 (unperturbed flow velocity). The parameter k 0 ≥ 0 gives the frictional (or weak) damping coefficient, while k 1 ≥ 0 measures square root type damping (discussed in Section 2.1).
We distinguish between panel and cantilever configurations below, with a key distinction being the allowance of rotational inertia. For the panel configurations, C and H, the dynamics are given by: For the cantilever configuration CF, we allow α ≥ 0: In both cases above, appropriate initial conditions are specified: Our numerical studies below will compare dynamics across all three configurations holding α = 0; however, in the discussion of theoretical results, we will address technical issues associated with inertia and damping specific to the CF configuration, and thus in said case we allow for α ≥ 0.
Remark 2.1. Though the clamped conditions have nice mathematical properties (e.g., H 2 0 displacements and preservation of regularity via extension by zero), the hinged boundary conditions are the easiest to work with numerically because the in vacuo modes are purely sinusoidal. Including a clamped or free condition gives rise to hyperbolic modes, as well as eigenvalues obtained via nontrivial transcendental equations. Additionally, the free end yields the aforementioned complications. On the other hand, we will see that it is easiest to destabilize the CF configuration, followed by H, with the C beam most resistant to the flutter instability.

Rotational Inertia and Damping
Standard aeroelasticity literature [18,19,46] omits the structural rotational inertia. A scaling argument is typically invoked for "thin" structures, yielding α = 0 [4,32]. Mathematically, the presence of α > 0 is regularizing and nontrivial: w t ∈ H 1 [14]. Thus, α > 0 is helpful, and results with α = 0 are considered stronger [47,15,14,13]. As such, the focus of the numerical study here is on α = 0. For the CF configuration the theoretical situation surrounding α is more complicated. The recent work [29] fully explores well-posedness and long-time behavior (global attractors [14]) in the CF configurations. We present these results below in Section 3.
To qualitatively study post-flutter dynamics, we here consider a simple nonlinear restoring force [48,36,37] described above. In the panel configurations C and H, the scalar nonlinearity can be treated as a locally Lipschitz perturbation, and the nonlinearity introduces no boundary terms. In the CF configuration the nonlinearity interacts with the free boundary condition (1.4), and nontrivial boundary traces appear in the analysis of trajectory differences. These boundary terms are only controlled with α > 0 (or by introducing boundary damping [36,37] or other velocity smoothing).
Another interesting issue pertains to the strength of appropriate structural damping mechanisms as a function of α ≥ 0. The term k 0 measures the strength of standard weak structural damping. For α = 0, this is an appropriate strength of damping for stabilization; when α > 0, it is clear from the standard energy expression (see Section 2.3) that one requires k 1 > 0 in order for the damping to control of kinetic energies. (Interesting spectral questions emerge for the less natural cases of α = 0 with k 1 > 0 and k 0 = 0, as well as α > 0 with k 1 = 0 and k 0 large-addressed below.) When k 1 > 0, we say that the damping is of the "square root" type, as ∂ 2 x is formally "half" the order of the principal stress operator B = ∂ 4 x [40]. This concept can be generalized to fractional powers [B] θ w t for θ ∈ [0, 1]-Kelvin-Voigt damping occurs at θ = 1; see the more detailed discussion in [29], and references [39,8,25].
Lastly, we point out the clear disparity between the inclusion of inertia α > 0, and the aerodynamic damping provided from the flow itself in (1.1): piston theory provides damping of the form βw t (with dissipative sign), yet, if one includes inertia, this "natural" weak damping is not adequate to stabilize energies. Thus, we have a mathematical incentive to take α = 0 if we wish to utilize the flow-provided damping to help bound or stabilize trajectories.

Notation
For a given domain D, its associated L 2 (D) will be denoted as || · || D (or simply || · || when the context is clear). Inner products in L 2 (D) are written (·, ·) D (or simply (·, ·) when the context is clear). We will also denote pertinent duality pairings as ·, · X×X , for a given Hilbert space X. The space H s (D) will denote the Sobolev space of order s, defined on a domain D, and H s 0 (D) denotes the closure of

Energies and State Space
The natural energy for linear beam dynamics is given by the sum of the potential and kinetic energies Enforcing that α = 0 for H and C configurations, and that α ≥ 0 for CF, we have that the dynamics evolve in the state space H, whose definition depends on the configuration: Owing to the conditional structure of the state space (depending on the configuration and, for CF, the value of α) we opt for a compact notation, using We can then cleanly write the state space across all cases as with a clear meaning in a given context. We will also critically utilize the following nonlinear energies associated to equations in (2.1) and (2.2): where the Π term represents the non-dissipative and nonlinear portion of the energy: As in the general case [14, Lemma 1.5.4], we can see that for any w ∈ H, 0 < η ≤ 2 and > 0, This yields the crucial fact that the positive nonlinear energy is dominant: for some c 0 , c 1 , C > 0 depending on b 1 , b 2 .

Solutions
We now discuss the pertinent notions of solution. We will distinguish between: (i) panel configurations (C and H), for which we consider only α = 0, and (ii) the cantilever CF, where we admit α ≥ 0.
• Weak solutions satisfy a variational formulation of (2.1). One of the key features of such solutions is that w tt is interpreted only distributionally.
• Strong solutions are weak solutions possessing additional regularity that permits a classical (pointwise) interpretation of the evolution (2.1). For CF with α > 0, boundary conditions are still interpreted weakly due to subtleties associated with inertial terms. (See [29] or [14].) • Generalized solutions are C([0, T ], H) limits of strong solutions (when they exist). These are weak solutions [14,13], but they admit smooth approximation, and thus inherit some properties of strong solutions holding with respect to the topology of H.
Precise definitions are now given. For weak solutions, the variational form of the problem is the same for all three C, H, CF, with the space H encapsulating the configuration and α ≥ 0.
where d/dt is taken in the sense of D (0, T ). Moreover, for any (χ, ψ) ∈ H Definition 2 (Panel strong solution). For C or H configurations, a weak solution to (2.1) is strong if it possesses the regularity: (w, w t ) ∈ L ∞ (0, T ; H 4 (0, L) × H) and d + dt + w t is right-continuous. More regular cantilever solutions are considered, but only for α > 0.
Definition 3 (Cantilever strong solution). For the CF configuration, we say a weak solution to (2.2) with α > 0 is strong if it possesses the regularity: (Such solutions satisfy the third-order boundary condition weakly-see [29].) Finally, the most convenient notion of solutions will be those induced by a semigroup flow:

Well-posedness
For all results below D, k 0 , β, U ≥ 0. We employ the energetic notations from Section 2.3.
. Generalized (and strong) solutions obey the energy identity for all t ≥ 0: is semigroup well-posed on H CF . Generalized (and strong) solutions obey the energy identity for all t ≥ 0: For such a linear problem, these results are well-established.
Similarly, for configuration CF with α > 0 and k 1 ≥ 0, the problem in (2.2) is (nonlinear) semigroup well-posed on H CF . Generalized (and strong) solutions obey the energy identity for all t ≥ 0: In the case of nonlinear clamped and hinged beams, the well-posedness has been known for sometime; consult the abstract theory in [9,14] and note the classical references (for both 1-D Krieger and 2-D Berger nonlinearities) [2,3,5,17,23]. In the case of the cantilevered beam CF with α ≥ 0, see the recent work [29] (which includes discussion of other cantilevered beams [36,37]).
The final case not addressed by the previous theorem is the nonlinear cantilever in the absence of rotational inertia. For this case [29] provides only an existence result.

Long-time Behavior of Trajectories
We can obtain the following global-in-time stability bounds for all generalized solution to (2.1) or (2.2) in the presence of piston-theoretic terms, β ≥ 0 and U ≥ 0: Theorem 3.4. Let D > 0 with β > 0 and U ≥ 0. In configurations C and H, we have the following global in time bound for generalized solutions to (2.1) for any amount of weak damping k 0 ≥ 0: In the case of configuration CF, with α > 0 (rotational inertia present) and k 1 > 0 (active square root damping), we have sup Remark 3.1. Formally, CF trajectories for (2.2) with no inertia, α = 0, do obey a global-in-time bound of the form (3.8), but the requisite multipliers (and associated integration by parts) are not permissible for the weak solutions guaranteed by Theorem 3.3. These global bounds are nontrivial, owing to the fact that the dynamics are non-gradient, as reflected by the energy-building terms under temporal integration in (3.6) and (3.7). They follow from the estimates in (2.6)-(2.7), taken together with a Lyapunov argument, and rely critically on the "good" structure of Π, reflected by the control of lower order terms in (2.6). In fact, more is true. For a dynamical system, a compact global attractor is a uniformly attracting (in the sense of the H topology), fully time-invariant set [14].
For configurations C and H with k 0 ≥ 0 the dynamical system (S(t), H) generated by generalized solution to (2.1) has a smooth and finite dimensional compact global attractor.
For configurations CF with α, k 1 > 0 and k 0 ≥ 0, the dynamical system (S α (t), H CF ) generated by generalized solution to (2.2) has a finite dimensional compact global attractor.
In the above theorem, smooth is taken to mean that trajectories on the attractor have (w, w t ) ∈ H ∩ H 4 (0, L) × H , where H indicates the configuration-dependent displacement space from (2.4); the attractor is bounded in this topology. Finite dimensional refers the fractal dimension of the attractorsee [14] for more details.
In this study, the attractor "contains" the essential post-flutter dynamics for (2.1) or (2.2). The aeroelastic model is non-conservative, thus the associated dynamical system has no strict Lyapunov function [14]. In this case the global attractor is not characterized as the unstable manifold of the equilibria set, so to obtain a compact global attractor one must show both asymptotic compactness and ultimately dissipativity of the dynamical system [14]. Since the infinite dimensional dynamical system associated to solutions has the property that trajectories converge-in a uniform sense-to an attractor, we effectively reduce the analysis of the dynamics to a "nice" finite dimensional set.
The existence and properties of the attractor for panel configurations are shown in [28] (though earlier proofs have certainly appeared for these configurations). This work also discusses the application of a recent tool: the theory of quasi-stability [10,14]. The theory provides the existence of so called exponential attractors [10,14,24]. A generalized fractal exponential attractor for the dynamics (S(t), H) is a forward invariant compact set A exp ⊂ H, with finite fractal dimension (possibly in a weaker topology), that attracts bounded sets with uniform exponential rate. We remove the word "generalized" if the fractal dimension of A exp is finite in H.
For configurations C and H with k 0 ≥ 0 the dynamical system (S(t), H) generated by generalized solution to (2.1) has a generalized fractal exponential attractor.
In the cantilever configuration CF, the result on the existence of a compact global attractor is much more recent [29].

Overview of Numerical Studies
Our numerical investigations of the piston-theoretic beam are multi-faceted; we seek to determine critical parameter values that lead to the onset of instability, and to understand essential qualitative properties of post-flutter dynamics. To be consistent with comparable engineering studies, and to be consistent across configurations, we neglect rotational inertia here.
In all of our experiments, we focus on the case α = 0 (and thus k 1 = 0).
We also take the static pressure p 0 (x) ≡ 0. In Section 5 we use physical parameters to emphasize the predictive capabilities of our methods. In Section 6, our primary interest is the qualitative properties of post-flutter nonlinear dynamics, and, as such, we choose mathematically convenient parameter values.
We utilize two main methods for these investigations: modal analysis (described in detail in Section 5.1) and finite differences. To analyze the onset of instabilities in the linear model, we use modal analysis as a predictive tool for finding critical parameters, and we use finite difference approximations to determine if these predictions are in agreement with the dynamics. Subsequently in Section 6, we peform finite difference simulations of the nonlinear model to investigate the post-flutter regime.
Finite difference simulations are accomplished via standard second-order difference formulas for spatial derivatives, with ghost values beyond the spatial boundary as necessary for the various boundary configurations. Temporal integration is accomplished using a stiff ODE solver. Our simulations utilize three basic types of initial data (configuration dependent): (i) modal initial displacement, (ii) polynomial initial displacement, and elementary initial velocity. Let us denotex = x/L.
• [nth Mode ID] For (i), we consider the in-vacuo mode shapes of a given configuration (detailed below in Section 5.1) as an initial displacement, with zero initial velocity profile. We restrict to mode numbers n = 1, ..., 5. For instance, in the case of CF conditions we have: where k 2 ≈ 4.6941 is the second Euler-Bernoulli cantilevered mode number and C 2 ≈ 1.0185.
• [Polynomial ID] For (ii), with the C and F configurations, we consider a polynomial of minimal degree satisfying the boundary conditions for the respective configuration in the experiment. For CF, we have We again take the initial velocity profile to be zero.
• [Elementary IV] For (iii), we take the initial displacement to be identically zero, and consider a simple initial velocity profile. In the case of C or H, we will take a symmetric quadratic function on [0, L]; for CF we simply take a linear function of positive slope on [0, L] vanishing at x = 0: In the sections to follow, we present our numerical results. In some sections we compare across boundary configurations, but for qualitative studies we often choose a single, particularly illustrative example in a given configuration. In addition, we sometimes consider how qualitative aspects of the simulations may or may not vary across differing initial conditions. Several visualizations involve energy plots (as a function of time), as well as plots of midpoint displacements (in time) for panel (C or H) configurations, or free end displacements for the cantilever (CF) configuration.

Linear Computational Analysis: Flutter Prediction
In this section we analyze the onset problem for flow-induced instability: namely, what combinations of parameters of interest-U, β for the flow, k 0 , L for the beam-yield stable or unstable asymptoticin-time dynamics. More specifically, our studies in this section we are concerned with convergence to stationary states or unbounded growth of displacements as t → ∞. (Note: when b 2 = 0, without the nonlinear restoring force, unstable trajectories will grow unboundedly.) For the readers convenience, we restate the system studied in this section here: For convenience, let us define the damping parameter k ≡ k 0 + β, which consists of both imposed damping k 0 ≥ 0, as well as piston-theoretic (flow) damping β > 0. The effects of the initial conditions are conservative (which is to say that, for β = k 0 = 0, the energy is conserved throughout deflection). With no imposed damping (k 0 = 0), the transient dynamics may or may not damp out; in this case, the only damping in the problem, then, is contributed by the flow in the form of βw t . Additionally, with U > 0, the non-dissipative piston term −βU w x on the RHS of (2.1) may induce instability (and thus oscillatory behavior in the linear energy E(t)). In this sense, stability can be seen as competition between the effect of the damping in βw t bleeding energy out of the system, and the non-dissipative flow effect U w x putting flow energy in. As U is increased the effect of the flow becomes more pronounced and the location of the spectrum is shifted. For all other parameters fixed (including β), a certain value U crit represents a bifurcation, and the dynamics become unstable and exponential growth (in time) of the solution is observed. With imposed damping in the model (k 0 > 0 and L, β fixed), there exists a U crit (k 0 ) such that for U < U crit , the transient dynamics are damped out, and for U > U crit instability will still be observed. We can think of U crit as an increasing function of the damping k 0 in the problem.

Modal Analysis
Modal analysis, often referred to in engineering literature, can in fact have various meanings. Generally, it refers to a Galerkin method whereby solutions to a structural or fluid-structural problem are approximated by in vacuo structural eigenfunctions ("modes") (e.g., [30,45,46,18,22]). Since the eigenfunctions of standard elasticity operators for a basis for the state space, a good well-posedness result for the full system justifies this type of approximation. This type of approximation can be dynamic, as in reducing an evolutionary PDE to a finite dimensional system of ODEs by truncation, or it can be stationary, reducing the problem of dynamic instability (for linear dynamics) to an algebraic equation. The latter is what we utilize and describe in detail.

Spatial Modes and Eigenvalues
Critical to any modal analysis are the in vacuo modes (eigenfunctions) associated with each configuration: C, H, CF. Since we are working with the Euler-Bernoulli beam (α = 0) for our simulations, the modes and associated eigenvalues can be computed in an elementary way. These functions are complete and orthonormal in L 2 (0, L), as well as complete and orthogonal in H 2 (0, L) (with respect to (·, ·) H , where H encodes the boundary conditions (2.4)).
So the modes in any configuration are linear combinations of the above, determined by invoking the boundary conditions; the boundary conditions also provide the specific values (κ n , ω n ).
For C and CF in Table 1, the mode numbers κ n L are obtained by numerically solving the characteristic equation. We have C : C n = −c n cos(κ n L) − cosh(κ n L) sin(κ n L) − sinh(κ n L) CF : C n = −c n cos(κ n L) + cosh(κ n L) sin(κ n L) + sinh(κ n L) .
The c n values are chosen to normalize the functions in the L 2 (0, L) sense.

Reduction to Perturbed Eigenvalue Problem
Let us consider the Galerkin procedure for the full linear equation beam equation: on (0, L), with boundary conditions according to the configuration being studied; recall that k = k 0 + β. We view the terms involving β, U, k 0 as perturbations and expand the solution via the in vacuo mode functions {s j } as w(t, x) = q j (t)s j (x). The {q j } represent smooth, time-dependent coefficients.
Plugging this representation into the (5.2), multiplying by s n , and integrating over (0, L) for each n we obtain: Orthonormality of the eigenfunctions can be invoked to produce diagonal terms, whereas the terms scaled by βU (∂ x s m , s n ) are off-diagonal and give rise to the instability of the ODE system.
Remark 5.1. To obtain a spectral simulation for the dynamics, one could simply truncate the system above for all 1 < m, n < N and solve the resulting N × N system for the unknown coefficients {q j (t)}, To simply determine the stability of the problem as a function of the given parameters, we invoke a standard engineering ansatz [46,21] (and references therein): assume simple harmonic motion in time of the beam according to some dominant (perturbed) frequency ω; we allow possible contribution from all functions s n for n = 1, 2, ..., N via coefficients labeled α n : (5.4) where N is a chosen dimensional truncation. We repeat by multiplying by mode functions s n , and then integration produces an eigenvalue problem in the perturbed frequency ω. With the off-diagonal entries (∂ x s m , s n ) in hand (for 1 ≤ m, n ≤ N with m = n), compute diagonal terms we create the matrix for 1 ≤ n, m ≤ N : for m = n. (5.5) As an example, for the configuration H, the entries in the matrix can be computed explicitly. Setting N = 4 yields the linear system For chosen parameter values of D, k 0 , β, U, L, we enforce the zero determinant condition for non-trivial solutions in ω: det A( ω) = 0, and solve for ω = [ ω 1 , ..., ω N ] T . The associated complex roots allows us to track the response of the natural modes to the perturbation terms-damping, (k 0 + β)w t and non-dissipative, βU w x .
• If Im( ω j ) > 0, then we call the configuration unstable, and the associated linear dynamics should exhibit exponential growth in time with rate Im( ω j ) and frequency Re( ω j ).
• Otherwise, the linear dynamics remain bounded in time.
Obviously, this method of determining stability has associated error-first, owing to the simple harmonic motion ansatz, and secondly, because of truncation. Testing the conclusions drawn from this type of modal stability determination against a full numerical scheme is one of the main numerical points of this consideration.
Remark 5.2. The emergence of an unstable perturbed eigenvalue (associated to Im( ω)) is precisely what is meant by self-excitation instability, which is commonly called flutter. In this way, we note that predicting the onset of flutter is linear. As mentioned above, for one to "observe" (numerically) the post-onset/post-flutter dynamics, a nonlinear elastic restoring force must be incorporated into the model to keep solutions bounded in time (as the underlying linear dynamics are driven by the aforementioned exponential growth).
Remark 5.3. From a practical point of view, one can make use of any reasonable N for modal calculations involving the H configuration-these mode functions are purely sinusoidal. In modal calculations involving CF or C, we must first solve for the k n , which is more numerically unstable for larger n. Additionally, these approximate k n are used in the numerical integration of cosh and sinh functions with large arguments, further propagating error. Therefore, for N ≥ 7 for C and CF, modal calculations are cumbersome and possibly inaccurate. Engineers have devised clever approximation of the mode functions is needed to circumvent this issue [20], though it does not concern us too much here as flutter is considered a "low-mode" phenomenon. As evidenced below, modal approximations perform quite well, even for N = 6.

Prediction of the Onset of Instabilities
We now predict the onset of instability for the model (2.2) as the U parameter increases, with other parameter values fixed and the nonlinear parameter b 2 = 0. Indeed, there exists some critical flow speed value, U crit , beyond which the system exhibits exponential growth in time.
In this section, we have adopted the non-dimensionalization procedure utilized in [45]. The parameters chosen here represent physical values for the prediction of the flutter phenomenon. Accordingly, we utilize these as the fixed baseline values throughout this section: • D = 23.9 (fixed across all simulations), • L = 300 (allowed to vary in the range 100 ≤ L ≤ 500, • β = 1.2 · 10 −4 (allowed to vary in the range 10 −5 ≤ β ≤ 10 −2 ). When employing finite difference simulations of (2.2), the task of identifying U crit usually requires several batches of simulations, usually in some bracketing approach. In general it is not practical, as qualitative observations of flutter are somewhat subjective and may require excessive long-time simulations (for instance, accurately determining U crit for L > 320). However, the modal approach outlined in Section 5.1.2 yields a direct, simple way to find U crit . For given parameter values, one can immediately analyze the perturbed eigenvalues associated with (5.5) to determine if an instability is present (Im(ω j ) > 0 for some j).
In Figure 1, we compare the results of a bracketing search for U crit with the modal approach for a range of L values across all three boundary configurations. The points are the U crit values found by a bisection-like search employing finite difference simulations, and the curves are the values predicted by the modal approach. We see that the modal approach is very consistent with the finite difference approach, and is substantially more convenient, for predicting U crit . The hierarchy of stability between the three configurations C, H, and CF is also revealed in Figure 1  Note that the modal approach also produces an ω crit associated to each U crit (not shown here), and the frequency of unstable oscillation is simply Re(ω crit ). Finally, as the modal analysis was such a good predictor of U crit for all configurations, the laborious finite difference simulations to search for U crit in the clamped-free case for L > 320 were not performed.
In Figure 2, we determine U crit as a function of β, with 10 −5 ≤ β ≤ 10 −2 . We observe monotonic linear decay in the U crit (β) value as β increases up to around 10 −3 , and then somewhat erratic behavior for larger β. Since β affects both the damping and destabilizing terms in the model (5.2), it is not obvious that its effects will remain monotonic on U crit as β increases further. Indeed, an analysis of the dispersion relation associated to (5.2) indicates that a critical transition occurs at β = 1. One could similarly extract critical values of L or β by fixing the other and varying U . Results maintain the same qualitative structure, for instance, with respect to the established stability hierarchy of boundary configurations and monotonic behavior of U crit .
Another relationship of interest is that of the damping coefficient k = k 0 + β on U crit . Note that in Figure 3 the region below a given curve represents stable combinations of k and U . The curves themselves represent the critical value U crit for any given damping coefficient k 0 . Conversely, given a U value, one can extract the necessary strength of damping to yield stability. Note that in cases where the boundary points are completely restricted (C and H), the relationship between U crit and k 0 is asymptotically linear, while for the CF configuration it is sublinear. The spurious behavior in the zoomed region on the right does not seem to be a numerical artifact, as the computations were done on a very fine scale. We conclude this section by empirically verifying the global nature of instability across initial conditions (for a given configuration). In Figure 4 we plot the dynamic energy as a function of time for the H configuration with U > U crit (i.e., U = 5 for L = 300, D = 23.9, β = 1.2 · 10 −4 , and k 0 = 0). The initial conditions are taken from among the first and second H modes (see Table 1), the polynomial initial displacement w 0 (x) =x 3 (1 −x) 3 , and the elementary initial velocity w 1 (x) =x(1 −x). We note that, in the above, all dynamic simulations reflect that U > U crit . Each simulation, in the absence of the nonlinear restoring force (b 2 = 0), shows the dynamics eventually growing exponentially. Moreover, it is clear from the log-scale plot-since all slope profiles are the same-that each simulation is growing in time by the same exponential factor, namely, the perturbed unstable eigenvalue. Also note, as we will observe in the non-linear case (b 2 > 0), there is a transient regime at the outset of the dynamics; however, with damping provided by the flow, βw t , the effect of the initial condition decays away after a short amount of time.

Nonlinear Computational Analysis: Qualitative Properties
In this section, we present a variety of simulations across configurations and parameter values to illustrate various aspects of the qualitative post-flutter behavior. In the simulations below, for the most part, we take the nonlinear parameter to be positive, b 2 > 0. This means, in particular, that the nonlinear restoring force is active. Thus, for any simulation with U > U crit , we will observe convergence to a limit cycle oscillation. We can then ask about the effect of any parameter (e.g., in-axis compression b 1 > 0) on the limit cycle.
For expediency of simulations, and also because the primary purpose of this section is qualitative, we allow parameters not of interest to be set to unity. This is to say that in Section 6, we take β = L = D = 1. Simulations below will consider varying b 1 , b 2 , k 0 and U .
For the readers convenience, we restate the system studied in this section here: (6.1) Implementation of the boundary condition for w xxx (L) in the CF configuration was computed by utilizing a 5-point stencil for w centered at x = L. The nonlinear terms in (6.1) involving b 2 ||w x || 2 were handled using the w x from the previous iterate (within the same time step) for the integral computation ||w x || = L 0 |w x | 2 dx.

Post-Flutter Examples
To ground our discussions, we begin by showing a collection of snapshots of the in vacuo linear (b 2 = 0) beam dynamics and corresponding post-flutter dynamics. In Figure 5  In the figures below, we provide a similar snapshot of the profile of a fluttering cantilever beam, taken with a polynomial initial displacement (see Section 4), rather than a modal initial displacement. We show approximately one period of the non-transient dynamics of a beam, after it has approached a limit cycle. On the left is the natural scale, and on the right, we decrease the y-scale to emphasize the nodal point. For comparison, we provide time snapshots of both the first and second cantilever modes in vacuo; note the similarity between the flutter profile and a superposition of the first and second in vacuo modes.

Energy Plots (Linear vs. Nonlinear)
In the following section we consider the C configuration. We fix the initial condition and primarily consider varying U around U crit , holding other parameters fixed. We focus on the behavior of the energies to illustrate the effect of including nonlinearity in the simulation; this is to say that, when b 2 > 0, the nonlinear restoring force is active and for U > U crit the trajectories remain bounded (see Section 3.3) and lock into a limit cycle. First, computed energies E(t) for the model (2.1) with no damping present (k 0 = −β = −1) and with no nonlinear or in-axis effects (b 1 = b 2 = 0) are shown in Figure 8. For this choice of k 0 = −1, we have U crit = 135.18. Energy profiles (log-scale) for various choices of U as multiples of U crit are given. Recall that a linear profile in the energy plot represents exponential growth or decay of the dynamics in the finite energy topology, with margin of in/stability depending on the slope of the profile. Note the exponential growth in energies for U > U crit while E(t) stays bounded for 0 < U < U crit , and is constant for U = 0. In the next simulation the damping parameter is set to k 0 = 0 so k = 1, and thus there is damping present from the flow (but no structural or imposed damping) in (2.2). This results in an approximate critical flow velocity of U crit = 135.9 (the presence of damping perturbs the critical flow velocity corresponding to the onset of instability by roughly 0.5%). Computations again were performed for various choices of U , with energies shown in Figure 9. Note the effect that the damping has on the decay of energies for all U < U crit (as E(t) → 0 when t → ∞)-each curve for U < U crit has a linear profile (or envelope) with negative slope, indicating exponential decay. Energy values for U ≥ U crit still grow exponentially. Finally, in Figure 10, the clamped, nonlinear beam (with b 1 = 0 and b 2 = 1) is considered and the nonlinear energy E (t) is shown. Note that E (t) (and hence E(t), due to (2.7)) now remains bounded for all supercritical (unstable) velocities U > U crit = 135.9, demonstrating the stability that is induced by nonlinear restoring force. From the point of view of trajectories, each remains bounded for all time, with global-in-time bound dependent upon U . For U < U crit , we note that the nonlinearity does not dramatically affect the stability observed in the linear case-consistent with the semilinear nature of the nonlinearity. To show the detailed difference between the energies E(t) and E(t), Figure 11 shows both quantities for k 0 = 0, b 1 = 0, and b 2 = 1, with U = 2U crit .

Convergence to Limit Cycles
In this section we demonstrate some properties of trajectories that converge to a limit cycle oscillation, as described above. Figure 12 shows the endpoint displacement a cantilevered beam CF with supercritical parameter values. It is clear that the envelope of the oscillations grows significantly and until around t = 1.5, and then undergoes an exponential decay to a nontrivial limit cycle behavior. In Figure 13, plots of the beam midpoint displacement are shown for a very large flow velocity (U = 5000), significant in-axis compression parameter (b = 50), and b 0 = 1 for selected values of the damping parameter k = k 0 + 1. The sensitivity of the dynamics to the damping parameter can be seen by noting the relative rate of initial decay of energy as k increases-note the quick decay to the limit cycle for larger values of k, though the limit cycle itself is preserved across damping values.

Uniformity of Limit Cycle
Next, we consider the CF configuration and provide tip profile at x = L = 1 for three different initial configurations: where κ 2 ≈ 4.6941 is the second Euler-Bernoulli cantilevered mode number (with L = 1) and C 2 = cos(κ 2 ) + cosh(κ 2 ) sin(κ 2 ) + sinh(κ 2 ) ≈ 1.0185; First, Figure 14 represents the nonlinear (b 2 = 1), in vacuo (β = k 0 = 0) tip displacements with the three initial configurations described above.  In these cases, we observe convergence to the "same" limit cycle for each of the initial conditions considered. Finally, note that, although we seem to see convergence to the same LCO, the transient regime and "time to convergence" are certainly affected by the choice of initial configuration.

Convergence to Nontrivial Equilibrium
For certain parameter combinations, the presence of excess in-plane compression leads to a "buckled" plate/beam configuration (a bifurcation, for the static problem). For b 2 = 1, the parameter b 1 = 50 is large enough for U = 100 to impart this behavior. The transient behavior decays more rapidly to the nontrivial steady state for larger values of k. A plot of the energy E(t) for different values of k is given in Figure 16. Regardless of the choice of k, all simulations for U = 100, b 2 = 1, and b 1 = 50 converge to the same nontrivial steady state and same linear energy level set. A plot of the nontrivial steady state w (the buckled beam displacement) is given in Figure 17. It is also possible to observe different choices of damping k resulting in convergence to different steady states. For U = 100, b 2 = 1, and b 1 = 100, the choices k = 1 and k = 2 actually produce nontrivial steady states that are negatives of one another. In Figure 18 the midpoint displacement is plotted for these two cases.

Non-simple Limit Cycle
It is possible to induce several different phenomena by manipulating parameter values for the nonlinear dynamics. In Figure 19 the midpoint displacement is shown for U = 5000, k = 100, b 1 = 5000, and b 2 = 1000. Note the initial transient dynamics are damped out quickly and the dynamics converges to a non-simple limit cycle.

Chaos with In-plane Tension
Here we consider the following initial configurations for the nonlinear H beam taken with substantial in-plane compression.
We note from [19] that this setup is a candidate to produce a deterministic, chaotic dynamical system. (In the reference [19], the system (2.2) is reduced down to a two modal system, and the qualitative properties-convergence to equilibria, convergence to a simple limit cycle, convergence to a non-simple limit cycle, and "chaos"-are mapped in the (b 1 , U ) parameter space.) We produce an interesting example below with large U . We note that for the given value of b 1 below, we are past the Euler buckling bifurcation point. First, from the four profiles corresponding to 0, we see no discernible notion of convergence of the dynamics. Secondly, we see no periodic behavior corresponding to the initial configuration with = 0. Additionally, it is not clear that the other configurations involving > 0 converge to a true limit cycle. In these cases, it is clear there is no uniformity of end behavior across initial configurations, as was the case in the previous section.
Remark 6.1. We continued simulations for smaller values of , and it was clear that no trend in magnitude or period of the oscillations emerged, even at very fine scales and long run-times.
The only consistent feature across the above simulations corresponds to their energy curves. All energy plots show the same features present in Figure 21; namely a transient period followed by an energy "bump" up to a higher level. The transition is smooth.

Effect of Magnitude of b 2
The effect of increasing the nonlinear parameter b 2 was studied at a supercritical flow velocity (U > U crit ) for beam. In Figure 22, energy profiles for several different choices of b 2 are given for the parameter k 0 = 0 and U = 150 (U crit = 135.97). Note that as b 2 increases, the energy plateau decreases (for an initial configuration fixed across all simulations).