A Predictive Model for Dry-Growth Icing on Composite Insulators under Natural Conditions

: Icing can adversely inﬂuence electric power system security. Two main issues are caused by icing: the overload of transmission lines, and the reduction in the insulation ability of the insulators. Most previous research has focused on the ﬂashover characteristics of ice-covered insulators, but research on the icing process of the insulator is seriously lacking. Considering the effect of icing shape, the outer airﬂow ﬁeld of an insulator was calculated and the local collision efﬁciencies of water droplets ( β 1 ) were investigated according to the Lagrange algorithm. The simulation showed that the values of β 1 on the insulator edge and rod are much higher than on the insulator surface, and both were signiﬁcantly inﬂuenced by the wind speed and median volume diameter (MVD) of the water droplets. Based on thermal balance equations, a dynamic dry-growth icing model was established. Using the natural icing conditions of Xuefeng Mountain (China) as an example, validation experiments were conducted on a composite insulator and the climate parameters measured by multi-cylinders were used to model the icing shape and mass. The results indicate that high wind speed and low temperature increase icing rate; the icing was mainly concentrated on the windward side and the greatest horizontal thickness was generally on the insulator edge. The dry-growth model had an average error lower than 25% for icing thickness and an average error lower than 20% for icing mass, which were affected by icing roughness.


Introduction
Icing is a normal phenomenon in cold regions, posing a serious problem for electrical systems. Two main issues are caused by icing: the overload of transmission lines, and the reduction in the insulation property of insulators [1][2][3][4]. To better understand the icing law and provide an appropriate guide for anti-icing and de-icing work, many scholars have researched this issue, resulting in considerable achievements.
Based on heat balance equations, Maeno and Makkonen [5] studied the process of wet-growth icing, and the results showed that the growth rate of an icicle increased with decreasing temperature and increasing wind speed. They also stated that the growth rate of an icicle was affected by the water supply, which could be determined using the cylindrical icing model. Szilder [6] proposed an approach using a three-dimensional (3D) random walk model to simulate the growth of an icicle, enabling the modelling of the ice formation due to the gravitational flow of liquid on objects with complex geometry. Similar theories were adopted to simulate ice density and roughness [7]. Considering the movement of supercooled water droplets in the air, Makkonen [8] proposed an ice growth model according to collision efficiency, sticking efficiency, and accretion efficiency on transmission lines. Makkonen [8] proposed that the growth rate of ice on transmission lines is related to environmental parameters including temperature, wind speed, liquid water content of air, and median volume diameter (MVD) of the water droplet distribution. The accretion efficiency of icing can be obtained by solving heat balance equations, whereas the collision efficiency of droplets can be calculated based on Finstad's fitting formula of transmission lines [9]. However, the introduced model [9] did not consider the effects caused by changing airflow fields and therefore could not accurately obtain the collision efficiency of a changing icing shape.
Fu et al. [10] computed time-dependent airflow and water droplet trajectory of overhead lines using boundary layer theory, which considered the effects caused by icing shape. Subsequently, the transmission line icing model was further improved using different methods like the finite element method [11,12].
However, different from a transmission line, an insulator has a more intricate outer shape, complicating the building of an insulator icing growth model. Most of the previous research has focused on the flashover properties of ice-covered insulators [13][14][15]. Few studies have reported the icing mechanism of insulators. Basically, two main insulator icing types exist: wet-growth icing and dry-growth icing. Wet-growth often appears with icicles growing, whereas dry-growth icing generally occurs in low-temperature and high wind speed environments, where the water droplets trapped on insulator surface can freeze quickly with no continual water film on the icing surface. In addition, the dry-growth icing shape is simple when no icicle grows. The collision efficiency (β 1 ) of composite insulator was calculated under different wind speeds and different MVDs [16,17]. Based on collision efficiency data, the authors introduced a simple rime growing model that can be used for icing thickness calculation under constant environment conditions, and the related tests were conducted in an artificial climate chamber. Nevertheless, the icing on the windward surface of the insulator was regarded as a uniform distribution in these studies, and the effects of changing icing shape on the airflow field and the icing growth rate were ignored. Notable differences between the tests results and actual icing results under natural conditions were also reported in these studies.
Moreover, Farzaneh and Chisholm [18] indicated that the resistance of the ice layer on the insulator is the key parameter for building a flashover model that can be used to predict the lowest flashover voltage. As most of the models are built based on semi-cylinder icing structure, the resistance can be easily obtained. Otherwise, the research is restricted when the icing shape becomes irregular and cannot be considered a semi-cylinder icing structure.
To explore the dry-growth icing characteristics of an insulator under natural conditions and numerically simulate icing shapes, we built a dynamic dry-growth icing model based on a composite insulator. Considering the influences of icing shape and climate factors, the model calculated the 3D airflow fields of the composite insulator with different icing degrees. Based on the Lagrange algorithm, the local collision efficiency (LCE) β 1 was renewed according to the different air flow field and water droplets trajectories. In considering heat, mass balance, and changes in environmental parameters, the ice thickness growth rate was obtained. The related insulator icing experiments were performed at the Xuefeng Mountain Natural Pollution and Icing Observation Station, located in Hunan Province, China, at an altitude of 1400 m [12].

Dry-Growth Icing Characteristics of Insulator
The insulator icing process can be divided into three parts: the collision process of water droplets on the insulator surface, the capture process of the collided water droplets, and the freezing process of the liquid water. The dry-growth icing shape and the icing growth rate are determined by the aforementioned three processes; otherwise, the three processes are affected by various factors including the environment temperature T a , wind speed U and direction, liquid water content of the air w, and the MVD of the water droplets. Moreover, the structure of the composite insulator shed also plays an important role in the icing process by affecting the icing shape and rate. According to long-term observation of insulator icing on Xuefeng Mountain, the dry-growth icing on the windward surface of insulator has been observed to accrete with multi-layer structures, especially for composite insulators. The ice accretes against the wind direction with different horizontal thickness h as shown in Figure 1a, which can easily and eventually cover the whole windward surface of the insulator and bridge the sheds. Inspired by the multi-layer icing phenomenon, we introduce a layer-divided icing calculation method, which considers the LCE, and the icing shapes of different layers are calculated according to its own LCE. According to long-term observation of insulator icing on Xuefeng Mountain, the dry-growth icing on the windward surface of insulator has been observed to accrete with multi-layer structures, especially for composite insulators. The ice accretes against the wind direction with different horizontal thickness h as shown in Figure 1a, which can easily and eventually cover the whole windward surface of the insulator and bridge the sheds. Inspired by the multi-layer icing phenomenon, we introduce a layer-divided icing calculation method, which considers the LCE, and the icing shapes of different layers are calculated according to its own LCE. As shown in Figure 1b, in concentrating on three main insulator positions: the edge, surface, and rod, the structure of 220 kV composite insulator type A was divided into a multi-layer structure (seven layers in Figure 1b for type A). The details of the parameters are listed in Table 1, where R1/R2 represents the radii of the sheds. During the icing process, the outer shape of the insulator changes with increasing icing mass. Hence, the external airflow field changes over time. To accurately simulate the icing shape, the outer airflow field should be renewed often and the collision efficiency β1 should also be recalculated.

External Flow Field Calculation
The basic structure of the insulator was built using AutoCAD, which then was imported into the ANSYS ICEM software for building the grid. With the help of the fluid simulation software FLUENT, the air flow field and the water droplet trajectories were obtained based on the finite element method. During the flow field calculation, the 3D Reynolds-averaged Navier-stokes equation (RANS) was used as the governing equation of the external incompressible viscous turbulent flow field. The standard k-ε turbulent model was used to close the RANS equations. When calculating the trajectories of water droplets, we assumed that the water droplets in the air were standard spheres, evenly distributed, and not be divided or conglomerated. Neglecting the Saffman lift, additional mass force and differential pressure force, the force analysis of water droplets only considered the viscous drag and gravity. The equation of motion for droplets can be expressed as: As shown in Figure 1b, in concentrating on three main insulator positions: the edge, surface, and rod, the structure of 220 kV composite insulator type A was divided into a multi-layer structure (seven layers in Figure 1b for type A). The details of the parameters are listed in Table 1, where R 1 /R 2 represents the radii of the sheds. According to long-term observation of insulator icing on Xuefeng Mountain, the dry-growth icing on the windward surface of insulator has been observed to accrete with multi-layer structures, especially for composite insulators. The ice accretes against the wind direction with different horizontal thickness h as shown in Figure 1a, which can easily and eventually cover the whole windward surface of the insulator and bridge the sheds. Inspired by the multi-layer icing phenomenon, we introduce a layer-divided icing calculation method, which considers the LCE, and the icing shapes of different layers are calculated according to its own LCE. As shown in Figure 1b, in concentrating on three main insulator positions: the edge, surface, and rod, the structure of 220 kV composite insulator type A was divided into a multi-layer structure (seven layers in Figure 1b for type A). The details of the parameters are listed in Table 1, where R1/R2 represents the radii of the sheds. During the icing process, the outer shape of the insulator changes with increasing icing mass. Hence, the external airflow field changes over time. To accurately simulate the icing shape, the outer airflow field should be renewed often and the collision efficiency β1 should also be recalculated.

External Flow Field Calculation
The basic structure of the insulator was built using AutoCAD, which then was imported into the ANSYS ICEM software for building the grid. With the help of the fluid simulation software FLUENT, the air flow field and the water droplet trajectories were obtained based on the finite element method. During the flow field calculation, the 3D Reynolds-averaged Navier-stokes equation (RANS) was used as the governing equation of the external incompressible viscous turbulent flow field. The standard k-ε turbulent model was used to close the RANS equations. When calculating the trajectories of water droplets, we assumed that the water droplets in the air were standard spheres, evenly distributed, and not be divided or conglomerated. Neglecting the Saffman lift, additional mass force and differential pressure force, the force analysis of water droplets only considered the viscous drag and gravity. The equation of motion for droplets can be expressed as: During the icing process, the outer shape of the insulator changes with increasing icing mass. Hence, the external airflow field changes over time. To accurately simulate the icing shape, the outer airflow field should be renewed often and the collision efficiency β 1 should also be recalculated.

External Flow Field Calculation
The basic structure of the insulator was built using AutoCAD, which then was imported into the ANSYS ICEM software for building the grid. With the help of the fluid simulation software FLUENT, the air flow field and the water droplet trajectories were obtained based on the finite element method. During the flow field calculation, the 3D Reynolds-averaged Navier-stokes equation (RANS) was used as the governing equation of the external incompressible viscous turbulent flow field. The standard k-ε turbulent model was used to close the RANS equations. When calculating the trajectories of water droplets, we assumed that the water droplets in the air were standard spheres, evenly distributed, and not be divided or conglomerated. Neglecting the Saffman lift, additional mass force and differential pressure force, the force analysis of water droplets only considered the viscous drag and gravity. The equation of motion for droplets can be expressed as: where the v a and v d are speeds of the airflow and water droplet, respectively; τ is the relaxation time of droplets and τ = (ρ d D w 2 )/(18µ); ρ d and D w are the density and the average diameter of droplets, respectively; f is function of viscous drag and f = C d Re w /24; C d is the drag coefficient; Re w is the Reynolds number of water droplets and Re w = ρ a × MVD (v a − v d )/µ; µ is the dynamic viscosity of air; ρ a is the density of air; and g is the gravitational acceleration [19,20]. As displayed in Figure 2a, the stream lines of water droplets have different bending degrees at the three different positions, which leads to different LCE. The LCE value directly determines the amount of water for icing, so a precise LCE is key to the real-time modeling of icing. With confirmed MVD and wind speed (U), the airflow field and droplet trajectories can be obtained.
where the va and vd are speeds of the airflow and water droplet, respectively; τ is the relaxation time of droplets and τ = (ρd Dw 2 )/(18μ); ρd and Dw are the density and the average diameter of droplets, respectively; f is function of viscous drag and f = CdRew/24; Cd is the drag coefficient; Rew is the Reynolds number of water droplets and Rew = ρa × MVD (va − vd)/μ; μ is the dynamic viscosity of air; ρa is the density of air; and g is the gravitational acceleration [19,20]. As displayed in Figure 2a, the stream lines of water droplets have different bending degrees at the three different positions, which leads to different LCE. The LCE value directly determines the amount of water for icing, so a precise LCE is key to the real-time modeling of icing. With confirmed MVD and wind speed (U), the airflow field and droplet trajectories can be obtained. In this paper, the speeds and relative positions of water droplets were used for calculating the LCE β1 according to Equation (2).
where A, B, C, and D represent the original locations of the four water droplets, A1, B1, C1, and D1, respectively; VA, VB, VC, and VD represent the collision locations and collision speeds of these four droplets, respectively, and S1 and S2 represent the areas of rectangle ABCD and A1B1C1D1, respectively. The ratio of S1 and S2 represents the dispersion degree of the water droplets: the smaller the S1/S2, the more the droplets disperse and the less the water droplets are trapped on one certain area around point P. In addition, the ratio of (VA + VB + VC + VD)/4U represents the collecting rate of droplets, which considers the change in velocity, so the collision efficiency β1 is more accurate. Also, the smaller the areas of S1 and S2, the more precise the calculation of β1 of point P. However, some mistakes may have occurred when the four points are not coplanar, as the algorithm of Equation (2) cannot be obtained, especially when the collided droplets are scattered, so we used three points to implement the algorithm above to improve reliability.
The layered divided structure of insulator type A has seven different diameters. The outermost layer has the maximum diameter of 188 mm (layer-1), whereas the diameter of the rod is only 28 mm (layer-7). The local collision efficiencies of every layer were calculated, and the data were integrated as shown in Figure 3a, where x represents the x-coordinate of the collision position. For every layer, the value of β1 varies with the change in collision position. Specifically, increases from both sides of the layer to the center achieves the maximum value on the center of the windward surface, which is In this paper, the speeds and relative positions of water droplets were used for calculating the LCE β 1 according to Equation (2).
where A, B, C, and D represent the original locations of the four water droplets, A 1 , B 1 , C 1 , and D 1 , respectively; V A , V B , V C , and V D represent the collision locations and collision speeds of these four droplets, respectively, and S 1 and S 2 represent the areas of rectangle ABCD and A 1 B 1 C 1 D 1 , respectively. The ratio of S 1 and S 2 represents the dispersion degree of the water droplets: the smaller the S 1 /S 2 , the more the droplets disperse and the less the water droplets are trapped on one certain area around point P. In addition, the ratio of (V A + V B + V C + V D )/4U represents the collecting rate of droplets, which considers the change in velocity, so the collision efficiency β 1 is more accurate. Also, the smaller the areas of S 1 and S 2 , the more precise the calculation of β 1 of point P. However, some mistakes may have occurred when the four points are not coplanar, as the algorithm of Equation (2) cannot be obtained, especially when the collided droplets are scattered, so we used three points to implement the algorithm above to improve reliability. The layered divided structure of insulator type A has seven different diameters. The outermost layer has the maximum diameter of 188 mm (layer-1), whereas the diameter of the rod is only 28 mm (layer-7). The local collision efficiencies of every layer were calculated, and the data were integrated as shown in Figure 3a, where x represents the x-coordinate of the collision position. For every layer, the value of β 1 varies with the change in collision position. Specifically, increases from both sides of the layer to the center achieves the maximum value on the center of the windward surface, which is the stagnation point (x = 0 mm). Furthermore, different layers have different collision efficiencies. Because of the inclined structure, the values of β 1 on layers-2-6 are much smaller than on the insulator edge and rod (layer-1 and layer-7). As a result, the ice first appears on the insulator windward edge and the rod (layer-1 and layer-7, respectively), which are shown in Figure 3b. Because of the inclined structure, the values of β1 on layers-2-6 are much smaller than on the insulator edge and rod (layer-1 and layer-7). As a result, the ice first appears on the insulator windward edge and the rod (layer-1 and layer-7, respectively), which are shown in Figure 3b. The layers close to the insulator edge are more likely to have a larger β1, but all these values are lower than 0.2 (β1(Layer-2) > β1(layer-6)). In addition, the collision area decreases with decreasing layer diameter, and a fan-shaped icing region is formed.

Influences of Environmental Factors
The wind speed U and the MVD of water droplets are the two main environmental factors that influence the collision characteristics of water droplets. To analyze the icing properties under different environmental conditions, the LCE of insulator type A under different wind speeds and different MVDs were calculated according to the method presented in Section 2.2.
Before encountering the obstacle insulator shed, the water droplets move along with the air flow at the same speed. However, when flowing around the insulator shed, these water droplets can deviate from the air flow stream and bump into the insulator surface, which is mainly determined by their inertia and the drag force of air.
As shown in Figure 4, with increasing wind speed and MVD, the LCE of the insulator edge (layer-1) increases differently: greater wind speed and MVD can both lead to a broader collision span and larger LCE, but an obvious saturation of LCE occurs when the wind speed is greater than 8 m/s, whereas no obvious saturation occurs when MVD increases from 30 to 110 μm. This finding can be explained by comparing the function of wind speed and MVD. Firstly, increasing wind speed can both increase the speed of water droplets and the air. Based on Equation (1), high-speed water droplets are more likely to move along the original orbit and bump into the insulator; however, highspeed airflow can provide a larger drag force that can promote the flowing of water droplets around the insulator. Therefore, the LCE increases minimally when these two factors reach a balance. Correspondingly, a larger MVD only increases the inertia of the water droplets, which is beneficial for water droplets' collision with the insulator. The layers close to the insulator edge are more likely to have a larger β 1 , but all these values are lower than 0.2 (β 1(Layer-2) > β 1(layer-6) ). In addition, the collision area decreases with decreasing layer diameter, and a fan-shaped icing region is formed.

Influences of Environmental Factors
The wind speed U and the MVD of water droplets are the two main environmental factors that influence the collision characteristics of water droplets. To analyze the icing properties under different environmental conditions, the LCE of insulator type A under different wind speeds and different MVDs were calculated according to the method presented in Section 2.2.
Before encountering the obstacle insulator shed, the water droplets move along with the air flow at the same speed. However, when flowing around the insulator shed, these water droplets can deviate from the air flow stream and bump into the insulator surface, which is mainly determined by their inertia and the drag force of air.
As shown in Figure 4, with increasing wind speed and MVD, the LCE of the insulator edge (layer-1) increases differently: greater wind speed and MVD can both lead to a broader collision span and larger LCE, but an obvious saturation of LCE occurs when the wind speed is greater than 8 m/s, whereas no obvious saturation occurs when MVD increases from 30 to 110 µm. This finding can be explained by comparing the function of wind speed and MVD. Firstly, increasing wind speed can both increase the speed of water droplets and the air. Based on Equation (1), high-speed water droplets are more likely to move along the original orbit and bump into the insulator; however, high-speed airflow can provide a larger drag force that can promote the flowing of water droplets around the insulator. Therefore, the LCE increases minimally when these two factors reach a balance. Correspondingly, a larger MVD only increases the inertia of the water droplets, which is beneficial for water droplets' collision with the insulator.

Icing Thickness and Shape on Insulator
Taking one small area of insulator surface as the target, such as element P shown in Figure 2b, and considering the mass balance during droplet collision, the icing growth rate can be expressed as: where M is the mass of accreted ice on the insulator surface; β2 and β3 are the sticking efficiency and accretion efficiency, respectively; w is the liquid water content; and Sp is the element area of the insulator surface. The dry-growth icing occurs under high wind speeds and low environmental temperature, and the water droplets impacting on insulator freeze quickly. Therefore, little bouncing droplets occur and the sticking efficiency is β2 ≈ 1. In addition, the accretion efficiency β3 can be obtained according to the heat balance equation: where qc is the convective thermal loss, qe is the heat loss due to evaporation, qr is the radiation loss, and ql is heat loss in warming the supercooled water to the icing temperature on the insulator surface. In detail, these four parts are calculated as: where Ta is the temperature of the environment, T is the temperature of the ice surface (T = 0 °C), hp is the convection coefficient of insulator surface, e(T) is the saturation vapor pressure over the ice surface at T, and the other parameters are shown in Table 2 [5].

Icing Thickness and Shape on Insulator
Taking one small area of insulator surface as the target, such as element P shown in Figure 2b, and considering the mass balance during droplet collision, the icing growth rate can be expressed as: where M is the mass of accreted ice on the insulator surface; β 2 and β 3 are the sticking efficiency and accretion efficiency, respectively; w is the liquid water content; and S p is the element area of the insulator surface. The dry-growth icing occurs under high wind speeds and low environmental temperature, and the water droplets impacting on insulator freeze quickly. Therefore, little bouncing droplets occur and the sticking efficiency is β 2 ≈ 1. In addition, the accretion efficiency β 3 can be obtained according to the heat balance equation: where q c is the convective thermal loss, q e is the heat loss due to evaporation, q r is the radiation loss, and q l is heat loss in warming the supercooled water to the icing temperature on the insulator surface. In detail, these four parts are calculated as: where T a is the temperature of the environment, T is the temperature of the ice surface (T = 0 • C), h p is the convection coefficient of insulator surface, e(T) is the saturation vapor pressure over the ice surface at T, and the other parameters are shown in Table 2 [5]. For the two parts on the equation's right side, q f is the latent heat released during freezing and q v represents the air flow heating, which is a small part that can be ignored. Therefore, the right side of Equation (4) can be expressed as: Substituting Equations (5)- (7) into Equation (4), β 3 can be expressed as: As shown in Figure 5, during wet-growth icing, a continuous water film is present on insulator surface, from where the unfrozen water flows out the insulator edge and provides the water supply for icicle growth. Due to the water film, the growing ice is vertical to the insulator surface as icing type 1. Otherwise, if the water film is thin and the water flow is not stable, the icicle cannot grow on the insulator edge, as with icing type 2, during which the accretion efficiency β 3 is also smaller than one and horizontal icing thickness accumulates along the insulator surface. Moreover, when water droplets colliding with the insulator freeze quickly into ice without water film forming and flowing, the icing can turn into dry-growth and accretes against the wind direction with different horizontal thickness h, like icing type 3. For the two parts on the equation's right side, qf is the latent heat released during freezing and qv represents the air flow heating, which is a small part that can be ignored. Therefore, the right side of Equation (4) can be expressed as: Substituting Equations (5)- (7) into Equation (4), β3 can be expressed as: As shown in Figure 5, during wet-growth icing, a continuous water film is present on insulator surface, from where the unfrozen water flows out the insulator edge and provides the water supply for icicle growth. Due to the water film, the growing ice is vertical to the insulator surface as icing type 1. Otherwise, if the water film is thin and the water flow is not stable, the icicle cannot grow on the insulator edge, as with icing type 2, during which the accretion efficiency β3 is also smaller than one and horizontal icing thickness accumulates along the insulator surface. Moreover, when water droplets colliding with the insulator freeze quickly into ice without water film forming and flowing, the icing can turn into dry-growth and accretes against the wind direction with different horizontal thickness h, like icing type 3. By this time, the accretion efficiency β3 ≈ 1, and M(1 − β3) mainly represents the evaporation loss. Therefore, type 2 and type 3 can both be considered as dry-growth icing without icicles. Based on mass balance, the icing thickness h in various positions can be calculated as: The ice density ρice (kg/m 3 ) can be calculated using Equation (10), which is an empirical formula determined by wind speed U (m/s), MVD (m), freezing temperature Tf (K) and environment temperature Ta (K) [21]. The constant in the denominator is an empirical constant based on the experiments.

Four Environment Parameters
The rotating multi-cylinder (RMC) method was developed by Langmuir and Blodgett [22], and which has been used in other studies [23,24]; however, the icing mass data could not be continuously By this time, the accretion efficiency β 3 ≈ 1, and M(1 − β 3 ) mainly represents the evaporation loss. Therefore, type 2 and type 3 can both be considered as dry-growth icing without icicles. Based on mass balance, the icing thickness h in various positions can be calculated as: The ice density ρ ice (kg/m 3 ) can be calculated using Equation (10), which is an empirical formula determined by wind speed U (m/s), MVD (m), freezing temperature T f (K) and environment temperature T a (K) [21]. The constant in the denominator is an empirical constant based on the experiments.

Four Environment Parameters
The rotating multi-cylinder (RMC) method was developed by Langmuir and Blodgett [22], and which has been used in other studies [23,24]; however, the icing mass data could not be continuously measured and the surveying work basically depended on manual measurements. To automatically obtain the environment parameters under natural conditions and in real-time, the RMC method was improved in this paper. As shown in Figure 6, five rotating cylinders with different diameters (d = 10, 15, 20, 25, and 30 mm, L = 230 mm) were used for monitoring the changing icing mass. In the experiments, the micro tension sensor was designed to connect with the motor and the icing cylinder, which can transmit the ice mass signal on every cylinder in real time. Also, the five cylinders were installed with enough space to ensure they did not affect other's air flow field during the icing period.
Energies 2017, 10, x 8 of 16 measured and the surveying work basically depended on manual measurements. To automatically obtain the environment parameters under natural conditions and in real-time, the RMC method was improved in this paper. As shown in Figure 6, five rotating cylinders with different diameters (d = 10, 15, 20, 25, and 30 mm, L = 230 mm) were used for monitoring the changing icing mass. In the experiments, the micro tension sensor was designed to connect with the motor and the icing cylinder, which can transmit the ice mass signal on every cylinder in real time. Also, the five cylinders were installed with enough space to ensure they did not affect other's air flow field during the icing period. Different than the insulator, the cylinder has a much simpler geometry and its collision efficiency β1 can be obtained according to Finstad's fitting formula of a conductor [9]. Moreover, the accretion efficiency β3 can be obtained according to Equation (8). Because β1 is determined by MVD, wind speed U, and diameter of the cylinder, β3 is determined by Ta, w, and U, based on Equation (3). The icing rates of cylinders can be expressed as a function of the four environmental parameters and the diameter of cylinder, which is shown as: where the dMi/dt (i = 1-4) is the growth rate of ice mass, di (i = 1-4) is the diameter of cylinder, and f is the function built based on Equation (3), where β1 and β3 will be substituted by environmental parameters: Using four cylinders as calculation samples and the remaining cylinder as the verification sample, the four climate parameters were back calculated by Equation (12), where Fi (i = 1-4) is the inverse function of f and is indirectly obtained using the differential evolution algorithm, which is one of the mature numerical calculate methods [25]. To ensure the accuracy of the climate parameters, the verification error of cylinder sample 5 had to be smaller than 10%, and the airflow field was updated when the icing significantly changed the outer shape of insulator. The detailed modelling flow chart is shown as Figure 7.
The term N is the maximum iterative steps, ∆t1 is the update time step of the icing structure, and ∆t2 is the update time step for calculating the climate parameters. Basically, ∆t2 is set to one minute, and ∆t1 is greater than ∆t2 to save computational time. The smaller the ∆t1, the more accurate the icing shape and thickness obtained. During the dry-growth icing, the horizontal ice thicknesses of different Different than the insulator, the cylinder has a much simpler geometry and its collision efficiency β 1 can be obtained according to Finstad's fitting formula of a conductor [9]. Moreover, the accretion efficiency β 3 can be obtained according to Equation (8). Because β 1 is determined by MVD, wind speed U, and diameter of the cylinder, β 3 is determined by T a , w, and U, based on Equation (3). The icing rates of cylinders can be expressed as a function of the four environmental parameters and the diameter of cylinder, which is shown as: where the dM i /dt (i = 1-4) is the growth rate of ice mass, d i (i = 1-4) is the diameter of cylinder, and f is the function built based on Equation (3), where β 1 and β 3 will be substituted by environmental parameters: Using four cylinders as calculation samples and the remaining cylinder as the verification sample, the four climate parameters were back calculated by Equation (12), where F i (i = 1-4) is the inverse function of f and is indirectly obtained using the differential evolution algorithm, which is one of the mature numerical calculate methods [25]. To ensure the accuracy of the climate parameters, the verification error of cylinder sample 5 had to be smaller than 10%, and the airflow field was updated when the icing significantly changed the outer shape of insulator. The detailed modelling flow chart is shown as Figure 7.
The term N is the maximum iterative steps, ∆t 1 is the update time step of the icing structure, and ∆t 2 is the update time step for calculating the climate parameters. Basically, ∆t 2 is set to one minute, and ∆t 1 is greater than ∆t 2 to save computational time. The smaller the ∆t 1 , the more accurate the icing shape and thickness obtained. During the dry-growth icing, the horizontal ice thicknesses of Energies 2018, 11, 1339 9 of 16 different layers are different from each other. The icing on each layer will be calculated separately, and the combination of the icing shape of the insulator means the integration of icing shape for every layer.
Energies 2017, 10, x 9 of 16 layers are different from each other. The icing on each layer will be calculated separately, and the combination of the icing shape of the insulator means the integration of icing shape for every layer.

Dry-Growth Icing Tests and Simulation
To verify the validity of the dry-growth icing model to ensure the icing shape and mass under changing environment conditions could be obtained, related icing tests were conducted at the Xuefeng Mountain Natural Pollution and Icing Observation Station, located in the typical icing region in Southern China. This area has an average annual precipitation of 1500 mm, maximum wind speed of 35 m/s, minimum air temperature of −15 °C, and a maximum icing thickness of 500 mm [26,27].
As shown in Figure 8, the five cylinders had uniform icing shapes and the icing mass increased over time. Notably, the icing rates were relatively small at 60 min or less and increased dramatically after 120 min. In addition, the bigger cylinder had heavier ice accretion because of its larger surface area. According to the flow chart in Figure 7, the icing mass of the four cylinders (No. 1-4) were imported into the calculation model, and the data of cylinder No. 5 were used for verification. The four climate parameters were obtained as shown in Figure 9. Over a four-hour icing period, the wind speed increased from 3.5 to 8.5 m/s, the temperature dropped over time, which corresponded to the increasing icing rates. The liquid water content of the air varied between 0.6 and 1 g/m 3 , and the MVD changed between 38 and 55 μm.

Dry-Growth Icing Tests and Simulation
To verify the validity of the dry-growth icing model to ensure the icing shape and mass under changing environment conditions could be obtained, related icing tests were conducted at the Xuefeng Mountain Natural Pollution and Icing Observation Station, located in the typical icing region in Southern China. This area has an average annual precipitation of 1500 mm, maximum wind speed of 35 m/s, minimum air temperature of −15 • C, and a maximum icing thickness of 500 mm [26,27].
As shown in Figure 8, the five cylinders had uniform icing shapes and the icing mass increased over time. Notably, the icing rates were relatively small at 60 min or less and increased dramatically after 120 min. In addition, the bigger cylinder had heavier ice accretion because of its larger surface area. According to the flow chart in Figure 7, the icing mass of the four cylinders (No. 1-4) were imported into the calculation model, and the data of cylinder No. 5 were used for verification. The four climate parameters were obtained as shown in Figure 9. Over a four-hour icing period, the wind speed increased from 3.5 to 8.5 m/s, the temperature dropped over time, which corresponded to the increasing icing rates. The liquid water content of the air varied between 0.6 and 1 g/m 3 , and the MVD changed between 38 and 55 µm.   As described in Section 3, the structure parameters of insulator type A were substituted into the dry-growth model. Based on the results obtained in Figure 9, considering the low icing rate and avoiding a massive calculation, the time step ∆t 1 was set to 60 min. During this period, the average environment parameters (Table 3) were used as the inputs of the model, and the grid was redrawn and the air and water flow fields were recalculated according to the changing icing shape. In the icing tests, shown in Figure 10, a tension senor with a proper measurement range of 40 kg, was connected to the insulator to survey the changing ice mass. The sampling rate of the tensor was set as once per second. The greatest ice thickness of the three positions (insulator surface, rod, and edge) was measured at 30 min intervals. This process improved the measurement precision by averaging the test values of the repeated measurements.
Energies 2017, 10, x 11 of 16 As described in Section 3, the structure parameters of insulator type A were substituted into the dry-growth model. Based on the results obtained in Figure 9, considering the low icing rate and avoiding a massive calculation, the time step ∆t1 was set to 60 min. During this period, the average environment parameters (Table 3) were used as the inputs of the model, and the grid was redrawn and the air and water flow fields were recalculated according to the changing icing shape. In the icing tests, shown in Figure 10, a tension senor with a proper measurement range of 40 kg, was connected to the insulator to survey the changing ice mass. The sampling rate of the tensor was set as once per second. The greatest ice thickness of the three positions (insulator surface, rod, and edge) was measured at 30 min intervals. This process improved the measurement precision by averaging the test values of the repeated measurements.  Figure 11 shows the comparative icing results of the four-hour test and simulation. Overall, of changeless wind direction during the test, the dry-growth icing had typical regional characteristics, with the icing mainly concentrated on the windward surface of insulator. However, the icing manifested with different characteristics at different stages.

Results and Discussion
(a)  Figure 11 shows the comparative icing results of the four-hour test and simulation. Overall, of changeless wind direction during the test, the dry-growth icing had typical regional characteristics, with the icing mainly concentrated on the windward surface of insulator. However, the icing manifested with different characteristics at different stages.

Results and Discussion
Energies 2017, 10, x 11 of 16 As described in Section 3, the structure parameters of insulator type A were substituted into the dry-growth model. Based on the results obtained in Figure 9, considering the low icing rate and avoiding a massive calculation, the time step ∆t1 was set to 60 min. During this period, the average environment parameters (Table 3) were used as the inputs of the model, and the grid was redrawn and the air and water flow fields were recalculated according to the changing icing shape. In the icing tests, shown in Figure 10, a tension senor with a proper measurement range of 40 kg, was connected to the insulator to survey the changing ice mass. The sampling rate of the tensor was set as once per second. The greatest ice thickness of the three positions (insulator surface, rod, and edge) was measured at 30 min intervals. This process improved the measurement precision by averaging the test values of the repeated measurements.  Figure 11 shows the comparative icing results of the four-hour test and simulation. Overall, of changeless wind direction during the test, the dry-growth icing had typical regional characteristics, with the icing mainly concentrated on the windward surface of insulator. However, the icing manifested with different characteristics at different stages.

Results and Discussion
(a) In the early icing period (0 ≤ t ≤ 60 min), icing first appeared on the insulator edge and rod, while only a little squamous ice accreted on the insulator windward surface. With the increase in icing, the layer-divided characteristic of icing became obvious; the icing outline changed to a stair-step shape, roughening the ice surface (60 ≤ t ≤ 180 min). In the later icing period (t > 180 min), the ice on the insulator rod and surface jutted out, the gap between sheds tended to be filled, and the insulator shed was bridged by continual rough ice.
As shown in Figure 12a, the changing ice shape of the vertical section along the insulator axis indicates that the icing rates at different positions of the insulator are different from each other. Specifically, the icing rates at the insulator edge and insulator rod were much greater than that at insulator surface, which corresponds to the field observation and the collision efficiency simulation results in Section 2. Two two phenomena should be noted, as shown in Figure 12b. Firstly, the icing rates were different during different icing periods; the icing rates at three positions on the insulator were all lower than 0.15 mm/min before 60 min, which dramatically increased after 120 min. A maximum icing thickness of 66 mm was recorded on the insulator rod after four hours, which can be explained by the changing ice shape and flow field. With increasing icing, the ice shape tends to be sharp and In the early icing period (0 ≤ t ≤ 60 min), icing first appeared on the insulator edge and rod, while only a little squamous ice accreted on the insulator windward surface. With the increase in icing, the layer-divided characteristic of icing became obvious; the icing outline changed to a stair-step shape, roughening the ice surface (60 ≤ t ≤ 180 min). In the later icing period (t > 180 min), the ice on the insulator rod and surface jutted out, the gap between sheds tended to be filled, and the insulator shed was bridged by continual rough ice.
As shown in Figure 12a, the changing ice shape of the vertical section along the insulator axis indicates that the icing rates at different positions of the insulator are different from each other. Specifically, the icing rates at the insulator edge and insulator rod were much greater than that at insulator surface, which corresponds to the field observation and the collision efficiency simulation results in Section 2.
(b) Figure 11. Icing shapes during a four-hour test and simulation. Icing was removed from the top layer for comparison. (a) The changing ice shapes on insulator sample (from left to right: t = 60, 120, 180, and 240 min); (b) Simulated ice shapes based on climate parameters and mesh reconstruction (from left to right: t = 0, 60, 120, 180, and 240 min).
In the early icing period (0 ≤ t ≤ 60 min), icing first appeared on the insulator edge and rod, while only a little squamous ice accreted on the insulator windward surface. With the increase in icing, the layer-divided characteristic of icing became obvious; the icing outline changed to a stair-step shape, roughening the ice surface (60 ≤ t ≤ 180 min). In the later icing period (t > 180 min), the ice on the insulator rod and surface jutted out, the gap between sheds tended to be filled, and the insulator shed was bridged by continual rough ice.
As shown in Figure 12a, the changing ice shape of the vertical section along the insulator axis indicates that the icing rates at different positions of the insulator are different from each other. Specifically, the icing rates at the insulator edge and insulator rod were much greater than that at insulator surface, which corresponds to the field observation and the collision efficiency simulation results in Section 2. Two two phenomena should be noted, as shown in Figure 12b. Firstly, the icing rates were different during different icing periods; the icing rates at three positions on the insulator were all lower than 0.15 mm/min before 60 min, which dramatically increased after 120 min. A maximum icing thickness of 66 mm was recorded on the insulator rod after four hours, which can be explained by the changing ice shape and flow field. With increasing icing, the ice shape tends to be sharp and Two two phenomena should be noted, as shown in Figure 12b. Firstly, the icing rates were different during different icing periods; the icing rates at three positions on the insulator were all lower than 0.15 mm/min before 60 min, which dramatically increased after 120 min. A maximum icing thickness of 66 mm was recorded on the insulator rod after four hours, which can be explained by the changing ice shape and flow field. With increasing icing, the ice shape tends to be sharp and rough, especially for the stagnation line of the insulator's windward edge and rod. The sharp and rough ice shape tended toward greater collision efficiency β 1 , which can easily capture the water droplets. As a consequence, the icing rate at a later icing stage will be greater than during the initial icing stage. Also, according to the changing climate parameters obtained in Figure 9, after 120 min, with a stable water content of 0.8 g/m 3 and a MVD close to 50 µm, the high wind speed of around eight m/s and the low temperature were both beneficial to icing. Secondly, the ice thickness simulation results were in good agreement with the test results before 60 min, with an average error less than 20%, which increased to 25% after 120 min. Moreover, assuming the insulator sheds with the same diameter had the same icing rate, by using finite element method, the icing mass could be calculated, which was compared with real-time monitored test results displayed in Figure 13. Similarly, the icing mass simulation results were all less than obtained in the test with an increasing average error of 20%, and the quality of icing on the whole insulator was 14.8 kg after four hours as shown in Table 4, which was even greater than mass of insulator itself (10.67 kg).
Energies 2017, 10, x 13 of 16 rough, especially for the stagnation line of the insulator's windward edge and rod. The sharp and rough ice shape tended toward greater collision efficiency β1, which can easily capture the water droplets. As a consequence, the icing rate at a later icing stage will be greater than during the initial icing stage. Also, according to the changing climate parameters obtained in Figure 9, after 120 min, with a stable water content of 0.8 g/m 3 and a MVD close to 50 μm, the high wind speed of around eight m/s and the low temperature were both beneficial to icing. Secondly, the ice thickness simulation results were in good agreement with the test results before 60 min, with an average error less than 20%, which increased to 25% after 120 min. Moreover, assuming the insulator sheds with the same diameter had the same icing rate, by using finite element method, the icing mass could be calculated, which was compared with real-time monitored test results displayed in Figure 13. Similarly, the icing mass simulation results were all less than obtained in the test with an increasing average error of 20%, and the quality of icing on the whole insulator was 14.8 kg after four hours as shown in Table 4, which was even greater than mass of insulator itself (10.67 kg).  These findings can be explained by the following two reasons. The simulation used a constant time step of one hour, which reduced the calculation but otherwise led to a larger deviation when the icing rate was large and the outer shape of insulator changed quickly. The greatly increased roughness can also increase simulation error when more sampling points are needed to build an accurate model. The mesh of the icing insulator should be extremely refined at the sharp icing edge, which, however, is hard to achieve with a reasonable computational cost. As displayed in Figure 14, in the initial stage of icing, smooth icing accumulated on the insulator edge and rod, which had a thickness that increased from two sides to the windward centre. The simulation was able to accurately model the icing shape. In the later stage of icing, the icing surface became extremely rough and sharp, and the simulation results deviated more.   These findings can be explained by the following two reasons. The simulation used a constant time step of one hour, which reduced the calculation but otherwise led to a larger deviation when the icing rate was large and the outer shape of insulator changed quickly. The greatly increased roughness can also increase simulation error when more sampling points are needed to build an accurate model. The mesh of the icing insulator should be extremely refined at the sharp icing edge, which, however, is hard to achieve with a reasonable computational cost. As displayed in Figure 14, in the initial stage of icing, smooth icing accumulated on the insulator edge and rod, which had a thickness that increased from two sides to the windward centre. The simulation was able to accurately model the icing shape. In the later stage of icing, the icing surface became extremely rough and sharp, and the simulation results deviated more.

Conclusions
The LCE β1 varies with the change in the collision position on the insulator. Specifically, increases from both sides of the shed toward the windward center creates the maximum β1 on the edge and rod. Higher wind speed and greater MVD both lead to a larger collision efficiency β1, but an obvious saturation of collision efficiency β1 occurs when U increases with MVD fixed, and the obvious saturation disappears when MVD increases with U fixed.
During the natural icing process, the stability of the flowing water film on the insulator surface was the key parameter determining icing type: wet-growth with icicle, wet-growth without icicle, or dry-growth. The latter two types have time-varying horizontal icing thicknesses that were simulated according to the change in LCE β1. Considering the redrawing of the grid and the change climate parameters, a dynamic dry-growth (including wet-growth without icicle) icing model was established based on heat and mass balance. The simulation of icing thickness and mass showed that high wind speeds and low temperatures promote the icing rate.
Verification tests were performed under natural conditions at Xuefeng Mountain. Taking advantage of the decentralized installation and automatic measurement of icing mass data, the RMC method was improved to measure climate parameters. The icing tests results proved that the drygrowth model has an average simulation error less than 25% for icing thickness and an average error less than 20% for icing mass, and the simulation accuracy decreased due to icing roughness.
Based on the distribution of LCE β1 on the insulator surface, an increase in the inclination angle of the insulator surface is beneficial for reducing ice accretion, and moderate increases in shed diameter and decreases in insulator rod diameter are also effective in slowing the ice-bridging process.
Author Contributions: X.J. and X.H. proposed the study concepts. X.H. and C.B. carried out the related experiments. X.H. built the predictive icing model, analyzed the test results, and wrote this paper. Z.Y. checked the language and revised the contents of this manuscript.

Conclusions
The LCE β 1 varies with the change in the collision position on the insulator. Specifically, increases from both sides of the shed toward the windward center creates the maximum β 1 on the edge and rod. Higher wind speed and greater MVD both lead to a larger collision efficiency β 1 , but an obvious saturation of collision efficiency β 1 occurs when U increases with MVD fixed, and the obvious saturation disappears when MVD increases with U fixed.
During the natural icing process, the stability of the flowing water film on the insulator surface was the key parameter determining icing type: wet-growth with icicle, wet-growth without icicle, or dry-growth. The latter two types have time-varying horizontal icing thicknesses that were simulated according to the change in LCE β 1 . Considering the redrawing of the grid and the change climate parameters, a dynamic dry-growth (including wet-growth without icicle) icing model was established based on heat and mass balance. The simulation of icing thickness and mass showed that high wind speeds and low temperatures promote the icing rate.
Verification tests were performed under natural conditions at Xuefeng Mountain. Taking advantage of the decentralized installation and automatic measurement of icing mass data, the RMC method was improved to measure climate parameters. The icing tests results proved that the dry-growth model has an average simulation error less than 25% for icing thickness and an average error less than 20% for icing mass, and the simulation accuracy decreased due to icing roughness.
Based on the distribution of LCE β 1 on the insulator surface, an increase in the inclination angle of the insulator surface is beneficial for reducing ice accretion, and moderate increases in shed diameter and decreases in insulator rod diameter are also effective in slowing the ice-bridging process.