Finite element analysis method of flow field of iced transmission line based on SA turbulence model

Abstract

The present invention discloses a finite element analysis method for flow field of ice-covered transmission lines based on the S-A turbulence model, which specifically includes the following steps: step one, determining the physical model and working conditions of ice-covered transmission lines; step two, establishing dimensionless boundary conditions and initial conditions of ice-covered transmission lines; step three, constructing a finite element grid for flow field of ice-covered transmission lines; step four, turbulent flow field control equations; step five, constructing a finite element numerical algorithm for solving the control equations, including time discretization, space discretization and numerical solution processes; step six, analyzing the numerical calculation results of ice-covered transmission lines. By transforming the coordinates along the streamline, a closed RANS equation based on the S-A model without convection terms is obtained, so that the RANS equation has an optimal solution when the Galerkin method is used for spatial discretization.

Description

Technical Field

[0001] The invention relates to the technical field of drilling slips, and in particular to a finite element analysis method of flow field of ice-covered transmission lines based on an SA turbulence model. Background Art

[0002] With the continuous development of high-voltage and ultra-high-voltage transmission lines in my country, the flow-induced effects caused by wind loads after ice coating have attracted more and more attention and research from many scholars. When the fluid flows around the surface of the structure, the fluid will give the structure an excitation, thereby inducing the vibration of the structure. The vibration forms of ice-covered transmission lines are mainly vortex-induced vibration and dancing. Vortex-induced vibration is manifested as the alternating shedding of vortices in the wake area after the wind flows around the transmission line, causing periodic vibration of the transmission line. When the aerodynamic coefficient formed by the wind flowing around the transmission line changes with the wind attack angle, triggering the Denharto dancing mechanism, the ice-covered transmission line will produce low-frequency, large-amplitude self-excited vibration, that is, dancing. Whether it is vortex-induced vibration or dancing, it is very important to study the flow morphology and aerodynamic coefficient of wind-covered ice transmission lines. Using wind tunnel tests to measure aerodynamic coefficients is one of the effective means, but the test cost is high and the study of flow field morphology is not detailed enough. In numerical algorithms, numerical simulation of flow based on the Reynolds time-averaged equation simulation method (RANS method) has become a research hotspot for many scholars.

[0003] Compared with the two-equation model of kw and k-ε, the RANS method of the one-equation Spalart-Allmaras (SA) model is more efficient and has been successfully applied to the simulation of flow around triangular prisms and square prisms. At present, the RANS equations are mainly solved by the characteristic line splitting method (CBS method) and the upwind finite element method (SUPG method), which are mainly first-order time discretization formats, resulting in certain differences between the numerical solutions and the experimental results, and the lift peak jump phenomenon (lift increases with the increase of wind attack angle, and suddenly decreases sharply at a certain critical wind attack angle) obtained in the calculation results of existing literature is not obvious. Summary of the invention

[0004] The object of the present invention is to provide a finite element analysis method for the flow field of ice-covered transmission lines based on the SA turbulence model and a method of using the same, so as to solve the problems raised in the above-mentioned background technology.

[0005] To achieve the above object, the present invention provides the following technical solution: a finite element analysis method for flow field of ice-covered transmission lines based on the SA turbulence model, wherein the finite element analysis method for flow field of ice-covered transmission lines specifically comprises the following steps:

[0006] Step 1: determine the physical model and operating conditions of ice-covered transmission lines;

[0007] Step 2: Establish dimensionless boundary conditions and initial conditions of ice-covered transmission lines;

[0008] Step 3: construct the finite element mesh of the flow field of ice-covered transmission lines;

[0009] Step 4: turbulent flow field control equation;

[0010] Step 5: construct a finite element numerical algorithm for solving the control equation, which includes time discretization, space discretization and numerical solution process;

[0011] Step six: Analysis of numerical calculation results of iced transmission lines.

[0012] Preferably, in the step 1, the conditions for establishing the physical model of the ice-covered transmission line are that the transmission line is of the LGJ-240 model with a diameter of 26.8 mm, and a crescent-shaped ice-covered model is adopted. Considering three ice-covered thickness conditions, the transmission line diameter is used as the characteristic length, and the model is dimensionless. The dimensionless dimensions are as follows: the transmission line diameter D=1, the small circle diameter d1=0.2633, and the ice-covered thickness H are 0.25, 0.5 and 1.0 respectively, wherein the physical model of H=0.5 is as follows Figure 1 As shown, the air density is 1.25kg / m 3 , the dynamic viscosity coefficient is 1.72×10 -5 pa.s, wind speed is 10m / s, corresponding Reynolds number Re=20000, wind direction is expressed by wind attack angle, numerical calculation is performed from 0° to 180° (Δα=10°), and encrypted calculation is performed where the local aerodynamic coefficient changes dramatically with wind attack angle.

[0013] Preferably, in step 2, dimensionless boundary conditions and initial conditions of the ice-covered transmission line are established, and the two-dimensional uniform incoming wind is considered. After dimensionless conversion, the inlet is the velocity inlet boundary condition, that is, the horizontal velocity u x is 1, the vertical velocity v y 0, auxiliary scalar =3; the coupling effect between wind and ice-covered transmission line is not considered, that is, the transmission line is fixed, which is a no-slip boundary condition, and the corresponding velocity and auxiliary scalar are all 0; the outlet is a pressure outlet, that is, the pressure p is 0, and the horizontal gradient of the other items is 0; the upper and lower boundaries are slip boundary conditions, that is, the vertical velocity is 0, and the horizontal velocity and the auxiliary scalar gradient along the vertical direction are 0, such as Figure 2 As shown, in the figure, x and y represent the horizontal and vertical directions respectively.

[0014] Preferably, in step 3, a finite element grid of the flow field of the ice-covered transmission line is constructed, the flow field area size is Ω[-16,36]×[-16,16], and linear quadrilateral units are used to mesh the space, such as Figure 3 shown.

[0015] Preferably, the turbulent flow field control equation in step 4 adopts an unsteady incompressible flow dimensionless RANS turbulence control equation closed by a one-equation SA model scalar equation;

[0016]

[0017] Among them, the subscripts i and j represent the subscripts of the space tensor, and the two-dimensional space coordinates are 1 and 2, such as x1 represents the horizontal x direction, x2 represents the vertical y direction, (x and y) represent the two-dimensional plane coordinates; u1 represents the horizontal speed u x , u2 represents the vertical velocity v y ; t represents time; p represents pressure; represents the turbulence auxiliary scalar, ν represents the kinematic viscosity, and ν t represents the eddy viscosity coefficient, which can be obtained by the auxiliary scalar field get:

[0018]

[0019] The parameter f v1 , c b1 , c b2 , σ,c w1 , f w is the process quantity in the SA model;

[0020] Introduce the following dimensionless expression:

[0021]

[0022] In the formula p * ,t * , and represent the dimensionless coordinates, velocity, pressure, time, turbulence auxiliary scalar and eddy viscosity coefficient, respectively, U ∞ , ρ, and D represent the flow field inlet velocity, density, and characteristic size (i.e., transmission line diameter), respectively;

[0023] For the convenience of expression, the “*” sign on the dimensionless expression is omitted, and the dimensionless control equation can be obtained:

[0024]

[0025] The expression of the process parameter is: g=r+c w2 (r 6 -r), Where d represents the minimum distance from the grid node to the solid wall, Ω ijrepresents vorticity; constant parameter: c b1 =0.1355,σ=2 / 3,c b2 =0.622,κ=0.41,c w1 =3.2391, c w2 =0.3, c w3 =2, c v1 =7.1;

[0026] Since equation (3) contains nonlinear convection terms, the optimal solution cannot be obtained when the standard Galerkin finite element method is used for spatial discretization. Therefore, the control equation without convection terms is obtained after coordinate transformation along the streamline:

[0027]

[0028] In the formula, x i ′ and s′ represent the coordinates and arc coordinates in the moving coordinate system respectively.

[0029] Preferably, the time in step 5 is discrete;

[0030] Now, by introducing the intermediate auxiliary speed, the momentum equation in equation (4) is split into the following two equations along the characteristic line:

[0031]

[0032] Where Δu i represents the speed correction term, The speed correction term Δu is ignored in the above formula. i Advanced small items

[0033] In order to improve the calculation accuracy, the third-order Runge-Kutta method along the uniform streamline is introduced to time discretize equation (5). Assuming that t is known n The velocity and pressure of the flow field at the given time step Δt = t n+1 -t n , t n+1 The intermediate auxiliary speed at the moment is:

[0034]

[0035] in, Indicates t n+1 The intermediate auxiliary speed at the moment, s represents t n+1 The position of the particle in the flow field at the moment, ds represents the position of the particle from t n Time to t n+1 The length of the trajectory that the flow field particle follows along the streamline at the moment, Indicates t n The speed of the moment; and Indicates tn Time and t n+1 / 2 The velocity slope at time, Indicates t n+1 The velocity slope at position s at time. The upper right superscript in the expression represents time, and the lower right subscript represents spatial position, reflecting the time and space attributes. n+1 The position at the time is defined as the static coordinate, and the other times are the moving coordinates. For the convenience of expression, the superscripts and subscripts on the right side of all variables under the static coordinates are omitted.

[0036] Introducing the definition of uniform streamline, from t n Time to t n+1 At the moment, assuming that the particle's speed is constant, the corresponding trajectory length ds along the streamline and the relationship between speed and time are as follows:

[0037]

[0038] In the formula, U n , o(Δt 2 ) represent t n Time to t n+1 Average speed at time t n The speed and second-order small quantity of the moment,

[0039] In order to overcome the difficulty of grid updating caused by coordinate transformation, the Taylor expansion along the uniform streamline is adopted to transform equations (7) to (9) into t n+1 Expression on the moment grid:

[0040]

[0041] Among them, t n+1 / 2 Prediction speed at the moment Using Euler method, we can get:

[0042]

[0043] In the formula, x k Respectively represent t n The speed and coordinates at the moment, the subscript k represents the space tensor subscript, the value is 1 and 2, 1 represents the horizontal direction, 2 represents the vertical direction, the definition is the same as the space tensor subscript i and j.

[0044] Substituting equations (13) to (15) into equation (7), we can obtain the intermediate auxiliary speed:

[0045]

[0046] Based on equations (6) and (17), as well as the coordinate transformation relationship, we can get tn+1 Speed ​​of the moment:

[0047]

[0048] Take the divergence on both sides of equation (18) and consider t n+1 The incompressible condition at time t, the pressure Poisson equation is as follows:

[0049]

[0050] The turbulent eddy viscosity coefficient in equations (16) and (17) is obtained by solving equation (4) to obtain the auxiliary scalar field And we can get that, Equation (4) can be discretized in time using the Euler method along the streamline:

[0051]

[0052] Preferably, the space discretization in step 5 is:

[0053] The Galerkin method is used to spatially discretize equations (16) to (20), and the variable u in the unit i , p and The same order interpolation functions are used to express them as follows:

[0054]

[0055] in, and They represent the velocity, pressure and turbulence scalar column vector of the unit node respectively; ψ represents the unit interpolation function,

[0056] The weight function adopts the unit interpolation function in equation (21). The unit interpolation function is multiplied by the left and right sides of equations (16), (17) and (20), and spatially integrated in the computational domain. The weak integral semi-implicit format is as follows:

[0057]

[0058] The expression of the unit matrix is:

[0059]

[0060] Similarly, the weak integral form after weighted integration of equations (19) and (20) is:

[0061]

[0062] The expression of the unit matrix is:

[0063] Preferably, the numerical solution process in step 5 is summarized as follows:

[0064] (1) Using the momentum equation (22), we can obtain t n+1 / 2 The speed of time

[0065] (2) Based on the Runge-Kutta method, equation (23), we can obtain t n+1 Intermediate auxiliary speed

[0066] (3) Using the pressure Poisson equation (26), we can obtain t n+1 The pressure of the moment n+1 ;

[0067] (4) By correcting equation (25) based on speed, we can obtain t n+1 The speed of time

[0068] (5) Based on the turbulence scalar equation (24), we can obtain t n+1 The turbulent scalar field at time

[0069] (6) Passing Calculate the turbulent eddy viscosity coefficient v t ;

[0070] (7) Return to the first step and continue the calculation at the next moment;

[0071] The relative error of the speed calculated at the two moments before and after is small enough (set to 10 in this paper). -5 ), stop the calculation, and finally obtain the flow field information around the ice-covered transmission line, including the velocity field and pressure field.

[0072] Preferably, in step 6, the analysis of the numerical calculation results of the ice-covered transmission line includes, through the finite element numerical algorithm of the present invention, finally obtaining the aerodynamic force coefficients acting on the ice-covered transmission line after being excited by the wind, including the drag coefficient C d , lift coefficient C l and torque coefficient C m ,

[0073]

[0074] Among them, F d 、F l and M represent the drag, lift and torque on the ice-covered transmission line respectively. The aerodynamic coefficients are obtained by integrating the pressure and viscous stress on the surface of the transmission line.

[0075] The numerical algorithm of the present invention can be further applied to solving the flow around a single ice-covered transmission line under different ice-covered shapes and different wind speeds, as well as solving the flow around a split line.

[0076] Compared with the prior art, the present invention has the following beneficial effects:

[0077] 1. Through the coordinate transformation along the streamline, the closed RANS equation based on the SA model without convection terms is obtained, so that the RANS equation has the optimal solution when the Galerkin method is used for spatial discretization.

[0078] 2. An intermediate auxiliary velocity term is introduced to split the closed RANS equation based on the SA model without convection terms, and the third-order Runge-Kutta method is used to time discretize it to improve the accuracy and stability of the numerical solution.

[0079] 3. The Taylor expansion along the uniform streamline is used to overcome the grid update problem caused by coordinate transformation. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] Figure 1 The crescent icing model and wind attack angle definition are shown, where Figure 1 (a) is the crescent icing model. Figure 1 (b) is defined as the wind attack angle;

[0081] Figure 2 Show the boundary conditions of flow around ice-covered transmission lines;

[0082] Figure 3 A grid diagram of flow around ice-covered transmission lines is shown;

[0083] Figure 4 The comparison of the aerodynamic coefficient curve with wind angle of attack is shown in Figure 2. Figure 4 (a) to Figure 4 (c) shows the variation of the average drag coefficient, lift coefficient and drag coefficient of a 0.25D crescent-shaped ice-covered transmission line with the wind attack angle; Figure 4 (d) to Figure 4 (f) shows the variation of the average drag coefficient, lift coefficient and drag coefficient of a 0.5D crescent-shaped ice-covered transmission line with the wind attack angle; Figure 4 (g) to Figure 4 (i) The average drag coefficient, lift coefficient and drag coefficient of a 1.0D crescent-shaped ice-covered transmission line are shown as a function of wind attack angle;

[0084] Figure 5 The aerodynamic coefficient time history curves for different ice thicknesses are shown in Figure 5 (a1) to Figure 5(i3) shows the aerodynamic coefficient time history curves corresponding to the three types of crescent-shaped ice-covered transmission lines of 0.25D, 0.5D and 1.0D when the wind attack angles are 0°, 20°, 30°, 60°, 90°, 150°, 160°, 170° and 180° respectively.

[0085] Figure 6 =20° shows the transient streamline diagram at different times in a cycle, where Figure 6 (a) to Figure 6 (e) shows the transient streamline diagram at different times within a cycle when t=t0, t=t0+1 / 4T, t=t0+1 / 2T, t=t0+3 / 4T, and t=t0+T respectively. DETAILED DESCRIPTION

[0086] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0087] See also Figure 1 - Figure 6 , an embodiment provided by the present invention: a finite element analysis method for flow field of ice-covered transmission lines based on SA turbulence model,

[0088] In the first embodiment, the finite element analysis method for flow field of ice-covered transmission lines specifically comprises the following steps:

[0089] Step 1: determine the physical model and operating conditions of ice-covered transmission lines;

[0090] Step 2: Establish dimensionless boundary conditions and initial conditions of ice-covered transmission lines;

[0091] Step 3: construct the finite element mesh of the flow field of ice-covered transmission lines;

[0092] Step 4: turbulent flow field control equation;

[0093] Step 5: construct a finite element numerical algorithm for solving the control equation, which includes time discretization, space discretization and numerical solution process;

[0094] Step six: Analysis of numerical calculation results of iced transmission lines.

[0095] In the step 1, the conditions for establishing the physical model of the ice-covered transmission line are that the transmission line is LGJ-240 model with a diameter of 26.8 mm, and a crescent-shaped ice-covered model is adopted. Considering three ice-covered thickness conditions, the transmission line diameter is used as the characteristic length, and the model is dimensionless. The dimensionless dimensions are as follows: the transmission line diameter D=1, the small circle diameter d1=0.2633, and the ice-covered thickness H are 0.25, 0.5 and 1.0 respectively, wherein the physical model of H=0.5 is as follows Figure 1 As shown, the air density is 1.25kg / m 3 , the dynamic viscosity coefficient is 1.72×10 -5 pa.s, wind speed is 10m / s, corresponding Reynolds number Re=20000, wind direction is expressed by wind attack angle, numerical calculation is performed from 0° to 180° (Δα=10°), and encrypted calculation is performed where the local aerodynamic coefficient changes dramatically with wind attack angle.

[0096] In step 2, the dimensionless boundary conditions and initial conditions of the ice-covered transmission line are established, and the two-dimensional uniform incoming wind is considered. After dimensionless conversion, the inlet is the velocity inlet boundary condition, that is, the horizontal velocity u x is 1, the vertical velocity v y 0, auxiliary scalar =3; the coupling effect between wind and ice-covered transmission line is not considered, that is, the transmission line is fixed, which is a no-slip boundary condition, and the corresponding velocity and auxiliary scalar are all 0; the outlet is a pressure outlet, that is, the pressure p is 0, and the horizontal gradient of the other items is 0; the upper and lower boundaries are slip boundary conditions, that is, the vertical velocity is 0, and the horizontal velocity and the auxiliary scalar gradient along the vertical direction are 0, such as Figure 2 As shown, in the figure, x and y represent the horizontal and vertical directions respectively.

[0097] In step 3, a finite element grid of the flow field of ice-covered transmission lines is constructed, the flow field area size is Ω[-16,36]×[-16,16], and linear quadrilateral units are used to divide the space into grids, such as Figure 3 shown.

[0098] The turbulent flow field control equation in step 4 adopts the unsteady incompressible flow dimensionless RANS turbulence control equation closed by the SA model scalar equation;

[0099]

[0100] Among them, the subscripts i and j represent the subscripts of the space tensor. When it is a two-dimensional space coordinate, the value is 1, 2, such as x1 represents the horizontal x direction, x2 represents the vertical y direction, (x and y) represents the two-dimensional plane coordinates; u1 represents the horizontal speed u x , u2 represents the vertical velocity v y ; t represents time; p represents pressure; represents the turbulence auxiliary scalar, ν represents the kinematic viscosity, and ν t represents the eddy viscosity coefficient, which can be obtained by the auxiliary scalar field get:

[0101]

[0102] The parameter f v1 , c b1 , c b2 , σ,c w1 , f w is the process quantity in the SA model;

[0103] Introduce the following dimensionless expression:

[0104]

[0105] In the formula p * ,t * , and represent the dimensionless coordinates, velocity, pressure, time, turbulence auxiliary scalar and eddy viscosity coefficient, respectively, U ∞ , ρ, and D represent the flow field inlet velocity, density, and characteristic size (i.e., transmission line diameter), respectively;

[0106] For the convenience of expression, the “*” sign on the dimensionless expression is omitted, and the dimensionless control equation can be obtained:

[0107]

[0108] The expression of the process parameter is:

[0109] Where d represents the minimum distance from the grid node to the solid wall, Ω ij represents vorticity; constant parameter: c b1 =0.1355,σ=2 / 3,c b2 =0.622,κ=0.41,c w1 =3.2391, c w2 =0.3, c w3 =2, c v1 =7.1;

[0110] Since equation (3) contains nonlinear convection terms, the optimal solution cannot be obtained when the standard Galerkin finite element method is used for spatial discretization. Therefore, the control equation without convection terms is obtained after coordinate transformation along the streamline:

[0111]

[0112] In the formula, x i ′ and s′ represent the coordinates and arc coordinates in the moving coordinate system respectively.

[0113] In the step 5, time is discrete;

[0114] Now, by introducing the intermediate auxiliary speed, the momentum equation in equation (4) is split into the following two equations along the characteristic line:

[0115]

[0116] Where Δu i represents the speed correction term, The speed correction term Δu is ignored in the above formula. i Advanced small items

[0117] In order to improve the calculation accuracy, the third-order Runge-Kutta method along the uniform streamline is introduced to time discretize equation (5). Assuming that t is known n The velocity and pressure of the flow field at the given time step Δt = t n+1 -t n , t n+1 The intermediate auxiliary speed at the moment is:

[0118]

[0119] in, Indicates t n+1 The intermediate auxiliary speed at the moment, s represents t n+1 The position of the particle in the flow field at the moment, ds represents the position of the particle from t n Time to t n+1 The length of the trajectory that the flow field particle follows along the streamline at the moment, Indicates t n The speed of the moment; and Indicates t n Time and t n+1 / 2 The velocity slope at time, Indicates t n+1 The velocity slope at position s at time. The upper right superscript in the expression represents time, and the lower right subscript represents spatial position, reflecting the time and space attributes. n+1 The position at the time is defined as the static coordinate, and the other times are the moving coordinates. For the convenience of expression, the superscripts and subscripts on the right side of all variables under the static coordinates are omitted.

[0120] Introducing the definition of uniform streamline, from t n Time to t n+1 At the moment, assuming that the particle's speed is constant, the corresponding trajectory length ds along the streamline and the relationship between speed and time are as follows:

[0121]

[0122]

[0123] In the formula, U n , o(Δt 2 ) represent t n Time to t n+1 Average speed at time t n The speed and second-order small quantity of the moment,

[0124] In order to overcome the difficulty of grid updating caused by coordinate transformation, the Taylor expansion along the uniform streamline is adopted to transform equations (7) to (9) into t n+1 Expression on the moment grid:

[0125]

[0126] Among them, t n+1 / 2 Prediction speed at the moment Using Euler method, we can get:

[0127]

[0128] In the formula, x k Respectively represent t n The speed and coordinates at the moment, the subscript k represents the space tensor subscript, the value is 1 and 2, 1 represents the horizontal direction, 2 represents the vertical direction, the definition is the same as the space tensor subscript i and j.

[0129] Substituting equations (13) to (15) into equation (7), we can obtain the intermediate auxiliary speed:

[0130]

[0131] Based on equations (6) and (17), as well as the coordinate transformation relationship, we can get t n+1 Speed ​​of the moment:

[0132]

[0133] Take the divergence on both sides of equation (18) and consider t n+1 The incompressible condition at time t, the pressure Poisson equation is as follows:

[0134]

[0135] The turbulent eddy viscosity coefficient in equations (16) and (17) is obtained by solving equation (4) to obtain the auxiliary scalar field And we can get that, Equation (4) can be discretized in time using the Euler method along the streamline:

[0136]

[0137] In step 5, the space is discretized:

[0138] The Galerkin method is used to spatially discretize equations (16) to (20), and the variable u in the unit i , p and The same order interpolation functions are used to express them as follows:

[0139]

[0140] in, and They represent the velocity, pressure and turbulence scalar column vector of the unit node respectively; ψ represents the unit interpolation function,

[0141] The weight function adopts the unit interpolation function in equation (21). The unit interpolation function is multiplied by the left and right sides of equations (16), (17) and (20), and spatially integrated in the computational domain. The weak integral semi-implicit format is as follows:

[0142]

[0143] The expression of the unit matrix is:

[0144]

[0145] Similarly, the weak integral form after weighted integration of equations (19) and (20) is:

[0146]

[0147]

[0148] The expression of the unit matrix is:

[0149] The numerical solution process in step 5 is summarized as follows:

[0150] (1) Using the momentum equation (22), we can obtain t n+1 / 2 The speed of time

[0151] (2) Based on the Runge-Kutta method, equation (23), we can obtain t n+1 Intermediate auxiliary speed

[0152] (3) Using the pressure Poisson equation (26), we can obtain t n+1 The pressure of the moment n+1 ;

[0153] (4) By correcting equation (25) based on speed, we can obtain t n+1 The speed of time

[0154] (5) Based on the turbulence scalar equation (24), we can obtain t n+1 The turbulent scalar field at time

[0155] (6) Passing Calculate the turbulent eddy viscosity coefficient v t ;

[0156] (7) Return to the first step and continue the calculation at the next moment;

[0157] The relative error of the speed calculated at the two moments before and after is small enough (set to 10 in this paper). -5 ), stop the calculation, and finally obtain the flow field information around the ice-covered transmission line, including the velocity field and pressure field.

[0158] In step 6, the analysis of the numerical calculation results of the ice-covered transmission line includes, through the finite element numerical algorithm of the present invention, finally obtaining the aerodynamic force coefficients acting on the ice-covered transmission line after being excited by the wind, including the drag coefficient C d , lift coefficient C l and torque coefficient C m ,

[0159]

[0160] Among them, F d 、F l and M represent the drag, lift and torque on the ice-covered transmission line respectively. The aerodynamic coefficients are obtained by integrating the pressure and viscous stress on the surface of the transmission line.

[0161] Figure 4The graphs are the relationship between the average drag coefficient, average lift coefficient and average torque coefficient of the transmission line under three ice thicknesses and the wind angle of attack, and the results are compared and analyzed with the existing wind tunnel test and numerical simulation results. From the results of the measured aerodynamic coefficient, the calculated results are basically consistent with the overall trend of the test and calculation results in the existing literature, and a lift coefficient spike jump phenomenon is presented, indicating that the SA model and algorithm adopted in the present invention can be used to numerically calculate the aerodynamic characteristics of the transmission line with different ice thicknesses. From a quantitative comparison, under the same wind angle of attack, when the wind angle of attack is between 0° and 40° and between 140° and 250°, the aerodynamic characteristics of the transmission line with different ice thicknesses can be numerically calculated. Between 180° and 180°, the calculated aerodynamic coefficients of the two ice thicknesses are basically consistent with the experimental results of Shimizu, and are better than the calculation results of Shinichi. It can be seen that within the small angle of attack range, better quantitative calculation results can be obtained. By comparing the calculation results of different ice thicknesses with the experimental results, the calculation results of ice thickness 1.0D are better than those of 0.5D, and the calculation results of 0.5D are better than those of 0.25D, indicating that the SA turbulence model has better accuracy in simulating the flow around the airfoil at small angles of attack.

[0162] Figure 5 The figure shows the comparison of the time history curves of the drag, lift and torque coefficients of the transmission line with different ice thickness under each typical wind attack angle. It can be seen from the figure that after the flow field around the ice-covered transmission line has evolved for a period of time, the drag coefficient, lift coefficient and torque coefficient acting on the ice-covered transmission line all show periodic changes, which further illustrates the periodic evolution of the wake of the ice-covered transmission line. At the same time, the amplitude fluctuation of the three aerodynamic coefficients can be obtained, which is helpful to study the vortex-induced vibration of the ice-covered transmission line.

[0163] Using the numerical calculation results, the streamline diagram under each working condition can be further obtained. Figure 5 When the wind attack angle is 20°, the upper tangent line of the ice-covered part is exactly in the horizontal direction. Figure 6 The transient streamline diagrams at five different moments in one cycle of this working condition are given. At this time, the boundary layer separation occurs when the ice is applied to the highest and lowest points in the vertical diameter direction of the transmission line from the end of the boundary layer flow. The vortices on the upper and lower sides of the wake area are generated, moved and disappeared periodically, making the aerodynamic coefficient fluctuate periodically, which will induce vortex-induced vibration of the transmission line. The symmetry of the wake vortex is destroyed relative to the angle of attack of 0°, making the lift period no longer twice the drag period.

[0164] The numerical algorithm of the present invention can be further applied to solving the flow around a single ice-covered transmission line under different ice-covered shapes and different wind speeds, as well as solving the flow around a split line.

[0165] It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above and that the invention can be implemented in other specific forms without departing from the spirit or essential features of the invention. Therefore, the embodiments should be considered exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description, and it is intended that all variations falling within the meaning and scope of the equivalent elements of the claims be included in the invention. Any reference numeral in a claim should not be considered as limiting the claim to which it relates.

Claims

1. Finite element analysis method of ice-covered transmission line flow field based on SA turbulence model, characterized in that: The finite element analysis method for flow field of ice-covered transmission lines specifically comprises the following steps: Step 1: determine the physical model and operating conditions of ice-covered transmission lines; Step 2: Establish dimensionless boundary conditions and initial conditions of ice-covered transmission lines; Step 3: construct the finite element mesh of the flow field of ice-covered transmission lines; Step 4: turbulent flow field control equation; Step 5: construct a finite element numerical algorithm for solving the control equation, which includes time discretization, space discretization and numerical solution process; Step 6: Analysis of numerical calculation results of ice-covered transmission lines; The turbulent flow field control equation in step 4 adopts the unsteady incompressible flow dimensionless RANS turbulence control equation closed by the SA model scalar equation; Among them, the subscripts i and j represent the spatial tensor subscripts, two-dimensional spatial coordinates, and the values ​​are 1 and 2, such as x1 represents the horizontal x direction, x2 represents the vertical y direction, and x and y represent the two-dimensional plane coordinates; u1 represents the horizontal speed u x , u2 represents the vertical velocity v y ; t represents time; p represents pressure; represents the turbulence auxiliary scalar, ν represents the kinematic viscosity, and ν t represents the eddy viscosity coefficient, which can be obtained by the auxiliary scalar field get: The parameter f v1 , c b1 , c b2 , σ,c w1 , f w is the process quantity in the SA model; Introduce the following dimensionless expression: In the formula p * ,t * , and represent the dimensionless coordinates, velocity, pressure, time, turbulence auxiliary scalar and eddy viscosity coefficient, respectively, U ∞ , ρ, and D represent the flow field inlet velocity, density, and characteristic size, i.e., the transmission line diameter, respectively; For the convenience of expression, the "*" sign on the dimensionless expression is omitted, and the dimensionless control equation can be obtained: The expression of the process parameter is: Where d represents the minimum distance from the grid node to the solid wall, Ω ij represents vorticity; constant parameter: c b1 =0.1355,σ=2 / 3,c b2 =0.622,κ=0.41,c w1 =3.2391, c w2 =0.3, c w3 =2, c v1 =7.1; Since equation (3) contains nonlinear convection terms, the optimal solution cannot be obtained when the standard Galerkin finite element method is used for spatial discretization. Therefore, the control equation without convection terms is obtained after coordinate transformation along the streamline: In the formula, x i ′ and s′ represent the coordinates and arc coordinates in the moving coordinate system respectively; In the step 5, time is discrete; Now, by introducing the intermediate auxiliary speed, the momentum equation in equation (4) is split into the following two equations along the characteristic line: Where Δu i represents the speed correction term, The speed correction term Δu is ignored in the above formula. i Advanced small items In order to improve the calculation accuracy, the third-order Runge-Kutta method along the uniform streamline is introduced to time discretize equation (5). Assuming that t is known n The velocity and pressure of the flow field at the given time step Δt = t n+1 -t n , t n+1 The intermediate auxiliary speed at the moment is: in, Indicates t n+1 The intermediate auxiliary speed at the moment, s represents t n+1 The position of the particle in the flow field at the moment, ds represents the position of the particle from t n Time to t n+1 The length of the trajectory that the flow field particle follows along the streamline at the moment, Indicates t n The speed of the moment; and Indicates t n Time and t n+1 / 2 The velocity slope at time, Indicates t n+1 The velocity slope at position s at time. The upper right superscript in the expression represents time, and the lower right subscript represents spatial position, reflecting the time and space attributes. n+1 The position at the time is defined as the static coordinate, and the other times are the moving coordinates. For the convenience of expression, the superscripts and subscripts on the right side of all variables under the static coordinates are omitted. Introducing the definition of uniform streamline, from t n Time to t n+1 At the moment, assuming that the particle's speed is constant, the corresponding trajectory length ds along the streamline and the relationship between speed and time are as follows: In the formula, U n 、o(Δt 2 ) represent t n Time to t n+1 Average speed at time t n The speed and second-order small quantity of the moment, In order to overcome the difficulty of grid updating caused by coordinate transformation, the Taylor expansion along the uniform streamline is adopted to transform equations (7) to (9) into t n+1 Expression on the moment grid: Among them, t n+1 / 2 Prediction speed at the moment Using Euler method, we can get: In the formula, x k Respectively represent t n The speed and coordinates at the moment, the subscript k represents the space tensor subscript, and the values ​​are 1 and 2, 1 represents the horizontal direction, and 2 represents the vertical direction, which is the same as the definition of the space tensor subscripts i and j; Substituting equations (13) to (15) into equation (7), we can obtain the intermediate auxiliary speed: Based on equations (6) and (17), as well as the coordinate transformation relationship, we can get t n+1 Speed ​​of the moment: Take the divergence on both sides of equation (18) and consider t n+1 The incompressible condition at time t, the pressure Poisson equation is as follows: The turbulent eddy viscosity coefficient in equations (16) and (17) is obtained by solving equation (4) to obtain the auxiliary scalar field ν. Equation (4) is time discretized using the Euler method along the streamline to obtain:

2. The finite element analysis method for flow field of ice-covered transmission lines based on SA turbulence model according to claim 1 is characterized in that: In the step 1, the conditions for establishing the physical model of the ice-covered transmission line are that the transmission line is of the LGJ-240 model with a diameter of 26.8 mm, and a crescent-shaped ice-covered model is adopted. Considering three ice-covered thickness conditions, the transmission line diameter is used as the characteristic length, and the model is dimensionless. The dimensionless dimensions are as follows: the transmission line diameter D=1, the small circle diameter d1=0.2633, and the ice-covered thickness H are 0.25, 0.5 and 1.0 respectively, wherein the air density of the physical model with H=0.5 is 1.25 kg / m 3 , the dynamic viscosity coefficient is 1.72×10 -5 pa.s, wind speed is 10m / s, corresponding Reynolds number Re=20000, wind direction is expressed by wind attack angle, numerical calculation is performed from 0° to 180°, and encrypted calculation is performed where the local aerodynamic coefficient changes dramatically with wind attack angle.

3. The finite element analysis method for flow field of ice-covered transmission lines based on SA turbulence model according to claim 1 is characterized in that: In step 2, the dimensionless boundary conditions and initial conditions of the ice-covered transmission line are established, and the two-dimensional uniform incoming wind is considered. After dimensionless conversion, the inlet is the velocity inlet boundary condition, that is, the horizontal velocity u x is 1, the vertical velocity v y 0, auxiliary scalar is 3; the coupling effect between wind and ice-covered transmission lines is not considered, that is, the transmission lines are fixed, which is a no-slip boundary condition, and the corresponding velocities and auxiliary scalars are all 0; the outlet is a pressure outlet, that is, the pressure p is 0, and the horizontal gradients of the other items are 0; the upper and lower boundaries are slip boundary conditions, that is, the vertical velocity is 0, and the horizontal velocity and the gradient of the auxiliary scalar along the vertical direction are 0, and x and y represent the horizontal and vertical directions, respectively.

4. The finite element analysis method for flow field of ice-covered transmission lines based on SA turbulence model according to claim 1 is characterized in that: In the step three, a finite element grid of the flow field of the ice-covered transmission line is constructed, the flow field area size is Ω[-16,36]×[-16,16], and linear quadrilateral units are used to divide the space into grids.

5. The finite element analysis method for flow field of ice-covered transmission lines based on SA turbulence model according to claim 1 is characterized in that: In step 5, the spatial discretization is performed by using the Galerkin method to discretize equations (16) to (20). The variable u in the unit i , p and The same order interpolation functions are used to express them as follows: in, and They represent the velocity, pressure and turbulence scalar column vector of the unit node respectively; ψ represents the unit interpolation function, The weight function adopts the unit interpolation function in equation (21). The unit interpolation function is multiplied by the left and right sides of equations (16), (17) and (20), and spatially integrated in the computational domain. The weak integral semi-implicit format is as follows: The expression of the unit matrix is: Similarly, the weak integral form after weighted integration of equations (19) and (20) is: The expression of the unit matrix is:

6. The finite element analysis method for flow field of ice-covered transmission lines based on SA turbulence model according to claim 1 is characterized in that: The numerical solution process in step 5 is summarized as follows: (1) Using momentum equation (22), we can obtain t n+1 / 2 The speed of time (2) Based on the Runge-Kutta method, equation (23), we can obtain t n+1 Intermediate auxiliary speed (3) Using the pressure Poisson equation (26), we can obtain t n+1 The pressure of the moment n+1 ; (4) By correcting equation (25) based on speed, we can obtain t n+1 The speed of time (5) Based on the turbulence scalar equation (24), we can obtain t n+1 The turbulent scalar field at time (6) Passing Calculate the turbulent eddy viscosity coefficient v t ; (7) Return to the first step and continue the calculation at the next moment; The relative error of the speed calculated at the two moments before and after is small enough, so it is set to 10 -5 , stop the calculation, and finally obtain the flow field information around the ice-covered transmission line, including velocity field and pressure field.

7. The finite element analysis method for flow field of ice-covered transmission lines based on SA turbulence model according to claim 1 is characterized in that: In step 6, the analysis of the numerical calculation results of the ice-covered transmission line includes, through the finite element numerical algorithm, finally obtaining the aerodynamic force coefficient acting on the ice-covered transmission line after being excited by the wind, including the drag coefficient C d , lift coefficient C l and torque coefficient C m , where F d 、F l and M represent the drag, lift and torque on the ice-covered transmission line respectively. The aerodynamic coefficients are obtained by integrating the pressure and viscous stress on the surface of the transmission line.

Images