Abstract
Dynamic instability of beams in complex structures caused by unsteady wind load has occurred more frequently. However, studies on the parametric resonance of beams are generally limited to harmonic loads, while arbitrary dynamic load is rarely involved. The critical frequency equation for simply supported Euler beams with uniform section under arbitrary axial dynamic forces is firstly derived in this paper based on the Mathieu-Hill equation. Dynamic instability regions with high precision are then calculated by a presented eigenvalue method. Further, the dynamically unstable state of beams under the wind force with any mean or fluctuating component is determined by load normalization, and the wind-induced parametric resonant response is computed by the Runge-Kutta approach. Finally, a measured wind load time-history is input into the dynamic system to indicate that the proposed methods are effective. This study presents a new method to determine the wind-induced dynamic stability of Euler beams. The beam would become dynamically unstable provided that the parametric point, denoting the relation between load properties and structural frequency, is located in the instability region, no matter whether the wind load component is large or not.
1. Introduction
The simply supported Euler beam, an important research object on the dynamic instability of elastic systems, is a common form of members in modern structures [1, 2]. Modern building structures are being constructed with greater height or span, resulting in smaller natural frequencies of beams in them, so parametric resonance of the beams may occur under unsteady strong wind forces. For example, many cases of roof collapse in windstorms caused by instability of significant beams have been reported [3, 4].
The dynamic instability of Euler beams has attracted wide attention of scholars in the early stage. Beliaev [5] analyzed the dynamic responses of simply supported Euler beams under axial load with the form of Main frequencies of parametric resonance were calculated, and the Mathieu-Hill equation for structural dynamic instability was derived. As a complement to Beliaev’s analysis, Krylov and Bogoliubov [6, 7] studied the influence of arbitrary boundary conditions on dynamic instability of Euler beams by the Galerkin method. To successfully solve the Mathieu-Hill equation for the problem of dynamic instability of Euler beams, Shtokalo [8] established the stability of linear differential equation with variable coefficients. McLachlan [9] discussed the theory and application of Mathieu equation. Based on these theories, Bolotin [10] built the critical frequency equation for structural dynamic instability and gave the approximate formula of the first three dynamic instability regions. With the development of finite element (FE) theories, Briseghella et al. [11] computed dynamic instability regions of Euler beams by FE method. He also compared the FE results with those obtained from the analytical method. Yang and Fu [12] discussed the dynamic instability of Euler beams with layered composite materials. Sochacki [13] concerned on simply supported Euler beams with additional elastic elements. Yan et al. [14] analyzed the effect of cracks on dynamic instability of Euler beams. However, in these previous studies for the engineering applications of the Mathieu-Hill equation, axial loads are all assumed as harmonic loads, whereas the wind load is arbitrary and cannot be simplified as a single harmonic wave in studying wind-induced parametric resonance. Therefore, it is necessary to establish an analytical method for computing the dynamic instability of Euler beams under arbitrary dynamic load.
In this research, the critical frequency equation of simply supported Euler beams with a uniform section under arbitrary dynamic load is firstly derived based on the Mathieu-Hill equation. An eigenvalue algorithm is then developed to compute dynamic instability regions with high precision. Further, the approaches for judging dynamically unstable state under any unsteady wind force and calculating wind-induced resonant responses are presented. Finally, a measured wind load time-history is adopted to describe the application of the proposed methods. At the same time, a discussion on the characteristics of wind-induced dynamic instability of Euler beams is also given.
2. Mathieu-Hill Equation for Arbitrary Dynamic Axial Load
The studied simply supported Euler beam with a uniform section is shown in Figure 1, whose length is , Young’s module is , and the moment of inertia of the section is . Under the action of axial unsteady wind force , the deflection of the beam can be denoted as [15] where is the time-varied deflection along the beam, is the th modal shape of the beam, and is the th modal coordinate.
The arbitrary force with the duration of can be transformed into Fourier series by even prolongation [16]: where and are the mean and fluctuating component of , respectively. and are the amplitude and angular frequency of the th harmonic wave, respectively, and .
If is measured by wind force time-history composed of discrete data from wind tunnel lab or in field, the maximum is .
Under the wind force, the Mathieu-Hill equation, corresponding to the th mode of the undamped Euler beam, can be written as [10] where is the th natural frequency of the beam loaded by and is the excitation parameter of the th harmonic wave for the th normal mode.
and can be gained by where is the th natural frequency of the unloaded beam and represents Euler’s critical buckling axial load of the beam.
The Mathieu-Hill equations corresponding to normal modes of different orders have the same form as long as and are calculated from the same order of mode, so the distributions of dynamic instability regions in different modes are identical. Consequently, the subscript in (6)-(7) can be omitted, and all calculations later will use the data corresponding to . The Mathieu-Hill equation is rewritten as
As is well known, the dynamic instability regions are surrounded by two solutions of the Mathieu-Hill equation with the same period, while the dynamic stability regions are encompassed by two solutions with different periods [10]. Hence, provided the condition that there exist periodic solutions of the period or () in (8) is found, the dynamic instability regions can be determined.
Periodic solutions of can be written in the form of Fourier series: where and are constant coefficients.
Substituting (9) into (8) and making the coefficients of terms of or be zero generate homogeneous system of linear equations with respect to and :
If and only if the coefficient determinant is equal to zero, (10) has nonzero solutions. Therefore, the condition that (8) has periodic solutions of is that the coefficient determinant of (10) is equal to zero. Considering both cases of and , the critical frequency equation corresponding to periodic solutions is obtained:
The critical frequency equation contains the relationship between the frequency spectrum characteristics (amplitude and frequency characteristics) of wind load and structural natural frequency. By solving the equation, the boundary of dynamic instability region , expressed by the ratio of critical load frequency and structural frequency, would be obtained. is called the frequency ratio, and is called the critical frequency ratio.
Similarly, in order to gain the periodic solution condition of the period , the periodic solution of period is written in the following form:
By substituting (12) into (8) and making the coefficients of terms of or be zero, the critical frequency equation corresponding to the periodic solutions of period is acquired:
By solving (11) with order determinant, the first odd dynamic instability regions are obtained. Besides, by solving (13) with order determinant and (14) with order determinant, the first even dynamic instability regions are achieved. Thus, the first 2 dynamic instability regions can be determined. Meanwhile, it is known from the above derivation that 2 should not exceed the maximum harmonic number .
3. Eigenvalue Algorithm for Computing Instability Boundary
Based on the critical frequency equation, the method of multiple scales or successive approximation is commonly used to determine the boundaries of dynamic instability regions [17–20]. However, the computation complexity or large time consuming of these methods generally increases rapidly with solution precision requirement. To ensure computation efficiency, less order of the determinant is usually selected, which results in less accurate instability boundaries. Here, a fast algorithm based on the eigenvalue problem is attempted to be established for acquiring instability boundaries with high precision.
Take the critical frequency equation (11); for example, matrices and can be defined as follows:
Thus, (11) can be expressed as
It is obvious that (16) converts the problem of solving dynamic instability regions to the eigenvalue problem [21, 22]. Since the eigenvalue problem is easy to be solved by means of a computer and the order of determinant has few effects on calculated amount, much more order of the determinant can be selected to meet high precision requirement. For matrices and of order, there are 2 eigenvalues , which generate boundary points of dynamic instability regions. By altering the value of , the spatial instability region in 2-dimensional coordinate system () can be constructed [23]. For the beam with known natural frequency, any unsteady wind force determines a point in the spatial coordinate system. If the point is located in an instability region, the beam would become dynamically unstable under such wind force and dynamic instability wouldnot occur otherwise.
Analogously, the instability boundaries according to the periodic solution of can be determined by solving the critical frequency equations (13) and (14) through the eigenvalue algorithm. For (13), the matrices and are The matrices and in (14) can be expressed as follows:
4. Calculation of Wind-Induced Dynamic Instability Region and Resonant Response
The value of is usually fairly large to meet precision requirement, so it is difficult to graphically display the 2-dimensional spatial instability region and prejudge the dynamic instability state. Here, a proportional load strategy is adopted to augment the mean or fluctuating component of wind force and set up the two- or three-dimensional dynamic instability regions [24, 25].
The wind force as expressed in (3), called the initial load here, is firstly normalized as follows: where is the maximum harmonic amplitude of and is the normalized wind force, called the basic load here.
According to (20) and (21), it is clear that the mean value and maximum harmonic amplitude of the basic load are 0 and 1, respectively. Each harmonic amplitude of the basic load is times less than that of the initial load .
Let denote the harmonic amplitude of the basic load, and then the basic load can be written as
Thus a new wind force can be generated by defining two parameters and to linearly increase the mean and fluctuating components of wind force based on the basic load: where is the newly generated wind force; and are the mean component and the maximum harmonic amplitude of , respectively; and are the ratios of and to Euler’s critical load of the beam. The reason why Euler’s critical load is adopted to denote the mean and fluctuating components is that the structural natural frequency and excitation parameter are both related to Euler’s critical load as shown in (7).
For the new load, its mean component of is and its fluctuating component or harmonic amplitude is times that of the basic load, while the harmonic frequency composition is identical to the initial load or the basic load.
is further simplified as follows: where denotes the harmonic amplitude of .
Meanwhile, the natural frequency and excitation parameter corresponding to can be computed by
The expression of excitation parameter can be further simplified by defining where represents the ratio of excitation parameter and harmonic amplitude of the basic load.
Thus, the excitation parameter of generated wind force can be calculated by harmonic amplitude of the basic load:
Moreover, the Mathieu-Hill equation consistent with is
That is
Therefore, various can be obtained only by adjusting or , and then the two-dimensional parametric plane or three-dimensional parametric space for dynamic instability regions can be graphically established according to the theory in Section 3 of this paper. For different wind forces, various distributions of instability regions will be drawn. Because the mean component and the peak harmonic amplitude are generally less than Euler’s critical load, there exist and .
When the parametric point is located in the instability region, the beam is deemed to be dynamically unstable. In order to verify such a conclusion, modal response of the beam can be calculated through the Mathieu-Hill equation (29). should keep increasing with time when the beam is unstable.
The Mathieu-Hill equation shown in (29) can be converted to the following form: Let Equation (29) is expressed in the vector type
For differential equations with the form of (32), the fourth-order Runge-Kutta method can be used to solve the unknowns [23, 26]. The response at the th iteration step can be obtained by the following iteration formula:
The time step length is set as s and initial disturbance of response is m in the latter computation.
5. Case Study
In order to describe the application of the above approaches, a real wind force time-history is taken from a wind tunnel test. The rigid model tests of a cylindrical reticulated shell, with the span of 103 m, the longitudinal length of 140 m, and the height of 40 m, were carried in a TJ-2 wind tunnel of Tongji University, China, in order to obtain the time-histories of wind pressure on the shell (Figure 2). Detailed information on the test can be referred to Zhou et al. [27]. Here, wind pressure at one node of the shell is extracted and converted to full scale, and the area of 16 m2 is employed to construct wind force on the beam, as shown in Figure 3.
The wind load data number is , and the sampling frequency is Hz, so the load duration is s. The load is simulated by the periodic load with period according to (3)–(5). Comparison of and the measured data is shown in Figure 4(a). To appreciate the two sets of data more clearly, a comparison of shorter load duration is also given in Figure 4(b). It can be seen that the original data are well fitted by the Fourier series. Mean component of the wind force is N and the maximum harmonic amplitude is N by (4), (5), and (20). The relationship between harmonic amplitude and corresponding harmonic frequency is shown in Figure 5, which indicates that larger amplitudes correspond to smaller frequencies. Moreover, the natural frequency of the beam is defined as Hz.
| (a) Load duration = |
| (b) Load duration = /10 |
Normalizing the initial load through (20)–(22) produces the basic load , whose magnitude and the relationship between harmonic amplitude and frequency are shown in Figure 6. Comparing with the initial load (Figures 4 and 5) shows that mean component and the maximum harmonic amplitude of the basic load are 0 and 1 m, respectively. Both the fluctuating components and harmonic amplitudes are reduced to times, while the harmonic frequencies are invariant.
(a) Load time-history
(b) Relationship between the amplitude and frequency
New wind load can be generated through (23)-(24) based on the basic load, and the dynamic instability regions can be established by substituting excitation parameters of the new load into (15)–(19). Figure 7 gives the left and right boundaries of the first dynamic instability region when is equal to 10, 40, and 80, respectively. It is obvious that the boundary values vary with and larger results in close boundaries with higher precision. Since the results under and are close, is set for 40 in the following calculations. Meanwhile, for convenience, only the first 100 data of wind force are selected for later computations ().
Bolotin [10] has pointed out that, if diagonal elements of the critical frequency equation are ignored, the boundary of the th instability region under arbitrary dynamic load can be calculated by the following approximate formula:
Equation (34) indicates that the th instability region only depends on the th harmonic wave. Figure 8 compares the first instability regions gained by (34) and the eigenvalue method, respectively. It is seen that two methods produce different results, especially under larger excitation parameters. Comparing other regions brings similar conclusions. This is because the influence of the th harmonic wave on the width of the th instability region is the order of [10], while the value of under wind load could be fairly large. Therefore, the interaction between different harmonic waves cannot be ignored for wind-induced dynamic instability analysis.
In order to further verify the precision of the eigenvalue method, points A(0.65, 0.8) and B(1.31, 0.8) on the parametric plane of Figure 8 are selected. Point A is located in the unstable region of the approximate formula while in the stable region of the eigenvalue method. Point B is located in the stable region of the approximate formula while in the unstable region of the eigenvalue method. Modal responses corresponding to points A and point B are calculated through (33), as shown in Figure 9. The modal response of point A is confined in a scope of about 0.01 m, while that of point B increases rapidly to about 5.0 m in the same duration. Hence, the structural system corresponding to point A is dynamically stable, but that corresponding to point B is dynamically unstable, which is consistent with the conclusion of the eigenvalue algorithm. This indicates that the eigenvalue algorithm provides more accurate dynamic instability region for arbitrary dynamic load.
(a) Point A
(b) Point B
Distribution of the first 40 dynamic instability regions (the first 20 odd regions corresponding to period and the first 20 even regions corresponding to period ) is shown as the shaded area in Figure 10. When the parametric point corresponding to the specific structure and wind force falls in the shaded area, structural dynamic instability occurs. In addition, according to the instability regions, the scope of critical instability frequency can be obtained for certain denoting the magnitude of wind force and the range of causing dynamic instability can also be evaluated when the frequency ratio is fixed.
When is assumed, we have according to (26). By s, we get Hz and . Thus, the wind force and the beam correspond to a parametric point C(0.043, 0.5), which is located in the shadow area of the parametric plane in Figure 11. Figure 11 is the partial enlarged view of Figure 10 in the scope of between 0.02 and 0.065. The modal response of point C is also shown in Figure 12. It is shown that the deflection of the beam increases rapidly in a very short time. Therefore, obvious dynamic instability is caused by this wind force.
The above analysis indicates that the dynamic instability of simply supported Euler beam with uniform section under any wind load can be determined according to the values of , , , and . For any wind load, its basic load can be constructed and the distribution of dynamic instability regions similar to Figure 10 can be drawn to determine the state of dynamic stability.
In order to study the influence of mean or fluctuating component on the scope of critical load frequency, parametric planes (, ) and (, ) can be built. Figure 13 provides the first 40 instability regions on the plane (, ) when , 0.5 and 0.99. It is demonstrated that the increase of generally widens the scope of critical frequency ratio under the same . That is, when mean component of wind load is fixed, increase of fluctuating component may transform the beam from dynamically stable state to unstable state. Moreover, even though the mean component is zero, increasing fluctuating component could also cause dynamic instability. In addition, when the mean component is close to Euler’s critical load, the beam would become dynamically unstable regardless of the value of or .
| (a) |
| (b) |
| (c) |
Figure 14 illustrates the first 40 instability regions on the plane (, ) when . It can be seen that the scope of critical frequency ratio is enlarged in general with the increase of under constant fluctuating components. Therefore, mean component growth would make the beam easier to be dynamically unstable.
The impact of or on the dynamic instability region is actually the impact of according to the solving process of dynamic instability regions. Increase (decrease) of or will enhance (reduce) the value of . The point (0.043, 0.35) in the shaded area with smaller than the point C is selected in Figure 11. Modal response of is shown in Figure 15. It can be found that the wind load corresponding to makes the beam become dynamically unstable, but the development trend of instability of is not as obvious as that of C. Therefore, larger mean or fluctuating component of wind load would produce more evident trend of instability.
Another two points (0.043, 0.13) and (0.043, 0.09) are also selected in Figure 11. is located in the unshaded area, while is located in the shaded area. Modal responses corresponding to the two points are shown in Figure 16. It is clear that, even though the mean or fluctuating component corresponding to point is less than that of , the beam is dynamically unstable under the wind load of , while it is stable under the load of . Consequently, when judging whether the structure is dynamically stable or not, the key does not lie in the magnitude of wind load but that the corresponding parametric point falls in dynamic instability regions or not. Because the shape of instability regions under complicated wind force is irregular, the beam might still become dynamically unstable even if the load component is very small. To determine the structural dynamic instability under any wind force, it is required to establish the parametric plane of instability regions and judge the position of the corresponding parametric point.
| (a) |
| (b) |
Furthermore, in order to express vividly the relationship between mean component, fluctuating component, and frequency ratio, three-dimensional parametric space (, , ) can be built. Figure 17 shows the first 4 spatial instability regions. According to the spatial instability regions, dynamic instability of simply supported Euler beam can be determined by any load parameters , and the frequency relationship of wind force and the beam.
6. Conclusions
Through establishing the Mathieu-Hill equation under arbitrary wind force and the computation approach of dynamic instability regions based on the eigenvalue algorithm, a new method for analyzing the dynamic stability of simply supported Euler beam of uniform section under unsteady wind force is established. By introducing the basic load, two-dimensional or three-dimensional parametric coordinate system can be built to determine dynamic instability of the beam and the dynamic instability is verified by calculating wind-induced resonant response. Analysis under the measured wind load shows that any unsteady wind load can be simulated by Fourier series for dynamic instability study and the eigenvalue method achieves sufficiently high precision in computing dynamic instability boundaries. When the parametric point corresponding to specific wind force and structure is located in the instability region, structural response augments divergently. Increase of mean or fluctuating component of wind force would generally widen the scope of instability region and produce more obvious development trend of dynamic instability. However, the magnitude of wind force is not the only factor deciding the dynamic instability state, and the key point is that the magnitude and frequency of wind load should meet a certain relationship with the natural frequency of the structure.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was jointly supported by the National Natural Science Foundation of China (Grant nos. 51222801, 51208126, and 11272095) and the Foundation for Distinguished Young Talents in Higher Education of Guangdong, China (Grant no. 2012LYM_0107).