Numerical Modeling of Ice Accumulation on Three-Dimensional Bridge Cables under Freezing Rain and Natural Wind Conditions

: In order to accurately predict the ice accumulation on bridge cables under two typical freezing rain conditions, rime and glaze ice, this paper proposes a numerical simulation framework based on the three-dimensional Messinger theory. Two technical challenges of determining the flow direction of unfrozen water and solving three-dimensional Messinger equations are solved in this research. Based on the outflow, mass was calculated according to the three-dimensional Messinger theory, and the flow direction of unfrozen water in each cell was determined by the resultant force of air shear stress and water film gravity. To solve the three-dimensional equations, an iterative method without finding the stagnation line was introduced. The final iced geometries were determined when the inflow mass ratio was satisfied with the converge criteria. Moreover, this modified numerical model was programmed and embedded into computational fluid software. For both two typical freezing rain conditions, the effects of temperature and wind speed on iced geometries were studied. The aerodynamic characteristics and galloping instability of bridge cables with different iced geometries were also investigated. These preliminary aerodynamic simulations can provide the basis for the wind-induced vibration analysis of the whole structure.


Introduction
Bridges are critical for connecting people and goods in the transportation system [1,2]. Freezing rain has a serious impact on daily human life, such as the destruction of infrastructure facilities, the reduction of agriculture outputs, and even the malfunction of aircraft [3,4]. In 2008, a massive freezing rain and snowstorm stuck southern China, causing severe damage to a vast number of buildings, bridges, and power units. Freezing raindrops are a kind of super-cooled water droplet that exhibit strong cohesion when they impact an object surface below 0 °C. The attached ice on a structure surface could change its aerodynamic characteristics and cause unexpected vibrations or damages. The wind resistance stability of cable after icing, especially under the action of strong wind, is important for the study of the potential impact of icing on bridges [5,6].
Cable-supported bridges are efficient types of bridges in mountainous areas [7,8]. However, rime ice and glaze ice climate weathers often occur in the mountains, which has a huge impact on the safety of these structures [9]. Installed structural health monitoring systems [10,11] could monitor the change of temperature to warn stakeholders [12]. The hazards of freezing rain to bridge cables include significantly increasing the self-weight of components and changing structural dynamic characteristics. In addition, icing leads to changes in the designed aerodynamic shape, which also significantly affects the structure's local wind vibration stability, and induces large vibrations such as cable galloping vibration [13,14]. Therefore, the mechanism of ice accumulation and final iced geometries are key issues in research on the aerodynamic performance of the cable system bridge after icing.
Previous research has investigated ice accumulation on power supply cables [15,16]. However, bridge cables cannot rotate, considering the stiffness of a bridge cable is far greater than that of a wireline cable [17]. Thus, these conclusions cannot be applied to the ice accumulation on bridge cables. Demartino et al. [18], Demartino and Ricciardelli [19], and Górski et al. [20] studied the ice accumulation phenomenon and the iced shapes of vertical and inclined cables. It was found that the runback water builds icicles, and a thin ice accumulation is formed. Koss et al. [21] conducted the experiment of ice accumulation on bridge cables and investigated the final iced shape and aerodynamics coefficients under different wind speed, temperature, and inclined angle. Since it is difficult to obtain a parametric model of iced shapes from experiments, numerical simulation plays an important role in the accurate prediction of ice accumulation geometries.
The common numerical model is based on the Messinger model, which has simplified assumptions on water film behavior [22]. The software FENSAP-ICE in Canada [23] considered air shear stress as the main driving force of water film flow. Myers et al. [24,25] developed a water film icing model based on lubrication theory and assumed the water film flow direction was affected by the shear stress, gravity, pressure gradients, and surface tension. Leng et al. [26] investigated that the dynamics of water film as affected by wind shear and gravity in the laboratory. While these studies are focused on in-cloud icing simulation of the airfoil [27], the results cannot be directly applied to freezing rain accumulation on bridge cable due to the significant difference in wind speed and temperature. However, the assumptions about water flow directions provide very important references.
The amount of icing and the shape of icing are two important issues in ice accumulation. Porcu [28] established a two-dimensional and three-dimensional random icing model of the cable wires, and reproduced the process of random ice formation on the circular protrusions. Myers [29] established a mathematical model of icing growth of numerous shapes, and showed that there existed a similar pattern in icing growth among these circular cross-sections. EB Lébatto [30] also conducted research on the impact of freezing rain icing on wire cables, which revealed that the ice accumulations were more easily formed in rime climate weather, and the final iced shapes varied for different wires. However, considering that the rigidity of cables in cable-supported bridges is much greater than that of power wires, and are difficult to be twisted due to their huge selfweight, these conclusions on power transmission wires cannot be simply applied to the cables of cable-supported bridges.
Recently, some researchers have focused on the effect of ice accumulation on bridge hangers [31] and circular cylinders [32,33]. Koss and Lund [34] revealed that iced cables in freezing rain weather have a potential for galloping instability, and the significant changes to aerodynamic cable forces should be paid more attention [35,36]. Xu [37] revealed the aerodynamic characteristics of pipelines under the glaze ice condition. The aerodynamic characteristics of iced bridge hangers were investigated by Guo et al. [38], who described the ice accumulation progress and the final thicknesses of two common types of iced shape. The cables on the cable-supported bridges are huge, and more complicated than circular cylinders on other types of bridges. Therefore, it is necessary to have a more realistic simulation of the ice accumulation of bridge cables to reveal the final iced shapes and their impacts on the bridge's wind-resistant performance.
Aiming at two typical freezing rains, this paper studied the mechanism and final geometries of the ice accumulation on bridge cables. Section 2 describes the numerical simulation framework of ice accumulation. Section 3 introduces the three-dimensional Messinger theory into the numerical models to simulate the ice accumulation progress of bridge cables, which could effectively handle the problem of unfrozen water flow direction. Section 4 contains the numerical setups of rime and glaze ice climate conditions. Section 5 to Section 7 show the simulated iced geometries on bridge cable in different climate conditions, and the aerodynamic forces, critical wind speed and wind attack angle of these iced cables are also studied. Finally, appropriate conclusions are drawn in Section 8.

Numerical Simulation Framework of Ice Accumulation
The numerical simulation of ice accumulation includes the following steps: (1) Building the geometric model of the object; (2) Choosing the appropriate numerical simulation method, including the multiphase flow model and the turbulence model; (3) Embedding the modified Messinger model into AYSYS FLUENT 19.0. After the solution of the flow field is completed, the ice-covered mass on the element surface can be obtained; (4) Based on the ice density and its relationship to temperature, the ice thickness can be obtained, and the iced geometry and surface grid can be reconstructed; (5) Determining whether the calculation has converged. If it is converged, the program will stop the simulation; otherwise, return to step (3).
Based on the conservation equations of the classical Messinger model and the software AYSYS FLUENT 19.0 [39], a three-dimensional simulation model that can simulate the ice accumulation process on bridge cable and analyze multiphase flow field was established. The flow direction of unfrozen water was determined by considering the air shear stress and gravity, and the three-dimensional Messinger equation was solved by the iteration method. Moreover, the whole process was programmed and embedded into the commercial software for ice accumulation analysis. According to this ice accumulation process, the ice thickness can be obtained, and the geometry of the bridge cable reconstructed. Furthermore, the ice-covered geometries of cables at different temperatures and wind speeds under rime and glaze conditions were analyzed individually.
The flow of the whole simulation process is outlined in Figure 1. Since the emphasis of this paper is to explain the specific numerical simulation method of rime-ice accumulation and glaze-ice accumulation on a three-dimensional element, the key steps (4) and (5) are described in detail. The Euler method was employed to achieve the task of calculating the two-phase flow using one set of mesh at the same time. The Euler method has a relatively computational efficiently and can simulate the ice accumulation problem of three-dimensional bridge element with complex geometry [40]. In order to simplify the analysis, the assumptions in the liquid phase were introduced to the Euler method as follows: (a) the droplets were treated as solid spheres; (b) the collision, deformation, and heat transfer inside the liquid phase were not considered; (c) the air turbulent impulse did not affect the movement of the water droplets.
Since the iced geometry of a bridge element continuously changes with increasing ice thickness, dense wall mesh makes it difficult to use the dynamic mesh method, so the k-ε turbulence model with a relatively large near-wall coefficient was selected.

The Governing Equation of Messinger Model
The essence of the Messinger model is a basic framework of ice accumulation simulation, and it assumes that water on the element surface exists in the form of water film. It can also comprehensively describe the relationship of mass conservation and energy conservation within a particular cell during the ice accumulation process, as shown in Figure  2. Since researchers can determine the unknown terms and the calculation methods of those terms, these two sets of conservation equations can be customized for different problems, making the Messinger model extremely applicable. However, for the three-dimensional ice accumulation problem, the unfrozen water flow direction on the element surface is very difficult to determine.
where denotes the rate of water droplets impingement, denotes the rate of water inflow, denotes the rate of water evaporation or sublimation depending on glaze ice case or rime ice case, denotes rate of water flowing out, and denotes the ice accumulation rate.

Energy Conservation Equation
The energy conservation is given by Equation (2): where E so denotes the freezing energy of liquid water, which includes the energy of liquid water itself and the latent heat released during the freezing process; E es denotes the energy required for liquid water to evaporation or sublimation, which also includes the release of latent heat during the phase changing process; E im denotes the energy from impinging droplets; E ou and E in denote the energy change in the control cell due to outflow and inflow, respectively; Q denotes the amount of friction heat between airflow and element surface; and Q denotes the convective heat [41].
For an arbitrarily surface cell A, if the temperature T and wind velocity v are already obtained, then the quantity , , and can be calculated. Therefore, the quantities needed to be solved are , , , and T (the temperature during ice accumulation process). In order to facilitate the solution, freezing fraction f is introduced, defined as the ratio of the mass of water frozen to ice in a control cell to total water mass of this cell.
Using Equation (3) substituting the m so in Equation (1): Thus, the Messinger model can be represented by Equations (2) and (4). The unknowns in these equations are m in , m ou , f and T, where T is the equilibrium temperature of the control cell after ice accumulation. Attention for the implicit relation of f and T, the temperature of the ice water mixture is 0 °C, which is means if 0 < f < 1, T = 1. There are three equations and four unknowns, which means the number of unknowns need to be further reduced for solving the equations.

Solution of the Three-Dimensional Messinger Model
Some scholars believe that a position only absorbs water droplets and never gets inflow from the surrounding positions, and this position is defined as "the stagnation point". As shown in Figure 3, point A is a stagnation point. In order to reduce the number of unknowns and solve Messinger equations, it is necessary to find this point at first and add the additional equation of m in =0; then, the m ou of the stagnation point can be solved and used as m in for adjacent cells. The unknown terms of those cells, once they have been solved, can again be treated as m in for the adjacent cells of those cells. Loop over all cells at the calculation domain and the solution will be obtained. However, for the three-dimensional problem, solving the Messinger equations has two major difficulties: (1) determining the outflow direction of unfrozen water, and (2) getting the position of the stagnation line. The solution schemes proposed in this paper are listed below.

Outflow Direction
It is assumed that the flow direction of the water film only relates to the resultant force F ⃗ of the air shear stress and the self-weight of water film. Assuming that a threedimensional control volume cell on the surface of a bridge cable is C0, and it is surrounded by other four control cells C1~C4, the normal vectors of the interfaces are A ⃗~A ⃗ which point to the direction away from C0, and the water film flow direction is shown in Figure  4. If the mass of unfrozen water in C0 at a certain time is m ou , the inflow mass of each adjacent cell is calculated according to Equation (5), where i is the identifier of surrounding control cell, and K i is defined by Equation (6).

Iterative Messinger Model
Since the stagnation point is difficult to determine in the three-dimensional problem, this paper employed an iterative Messinger model to solve the conservation equations. The flowchart of the iteration process is shown in Figure 5. Firstly, assuming all control cells are stagnation points, after introducing the additional condition of m in = 0, every cell can get an initial solution m ou through Messinger equations. According to the determining scheme of outflow direction proposed in Section 3.2.1, updating the inflow mass of each surrounding cell, the calculation is based on Equation (7) to determine whether it is converged. If it is not converged, then the undated m in can be used to solve the equations and to continue the steps above until it is converged.
In this section, the numerical simulation method of the ice accumulation process under freezing rain conditions was studied. The Messinger model was modified for solving the three-dimensional problems, which effectively considered the unfrozen water flow directions and solved the three-dimensional Messinger equations by the iteration method. The whole process is suitable for the simulation of ice accumulation on bridge cables under freezing rain conditions.

Numerical Simulation Setups
Based on the meteorological data of the Yunnan-Guizhou Plateau in China, the distribution law and duration of freezing rain were obtained, and information such as temperature, rainfall intensity and wind speed during the freezing rain period were extracted. This measured environmental information could provide a solid basis for the realistic numerical simulations of rain and ice accumulation on cables.
Rime and glaze are two common types of ice accumulation for freezing rain conditions. Rime is built up of ice (in the form of small white crystals) when small water droplets freeze upon contact with a cold surface. Glaze ice forms when a relatively large amount of liquid water (a raindrop, or several) freezes while in contact with a surface. Although they are the result of freezing rain, or freezing after rain or heavy fog, these two climate conditions and icing simulations should be individually studied due to their variations in water droplet diameter and ambient temperature.

Rime Climate Condition
The diameter of a rime water droplet is small, about 0.01~0.08mm. Therefore, the water droplets are often suspended in the air with a lower temperature. Ice accretion under rime weather has a milky white opaque body with bubble voids, loose texture, and an undulating surface with no fixed shape. The boundary conditions should be precisely simulated according to these characteristics. The left entrance is set to be the mixed phase's velocity inlet, and the right side is defined as the pressure outlet of the mixed-phase, while the other four sides are all set to be symmetrical. All the boundary settings for rime ice simulation are shown in Figure 6. The structural hexahedral mesh was employed to improve computational efficiency and quality, and the mesh quantity was densified only near the cable surface. The total mesh number was 1.3 million. The minimum grid volume was 4.4 × 10 −6 m 3 , and the maximum grid volume was 9.1 × 10 −4 m 3 . The mesh layout is shown in Figure 7. To ensure the accuracy of simulation, the time step was set to be 0.01 s. Based on this small time step, it was too difficult to simulate the whole process of ice accumulation, which normally takes several hours. Therefore, to speed up the simulation process, the mesh deformation obtained in each step was multiplied by 1000 for the rime ice case, and the deformation of each grid was limited below 0.028 mm per step to ensure the negative volume mesh did not appear. Finally, there was a 1000-time step for simulating the rime ice accumulation, which corresponded to 10,000 s (2.8 h) in the actual climatic condition. The ice accretion mass and thickness were generally logarithmically related with the time symmetry symmetry step. The ice accretion grew rapidly at the beginning, slowing down over time, and the growth was not significant at the set-up step.

Glaze Climate Condition
The meteorological conditions of glaze weather are significantly different from those of rime weather. The water content of the air in glaze weather is about 2.5 × 10 −7 , but the diameter of the water droplets is larger, generally in the range of 0.1~0.5 mm. Ice accretion under glaze conditions has a transparent glass body and a hard texture which is not easy to break, and a strong adhesion to the object surface. Therefore, it can cause major damage to bridge cables. In the glaze ice case, since the two-phase flow was coupled at the entrance and very complex to express, the left side was only set as the air velocity entrance to simulate the natural wind and the upper surface was set as the velocity inlet of the water droplets to simulate the rainfall. Moreover, the right side was the mixed-phase pressure outlet, and the other three sides were set to be symmetrical. The boundary conditions for glaze ice are shown in Figure 8. The mesh settings and qualities were nearly same as that in the rime climate conditions, with a total mesh number of 1.3 million, a minimum grid volume of 4.4 × 10 −6 m 3 , and a maximum grid volume of 9.1 × 10 −4 m 3 . The mesh layout is shown in Figure 9. The calculation time step was set to be 0.01 s. However, for the glaze ice case, the deformation of each analysis step was multiplied by 200 and restricted so as not to exceed 0.08 mm. The total number of time steps was 5000, which corresponded to 10,000 s (2.8 h) for actual time. Similar to the rime condition, the ice accretion of the glaze was generally logarithmically related to the time step. The growth of ice accretion was not significant at the setup step.

Validation of Numerical Method
In order to verify the proposed three-dimensional ice accretion model, the numerical results were compared with the cylindrical icing experiment results by Koss et al. [21].
The experimental temperatures were −5 °C, −3 °C and −1 °C, and the wind speed was 30 m/s, which was similar to the rime climate condition. The cylinder diameter was 3.81 cm, and 50 cm in length. The numerical simulation was under the same condition for the calculation, and Figure 10 shows the results of ice accretion on circular cylinders after 30 min. When the temperature was −3 °C and −1 °C, the numerical simulation results were generally consistent with the icing experiment. When the temperature was −5 °C, the numerical results had little differences with the experiment. This was due to the certain deviation of convection heat dissipation coefficient in the simulation and the actual results. Considering that the temperature of freezing rain weather is generally between −0 °C and −5 °C, it is reasonable to use this method to simulate the ice accretion on bridge cables.

Flow Field in the Rime Climate
To illustrate the characteristics of the flow field in the simulated rime environment, this paper shows a simulated flow field with a wind speed of 2 m/s, a droplet diameter of 10 μm, and water content of 5 × 10 −6 in the air. The water content of the whole calculation domain was initialed to be zero at the beginning. Thus, only the entrance had a liquid phase when the calculation started, as seen by the water content shown in Figure 11a. As time goes on, the water droplets passed through the location of the bridge cable, and finally left the computational domain from the pressure outlet in the right face. Figure 11b gives the liquid phase distribution at 20 s after the calculation. At this time, the water content at the windward side was uniform. However, on the leeward side, the water content had some disturbances due to air turbulence. Since the stable stage of the flow field around the cable was achieved, this moment was used as the initial stage for the simulation of rime ice accumulation on bridge cables.

Rime Iced Geometries of Three-Dimensional Bridge Cables
Since the wind speed and environment temperature determine the final iced geometries of bridge cables, the simulation process was carried out under three different wind speeds (v = 1 m/s, 2 m/s, and 5 m/s), and each wind speed was considered for the four common rime temperatures, which were −0.1 °C, −0.5 °C, −1 °C and −2 °C. The diameter of cable was 16 cm. Figure 12 shows the relationship between ice thicknesses and temperature at the wind speed of 1m/s. It is obvious that the iced geometries had no remarkable differences for the three temperatures below −0.5 °C, and the icing was mainly on the windward side (−75° < θ < 75°). Moreover, the ice thickness of −0.1 °C was relatively uniform at each position of the cross-section. Due to relatively lower wind speed, the water droplets froze in the windward region when the temperature was below −0.5 °C. However, the water droplets were gradually frozen in the higher temperatures when blown to the leeward side, instead of completely freezing on the windward side.  Figure 13 gives the final iced geometries at various temperatures under the wind speed of 2 m/s. As the wind speed increased, more water droplets impinged the cable surface per unit time, making the ice thicknesses larger than that of the lower wind speed at all temperatures except for −0.1 °C. Also caused by increasing wind speed, more droplets with high kinetic energy were received on the cable surface, but they were difficult to freeze for relatively higher temperatures. This resulted in the evaporation or shedding of Icing thickness(mm)

Iced Geometries at v = 2 m/s
Overlapping when temperatures below -0.5 °C the raindrops, which can be concluded from the −0.1 °C case; the ice accumulation part only concentrated on the windward side from −30° to 30°. For the temperatures of −1 °C and -2 °C, the iced geometries were similar to that of lower wind speed, and the difference was the thickness. As for the temperature of −0.5 °C, the droplets could not fully freeze at the windward side and were gradually frozen while flowing along the cable surface.

Iced Geometries at v = 5 m/s
The ice thicknesses at a wind speed of 5 m/s with different temperatures are shown in Figure 14. The remarkable difference, when compared to previous wind speed cases, was the distribution at different temperatures having different patterns. More water droplets impacted the cable surface with the rising wind speed, and the total ice thicknesses all increased under −0.5 °C, −1 °C, and −2 °C. The maximum thickness at −2 °C was above 50 mm, and due to the low temperature, the icing zone was still concentrated on the windward part. However, for −0.5 °C and −1 °C, the ice accumulation zone was around the cable surface, formed while the droplets were blown back to the leeward side. However, ice accumulation could hardly be found at the highest temperature −0.1 °C.

Aerodynamic Characteristics
The icing changed the original shape of the bridge cables and affected their original aerodynamic characteristics, which could result in aerodynamic instability. Based on the obtained iced shape of the bridge cables, three typical iced shape were selected to investigate their aerodynamic characteristics: R1 (−0.5 °C, 2 m/s), R2 (−1 °C, 2 m/s), and R3 (−2 °C, 2 m/s). Their aerodynamic coefficients were obtained, and the critical wind speed and wind attack angle of the cable galloping were further studied.

Aerodynamic Coefficients
A two-dimensional model was used to analyze the aerodynamic coefficients of the bridge cables under rime ice conditions. The calculation domain was set as 16 m long and 8 m wide. The left boundary was the velocity inlet, the incoming speed was 20 m/s, and the turbulence intensity was taken as 0.2%. The right boundary was the pressure outlet, the upper and lower boundaries were symmetrical, and the bridge cable was simulated by the non-slip wall boundary condition. The wind attack angle varied from 0° to 180°, and the step length of the non-steady-state analysis was set as 0.0005 s.
The comparison of the drag coefficients of the three models is shown in Figure 15a. It can be seen that as the thickness of the icing increased, the maximum value of the drag coefficient became larger. The comparison of lift coefficients is shown in Figure 15b, and these three models hold negative slopes. With the decrease in ambient temperature, the fluctuation of the lift coefficient curve increased.

Critical Wind Speed and Wind Attack Angle of Galloping
To determine the critical wind speed and wind attack angle of galloping, the bridge cable was simulated to have a mass of 100 kg/m, a frequency of 0.4 Hz, a damping ratio of 0.01, and a diameter of 16 cm. For the aforementioned three bridge cable models in the rime condition, the variation of Den Hartog coefficients with the wind attack angle are illustrated in Figure 16. Figure 16a shows that there was no galloping risk for the R1 model, because the Den Hartog coefficient was greater than 0. However, when the wind attack angle was 38° and 106°, the Den Hartog coefficient had a minimum value. Figure  16b shows that the R2 model had possible galloping vibrations at the wind attack angles of 16°, 66°, 115°, and 180°. The calculated critical wind speeds for these four wind attack angles are shown in Table 1. The R3 model had the risk of galloping vibration at three wind speeds (35°, 86°, and 156°). It can be seen that when the wind attack angle was 35°, the critical wind speed was relatively low. This high potential of galloping vibrations should be paid more attention. The corresponding results are shown in Figure16c and Table 2.

Simulation of the Glaze Climate
Since the setting of the boundary conditions, the gaseous phase and the liquid phase entered the flow field from different surfaces, which was inconsistent with the fact that the incoming flow was rich in super-cooled droplets. Therefore, in order to eliminate this deviation, the cable was located far away from both inlets where the two-phase flow was fully mixed. Initially, the air moisture content of the whole calculation domain was zero. When the water droplets entered the calculation domain, only the upper part of the calculation domain had water droplets, as shown in Figure 17a. Since the left side was the velocity inlet of air, the moisture content near the left entrance was always 0. Droplets were continuously falling and accumulating on the ground and eventually were blown out from the pressure outlet, as shown in Figure 17b. Because of the falling variation rate of raindrops with different sizes, 20 s after the rainfall start was set to be the initial time of glaze ice case, which was the stable stage of rainfall at this time. From the beginning stage of rainfall to its stabilization, the velocities of wind and water droplets were constantly changing. The liquid phase velocity field is shown in Figure 18. Due to the wind blowing from the left side, the water content near the left surface was zero. Water droplets received their energy from the wind to flow to the right; hence, the velocity gradually increased from left to right. The airflow field after the stable rainfall process is shown in Figure 19. It shows that the wind speed near the upper side was zero, since the inlet of the liquid phase acted as a wall for the gaseous phase. Due to the cable being located away from the liquid phase inlet, the flow field around it was not affected. Figure 19. Velocity of airflow field at stable stage.

Glaze Iced Geometries of Three-Dimensional Bridge Cables
The glaze ice simulation was also carried out under three different wind speeds (v = 1 m/s, 2 m/s, and 5 m/s), and each wind speed was considered for the four common rime temperatures, which are −0.1 °C, −0.5 °C, −1 °C, and −2 °C. Since the glaze ice appears at relatively higher temperatures, the simulation temperatures were slightly different from the rime ice case. Figures 20 and 21 show the three-dimensional iced geometries on the cable surface when the wind speed was 1m/s. At temperature conditions except −0.1 °C, ice accumulation mainly happened in −30° to 180°, the water droplet receiving zone. At the temperature of −0.1 °C, the icing zone was from −180° to 30°. This is because when the temperature was high the water could not fully freeze in the droplet receiving area, and instead they gradually froze when flowing back to the downstream side. More raindrops quickly froze at the droplet receiving zone with the drop in temperature, and the maximum ice thickness also increased. give the general situation of the ice accumulation on the cable surface when the wind speed was 2 m/s. Since the wind speed increased, the surface received more water droplets per unit of time. At the temperature of −1 °C and −2 °C, the icing range was wider, and ice was thicker than that of lower wind speed. When the temperature was high, only a small part of the cable surface was covered by ice, because the water on the cable was difficult to freeze and accumulate, due to high wind speed and high temperature. Icing thickness(mm)

Iced Geometries at v = 5 m/s
The iced geometries at different temperatures are shown in Figures 24 and 25 for when the wind speed rose to 5 m/s. At the temperature of −0.1 °C and −0.5 °C, almost no ice accumulation could be found on the cable surface. While the total ice thickness was larger than that of lower wind speed, this was due to the increase in wind speed and more water droplets impacting the cable surface. When the temperature dropped to −1 °C and −2 °C, the icing range was wider than that of 2 m/s, but the ice thickness was smaller. This was due to the complex coupling relationship between the wind speed, the kinetic energy of the water droplets, the cooling rate, and the receiving rate of the water droplets.

Aerodynamic Characteristics
The numerical model setups used in the glaze ice condition were the same as those in rime ice conditions. Iced shapes were studied under three different glaze ice conditions, that is, G1 (−2 °C, 1 m/s), G2 (−2 °C, 1 m/s), and G3 (−2 °C, 1 m/s). The aerodynamic coefficients, critical wind speed, and wind attack angle of the cable galloping are described below.

Aerodynamic Coefficients
The average aerodynamic coefficients of three models are shown in Figure 26a. The drag coefficient of the G1 model had a convex shape in the range of 80° to 160°, and the G2 and G3 models had two convex shapes in the range of 120° to 280°. The three models all had concave shapes in the range of 200° to 240°. The comparison of the lift coefficients of the three models is shown in Figure 26b. The lift coefficient of the G1 model decreased rapidly in the range of 60° to 120°. In the G2 model, the lift coefficient decreased rapidly

Critical Wind Speed and Wind Attack Angle of Galloping
The variation of the Den Hartog coefficient of each icing model with the angle of attack of the wind is shown in Figure 27. Each model had the potential for galloping variations. Table 3 lists three critical wind speeds of the G1 models, and indicates that the galloping vibrations may occur with the lower critical wind speed at the wind attack angle of 66° and 226°. The G2 model could have galloping vibrations at the wind speed of 17.5 m/s when the wind attack angle is 346°, which is described in the Table 4. Under the low wind speed at 28 m/s and 21.8 m/s, the G3 model (Table 5)

Rime Iced Geometries under Different Temperature and Wind Speed
As illustrated in Table 6, under the temperature of −0.1 °C, the ice accretion thickness decreased with the increase in wind speed. When the wind speed was 1 m/s, most of the positions of the cable were covered with ice, and when the wind speed increased to 2 m/s and 5 m/s, only a small part of the cable surface was covered with ice. This was due to the fact that too many water droplets with high kinetic energy were received by the cable surface at the same time, resulting in evaporation or shedding of the water droplets at this high temperature. Under the temperature of −0.5 °C, the ice-covered range on the cable surface increased with the increase in wind speed. This was due to the increase in water droplet rate. When the wind speed was 5 m/s, the whole cable surface was covered with ice.
When the temperature dropped to −1 °C, the ice cover pattern at the wind speed of 1 m/s was similar with that at a temperature of −0.5 °C. The ice cover was only affected by the distribution of water droplets, and when the wind speed increased the icing thickness increased, while the ice cover range decreased. This was because the distance supercooled moisture can flow before freezing was shortened.
When the temperature continued decrease to −2 °C, there was no change in the ice cover pattern on the cable surface when the wind speed was 1 m/s. The ice cover pattern on the surface of the cable was not smooth when the wind speed was 5 m/s, which was due to the turbulence caused by the faster wind speed.

Glaze Iced Geometries under Different Temperature and Wind Speed
As illustrated in Table 7, under various wind speeds the ice accretion thickness increased as the temperature decreased. The most significant difference with the rime condition was that the final iced geometry of rime was generally symmetric, while the glaze iced geometries were various at different angles. This was due to the water phase inlets being different under these two climate conditions: water droplets are small and suspended in the air under rime conditions, while the diameter of water droplets for glaze condition is larger. For the same reason, the ice covering range of the glaze condition was more concentrated, while the rime iced covering was dispersed.

Cable Galloping under Glaze and Rime Conditions
Cable galloping is a complex phenomenon of fluid-structure interaction. Galloping occurs when a cable presents an asymmetric profile to a moderate to high wind flow. The above sections studied cable galloping issues for different rime and glaze conditions. Table 8 summarizes the minimal critical wind speed and associated wind attack angles for each model. It can be seen that cables in the rime conditions had lower critical wind speed than in glaze conditions. The amplitude and wind speed for the wind attack angle 66° for model R3 are shown in Figure 28. When the wind speed was 12.5 m/s, the amplitude of the cable could reach 0.5 m, and when the wind speed reached 20 m/s, the cable amplitude increased to about 1 m/s. With the wind speed increasing, the galloping amplitude linearly increased. This could lead to a collision between the cables. The time for stabilizing the galloping vibration was 240 s at 20 m/s, 160 s at 30 m/s and 100 s at 40 m/s. The amplitude and wind speed for the wind attack angle 346° for model G2 are shown in Figure 29. When the wind speed was 25 m/s, the amplitude of the vibration was about 0.7 m, and the time required for stabilizing the galloping vibration was about 640 s. When the wind speed was 30 m/s, the amplitude of the vibration was about 1.1 m, and the time for stabilizing the galloping vibration was about 400 s. When the wind speed was 40 m/s, the amplitude of the vibration was about 1.7 m, and the time for stabilizing the galloping vibration was about 220 s.

Conclusions
This paper proposes a numerical simulation framework based on the three-dimensional Messinger theory to simulate the ice accumulation on bridge cables under two typical freezing rain climate conditions: rime ice and glaze ice. Based on the measured meteorological data of Yunnan-Guizhou Plateau in China, the realistic climate conditions of freezing rains were achieved. The impact of wind speed and temperature on the shape of ice accumulation in rime and glaze climate was studied. The associated aerodynamic characteristics were also investigated in this paper. Several conclusions can be drawn: (1) To accurately predict the final iced shape, a three-dimensional ice accumulation model based on Messinger theory was proposed, which was programmed and embedded into computational fluid software. The proposed model can consider factors such as wind speed, temperature, air moisture content, and raindrop falling speed. By controlling these factors, the possible iced shapes of bridge cables under different weather conditions can be obtained. (2) In rime weather, within the range of the diameter of the droplets, the increase in the diameter of the droplets did not have a significant impact on the icing. Wind speed had two effects on the formation of icing: when the temperature was higher, excessive wind speed reduced the moisture on the surface of the stay cable; when the temperature was lower, the increase in wind speed brought more water to the cable and increased icing. (3) In glaze weather, the decrease in temperature leads to the formation of ice coating, and affects the final iced geometries by changing the amount of free water on the surface of the bridge cable. The effect of wind speed was similar to that of rime weather. The impact of the diameter of raindrops on the shape of the icing was different to the rime weather. This is because large raindrops are more difficult to flow around under the influence of wind, and large raindrops fall faster than small ones. (4) For the aerodynamic characteristics of bridge cables with different iced geometries, the lift coefficients and drag coefficients were obtained, and the critical wind speed and wind attack angle of galloping were determined by Den Hartog criterion. These preliminary aerodynamic simulations under the typical freezing rain conditions can provide the basis for the wind-induced vibration analysis of the whole structure.
This study mainly proposes a numerical method for simulating the ice accretion on cables or circular cylinders in rime or glaze conditions under a certain inclination angle. Future work will focus on exploring other complex geometries under various inclination angles of engineering interest.  Data Availability Statement: The data of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.