CN114528781A - CFD numerical simulation method combining actual measurement wind-induced snow drift - Google Patents

Abstract

A CFD numerical simulation method combining actual measurement wind-induced snow drift relates to a CFD numerical simulation method of wind-induced snow drift. The method aims to solve the technical problems that relative motion between wind and snow is not considered and accuracy is poor in the existing numerical simulation method. The method comprises the following steps: modeling based on an initial snow distribution form measured by field actual measurement and actually measured snow particle data, considering the relative slip speed of wind phase and snow phase, adopting programming software Visual Studio to write UDF and introduce finite element numerical simulation analysis software Ansys Fluent to simulate the uneven snow distribution of the structure, and utilizing software Tecplot to carry out post-processing on simulated data, thereby establishing a CFD numerical simulation method which combines the actually measured wind-induced snow drift considering the relative movement of wind and snow. The invention improves the accuracy of the numerical simulation method, and can be used in the field of building structure design in cold regions.

Description

A CFD Numerical Simulation Method Combined with Measured Wind-induced Snow Drift

technical field

The invention belongs to the field of civil engineering and relates to a CFD numerical simulation method of wind-induced snow drift.

Background technique

In recent years, blizzards in cold areas have brought great inconvenience to people's travel, communication, production and life, and also greatly threatened the safety of houses and public buildings. Among them, wind-induced snow drifts can easily lead to excessive loads at local locations of buildings. However, uneven snow loads on such buildings tend to have additional effects on the structure, making it more susceptible to instability. Therefore, the research and prevention of wind and snow has become more and more important.

In the early research on wind and snow disasters, limited by the computer conditions, scholars at home and abroad focused on actual measurement and experimental research. With the continuous progress of society and the rapid development of science and technology, numerical simulation based on computational fluid dynamics Research is gradually being adopted by the general public. The current numerical simulation study of wind-induced snow drift adopts the assumption of one-way coupling, and the parameter setting adopts the empirical value. In practice, the deposition and erosion will affect the snow boundary conditions, thereby further affecting the flow field distribution, and this method is essentially the calculation of the snow phase diffusion in a stable wind field, which is very different from the actual situation and reduces the the accuracy of the simulation.

SUMMARY OF THE INVENTION

The present invention aims to solve the technical problems that the relative motion between wind and snow is not considered and the accuracy is poor in the existing numerical simulation method, and provides a CFD numerical simulation method combined with the measured wind-induced snow drift. The method considers the relative motion of wind and snow, and optimizes a series of processes from the modeling of wind and snow two-phase flow, to parameter setting, and simulation calculation, and obtains the CFD simulation value of wind-induced snow drift.

The CFD numerical simulation method of the present invention combined with the measured wind-induced snow drift is carried out according to the following steps:

Step 1: Collect the basic information required for numerical simulation: collect the geometric characteristics of the building under study, the geographic location of the building, the topographic parameters, the average wind speed V ₀ at the reference height H ₀ , the wind profile specification, the turbulent kinetic energy specification, Turbulence intensity specification, integrated turbulence scale specification, turbulent energy dissipation rate specification; meanwhile, the snow particle density, snow particle diameter d _s , angle of repose, wind speed, wind direction and the initial snow depth of the building are measured. , calculate the threshold friction velocity, the snow particle velocity of jumping, the ground roughness, the building wall roughness, the snow surface roughness and the snow particle settling velocity w _f ;

Step 2: Establish a grid model required for numerical simulation: use the building geometric information and the initial snow distribution pattern collected in step 1 to establish a model and divide the grid to obtain a grid file;

Step 3: Write the air term boundary conditions for numerical simulation: use the geographic location information of the building location collected in Step 1, the average wind speed at the reference height, the wind profile, the turbulent kinetic energy, the turbulent intensity, the integrated turbulent scale, and the turbulent dissipation rate. The boundary conditions of the air phase are written in the normative information, and then based on the C language and the Fluent secondary development interface in Visual Studio, the exponential wind profile, turbulent kinetic energy, and turbulent dissipation rate programs are written, and the C language file is saved. Step 2: Grid model files are placed in the same folder;

Step 4: Write the snow boundary conditions for numerical simulation: use the snow particle density collected in step 1, the gravitational acceleration at the location of the building, and the threshold friction speed, and use the empirical formula of snow concentration and the snow movement speed used in the formula. , Snow particle friction velocity, jump and critical height information of the suspension layer, write the snow phase volume fraction formula based on C language and Fluent secondary development interface in Visual Studio, as the inlet boundary condition of the snow term, save the C language file , and put it in the same folder as the model grid file in step 2 and the C language file in step 3;

Step 5: Program access required for numerical simulation. The specific operations are as follows: Open the Ansys-Fluent interface and select 2D/3D, double precision, and parallelism. Under General Options-Working Directory, select the grid file containing step 2 and the C in step 3. Language file, the folder directory of the C language file in step 4, select the Set upCompilation Environment for UDF interface under the environment to open the software; select UserDefined–Functions–Compiled under the Ansys-Fluent main menu, and select the C language file in steps 3 and 4. Load and compile to get the inlet velocity and volume fraction of the air and snow phases;

Step 6: Using Ansys-Fluent software, first set the model, material, and solver parameters, and then import the grid model obtained in step 2, the inlet velocity and volume fraction of the air phase and snow phase obtained in step 5, and carry out two-phase wind and snow. Flow simulation to obtain convergence files for wall shear forces, air and snow velocities and volume fractions;

Step 7: Use the convergence file calculated in Step 6 to perform post-processing of the wall shear force and velocity: Save the convergence file calculated by Ansys-Fluent in the format supported by Tecplot, use Tecplot for post-processing calculation, and calculate by extracting the calculation The components of the velocity, pressure, and wall shear force in the x, y, and z directions in the domain three-dimensional space. Use the Tecplot custom function Specify Equations to make the cloud map and equivalent value of the pressure, velocity, and wall shear force within the calculation domain. line graph;

把壁面剪切力转换为壁面的摩擦速度，其中τ为壁面剪切力，ρ为空气密度；对壁面摩擦速度与阈值磨擦速度进行比较，当壁面摩擦速度大于阈值磨擦速度时，积雪发生侵蚀；当壁面摩擦速度小于阈值磨擦速度时，积雪发生沉积；分别计算侵蚀量、沉积量，并且计算积雪改变厚度，最后将积雪改变厚度与建筑物的初始积雪深度进行叠加，得到积雪的最终分布形态，完成结合实测的风致雪漂移的CFD数值模拟。Step 8: Post-processing calculation of numerical simulation of snow deposition erosion: Based on the wall shear force calculated in Step 7, extract the wall shear force on the bottom surface of the model or the horizontal and vertical midline of the bottom surface through Tecplot, and use the formula

Convert the wall shear force to the wall friction velocity, where τ is the wall shear force, ρ is the air density; compare the wall friction velocity with the threshold friction velocity, when the wall friction velocity is greater than the threshold friction velocity, the snow will erode ; When the friction velocity of the wall is less than the threshold friction velocity, the snow is deposited; the erosion amount and deposition amount are calculated respectively, and the snow change thickness is calculated. Finally, the snow change thickness and the building's initial snow depth are superimposed to obtain the The final distribution of snow, complete the CFD numerical simulation combined with the measured wind-induced snow drift.

Further, the solution process of the grid file described in the step 2 is: in Ansys-Icem, the model calculation domain is determined according to the geometric characteristics of the building and the calculation domain setting satisfies that the blocking rate is less than 3%, and then establishes according to the model in Icem. For the mesh model, in order to ensure the smooth convergence of the calculation and the accuracy of the simulation, set a reasonable mesh size according to the size of the model, control the mesh skewness (Skewness) of triangles and tetrahedra not greater than 0.95, and control the mesh aspect ratio (Aspect Ratio). , control the mesh squashing degree (Squish) between 0-1, and finally output the solved mesh file in Ansys-Fluent format.

Further, the reasonable grid size refers to: the grid model controls the structured grid quality to be 0.6-1, and the unstructured grid quality is 0.3-1; the grid aspect ratio is: the grid in the flow core area. The aspect ratio is less than 5:1, and the aspect ratio of the boundary layer mesh is less than 10:1.

Further, the snow concentration empirical formula described in step 4 is the snow concentration empirical formula proposed by Pomeroy & Gray in 1990 edition and Pomeroy & Male in 1992 edition.

Further, the simulation of the two-phase flow of wind and snow described in step 6, the specific process is as follows:

1) Using Ansys-Fluent software, select Transient Transient in the Time column under General;

2) Select the Mixture model in the Multiphase column under Models, select the Realizable k-epsilon model under the k-epsilon two-equation model in the Viscous column under Models, and select Phase in the Multiphase column under Models Interaction-slip, set the sliding speed algebraic equation;

3) Select Fluid under Material Materials, add air phase, set Viscosity viscosity, Density density and create a new material snow phase, set snow phase Viscosity viscosity to be consistent with air and set Density density and particle size according to step 1;

4) Under the volume grid Cell Zone Conditions, set the required fluid fluid and solid solid in the grid computational domain, and set the boundary conditions in the Boundary Conditions: the inlet boundary is the velocity inlet, the two sides and top of the computational domain are symmetrical boundaries, The ground and model surfaces are non-slip walls, and the exit boundary is a free exit;

5) In the Boundary Conditions-inlet (inlet), access the air phase boundary conditions of step 3: set the wind profile, turbulent kinetic energy and dissipation rate, and access the Volume Fraction in the Boundary Conditions-inlet (inlet)-snow snow phase-Multiphase The volume fraction of the snow phase in step 4; the Fluent-Mixture model is calculated by the finite volume method. Since there are only two phases of the wind and snow flow, and the sum of the integral fractions of the two phases is 1, only the volume fraction of the snow phase needs to be given, and the air phase is given by Software calculation without setting;

6) In the selection method Methods under the solver setting Solution, select the SIMPLE algorithm;

7) In the Select Monitors-Residual under the Solver setting Solution, set the Continuity Equation, Momentum Equation and the Convergence Accuracy of k-epsilon, where the Convergence Accuracy of k-epsilon is set to 1e- ⁶ Accuracy;

8) In the initialization Initialization under the solver setting Solution, select the standard initialization Standard Initialization, and initialize at the entrance;

9) Finally, set the calculation time step and the number of calculation steps in Run Calculation, start the calculation, and obtain the convergence files of the wall shear force, the velocity and volume fraction of the air phase and the snow phase.

其中A_ero为常系数，u_*为入口摩擦速度，u_*t为阈值速度。Further, in step 8, the calculation formula of the erosion amount is:

where A _ero is the constant coefficient, u _* is the inlet friction velocity, and u _*t is the threshold velocity.

其中w_f中为雪的沉降速度，φ_s为雪的质量浓度，φ_s＝ρ_sf，f为雪的体积分数，ρ_s为雪颗粒密度。Further, in step 8, the calculation formula of the deposition amount is:

where w _f is the snow settling velocity, φ _s is the snow mass concentration, φ _s = ρ _s f, f is the snow volume fraction, and ρ _s is the snow particle density.

其中q_s为q_ero或q_dep，γ为积雪的最大体积分数，ρ_s为雪颗粒密度。Furthermore, in step 8, the formula for calculating the thickness of snow change is:

where q _s is q _er o or q _dep , γ is the maximum volume fraction of snow cover, and ρ _s is the snow particle density.

The method of the invention uses Icem CFD software to establish a model, applies the measured cube and surrounding snow distribution in the establishment of the model, considers the initial snow distribution of wind-induced snow drift, and uses Ansys Fluent to carry out numerical simulation analysis , This method is based on the two-phase flow theory of the Euler-Euler model, combined with the numerical calculation and analysis program of the improved turbulence model using Visual Studio, and the measured wind speed, wind direction, and physical properties of snow particles (shape, density, size, etc. ), the threshold speed of wind and snow drift, the falling speed of snow particles, the movement mode of wind and snow drift, the angle of repose and other data are used in the simulation to simulate the snow distribution of the building, and complete the CFD value combined with the measured wind-induced snow drift. simulation. The simulation results were compared with the classical experiment of wind-induced snow drift, and with the method that was not combined with the actual measurement, to verify the accuracy of the simulation method. The invention considers the relative motion between wind and snow and combines it with actual measurement, thereby improving the accuracy of numerical simulation of wind-induced snow drift.

The overall flow diagram of the CFD numerical simulation method combined with the measured wind-induced snow drift of the present invention is shown in Figure 1. The numerical simulation method of the present invention combined with the measured wind-induced snow drift has the following advantages:

1. The relative slip velocity between the wind phase and the snow phase is considered, and the effect between the two phases is considered by solving the relative velocity between the two phases.

2. The snow distribution pattern of wind-blown snow obtained from the actual measurement is used in the establishment of the model, and the influence of snow deposition and erosion on the flow field distribution in the calculation domain is considered.

3. Use experimental tools to measure various data based on local conditions and apply them in numerical simulation, so that the simulation results are more in line with reality and more accurate.

The method of the present invention helps to recognize the law of wind and snow movement, helps prevent the collapse of buildings, and reduces casualties and economic losses. This method can be used in the field of building structure design in cold regions.

Description of drawings

1 is a schematic diagram of the overall flow of the CFD numerical simulation method of the present invention in conjunction with the measured wind-induced snow drift;

Fig. 2 is the grid division of embodiment 1 and the figure of computation domain;

Fig. 3 is the snow distribution form figure that obtains through step 7 of embodiment 1;

Fig. 4 shows the snow flux per unit time obtained in Example 1, the snow flux per unit time obtained by the classical experiment, and the snow flux per unit time obtained by the method not combined with actual measurement in Example 2 Quantitative comparison chart.

Detailed ways

Verify the beneficial effects of the present invention with the following examples:

Embodiment 1: This embodiment is a numerical simulation method combined with the measured wind-induced snow drift, and is carried out according to the following steps:

z₀与v₀分别表示参考点的高度(1m)及对应的参考平均风速(5m/s)，z表示任一高度；α表示地貌参数，对于B类地貌，α取值为0.16；湍动能采用公式

湍能耗散率采用公式

其中，C_μ为经验常数，取0.09；I和l分别为湍流强度和湍流长度尺度，其中：Step 1: Collect the basic information required for numerical simulation: Collect the geometric characteristics of the building to be studied: a cube model with a side length of 1m, the geographic location of the building, and the geomorphological parameters of 0.16. The average wind speed at a reference height of 10m is 5m/s , the wind profile adopts the formula

z ₀ and v ₀ represent the height of the reference point (1m) and the corresponding reference average wind speed (5m/s), respectively, z represents any height; α represents the landform parameter, for B-type landform, α is 0.16; turbulent kinetic energy using the formula

Turbulent energy dissipation rate using the formula

Among them, C _μ is an empirical constant, which is taken as 0.09; I and l are the turbulent intensity and turbulent length scales, respectively, where:

积分湍流尺度采用公式

其中，z_G、z_b、α'、均为地貌参数，在日本AIJ规范中，B类地貌对应Ⅱ类平地，以上三个地貌参数取值分别为5m、350m和0.15。The turbulent kinetic energy intensity adopts the formula

The integral turbulence scale adopts the formula

Among them, z _G , z _b , α', are all geomorphological parameters. In the Japanese AIJ code, the B-type landform corresponds to the II-level flat land, and the values of the above three geomorphic parameters are 5m, 350m and 0.15 respectively.

The local snow particle diameter is 0.00015m, the angle of repose is 45°, the snow density is 145kg/m ³ , the calculated threshold friction speed is 0.2m/s, The jumping snow particle velocity is 0.2m/s, the ground roughness coefficient is 0.25, and the snow particle settling velocity is 0.2m/s.

Step 2: Establish the grid model required for numerical simulation: In Ansys-Icem CFD, the model is established and the grid is divided according to the building geometric information and the initial snow distribution shape collected in step 1, and the model calculation is determined. The length, width and height are 20m, 10m and 8m respectively, and the blockage rate of the computational domain is set to 1.25%. Then, a grid model is established based on the model in Icem CFD, and the number of grids is 3.6 million according to the size of the model. The mesh model controls the structured mesh quality at 0.8, the unstructured mesh quality at 0.8, and finally outputs the solved mesh file in Ansys-Fluent format. The graph of mesh division and computational domain is shown in Fig. 2; in Fig. 2, 1 is the symmetric boundary, 2 is the free development boundary, 3 is the velocity inlet, and 4 is the non-slip wall.

Step 3: Write the boundary conditions of the air term for numerical simulation: use the geographic location information of the building collected in Step 1, the average wind speed at the reference height, the wind profile, the turbulent kinetic energy, the turbulent intensity, the integral turbulence scale, and the turbulent dissipation rate specification Write the boundary conditions of the air phase based on the C language and the Fluent secondary development interface in Visual Studio, write the exponential wind profile, turbulent kinetic energy, and turbulent dissipation rate, save the C language file, and compare it with the second step. Lattice model files are placed in the same folder;

其中ρ_s表示雪颗粒密度，g为重力加速度，取9.81m/s²，u_p为跃移层雪粒平均速度，u_p＝2.8u_*t，u_*为入口摩擦速度；根据Beyers等实测拟合得到的摩擦速度与10m高度处风速u₁₀的关系确定，

雪相运动速度取5m/s、雪颗粒摩擦速度取0.15m/s、u_*t为阈值速度；取0.15m/s；h_s为跃移层和悬移层的临界高度，

雪相运动速度5m/s、雪颗粒摩擦速度0.15m/s、在Visual Studio中以C语言、Fluent二次开发接口为基础编写雪相体积分数公式作为雪项的入口边界条件，保存C语言文件，并将其与步骤二模型网格文件、步骤三的C语言文件放在同一文件夹中；Step 4: Write the snow boundary conditions for numerical simulation: using the density of snow particles collected in step 1 of 145kg/m ³ , the local acceleration of gravity 9.81m/s ² , and the threshold friction velocity of 0.2m/s, Pomeroy & Gray (1990), The empirical formula for snow concentration proposed by Pomeroy & Male (1992) is the formula

where ρ _s is the density of snow particles, g is the acceleration of gravity, which is 9.81m/s ² , u _p is the average speed of snow particles in the transition layer, u _p = 2.8u _*t , u _* is the inlet friction speed; according to the actual measurement of Beyers et al. The relationship between the friction speed obtained by fitting and the wind speed u ₁₀ at a height of 10m is determined,

The speed of snow phase movement is 5m/s, the friction speed of snow particles is 0.15m/s, u _*t is the threshold speed; 0.15m/s is taken; h _s is the critical height of the transition layer and the suspension layer,

Snow phase motion speed is 5m/s, snow particle friction speed is 0.15m/s, and the snow phase volume fraction formula is written in Visual Studio based on C language and Fluent secondary development interface as the inlet boundary condition of the snow term, and the C language file is saved , and put it in the same folder as the model grid file in step 2 and the C language file in step 3;

Step 5: Program access required for numerical simulation. The specific operations are as follows: Open the Ansys-Fluent interface and select 3D, double precision, and parallel. Under General Options-Working Directory, select the folder directory containing the grid file in step 2 and the custom program in step 3 and 4, and select the Set up Compilation Environmentfor UDF interface under the environment to open the software. Select the C language files in steps 3 and 4 under User-Defined-Functions-Compiled in the Ansys-Fluent main menu to load and compile, and obtain the inlet velocity and volume fraction of the air phase and snow phase;

Step 6: Using Ansys-Fluent software, first set the model, material, and solver parameters, and then import the grid model obtained in step 2, the inlet velocity and volume fraction of the air phase and snow phase obtained in step 5, and carry out two-phase wind and snow. Flow simulation to obtain the convergence files of wall shear force, air and snow phase velocities and volume fractions, as follows:

1) Select transient Transient in the Time column under General;

2) Select the Mixture model in the Multiphase column under Models, select the Realizable k-epsilon model under the k-epsilon two-equation model in the Viscous column under Models, and select the Multiphase multiphase model under Models. Select Phase Interaction-slip and set the relative sliding speed to 0.25m/s;

3) Under the material Materials, select Fluid, add the air phase, set the Viscosity viscosity and Density density, and create a new material snow phase, set the snow phase Viscosity viscosity to be the same as air, and set the Density density to 145kg/m ³ according to step 1 (by Snow Fork Snow characteristic analyzer), particle size 0.00015m (measured by sample sieve);

4) Set the required fluid fluid and solid solid in the grid computing domain under the volume grid Cell Zone Conditions, and set the boundary conditions in the Boundary Conditions: the inlet boundary is the velocity inlet, the sides and top of the computational domain are symmetrical boundaries, the ground and the model The surface is a non-slip wall, and the outlet boundary is a free outlet;

5) In the Boundary Conditions-inlet (inlet), access the air phase boundary conditions of step 3: set the wind profile, turbulent kinetic energy and dissipation rate, and access the Volume Fraction in the Boundary Conditions-inlet (inlet)-snow snow phase-Multiphase The volume fraction of the snow phase in step 4; the Fluent-Mixture model is calculated by the finite volume method. Since there are only two phases of the wind and snow flow, and the sum of the integral fractions of the two phases is 1, only the volume fraction of the snow phase needs to be given, and the air phase is given by Software calculation requires no setup.

6) Select the SIMPLE algorithm in Methods under the solver setting Solution.

7) Under the solver setting Solution, select Monitors-Residual to set the convergence precision of the continuity equation, the momentum equation and the k-epsilon. In this embodiment, the convergence precision is set to 1e- ⁶ precision.

8) Select the standard initialization StandardInitialization in the initialization initialization under the solver setting Solution to initialize at the entrance.

9) Final calculation Set the number of calculation steps in Run Calculation to 1000 steps and the accuracy to 10 ^-6 , start the calculation, and obtain the convergence files of wall shear force, air phase and snow phase velocity and volume fraction.

Step 7: Use the convergence file calculated in Step 6 to perform post-processing of the wall shear force and velocity: Save the convergence file calculated by Ansys-Fluent in the format supported by Tecplot, use Tecplot for post-processing calculation, and calculate by extracting the calculation The components of the velocity, pressure, and wall shear force in the x, y, and z directions in the domain three-dimensional space. Use the Tecplot custom function Specify Equations to make the cloud map and equivalent value of the pressure, velocity, and wall shear force within the calculation domain. Line diagram; the distribution of snow cover is shown in Figure 3.

计算侵蚀量，其中A_ero为常系数，取7.0×10^-4；当壁面摩擦速度小于阈值磨擦速度时沉积，采用公式

计算沉积量，其中w_f为雪的沉降速度，φ_s为雪的质量浓度，φ_s＝ρ_sf，f为雪的体积分数，ρ_s为雪颗粒密度；再采用公式

计算积雪改变厚度，其中q_s为q_ero或q_dep，γ为积雪的最大体积分数，取0.62；最后将积雪改变厚度与步骤一的建筑物的初始积雪深度进行叠加，得到积雪的最终分布形态，完成结合实测的风致雪漂移的CFD数值模拟。Step 8: Post-processing calculation for numerical simulation of snow deposition erosion: Based on the wall shear force calculated in step 7: Based on the wall shear force calculated in step 7, extract the wall shear force on the bottom of the model through Tecplot, and use a custom function Convert the wall shear stress to the friction velocity of the wall, and compare the wall friction velocity with the threshold friction velocity. When the wall friction velocity is greater than the threshold friction velocity, it is erosion, using the formula

Calculate the amount of erosion, where A _ero is a constant coefficient, taking 7.0×10 ^-4 ; when the wall friction velocity is less than the threshold friction velocity, the formula is used for deposition

Calculate the amount of deposition, where w _f is the sedimentation velocity of snow, φ _s is the mass concentration of snow, φ _s = ρ _s f, f is the volume fraction of snow, and ρ _s is the density of snow particles; then use the formula

Calculate the changed thickness of snow cover, where q _s is q _ero or q _dep , γ is the maximum volume fraction of snow cover, which is taken as 0.62; finally, superimpose the changed thickness of snow cover with the initial snow cover depth of the building in step 1 to obtain the The final distribution of snow, complete the CFD numerical simulation combined with the measured wind-induced snow drift.

Embodiment 2: This embodiment simulates the numerical value of wind-induced snow drift without combining actual measurement and without considering the relative slip velocity of wind phase and snow phase. The specific method is carried out according to the following steps:

z与v₀分别表示参考点的高度(1m)及对应的参考平均风速(5m/s)；α表示地貌参数，对于B类地貌，α取值为0.16；湍动能采用公式

湍能耗散率采用公式

其中，C_μ为经验常数，取0.09；I和l分别为湍流强度和湍流长度尺度，湍动能强度采用公式

积分湍流尺度采用公式

其中，z_G、z_b、α'、均为地貌参数，在日本AIJ规范中，B类地貌对应Ⅱ类平地，以上三个地貌参数取值分别为5m、350m和0.15。Step 1: Collect the basic information required for numerical simulation: the geometric characteristics of the building under study: a cube model with a side length of 1m, the geographic location of the building where the geomorphological parameters are 0.16, the average wind speed at a reference height of 10m is 5m/s, The wind profile adopts the formula

z and v ₀ represent the height of the reference point (1m) and the corresponding reference average wind speed (5m/s), respectively; α represents the geomorphic parameter, for Type B landform, the value of α is 0.16; the turbulent kinetic energy adopts the formula

Turbulent energy dissipation rate using the formula

Among them, C _μ is an empirical constant, which is taken as 0.09; I and l are the turbulent flow intensity and turbulent flow length scale, respectively, and the turbulent kinetic energy intensity adopts the formula

The integral turbulence scale adopts the formula

Among them, z _G , z _b , α', are all geomorphological parameters. In the Japanese AIJ code, the B-type landform corresponds to the II-level flat land, and the values of the above three geomorphic parameters are 5m, 350m and 0.15 respectively.

Relevant numerical values are adopted through experience: the diameter of snow particles is 0.00015m, the angle of repose is not considered, the snow density is 150kg/m ³ , the calculated threshold friction speed is 0.2m/s, the snow particle speed is 0.2m/s, and the ground roughness coefficient is 0.25 , Snow particle settling velocity 0.2m/s.

Step 2: Establish the grid model required for numerical simulation: In Ansys-Icem CFD, the model is established and the grid is divided according to the building geometric information and the initial snow distribution shape collected in step 1, and the model calculation is determined. The length, width and height are 20m, 10m and 8m respectively, and the blockage rate of the computational domain is set to 1.25%. Then, a grid model is established based on the model in Icem CFD, and the number of grids is 3.6 million according to the size of the model. The mesh model controls the structured mesh quality at 0.8, the unstructured mesh quality at 0.8, and finally outputs the solved mesh file in Ansys-Fluent format.

Step 3: Write the boundary conditions of the air term for numerical simulation: use the geographic location information of the building collected in Step 1, the average wind speed at the reference height, the wind profile, the turbulent kinetic energy, the turbulent intensity, the integral turbulence scale, and the turbulent dissipation rate specification Write the boundary conditions of the air phase based on the existing technology in Visual Studio, based on the C language and the Fluent secondary development interface, write the exponential wind profile, turbulent kinetic energy, and turbulent dissipation rate, save the C language file, and put it Put it in the same folder as the step 2 mesh model file;

其中ρ_s表示雪颗粒密度，g为重力加速度，取9.81m/s²，u_p为跃移层雪粒平均速度，u_p＝2.8u_*t，u_*为入口摩擦速度，根据Beyers等实测拟合得到的摩擦速度与10m高度处风速u₁₀的关系确定，

雪相运动速度取5m/s、雪颗粒摩擦速度取0.15m/s、u_*t为阈值速度；取0.15m/s；h_s为跃移层和悬移层的临界高度，

雪相运动速度5m/s、雪颗粒摩擦速度0.15m/s、在Visual Studio中以C语言、Fluent二次开发接口为基础编写雪相体积分数公式作为雪项的入口边界条件，保存C语言文件，并将其与步骤二模型网格文件、步骤三的C语言文件放在同一文件夹中；Step 4: Write the snow phase boundary conditions for numerical simulation: using the snow particle density of 150kg/m ³ collected in step 1, the local gravitational acceleration 9.81m/s ² , and the threshold friction velocity of 0.2m/s, using Pomeroy & Gray (1990), The empirical formula for snow concentration proposed by Pomeroy & Male (1992) is the formula

where ρ _s is the density of snow particles, g is the acceleration of gravity, 9.81m/s ² is taken, u _p is the average speed of snow particles in the transition layer, u _p =2.8u _*t , u _* is the friction velocity at the entrance, according to the actual measurement of Beyers et al. The relationship between the friction speed obtained by fitting and the wind speed u ₁₀ at a height of 10m is determined,

The speed of snow phase movement is 5m/s, the friction speed of snow particles is 0.15m/s, u _*t is the threshold speed; 0.15m/s is taken; h _s is the critical height of the transition layer and the suspension layer,

Snow phase motion speed is 5m/s, snow particle friction speed is 0.15m/s, and the snow phase volume fraction formula is written in Visual Studio based on C language and Fluent secondary development interface as the inlet boundary condition of the snow term, and the C language file is saved , and put it in the same folder as the model grid file in step 2 and the C language file in step 3;

Step 5: Program access required for numerical simulation. The specific operations are as follows: save the grid file in step 2 and the custom programs in steps 3 and 4 in the same folder, open the Ansys-Fluent interface and select 3D, double precision, and parallel. Select the above folder directory under GeneralOptions-Working Directory, and select the Set up CompilationEnvironment for UDF interface under the environment to open the software. Select User-Defined-Functions under the Ansys-Fluent main menu

- Under Compiled, select the C language files in steps 3 and 4 to load and compile, and obtain the inlet velocity and volume fraction of the air phase and the snow phase;

Step 6: Set the model, material and solver parameters of the numerical simulation software in detail, use Ansys-Fluent to import the grid model, and simulate the two-phase flow of wind and snow based on the secondary development programming imported in Step 5, as follows:

1) Select transient Transient in the Time column under General;

2) Select the Mixture model in the Multiphase column under Models, select the Realizable k-epsilon model under the k-epsilon two-equation model in the Viscous column under Models, and select the Multiphase multiphase model under Models. When Phase Interaction-slip is selected, the relative slip velocity of the wind phase and the snow phase is not considered at this time;

3) Under the material Materials, select Fluid, add the air phase, set the Viscosity viscosity and Density density, and create a new material snow phase, set the snow phase Viscosity viscosity to be the same as air, and set the Density density to 150kg/m ³ and the particle size to 0.00015 according to step 1. m;

4) Set the required fluid fluid and solid solid in the grid computing domain under the volume grid Cell Zone Conditions, and set the boundary conditions in the Boundary Conditions: the inlet boundary is the velocity inlet, the sides and top of the computational domain are symmetrical boundaries, the ground and the model The surface is a non-slip wall, and the outlet boundary is a free outlet;

5) In the Boundary Conditions-inlet (inlet), access the air phase boundary conditions of step 3: set the wind profile, turbulent kinetic energy and dissipation rate, and access the Volume Fraction in the Boundary Conditions-inlet (inlet)-snow snow phase-Multiphase The volume fraction of the snow phase in step 4; the Fluent-Mixture model is calculated by the finite volume method. Since there are only two phases of the wind and snow flow, and the sum of the integral fractions of the two phases is 1, only the volume fraction of the snow phase needs to be given, and the air phase is given by Software calculation without setting;

6) Select the SIMPLE algorithm in Methods under the solver setting Solution.

7) Under the solver setting Solution, select Monitors-Residual to set the convergence precision of the continuity equation, the momentum equation and the k-epsilon. In this embodiment, the convergence precision is set to 1e- ⁶ precision.

8) Select the standard initialization StandardInitialization in the initialization initialization under the solver setting Solution to initialize at the entrance.

9) Final calculation Set the number of calculation steps in Run Calculation to 1000 steps, and the accuracy is 10 ^-6 , start the calculation, and obtain the convergence files of the wall shear force, the velocity and volume fraction of the air phase and the snow phase.

Step 7: Use the convergence file calculated in Step 6 to post-process the wall shear force and velocity: After the Ansys-Fluent calculation converges, the file is saved in the Tecplot supported format. Post-processing calculation with Tecplot, by extracting the calculated components of x, y, z velocity, pressure, and wall shear force, and use Tecplot's custom function Specify Equations to make pressure, velocity, and wall shear within the computational domain. Force cloud map, contour map.

其中A_ero为常系数，取7.0×10^-4；当壁面摩擦速度小于阈值磨擦速度时沉积，采用公式

其中w_f中为雪的沉降速度，φ_s为雪的质量浓度(φ_s＝ρ_sf，f为雪的体积分数，ρ_s为雪颗粒密度)计算侵蚀量、沉积量；基于步骤一的初始积雪深度，采用公式

计算积雪改变厚度、叠加积雪改变量(其中，统一用q_s代替q_ero和q_dep，γ为积雪的最大体积分数，取0.62)，得到积雪的最终分布形态，完成风致雪漂移的数值模拟。Step 8: Post-processing calculation for numerical simulation of snow deposition erosion: Based on the wall shear force calculated in Step 7, extract the wall shear force on the bottom of the model through Tecplot, and convert the wall shear stress into the friction velocity of the wall through a custom function , the wall friction velocity is compared with the threshold friction velocity. When the wall friction velocity is greater than the threshold friction velocity, the erosion is carried out using the formula

where A _ero is a constant coefficient, which is taken as 7.0×10 ^-4 ; when the wall friction velocity is less than the threshold friction velocity, the deposition is made using the formula

Wherein w _f is the sedimentation velocity of snow, φ _s is the mass concentration of snow (φ _s = ρ _s f, f is the volume fraction of snow, ρ _s is the density of snow particles) to calculate the amount of erosion and deposition; Initial snow depth, using the formula

Calculate the change thickness of snow cover and the change amount of superimposed snow cover (where q _s is used to replace q _ero and q _dep , and γ is the maximum volume fraction of snow cover, which is taken as 0.62) to obtain the final distribution form of snow cover, and complete the wind-induced snow drift numerical simulation.

The snow flux per unit time shown by the simulation results of Example 1, and the snow accumulation per unit time obtained in Example 2 without combining with the actual measurement and ignoring the relative slip velocities of the wind and snow phases The snow flux and the snow flux curve per unit time of the classic experiment are compared, as shown in Figure 4. It can be seen from Figure 4 that Example 2 is not combined with the actual measurement, and does not consider the wind phase and snow phase. The snow flux per unit time obtained under the condition of relative slip velocity is quite different from the curve obtained by the classical experimental method, and the simulation results of Example 1 are in good agreement with the results of the classical experiment. The accuracy of the simulation method of 1 is high.

Claims (8)

1. A CFD numerical simulation method combined with actually measured wind-induced snow drift is characterized by comprising the following steps of:

the method comprises the following steps: collecting basic information required for numerical simulation: collecting the geometric characteristics of the building under study, the geographical position of the building, the geomorphic parameters, the reference height H₀Average wind speed V₀Wind profile specification, turbulence energy specification, turbulence intensity specification, integral turbulence scale specification and turbulence energy dissipation rate specification; simultaneously measuring the snow particle density and the snow particle diameter d of the accumulated snow at the location of the building_sAngle of repose, wind speed, wind direction and initial snow depth of building, calculating threshold friction speed, speed of jumping snow particles, and ground roughnessRoughness of building wall surface, roughness of snow surface and sedimentation velocity w of snow particles_f；

Step two: establishing a grid model required by numerical simulation: building a model and dividing grids by using the geometric information of the building and the initial snow distribution form collected in the step one to obtain a grid file;

step three: writing air term boundary conditions for numerical simulation: compiling boundary conditions of air phase by using the geographical position information of the place where the building is located, the reference height average wind speed, the wind profile, the turbulence energy, the turbulence intensity, the integral turbulence scale and the turbulence dissipation rate specification information collected in the step one, compiling exponential type wind profile, turbulence energy and turbulence dissipation rate programs based on the C language and Fluent secondary development interfaces in Visual Studio, storing a C language file, and placing the C language file and the grid model file in the step two in the same folder;

step four: writing snow phase boundary conditions for numerical simulation: writing a snow phase volume fraction formula in Visual Studio based on a C language and Fluent secondary development interfaces as an entrance boundary condition of a snow item by using the snow particle density collected in the step one, the gravity acceleration of the place where the building is located and the threshold friction speed and the snow phase motion speed, the snow particle friction speed, the jump and the critical height information of a suspension layer adopted in a snow concentration empirical formula and a formula thereof, storing a C language file, and putting the C language file, the model grid file in the step two and the C language file in the step three in the same folder;

step five: numerically simulating the required program access, specifically operating as follows: opening an Ansys-Fluent interface to select 2D/3D, double precision and parallelism, selecting a folder Directory containing a grid file in the second step, a C language file in the third step and a C language file in the fourth step under a General Options-Working Directory, and selecting Set up compatibility Environment for UDF interface opening software under the Environment; selecting User Defined-Functions-completed under the Ansys-Fluent main menu, and selecting the C language files in the third step and the fourth step for loading and compiling to obtain the inlet speeds and volume fractions of the air phase and the snow phase;

step six: adopting Ansys-Fluent software, firstly setting model, material and solver parameters, then introducing the grid model obtained in the step two and the inlet speed and volume fraction of the air phase and the snow phase obtained in the step five, and carrying out wind-snow two-phase flow simulation to obtain a convergence file of wall surface shearing force, the speed and volume fraction of the air phase and the snow phase;

step seven: and (5) performing post-processing on the wall surface shearing force and the wall surface shearing force by using the convergence file calculated in the step six: storing the convergence file calculated by Ansys-Fluent in a Tecplot support format, performing post-processing calculation by using the Tecplot, and making a cloud map and a contour map of pressure, speed and wall shear force in a calculation domain range by extracting components of speed, pressure and wall shear force in x, y and z directions in a three-dimensional space of a calculation domain and using a Tecplot custom function specificity Equations;

step eight: post-processing calculation of numerically simulated snow deposit erosion: extracting the wall shear force on the bottom surface of the model or the horizontal and vertical center lines of the bottom surface through Tecplot based on the wall shear force obtained by the calculation in the step seven, and obtaining the wall shear force through a formula

Converting the wall surface shearing force into the wall surface friction speed, wherein tau is the wall surface shearing force, and rho is the air density; comparing the wall surface friction speed with a threshold friction speed, and when the wall surface friction speed is greater than the threshold friction speed, corroding the accumulated snow; when the wall surface friction speed is lower than the threshold friction speed, the accumulated snow is deposited; and respectively calculating erosion amount and deposition amount, calculating the snow accumulation change thickness, finally superposing the snow accumulation change amount and the initial snow accumulation depth of the building to obtain the final distribution form of the snow accumulation, and completing the CFD numerical simulation combining the actually measured wind-induced snow drift.

2. A CFD numerical simulation method in combination with measured wind-induced snow drift according to claim 1, characterized in that the solving step of the grid file in step two is: determining a model calculation domain according to the geometric characteristics of a building in the Ansys-Ican, setting the calculation domain to meet the requirement that the blocking rate is less than 3%, then establishing a mesh model in the Ican according to the model, setting reasonable mesh size according to the size of the model, regulating and controlling the mesh skewness to be that a triangle and a tetrahedron are not more than 0.95, regulating and controlling the aspect ratio of the mesh, regulating and controlling the mesh flattening degree to be between 0 and 1, and finally outputting a solved mesh file in an Ansys-Fluent format.

3. The method according to claim 2, wherein the reasonable grid size is: the grid model controls the quality of the structured grid to be 0.6-1, and the quality of the unstructured grid to be 0.3-1.

4. A CFD numerical simulation method incorporating measured wind induced snow drift according to claim 2, characterised in that the grid aspect ratio is: the grid aspect ratio in the flow core area is less than 5:1, and the boundary layer grid aspect ratio is less than 10: 1.

5. A CFD numerical simulation method combining measured wind-induced snow drift according to claim 1, 2, 3 or 4, characterized in that the wind-snow two-phase flow simulation in the sixth step is as follows:

1) adopting Ansys-Fluent software to select Transient in the column of Time under General;

2) selecting a texture model in one column of the multi-stream multi-Phase under the model Models, selecting a Rearizable k-epsilon model in a k-epsilon two-equation model in one column of the Viscous under the model Models, selecting a Phase Interaction-slip in one column of the multi-stream multi-Phase under the model Models, and setting a slip speed algebraic equation;

3) selecting Fluid under material Materials, adding air phase, setting Viscosity Viscosity and Density Density, creating material snow phase, setting the Viscosity Viscosity to be consistent with air phase, and setting Density Density and particle size according to the first step;

4) setting the required fluid and Solid solution in the grid calculation domain under the volume grid Cell Zone Conditions, and setting the Boundary Conditions in Boundary Conditions: the entrance boundary is a speed entrance, the two sides and the top of the calculation domain are symmetrical boundaries, the ground and the model surface are non-slip wall surfaces, and the exit boundary is a free exit;

5) air phase Boundary Conditions at Boundary Conditions-inlet access step three: setting a wind section, turbulent kinetic energy and dissipation rate, and accessing the Volume Fraction of the snow phase in the fourth step under Boundary Conditions-inlets-snow phase-Multiphase; the Fluent-texture model is calculated by adopting a finite volume method, and as the wind and snow flows only have two phases and the integral fraction sum of the two phases is 1, the volume fraction of the snow phase is only required to be given, and the air phase is calculated by software without setting;

6) selecting a SIMPLE algorithm in a selection method under the Solution set by a solver;

7) in the selection monitor-resolution under the solver setting Solution, the continuity equation, the momentum equation and the convergence accuracy of k-epsilon are set, wherein the convergence accuracy of k-epsilon is set to 1e^-6Precision;

8) in the Initialization under the Solution set by the solver, the Standard Initialization is selected and initialized at the entrance;

9) and finally, setting the calculated time step length and the calculated step number in Run Calculation, and starting Calculation to obtain a convergence file of the wall shear force, the air phase and snow phase speed and the volume fraction.

6. A CFD numerical simulation method in combination with measured wind-induced snow drift according to claim 1, 2, 3 or 4, characterized in that in step eight, the erosion amount is calculated by the formula:

wherein A is_eroIs a constant coefficient of u_*As inlet friction speed, u_*tIs the threshold speed.

7. A CFD numerical simulation method in combination with measured wind-induced snow drift according to claim 1, 2, 3 or 4, characterised in thatIn step eight, the calculation formula of the deposition amount is as follows:

wherein w_fThe speed of sedimentation of the snow in the middle, phi_sIs the mass concentration of snow, phi_s＝ρ_sf, f is the volume fraction of snow, ρ_sIs the snow particle density.

8. A CFD numerical simulation method in combination with measured wind-induced snow drift according to claim 1, 2, 3 or 4, characterized in that in step eight, the calculation formula of the snow change thickness is:

wherein q is_sIs q is_eroOr q_depγ is the maximum volume fraction of accumulated snow, ρ_sIs the snow particle density.

Images