A destination-preserving model for simulating Wardrop equilibria in traffic flow on networks

1.   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,10,20]. The onset of flutter, for particular flow parameters, represents a bifurcation of the system ^[24], and the qualitative properties of the post-flutter dynamics can be analyzed from an infinite dimensional/control-theoretic point of view ^[10,14,25].

Beam (or plate) flutter is a topic of great interest in the engineering literature ^{[1,7,16,17,35,41,42,43,46]}, as well as, more recently, the mathematical literature ^{[10,26,32,33]} (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 piston-theoretic flow terms that, although simplified, demonstrate good correspondence of the dynamics with what is empirically known about the onset of flutter and post-flutter behaviors ^[1,42]. We will demonstrate the remarkable ability of these simplified equations to capture a robust set of un/stable behaviors associated with structures in axial flow. We also provide modern, rigorous context for these numerical studies by discussing recent theoretical results for the nonlinear dynamical systems associated to the PDE models given here.

Our studies primarily take the form of computational stability analyses for a beam distributed on $ x \in [0, L] $, immersed in an inviscid potential flow of sufficiently high Mach number ^[1]; we consider three beam configurations: clamped-clamped panel (C), hinged-hinged panel (H), and clamped-free cantilever (CF). The fluid flows with steady velocity $ U \ge 0 $ in the $ x $-direction (along the length of the beam), in contrast with so-called normal flow; thus, $ x = 0 $ is the leading edge, and $ x = L $ is the trailing edge. The cantilever CF is taken to be a flag, clamped at the leading edge, and free at the trailing edge.^* The baseline linear structural model we consider is the standard Euler-Bernoulli beam, though for the CF configuration we mention the so called Rayleigh beam (allowing for rotational inertia).^† Additionally, we consider extensible nonlinear beams in sense of Woinowsky-Krieger ^[45]—see ^[27] for more details. For the flow theory, we employ the linear piston theory, which reduces the flow dynamics to a simple non-conservative right hand side, written through the material derivative of structural variable ^[1,25].

^*A recent configuration of interest is the so called inverted flag, where the free edge leads; see ^[29,39].

^†See ^[4] for a nice description and comparison of four fundamental beam theories.

As with all flutter problems ^[20,35], and as we will demonstrate in the sequel, the onset of instability can be studied via linear theory—typically as an eigenvalue problem (see e.g., ^[42,43]). Indeed, if one can model the dynamic load across the structure, the issue is to determine under what conditions the non-conservative flow loading destabilizes structural (eigen)modes. However, to study the qualitative properties of the dynamics in the post-flutter regime, analysis will require some physical nonlinear restoring force to keep trajectories bounded in time ^{[10,14,16,17]}. Indeed, for the post-flutter dynamics of a given system, the theory of nonlinear dynamical systems is often used, e.g., compact global attractors ^[10,12] as well as exponential attractors ^[9,26].

1.1.   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: (ⅰ) a "modal" analysis, and (ⅱ) 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., ^[24,28,43]), we provide comparisons of the dynamics across the configurations. The inclusion of structural nonlinearity distinguishes this work from ^{[28,40,42,43]}, 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 ^[10,26], and those in ^[27] for cantilevers with piston-theoretic terms.

Our analysis is done in the spirit of the impressive (by now classical) references ^[24,25] 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 ^[24] performs a tremendously thorough ODE analysis for a truncated—low mode—system, while ^[25] provides an infinite dimensional framework for discussions instability. We also mention the classical reference ^[18], 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 ^[20]. In the discussion of piston theory, we mention the modern engineering references ^[42,43], classic engineering references ^[1,7], and the recent mathematical-minded surveys ^[10]. Other mathematical references addressing piston-theoretic models include ^[6,11,12,26]. The article ^[26]: (ⅰ) performs numerical simulations (which guide the simulations herein); (ⅱ) 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 ^[27] follows the general theory/numerics approach in ^[26], but focuses specifically on the CF configuration and the effect of rotational inertia in the beam.

1.2.   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 to 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 ^{[12,13,14,44]}).

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 ^[1,43]. The pressure can then be written in terms of the down-wash of the fluid $ W = (\partial_t+U\partial_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] in $ W $

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 $ \beta > 0 $ is a fluid density parameter that typically depends on $ U $, though numerically we consider $ \beta $ decoupled from $ U $.^‡

^‡See ^[43] for a discussion of the flow non-dimensionalization, and further discussion of characteristic parameter values.

For the structural model, we assume a slightly extensible beam ^[31]. Taking the standard Kirchhoff-type hypotheses ^[4,30] 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 ^[30]. From this, we obtain the beam analog of the full von Karman plate equations ^[31], 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 $ \alpha \ge 0 $, and $ D_1, D_2 \ge 0 $ are elastic coefficients.^§

^§The coefficients $ \alpha $, $ D_1 $, $ D_2 $ are related as $ \alpha D_1 = D_2 $ in the physical derivation of the beam equation—see ^[31].

The boundary conditions $ BC(u, w) $ are determined by the configuration to be studied.

Remark 1. At present, all rigorous mathematical analyses of the conservative (1.3) model (or the 2-D analog) require $ \alpha > 0 $ ^[31,27]. However, when considering purely transverse $ w $ dynamics of thin beams and plates, $ \alpha = 0 $ is almost always taken (see the discussion in [30,Ch.1]).

The beam model under consideration here is a simplification of the above system (1.3): take $ u_{tt} \approx 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 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 \in \mathbb R $ measures in-plane stretching ($ b_1 < 0 $) or compression ($ b_1 > 0 $) at equilibrium (pre-stressing); $ b_2 \ge 0 $ gives the strength of the effect of stretching on bending (i.e., the nonlinear restoring force). We delay a detailed discussion of the inertial parameter $ \alpha \ge 0 $, as well as associated structural damping in the model, until the next section.

Remark 2. This model has been historically referred to as the Woinowsky-Krieger beam ^[34,45], but has also been referred to as Kirchhoff or Berger—the latter used owing to the plate equation of the same nonlinear structure ^[26,11]. Classic references ^{[2,3,5,15,22]} 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)\approx 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.5) is mathematically non-trivial, as well as numerically challenging—^[27] addresses this condition, and we elaborate on it below in Section 6.

Remark 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.5) 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 ^[21,41]. 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.3) or the recent inextensible beam models ^{[21,38,41,46]} 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: (ⅰ) it is only 1-D, not accounting for in-axis dynamics; (ⅱ) nonlinear effects studied here are strictly due to extensibility, and we do not consider more sophisticated nonlinear models; (ⅲ) 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 ^[26,27]. Our studies also serve as a baseline for future studies of more complex models, testing conjectured similarities between these dynamics in certain regimes.

2.   PDE model and configurations

Recalling from Section 1.2: $ D > 0 $ is the elastic stiffness parameter, while the $ b_1 \in \mathbb R $ corresponds to in-axis forces and $ b_2 \ge 0 $ to the nonlinear restoring force; the function $ p_0(x) $ represents static pressure on the surface of the beam, with flow parameters $ \beta > 0 $ (density) and $ U\ge 0 $ (unperturbed flow velocity). The parameter $ k_0 \ge 0 $ gives the frictional (or weak) damping coefficient, while $ k_1 \ge 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 $ \mathbf{CF} $, we allow $ \alpha \ge 0 $:

In both cases above, appropriate initial conditions are specified: $ w(t = 0) = w_0, \; \; w_t(t = 0) = w_1 $, and we take $ p(x, t) = p_0(x)-\beta[w_t+Uw_x] $.

Our numerical studies below will compare dynamics across all three configurations holding $ \alpha = 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 $ \alpha \ge 0 $.

Remark 4. Though the clamped conditions have nice mathematical properties (e.g., $ H_0^2 $ 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.

2.1.   Rotational inertia and damping

Standard aeroelasticity literature ^{[16,17,18,43]} omits structural rotational inertia. A scaling argument is typically invoked for "thin" structures, yielding $ \alpha = 0 $ ^[4,30]. Mathematically, the presence of $ \alpha > 0 $ is regularizing and nontrivial: $ w_t \in H^1 $ ^[12]. Thus, $ \alpha > 0 $ is helpful, so results with $ \alpha = 0 $ are considered stronger ^{[11,12,13,44]}. As such, the focus of the numerical study here is on $ \alpha = 0 $. For the CF configuration the theoretical situation surrounding $ \alpha $ is more complicated. The recent work ^[27] fully explores well-posedness and long-time behavior (global attractors ^[12]) 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 ^[34,45] 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.5), and non-trivial boundary traces appear in the analysis of trajectory differences. These boundary terms are only controlled with $ \alpha > 0 $ (or by introducing boundary damping ^[34] or other velocity smoothing).

Another interesting issue pertains to the strength of appropriate structural damping mechanisms as a function of $ \alpha \ge 0 $. The term $ k_0 $ measures the strength of standard weak structural damping. For $ \alpha = 0 $, this is an appropriate strength of damping for stabilization; when $ \alpha > 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 $ \alpha = 0 $ with $ k_1 > 0 $ and $ k_0 = 0 $, as well as $ \alpha > 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 $ \partial_x^2 $ is formally half the order of the principal stress operator $ B = \partial_x^4 $ ^[37]. This concept can be generalized to fractional powers $ [B]^\theta w_t $ for $ \theta\in[0, 1] $—Kelvin-Voigt damping occurs at $ \theta = 1 $; see the more detailed discussion in ^[27], and references ^[8,23,36].

Lastly, we point out the clear disparity between the inclusion of inertia $ \alpha > 0 $, and the aerodynamic damping provided from the flow itself in (1.2): Piston theory provides damping of the form $ \beta 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 $ \alpha = 0 $ if we wish to utilize the flow-provided damping to help bound or stabilize trajectories.

2.2.   Notation

For a given domain $ D $, its associated $ L^{2}(D) $ norm will be denoted as $ ||\cdot ||_D $ (or simply $ ||\cdot|| $ when the context is clear). Inner products in $ L^{2}(D) $ are written $ (\cdot, \cdot)_{D} $ (or simply $ (\cdot, \cdot) $ when the context is clear). We will also denote pertinent duality pairings as $ \left\langle \cdot, \cdot \right\rangle _{X\times X^{\prime }} $, 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_{0}^{s}(D) $ denotes the closure of $ C_{0}^{\infty }(D) $ in the $ H^{s}(D) $-norm $ \Vert \cdot \Vert _{H^{s}(D)} $ or $ \Vert \cdot \Vert _{s} $.

2.3.   Energies and state space

The natural energy for linear beam dynamics is given by the sum of the potential and kinetic energies

Enforcing that $ \alpha = 0 $ for H and C configurations, and that $ \alpha \ge 0 $ for CF, we have that the dynamics evolve in the state space $ \mathcal H $, whose definition depends on the configuration:

where

Owing to the conditional structure of the state space (depending on the configuration and, for CF, the value of $ \alpha $) we opt for a compact notation, using

with corresponding inner-product $ (u, w)_H\equiv D(u_{xx}, w_{xx}) $ (equivalent to the standard $ H^2(0, L) $ inner product), and

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 $ \Pi $ term represents the non-dissipative and nonlinear portion of the energy:

As in the general case [12,Lemma 1.5.4], we can see that for any $ w \in H $, $ 0 < \eta\leq 2 $ and $ \epsilon > 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 $.

3.   Theory and discussion

3.1.   Solutions

We now discuss the pertinent notions of solution. We will distinguish between: (ⅰ) panel configurations ($ \mathbf C $ and $ \mathbf H $), for which we consider only $ \alpha = 0 $, and (ⅱ) the cantilever CF, where we admit $ \alpha \ge 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 (point-wise) interpretation of the evolution (2.1). For CF with $ \alpha > 0 $, boundary conditions are still interpreted weakly due to subtleties associated with inertial terms. (See ^[27] or ^[12].)

● Generalized solutions are $ C([0, T], \mathcal H) $ limits of strong solutions (when they exist). These are weak solutions ^[11,12], but they admit smooth approximation, and thus inherit some properties of strong solutions holding with respect to the topology of $ \mathcal 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 $ \mathcal H $ encapsulating the configuration and $ \alpha \ge 0 $.

Definition 3.1 (Weak Solutions). We say $ w \in H^1(0, T; L^2(0, L)) $ is a weak solution on $ [0, T] $ if:

and, for every $ \phi \in H $, we have

where $ d/dt $ is taken in the sense of $ \mathcal{D}'(0, T) $. Moreover, for any $ (\chi, \psi) \in \mathcal H $

Definition 3.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)\in L^\infty(0, T; H^4(0, L)\times H) $ and $ \frac{d^+}{dt^+}w_t $ is right-continuous.

More regular cantilever solutions are considered, but only for $ \alpha > 0 $.

Definition 3.3 (Cantilever strong solution). For the CF configuration, we say a weak solution to (2.2) with $ \alpha > 0 $ is strong if it possesses the regularity: $ w\in L^\infty(0, T; W) $, where

and $ w_t\in L^\infty(0, T; H^2_*(0, L)) $. In addition $ \frac{d^+}{dt^+}w_t $ is right-continuous and $ L^\infty(0, T; H^1_*(0, L)) $ with $ H^1_*(0, L) \equiv \{w \in H^1(0, L)\; :\; w(0) = 0\} $. (Such solutions satisfy the third-order boundary condition weakly—see ^[27].)

Finally, the most convenient notion of solutions will be those induced by a semigroup flow:

Definition 3.4 (Semigroup well-posedness). We say (2.1) or (2.2) is semigroup well-posed if there is a family of locally Lipschitz operators $ t\mapsto \mathcal S(t) $ on $ \mathcal H $, such that: For any $ y_0 = (w_0, w_1)\in \mathcal H $ the function $ t\mapsto S(t)y_0 $ (ⅰ) is in $ C([0, T]; \mathcal H) $, (ⅱ) is a weak solution to (2.1), and (ⅲ) is a strong $ C([0, T]; \mathcal{ H}) $ limit of strong solutions to (2.1) on every $ [0, T] $, $ T > 0 $. In particular, solutions are unique and depend continuously in $ C([0, T]; \mathcal H) $ on the initial data from $ \mathcal H $. Thus $ t\mapsto S(t)y_0 $ is a generalized solution. Furthermore, there exists some subset of $ \mathcal H $, denoted $ \mathcal D(\mathbb A) $ (with $ \mathbb A $ the nonlinear generator of $ \mathcal S(t) $), invariant under the flow, such that all solutions originating there are strong solutions.

3.2.   Well-posedness

For all results below $ D, k_0, \beta, U \ge 0 $. We employ the energetic notations from Section 2.3.

Theorem 3.5 (Linear semigroup well-posedness). For configurations C and H, the linear problem in (2.1) (with $ b_1 = b_2 = 0 $) is semigroup well-posed on $ \mathcal H_{\mathbf C} $ or $ \mathcal H_{\mathbf H} $ (resp.). Generalized (and strong) solutions obey the energy identity for all $ t \ge 0 $:

Similarly, for configuration CF with $ \alpha, k_1 \ge 0 $, the linear problem in (2.2) (with $ b_1 = b_2 = 0 $) is semigroup well-posed on $ \mathcal H_{\mathbf{CF}} $. Generalized (and strong) solutions obey the energy identity for all $ t \ge 0 $:

For such a linear problem, these results are well-established.

Theorem 3.6 (Nonlinear semigroup well-posedness). Let $ b_1 \in \mathbb R $ and $ b_2 > 0 $. For configurations C and H, the problem in (2.1) is semigroup well-posed on $ \mathcal H_{\mathbf C} $ or $ \mathcal H_{\mathbf H} $ (resp.). Generalized (and strong) solutions obey the energy identity for all $ t \ge 0 $:

Similarly, for configuratio CF with $ \alpha > 0 $ and $ k_1 \ge 0 $, the problem in (2.2) is (nonlinear) semigroup well-posed on $ \mathcal H_{\mathbf{CF}} $. Generalized (and strong) solutions obey the energy identity for all $ t \ge 0 $:

In the case of nonlinear clamped and hinged beams, the well-posedness has been known for sometime; consult the abstract theory in ^[12] and note the classical references (for both 1-D Krieger and 2-D Berger nonlinearities) ^{[2,3,5,15,22]}. In the case of the cantilevered beam CF with $ \alpha \ge 0 $, see the recent work ^[27] (which includes discussion of other cantilevered beams ^[34]).

The final case not addressed by the previous theorem is the nonlinear cantilever in the absence of rotational inertia. For this case ^[27] provides only an existence result.

Theorem 3.7. Let $ b_1 \in \mathbb R $ with $ b_2 \ge 0 $ and $ k_1 \ge 0 $. For (2.2) with $ \alpha = 0 $ in the CF configuration, weak solutions exist on $ [0, T] $ for any $ T > 0 $.

Without boundary damping (e.g., ^[34]), existence of strong solutions for CF with $ \alpha = 0 $ is open.

3.3.   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, $ \beta \ge 0 $ and $ U \ge 0 $:

Theorem 3.8. Let $ D > 0 $ with $ \beta > 0 $ and $ U \ge 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 \ge 0 $:

In the case of configuration CF, with $ \alpha > 0 $ (rotational inertia present) and $ k_1 > 0 $ (active square root damping), we have

Remark 5. Formally, CF trajectories for (2.2) with no inertia, $ \alpha = 0 $, do obey a global-in-time bound of the form (3.7), but the requisite multipliers (and associated integration by parts) are not permissible for the weak solutions guaranteed by Theorem 3.7.

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.5) and (3.6). They follow from the estimates in (2.6)–(2.7), taken together with a Lyapunov argument, and rely critically on the "good" structure of $ \Pi $, 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 $ \mathcal H $ topology), fully time-invariant set ^[12].

Theorem 3.9. Let $ D > 0 $ with $ \beta > 0 $ and $ U\ge 0 $. Also let $ b_2 > 0 $ with $ b_1 \in \mathbb R $.

For configurations C and H with $ k_0 \ge 0 $ the dynamical system $ (S(t), \mathcal H) $ generated by generalized solution to (2.1) has a smooth and finite dimensional compact global attractor.

For configurations CF with $ \alpha, k_1 > 0 $ and $ k_0 \ge 0 $, the dynamical system $ (S_{\alpha}(t), \mathcal H_{\mathbf{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)\in \mathcal H \cap \big(H^4(0, L) \times H\big) $, 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 attractor—see ^[12] 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 ^[12]. 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 ^[12]. 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 ^[26] (though earlier proofs have certainly appeared for these configurations). This work also discusses the application of a recent tool: the theory of quasi-stability ^[9,12]. The theory provides the existence of so called exponential attractors ^[9,12]. A generalized fractal exponential attractor for the dynamics $ (S(t), \mathcal H) $ is a forward invariant compact set $ A_{\text{exp}} \subset \mathcal 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_{\text{exp}} $ is finite in $ \mathcal H $.

Theorem 3.10. Let $ D > 0 $ with $ \beta > 0 $ and $ U\ge 0 $. Also let $ b_2 > 0 $ with $ b_1 \in \mathbb R $.

For configurations C and H with $ k_0 \ge 0 $ the dynamical system $ (S(t), \mathcal H) $ generated by generalized solution to (2.1) has a generalized fractal exponential attractor.

If the imposed damping is sufficiently large, i.e., $ k_0 > k_*(U, \beta, b_i, D, L, \beta) > 0 $, the dynamical system $ (S(t), \mathcal H) $ generated by generalized solutions to (2.1) has a fractal exponential attractor.

In the cantilever configuration CF, the result on the existence of a compact global attractor is much more recent ^[27].

4.   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.

We also take the static pressure $ p_0(x) \equiv 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 perform 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): (ⅰ) modal initial displacement, (ⅱ) polynomial initial displacement, and (ⅲ) elementary initial velocity. Let us denote $ \hat x = x/L $.

● [$ n $th Mode ID] For (ⅰ), 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\approx 4.6941 $ is the second Euler-Bernoulli cantilevered mode number and $ \mathcal C_2 \approx 1.0185 $.

● [Polynomial ID] For (ⅱ), with the C and CF 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 Ⅳ] For (ⅲ), 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 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.

5.   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, \beta $ for the flow, $ k_0, L $ for the beam—yield stable or unstable asymptotic-in-time dynamics. More specifically, our studies in this section are concerned with convergence to stationary states or unbounded growth of displacements as $ t\to \infty $ (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\equiv k_0+\beta $, which consists of both imposed damping $ k_0 \ge 0 $, as well as piston-theoretic (flow) damping $ \beta > 0 $.

The effects of the initial conditions are conservative (which is to say that, for $ \beta = 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 $ \beta w_t $. Additionally, with $ U > 0 $, the non-dissipative piston term $ -\beta Uw_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 $ \beta w_t $ bleeding energy out of the system, and the non-dissipative flow effect $ Uw_x $ inputting flow energy. 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 $ \beta $), a certain value $ U_{\text{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 $, $ \beta $ fixed), there exists a $ U_{\text{crit}}(k_0) $ such that for $ U < U_{\text{crit}} $, the transient dynamics are damped out, and for $ U > U_{\text{crit}} $ instability will still be observed. We can think of $ U_{\text{crit}} $ as an increasing function of the damping $ k_0 $ in the problem.

5.1.   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., ^{[16,17,21,28,42,43]}). Since the eigenfunctions of standard elasticity operators form 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.

5.1.1.   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 ($ \alpha = 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 $ (\cdot, \cdot)_H $, where $ H $ encodes the boundary conditions (5.1)).

For the in vacuo dynamics, we have: $ w_{tt}+D\partial_x^4 w = 0. $ Separating variables yields the dispersion relation: $ \kappa^4_n = \dfrac{\omega_n^2}{D}, $ where we have labeled temporal frequencies $ \omega_n $ and associated mode numbers $ \kappa_n $, $ n = 1, 2, .... $ The spatial problem is then $ \partial_x^4s+\kappa_n^4 s = 0 $, giving the fundamental set^¶:

^¶This choice of linear combinations is particularly convenient.

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 $ (\kappa_n, \omega_n) $.

For C and CF in Table 1, the mode numbers $ \kappa_nL $ are obtained by numerically solving the characteristic equation. We have

The $ c_n $ values are chosen to normalize the functions in the $ L^2(0, L) $ sense.

5.1.2.   Reduction to perturbed eigenvalue problem

Let us consider the Galerkin procedure for the full linear beam equation:

on $ (0, L) $, with boundary conditions according to the configuration being studied; recall that $ k = k_0+\beta $. We view the terms involving $ \beta, U, k_0 $ as perturbations and expand the solution via the in vacuo mode functions $ \{s_j\} $ as $ w(t, x) = \sum 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 (summing over $ m $):

Orthonormality of the eigenfunctions can be invoked to produce diagonal terms, whereas the terms scaled by $ \beta U(\partial_xs_m, s_n) $ are off-diagonal and give rise to the instability of the ODE system.

Remark 6. 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\times N $ system for the unknown coefficients $ \{q_j(t)\} $, yielding the approximate solution $ w_N(x, t) = \sum\limits_{j = 1}^Nq_j(t)s_j(x). $

To simply determine the stability of the problem as a function of the given parameters, we invoke a standard engineering ansatz ^[20,43] (and references therein): assume simple harmonic motion in time of the beam according to some dominant (perturbed) frequency $ \widetilde\omega $; we allow possible contribution from all functions $ s_n $ for $ n = 1, 2, ..., N $ via coefficients labeled $ \alpha_n $:

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 $ \widetilde{\omega} $. With the off-diagonal entries $ (\partial_x s_m, s_n) $ in hand (for $ 1 \le m, n \le N $ with $ m \neq n $), compute diagonal terms

we create the matrix for $ 1 \le n, m \le N $:

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, \beta, U, L $, we enforce the zero determinant condition for non-trivial solutions in $ \widetilde{\omega} $:

and solve for $ \widetilde \omega = [\widetilde \omega_{1}, ..., \widetilde \omega_N]^T $. The associated complex roots allows us to track the response of the natural modes to the perturbation terms—damping, $ (k_0+\beta)w_t $ and non-dissipative, $ \beta Uw_x $.

● If $ Im(\widetilde{\omega}_j) > 0 $, then we call the configuration unstable, and the associated linear dynamics should exhibit exponential growth in time with rate $ Im(\widetilde{\omega}_j) $ and frequency $ Re(\widetilde{\omega}_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 7. The emergence of an unstable perturbed eigenvalue (associated to $ Im(\widetilde \omega) $) 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 8. 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\ge 7 $ for C and CF, modal calculations are cumbersome and possibly inaccurate. Engineers have devised clever approximations of the mode functions to circumvent this issue ^[19], 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 $.

5.2.   Prediction of the onset of instabilities

We now predict the onset of instability for the model (5.1) 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_{\text{crit}} $, beyond which the system exhibits exponential growth in time.

In this section, we have adopted the non-dimensionalization procedure utilized in ^[42]. 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 \le L \le 500 $,

● $ \beta = 1.2\cdot 10^{-4} $ (allowed to vary in the range $ 10^{-5} \le \beta \le 10^{-2} $).

When employing finite difference simulations of (5.1), the task of identifying $ U_{\text{crit}} $ usually requires several batches of simulations, usually in some bracketing approach. In general this is not practical, as qualitative observations of flutter are somewhat subjective and may require excessive long-time simulations (for instance, accurately determining $ U_{\text{crit}} $ for $ L > 320 $). However, the modal approach outlined in Section 5.1.2 yields a direct, simple way to find $ U_{\text{crit}} $. For given parameter values, one can immediately analyze the perturbed eigenvalues associated with (5.5) to determine if an instability is present ($ Im(\tilde{\omega}_j) > 0 $ for some $ j $).

In Figure 1, we compare the results of a bracketing search for $ U_{\text{crit}} $ with the modal approach for a range of $ L $ values across all three boundary configurations. The points are the $ U_{\text{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_{\text{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 $ \omega_{\text{crit}} $ associated to each $ U_{\text{crit}} $ (not shown here), and the frequency of unstable oscillation is simply $ Re(\omega_{\text{crit}}) $. Finally, as the modal analysis is such a good predictor of $ U_{\text{crit}} $ for all configurations, the laborious finite difference simulations to search for $ U_{\text{crit}} $ in the clamped-free case for $ L > 320 $ were not performed.

In Figure 2, we determine $ U_{\text{crit}} $ as a function of $ \beta $, with $ 10^{-5}\le \beta\le 10^{-2} $. We observe monotonic linear decay in the $ U_{\text{crit}}(\beta) $ value as $ \beta $ increases up to around $ 10^{-3} $, and then somewhat erratic behavior for larger $ \beta $. Since $ \beta $ affects both the damping and destabilizing terms in the model (5.2), it is not obvious that its effects will remain monotonic on $ U_{\text{crit}} $ as $ \beta $ increases further. Indeed, an analysis of the dispersion relation associated to (5.2) indicates that a critical transition occurs at $ \beta = 1 $.

One could similarly extract critical values of $ L $ or $ \beta $ 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_{\text{crit}} $.

Another relationship of interest is that of the damping coefficient $ k = k_0+\beta $ on $ U_{\text{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_{\text{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_{\text{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_{\text{crit}} $ (i.e., $ U = 5 $ for $ L = 300 $, $ D = 23.9 $, $ \beta = 1.2\cdot 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) = \hat x^3(1-\hat x)^3 $, and the elementary initial velocity $ w_1(x) = \hat x (1-\hat x) $.

We note that, in the above, all dynamic simulations reflect that $ U > U_{\text{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, $ \beta w_t $, the effect of the initial condition decays away after a short amount of time.

6.   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_{\text{crit}} $, we will observe dynamic end behavior. We can then ask about the effect of any parameter (e.g., in-axis compression $ b_1 > 0 $) on the non-transient behavior.

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 $ \beta = 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:

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||^2 = \int_0^L|w_x|^2\, dx $.

6.1.   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 on the left we have the in vacuo ($ k_0 = \beta = b_1 = b_2 = U = 0 $) dynamics with the first H mode as an initial displacement (resulting in point-wise harmonic motion). To the right, we have a post-flutter—$ U > U_{\text{crit}} = 360 $ and $ b_2 = 1 $—sequence of snapshots that are demonstrative of a typical H fluttering beam (Note that for $ U > 0 $, the flow is from left to right).

In the Figures 6 and 7 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.

6.2.   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_{\text{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_{\text{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 = -\beta = -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_{\text{crit}} = 135.18 $. Energy profiles (log-scale) for various choices of $ U $ as multiples of $ U_{\text{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_{\text{crit}} $ while $ E(t) $ stays bounded for $ 0 < U < U_{\text{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_{\text{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_{\text{crit}} $ (as $ E(t)\to 0 $ when $ t\to\infty $)—each curve for $ U < U_{\text{crit}} $ has a linear profile (or envelope) with negative slope, indicating exponential decay. Energy values for $ U\ge U_{\text{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 $ \mathcal{E}(t) $ is shown. Note that $ \mathcal{E}(t) $ (and hence $ E(t) $, due to (2.7)) now remains bounded for all supercritical (unstable) velocities $ U > U_{\text{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_{\text{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 $ \mathcal{E}(t) $, Figure 11 shows both quantities for $ k_0 = 0 $, $ b_1 = 0 $, and $ b_2 = 1 $, with $ U = 2U_{\text{crit}} $.

6.3.   Blow up for nonlinear CF with linear boundary conditions

We recall that, for conservation of energy to hold (or, equivalently, to derive the equations of motion through Hamilton's principle) in the CF configuration, the boundary conditions (1.5) must be altered; in the panel configurations C and H, standard linear boundary conditions (2.1) are correct. (For more discussion, see ^[27].)

A natural question, then, for the nonlinear cantilever ($ b_2 > 0 $) is to ask about the stability of the (2.2) model in the presence of standard linear clamped-free boundary conditions (as in (6.1)). We emphasize that this combination is non-physical, as there is no conservation of energy associated to the dynamics. Below, in Figure 12, we achieve arbitrary growth (in time) of displacements (or energies) for finite difference solutions to:

for $ c $ sufficiently large.

Specifically, if $ c = 12 $, the dynamics exhibit periodic behavior; when $ c = 13 $, the dynamics grow with exponential rate. When the boundary conditions are adjusted appropriately (i.e., $ \partial_x^3w(1) = ||w_x||^2w_x(1) $), the nonlinear energy $ \mathcal E(t) $ is nearly perfectly conserved (Figure 12).

6.4.   Convergence to limit cycles

In this section we demonstrate some properties of trajectories that converge to a limit cycle oscillation, as described above.

Figure 13 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 decay to a nontrivial limit cycle behavior.

In Figure 14, 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.

6.5.   Uniformity of limit cycle

Next, we consider the CF configuration and provide tip profile at $ x = L = 1 $ for three different initial configurations:

● [2nd Mode ID] $ w(0, x) = s_2(x) = [\cos(\kappa_2x)-\cosh(\kappa_2x)]-\mathcal C_2[\sin(\kappa_2x)-\sinh(\kappa_2x)], \; w_t(0, x) = 0 $, where $ \kappa_2\approx 4.6941 $ is the second Euler-Bernoulli cantilevered mode number (with $ L = 1 $) and $ \mathcal C_2 = \left[\dfrac{\cos(\kappa_2)+\cosh(\kappa_2)}{\sin(\kappa_2)+\sinh(\kappa_2)}\right]\approx 1.0185 $;

● [Polynomial ID] $ w(0, x) = -4x^5+15x^4-20x^3+10x^2 $, $ w_t(0, x) = 0 $;

● [Linear Ⅳ] $ w(0, x) = 0 $, $ \; w_t(0, x) = x $.

First, Figure 15 represents the nonlinear ($ b_2 = 1 $), in vacuo ($ \beta = k_0 = 0 $) tip displacements with the three initial configurations described above. Figure 16 represents fluttering dynamics across the three different initial configurations.

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.

6.6.   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 17. 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 18.

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 roughly negatives of one another. In Figure 19 the midpoint displacement is plotted for these two cases.

6.7.   Non-simple limit cycle

It is possible to induce several different phenomena by manipulating parameter values for the nonlinear dynamics. In Figure 20 the midpoint displacement is shown for $ U = 5000 $, $ k = 101 $, $ b_1 = 5000 $, and $ b_2 = 5000 $. Note the initial transient dynamics are damped out quickly and the dynamics converges to a non-simple limit cycle.

6.8.   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 ^[18] that this setup is a candidate to produce a deterministic, chaotic dynamical system. (In the reference ^[18], the system (2.2) is reduced down to a two mode 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 (Figure 21).

First, from the four profiles corresponding to $ \epsilon \searrow 0 $, we see no discernible notion of convergence of the dynamics. Secondly, we see no periodic behavior corresponding to the initial configuration with $ \epsilon = 0 $. Additionally, it is not clear that the other configurations involving $ \epsilon > 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 9. We continued simulations for smaller values of $ \epsilon $, 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 22; namely a transient period followed by an energy "bump" up to a higher level. The transition is smooth.

6.9.   Effect of magnitude of $ b_2 $

The effect of increasing the nonlinear parameter $ b_2 $ was studied at a supercritical flow velocity ($ U > U_{\text{crit}} $) for the beam. In Figure 23, energy profiles for several different choices of $ b_2 $ are given for the parameter $ k_0 = 0 $ and $ U = 150 $ ($ U_{\text{crit}} = 135.97 $). Note that as $ b_2 $ increases, the energy plateau decreases (for an initial configuration fixed across all simulations).

7.   Synopsis of main numerical observations

Here we briefly outline the main observations in Sections 5 and 6 above, and succinctly connect to the theory in Section 3.

Modal approximation works well

In Section 5.1 we describe the "modal" approach used commonly in engineering to extract the flutter point (or point of instability, with respect to critical parameters). We demonstrate in Section 5.2 that—for physical parameters—the modal approach faithfully replicates the onset of instability with various parameters obtained empirically (through direct numerical simulations of the PDE dynamics). This approximation is very good across configurations, and across large ranges of relevant parameters.

Extensible nonlinearity controls non-conservative effects

In various simulations above we have confirmed what we established theoretically: the extensible nonlinearity featured here in (2.1) ($ b_2 > 0 $) is sufficient to control the non-dissipative (non-conservative) effects introduced by the flow through piston theory (1.2). This is independent of configuration, of course making the necessary change to the boundary condition in the CF configuration—see Section 6.3. When the damping due to piston-theory is taken into account, the overall bound on trajectories, as the transient dynamics decay, seems independent of the initial data (i.e., its size and type), and the magnitude of the non-transient behavior (typically an LCO) depends on the intrinsic parameters, such as $ b_2, U, \beta, L $, etc.

Onset of instability is linear

Through our simulations, we have observed that below the point of onset (for instance, when $ U < U_{\text{crit}} $ in a particular configuation) the dynamics are effectively not altered by the presence of the (higher order) semilinear effect of extensibility. This means that the emergence of instability (bifurcation) is not really affected by the value of $ b_2 $. In this way, we may claim that the onset of flutter is a purely linear phenomenon, regarding the manner in which $ \beta Uw_x $ perturbs the in vacuo eigenmodes of the system; however, to view the flutter dynamics without seeing energies grow exponentially in time (as per the previous point in this section), one must include nonlinear restoring force $ b_2 > 0 $.

Instability in each critical parameter

As can be expected, we can de/stabilize the dynamics through a myriad of parameters. Unsurprisingly, the flow parameters $ \beta $ and $ U $ are the principal candidates for the onset of instability. However, within certain limits, we can utilize $ L $ and $ b_1 $ (for C and H) to affect the onset of instability. On the other hand, if we impose damping of the form $ k_0 $, we can re-stabilize eigenvalues that have become unstable due to any of the aforementioned parameters. This is to say that damping of the form $ k_0w_t $ (for $ k_0 $ sufficiently large) is "strong" enough to stabilize flutter. For a fixed $ U > U_{\text{crit}} $, increasing $ k_0 $ provides two regimes: Initially, it increases the rate of convergence to the stable LCO, but, for $ k_0 $ sufficiently large, the flutter dynamics are eventually damped out; see Figure 14.

Relative stability across configurations

As noted from various plots above, principally Figure 1, the configurations have a relative hierarchy of stability: the clamped panel C is the most stable; the hinged panel H is less stable; the cantilever CF is least stable. This demonstrates the extent to which the boundary conditions (which serve to distinguish the in vacuo eigenfunctions) affect the inherent stability of each configuration in axial flow. It is important to emphasize, then, that across all possible configurations one could consider with the boundary conditions presented here (including all left/right combinations), the CF configuration is the easiest to make flutter.

Uniformity of end behavior

We have seen that, in most cases with no pre-stressing of the structure (when $ b_1 = 0 $), a very strong uniformity of end-behavior with respect to size and type (displacement, velocity) of initial data. This is true in the fully linear case, when trajectories are stable or when they are unstable, e.g., Figure 4; in the latter case, we note the same rate of exponential growth in time of trajectories independent of initial condition. In the nonlinear case—where interesting behaviors are found—we typically observe convergence to the same LCO for fixed parameters and varying initial conditions—see Figure 14.

Reproduction of qualitative behaviors

It is worth emphasizing here that this very simplified model (extensible beam with piston theoretic RHS) for a very complex phenomenon (axial flow flutter) can replicate all expected qualitative behaviors from the engineering point of view. This includes convergence to simple LCOs (Section 6.4), non-trivial LCOs (Figure 20), and non-trivial (buckled) steady states (Figures 17 and 18). We have even demonstrated a compelling case for deterministic chaos in Section 6.8, by playing large flow values $ U $ against large in-plane pre-stressing $ b_1 > > 0 $ for the hinged panel H.

Conflict of interest

The authors declare that there is no conflict of interest.

A.   Appendix: Fundamental frequencies