Calculation method for predicting icing growth of wire

Abstract

The invention discloses a calculation method for predicting conductor icing growth, which is characterized by comprising the following steps of: 1, dividing a flow field calculation domain into a potential flow region and a near-wall region, and further obtaining air velocity distribution of all points to be solved in the calculation domain; step 2, searching a limit collision point on the surface of the zero-collision conductor by combining a water drop motion control equation and a dichotomy iterative approximation mode; 3, substituting the water drop collision coefficient and the air velocity distribution into the thermodynamic model, and establishing a mass conservation model and an energy conservation equation in a conductor icing surface unit body based on the thermodynamic model in a Messinger infinitesimal control body; 4, the icing thickness corresponding to each infinitesimal is calculated, and coordinates of a new ice shape boundary are obtained in the icing increasing direction; and possible electric power system accidents can be predicted in advance, so that the stable operation of the electric power system can be ensured, and unnecessary social and economic losses can be reduced.

Description

Technical Field

[0001] The invention belongs to the technical field of transmission line icing monitoring, and in particular relates to a calculation method for predicting the growth of conductor icing. Background Art

[0002] In recent years, transmission line icing disasters have frequently occurred around the world. This natural phenomenon has become a major challenge affecting the safe and stable operation of power systems. These disasters not only directly threaten the reliability of the power grid, but also frequently cause large-scale power outages due to their strong suddenness and wide impact range, causing immeasurable losses to residents' lives, industrial production and socio-economic development. In view of this, we deeply explore the physical process of conductor icing and develop an accurate and efficient numerical calculation method for conductor icing. This provides a scientific basis for the anti-icing and de-icing design of transmission lines, predicts the evolution trend of conductor icing under different climatic conditions, and thus takes countermeasures in advance to effectively reduce the risk of power outages caused by natural icing. Summary of the invention

[0003] The technical problem to be solved by the present invention is: to provide a calculation method for predicting the growth of conductor ice coverage, so as to effectively predict ice coverage and make early predictions of possible power system accidents, thereby ensuring the safe and stable operation of the power system and reducing unnecessary social and economic losses.

[0004] Technical solution of the present invention:

[0005] A calculation method for predicting the growth of ice covering on a conductor, the method comprising:

[0006] Step 1: Divide the flow field calculation domain into the potential flow area and the near-wall area, and then obtain the air velocity distribution of all the points to be calculated in the calculation domain;

[0007] Step 2, combining the water drop motion control equation and the binary division iterative approximation method to find the limit collision point on the zero-collision wire surface;

[0008] Step 3: Substitute the water drop collision coefficient and air velocity distribution into the thermodynamic model, and establish the mass conservation model and energy conservation equation of the conductor ice-covered surface unit body based on the thermodynamic model in the Messinger micro-element control body;

[0009] Step 4: Calculate the ice thickness corresponding to each microelement, and then obtain the coordinates of the new ice shape boundary along the ice growth direction.

[0010] Methods for calculating the air velocity distribution at all points in the domain include:

[0011] Potential flow field solution:

[0012]

[0013] is the velocity potential function of the source point Q; q is the normal derivative of the velocity potential function of the object boundary; n Q Represents the unit normal vector of point Q; n x is the component of the unit normal vector on the boundary in the x-axis direction; v is the incoming flow velocity at infinity, and the partial derivative of the total potential flow at the point to be solved with respect to the coordinate axis is the velocity component of each point to be solved;

[0014] By calculating the value of the Reynolds number R of the boundary layer momentum thickness, the position of the boundary layer flow transitioning from laminar flow to turbulent flow can be determined; the Karman boundary layer momentum integral relationship is obtained:

[0015]

[0016] On the outer boundary of the boundary layer, the fluid flows in an ideal potential manner, and then:

[0017]

[0018] Where δ is the thickness of the boundary layer; p is the pressure, and Ue is the boundary layer velocity;

[0019] In the near-wall region, Vh is the velocity at a vertical distance h from the wall. The boundary restriction conditions include: when the wall velocity is 0, that is, h = 0, Vh = 0; on the outer boundary of the boundary layer, the velocity is the flow velocity in the potential flow zone; on the outer boundary of the boundary layer, there is no velocity gradient at this time, that is, the velocity gradient is 0; combining the above conditions and using the Bolhausen method, the velocity distribution in the boundary layer is approximated by a fourth-order polynomial, and the following is obtained:

[0020]

[0021] Where ds is the infinitesimal unit length on the boundary layer;

[0022] The expression of friction shear stress τ0 is:

[0023]

[0024] The expression of dimensionless pressure gradient coefficient B is:

[0025]

[0026] Laminar boundary layer momentum thickness θ l The expression is:

[0027]

[0028] The velocity in the turbulent zone changes approximately exponentially, and the empirical expression is:

[0029]

[0030] The friction shear stress and friction coefficient are:

[0031]

[0032] Momentum thickness of the boundary layer in the turbulent zone θ t The expression is

[0033]

[0034] In the formula, str represents the position of the boundary between laminar flow and turbulent flow, is the momentum thickness of the boundary layer on the conductor surface; μ is the air viscosity coefficient, ρ is the fluid density; U e is the boundary layer velocity, and the air velocity distribution in the near-wall layer is obtained after calculation.

[0035] Methods for finding the limit collision point on the zero-collision wire surface include:

[0036] Assume that when the initial points are i and j, the trajectories collide with the wires, and the collision angles are

[0037] θ i ,θ j ; When the initial point is k, the trajectory bypasses the wire without collision. If the trajectory is tangent to the wire, this trajectory is what we are looking for, and the corresponding limit collision angle is θ max； Finding the limiting trajectory of the water droplet with zero collision with the wire is a gradually approaching numerical calculation process; the general principle for finding the limiting trajectory of the upper zero collision wire is: assuming that the motion trajectory calculated this time collides with the wire, the initial coordinate point of the water droplet emission will move upwards during the next iteration; conversely, assuming that the motion trajectory calculated this time does not collide with the wire, the initial coordinate point of the water droplet emission will move downwards during the next iteration.

[0038] The binary search iterative method is used to solve the coordinates of the limit point on the surface of the water droplet zero collision wire, which includes:

[0039] Step 2.1, select the initial coordinate points of the two trajectories, the initial y-direction coordinate of one trajectory is 0, that is, on the x-axis, recorded as Y0, and the initial y-direction coordinate of the other trajectory is 1.1 times the radius, recorded as Y1. Except for the case where the collision coefficient is 0, the water droplets emitted from the Y0 position will definitely collide with the wire, in this case, set flag_0 = -1, and the water droplets emitted from the Y1 position will definitely bypass the upper surface of the wire, in this case, set flag_1 = 1;

[0040] Step 2.2, determine whether the water droplet motion trajectory emitted from the position with the ordinate (Y0+Y1) / 2 collides with the surface of the ice-covered wire. If a collision occurs, the midpoint flag is immediately set to -1, that is, flag_m=-1; since flag_m×flag_1<0, the water droplet emission position corresponding to the upper limit trajectory should be between (Y0+Y1) / 2 and Y1, and Y0 needs to be reassigned, that is, Y0=(Y0+Y1) / 2; assuming that no collision occurs, the midpoint flag is immediately set to 1, that is, flag_m=1; since flag_m×flag_0<0, the water droplet emission position corresponding to the upper limit trajectory should be between Y0 and (Y0+Y1) / 2, and Y1 needs to be reassigned, that is, Y1=(Y0+Y1) / 2;

[0041] Step 2.3, iterative stop condition for searching zero-collision wire trajectory: when the condition |Y0-Y1|≤10 -4 If it is not true, the iteration calculation will continue. Otherwise, when the above conditions are met, the iteration will stop.

[0042] Step 2.4, calculate the straight-line distance from the coordinate point on the trajectory to the arc microelement, that is, find the minimum distance and the corresponding microelement unit, which is the collision limit area.

[0043] The mass conservation model and energy conservation equations for the unit body of the conductor ice-covered surface are established as follows:

[0044] The mass and energy balance in the wire surface microelement control body is obtained as follows:

[0045]

[0046] In the formula, m im is the mass of the water droplet that hits the microelement; m in and m out are the mass of liquid water flowing into and out of the current control volume microelement respectively; m eva is the mass of liquid water evaporated; m ice Q is the final mass of frozen water droplets; f Energy released by freezing of water droplets; Q h is the energy generated by the collision of water droplets; Q r is the heat of friction between air and wall; Q e is the evaporation energy of liquid water; Qc is the air convection heat transfer; Q in , Q out are the energies flowing into and out of the microelement of the current control volume respectively; the specific expressions are as follows:

[0047]

[0048] The energy released when a water drop freezes is:

[0049] Q f =m ice [c i (T f -T s )-L f ]

[0050] The energy brought by the water droplet hitting the microelement surface is:

[0051]

[0052] The heat loss due to friction between air and wall is:

[0053]

[0054] The energy lost by evaporation of liquid water is:

[0055] Q e =m eva [C d (T s -T f )+L e ]

[0056] The natural convection heat transfer is:

[0057] Q c =h c A(T s -T)

[0058] Where, Ts is the surface temperature of the ice-covered conductor; T f is the freezing temperature of water droplets; T is the ambient temperature, K, L f , L e C is the latent heat of ice melting and the latent heat of ice sublimation and evaporation respectively; a , C d , C i are the specific heats of air, water and ice respectively; h c is the convective heat transfer coefficient; U e is the boundary layer velocity; V c The water drop collision speed; A is the microelement length;

[0059] r is the recovery coefficient of kinetic heating, expressed as:

[0060]

[0061] Where Pr is the Prandtl number, n0 = 1 / 2 in the laminar region; n0 = 1 / 3 in the turbulent region; V c is the wind speed; m im Expressed as:

[0062] m im =β·LWC·V c ·A

[0063] β represents the local collision coefficient of the infinitesimal; LWC is the liquid water content, and the freezing coefficient f is:

[0064]

[0065] In the formula, L ew is the Lewis number, with a value of 1; R v is the water vapor volume constant, with a value of 461 J / (Kg·k); Ф is the relative humidity of the atmosphere, usually taking a value of 100% in the icing environment; P sat (Ts) and P sat (T) are the saturated water vapor pressures corresponding to the wire surface temperature T s and the ambient temperature T respectively; the unknown parameters are m in , f, and Ts, and an iterative method is used to solve them.

[0066] The method of solving by iteration includes: initially assuming that the wire surface temperature Ts is 0. If the calculation result is 0 < f < 1, then Ts = 273.15 K; if the calculation result is f < 0, then let f = 0 and then solve for Ts inversely; if the calculation result is f > 1, then let f = 1 and then solve for Ts inversely; the numerical calculation starts from the stagnation point and calculates in both the upper and lower directions.

[0067] When calculating the ice thickness corresponding to each infinitesimal and then obtaining the coordinates of the new ice shape boundary along the ice growth direction, the least squares method is used to fit the discrete points of the new boundary directionally to obtain the boundary curve function; the boundary curve expression will be used to update the boundary conditions for the iterative calculation of the next time step; the value of the time step is adjusted according to the actual icing degree and the ice growth rate; specifically including: when the icing degree increases by about 3.4%, it is regarded as the completion of one iterative calculation, and the time step is calculated and determined in combination with the ice growth rate under this condition. The ice thickness corresponding to the icing degree of 3.4% is 0.034D, and the optimal time step is determined in combination with the ice growth rate to determine the optimal time step.

[0068] Advantages of the present invention:

[0069] First, according to the wire diameter and wind speed conditions, the present invention calculates the distribution of the air flow field velocity, and then combines the median volume diameter of the water droplets to obtain the distribution law of the collision coefficient of the water droplets on the wire surface and the velocity vector of the water droplets colliding with the wire surface. Substituting the above calculation results and the liquid water content in the air into the thermodynamic model, the convective heat transfer coefficient, the local freezing coefficient, and the infinitesimal freezing mass are calculated, and the ice shape of the wire icing growth is calculated. This method can effectively predict icing, make an advance prediction of possible power system accidents, and thus ensure the stable operation of the power system, and further reduce unnecessary social and economic losses. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 It is a schematic diagram of a specific implementation process of the present invention;

[0071] Figure 2 A schematic diagram of solving the trajectory of a water droplet zero-collision wire surface in a specific embodiment of the present invention;

[0072] Figure 3 It is a schematic diagram of a method for finding the upper and lower limit points of a water droplet zero-collision wire surface in a specific embodiment of the present invention;

[0073] Figure 4 It is a schematic diagram of the mass and energy conservation model of the ice-covered surface control body in a specific embodiment of the present invention;

[0074] Figure 5 It is a schematic diagram of a method for determining an adaptive time step according to a specific embodiment of the present invention. DETAILED DESCRIPTION

[0075] The present invention provides a method for calculating the growth of ice on a conductor, comprising the following steps:

[0076] S1: Divide the flow field calculation domain into the potential flow area and the near-wall area, and then obtain the air velocity distribution of all the points to be calculated in the calculation domain. The following formula is the potential flow field solution method:

[0077]

[0078] is the velocity potential function of the source point Q; q is the normal derivative of the velocity potential function of the object boundary; n Q represents the unit normal vector of point Q; n x is the component of the unit normal vector on the boundary in the x-axis direction; v is the velocity of the incoming flow at infinity. The partial derivative of the total potential flow at the point to be solved with respect to the coordinate axis is the velocity component of each point to be solved.

[0079] By calculating the value of the Reynolds number R of the boundary layer momentum thickness, the position of the dividing point where the boundary layer flow transitions from laminar flow to turbulent flow can be determined.

[0080] According to the Karman boundary layer momentum integral relationship:

[0081]

[0082] On the outer boundary of the boundary layer, the fluid flows in an ideal potential manner, and then:

[0083]

[0084] Where δ is the thickness of the boundary layer; p is the pressure, and Ue is the boundary layer velocity.

[0085] In the near-wall region, Vh is the velocity at a vertical distance h from the wall. The boundary restriction conditions are as follows: a) When the wall velocity is 0, that is, h = 0, Vh = 0; b) On the outer boundary of the boundary layer, the velocity is the flow velocity in the potential flow zone; c) On the outer boundary of the boundary layer, there is no velocity gradient at this time, that is, the velocity gradient is 0; Combining the above conditions and using the Pohlhausen method, the velocity distribution in the boundary layer is approximated by a fourth-order polynomial, and we can get:

[0086]

[0087] Where ds is the infinitesimal unit length on the boundary layer.

[0088] The expression of friction shear stress τ0 is:

[0089]

[0090] The expression of dimensionless pressure gradient coefficient B is:

[0091]

[0092] Laminar boundary layer momentum thickness θ l The expression is:

[0093]

[0094] The velocity in the turbulent zone changes approximately exponentially, and the empirical expression is:

[0095]

[0096] The friction shear stress and friction coefficient are:

[0097]

[0098] Momentum thickness of the boundary layer in the turbulent zone θ t The expression is

[0099]

[0100] In the formula, str represents the position of the boundary between laminar flow and turbulent flow, is the momentum thickness of the boundary layer on the conductor surface; μ is the air viscosity coefficient, ρ is the fluid density; U e is the boundary layer velocity. The air velocity distribution in the near-wall layer is obtained after calculation.

[0101] S2: Combine the water drop motion control equation and the bisection iterative approximation method to find the limit collision point on the zero-collision wire surface. Figure 2 As shown in the figure, assuming that when the initial points are i and j, their trajectories collide with the wire, and the collision angles are θ i ,θj When the initial point is k, its trajectory bypasses the wire without collision. If the trajectory is exactly tangent to the wire, this trajectory is what we are looking for (as shown by the red dotted line), and its corresponding limit collision angle is θ max Finding the limit trajectory of the water droplet with zero collision with the wire is a gradually approaching numerical calculation process. Taking the upper trajectory as an example, the general principle of finding the limit trajectory of the upper zero collision wire is: assuming that the motion trajectory calculated this time collides with the wire, the initial coordinate point of the water droplet emission should move up in the next iteration. Conversely, assuming that the motion trajectory calculated this time does not collide with the wire, the position of the initial coordinate point of the water droplet emission should move down in the next iteration.

[0102] For the sake of convenience and speed, the present invention adopts a binary search iterative method to solve the coordinates of the limit point on the surface of the water drop zero collision wire. Taking the solution of the coordinates of the limit point of the upper zero collision wire as an example, the specific ideas are as follows:

[0103] (1) Select the initial coordinate points of the two trajectories. The initial y-direction coordinate of one of the trajectories is 0, that is, on the x-axis, and for convenience, it is recorded as Y0. The initial y-direction coordinate of the other trajectory is 1.1 times the radius, recorded as Y1. Except for the case where the collision coefficient is 0, the water droplets emitted from the Y0 position will definitely collide with the wire, in which case flag_0 = -1 is set. The water droplets emitted from the Y1 position will definitely bypass the upper surface of the wire, in which case flag_1 = 1 is set.

[0104] (2) Determine whether the trajectory of the water droplet emitted from the position with the ordinate (Y0+Y1) / 2 collides with the surface of the ice-covered conductor. Assuming a collision, the midpoint flag is immediately set to -1, that is, flag_m = -1; since flag_m×flag_1<0, the water droplet emission position corresponding to the upper limit trajectory should be between (Y0+Y1) / 2 and Y1. Y0 needs to be reassigned, that is, Y0 = (Y0+Y1) / 2. Assuming no collision occurs, the midpoint flag is immediately set to 1, that is, flag_m = 1; since flag_m×flag_0<0, the water droplet emission position corresponding to the upper limit trajectory should be between Y0 and (Y0+Y1) / 2. Y1 needs to be reassigned, that is, Y1 = (Y0+Y1) / 2.

[0105] (3) Iterative stop condition for searching for zero-collision wire trajectories: When the condition |Y0-Y1|≤10 -4 If it is not true, continue to repeat step 2 for iterative calculation. Otherwise, when the above conditions are met, the iteration stops.

[0106] (4) Since the coordinates of all micro-elements on the wire boundary are known, the straight-line distance from the coordinate point on the trajectory to the arc micro-element can be calculated, and the minimum distance and its corresponding micro-element unit can be found. This micro-element unit is the collision limit area. Since the micro-element division interval is originally very small, the specific coordinates of the extreme point can be replaced by the coordinate value of the midpoint of the micro-element.

[0107] S3: Substitute the water drop collision coefficient, air velocity distribution and other results calculated in step 1 and step 2 into the thermodynamic model calculation. Based on the thermodynamic model in the Messinger micro-element control body, the mass conservation model and energy conservation equation in the unit body of the conductor ice-covered surface can be established, such as Figure 4 shown.

[0108] The mass and energy balance in the wire surface microelement control body can be obtained:

[0109]

[0110] In the formula, m im is the mass of the water droplet that hits the microelement; m in and m out are the mass of liquid water flowing into and out of the current control volume microelement respectively; m eva is the mass of liquid water evaporated; m ice is the final mass of the frozen water drop. The energy released by the freezing of the water drop is Q f ; Energy generated by water drop collision Q h ; Heat Q caused by friction between air and wall r ; Liquid water evaporation energy Q e ; Air convection heat transfer Qc; Q in , Q out are the energies flowing into and out of the current control volume element respectively.

[0111] The specific expressions corresponding to the relevant parameters are as follows:

[0112]

[0113] The energy released when a water drop freezes is:

[0114] Q f =m ice [c i (T f -T s )-L f ]

[0115] The energy brought by the water droplet hitting the microelement surface is:

[0116]

[0117] The heat loss due to friction between air and wall is:

[0118]

[0119] The energy lost by evaporation of liquid water is:

[0120] Q e =m eva [C d (T s -T f )+L e ]

[0121] The natural convection heat transfer is:

[0122] Q c =h c A(T s -T)

[0123] Where, Ts is the surface temperature of the ice-covered conductor, K; T f is the freezing temperature of water droplets, which is 273.15K; T is the ambient temperature, K, L f , L e C is the latent heat of ice melting and the latent heat of ice sublimation and evaporation respectively; a , C d , C i are the specific heats of air, water and ice respectively; h c is the convective heat transfer coefficient; U e is the boundary layer velocity; V c The water drop collision speed; A is the infinitesimal length.

[0124] r is the recovery coefficient of kinetic heating, which can be expressed as:

[0125]

[0126] Where Pr is the Prandtl number, which is usually 0.72. In laminar flow n0 = 1 / 2; in turbulent flow n0 = 1 / 3; V c is the wind speed. m im It can be expressed as:

[0127] m im =β·LWC·V c ·A

[0128] β represents the local collision coefficient of the microelement; LWC is the liquid water content. The freezing coefficient f is:

[0129]

[0130] Where, L ew is the Lewis number, which takes the value of 1; R vis the volume constant of water vapor, which can be taken as 461 J / (Kg·k); Ф is the relative humidity of the atmosphere, usually taken as 100% in the icing environment; P sat (Ts) and P sat (T) are the saturated water vapor pressures corresponding to the surface temperature T s of the wire and the ambient temperature T respectively. The unknown parameters are m in , f and Ts. Here, an iterative method is used to solve. Specifically: initially assume that the surface temperature Ts of the wire is 0. If the calculation result is 0 < f < 1, then Ts = 273.15K; if the calculation result is f < 0, then f needs to be set to 0, and then Ts is solved inversely. If the calculation result is f > 1, then f needs to be set to 1, and then Ts is solved inversely. The numerical calculation starts from the stagnation point and calculates in both the upper and lower directions. Therefore, for the first microelement, m in = 0, Q in = 0.

[0131] S4: Calculate the ice thickness corresponding to each microelement, and then obtain the coordinates of the new ice shape boundary along the ice growth direction.

[0132] Use the least squares method to fit the discrete points of the new boundary directionally to obtain the boundary curve function. This boundary curve expression will be used to update the boundary conditions for the iterative calculation of the next time step. If the time step is too large, it will cause a large deviation in the ice shape of the wire. If the value is too small, it will lead to too many iterations, reduce the operation efficiency, and increase the possibility of error because the ice thickness increment is small each time. Considering the accuracy and calculation efficiency of the numerical calculation, the value of the time step should be adjusted according to the actual icing degree and ice growth rate. The present invention suggests that when the icing degree increases by about 3.4%, it can be regarded as the completion of one iterative calculation, and the time step is calculated and determined in combination with the ice growth rate under this condition. The specific determination method is as Figure 5 shown. The ice thickness corresponding to the icing degree of 3.4% is 0.034D, and the optimal time step can be determined in combination with the ice growth rate .

Claims

1. A calculation method for predicting the growth of ice covering on a conductor, characterized in that: The method comprises: Step 1: Divide the flow field calculation domain into the potential flow area and the near-wall area, and then obtain the air velocity distribution of all the points to be calculated in the calculation domain; Step 2, combining the water drop motion control equation and the binary division iterative approximation method to find the limit collision point on the zero-collision wire surface; Step 3: Substitute the water drop collision coefficient and air velocity distribution into the thermodynamic model, and establish the mass conservation model and energy conservation equation of the conductor ice-covered surface unit body based on the thermodynamic model in the Messinger micro-element control body; Step 4: Calculate the ice thickness corresponding to each microelement, and then obtain the coordinates of the new ice shape boundary along the ice growth direction.

2. A calculation method for predicting the growth of conductor ice cover according to claim 1, characterized in that: Methods for calculating the air velocity distribution at all points in the domain include: Potential flow field solution: is the velocity potential function of the source point Q; q is the normal derivative of the velocity potential function of the object boundary; n Q Represents the unit normal vector of point Q; n x is the component of the unit normal vector on the boundary in the x-axis direction; v is the incoming flow velocity at infinity, and the partial derivative of the total potential flow at the point to be solved with respect to the coordinate axis is the velocity component of each point to be solved; By calculating the value of the Reynolds number R of the boundary layer momentum thickness, the position of the boundary layer flow transitioning from laminar flow to turbulent flow can be determined; from the Karman boundary layer momentum integral relationship, we get: On the outer boundary of the boundary layer, the fluid flows in an ideal potential manner, and then: Where δ is the boundary layer thickness; p is the pressure, and Ue is the boundary layer velocity; In the near-wall region, Vh is the velocity at a vertical distance h from the wall. The boundary restriction conditions include: when the wall velocity is 0, that is, h = 0, Vh = 0; on the outer boundary of the boundary layer, the velocity is the flow velocity in the potential flow zone; on the outer boundary of the boundary layer, there is no velocity gradient at this time, that is, the velocity gradient is 0; combining the above conditions and using the Bolhausen method, the velocity distribution in the boundary layer is approximated by a fourth-order polynomial, and the following is obtained: Where ds is the infinitesimal unit length on the boundary layer; The expression of friction shear stress τ0 is: The expression of dimensionless pressure gradient coefficient B is: Laminar boundary layer momentum thickness θ l The expression is: The velocity in the turbulent zone changes approximately exponentially, and the empirical expression is: The friction shear stress and friction coefficient are: Momentum thickness of the boundary layer in the turbulent zone θ t The expression is In the formula, str represents the position of the boundary between laminar flow and turbulent flow, is the momentum thickness of the boundary layer on the conductor surface; μ is the air viscosity coefficient, ρ is the fluid density; U e is the boundary layer velocity, and the air velocity distribution in the near-wall layer is obtained after calculation.

3. A calculation method for predicting the growth of ice covering on a conductor according to claim 1, characterized in that: Methods for finding the limit collision point on the zero-collision wire surface include: Assume that when the initial points are i and j, the trajectories collide with the wires, and the collision angles are θ i ,θ j ; When the initial point is k, the trajectory bypasses the wire without collision. If the trajectory is tangent to the wire, this trajectory is what we are looking for, and the corresponding limit collision angle is θ max； Finding the limiting trajectory of the water droplet with zero collision with the wire is a gradually approaching numerical calculation process; the general principle for finding the limiting trajectory of the upper zero collision wire is: assuming that the motion trajectory calculated this time collides with the wire, the initial coordinate point of the water droplet emission will move upwards during the next iteration; conversely, assuming that the motion trajectory calculated this time does not collide with the wire, the initial coordinate point of the water droplet emission will move downwards during the next iteration.

4. A calculation method for predicting the growth of ice covering on a conductor according to claim 3, characterized in that: The binary search iterative method is used to solve the coordinates of the limit point on the surface of the water droplet zero collision wire, which includes: Step 2.1, select the initial coordinate points of the two trajectories, the initial y-direction coordinate of one trajectory is 0, that is, on the x-axis, recorded as Y0, and the initial y-direction coordinate of the other trajectory is 1.1 times the radius, recorded as Y1. Except for the case where the collision coefficient is 0, the water droplets emitted from the Y0 position will definitely collide with the wire, in this case, set flag_0 = -1, and the water droplets emitted from the Y1 position will definitely bypass the upper surface of the wire, in this case, set flag_1 = 1; Step 2.2, determine whether the water droplet motion trajectory emitted from the position with the ordinate (Y0+Y1) / 2 collides with the surface of the ice-covered wire. If a collision occurs, the midpoint flag is immediately set to -1, that is, flag_m=-1; since flag_m×flag_1<0, the water droplet emission position corresponding to the upper limit trajectory should be between (Y0+Y1) / 2 and Y1, and Y0 needs to be reassigned, that is, Y0=(Y0+Y1) / 2; assuming that no collision occurs, the midpoint flag is immediately set to 1, that is, flag_m=1; since flag_m×flag_0<0, the water droplet emission position corresponding to the upper limit trajectory should be between Y0 and (Y0+Y1) / 2, and Y1 needs to be reassigned, that is, Y1=(Y0+Y1) / 2; Step 2.3, iterative stop condition for searching zero-collision wire trajectory: when the condition |Y0-Y1|≤10 -4 If it is not true, the iteration calculation will continue. Otherwise, when the above conditions are met, the iteration will stop. Step 2.4: Calculate the straight-line distance from the coordinate points on the trajectory to the arc microelement, that is, find the minimum distance and the corresponding microelement unit, and this microelement unit is the collision limit area.

5. The calculation method for predicting the growth of ice covering on a conductor according to claim 1 is characterized in that: Establish the mass conservation model and energy conservation equation in the unit volume of the conductor icing surface, including: From the mass and energy balance in the microelement control volume on the conductor surface: In the formula, m im is the mass of the water droplet that hits the microelement; m in and m out are the mass of liquid water flowing into and out of the current control volume microelement respectively; m eva is the mass of liquid water evaporated; m ice Q is the final mass of frozen water droplets; f Energy released by freezing of water droplets; Q h is the energy generated by the collision of water droplets; Q r is the heat of friction between air and wall; Q e is the evaporation energy of liquid water; Qc is the air convection heat transfer; Q in , Q out are the energies flowing into and out of the microelement of the current control volume respectively; the specific expressions are as follows: The energy released by the freezing of water droplets is: Q f =m ice [c i (T f -T s )-L f ] The energy brought by the impact of water droplets on the microelement surface is: The heat loss due to the friction between air and the wall is: The energy loss due to the evaporation of liquid water is: Q e =m eva [C d (T s -T f )+L e ] The natural convection heat transfer is: Q c =h c A(T s -T) Where, Ts is the surface temperature of the ice-covered conductor; T f is the freezing temperature of water droplets; T is the ambient temperature, K, L f , L e C is the latent heat of ice melting and the latent heat of ice sublimation and evaporation respectively; a , C d , C i are the specific heats of air, water and ice respectively; h c is the convective heat transfer coefficient; U e is the boundary layer velocity; V c The water drop collision speed; A is the microelement length; r is the recovery coefficient of kinetic heating, expressed as: Where Pr is the Prandtl number, n0 = 1 / 2 in the laminar region; n0 = 1 / 3 in the turbulent region; V c is the wind speed; m im Expressed as: m im =β·LUCK·V c A β represents the local collision coefficient of the microelement; LWC is the liquid water content, and the freezing coefficient f is: Where, L ew is the Lewis number, which takes the value of 1; R v is the water vapor volume constant, which is 461 J / (Kg·k); Ф is the relative humidity of the atmosphere, which is usually 100% in an icing environment; P sat (Ts) and P sat (T) are the surface temperature of the conductor T s The saturated water vapor pressure corresponding to the ambient temperature T; the unknown parameters are m in , f and Ts are solved by iterative method.

6. A calculation method for predicting the growth of ice covering on a conductor according to claim 5, characterized in that: The method of solving by iteration includes: initially assume that the conductor surface temperature Ts is 0. If the calculation result is 0 < f < 1, then Ts = 273.15K; if the calculation result is f < 0, then let f = 0 and then solve for Ts inversely; if the calculation result is f > 1, then let f = 1 and then solve for Ts inversely; the numerical calculation starts from the stagnation point and calculates towards the upper and lower ends respectively.

7. The calculation method for predicting the growth of ice covering on a conductor according to claim 1, characterized in that: When calculating the ice thickness corresponding to each microelement and then obtaining the coordinates of the new ice shape boundary along the ice growth direction, the least squares method is used to directionally fit the discrete points of the new boundary to obtain the boundary curve function; the boundary curve expression will be used to update the boundary conditions for the iterative calculation of the next time step; the value of the time step is adjusted according to the actual ice coverage and ice growth rate; specifically, when the ice coverage increases by about 3.4%, it is considered that an iterative calculation is completed, and the time step is determined by combining the ice growth rate under this condition. The ice thickness corresponding to the ice coverage of 3.4% is 0.034D, and combined with the ice coverage rate Determine the optimal time step.

Images