Limit Cycle Bifurcations In The In-plane Galloping Of Iced

Preparing to load PDF file. please wait...

0 of 0
Limit Cycle Bifurcations In The In-plane Galloping Of Iced

Transcript Of Limit Cycle Bifurcations In The In-plane Galloping Of Iced

Journal of Applied Analysis and Computation Volume 10, Number 4, August 2020, 1355–1374

Website: DOI:10.11948/20190202

Peng Liu1,2,3,†, Anqi Zhou3,4, Bing Huo5 and Xijun Liu3,4
Abstract In this paper, we establish a mathematical model to describe inplane galloping of iced transmission line with geometrical and aerodynamical nonlinearities using Hamilton principle. After Galerkin Discretization, we get a two-dimensional ordinary differential equations system, further, a near Hamiltonian system is obtained by rescaling. By calculating the coefficients of the first order Melnikov function or the Abelian integral of the near-Hamiltonian system, the number of limit cycles and their locations are obtained. We demonstrate that this model can have at least 3 limit cycles in some wind velocity. Moreover, some numerical simulations are conducted to verify the theoretical results.
Keywords Galloping, limit cycle, bifurcation, Melnikov function.
MSC(2010) 37G15, 34C07.
1. Introduction
Galloping of iced transmission line is a motion involving low frequency, large amplitude and self-excited properties due to the instability of the aerodynamic forces caused by wind flow acting on the non-circular section. The classic mechanisms in galloping researches include vertical galloping mechanism by Den. Hartog [2], torsional galloping mechanism by Nigol and Buchan [17, 18], and eccentric inertia mechanism proposed by Yu et al. [22].
On the basis of above mechanism, many other galloping models were proposed. Fuhao Liu et al. established a two-degree-of-freedom model of iced, electrical quad
†the corresponding author. Email address:[email protected](P. Liu) 1Shanxi Key Laboratory of Material Strength & Structural Impact, Taiyuan
University of Technology, 79 Yingze West Street, 030024 Taiyuan, China 2National Demonstration Center for Experimental Mechanics Education,
Taiyuan University of Technology, 79 Yingze West Street, 030024 Taiyuan, China 3School of Mechanical Engineering, Tianjin University, 135 Yaguan Road, 300354 Tianjin, China 4Tianjin Key Lab Nonlinear Dynamics and Chaos Control, Tianjin University, 135 Yaguan Road, 300354 Tianjin, China 5College of Mechanical Engineering, Tianjin University of Science and Technology, 1038 Dagu South Road, 300222 Tianjin, China ∗The authors were supported by National Natural Science Foundation of China (51808389), Project of Shanxi Youth Foundation (201901D211067), Project of Tianjin Youth Foundation (18JCQNJC08000), and the Scientific and Technological Innovation Programs of Higher Education Institutions in Shanxi (2019L0254).


P. Liu, A. Zhou, B. Huo & X. Liu

bundle conductor [10]. They analyzed co-dimension-2 bifurcation in galloping using centre manifold and invertible linear transformation. Luongo et al. proposed a new model based on curved-beam theory [13,14]. Branches of periodic solutions were obtained by employing complementary solution methods, and their stability evaluated as functions of wind velocity. Lei hou et al built a two-degree-of-freedom differential equations containing lateral and torsional motion by Lagrange equation [6]. They investigated the existence of chaotic motion in galloping system and found that there coexisted three resonance patterns in system, the mutual transformation of which may lead to chaotic motion. Zhang et al. [23] developed a useful tool for a bundle conductor of an electrical transmission line by using a three-degree-of-freedom hybrid model. This model accommodated interactions of the vertical, horizontal and torsional motions, non-linear aerodynamic loads, a non-uniform ice geometry and a variation of the wind along a span. Our team proposed a three-degree-of-freedom model to describe the nonlinear interactions between the in-plane, out-of-plane and torsional vibration [11]. Eigenvalue analysis was applied to determine the switch of different mode such as mono-modal, bi-modal and multi-modal galloping. Pierre et al. [15] explored a new approach to analyze galloping for a minimal aspect ratio of accretion shapes, by including in the quasi-static analysis the effect of the torsional vibration on the lift force. Zhaohong Qin applied Lyapunov stability theory get the critical wind speed, and obtained the effects of structural parameters to the critical wind speed [19]. Bin Liu [8] built a continuum model of the conductor galloping with D-shape ice and the aerodynamic forces were described by using the quasi steady hypothesis, where the aerodynamic coefficients were expanded by the polynomial curves with a third order and a ninth order respectively. Zhimiao Yan [21] formulated a nonlinear galloping model with regard of bending, rotation and eccentricity of cross section and the bifurcation and stability of the two cases (1:1 resonance and 2:1 resonance) were analyzed. In view of extensive literature review, we find that most of papers in galloping area focus and study equilibrium stability, internal resonance vibration, center manifold, bifurcation analysis, and even chaos, however, the limit cycle bifurcations are seldom investigated in galloping of ice transmission line. Based on our previous research, limit cycle bifurcation may be a useful tool to explore the galloping problem. Hence, the objective of the present paper is to study limit cycle bifurcations in galloping of iced transmission line. Meanwhile, our work may benefit the application of limit cycle bifurcations in engineering problem.
Limit cycle is a common phenomenon in science and engineering, and was first discovered by Poincar´e. Limit cycles are generated through bifurcations in many different ways, i.e., through Hopf bifurcation from a center or a focus, via Poincar´e bifurcation from closed orbits, or separatrix from homoclinic or heteroclinic orbits [4]. The most well-known 23 problems related limit cycles were presented by Hilbert [5], among which Hilbert’s 16th problem is still open up to now. The second part of Hilbert’s 16th problem is to find the maximal number of limit cycles and their relative locations in the planar polynomial system of degree n,

x˙ = Pn(x, y), y˙ = Qn(x, y),
where Pn(x, y) and Qn(x, y) are nth-degree polynomials in x and y. Although many research results have been obtained, the precise number of limit cycles is not obtained. Then the weakened Hilbert’s 16th problem was presented by Arnold [1] to help understand and attack the problem. The weakened problem is to find the

Limit cycle bifurcations in galloping


maximum number of isolated zeros of the Abelian integral or Melnikov function,

M (h, δ) =

Qn(x, y)dx − Pn(x, y)dy,

H (x,y)=h

where H(x, y) is a polynomial in x and y. A contrast on galloping models with different degree-of-freedoms was conducted
by Bin Liu [9], the result showed that the one-degree-of-freedom model can also present the trend of galloping amplitude change from the point view of qualitative analysis. Hence, in this paper, we propose an one-degree-of-freedom model, which describes the in-plane motion.
The rest of paper is organized as follows. Based on Hamilton principle, the mathematical model is derived in Section 2 to describe the in-plane nonlinear vibration. Galerkin method is employed to spatially disperse the established model to get an ordinary equation system. Then using rescaling, we get a near-Hamiltonian system. The unperturbed Hamiltonian system has two centers and one saddle. Some preliminaries on limit cycle bifurcations are given in Section 3. In Section 4, we get the expansions of Melnikov functions of near-Hamiltonian system, then apply the theorem in Section 3 to obtain the number of limit cycles in present model. To verify the theoretical results, we conduct some numerical simulations in Section 5. Finally, a conclusion is drawn in Section 6.

2. Formulation of the system
The transmission line is modeled as a body made of flexible cable with length l, before modeling, we have following assumptions:
(i) The thin ice accretion is in crescent shape, and assumed to be uniform along the transmission line. The transmission line is supposed to bear tension only but not resist compression and bending moment;
(ii) The transmission tower is considered to be rigid so that the ends of the transmission line are fixed, and the sag-to-span ratio of transmission line is small, the expression of transmission line satisfies catenary equation [16];
(iii) We only consider the in-plane dynamic, which is on the basis of Den. Hartog galloping mechanism, while the axial vibration, out-of-plane vibration, and torsional vibration is ignored;
(iv) Based on the quasi-steady analysis, wind forces are evaluated at constant wind speed and angle of attack, and these wind forces are employed to galloping system.
Under above assumptions and considering the geometrical and aerodynamical nonlinearities, we establish the galloping model, the schematic of the transmission line is shown in Figure 1 (a). The transmission line is placed on the initial configuration Γ0 under the gravity. Γ¯ is the reference configuration, which swings the angle φ from Γ0 due to the static aerodynamic forces act on the transmission line at the time t = 0. Γ denotes the displaced configuration at the time t > 0, and v denotes the in-plane dynamic displacements at the time t. A wind flow blows normally to the plane where Γ locates, and the aerodynamic forces produced by wind flow act on the crescent-shaped uniformly along the transmission line.


P. Liu, A. Zhou, B. Huo & X. Liu





¡(x0, y0+v) A'




B ds










Figure 1. The theoretical model of transmission line: (a) configuration, and (b) schematic diagram of dynamic displacement.

A segment of the transmission line of an infinitesimal length dx is considered and its dynamic displacement is illustrated in Figure 1 (b), from which we have

# »

#» dy0 #»

AB = dx i + dx j ,


where y0 is the catenary equation [16] as follows


y0 = − 2T0 sinh mgx sinh mg(l − x) ,





where T0 is the initial tension of the transmission line, m is the mass per unit length of the iced transmission line, and g is gravitational acceleration constant , hereinafter, we denote y0x = ddyx0 .
The arc length of the segment in the reference position is obtained in the fol-
lowing equation

# » ds0 = AB =

1 + y2 dx ≈ 1 + 1 y2 − 1 y4 .


2 0x 8 0x

The vector of the deformed segment is expressed as

A# ′B»′ = dx #i» + (y0x + vx) dx #j»,

(2.3) (2.4)

Limit cycle bifurcations in galloping

y FL









Figure 2. section.

The schematic of velocity vector and the aerodynamics forces acting on the conductor

where vx = ddxv , similarly, the arc length of the deformed segment is written ds = A# ′B»′ = 1 + (y0x + vx)2dx.


Due to the small sag-to-span ratio, we use Taylor expansion up to the third order to get the strain expression of segment, given by

ε = dsd−sd0s0 =

√ 1+(y0x+vx)2 − 1


≈ (y0x − y03x)vx + ( 12 − y02x)vx2 − 12 y0xvx3.


Then the potential energy of the transmission line is given by



V = 0 T0 + 2 EAε (ds − ds0)


where EA is the in-plane stiffness. The schematic of velocity vector and the aerodynamics forces acting on the conductor section is shown in Figure 2. Then the kinetic energy of the iced transmission line is expressed as

T = l 1 mv˙ 2ds0, 02
and the virtual work performed by all forces is given


W = Fyvds0,


where Fy is the aerodynamic force acting on the in-plane direction of the transmission line [22], represented as

Fy = FL cos α1 − FD sin α1,



P. Liu, A. Zhou, B. Huo & X. Liu

in which FL and FD are lift and drag forces due to wind, α1 is angle related to vertical velocity. The expression of FL and FD are:

FL = 12 ρUr2DCL, FD = 12 ρUr2DCD,


where CL and CD are lift and drag coefficients, ρ is the air density, D is th√e diameter of the bare conductor, Ur is the relative wind velocity given by Ur = U 2 + v˙2, where U is the mean wind velocity and treated as a bifurcation parameter. Hence,
Fy in (2.10) can be expressed as

Fy = 12 ρUr2DCy,


in which Cy is the pure aerodynamic coefficient acting on the in-plane of the trans-

mission. The coefficient Cy can be written as a function of the dynamical angle of

attack α, i.e.,

Cy = ry1 α + ry2 α2 + ry3 α3,


where coefficients ryi (i = 1, 2, 3) are obtained by curve-fitting from the experimental data in [23] and shown in Table 1, α is determined by α = θ0 − α1, where θ0 is the initial equilibrium angle of ice. For small angles, α can be expressed as α = θ0 − Uv˙ .
Now, applying the Hamilton principle

δ (T − V + W )dt = 0,


we get the partial differential equation (PDE) of galloping of iced transmission line

mv¨ 1 + 21 y02x − 18 y04x + 34 EA(15y02x − 2)vx2vxx + 3T0 − 136 EAy0x 16 − 40y02x − y04x vxvxx + − 116 EA 16y02x + 9y04x + 12 T0 3y02x − 2 vxx −Fy 1 + 21 y02x − 18 y04x = 0.


where the dot denotes differentiation with respect to t. In order to get the ordinary differential equation (ODE) of the model, we use
Galerkin Discretization method

π v(x, t) = V (t) sin( x),


where V (t) is the dynamic displacements for the in-plane mode, and sin( πl x) is the trail function.
Then substituting (2.16) into (2.15), and applying Galerkin procedure to get the following ODE model

V¨ + a1V + 2ξωV˙ = a2V 2 + a3V 3 + b1V˙ + b2V˙ 2 + b3V˙ 3 + b4V˙ 4 + b5V˙ 5, (2.17)

where ω and ξ is the natural frequency and the damping ratio of the in-plane vibration, respectively, ai (i = 1, 2, 3) and bi (i = 1, . . . , 5) are the integral constants

Limit cycle bifurcations in galloping


Table 1. Conductor data and curve-fitting coefficients of Cy [23]

Parameter Definition



Mass per unit length



Cross-sectional area of iced transmission line



Diameter of bare conductor



Elastic modulus

4.78 × 1010


Length of the transmission line



Air density



Initial tension



Initial equilibrium angle of ice



Damping ratio



Fitting coefficient of Cy


Fitting coefficient of Cy


Fitting coefficient of Cy

-0.1667 -4.0547 8.3581

kg/m mm2
m N/m2
m kg/m3
N degree

related to physical parameters (see in Table 1) of transmission line, for reading
convenience, the expression of ai (i = 1, 2, 3) and bi (i = 1, . . . , 5) are given in Appendix. From the expression of bi (i = 1, . . . , 5) in Appendix, we can find that b1, b3, b4, b5 are corresponding to wind velocity U , except for b2.
Now we introduce a transformation y = x˙ to the model (2.17) to get a second-
order non-autonomous system, determined by

  x˙ = y,  y˙ = −a1x + a2x2 + a3x3 − 2ξωy  +b1y + b2y2 + b3y3 + b4y4 + b5y5.


For convenience of analysis, we simplify the system (2.18) by the following rescalings:

x → µ1x, y → µ2y, t → µ3t,


where µ1 = aa12 , µ2 = aa1322 , and µ3 = √a1. By the above, the model (2.18) is transformed to the following equivalent one:

  x˙ = y,  y˙ = −x + x2 + ax3 + b∗01y + b∗02y2 + b∗03y3 + b∗04y4 + b∗05y5,



a = a3µ31 ,
µ2 µ3

b∗ = b3µ22 ,



b∗01 = µb13 − 2ξ, b∗02 = b2µµ32 ,

b∗ = b4µ32 ,



b∗ = b5µ42 .




Further, if we consider low velocity wind (U ≤ 20m/s), then, according to (2.21) and Table 1, the coefficients b∗0i (i = 1, ..., 5) can be written as b∗0i = ε′b0i (i = 1, ..., 5), ε′ ≪ 1, then the system (2.20) becomes

  x˙ = y,  y˙ = −x + x2 + ax3 + ε′ b01y + b02y2 + b03y3 + b04y4 + b05y5 ,



P. Liu, A. Zhou, B. Huo & X. Liu

xc xs

Figure 3. The phase portrait of system (2.23)

which is a near-Hamiltonian system. When ε′ = 0, the system becomes a planar Hamiltonian system as follows

  x˙ = y,  y˙ = −x + x2 + ax3,


with the Hamiltonian function

H(x, y) = 1 y2 + 1 x2 − 1 x3 − a x4. 2234


Taking parameter values from Table 1, the unperturbed system (2.23) has two
elementary centers (0, 0) and (xc, 0), a double homoclinic loop L = L1 ∪ L2 passing through a hyperbolic saddle (xs, 0), where xc and xs are the roots of −x+x2 +ax3 = 0, and L1 ≡ L|xxs . The phase portrait of system (2.23) is shown in Figure 3.
Let hc1 ≡ H(0, 0) = 0, hc2 ≡ H(xc, 0) = 0.039412571 and hs ≡ H(xs, 0) = 0.2838521430, we can easily get 0 ≤ hc1 < hc2 < hs. In unperturbed system (2.23), there are three families of periodic orbits denoted by Lh, L1h and L2h, respectively, where

Lh : H(x, y) = h, h ∈ (hs, +∞); L1h : H(x, y) = h, h ∈ (hc1 , hs); L2h : H(x, y) = h, h ∈ (hc2 , hs).

Limit cycle bifurcations in galloping


3. Preliminaries

To investigate the limit cycles bifurcation in near-Hamiltonian system (2.22), we give some preliminaries in this section. Consider a C∞ system of form

x˙ = Hy + εp(x, y, δ), y˙ = −Hx + εq(x, y, δ),



H(x, y),

p(x, y, δ) =

aij xiyj ,


q(x, y, δ) =

bij xiyj ,


are C∞ functions, ε ≥ 0 is small and δ ∈ D ⊂ Rn is a vector parameter with D compact. When ε = 0, the unperturbed system of (3.1) is

x˙ = Hy, y˙ = −Hx.


Suppose the system (3.2) has a singular point, without loss of generality, we can assume it is at origin. Then, if the singular point at the origin is elementary, we may assume

H(x, y) = 1 y2 + h20x2 +

hij xiyj ,



h20 ̸= 0.

We assume the above system has a family of periodic orbits given by Lh = {H(x, y) = h, h ∈ (α, β)}. The first order Melnikov function or the Abelian integral of system (3.1) is expressed as follows

M (h, δ) = qdx − pdy =

(px + qy)dxdy, h ∈ (α, β),
H ≤h


the number of isolated zero roots of which relates to the number of limit cycles of

sytem (3.1). For convenience, we denote px + qy =

cij xiyj .


Now we consider two different cases of singular points and give the corresponding

expansion of the above Melnikov function.

• Case 1: We suppose the system (3.2) has an elementary center at the origin, and without loss of generality, we may suppose the expression of Hamiltonian function near the orgin is given by

H(x, y) = ω (x2 + y2) +

hij xiyj ,



ω > 0.


Then we have


P. Liu, A. Zhou, B. Huo & X. Liu

Theorem 3.1 ( [3]). Let (3.4) hold. Then M (h, δ) is C∞ in 0 ≤ h ≪ 1 with

M (h, δ) = h l≥0 bl(δ)hl+1 + O(hl+1)


formally for 0 ≤ h ≪ 1. The expression of the coefficient bl is given in [7], as fol lows

b0 = 2πc00,


= c00π

15 2

h230 + h203

+ 23

h221 + h212 + 2h30h12 + 2h21h03

−c10π (h12 + 3h30) − c01π (h21 + 3h03) + c20π + c02π,



where c00, c10, c01, c20 and c02 are determined by the coefficient of px + qy = cij xiyj .
• Case 2: Suppose the system (3.2) has a double homoclinic loop L = L1 ∪ L2 passing through the origin which is a hyperbolic saddle, where L1 = L|x≤0 and L2 = L|x≥0. We give notation Lh to denote a family of periodic orbits defined by H(x, y) = h for 0 < h ≪ 1 and Li(h) (i = 1, 2) to denote two families of periodic orbits determined by H(x, y) = h for 0 < −h ≪ 1. Then, without loss of generality, we may assume

H(x, y) = λ (y2 − x2) +

hij xiyj ,



λ > 0.


Then we have Theorem 3.2 ( [3]). Let (3.7) hold. Then

M (h, δ) = c0(δ) + 2c1(δ)h ln |h| + c2(δ)h + 2c3(δ)h2 ln |h| + O(h2), Mi(h, δ) = c0i(δ) + c1(δ)h ln |h| + c2i(δ)h + c3(δ)h2 ln |h| + O(h2),
i = 1, 2,



c0i(δ) = Li qdx − pdy, i = 1, 2,

c0(δ) = c01(δ) + c02(δ),




+b01 λ


c2i(δ) = Li (px + qy − a10 − b01)dt + bic1,

c2(δ) = c21(δ) + c22(δ),

c3(δ) = . . .

for some constants b, b1 and b2.

i = 1, 2,
ModelExpressionLimit Cycle BifurcationsLiuLimit Cycles