# Dynamic Stability of Euler Beams under Axial Unsteady Wind Force

Authors: You-Qin Huang; Han-Wen Lu; Ji-Yang Fu; Ai-Rong Liu; Ming Gu
Journal: Mathematical Problems in Engineering
 (2014)

---

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


---

## Body

## 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
(1)
P
(
t
)
=
P
0
+
P
t
cos
⁡
θ
t
.
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 Figure1, whose length is l, Young’s module is E, and the moment of inertia of the section is I. Under the action of axial unsteady wind force f
(
t
), the deflection of the beam can be denoted as [15]
(2)
υ
(
x
,
t
)
=
∑
i
=
1
∞
q
i
(
t
)
sin
i
π
x
l
,
where υ
(
x
,
t
) is the time-varied deflection along the beam, sin
⁡
(
i
π
x
/
l
) is the ith modal shape of the beam, and q
i
(
t
) is the ith modal coordinate.Figure 1
The simply supported Euler beam under wind force.The arbitrary forcef
(
t
) with the duration of t
max
⁡ can be transformed into Fourier series by even prolongation [16]:
(3)
f
(
t
)
=
f
0
+
∑
k
=
1
∞
(
f
k
cos
⁡
k
θ
t
)
,
t
∈
[
0
,
t
max
⁡
]
,
(4)
f
0
=
1
t
max
⁡
∫
0
t
max
⁡
f
(
t
)
d
t
,
(5)
f
k
=
2
t
max
⁡
∫
0
t
max
⁡
f
(
t
)
cos
⁡
k
θ
t
d
t
,
where f
0 and ∑
k
=
1
∞
(
f
k
cos
⁡
k
θ
t
) are the mean and fluctuating component of f
(
t
), respectively. f
k and k
θ are the amplitude and angular frequency of the kth harmonic wave, respectively, and θ
=
π
/
t
max
⁡.Iff
(
t
) is measured by wind force time-history composed of N discrete data from wind tunnel lab or in field, the maximum k is N
-
1.Under the wind force, the Mathieu-Hill equation, corresponding to theith mode of the undamped Euler beam, can be written as [10]
(6)
q
¨
i
+
Ω
i
2
[
1
-
∑
k
=
1
∞
(
2
μ
i
,
k
cos
⁡
k
θ
t
)
]
q
i
=
0
,
where Ω
i is the ith natural frequency of the beam loaded by f
0 and μ
i
,
k is the excitation parameter of the kth harmonic wave for the ith normal mode.Ω
i and μ
i
,
k can be gained by
(7)
Ω
i
=
ω
i
1
-
f
0
f
*
i
,
μ
i
,
k
=
f
k
2
(
f
*
i
-
f
0
)
,
f
*
i
=
i
2
f
*
=
i
2
π
2
E
I
l
2
,
where ω
i is the ith natural frequency of the unloaded beam and f
* 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μ
i
,
k and Ω
i are calculated from the same order of mode, so the distributions of dynamic instability regions in different modes are identical. Consequently, the subscript i in (6)-(7) can be omitted, and all calculations later will use the data corresponding to i
=
1. The Mathieu-Hill equation is rewritten as
(8)
q
¨
+
Ω
2
[
1
-
∑
k
=
1
∞
(
2
μ
k
cos
⁡
k
θ
t
)
]
q
=
0
.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 T or 2
T (T
=
2
t
max
⁡) in (8) is found, the dynamic instability regions can be determined.Periodic solutions of2
T can be written in the form of Fourier series:
(9)
q
(
t
)
=
∑
n
=
1,3
,
5
∞
(
a
n
sin
n
θ
t
2
+
b
n
cos
⁡
n
θ
t
2
)
,
where a
n and b
n are constant coefficients.Substituting (9) into (8) and making the coefficients of terms of sin
⁡
(
n
θ
t
/
2
) or cos
⁡
⁡
(
n
θ
t
/
2
) be zero generate homogeneous system of linear equations with respect to a
n and b
n:
(10)
(
1
+
μ
1
-
θ
2
4
Ω
2
)
a
1
-
(
μ
1
-
μ
2
)
a
3
-
(
μ
2
-
μ
3
)
a
5
+
⋯
=
0
,
-
(
μ
1
-
μ
2
)
a
1
+
(
1
+
μ
3
-
9
θ
2
4
Ω
2
)
a
3
-
(
μ
1
-
μ
4
)
a
5
+
⋯
=
0
,
-
(
μ
2
-
μ
3
)
a
1
-
(
μ
1
-
μ
4
)
a
3
+
(
1
+
μ
5
-
25
θ
2
4
Ω
2
)
a
5
+
⋯
=
0
,
⋮
(
1
-
μ
1
-
θ
2
4
Ω
2
)
b
1
-
(
μ
1
+
μ
2
)
b
3
-
(
μ
2
+
μ
3
)
b
5
+
⋯
=
0
,
-
(
μ
1
+
μ
2
)
b
1
+
(
1
-
μ
3
-
9
θ
2
4
Ω
2
)
b
3
-
(
μ
1
+
μ
4
)
b
5
+
⋯
=
0
,
-
(
μ
2
+
μ
3
)
b
1
-
(
μ
1
+
μ
4
)
b
3
+
(
1
-
μ
5
-
25
θ
2
4
Ω
2
)
b
5
+
⋯
=
0
,
⋮If and only if the coefficient determinant is equal to zero, (10) has nonzero solutions. Therefore, the condition that (8) has periodic solutions of 2
T is that the coefficient determinant of (10) is equal to zero. Considering both cases of a
n and b
n, the critical frequency equation corresponding to 2
T periodic solutions is obtained:(11)
|
1
±
μ
1
-
θ
2
4
Ω
2
-
μ
1
±
μ
2
1
±
μ
3
-
9
θ
2
4
Ω
2
symmetry
-
μ
2
±
μ
3
-
μ
1
±
μ
4
1
±
μ
5
-
25
θ
2
4
Ω
2
⋯
⋯
⋯
⋯
⋯
⋯
-
μ
i
-
1
±
μ
i
⋯
-
μ
i
-
j
±
μ
i
+
j
-
1
⋯
1
±
μ
2
i
-
1
-
(
2
i
-
1
)
2
θ
2
4
Ω
2
⋯
⋯
⋯
⋯
⋯
⋯
⋯
|
=
0
.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θ
*
/
2
Ω, expressed by the ratio of critical load frequency and structural frequency, would be obtained. θ
/
2
Ω is called the frequency ratio, and θ
*
/
2
Ω is called the critical frequency ratio.Similarly, in order to gain the periodic solution condition of the periodT, the periodic solution of period T is written in the following form:
(12)
q
(
t
)
=
b
0
+
∑
n
=
2,4
,
6
∞
(
a
n
sin
n
θ
t
2
+
b
n
cos
⁡
n
θ
t
2
)
.By substituting (12) into (8) and making the coefficients of terms of sin
⁡
(
n
θ
t
/
2
) or cos
⁡
⁡
(
n
θ
t
/
2
) be zero, the critical frequency equation corresponding to the periodic solutions of period T is acquired:(13)
|
1
+
μ
2
-
4
θ
2
4
Ω
2
-
μ
1
+
μ
3
1
+
μ
4
-
16
θ
2
4
Ω
2
symmetry
-
μ
2
+
μ
4
-
μ
1
+
μ
5
1
+
μ
6
-
36
θ
2
4
Ω
2
…
…
…
…
…
…
-
μ
i
-
1
+
μ
i
+
1
…
-
μ
i
-
j
+
μ
i
+
j
…
1
+
μ
2
i
-
(
2
i
)
2
θ
2
4
Ω
2
…
…
…
…
…
…
…
|
=
0
,
(14)
|
1
-
μ
1
-
μ
2
…
-
μ
i
-
1
…
-
2
μ
1
1
-
μ
2
-
4
θ
2
4
Ω
2
-
2
μ
2
-
μ
1
-
μ
3
1
-
μ
4
-
16
θ
2
4
Ω
2
symmetry
…
…
…
…
…
…
-
2
μ
i
-
1
…
-
μ
i
-
j
-
μ
i
+
j
-
2
…
1
-
μ
2
i
-
2
-
(
2
i
-
2
)
2
θ
2
4
Ω
2
…
…
…
…
…
…
…
|
=
0
.By solving (11) with m order determinant, the first m odd dynamic instability regions are obtained. Besides, by solving (13) with m order determinant and (14) with m
+
1 order determinant, the first m even dynamic instability regions are achieved. Thus, the first 2m dynamic instability regions can be determined. Meanwhile, it is known from the above derivation that 2m should not exceed the maximum harmonic number N
-
1.


## 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 A and B can be defined as follows:
(15)
A
=
[
1
±
μ
1
-
μ
1
±
μ
2
1
±
μ
3
symmetry
-
μ
2
±
μ
3
-
μ
1
±
μ
4
1
±
μ
5
⋯
⋯
⋯
⋯
⋯
⋯
-
μ
i
-
1
±
μ
i
⋯
-
μ
i
-
j
±
μ
i
+
j
-
1
⋯
1
±
μ
2
i
-
1
⋯
⋯
⋯
⋯
⋯
⋯
⋯
]
,
B
=
[
1
9
25
⋱
(
2
i
-
1
)
2
⋱
]
.Thus, (11) can be expressed as
(16)
|
A
-
(
θ
2
Ω
)
2
B
|
=
0
.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 A and B of m order, there are 2m eigenvalues θ
*
/
2
Ω, which generate boundary points of m dynamic instability regions. By altering the value of μ
1
,
μ
2
,
…
,
μ
2
m
-
1, the spatial instability region in 2m-dimensional coordinate system (θ
/
2
Ω
,
μ
1
,
μ
2
,
…
,
μ
2
m
-
1) 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 ofT can be determined by solving the critical frequency equations (13) and (14) through the eigenvalue algorithm. For (13), the matrices A and B are
(17)
A
=
[
1
+
μ
2
-
μ
1
+
μ
3
-
μ
2
+
μ
4
…
-
μ
i
-
1
+
μ
i
+
1
…
-
μ
1
+
μ
3
1
+
μ
4
-
μ
1
+
μ
5
…
…
…
-
μ
2
+
μ
4
-
μ
1
+
μ
5
1
+
μ
6
…
-
μ
i
-
j
+
μ
i
+
j
…
…
…
…
…
…
…
-
μ
i
-
1
+
μ
i
+
1
…
-
μ
i
-
j
+
μ
i
+
j
…
1
+
μ
2
i
…
…
…
…
…
…
…
]
,
B
=
[
4
16
36
⋱
(
2
i
)
2
⋱
]
.
The matrices A and B in (14) can be expressed as follows:
(18)
A
=
[
1
-
μ
1
-
μ
2
…
-
μ
i
-
1
…
-
2
μ
1
1
-
μ
2
-
μ
1
-
μ
3
…
…
…
-
2
μ
2
-
μ
1
-
μ
3
1
-
μ
4
…
-
μ
i
-
j
-
μ
i
+
j
-
2
…
…
…
…
…
…
…
-
2
μ
i
-
1
…
-
μ
i
-
j
-
μ
i
+
j
-
2
…
1
-
μ
2
i
-
2
…
…
…
…
…
…
…
]
,
(19)
B
=
[
0
4
16
⋱
(
2
i
-
2
)
2
⋱
]
.


## 4. Calculation of Wind-Induced Dynamic Instability Region and Resonant Response

The value ofm is usually fairly large to meet precision requirement, so it is difficult to graphically display the 2m-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 forcef
(
t
) as expressed in (3), called the initial load here, is firstly normalized as follows:
(20)
f
-
=
max
⁡
k
|
f
k
|
,
(21)
f
′′
(
t
)
=
f
(
t
)
-
f
0
f
-
=
∑
k
=
1
∞
f
k
f
-
cos
⁡
k
θ
t
,
where f
- is the maximum harmonic amplitude of f
(
t
) and f
′
′
(
t
) 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 f
- times less than that of the initial load f
(
t
).Letf
k
′
′ denote the harmonic amplitude of the basic load, and then the basic load can be written as
(22)
f
′
′
(
t
)
=
∑
k
=
1
∞
f
k
′
′
cos
⁡
k
θ
t
.Thus a new wind forcef
′
(
t
) can be generated by defining two parameters α and β to linearly increase the mean and fluctuating components of wind force based on the basic load:
(23)
f
′
(
t
)
=
f
0
′
+
f
k
m
′
f
′
′
(
t
)
=
α
f
*
+
β
f
*
f
′
′
(
t
)
=
α
f
*
+
∑
k
=
1
∞
β
f
*
f
k
′
′
cos
⁡
k
θ
t
=
α
f
*
+
∑
k
=
1
∞
β
f
*
f
-
f
k
cos
⁡
k
θ
t
,
α
=
f
0
′
f
*
,
β
=
f
k
m
′
f
*
,
where f
′
(
t
) is the newly generated wind force; f
0
′ and f
k
m
′ are the mean component and the maximum harmonic amplitude of f
′
(
t
), respectively; α and β are the ratios of f
0
′ and f
k
m
′ 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α
f
* and its fluctuating component or harmonic amplitude is β
f
* times that of the basic load, while the harmonic frequency composition is identical to the initial load or the basic load.f
′
(
t
) is further simplified as follows:
(24)
f
′
(
t
)
=
f
0
′
+
∑
k
=
1
∞
f
k
′
cos
⁡
k
θ
t
,
where f
k
′ denotes the harmonic amplitude of f
′
(
t
).Meanwhile, the natural frequency and excitation parameter corresponding tof
′
(
t
) can be computed by
(25)
Ω
′
=
ω
1
-
f
0
′
f
*
=
ω
1
-
α
f
*
f
*
=
ω
1
-
α
,
μ
k
′
=
f
k
′
2
(
f
*
-
f
0
′
)
=
(
β
f
*
/
f
-
)
f
k
2
(
f
*
-
α
f
*
)
=
(
β
/
f
-
)
f
k
2
(
1
-
α
)
=
β
2
(
1
-
α
)
f
k
′
′
.The expression of excitation parameter can be further simplified by defining(26)
μ
′
=
β
2
(
1
-
α
)
,
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:(27)
μ
k
′
=
μ
′
f
k
′
′
.Moreover, the Mathieu-Hill equation consistent withf
′
(
t
) is
(28)
q
¨
+
Ω
2
[
1
-
2
μ
′
∑
k
=
1
∞
(
f
k
′
′
cos
⁡
k
θ
t
)
]
q
=
0
.That is(29)
q
¨
+
Ω
2
[
1
-
2
μ
′
f
′
′
(
t
)
]
q
=
0
.Therefore, variousμ
′ can be obtained only by adjusting α or β, and then the two-dimensional parametric plane (
θ
/
2
Ω
,
μ
′
) or three-dimensional parametric space (
θ
/
2
Ω
,
α
,
β
) 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 α
∈
[
0,1
] and β
∈
[
0,1
].When the parametric point(
θ
/
2
Ω
,
μ
′
) is located in the instability region, the beam is deemed to be dynamically unstable. In order to verify such a conclusion, modal response q
(
t
) of the beam can be calculated through the Mathieu-Hill equation (29). q
(
t
) should keep increasing with time when the beam is unstable.The Mathieu-Hill equation shown in (29) can be converted to the following form:
(30)
x
˙
=
-
Ω
2
[
1
-
2
μ
′
f
′
′
(
t
)
]
q
,
q
˙
=
x
.
Let
(31)
y
=
{
x
q
}
.
Equation (29) is expressed in the vector type
(32)
y
˙
=
f
(
y
,
t
)
.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 q
(
t
) at the nth iteration step can be obtained by the following iteration formula:
(33)
y
n
=
y
n
-
1
+
h
6
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
,
k
1
=
f
(
t
n
-
1
,
y
n
-
1
)
,
k
2
=
f
(
t
n
-
1
+
h
2
,
y
n
-
1
+
h
2
k
1
)
,
k
3
=
f
(
t
n
-
1
+
h
2
,
y
n
-
1
+
h
2
k
2
)
,
k
4
=
f
(
t
n
-
1
+
h
,
y
n
-
1
+
h
k
3
)
.The time step length is set ash
=
0.002 s and initial disturbance of response is q
(
t
=
0
)
=
0.001 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 (Figure2). 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 f
(
t
) on the beam, as shown in Figure 3.Figure 2
Rigid model in the wind tunnel tests.Figure 3
Measured wind force acting on the beam.The wind load data number isN
=
6000, and the sampling frequency is F
s
=
9.58 Hz, so the load duration is t
max
⁡
=
N
/
F
s
=
626.3 s. The load is simulated by the periodic load f
(
t
) with period 2
t
max
⁡ according to (3)–(5). Comparison of f
(
t
) 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 f
0
=
4.19
e
3 N and the maximum harmonic amplitude is f
-
=
288.87 N by (4), (5), and (20). The relationship between harmonic amplitude f
k and corresponding harmonic frequency k
/
(
2
t
max
⁡
) is shown in Figure 5, which indicates that larger amplitudes correspond to smaller frequencies. Moreover, the natural frequency of the beam is defined as ω
=
5 Hz.Comparison of measured and fitted wind forces.
(a)
Load duration =t
max
⁡
(b)
Load duration =t
max
⁡/10Figure 5
Relationship between the harmonic amplitude and frequency in the fitted wind force.Normalizing the initial load through (20)–(22) produces the basic load f
′
′
(
t
), 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 f
- times, while the harmonic frequencies are invariant.Time-history and frequency feature of the basic load.
(a)
Load time-history
(b)
Relationship between the amplitude and frequencyNew 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 m is equal to 10, 40, and 80, respectively. It is obvious that the boundary values vary with m and larger m results in close boundaries with higher precision. Since the results under m
=
40 and m
=
80 are close, m is set for 40 in the following calculations. Meanwhile, for convenience, only the first 100 data of wind force are selected for later computations (N
=
100).Figure 7
Comparison of the first instability regions under differentm’s.Bolotin [10] has pointed out that, if diagonal elements of the critical frequency equation are ignored, the boundary of the kth instability region under arbitrary dynamic load can be calculated by the following approximate formula:
(34)
θ
*
2
Ω
≈
1
k
1
±
μ
k
.Equation (34) indicates that the kth instability region only depends on the kth 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 ith harmonic wave on the width of the kth instability region is the order of (
μ
k
±
μ
i
)
2 [10], while the value of (
μ
k
±
μ
i
)
2 under wind load could be fairly large. Therefore, the interaction between different harmonic waves cannot be ignored for wind-induced dynamic instability analysis.Figure 8
Comparison of the first instability regions by the approximate formula and eigenvalue method.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 Figure8 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 q
(
t
) 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.Modal responseq
(
t
) of the parametric points A and B.
(a)
Point A
(b)
Point BDistribution of the first 40 dynamic instability regions (the first 20 odd regions corresponding to period2
T and the first 20 even regions corresponding to period T) is shown as the shaded area in Figure 10. When the parametric point (
θ
/
2
Ω
,
μ
′
) 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 θ
/
2
Ω is fixed.Figure 10
Distribution of the first 40 instability regions on the parametric plane(
θ
/
2
Ω
,
μ
′
).Whenα
=
β
=
0.5 is assumed, we have μ
′
=
0.5 according to (26). By t
max
⁡
=
100
/
9.58
=
10.44 s, we get θ
=
π
/
t
max
⁡
=
0.3 Hz and θ
/
2
Ω
=
0.043. 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 θ
/
2
Ω between 0.02 and 0.065. The modal response q
(
t
) 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.Figure 11
Location of the parametric points on the parametric plane(
θ
/
2
Ω
,
μ
′
).Figure 12
Modal responseq
(
t
) of the parametric point C.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α, β, t
max
⁡, 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 (θ
/
2
Ω, α) and (θ
/
2
Ω, β) can be built. Figure 13 provides the first 40 instability regions on the plane (θ
/
2
Ω, β) when α
=
0, 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 θ.The first 40 instability regions on the parametric plane(
θ
/
2
Ω
,
β
) under different α’s.
(a)
α
=
0
(b)
α
=
0.5
(c)
α
=
0.99Figure14 illustrates the first 40 instability regions on the plane (θ
/
2
Ω, α) when β
=
0.5. 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.Figure 14
The first 40 instability regions on the parametric plane(
θ
/
2
Ω
,
α
) when β
=
0.5.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 C
1 (0.043, 0.35) in the shaded area with smaller μ
′ than the point C is selected in Figure 11. Modal response of C
1 is shown in Figure 15. It can be found that the wind load corresponding to C
1 makes the beam become dynamically unstable, but the development trend of instability of C
1 is not as obvious as that of C. Therefore, larger mean or fluctuating component of wind load would produce more evident trend of instability.Figure 15
Modal responseq
(
t
) of the parametric point C
1.Another two pointsC
2(0.043, 0.13) and C
3(0.043, 0.09) are also selected in Figure 11. C
2 is located in the unshaded area, while C
3 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 C
3 is less than that of C
2, the beam is dynamically unstable under the wind load of C
3, while it is stable under the load of C
2. 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.Modal responseq
(
t
) of the parametric points C
2 and C
3.
(a)
C
2
(b)
C
3Furthermore, in order to express vividly the relationship between mean component, fluctuating component, and frequency ratio, three-dimensional parametric space (θ
/
2
Ω, α, β) 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 θ
/
2
Ω of wind force and the beam.Figure 17
Distribution of the first 4 instability regions on the parametric space(
θ
/
2
Ω
,
α
,
β
).


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



---

*Source: 434868-2014-03-23.xml*