LIMIT CYCLE BIFURCATIONS IN THE IN-PLANE GALLOPING OF ICED TRANSMISSION LINE∗

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.


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 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é. Limit cycles are generated through bifurcations in many different ways, i.e., through Hopf bifurcation from a center or a focus, via Poincaré 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, where P n (x, y) and Q n (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 maximum number of isolated zeros of the Abelian integral or Melnikov function, 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.

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.  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 where y 0 is the catenary equation [16] as follows where T 0 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 y 0x = dy0 dx . The arc length of the segment in the reference position is obtained in the following equation 3) The vector of the deformed segment is expressed as where v x = dv dx , similarly, the arc length of the deformed segment is written 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 Then the potential energy of the transmission line is given by 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 and the virtual work performed by all forces is given where F y is the aerodynamic force acting on the in-plane direction of the transmission line [22], represented as in which F L and F D are lift and drag forces due to wind, α 1 is angle related to vertical velocity. The expression of F L and F D are: where C L and C D are lift and drag coefficients, ρ is the air density, D is the diameter of the bare conductor, U r is the relative wind velocity given by U r = √ U 2 +v 2 , where U is the mean wind velocity and treated as a bifurcation parameter. Hence, F y in (2.10) can be expressed as 12) in which C y is the pure aerodynamic coefficient acting on the in-plane of the transmission. The coefficient C y can be written as a function of the dynamical angle of attack α, i.e., where coefficients r yi (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 we get the partial differential equation (PDE) of galloping of iced transmission line where the dot denotes differentiation with respect to t. In order to get the ordinary differential equation (ODE) of the model, we use 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 where ω and ξ is the natural frequency and the damping ratio of the in-plane vibration, respectively, a i (i = 1, 2, 3) and b i (i = 1, . . . , 5) are the integral constants Initial equilibrium angle of ice 40 degree ξ Damping ratio 0.0515 ry 1 Fitting coefficient of Cy -0.1667 ry 2 Fitting coefficient of Cy -4.0547 ry 3 Fitting coefficient of Cy 8.3581 related to physical parameters (see in Table 1) of transmission line, for reading convenience, the expression of Now we introduce a transformation y =ẋ to the model (2.17) to get a secondorder non-autonomous system, determined by For convenience of analysis, we simplify the system (2.18) by the following rescalings: a2 , and µ 3 = √ a 1 .

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 forṁ where are C ∞ functions, ε ≥ 0 is small and δ ∈ D ⊂ R n is a vector parameter with D compact. When ε = 0, the unperturbed system of (3.1) iṡ (3.2) 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 We assume the above system has a family of periodic orbits given by L h = {H(x, y) = h, h ∈ (α, β)}. The first order Melnikov function or the Abelian integral of system (3.1) is expressed as follows the number of isolated zero roots of which relates to the number of limit cycles of sytem (3.1). For convenience, we denote p x + q y = i+j≥0 c ij x i y j . 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 Then we have
• Case 2: Suppose the system (3.2) has a double homoclinic loop L = L 1 ∪ L 2 passing through the origin which is a hyperbolic saddle, where L 1 = L| x≤0 and L 2 = L| x≥0 . We give notation L h to denote a family of periodic orbits defined by H(x, y) = h for 0 < h ≪ 1 and L i (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 Then we have

Expansions of Melnikov functions of system (2.22)
In this section, we give the expansion of the related Melnikov function corresponding to system (2.22). First, we recall the results obtained in Section 2, the unperturbed system (  Here we need to note that we only give the formulas of b 1,0 and b 1,1 , since we have mentioned in section 2 that we only have one bifurcation parameter U in this system. To obtain the coefficients b 2,j (δ) in (4.1), we introduce a transformation and a time scaling t → √μ The Hamiltonian function of system (4.5)| ε=0 is For the coefficient c 1 (δ) in (4.2) and (4.3), we make a change of variables of the form in system (2.22) (4.5) The corresponding Hamiltonian function of system (4.5)| ε=0 is (1 + 3ax s ) andh 40 = − a 4μ 2 . Now using the formula in Theorem 3.2, we obtain To get the coefficients c 0 (δ) and c 0i (δ), i = 1, 2 in (4.2) and (4.3), we apply the formulas in Theorem 3.2 to obtain whereũ − is the intersection point between homoclinic loop L 1 andũ-axis,ũ + is the intersection point between homoclinic loop L 2 andũ-axis, and f (ũ) = 1 − 2h 30ũ − 2h 40ũ 2 .

Number of limit cycles of system (2.22)
So far, we have obtained all coefficients in (4.1), (4.2) and (4.3). Now we use these coefficients to study limit cycles. For convenience, we calculate I 1j and I 2j , j = 1, . . . , 5 at the beginning, to get , substituting the parameter values in Table 1 into which, we obtainŪ = 3.787861482.  5 and h 6 = h s + 100 we have Applying the method in [3], we obtain that there exist 1 limit cycle between homoclinic loop L 1 and center (0, 0), and one limit cycle surrounding the double homoclinic loop L.
Since b 2,0 = b 1,0 = 0, similar as [3], taking b 1,0 or b 2,0 as a free parameter, we can know that there exists one more limit cycle near center (0, 0) or (x c , 0). Hence, we have following theorem. Theorem 4.1. There exists some U nearŪ = 3.787861482 such that the system (2.22) has at least 3 limit cycles: 1 limit cycle is near one of the center, 1 limit cycle is between homoclinic loop L 1 and center (0, 0), and last limit cycle surrounds the double homoclinic loop L. Now solving the equation c 01 = 0 yieldsŨ 1 = 19.56011461, substituting which into coefficients b 1,0 , b 2,0 , c 1 , c 02 and c 0  which means there exists one limit cycle near the homoclinic loop L 1 . To examine the limit cycle between homoclinic loop L 1 and center (0, 0), between homoclinic loop L 2 and center (x c , 0), and surrounding double homoclinic loop L, taking U = U 1 = 19.56011461 and some ε 1 , ε 2 , ε 3 , ε 4 , ε 5 positive and sufficiently small, and  and take some ε 1 , ε 2 , ε 3 , ε 4 , ε 5 positive and sufficiently small, and Employing the method in [3], we have following theorem.

Theorem 4.3.
There exists some U nearŨ 2 = 4.08330378 such that the system (2.22) has at least 3 limit cycles: 1 limit cycle is near homoclinic loop L 2 , 1 limit cycle is between homoclinic loop L 1 and center (0, 0), and 1 limit cycle surrounds the double homoclinic loop L.  Applying the method in [3], we have following theorem.

Theorem 4.4.
There exists some U nearŨ 3 = 181.8652959 such that the system (2.22) has at least 2 limit cycles: 1 limit cycle near double homoclinic loop L, and the limit cycle locates between homoclinic loop L 1 and center (0, 0).

Numerical simulation
In this section, we present some numerical examples and simulations to demonstrate the theoretical results obtain in previous section, with physical parameter values taken from Table 1. As we mentioned in Section 2, we only consider the low wind velocity U ≤ 20m/s. Taking U = 3.9 and the initial value x(0) = 0.08, y(0) = 0, we give the time evolution of x and the phase portrait of system (2.22) in Figure 4. Figure 4 shows that there indeed exists a limit cycle near center (0, 0), and we can demonstrate the limit cycle is stable based on the time evolution and phase portrait. The amplitude of limit cycle is approximately 0.06, timed by dimensionless scale µ 1 in (2.19), which is equal to about 0.22 m in reality.
When we select U = 3.9 and the initial value at x(0) = 2.95, y(0) = 0, the simulation of system (2.22) is shown in Figure 5. It shows a good agreement with  theoretical prediction, confirming that there is one limit cycle near center (x c , 0), similarly, from the time evolution and phase portrait, we can determine that it is stable. The amplitude of limit cycle is also approximately 0.06, multiplied by dimensionless scale µ 1 in (2.19), which is equal to about 0.22 m. Further, we conduct a numerical simulation on original system (2.17). When U = 3.9, the steady-state time history is shown in Figure 6, from which we can see that the galloping amplitude is nearly 0.22 m. It is in agreement with the theoretical prediction.
At last, taking U = 19.55 and initial value x(0) = −0.62, y(0) = 0, the simulation of system (2.22) is shown in Figure 7. The simulation results indicates the limit cycle near the double homoclinic loop L is unstable, and it will converge to a stable limit cycle surrounding L. The maximal amplitude of limit cycle is about 2.5, which is about 9.17 m in reality. From the steady-state time history of system (2.17) as shown in Figure 8, we can get that the vibration amplitude is about 10 m. The authors conjecture that the reason of the difference between the theoretical and numerical result is induced by high wind velocity.

Conclusion and discussion
The in-plane galloping of iced transmission line is researched in this paper. Considering the geometrical and aerodynamical nonlinearities, we get an one-degree-offreedom model. After rescaling, we get a near-Hamiltonian system, on which we use Melnikov function expansions to study the limit cycles. We find that when the wind velocity is near U = 3.787861482, there are at least 3 limit cycles in this model, 1 limit cycle is near one of the center, 1 limit cycle is between homoclinic loop L 1 and center (0, 0), and last limit cycle surrounds the double homoclinic loop L. By employing numerical simulation, we illustrate the theoretical results and determine the stability of these limit cycles.
The results obtained in this paper show that under some engineering and environmental condition, the iced transmission line system may behave stable limit cycle motion, which is harm to the safety of transmission system. Since the motivation of galloping study is to find the mechanism and to make some suggestions on anti-galloping, in future work, we may select other physical parameters as bifurcation parameters to study the limit cycles and to get the optimal physical parameter combinations for avoiding galloping.
The main contribution of this paper is applying Han's method to mechanics, especially to determine the oscillation of in-plane galloping of iced transmission line. The method we used in this paper is also applied to many other disciplines, such as global bifurcations for a generalized codimension-4 Duffing-van der pol equation [12], periodic waves in a population model with density-dependent migrations and Allee Effect [20], and so on.