2.1. Experimental facility

Flag flutter experiments were performed in a suction-type wind tunnel at Auburn University, which has a 244 cm $\times$ 61 cm $\times$ 61 cm square test section (please refer to § 2.6 for flow characterization). At approximately three quarters of the test section’s length from its entrance, a sample holder spanned the test section vertically. The vertical configuration was chosen to minimize the effects of gravity through tensioning or bending. The sample holder consisted of two machined aluminium strips, the longer one of which was anchored to the test section at both ends whereas the shorter one was affixed to the longer one by means of eight flat-head screws. Together, the strips formed a streamlined body based on the SD8020-015-88 aerofoil with a chord length of 42 mm and a maximum thickness of 6 mm plus the thickness of whatever flag was clamped between them (figure 2). Spacers were used in the front part of the holder to keep both strips parallel. Through a tap test on the mast equipped with a 0.1 mm spacer, the lowest natural frequency and the associated damping ratio were computed to be 40 Hz and 0.011, respectively. The sample holder was sufficiently straight (Appendix A.1).

Figure 2. The sample holder with a sample in the wind tunnel. Picture mirrored for consistency with conventional flow direction from left to right: $L$ length; $H$ height.

2.2. Samples

Flag samples were cut from various materials, with the main direction of curvature from manufacturing aligned axially so as to maximize transverse flatness. All had a free length of $L = 148.3$ mm (5.84”) and a height of $H = 152.4$ mm (6”), which makes $H^*=1.03$ . Based on the resulting ratio of normal wall distance to flag length of 2.05 and the numerical results of Alben (Reference Alben2015), we do not expect the lateral test section walls to affect the stability and flutter characteristics of our flags. Based on Hiroaki & Watanabe (Reference Hiroaki and Watanabe2024), also the vertical separation of the flags to the top and bottom of the test section is sufficient to exclude confinement effects.

The samples were all laser printed with a chequerboard pattern of 4.064 mm (0.16”) pitch, which resulted in a grid of $37\times 38$ trackable corner features. With this pitch we expect to capture the range of spatial frequencies relevant for hysteresis, while the random uncertainty associated with stereo motion tracking does not cause too much noise in the spatial derivatives obtained by finite differencing. Table 1 lists the mechanical properties of all samples along with the number of data sets acquired. The plate bending stiffness of printed sheet stock was estimated from the tip displacement of strips cantilevered in a fixture, bending under their own weight. The fixture was flipped once to eliminate the influence of curvature inherent to the strips. Euler–Bernoulli beam theory provided the relationship between displacement and plate stiffness. By choosing appropriate strip lengths, care was taken to limit the tip slope. At rest, the tip of the flags deviated from the axial direction by angles of $0.5^\circ - 32^\circ$ due to inherent longitudinal curvature and gravity. The air flow eliminated any large curvature before the onset of oscillations, however, a possible impact on the results is discussed in § 3.4.

Table 1. Properties of flag samples: $h$ , thickness; $m$ , area density; $D$ , plate bending stiffness; $\mu$ , mass ratio; $H$ , flag height.

2.3. Motion tracking methods

Two complementary methods of motion tracking were employed, surface tracking from one side and edge tracking from below. For surface tracking, two VEO4K 990L cameras (Vision Research Inc.) were equipped with Nikkor 50 mm lenses set to an f-number of 16. They were vertically arranged as a stereo camera pair, the optical centres $36.23\pm 0.12$ cm apart and the optical axes converging at an angle of $35.6^\circ\pm 0.7^\circ$ . Vertical rather than horizontal arrangement was chosen for a better view during flutter. For edge tracking, a VEO 640L camera (Vision Research Inc.) was equipped with a Scheimpflug adapter and a Nikkor 50 mm lens set to an f-number of 2.8. The resolution was approximately 15 pixels mm^–1 and 10–13 pixels mm^–1 for surface tracking and for edge tracking, respectively.

Multiple high-power light-emitting diode (LED) arrays (MultiLED QX, GS Vitec GmbH) illuminated the flag for surface tracking. With transparent flags, a white surface behind the flag was illuminated instead. For edge tracking a single LED array was combined with two slits to illuminate only the lower edge. The LEDs were computer-controlled and turned on only during imaging to limit light-induced heating.

For calibration and data processing, MATLAB’s Computer Vision Toolbox was used. For the calibration of the stereo camera pair, laser-printed chequerboard patterns were applied to a precision-machined aluminium plate and waved in the volumetric region of interest (details in Appendix A.2). Also the edge tracking camera was calibrated by waving a chequerboard plate, and the plane that the flag’s edge would move in was calibrated with another chequerboard plate which was clamped to the sample holder.

The surface tracking images were processed by detecting the chequerboard pattern printed on all samples, using MATLAB’s gradient-based function detectCheckerboardPoints, and triangulation of corresponding features. The edge tracking images were processed with a MATLAB script that followed the path of highest image intensity. Local surface curvature was approximated by means of the first and second derivatives of the surface’s $y$ -coordinate along the tracking pattern’s principal directions. To reduce noise in figures of this document, the derivatives were evaluated on smoothing splines, whereas for Supplementary movies, they were computed directly. Gaussian curvature was computed as the product of both curvatures in this document, whereas for the movies, it was computed from the second fundamental form, which is independent of the coordinate system. All derivatives were computed via finite differencing. Uncertainties are detailed in Appendix A.2. To quantify a flag’s local bending stiffness, a polygonal cross-section with the flag’s known thickness was reconstructed. The smaller principal second area moment, $I$ , of the cross-section relevant for bending stiffness was calculated and normalized by $I_{\textit{flat}}$ , the value for a flat rectangular cross-section. A higher $I/I_{\textit{flat}}$ indicates greater bending resistance, with $I/I_{\textit{flat}} = 1$ denoting minimal stiffness.

2.4. Wake flow visualization

The beam of a continuous green laser (RN-532nm-G1000, Beijing Ranbond Technology Co., Ltd.) was expanded to a sheet through two planoconcave cylindrical lenses. The sheet was aligned with a horizontally mounted chequerboard calibration plate, which was placed approximately 6 mm above the test section’s centre line. As a tracer, fog from a fog machine (Z-800II, Antari Lighting & Effects Ltd.) operated with water-based liquid (Fog Worx, Sanco Industries, Inc.) was made to be ingested by the wind tunnel. A VEO 640L high speed camera (Vision Research Inc.) equipped with a Scheimpflug adapter and a 50 mm Nikkor lens (f-number 1.8) was pointed at the plate at an angle from below and calibrated by the same method as the edge tracking camera. The images were digitally reprojected to obtain a normal perspective. The observed area measured approximately 32 cm ( $\approx 2.2 L$ ) streamwise and 25 cm laterally, with a resolution of 5–7 pixels mm^–1. The flag’s tip was visible on the images in most conditions.

2.5. Instrumentation, data acquisition and signal generation

The sample holder was equipped with a set of four strain gauges at each end that chiefly responded to lateral bending. In combination with two DMD4059 strain gauge to direct current isolated transmitters (Omega Engineering Inc.) they recorded the frequency content of structural vibrations. This helped monitor resonance between the structural vibrations due to the fan motor’s drive (0–60 Hz), the sample holder and the sample. For the outcomes see Appendix A.1. The airspeed, $U$ , was measured with a pitot–static probe and a 1 INCH-G-4V pressure transducer (All Sensors Corporation). The tube was mounted 1.4 m upstream of the sample holder, 11 cm below the ceiling and 10 cm from the lateral wall. All data acquisition and signal generation was handled through a USB-6341 data acquisition device (National Instruments Corp.). The wind tunnel’s fan motor was powered by an adjustable speed drive (Toshiba H7), controlled by an analogue signal.

2.6. Flow characterization and experimental conditions

The air density was computed from pressure, temperature and relative humidity. The turbulence intensity at the location of the sample holder was characterized through hot-wire anemometry, with the sample holder removed but the pitot–static probe in place. Right on the centre line, 10 cm above it, and 10 cm below it, the turbulence intensity was 0.5 % at 1 m s^–1, 0.3 % at 5 m s^–1 and 0.25 % at 20 m s^–1. Table 2 informs about the geometric blockage of the test section. At the critical airspeeds (between 3.2 m s^–1 and 13 m s^–1), the Reynolds number based on the free flag length, $\textit{Re}_L$ , was in the range of 30 000 and 120 000, and between 8600 and 34 000 based on the chord length of the sample holder. From wake flow visualization during the non-fluttering state, evidently the flow just upstream of the flags’ tip was fully laminar up to $U=3.6\, \mathrm{m\,s^{-1}}$ . For $3.75\, \mathrm{m\,s^{-1}}\leqslant U\leqslant 6\, \mathrm{m\,s^{-1}}$ , patches of disturbed smoke were visible or one side was turbulent and the other side laminar. For $U\geqslant 7.5\, \mathrm{m\,s^{-1}}$ , the flow on both sides was turbulent. Examples are presented in figure 3.

Table 2. Blockage ratios caused by the sample holder and the fluttering flag from axial projection, at selected conditions.

Figure 3. Examples of wake visualization for flow characterization. (a) Fully laminar at an airspeed of $U=3.25\, \mathrm{m\,s^{-1}}$ (TerraSlate sample TE1), (b) fully laminar at $U=3.6\, \mathrm{m\,s^{-1}}$ (transparency sample TP2), (c) turbulent on one side at $U=4.35\, \mathrm{m\,s^{-1}}$ (copy paper sample CP1), (d) fully turbulent at $U=7.5\, \mathrm{m\,s^{-1}}$ (stencil sheet sample SS2). The vertical dark lines originate from scratches on the acrylic test section wall. Here $L$ is flag length. Sample identifiers refer to table 1.

Once a sample was clamped, the critical speeds were first explored through manual speed control. Then, an automated sequence of upsweep and downsweep with low acceleration rates determined the critical speeds more precisely. With those speeds, the stepwise velocity sweep for edge tracking was executed. At every step, we recorded a time series corresponding to roughly five flapping cycles based on the nominal flapping frequency determined before, with 60 images per cycle, which resulted in 600–1500 f.p.s. (frames per second). This was complemented by manually prescribing selected airspeeds for stereo tracking. We targeted recording a time series corresponding to at least five flapping cycles, with between 37 and 60 images per cycle, which resulted in 600–925 f.p.s. This test procedure was executed between one and three times for each sample, as indicated in table 1. The samples were removed from the sample holder between executions, and executions of the same sample were on different days. We refer to the data obtained from one execution as one data set. For at least one sample per material, also wake visualization was done at the same conditions as stereo tracking.

3.1. Flutter motion

Figure 4 illustrates typical motion tracking results. Most samples exhibited flutter in the lowest structural mode with only one neck (Eloy et al. Reference Eloy, Lagrange, Souilliez and Schouveiler2008), as observed in figure 4(a). In contrast, copy paper samples exhibited two-neck flutter, illustrated in figure 4(c), as anticipated from the mass ratios reported in table 1. The motion contains elements of both standing and travelling waves, with the tip moving along a figure-eight path, mostly pointing away from its travel direction, shown in figures 4(a) and 4(c). Figures 4(b) and 4(d) display local curvature envelopes, highlighting mode shape differences via the number of waists. Supplementary movies 1–4 provide visualizations of edge tracking data. Figure 4(e) presents a three-dimensional (3-D) reconstruction of an instantaneous surface with the lower edge locations superimposed for context. Our method reconstructs the entire surface, excluding the outermost thin rectangles. Supplementary movies 5 and 6 animate the 3-D tracking for samples CS1 and TE2, respectively.

Figure 4. Visualization of motion tracking data from (a–d) edge tracking and (e) from surface tracking. (a,c) Superposition of edge shapes with five emphasized instances and highlighted path of the free end. (b,d) Envelopes of local curvature computed from the data of (a,c). (e) Reconstruction of the surface in one instance and superposition of all locations of the lower edge. (a,b,e) Cardstock sample CS1, (c,d) copy paper sample CP2 (refer to table 1).

Figure 5 illustrates the tip amplitude variation through a stepwise velocity sweep. In the initial upsweep from $U=0\ \mathrm{m\,s^{-1}}$ , vibrations exhibit small amplitudes, sometimes appearing irregular or expressing a dominant frequency. For TE2, flutter initiates suddenly at a critical velocity, whereas SS1 shows small-amplitude LCO before flutter sets in, a pattern seen in most datasets. After transition to flutter, amplitudes change with increased velocity. During downsweep, amplitudes adhere to the upper hysteresis branch until flutter ceases. For TE2, amplitude changes are less smooth compared with SS1, possibly due to kinks that travel diagonally across the flag once per motion cycle, introducing irregularity in the motion (see Supplementary movies 2 and 6). Each data set displays a distinct amplitude plot shape (refer to Supplementary figure S1). Modes during small-amplitude LCOs differ markedly from flutter, featuring transverse bending, longitudinal twisting and asymmetry, except in SS1 and SS2, where similarity exists, yet accompanied by snap-through events (see Supplementary movie 7).

Figure 5. Evolution of tip amplitude throughout a stepwise velocity sweep of TerraSlate sample TE2 and stencil sheet sample SS1 (refer to table 1).

To further elucidate the flutter kinematics and form conjectures about how the flag exchanges energy with the flow, crucial for our ODE model, figure 6 illustrates the motion depicted in figure 4(a,c) for a half-cycle of flutter with superimposed velocity vectors. Drawing from Connell & Yue (Reference Connell and Yue2007), we consider the pressure field as a potential energy reservoir exchanged with the flag’s kinetic energy. Although we lack pressure data, the flag segment’s local curvature and travel direction permit qualitative insights into the adjacent pressure field. In one-neck flutter (figure 6 a–f), the tip’s figure-eight path correlates with a curvature reversal over a significant flag portion (figure 6 a–c) before motion reverses (figure 6 c,d). In two-neck flutter (figure 6 g–l), the first section of the flag often bends oppositely to the second section. Thus, in both modes, large portions of the flag predominantly move towards the convex side. If convex and concave areas experience negative and positive pressures, respectively (Thoma Reference Thoma1939), the pressure field mostly exerts positive work on the flag. For the purpose of modelling in § 4, we will assume that fluid forcing on the flag works all the time from one reversal to the other, and always in the direction of motion.

Figure 6. Kinematics of cardstock sample CS1 (a–f) and copy paper sample CP2 (g–l) over a half-cycle of flutter from edge tracking. The figure-eight tip path is indicated. Velocity vectors point into the local direction of motion. The vectors’ scale is normalized per data set. For sample properties, refer to table 1.

Taneda (Reference Taneda1968) noted that the transverse oscillations are phase-shifted at the flag’s rearmost sections, causing waves to move backward. Near the tip, an instantaneous upward deflection occurs with significant upward bending; during reversal (figure 6 c,i), the velocity is directed towards the flag’s concave side with a major upstream component. This indicates the flag loses energy to the flow as it opposes the pressure field. However, the flag may recover some energy from the pressure field on adjacent sections moving in the opposite direction (figure 6 c,i) and on the tip itself after reversal (figure 6 d,j). These observations will be revisited when introducing the ODE-based model in § 4.

3.2. Critical airspeed and repeatability

Figure 7 presents the reduced airspeeds for both flutter onset ( $U_c^*$ ) and stop ( $U_d^*$ ) across all data sets, which are approximately arranged by increasing mass ratio. Data were gathered from continuous and stepwise velocity sweeps, and 3-D tracking. Consistent with Eloy et al. (Reference Eloy, Kofman and Schouveiler2012), the stop speed (black) shows greater consistency across similar experiments than the onset speed (grey). The lighter copy paper samples, CP1, CP2 and CP3, exhibit a higher stop speed $U_d^*$ due to experiencing two-neck flutter.

Figure 7. Reduced airspeed $U^*$ of flutter onset ( $U_c^*$ , grey symbols) and of flutter stop ( $U_d^*$ , black symbols) for all samples. Thick vertical lines demarcate uncertainty from stepwise speed increments during edge tracking. Triangles and horizontal lines are from continuous velocity sweeps. Vertically stacked dots with connecting lines mark the intervals in which transition happened during 3-D tracking. Data sets are roughly sorted by increasing mass ratio. Data set labels refer to table 1. Two onset speeds from continuous sweeping marked with an asterisk had frequency matching between fan motor drive and flag vibrations (refer to Appendix A.1). Here $U$ is airspeed, $L$ flag length, $m$ area-specific mass, $D$ plate bending stiffness.

Most data sets exhibit consistent internal agreement across repetitions. However, copy paper and cardstock samples display notable variability across data sets, even when the same sample is tested. This variability is mainly attributed to humidity imbalances between the air and the sample, causing mild deformation during experiments. Factors such as minor clamping differences and possible plastic deformations during handling and storage may also contribute. Despite this, conclusions regarding the correlation between shape and flutter hysteresis behaviour remain valid since the 3-D flag shape is assessed each time. The highest repeatability across data sets was observed with the thickest plastic samples, SS1 and SS2.

3.3. Analysis of cross-sections throughout a velocity sweep

Figure 8 illustrates how the normalized tip amplitude, $A^* \equiv A/L$ , varies with the reduced airspeed, $U^* \equiv UL\sqrt {m/D}$ , during a velocity sweep employing edge tracking. The sample, CP2, exhibits an 82 % relative hysteresis loop extent, defined as $\varXi \equiv (U_c - U_d) / U_d$ . Locations (a) to (d) on this loop, close to transitions between static and fluttering phases, were captured using stereo motion tracking, with 230 frames per location at 920 f.p.s. – translating to six and three flapping periods for states (c) and (d), respectively. In state (c), due to instantaneously shallow view angles at the flapping tip, approximately 8 % of frames could not be reconstructed, but this is unlikely to significantly affect the findings. Figures 8(a) to 8(d) display statistics for the normalized stiffness, $I/I_{\textit{flat}}$ (defined in § 2.3), along the sample. The shaded grey area reflects the full time series’ envelope, and the lines depict example instantaneous distributions.

Figure 8. Reduced amplitude, $A^*$ , from edge tracking as a function of reduced airspeed, $U^*$ , throughout a velocity sweep. (a–d) Envelope (shaded area) of normalized cross-sectional stiffness, $I/I_{\textit{flat}}$ , from surface tracking and instantaneous distributions (lines) at three instances. Data of sample CP2 from table 1.

At point (a), no periodic motion is observed, only vibrations with varying amplitudes around 0.03 mm and approximately 45 Hz. The large transverse curvature, corresponding to a cross-sectional tip thickness of approximately 3 mm, results in significantly greater stiffness compared with a flat sheet of copy paper. While the sample holder enforces a straight clamped edge, the stiffness increases monotonically towards the free end. At point (b), aperiodic motion is noted but flutter remains suppressed. The wider envelope of $I/I_{\textit{flat}}$ signifies vibrations, yet retains the same shape as in (a). Point (c) is on the upper hysteresis branch, at the same velocity as (b). Here, flutter coexists with local flag flattening, allowing the envelope to approach 1.0 along much of its length, although flattening may occur momentarily only at selected points. For point (d), on the upper branch near the flutter stop, the envelope is higher than at (c), with local $I/I_{\textit{flat}}$ akin to the motionless state, (a). Complete flattening occurs except along the last 10 % of the length. Similar observations are applicable to all data sets except those lacking significant hysteresis or using stencil sheets (explained further below). These findings support the hypothesis that stiffness-inducing flatness defects are diminished during flutter but not on the lower hysteresis branch, likely contributing to hysteresis as explained by Eloy et al. (Reference Eloy, Kofman and Schouveiler2012). Supplementary movies 8 and 9 provide animations of points (b) and (d) of figure 8, respectively, also illustrating several cross-sections.

The stencil sheet samples with $h=0.26\, \mathrm{mm}$ , the thickest and heaviest in the study, exhibit a small-amplitude LCO with $A$ between $1.9-3.7\, \mathrm{mm}$ before flutter occurs (e.g. see sample SS1 in figure 5). This phenomenon appears unaffected by transverse curvature as the $I/I_{\textit{flat}}$ envelope reduces to 1.0 like in figure 8(c). This was observed consistently across all three data sets using this material (see table 1). It remains unclear if the same physics apply here as for other samples or if our $I/I_{\textit{flat}}$ criterion is inadequate. The sample holder’s finite rigidity is likely not a factor, given the flag root’s minor lateral vibration, approximately 0.1 mm in amplitude. We opted not to investigate this case further but note that Tang et al. (Reference Tang, Yamamoto and Dowell2003) observed a similar LCO state in an aluminium panel with a low height-to-thickness ratio of $H/h\approx 325$ and $H^*$ of 0.48 (for reference compare with our table 1), and suggested aerodynamic nonlinearities might cause distinct small- and large-amplitude LCO states. Conversely, Hiroaki & Watanabe (Reference Hiroaki and Watanabe2024) found no hysteresis for flags with $0.3\leqslant H^*\leqslant 1.17$ , $120\leqslant H/h \leqslant 470$ , $\mu \approx 2.9$ .

3.4. Correlation between flatness defects and onset delay

To explore the relationship between increased stiffness due to flatness defects and onset delay, we analyse two statistical metrics of $I/I_{\textit{flat}}$ in representative non-fluttering states (e.g. point (a) in figure 8), derived from time series surface tracking data.

1. (i) The median across all locations and instances, noted for its resistance to outliers.

2. (ii) The temporal median of the weighted sum of $I/I_{\textit{flat}}$ , using local curvature variation during flutter as weights. This variation is the envelope height of curvature, shown in figure 4(b,d). Although edge tracking data appear in the figure, here we use surface tracking.

The second metric assumes that stiffness of locations experiencing higher bending in each sample’s mode shape (two necks in copy paper, one neck in others) has more effect. Since non-fluttering states exhibit minor motion, using temporal median, average or minimum is not impactful. Figure 9 plots these metrics against $\varXi$ and $U_c^*-U_d^*$ as hysteresis metrics. Critical airspeeds are derived from 3-D or edge tracking (see figure 7), which explains high hysteresis uncertainties in some data points. All four combinations of these metrics indicate a positive correlation between the width of the hysteresis loop and the increase in stiffness due to transverse flatness defects.

Figure 9. Correlation of flatness defects and hysteresis. Two metrics of normalized cross-sectional stiffness, $I/I_{\textit{flat}}$ , in a representative non-fluttering state versus (a) $\varXi \equiv (U_c-U_d)/U_d$ and (b) reduced velocity difference, $U_c^*-U_d^*$ . Marker size indicates the deviation of a tangent to the flag tip from the wind tunnel axis in quiescent air, with some data points labelled for reference. Horizontal lines demarcate uncertainty in hysteresis (only one per data set). Two data sets marked with an asterisk had frequency matching between fan motor drive and flag vibrations before onset (refer to Appendix A.1). Here $U_c$ is onset speed, $U_d$ stop speed, $U^*$ reduced airspeed as defined in (1.1).

Figure 10. Straightening of an initially curved flag. (a) Quiescent air, (b) airspeed close to the speed where flutter stops (not fluttering). Above the graphs are the vertical projection of the shape of the flag’s upper edge (black) and lower edge (grey) with the ideally straight shape (dash–dotted). The left-hand end is clamped. The graphs are the envelope of normalized cross-sectional stiffness, $I/I_{\textit{flat}}$ , from surface tracking. Data of transparency sample TP2 (refer to table 1).

Four data points with low hysteresis but relatively high $I/I_{\textit{flat}}$ exhibit a flag shape with pronounced longitudinal curvature in still air, illustrated by larger markers in figure 9. Upon increasing the airspeed from $U=0\, \mathrm{m\,s^{-1}}$ to the non-fluttering state, the flag shape changes, as depicted in figure 10. While $I/I_{\textit{flat}}$ remains largely consistent in magnitude during this alignment in all four cases, the resulting S-shape might be related to low hysteresis. Alternative flatness metrics, potentially better than $I/I_{\textit{flat}}$ , could enhance the correlation by considering the entire 3-D surface rather than cross-sectional views.

Figure 11 illustrates how flatness defects affect similar samples by showing the tip amplitude evolution (above the coloured panels) for samples CS1 ( $\varXi \approx 0.02$ ) and CS2 ( $\varXi \approx 0.37$ ), along with their longitudinal, transverse and Gaussian curvatures. This is observed in two states: still near stop speed (figure 11 a,c) and during flutter at similar amplitudes (figure 11 b,d). In CS1’s ‘still’ state (figure 11 a), the tip vibrates with an amplitude of 0.2 mm at 14 Hz, while CS2 ( figure 11 c) shows aperiodic motion at 0.1 mm. During flutter, the longitudinal curvature greatly exceeds transverse curvature, displaying a smooth banded pattern. Consequently, Gaussian curvature – the product of the principal curvatures – is similar to the transverse curvature masked by multiplication with the longitudinal curvature. Here CS1 is relatively flat before flutter, shown by low transverse curvature (left-most column, central subpanel). The transition between states involves gradual amplitude changes and minor hysteresis. Conversely, CS2 exhibits significant transverse curvature (figure 11 c), which disappears during flutter except at the lateral edges (figure 11 d). Differences in transverse curvature account for differing onset speeds, and the reduction of curvature during flutter explains similar stop speeds. Supplementary movie 7 shows animations just before flutter onset of CS1 and CS2, and movie 5 depicts CS1’s hysteresis loop.

The motion tracking data showed small areas of fluctuating transverse and Gaussian curvature in the centre of the sample during flutter, likely from measurement error (noting the $0.3\ \mathrm{rad\,m^{-1}}$ uncertainty described in Appendix A.2). However, negative Gaussian curvature consistently appeared at the lateral edges when longitudinal bending was large, in all samples. This suggests transverse bending towards the large-scale convex side, possibly due to Poisson’s effect, where the concave side’s compression and convex side’s stretching are slightly offset by transverse bending at the edge. Therefore, it may not be completely accurate to claim that Gaussian curvature incurs a prohibitively high energetic cost in flag bending, as Eloy et al. (Reference Eloy, Kofman and Schouveiler2012) argued.

Figure 11. Reduced tip motion amplitude, $A^*\equiv A/L$ , of cardstock samples CS1 and CS2 as a function of reduced airspeed, $U^*$ (defined in 1.1) (above the coloured panels) together with metrics of instantaneous curvature for CS1 (a,b) and CS2 (c,d) in still (a,c) and flutter (b,d) states. The sample holder is located at the left-hand edge of each tile. From top to bottom are the vertical projection of the shape of the flag averaged along the span, longitudinal curvature, transverse curvature, Gaussian curvature. Here $A$ is amplitude, $L$ flag length. Sample identifiers refer to table 1.

3.5. Wake flow visualization

Figure 12. Flow visualization of wake structure. (a–c) copy paper sample CP1; (d–f) cardstock sample CS2. (a,d) Here $U$ just below onset speed, not fluttering; (b,e) nominally same $U$ as (a,d), fluttering; (c,f) $U$ just above stop speed, fluttering. Here $L$ is flag length, $U$ airspeed, $\textit{Str}_A\equiv 2fA/U$ , $\textit{Str}_L\equiv fL/U$ , $A$ tip amplitude, $f$ flutter frequency. Within each row, the instant (or phase angle) within the motion cycle is approximately the same. Sample identifiers refer to table 1. See also Supplementary movie 10.

Figure 12 presents a series of wake flow visualizations of a light copy paper sample (figure 12 a–c) and a heavier cardstock sample (figure 12 d–f). Figures 12(a,b) and 12(d,e) show airspeeds just below $U_c$ , contrasting the lower and upper hysteresis branches. Although frequency or $\textit{Str}_L \equiv fL/U$ is similar, the amplitude, characterized by $\textit{Str}_A \equiv 2fA/U$ , differs significantly. Figures 12(c) and 12(f) display flutter occurring at airspeeds slightly above $U_d$ , with figure 12(f) matching the conditions in figure 11(d). On the lower branch (figure 12 a,d), small-amplitude vibrations do not cause vorticity roll-up in the near wake, unlike during flutter where vortices roll up. The value of $\textit{Str}_A$ indicates the tip displacement to advection distance ratio per period of motion, with higher $\textit{Str}_A$ typically leading to tighter wake curves and discrete vortex formations. However, this is not always true across different mass ratios. For example, figure 12(c) shows a dominant vortex, while figure 12(e) with higher $\textit{Str}_A$ has multiple weak vortices, suggesting that the mode shape variance may play a role. These findings suggest the wake may contribute to hysteresis, as discrete vortices, thought to amplify flag motion, are absent in small-amplitude LCO (figure 12 a,d). However, during flutter close to $U_d$ , as in figure 12(f) no dominant vortex formation was observed with any but the lightest material, as in figure 12(c). Consequently, pressure feedback from the vortical wake is unlikely to sustain flutter in heavier flags, correlating with the assertion of Hiroaki & Watanabe (Reference Hiroaki and Watanabe2024) that vortex formation is less important as a source of aerodynamic nonlinearity that contributes to bistability. See Supplementary movie 10 for a moving version of figure 12.

As a secondary observation we note that sometimes the boundary layer state on both sides of the flag appears to alternate between laminar and turbulent, as in figure 12(e,f). Namely, the advancing side (above the flag tip in the figure) has laminar flow, while the retreating side (below the tip) is turbulent. This may be related to the pressure gradient along the flag, which is adverse on the retreating side which has relatively high pressure towards the tip and favourable on the advancing side which has relatively low pressure towards the tip.

5.1. The qualitative behaviour of a reference configuration

Figure 14. Response of the reference configuration of the model in terms of limit cycle amplitude, $A$ , (main panel) and corresponding time histories (a–h). Here $\tau$ is dimensionless time, $Y$ position of the oscillator, $U$ flow speed, $U_c$ critical speed of flutter onset, $U_d$ critical speed of flutter stop, $U_n$ single critical speed of the notchless oscillator. In the main panel, small dots mark the data points. When the notch term is (de-)activated, the amplitude follows the (grey dashed) black line.

Figure 14 illustrates the model’s reaction in the reference configuration to a stepwise increase and subsequent decrease in velocity. The labelled points on the central amplitude plot correspond to time-domain satellite plots. With $\chi$ at zero instead of 0.5, the amplitude response, indicated by the grey dashed line, shows no hysteresis.

At $U=0$ (point (a)), low amplitude oscillations at the forcing frequency $\omega _e$ (refer to (4.1)) are noted. As $U$ rises, oscillations at a frequency of approximately $\sqrt {(1+\chi )}$ emerge and grow, as seen in points (b) and (c). Until just before the critical flutter speed $U_c$ , the amplitude stays under $10^{-4}$ . Slightly exceeding $U_c$ , the amplitude begins a gradual increase until it matches $Y_{\textit{ref}}$ , where it escalates sharply, leading to significant flutter oscillations, shown in figure 14(d). With further increase in $U$ , the amplitude continues to increase, seen in figure 14(e). During downsweep, for $U\gt U_c$ , the amplitude mirrors the upsweep behaviour. For $U_d\lt U\lt U_c$ , it remains high without any discontinuity at $U_c$ , shown in figures 14(f) and 14(g). In this velocity range, the system is bistable, which manifests as hysteresis. Only when $U$ decreases below the stop speed $U_d$ , the model returns to the non-fluttering state with low-amplitude oscillations, seen in figure 14(h). During this transition, the amplitude decreases rapidly once it reaches $Y_{\textit{ref}}$ .

In the fluttering state, the amplitude increases roughly as $(U-U_n)^{0.6}$ (refer to figure 14 inset), where $U_n$ is the critical speed for a notchless oscillator, and the stiffening notch term has minimal effect (see the grey dashed line for comparison). This parallels the findings of Eloy et al. (Reference Eloy, Kofman and Schouveiler2012), though they used $\sqrt {U^*-U_c^*}$ (or maybe $\sqrt {U^*-U_d^*}$ as appears from their figure 2d and table 2) for curve fitting. Approaching $U_d$ on the upper branch, the amplitude is smaller compared with the notchless oscillator. As observed in experiments, transitions between states are marked by discrete amplitude jumps (figure 8). The amplitude plot resembles sample SS1 (figure 5), but in the ODE model, lower branch amplitudes do not reach macroscopic levels. Also, see Eloy et al. (Reference Eloy, Lagrange, Souilliez and Schouveiler2008), figure 2, for a similar experimental hysteresis loop. Note that with $\chi =0$ , a continuous amplitude increase is observed for $U\lt U_n$ , evident on a logarithmic scale with $A=2.2\times 10^{-5}$ at $U=0.28$ and $A=4.0\times 10^{-5}$ at $U=0.3$ . Similarly, with $\chi =0.5$ , the amplitude continuously rises for $U\lt U_c$ .

5.2. The mechanism of hysteresis

To elucidate the model’s behaviour (4.1), we analyse the energy balance between the work done by the fluid forcing term, $\textrm {sgn}(\dot {Y})U|Y|$ , and the dissipation term, $2\zeta \dot {Y}$ , over one oscillation cycle. The fluid forcing energy supplied per cycle is given by

(5.1) \begin{equation} E_{\textit{in}}={\int _{-A}^{A}U|Y|\,{\rm d}Y}+{\int _{A}^{-A}-U|Y|\,{\rm d}Y}=2UA^2, \end{equation}

where the negative sign in the second integrand arises because $\textrm {sgn}(\dot {Y})=-1$ . The energy dissipated by damping lacks a general analytical form due to the unknown non-harmonic shape of $Y$ , discussed further in § 5.4. Nonetheless, for small $\zeta$ , approximating $Y$ as harmonic with angular frequency $\omega$ , $Y=A\sin {(\omega \tau )}$ , is feasible. The dissipated energy, calculated from the damping force and velocity product over a period, is

(5.2) \begin{equation} E_d=\int _0^{2\pi /\omega }2\zeta \dot {Y}\dot {Y}{\rm d}\tau =2\pi \zeta \omega A^2. \end{equation}

These energy expressions are scaled by $m\omega _0^2$ , where $m$ is the dimensional mass.

Both $E_{\textit{in}}$ and $E_d$ scale with $A^2$ . However, $E_{\textit{in}}$ remains frequency-independent, whereas $E_d$ is proportional to frequency. This distinction clarifies the oscillator’s stable points. The frequency can be approximated by the natural frequency of the left-hand side of the ODE, considering $A\ll Y_{\textit{ref}}$ (first scenario) or $A\gg Y_{\textit{ref}}$ (second scenario). In the first scenario, oscillation is damped with a constant augmented stiffness of $1+\chi$ , having a non-damped natural frequency of $\sqrt {1+\chi }$ . Utilizing the ratio between damped and non-damped natural frequencies, $\sqrt {1-\zeta ^2}$ (Fowles Reference Fowles1986), the frequency becomes

(5.3) \begin{equation} \omega =\sqrt {(1+\chi )(1-\zeta ^2)}. \end{equation}

In the second scenario, the notch term is negligible as it influences only a minor portion of the cycle (see figure 14 for $U\gt U_c$ ), and a standard result for the Duffing equation (Jordan & Smith Reference Jordan and Smith2007) can be repurposed to obtain the natural frequency as

(5.4) \begin{equation} \omega =\sqrt {1-\frac {3}{4}\gamma A^2}. \end{equation}

Here, the damping term is neglected.

The critical onset speed, $U_c$ , can be deduced by balancing the energies, $E_d = E_{\textit{in}}$ ((5.1), (5.2)), along with the small-amplitude frequency formula for $\omega$ (5.3):

(5.5) \begin{equation} U_c = \pi \zeta \sqrt {(1+\chi )(1-\zeta ^2)}. \end{equation}

For $U\lt U_c$ , $E_{\textit{in}} - E_d$ remains negative, leading to damping of oscillations in absence of harmonic excitation, $\lim _{\tau \rightarrow \infty }A=0$ . As $\zeta$ increases, this $U_c$ estimate deviates from numerical results due to overlooked departure from harmonic motion caused by damping-fluid forcing interaction (explained in § 5.4). On the other hand, for $U\gt U_c$ , energy accumulates in the oscillator, resulting in large-amplitude oscillations reaching a limit cycle. The LCO amplitude can be estimated using these equations, employing the large-amplitude frequency approximation (5.4) in the energy balance, $E_d = E_{\textit{in}}$ ((5.1), (5.2)):

(5.6) \begin{equation} A = \sqrt {-\frac {4}{3\gamma }\left (\left (\frac {U}{\pi \zeta }\right )^2 - 1\right )}. \end{equation}

The ODE model suggests amplitudes can become very large if $U$ is sufficiently high, although in reality, the system’s amplitude is limited by the flag’s length and decreases beyond a particular speed, as depicted in figures 5 and 8 or in Eloy et al. (Reference Eloy, Kofman and Schouveiler2012).

For $\chi \gt 0$ , the transition from large to small $A$ as $U$ decreases below $U_d$ cannot be explained by straightforward analytical expressions due to the interplay between the notch term and the $\gamma$ -term, preventing estimations of $\omega$ and $E_d$ . As $A$ approaches $Y_{\textit{ref}}$ , the stiff region of the notch that accelerates the motion becomes more important, which causes an additional energy loss, which at $U_d$ leads to the transition to smaller $A$ . Equation (5.5) with $\chi =0$ , representing the critical speed of the notchless oscillator, can be used to approximate $U_d$ .

Figure 15. Illustration of the mechanism of bistability and hysteresis. Amplitude $A$ $(a_{1} -c_{1})$ and frequency $\omega$ ( $(a_{2} -c_{2})$ , from fast Fourier transform) as a function of velocity $U$ from a stepwise velocity sweep $(a_{1} ,a_{2})$ ; as a function of time $\tau /(2\pi )$ from onset transition $(b_{1}, b_{2})$ and from stop transition $(c_{1}, c_{2})$ . Solid lines with dots, numerical model; dashed lines, (5.6); dash–dotted lines, (5.3); dotted lines, (5.4) based on numerical amplitude. Theoretical onset velocity $U_c$ , (5.5).

To demonstrate the mechanism of bistability and hysteresis, figure 15 displays the reference configuration’s amplitude and frequency during a velocity sweep (figure 15 $a_{1,2}$ ) and transitions between states (figure 15 $b_{1,2}, c_{1,2}$ ). The analytical expressions (5.3) to (5.6) align fairly with numerical results. As indicated in figure 15( $a_{1}$ ), (5.6) slightly underestimates $A$ , due to the harmonic approximation overestimating energy loss, leading to a power equilibrium at lower amplitudes than the numerical model. While sweeping up from rest with $U\lt U_c$ , the frequency and dissipation rate remain constant (figure 15 $a_{1,2}$ ) and so high that amplitude growth is inhibited. At $U_c$ , amplitude increases as the fluid forcing provides excess energy, initiating flutter. During this transition (figure 15 $b_{1,2}$ ), frequency briefly drops while the amplitude rises sharply, then stabilizes; from here, frequency keeps growing with $U$ . During downward sweep, the frequency drops below that of small vibrations (possible because $A\gg Y_{\textit{ref}}$ ), enabling sustained flutter. Eventually, the frequency returns to its initial value when flutter ends. During this stop transition (figure 15 $c_{1,2}$ ), oscillations undergo several low-frequency cycles before a frequency increase, influenced by the notch term, causes a rapid amplitude drop, ending the transition.

The experimental data sets that allow for a frequency measurement just prior to flutter typically show an increase in frequency at onset, resembling the model’s behaviour, though occasionally a slight decrease is noted. When frequency can be measured over most of the lower hysteresis branch, it generally rises with $U$ (as opposed to being constant in the model). During downsweep, the frequency does not usually drop below the vibration frequency, except in four data sets (CP1, TE2, TP1, TP2) where it does not monotonically rise with $U$ . This represents a subtle distinction between the model and the physical system. The amplitude evolution during onset transition matches results from more complex numerical models (see figure 3a of Michelin et al. (Reference Michelin, Smith and Glover2008) and figure 7 of Accardo et al. (Reference Accardo, Eugeni, Mastroddi and Romano2013)). As $U$ approaches the critical speed, steady state is reached more slowly, as noted by Tang & Païdoussis (Reference Tang and Païdoussis2007) for their model. Fluctuations in $\omega$ in figure 15( $c_{2}$ ) for $100\lt \tau /(2\pi )$ are artefacts resulting from computing $\omega$ from zero crossings in presence of the excitation frequency $\omega _e$ (refer to figure 14 c). The fluctuations in $\omega$ seen in figure 15( $a_{2}$ ) are due to the limited frequency resolution of the fast Fourier transform in combination with a simple peak-finding method.

To summarize, the bistability for $U_d\lt U\lt U_c$ stems from the relative importance of the notch term within the instantaneous range $[-A,A]$ . On the lower hysteresis branch, where $A$ is small, the notch enhances stiffness across the period, increasing frequency and damping, which stabilizes the vibration. On the upper branch, where $A$ is large, primarily the hardening spring term $-\gamma Y^2$ influences frequency and hence the damping, while the notch’s effect on damping is insufficient to capture the oscillator. Interestingly, the notch does not directly cause bistability by stiffening but rather through enhanced damping. This indirect damping might also affect the physical flag by inducing an energy loss that inhibits flutter onset. However, our data does not show that frequency is lower on the upper branch or that otherwise alike samples with larger flatness defects vibrate at higher frequencies before onset, implying energy loss must have a different origin. Possibly, it may result from material or aerodynamic damping while the flag vibrates differently than its flutter mode, as noted in § 3.1. Similarly, the hardening term limits amplitude not by opposing motion directly but through frequency increases and hence dissipation that comes with hardening. This could happen in a fluttering flag, where increase in frequency and in strength of vortex shedding causes energy loss.

5.3. The effects of parameter variations

We explore the effects of parameter variations in figure 16, where critical speeds are plotted as a function of pairs of two parameters. Data points from the reference configuration are identified for orientation. In figure 16(a), the variation of the height of the stiffening notch, $\chi$ , is investigated for three damping ratios, $\zeta$ . As indicated in (5.5), $U_c$ increases with increasing $\chi$ or $\zeta$ and decreases inversely. In contrast, $U_d$ remains consistent across $\chi$ variations, whereas $\zeta$ influences it in approximately proportional amounts to $U_c$ . When $\chi =0$ , $U_c = U_d$ , showing no hysteresis. The influence of $\chi$ aligns with the experimental finding that stiffening a flag’s straight state via flatness defects raises $U_c$ , but negligibly affects $U_d$ (figures 7 and 9; see also Eloy et al. (Reference Eloy, Kofman and Schouveiler2012)). The rise in critical speeds with higher $\zeta$ mirrors the $U_c$ increase seen when the Reynolds number drops in stability studies and numerical models (Connell & Yue Reference Connell and Yue2007), as kinetic energy dissipation pushes an oscillating system towards a static state.

Figure 16. Onset speed, $U_c$ , and stop speed, $U_d$ , as a function of (a) notch height, $\chi$ , for multiple damping ratios, $\zeta$ , (b) notch half-width, $Y_{\textit{ref}}$ , for multiple values of the quadratic stiffening coefficient, $\gamma$ , (c) notch shape parameter, $\alpha$ , for multiple $\gamma$ . Based on the reference configuration defined in § 4. Here $U$ is flow speed.

In figure 16(b), the notch half-width $Y_{\textit{ref}}$ is varied across different hardening coefficients $\gamma$ . At low $|\gamma |$ , $U_c$ remains constant regardless of notch width alterations, but for larger $|\gamma |$ , $U_c$ increases as $Y_{\textit{ref}}$ is increased. Meanwhile, $U_d$ increases consistently with increasing $Y_{\textit{ref}}$ , and this is more pronounced with higher $|\gamma |$ . Large values for both parameters lead to the complete disappearance of hysteresis. These phenomena are explained by the shape of the stiffness term $1+\chi N(\alpha ,\beta ,Y/Y_{\textit{ref}})-\gamma Y^2$ . At lower $Y_{\textit{ref}}$ , the notch is much narrower than the valley formed by the hardening spring term, bounded vertically by $1+\chi$ (light grey curve in figure 17 a; see figure 13 b for labels). The notch width notably affects frequency only at small amplitudes, lightly influencing medium-amplitude oscillations just above $U_d$ , thereby altering $U_d$ itself. In contrast, when the width of the notch approaches that of the valley, the notch blends in with the valley, losing its distinct peak (black curve in figure 17 a). Depending on the shoulder slope from $\alpha$ , the notch’s top mirrors the underlying parabola $-\gamma Y^2$ , forming two local extrema that exceed $1+\chi$ (medium grey curve in figure 17 a). The more the notch blends in and loses its peak, the more the effects on onset and stop become indistinguishable, eventually eliminating hysteresis. The extent to which these mechanisms occur in physical flag flutter remains uncertain; for instance, $Y_{\textit{ref}}$ might relate to the wavenumber of transverse curvature defects in a flag.

Figure 17. Effects of large notch half-width, $Y_{\textit{ref}}$ , and/or large notch shape parameter, $\alpha$ . (a) Stiffness term of (4.1) for three parameter sets. (b) Limit cycle amplitude of a velocity sweep with $\alpha =64$ and $Y_{\textit{ref}}=0.6$ whereby $U$ reverses once before large-scale motion sets on, and once after. Here $Y$ is oscillator position, $A$ is the amplitude of $Y$ , $U$ is flow speed. Refer to § 4 for definitions.

In figure 16(c), the notch shape parameter $\alpha$ is modified for different $\gamma$ values. Changes in $\gamma$ mainly impact $U_d$ while $U_c$ remains unaffected, as $\gamma$ is absent from (5.5). However, $U_c$ can be influenced if $\alpha$ is large enough, leading to two peaks in the stiffness term, which increase the onset speed. Here, the amplitude grows to macroscopic levels of the order of $Y_{\textit{ref}}$ as the oscillator traverses the dimple at the top of the notch (see figure 17 b) before transitioning to large-amplitude motion. Lowering $U$ prior to this transition avoids hysteresis (dashed line). This is analogous to wind tunnel findings, where two flutter states can occur, showing a small-amplitude LCO before (large-amplitude) flutter onset, as observed in stencil sheet samples (e.g. SS1 in figure 5) and in the plate flutter case of Tang et al. (Reference Tang, Yamamoto and Dowell2003).

5.4. Oscillation dynamics

Figure 18. Phase portraits of (a) the ODE model (reference configuration of § 4, flow speed $U=0.4$ ), a van der Pol oscillator (nonlinear damping parameter $\varepsilon =0.2$ ), and copy paper sample CP3 (refer to table 1) at a flow speed slightly beyond the critical speed ( $U=6.61\ \mathrm{m\,s^{-1}}$ , reduced airspeed $U^*=19.3$ ), (b) the ODE model with different amounts of damping, $\zeta$ , at a flow speed slightly beyond the critical speed (notch height $\chi =0$ , quadratic stiffening coefficient $\gamma =-0.3$ , excitation amplitude $E=0$ ). The circular path of harmonic oscillation is given for reference. For normalization, the $Y$ -amplitude was divided by the peak amplitude, $A$ , and the velocity, $\dot {Y}$ , was divided by $A\omega$ . Here $Y$ is oscillator position, $\omega$ angular frequency. Note that progressive time corresponds to clockwise rotation.

Following our prior analysis focused on amplitudes and frequencies, we conclude this section by contrasting the ODE model’s motion, characterized in terms of displacement, $Y$ , and velocity, $\dot {Y}$ , with experimental data and other common oscillator models. This complimentary time-domain viewpoint enhances the understanding of the model’s behaviour and highlights both its similarities and differences to other oscillators. Figure 18(a) showcases four phase portraits.

1. (i) Harmonic motion serving as a baseline, corresponding to a sinusoidal waveform with a circular phase portrait.

2. (ii) The ODE model in reference configuration at $U=0.4$ , exhibiting non-harmonic motion due to non-constant stiffness and external forcing (refer to § 5.2).

3. (iii) A weakly nonlinear van der Pol oscillator characterized by the equation

   (5.7) \begin{equation} \ddot {Y}-\varepsilon \big(1-Y^2\big)\dot {Y}+Y=0. \end{equation}
   Here, damping is contingent upon the state variable $Y$ , switching between energy dissipation and addition for $|Y|\gt 1$ and $|Y|\lt 1$ , respectively. This oscillator is utilized in VIV wake models (Facchinetti et al. Reference Facchinetti, de Langre and Biolley2004). Those models were an inspiration for our own model, yet this oscillator is not directly applicable to flag flutter due to its self-exciting nature. We use $\varepsilon =0.2$ to match the experimental phase portrait.

4. (iv) The transverse tip displacement during sample CP3’s flutter. When other locations along the rear-half of the flag are chosen instead, the shape of the curve is similar and a little more circular at some of them.

In the experimental motion cycle, $\dot {Y}$ peaks after $Y$ crosses zero. The weakly nonlinear van der Pol oscillator behaves qualitatively similar, whereas in the ODE model, peak velocity occurs before $Y$ crosses zero. While both the van der Pol and ODE models act as relaxation oscillators, a notable distinction is that the van der Pol oscillator’s energy input peaks when $-1\lt Y\lt 1$ , while damping reduces motion for $|Y|\gt 1$ . Conversely, the ODE model receives energy at the reversal point, with constant damping. The fluid forcing term $\textrm {sgn}(\dot {Y})U|Y|$ , if moved to the left-hand side of (4.1), can be interpreted as a spring term, adding $UY^2$ in potential energy at every sign change in $\dot {Y}$ . Despite these differences, all motions remain close to harmonic, and the motion patterns of the ODE model and the physical flag can still be considered analogous. For future work, the fluid forcing term in the ODE model could be adjusted to better match the experimental phase relationship.

Figure 18(b) depicts how increased damping affects the system. The peak velocity occurs before $Y$ crosses zero, regardless of damping $\zeta$ . A higher $\zeta$ results in a departure from harmonic motion, with an asymmetry where the initial acceleration is followed by an increasingly slow phase. For large $\zeta$ , the motion transitions to a slow, creeping behaviour, with kinetic energy considerably lower than potential energy (note that due to the normalization by $A\omega$ , this is not directly apparent in the phase portraits). Thus, with large $\zeta$ , the model fails to accurately describe the motion of a fluttering flag.