Nonlinear modelling of aerodynamic self-excited forces: An experimental study Journal of Wind Engineering & Industrial Aerodynamics

The bridge aerodynamics research community is currently discussing several nonlinear wind load models for bridge decks, but no de ﬁ nite conclusion on which model is superior to the others is currently available. In this paper, we use experimental data for a double-deck section model tested in an advanced forced vibration rig to study the observed nonlinearities and to gain insight into what characteristics the nonlinear load model should be capable of modelling. Single harmonic horizontal, vertical and pitching motion; combined motion; and stochastic motion are considered. This approach allows the investigation of a more extensive range of nonlinear behaviours than regular wind tunnel testing. The typical nonlinear characteristics observed are mean drift, deviation from superposition and harmonic distortion. Further, we introduce a simple response-surface model for force prediction using polynomial combinations of the inputs and its derivatives. The model helps to gain further insight into the nonlinearity of the problem at hand and to select which re ﬁ ned modelling approach can be used in future work.


Introduction
The Norwegian Public Roads Administration (NPRA) plans to build a ferry free road on the west coast of Norway that will stretch from Kristiansand in the south to Trondheim in the north. The plans include several long-span bridges, since there are many wide and deep fjords along the route. The Sulaford is one of the significant challenges, and the plans include a suspension bridge with a 2800 m main span as one of the alternatives to cross the fjord (Multiconsult and Brancaleoni, 2015).
Multiple authors have suggested multi-deck cross-sections for longspan suspension bridges due to their superior aerodynamic performance. (Ogawa et al., 2002;Sato et al., 2000Sato et al., , 2002. (G. Diana, Fiammenghi, et al., 2013) states that for bridge spans longer than 2000 m, double or triple bridge deck sections are necessary (Larsen, 2008). reported that twin-decks outperform single deck cross-sections in terms of flutter, but exhibit more energetic vortex-induced vibrations, which is also supported by (W. L. Chen et al., 2014;Kwok et al., 2012). Closed box single deck girders have been the subjects of extensive research over the last three decades, leading to well-established linear load models (C. G Bucher, 1987;Caracoglia and Jones, 2003;Costa and Borri, 2006;Nowicki and Flaga, 2011;Øiseth et al., 2010;Sarkar et al., 1992Sarkar et al., , 2009). For the aeroelasticity, the state-of-the-art linear model is the aerodynamic derivatives model. The development of the theory stems from Theodorsen's model for an ideal flat plate in a potential flow abiding the Kutta condition (Theodorsen, 1934). (Scanland and Tomko, 1971) generalized the idea to include bluff bodies. Other relevant models that are widely accepted are the linearised quasi-steady theory, and the full quasi-steady theory (Kovacs et al., 1990;Miyata and Sato, 1995). Scanlan's model has been tested for multi-box girders (Andersen et al., 2016;G. Diana et al., 2004;Yang et al., 2015), but both CFD simulations and wind tunnel tests show that the multi-deck is prone to nonlinear aerodynamic behaviour (G. Diana et al., 2004;Skyvulstad et al., 2017;Zhang et al., 2017;Zhou et al., 2019;Zhou et al., 2018). Double-deck sections can also be more prone to exhibit nonlinear behaviour than single-deck sections (Skyvulstad et al., 2017). The nonlinearities of the multi-deck cross-sections might stem from the complex fluid flow due to the flow separation on the upstream deck and the highly turbulent flow hitting the downstream deck and its interaction, giving higher-order fluid memory in the system. (Argentini et al., 2012;Kwok et al., 2012;Xiang and Ge, 2002).
Identification of nonlinear self-excited forces for single-deck crosssections has been experimentally investigated in (Guangzhong Gao et al., 2018) for single degree torsional flutter, utilizing a free-vibration rig. The centre part of the sectional model was a force balance allowing extracting accurate self-excited forces during free vibration (Guangzhong Gao, Zhu, Wang, et al., 2020). further tested another single deck cross-section that was susceptible to two-degree of freedom flutter with post flutter limit cycle oscillations (Zeng shun Chen et al., 2017Chen et al., , 2018. used a forced vibration rig to identify nonlinear unsteady galloping forces on prisms and rectangular cylinders.
There exist several nonlinear wind load models for bridge decks. The choice of model depends on the problem at hand, what types of nonlinear properties the system exhibits, whether or not fluid memory is included, which parameters are important, what type of experimental data exists, how much experimental data is needed, and more. The ease of use for bridge engineers and the physical interpretations of the model parameters are also relevant. Multiple nonlinear load models for self-excited and buffeting forces have been suggested in the bridge aerodynamics research community.
(Giorgio Diana et al., 1993) proposed a corrected quasi-steady theory that implements the fluid memory effect into a nonlinear framework of the quasi-steady theory (G. Diana et al., 2008). suggest modelling nonlinear aerodynamic loads as a function of the effective angle of attack, which depends on the structural motion and the turbulence. The researchers used rheological models to model nonlinearities. In later research (G. Diana et al., 2010), used a polynomial together with the rheological model to model the nonlinearities (G. Diana, Rocchi, et al., 2013). stated that merging all the contributions into one unique variable was not possible for the multi-deck cross-section considered. Not strictly relying on the angle of attack adds flexibility to the model, but more coefficients need to be determined.
Some models separate the high-and low-frequency part of the turbulence. Linear models are often used for the high-frequency component while considering the low-frequency component as a slowly varying change of the mean angle of attack. The superposition of the two parts provides the total load. These models follow the findings from (Bocciolone et al., 1992), which show that there is low-frequency turbulence at the bridge site that creates a slow but considerable variation in the effective angle of attack (Fenerci and Øiseth, 2018). also reports such an effect. Several contributions suggest using nonlinear quasi-steady models without memory for the low-frequency part and using standard linear models based on rational or indicial functions for the high-frequency part (X. Chen and Kareem, 2001;G. Diana et al., 1995). There are also hybrid (X. Chen and Kareem, 2003), modified hybrid (Wu and Kareem, 2013b), and improved band superposition (G. Diana, Rocchi, et al., 2013) models that utilize similar principles of high-and low-frequency separation.
The Volterra series is another more general nonlinear model that utilizes the summation of higher-order convolutions to include the higher-order memory of the system (Schetzen, 1980;Volterra, 1959). In the frequency domain, the Volterra series is expressed as generalized frequency response functions. Some authors have used the load model for bridge aerodynamics in the frequency domain (Carassale et al., 2014;Carassale and Kareem, 2010) and time-domain (Deno€ el and Carassale, 2015;Wu & Kareem, 2013a, 2013c, 2015a, 2015b. In the contributions mentioned above, the data for calibrating and validating the models mainly come from CFD and analytical expressions, while the use of experimental data is limited. (Wu and Kareem, 2011) suggests using artificial neural networks (ANN) to model nonlinear aerodynamic forces. They used variables related to motions and turbulence as independent inputs into the network without combining them into an effective angle of attack. (Zhou et al., 2018;Zhou et al., 2019;Zhou et al., 2019) suggested the nonlinear aerodynamic force model (NAFM) and the general nonlinear aerodynamic force model (GNAFM). The models were calibrated using data from CFD simulations of bridge decks in harmonic motion. Data from free vibration tests and full aeroelastic model tests were used to validate the models. The model consists of a system of nonlinear differential equations that take static, motion-induced, buffeting, and unsteady effects into account.
(Guang zhong Gao and Zhu, 2017; Guangzhong Gao and Zhu, 2015) Developed a nonlinear self-excited force model utilizing a polynomial combination of the inputs and an energy equivalent principle on the specific problem to reduce the number of necessary polynomial terms. This paper presents a detailed investigation of the nonlinear aerodynamic characteristics of a double-deck section with curved undersides. The section exhibits nonlinear aerodynamic behaviour at relatively small motion amplitudes, which makes it an interesting case study for nonlinear force modelling. As outlined above, there are many nonlinear aerodynamic load models for bridge decks. However, experimental data for validating the models for general motion are rare and experimental data from combined vertical and pitching motions are also very scarce. CFD simulations, analytical expressions of idealized models, or wind tunnel tests that apply single harmonic motion are thus the basis for the presented modelling approaches. To gain further knowledge based on experimental data, we have used a single harmonic motion as well as stochastic combined vertical and pitching motions. The experimental campaign gives insight into the nonlinear characteristics that are important for the advanced models to capture and also sheds some light on which effects need to be included in the modelling. A simple response surface model for modelling of nonlinear aerodynamic forces is also proposed and validated against the experimental data, giving more insight into the nonlinearity at hand.

Nonlinear models
In this chapter, some nonlinear models are reviewed. These models typically have no/weak memory with fairly basic input parameter combinations. For a more in-depth overview of nonlinear models, refer to the cited paper in the introduction.

Quasi-steady theory and linearised quasi-steady theory
It is common to use the quasi-steady theory for nonlinear modelling of aerodynamic forces for bridge decks. Fig. 1 illustrates the central concept and also shows the variables involved. Here, Vis mean wind speed, u, and w represent the horizontal and vertical turbulence components. V rel represents the relative wind speed; βis the wind angle of attack; while B and D symbolize the width and height of the cross-section, respectively. Configuration (a) in the figure is the static equilibrium condition. Configuration (b) is the instantaneous displaced configuration defined by r n , where n 2 fx; z; θg denotes displacements in the horizontal and vertical direction, as well as the pitching rotation. The overbar denotes the mean value, and the dot indicates the derivative with respect to time. F n , where n 2 fx;z;θ;L;D;Mg, denotes the forces in the global directions x, z, and θ and the local relative wind directions L, D, and M. The quasi-steady theory assumes that the unsteady effects are negligible. This assumption means that the forces acting on the cross-section are only dependent on the instantaneous angle of attack and the relative velocity and that it is, thus, not necessary to consider the motion history. These assumptions are fulfilled for cases when the time it takes for a parcel of air to pass the section is short compared to the oscillation period of the motion. One can then use the static coefficients as a function of the angle of attack to determine the wind load. A much-used linearisation of the quasi-steady theory is deducted by linearising the static coefficients around a mean angle of attack and neglecting the higher-order terms of turbulence and motion as follows: (1) (4) where α and α f denote the mean and fluctuating part of the angle of incidence, respectively. The full quasi-steady theory is a nonlinear description without memory that utilizes the static coefficients as the basis for the nonlinearity. The formulation also includes cross-terms because the expanded expressions include terms such as r z r θ ,u _ r x and similar terms. The force coefficients are also functions of the angle of attack α, which enhance the nonlinear modelling capabilities.

Hysteretic load models
The quasi-steady theory does not include the time derivative of the pitching motion, which is a severe weakness since aerodynamic damping for the pitching motion is essential. It is possible to include this effect by using another point as a reference point. This has, for instance, been done by (G. Diana et al., 2008), where the following definition of the angle of attack is used.
Here, ψ is the dynamic angle of attack, and B* is a constant that may be different for each force component, but which often is taken as the width of the section. Using the dynamic angle of attack simplifies the relationship between the motion, turbulence, and load, as follows: Different nonlinear functions are applied to map the effective angle of attack to the wind-induced forces (G. Diana et al., 2008). suggests using a rheological model (Wu and Kareem, 2013b), used a polynomial, and (G. Diana et al., 2010) applies a combination of the polynomial and rheological models. The equation below presents a polynomial model.
2 þ a 6 ψ 3 þ a 7 ψ 2 _ ψ::: The hysteretic load model with polynomial mapping is a nonlinear load model with a weak memory since both the input variables and their first and sometimes the second derivatives with respect to time are included in the modelling. This provides the opportunity to model a phase lag between the inputs and forces and nonlinear hysteresis loops.
The formulation also includes higher-order cross-terms since expanding the expression above reveals products of the inputs, such as _ r z _ r θ , _ r z _ r 2 θ and similar products. (Wu and Kareem, 2011) use an artificial neural network as a nonlinear mapping of the inputs to the wind forces. The authors included a 2DOF dynamic system into the neural network, which makes the network predict the acceleration at t nþ1 . They applied a bipolar sigmoidal network up to the third-order, including all cross-terms. The acceleration of the vertical motion can be written as follows: € r nþ1 z ¼ ℚ L € r n z ; _ r n z ; r n z ; € r n θ ; _ r n θ ; r n θ ; V n ; u n ; w n ; V nþ1 ; u nþ1 ; w nþ1

Artificial neural network
where ℚ L is the nonlinear function described by the neural network that maps the inputs to the forces. The neural network model also has some memory since the variables and their time derivatives are used as inputs.
The cross-terms are directly implemented in the model up to the order of the neural network. The summary presented above, together with the nonlinear load models presented in the introduction illustrates that there are a vast number of nonlinear models available in the literature and that the models have somewhat different modelling capabilities. The linear models for bridge aerodynamics are well established, and there are standardized procedures for wind tunnel testing and the modelling of the forces. Nonlinear modelling is, however, a much more challenging field since deciding which models to use and what type of nonlinear function one should use to model the forces are not straightforward decisions. This paper, therefore, provides experimental results that might shed some light on which terms are important to include in the modelling and which might be disregarded.

Nonlinear self-excited force model
(Guang zhong Gao and Zhu, 2017;Guangzhong Gao and Zhu, 2015) Developed an energy equivalent principle (EEP) for identifying the nonlinear self-excited forces. The general form of the force model is a polynomial combination of the inputs, similar to the polynomial function shown in equation (7), but with motion/velocity input and not the dynamic angle of attack. The model is used in several works, and often specialized to handle a single type of limit cycle oscillation event such as VIV-Galloping(Guangzhong Gao, Zhu, Li, et al., 2020), galloping instability (Zengshun Chen et al., 2020) and torsional flutter (Guangzhong  Gao et al., 2018). The specialization utilizes the energy equivalent principle to find the elements of the model who contribute to the important energy balance of the system, disregarding the other elements that only contribute to the direct force or stiffness. This approach simplifies the identification problem significantly. An example of the specialized model for nonlinear self-excited pitching moment of torsional flutter LCO is given as (Guangzhong Gao et al., 2018): Where, the A's, are unknown reduced frequency-dependent coefficients. The NSEF-model, in its full form, is a general model with weak memory due to inclusion of velocity and displacement terms. In the reduced form shown, memory is also included.

Experimental setup
The wind tunnel experiments were carried out in the Fluid Mechanics Laboratory at the Norwegian University of Science and Technology (NTNU). The wind tunnel is a closed-circuit, single fan wind tunnel with a maximum wind speed of 30 m/s. The test section is approximately 11 m long, 2 m high, and 2.7 m wide, and the flow is nearly laminar with a uniform velocity profile at the inlet and with a turbulence intensity of approximately 0.2% (Adaramola, M. S. Krogstad, 2009). A pitot probe 4 m upstream of the cross-section measures the mean wind speed. Forced vibration tests or free vibrations tests are commonly used to characterize the aerodynamic properties of bridge decks. A particular advantage with the forced vibration tests is that it is straightforward to directly measure the forces that act on the section model for any applied motion. It is, however, a disadvantage that there is no proper fluid-structure interaction since the movement is prescribed and not affected by the aerodynamic forces. Free vibration tests capture the fluid-structure interaction effect, but it is more difficult to determine the aerodynamic forces. Forced vibration tests are used to obtain the results presented in this paper.
The forced vibration rig consists of two 3DOF actuators on each side of the wind tunnel. The actuators consist of two linear motion slides driven by ball screws that move the section model both vertically and horizontally. Torsional motion is ensured by zero-backlash shafts couplings connected to the ball screw and a servo motor with a planetary gear. The 6 servo motors are controlled by a multi-axis modular control system (MC4U from ASC Motion Control) support the section model. The vertical and horizontal axes can travel AE10 cm, and the torsional axes can travel AE90 . In-between the motors and the cross-section are two 6DOF Gamma load cells from ATI Industrial Automation, which enable the high precision force and torque measurements. A full description of the setup is given in (Siedziako, 2018;Siedziako et al., 2017). Fig. 2, Fig. 3 and Fig. 4 shows photos/illustrations of the rig. The rig can force the section model to move in an arbitrary direction in 3 degrees of freedom, with high precision. This makes the rig well suited for investigating the nonlinear properties of the bridge cross-section by applying an advanced bridge motion. Fig. 2 shows a drawing of a 1:50 scale section model of the crosssection suggested for the Sulafjorden crossing. The section is unsymmetrical and has curved undersides. The shape of the section model is milled from Divinycell H60 material, while two aluminium pipes with diameters of 40 mm and wall thicknesses of 1 mm provide the stiffness. The surface is foliated for a smooth finish. All tests are performed with the smallest deck upstream, without details such as handrails and windscreens. A total of four types of wind tunnel tests were carried out, as follows: (i) a slowly varying pitching motion is used to obtain the static force coefficients at two mean wind velocities; (ii) single harmonic motions tests are performed at two mean wind velocities, and several frequencies are considered; (iii) stochastic broad-banded motions are used to test the performances of the models for a more complicated motion pattern, and these tests are also performed at two mean wind velocities; and (iv) finally, the performance of the section concerning   vortex shedding is tested by keeping the section still while slowly increasing the mean wind velocity. Table 1 provides a summary of the tests.

Wind tunnel results
This chapter presents the experimental results from the wind tunnel tests. First, the static coefficients are presented together with a detailed discussion of their nonlinear characteristics. Further, tests at varying mean wind velocities are presented to investigate if the aerodynamic characteristics depends on the Reynolds number. Then, the experimental results considering single harmonic motions are given to obtain aerodynamic hysteresis and aerodynamic derivatives. Finally, the stochastic motion is considered to gain further insight into the observed nonlinearities and to study the importance of cross-terms of the input motions in the modelling.

Static force coefficients
Fig. 5 shows the static force coefficients obtained for the mean wind velocities of 6 and 10 m/s. The results for the two velocities are slightly different, which indicates that the forces are sensitive to the Reynolds number. The drag coefficient is, as expected, nonlinear since either pitching up or down gives higher drag and should result in a quadratic nonlinearity. The plot also shows that a higher-order nonlinearity is present since there is a bump in the curve at approximately 2 , which cannot be modelled by a quadratic polynomial. It is also interesting to observe that the bump at approximately 2 is present for both wind speeds. Still, the nonlinear characteristics are a bit different since the curves have different trends. The results for the lift coefficient are more consistent in the sense that the curves have an almost linear trend and that they are similar for both the velocities. The results do, however, show some discrepancies in the range of þ -2 , where the curves are slightly different and also show some indications of nonlinear behaviour. The results for the pitching moment coefficient are more interesting since apparent differences are present for the two velocities and nonlinearities are present. It is interesting to observe that the aerodynamic characteristics change significantly at þ1.5 at 10 m/s and that this change is not visible at 6 m/s. It is a well-known fact that the flow structure around a curved body is highly dependent on the Reynolds number. This dependency is problematic in wind tunnel testing since it is impossible to achieve the same Reynolds number in the model-and full-scale, which introduces uncertainties in the modelling of the wind loads (Kwok et al., 2012). reported that the aerodynamic properties of the cross-section of the Stone Cutters Bridge are sensitive to changes in the Reynolds number. This problem was avoided by performing all comparable wind tunnel tests at one representative wind velocity to ensure that the results corresponded to each other. Fig. 6 shows static coefficients as a function of wind speed. The overall shapes of the curves are similar to other curved cross-sections in the sense that the curves are relatively constant in the beginning and drop as the Reynolds number increases. Parts of the curve is not shown since significant vibrations of the section model due to vortex shedding were observed in the areas with slightly darker backgrounds. The first, second and third natural frequencies of the section model are 10.5 Hz, 42 Hz and 73 Hz, respectively. A lowpass filter with a 5 Hz cut off frequency is applied for the static test results to remove high-frequency noise.

Single harmonic motion
Measuring the aerodynamic forces when the section model is in single harmonic motion is standard procedure when investigating the aerodynamic properties of bridge decks. The experimental setup measures forces on each end of the section model with load cells. The measurements consist of self-excited forces, inertia forces due to the mass of the section model and also inertia and small damping forces due to the surrounding air. The contribution from the inertia forces thus needs to be removed to obtain the self-exited forces. Removing the inertia forces is a challenging task. The method suggested by (Han et al., 2014) where one perform the same test in-wind and in still-air and assume that the self-excited forces can be obtained as the difference of the measured time series have been used in this paper. This approach is an approximation since the still-air tests are not done in a vacuum and will thus be affected by the surrounding air. This procedure is, however, the most robust way of removing the inertia forces and the inaccuracies introduced are not significant. Fig. 7 shows the force coefficients as functions of the dynamic angle of attack defined in Eq. (5), where B * is chosen as B ¼ 0.74 m. The results from low-frequency motion should collapse onto the static coefficient, as stated in the full quasi-steady theory. Furthermore, linear aerodynamic models such as, for instance, aerodynamic derivatives give an inclined elliptical-shaped hysteresis. A hysteresis that deviates from this indicates nonlinear behaviour. The centre of the hysteresis does not correspond to the static coefficients for some of the cases. This can be caused by a quadratic (even order) nonlinearity and is discussed later in the paper. Studying the aerodynamic hysteresis plotted against displacement or rotation is a good method for analyzing the energy transfer in a system since the area inside the hysteresis is proportional to the energy generated or absorbed in each cycle. The direction of the cycle indicates if energy is absorbed (clockwise) or dissipated (counterclockwise). This is generally not the case for a hysteresis plotted against the dynamic angle of attack since the dynamic angle of attack includes the pitching motion, the derivative of the pitching motion and possibly the vertical velocity of the section. This will change the overall shape of the hysteresis. This is illustrated later in the paper.
Moment: Fig. 7f) shows that the moment coefficient obtained from the vertical motion tests at 6 m/s follows a slightly nonlinear static coefficient for the 0.25 Hz motion, while it is changes to a more elliptical shape as the frequency increases. This shows the interchange between the quasi-steady (nonlinear memoryless) and unsteady (linear with memory) regions, which are well documented in the literature, and is also one of the basic premises of the band superposition and hybrid models. It seems that the pitching moment at 6 m/s is well behaved for all the tests considered. The same conclusion for the pitching motion at 6 m/ s can also be drawn but is not shown here.
On the other hand, the data at 10 m/s in Fig. 7d) and e) show explicit nonlinear behaviour for both the low-and high-frequency motion. The shapes of the hysteresis are considerably different for the vertical motion and the pitching motion. A similar observation is made in (G. Diana, Rocchi, et al., 2013). The experimental results for this particular section illustrate that it is not possible to model the forces as a function of the Drag: Fig. 7a shows that the drag coefficient is nonlinear and that the drag hysteresis twists, which are also observed in the literature (G. Diana et al., 2010;Zhou et al., 2018). It is also observed that the hysteresis curve does not evolve around the static coefficients and deviates quite a bit from them. This is discussed further later in the paper, but it is important to note that the self-excited drag forces are, in general, small compared to the mean drag force and the inertia forces that need to be removed. This makes the experimental results for the self-excited drag forces slightly more uncertain than the others.
Lift: Fig. 7c shows the lift coefficient obtained at the mean wind velocity of 6 m/s. The hysteresis rotates and becomes larger when the frequency is increasing, which is similar to the curve for the pitching moment. The static lift coefficient is nearly linear, and the low-frequency motion is also nearly linear, indicating that the static coefficient is a good indicator of whether nonlinear characteristics can be expected. Fig. 8 shows the drag, lift and pitching moment coefficient hysteresis loops for the vertical motion at 6 and 10 m/s. The frequency of the single harmonic motion at 6 m/s is 1.1 Hz, while the frequency at 10 m/s is 1.7 Hz, making the reduced velocity for both cases close to 5.5. It is a general hypothesis that the aerodynamic forces can be modelled using reduced velocity and, thus, reduced frequency. The results presented illustrate that this is not the case for the particular section considered. The results in Fig. 6 shows that the static force coefficients are sensitive to changes in the Reynolds number. This dependency is probably the main reason that the two hysteresis phenomena are so different. In conclusion, the Reynolds-dependency also seems to affect the dynamic flow around the cross-section and the associated forces.
The aerodynamic derivatives are also functions of the reduced velocity. It is, thus, interesting to study the aerodynamic derivatives of the section and see if any differences can be observed. Fig. 9 presents the aerodynamic derivatives from the standard forced vibration tests extracted using the Han method (Han et al., 2014). The results are shown in the Zasso convention (Zasso, 1996). The blue circles are 10 m/s tests, while the red circles are the 6 m/s tests. It is clear that for some of the aerodynamic derivatives, the results cannot be characterized by the reduced velocity alone since different trends for the tests at 6 and 10 m/s are observed. It should be noted that modelling self-excited forces using aerodynamic derivatives yields inaccurate results for the sections that  exhibit nonlinear aerodynamic behaviour. The observations made are, nevertheless, a clear indication of that Reynolds number dependency is an issue in this setting also. Fig. 8 also shows that the mean forces for a dynamic motion are different from the static values at a zero angle of attack. Fig. 10 shows a time series of the measured motion-induced pitching moment at 10 m/s. The figure shows the results for four single harmonic motions with different frequencies. The force measured between the harmonic motions is the static force, and it is observable that the mean value of the motioninduces forces does not coincide with the static forces. This difference is a nonlinear effect that is not taken into account in a linear framework such as, for instance, modelling motion-induced forces using aerodynamic derivatives. Fig. 10 also shows the same force time series for drag and lift. It is observed that the effects are considerably less for the more linear lift force but are very prominent for the drag force. If quadratic nonlinearities are present, a mean drift is expected. Fig. 11 shows the change in the mean value for all the forces for all the harmonic tests performed. Two mean wind velocities, i.e., 10 m/s, and 6 m/s, and a frequency range between 0.25 Hz and 2.5 Hz are given. The figure unveils a more complicated situation than in Fig. 10, where it  velocities. The blue lines shows the results at V red ¼ V/fB ¼ 7.4 (6 m/s, vertical motion, 1.1 Hz) and V red ¼ 7.9(10 m/s, vertical motion, 1.7 Hz) but with different Reynolds numbers, 2.1e4(6 m/s) and 3.5e4(10 m/s). Mean hyst denotes the mean over the entire hysteresis, and Static indicates the static value at an angle of attack equal 0 . (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) seems that the mean drift is constant. Some general comments can be made. The mean drift for the 6 m/s tests tends to trend towards the 10 m/ s static coefficient, and vice versa. The effect of a rotation of 2 seems to be more frequency-independent than vertical motion, and this could be due to an effective angle of attack dependence. The drift is, in general, larger for higher wind velocities. Note that the static coefficient is not the same for a different test at the same wind-speed because the static coefficients are found from the start of that actual test to minimize the uncertainties in the experiment. The differences in the moment and lift are small but are more significant for drag. This could be due to the small angular offset of the mean angle of attack in the incoming flow and/or the section model alignment, slightly different wind velocity or other sources. Nevertheless, the errors are rather small and do not alter the conclusion regarding the mean drift. Fig. 12 shows an in-depth view of the force coefficients measured at 10 m/s for a single harmonic pitching motion at 2 Hz, with an amplitude of 2 . Fig. 12a) shows the time series of all 16 cycles in light grey and the mean of the cycles in solid red, and it is observed that the drag coefficient  shows less repeatability compared to the lift and moment coefficients. The hysteresis shows that there are higher-order components and, thus, significant nonlinearities in all tree load components since the hysteresis do not have elliptical shapes. The dynamic coefficients are plotted as functions of the dynamic angle of attack. Fig. 12d) shows the same experimental results, but here the dynamic coefficient is plotted against the pitching motion. Comparing Fig. 12a and d, it is observed that even though it is the same experiment that is shown, the two hysteresis have different shapes. This is obviously because one set of hysteresis is plotted as function of the angle of attack while the other is plotted as a function of the pitching motion. To be able to evaluate the energy transfer in the system, it is most convenient to plot the hysteresis as a function of the motion and not the angle of attack since the area of the hysteresis is then equal to the work done. Energy dissipation is seen for drag and moment (counterclockwise direction of the plot, not shown), and an interplay between generation and dissipation of energy is seen for the lift force.
The Fourier amplitudes in Fig. 12c) also show the higher-order harmonics clearly. The higher-order harmonics unveil the integer multiples of the fundamental motion frequencies. Fig. 12b) shows time-domain measurements of the motion-induced forces, and the deviation from the sinusoidal shape of the motion-induced forces indicates nonlinearities. In conclusion, all three load components show clear signs of nonlinearities.

Stochastic motion
Wind-induced vibrations range from near harmonic motion due to ViV lock-in to stochastic motion due to buffeting response in highly turbulent wind. Traditional wind tunnel testing is based on single harmonic motions. In the linear range, superposition of single harmonic motions is satisfactory for considering all types of motion, but in the nonlinear range, the superposition principle does not hold. Therefore, it is relevant to test the cross-section in stochastic motion rather than only doing simple single harmonic tests.
The Cholesky decomposition approach (Shinozuka, 1971) is used to simulate vertical and pitching motion, with the assumption that the spectral densities are constant between 0 and 2.5/3.5 Hz and zero elsewhere for both the vertical and pitching motion and that the vertical and pitching motion are uncorrelated.
Identifying nonlinearities in a stochastic time-series can be challenging, but numerous approaches and tools exist. One of them is to examine the Fourier spectrum of the input and the output and search for harmonic distortion. Fig. 13 shows the Fourier amplitudes of the input and outputs. The testing rig is precise enough that the cut off frequency of 3.5 Hz is strictly held. All energy over this limit comes from harmonic distortion. It is, however, difficult to say much about possible harmonic distortion in the frequency range of the stochastic broad-banded input motion. Higher-order spectra can be used for this purpose (Hosseini et al., 2011); the higher-order spectrum is the Fourier transforms of the higher-order cumulants (Kim and Powers, 1979;Udea and Dowell, 1984). The higher-order spectrum can detect higher-order nonlinearities through the inclusion of phase information (Nikias and Raghuveer, 1987), the bispectrum can detect quadratic nonlinearities, the trispectrum can detect cubic nonlinearities, and so on. We chose to use the Fig. 11. Mean value of the static part and the harmonic part of the consecutive harmonic motion series. Vert ¼ vertical motion, Pit ¼ pitching motion. The X-axis names indicate the amplitude, varamp denotes amplitudes equal {95,80,50,35,28,23,20,16 mm} for the 0.5-2.5 Hz harmonics.
cross-bicoherence between the input and output to detect nonlinearities (Hagihira et al., 2001;Hayashi et al., 2007)  transform of the kth window of signal b and a, respectively. * denotes a complex conjugate, and |x| denotes the absolute value. The crossbicoherence is always between 0 and 1. Fig. 14 and Fig. 15 show the cross-bicoherence spectrum between the motion and the forces at 6 and 10 m/s, respectively. The bicoherence is set to zero if the bispectrum is less than 0.1% of its maximum value to make the figures easier to interpret. Table 2 gives the maximum value and the mean value of the bicoherence. The quadratic distortion of the pitching moment at 10 m/s is considerable. The drag force has a considerable amount of energy that comes from the quadratic coupling. The lift force has less quadratic distortion for the considered motions, which indicates a more linear behaviour. The spectral density in Fig. 13 confirms this observation since there is considerably less frequency content above 2.5 Hz for the drag force compared to the others. Some articles discuss whether it is necessary to include cross-terms, such as r z r θ and r z r 2 θ , or if it is sufficient to only consider direct terms, such as r 2 z and r 3 z , in the modelling (Wu et al., 2013;Wu and Kareem, 2011). It is generally necessary to include cross-terms for a nonlinear system since the principle of superposition is not valid. The validity of the assumption of superposition can be investigated experimentally. If we, for instance, assume that we can model the lift force due to vertical and pitching motion separately, the following expression applies.
where F z ðtÞ is the lift force and f ð⋯Þ symbolizes a nonlinear function with or without memory. The stochastic motion is used to test whether Eq.
(11) is a fair approximation. A total of three stochastic motions involving both the vertical and pitching motion are first generated. Each realization is tested in the wind tunnel three times. Only the vertical motion is applied in the first test, only the pitching motion is applied in the second  test, while the third test is the sum of the vertical and pitching motion from the first and second test. Fig. 16 shows a 20 s window of one of the realizations tested at 10 m/s in the wind tunnel. The top subplot shows the motions while the other subfigures show the drag, lift, and pitching moment coefficients. The blue lines are the measured self-excited forces for the combined motion F z ðtÞ ¼ f zθ ðr z ;r θ Þ, while the red lines are the selfexcited forces obtained by adding the contribution from only the vertical motion and only the pitching motion F z ðtÞ ¼ f z ðr z Þ þ f θ ðr θ Þ. It is observed that the approximation suggested by Eq. (11) is almost perfect for lift, is acceptable for the pitching moment, while it is a bit poorer for the drag coefficient. Note that the drag forces have a much lower signal to noise ratio than the other force components since the self-excited drag is small. Table 3 shows the normalized mean square error for all the combinations tested. The normalized mean square error is 1 if the time series are equal, while it is À ∞if the time series do not match. It should be noted, as indicated in the table, the same motion has been used for the tests at 6 m/ s, and 10 m/s. In other words, the changes in the performance when considering the same stochastic motion are caused by different wind velocities. The bicoherence, static coefficients and hysteresis loops presented earlier show that the pitching moment exhibits more nonlinear behaviour at 10 m/s compared to at 6 m/s. The results also show that the lift force exhibit a fairly linear behaviour at both wind velocities and that the drag force measured at 6 m/s seems more nonlinear than the one obtained at 10 m/s. Table 3 supports these observations and gives a clear indication that the deviation from the superposed results is a distinct feature for the cases where nonlinear behaviours were observed earlier. If the forces generated by vertical and pitching motion cannot be superposed, it means that the nonlinear model needs to include the products of the input motions, such as _ r z r 2 θ . Fig. 16 shows a plot of the force coefficients obtained for the combined vertical and pitching motions and the superposed results from vertical and pitching movements tested separately. The figure also shows the individual motions. Considering the pitching moment coefficient, we see that the deviation is most significant at the peaks of the time series. The results also show that the error is at its largest when both motions are large.
This observation is a clear sign of the interaction between the two motions and can only be modelled by the cross-terms in the model. For drag, the peaks of the combined motion are almost always higher than from the superposed results. This finding makes sense due to the quadratic nature of the drag force. The superposition principle is clearly a better approximation for the lift force; however, deviations in the peaks are visible and are most significant when both the vertical and pitching motions are large at the same time.
It is also interesting to note that the NMSE results from Table 3 are similar for the different motions tested at 6 m/s, which indicates that the duration of the stochastic motion is sufficient and that the results are consistent.

Computational results
Predicting self-excited forces using several of the models presented is considered to be out of the scope of this paper. A rather simple responsesurface model is instead proposed and used to model the self-excited forces. The simplicity of the model makes it a good candidate for evaluating the necessary characteristics of nonlinear load models because the model can easily be expanded, and the contribution of the different components can be deducted.

Response surface
The response surface proposed is a polynomial combination of variables describing the motion-induced forces in laminar flow, as follows: Note that in this description, all the displacement, velocity and acceleration terms are included, and no terms that might be of less importance are disregarded. Since the velocity and acceleration depend on the motion history, the memory effect is partly taken into account. It is, however, clear that this approach is not capable of modelling fluid memory as accurate as by using convolution integrals.
The self-excited force F i ðxÞ expressed in terms of a response surface reads as follows: x 4 þ ::: β 11 x 2 1 þ β 12 x 1 x 2 þ β 13 x 1 x 3 þ ::: β 111 x 3 1 þ β 112 x 1 x 1 x 2 þ β 113 x 1 x 1 x 3 þ ::: :::þ (13) Fig. 15. Cross-Bicoherence between the motion and force, 10 m/s tests. The motion region is between 0 and 3.5 Hz. where β n are coefficients that can be obtained by fitting the model to the experimental data. For response surface models with low orders, a physical interpretation of the parameters is possible. Still, it may become difficult to consider the higher-order terms as physical quantities since several combinations of parameters may provide an almost equally good fit to the experimental data. Since the model is a part of a class of response-surface models, standardized methods and software for identification of the model coefficients are readily available. The model is also linear in its coefficients, making it possible to use linear regression. The authors used the MATLAB function fitlm to obtain the results presented in this paper. Both pitching and vertical motions are used as inputs in Eq (13). Simpler models can be developed for galloping, vortex-induced vibration and torsional flutter since only vertical or torsional motions are relevant. Simplifications can also be introduced when considering limit cycle oscillations where the vertical and torsional motions are related (Guang zhong Gao and Zhu, 2017). presents a model for transverse instabilities of slender rectangular prisms where such assumptions have been successfully introduced.
Due to the strong Reynolds number dependency seen in Fig. 6, the coefficient of the response surface model will be calibrated on single wind speeds, and therefore in effect become Reynolds number dependent. At a given Reynolds number the coefficients do not change, and the model can be used with any arbitrary motions as long as the motion is within the range of the training data.

Force predictions
The response surface presented above can be challenging to use since there are many unknown coefficients when applying a high-order model. It is also an open question if it is necessary to include all the terms in the modelling. It is recommended that one starts with a low order model and adds terms gradually. It is also vital that the performance of the model is investigated considering a different dataset than the one used to obtain   Fig. 17 shows the normalized mean square errors (NMSE) of different response-surface models when considering the stochastic motions discussed earlier in the paper. The results for the pure vertical (V) and pitching (P) motions, as well as the combined vertical and pitching (PV) motions are presented for the mean wind velocity values of 6 and 10 m/s. The acceleration terms in Eq. (12) have been omitted for the drag and lift forces, while they are included when predicting the pitching moment.
The results for the lift coefficient correspond very well with the validation data, and the higher-order models perform only marginally better.
The bicoherence, hysteresis loops and static coefficients presented earlier clearly show that the self-exited drag is nonlinear. Fig. 17 also support this observation since the second-order model fits significantly better than the linear model for all cases except for the drag force caused by the vertical motion at 10 m/s. It seems that the full second-order model, including cross-terms, is the best model when balancing low order and high performance. For the combined pitching and vertical motions, the same conclusion can be drawn since the full second-order model performs the best.
The results for the pitching moment coefficients in Fig. 17 are presented both with and without the acceleration terms. The models with acceleration include the x 5 and x 6 -terms in Eq. (12), while these, as pointed out earlier, have been disregarder for all the other models presented since they had a very minor influence on the results for the drag and lift. The terms are, however, of crucial importance for the modelling of the pitching moment, and it is likely that this is because the memory effect is more important for the pitching moment and that the acceleration terms try to compensate for this. The coefficient b 5 corresponds to an added mass moment of inertia and it is, therefore, interesting to study its magnitude. If we consider the section as two ideal flat plates next to each other with no aerodynamic interference, it is possible to obtain an added mass moment of inertia based on the theoretical values for the ideal flat plate. This is clearly an oversimplification, but the added mass corresponding to the b 5 coefficient is approximately 10 times larger, which indicates that the obtained coefficient is slightly too high to have a physical meaning. This observation may indicate that the model is trying to compensate for effects that are not included in the model, for instance, proper fluid memory. It is also interesting to note that the theoretical and calculated added mass for the lift were approximately equal for the first order model when the acceleration terms are included in the modelling. If the importance is that the acceleration terms where due to higher-order nonlinearity in the system, the fit would be better for higher-order models, but this is not observed in the results. Nevertheless, it is concluded that the acceleration terms are essential for the prediction of the pitching moment coefficient. In the following, only the models Fig. 17. NMSE values considering various response surface models. Lin ¼ linear direct terms. Lin þ diag ¼ linear model plus second-order direct terms. Lin þ cross ¼ linear model plus products of inputs. Full 2nd, 3rd and 5th represents full second third and fifth order response surfaces. T ¼ pitching motion, V ¼ vertical motion, PV ¼ vertical and pitching motion, WA ¼ with acceleration, and WOA ¼ without acceleration.
including the acceleration terms are considered. The results for the pitching moment coefficient at 6 m/s indicate linear behaviour because no higher-order model performs significantly better than the linear one. The earlier presented data also support this observation. For the 10 m/s tests, the full second-order model seems to be the model with the lowest order and good performance for both the vertical and pitching motions. It seems that the pitching moment due to the pitching motion is more nonlinear then the pitching moment due to vertical motion. It is interesting to observe that a model with full 2nd-order fits better than 1storder with second order direct or first-order cross-terms for the test at 10 m/s with the combined pitching and vertical motions (PV). Fig. 18 shows a comparison of the measured pitching moment coefficient and predictions using different response surfaces. The linear model without acceleration has less amplitude than the measured one, and it lags a bit behind. The lag is hardly visible in the plot but can be important for motion-induced damping of the final structure. The secondorder model with acceleration fits best but still struggles to capture the peaks in the experimental data.
Nonlinear models are, in general, more difficult to identify than linear models since the principle of superposition is not valid, which in turn, makes the type of motion used for training more important. It is, therefore, interesting to study how the model performs for types of motion other than those used in the identification process. Fig. 19 shows the hysteresis loop predicted using the full second and third-order response surfaces presented above. The single harmonic motion is a pitching motion at 1.1 Hz with 2 amplitude, and the mean wind velocity is 6 m/s. The models predict the measured self-excited forces with decent accuracy, but some discrepancies are clearly present. The results are best for the lift force and pitching moment, which is probably explained by their close to linear behaviour. The prediction for the drag coefficient is slightly inaccurate, and it is, therefore, interesting to observe if better accuracy can be obtained if a stochastic motion that is closer to a single harmonic motion is applied. A motion corresponding to a single degree of freedom system with a damping ratio of 1.66 (overdamped) driven by white noise is therefore applied to identify an alternative response surface to investigate if the results are sensitive to the motion applied. The predictions obtained using the resulting response surface are also shown in Fig. 19. The results show that the performances are very similar for the lift and the moment coefficients and that the response surface obtained when using a narrow-banded motion actually gives results that are closer to the predictions of the other response surfaces than the measured data. A significant difference is observed in the prediction of the drag coefficient, but the results did not improve significantly.

Concluding remarks
A series of experimental tests have been conducted to gain further insight into the nonlinear characteristics of a double-deck section with curved undersides. A simple response surface model has been proposed and validated against experimental data. The following observations are made.
Cross-sections having aerodynamic characteristics that are sensitive to changes in the Reynolds number can exhibit nonlinear aerodynamic behaviour for some wind speeds and exhibit linear behaviour for others. A significant change in the nonlinear characteristics as a function of the mean wind speed is difficult to include in any model based on scaled model tests and pose a significant challenge for the research community. One of the main challenges in nonlinear modelling is to select a nonlinear function that is capable of modelling the nonlinear characteristics and at the same time ensuring that a robust identification of the model coefficients can be carried out. It is, therefore, important to know which effects must be included in the modelling and which can be disregarded. It is generally necessary to include cross-terms, such as _ r z r 2 θ , for a nonlinear system. Designated wind tunnel tests were carried out to investigate the importance of the cross-terms. It can be concluded that the cross-terms are essential when the force components exhibit significant nonlinear behaviour. This implies that it is, in general, not possible to model the force component separately by applying the superposition principle. This makes both the wind tunnel testing and model identification more complicated. The simple response surface model proposed is capable of partly modelling a memory effect when including displacement, velocity and acceleration terms as inputs. This memory effect was essential to include in the modelling of the pitching moment for the particular section tested. The Volterra series is the natural extension of the response surface model applied in this paper since many of the same terms can be used as inputs. It is then sufficient to use only displacements as inputs since the convolutions model the memory effect. Multiple input single output second-order Volterra models seem to be attractive alternatives since good performance was achieved using a quadratic response surface, including cross-terms, in this paper. The response surface model shows promise in the modelling of weakly nonlinear self-excited forces for low amplitude stochastic random Fig. 18. Timeseries plot of the measured and predicted pitching moments for a 10 m/s, pitching stochastic motion with a frequency content between 0 and 3.5 Hz. The 1st-order model without acceleration lags 0.12 s behind. The plot is part of a 600 s test. WA ¼ with acceleration and WOA ¼ without acceleration. motions. The model is validated using time histories of measured forces for motons that were not used to obtain the model coefficients. Random motions corresponding to rectangular auto-and cross-spectral densities and as well as small-amplitude single harmonic motion data was applied. The performance of the model concerning aeroelastic phenomenon resulting in limit cycle oscillations remains an open question and should be studied in detail in future work. Since it might be possible to disregard some of the terms in the model, it is also recommended that the model is carefully reviewed before it is used to model nonlinear aeroelastic pehnomena that results in limit cycle oscilations.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Fig. 19. Measured vs predicted dynamic coefficients from single harmonic motion. 6 m/s, 1.1 Hz single harmonic pitching motion. Predicted forces are from the response surface models calibrated on 6 m/s pitching stochastic and narrow motion and data. The hysteresis loops are plotted against the dynamic angle of attack in accordance with Fig. 7. WA ¼ with acceleration.