A method for simulating dust dispersion in an open pit coal mine

Abstract

This invention discloses a method for simulating dust diffusion in open-pit coal mines, belonging to the field of dust diffusion simulation in mining areas. The method includes the following steps: S1, deploying a multi-scale meteorological dust collaborative monitoring network in the mining area and acquiring a standardized time-series monitoring dataset; S2, generating hourly three-dimensional meteorological parameter fields for the entire mining area with horizontal and vertical stratification; S3, constructing a full-element dynamic CFD boundary condition system covering the flow field and thermal field; S4, conducting transient gas-solid two-phase dust diffusion coupled simulation, outputting hourly particle size dust concentration fields and turbulent kinetic energy fields for the entire mining area; S5, using an RNG k-ε turbulence model coupled with a discrete phase model to conduct transient gas-solid two-phase dust diffusion coupled simulation, outputting hourly particle size dust concentration fields and turbulent kinetic energy fields for the entire mining area; S6, updating and iterating. Using the above-mentioned open-pit coal mine dust diffusion simulation method, the hourly stratified spatiotemporal distribution of dust in the mining area under day and night conditions, mining dynamics, water spraying for dust suppression, and multi-process coupling can be accurately reproduced.

Description

Technical Field

[0001] This invention relates to the field of dust diffusion simulation technology in mining areas, and more particularly to a method for simulating dust diffusion in open-pit coal mines. Background Technology

[0002] Open-pit coal mines continuously generate large amounts of dust throughout the entire process of mining, transportation, blasting, and waste disposal, which can easily lead to pneumoconiosis among workers and excessive air pollution in the surrounding area. Simultaneously, due to the large elevation differences and fragmented terrain of the mine, diurnal temperature variations create thermal convection / nighttime inversion layers, resulting in dramatic spatiotemporal changes in wind and temperature fields and highly complex dust diffusion patterns. Traditional empirical estimations relying on fixed dust generation coefficients to simply calculate the total amount of dust in a region cannot meet the needs of refined dust suppression design and environmental impact assessment. Summary of the Invention

[0003] The purpose of this invention is to provide a method for simulating dust diffusion in open-pit coal mines, thereby solving the aforementioned technical problems.

[0004] To achieve the above objectives, the present invention provides a method for simulating dust diffusion in open-pit coal mines, comprising the following steps: S1. Deploy a multi-scale meteorological dust collaborative monitoring network in the mining area and obtain a standardized time-series monitoring dataset: Set up a 10-meter-high high-level benchmark meteorological station at the top of the highest spoil heap upstream of the prevailing wind direction, and deploy ground meteorological monitoring points and dust monitoring points within the mining area and within a set range of the mining area to collect meteorological parameters of all elements and dust concentration data of different particle sizes. After outlier removal, consistency verification and standardization processing, a standardized time-series monitoring dataset is generated. S2. Based on the standardized time-series monitoring dataset output in step S1, and combined with the three types of terrain feature parameters extracted from the mine digital elevation model, namely elevation difference, slope aspect, and terrain shading, a three-dimensional inverse distance weighted interpolation model with terrain correction factor is constructed. Vertical constraints are then applied using the atmospheric boundary layer wind power law profile to generate a three-dimensional meteorological parameter field of the entire mining area with horizontal and vertical layers for each hour. S3. Decompose the hourly three-dimensional meteorological parameter field of the entire mining area output in step S2 into three types of boundary parameters: inlet dynamic wind profile, wall grid net radiative heat flux, and outlet static pressure. Automatically load the boundary parameters hourly through user-defined functions to construct a full-element dynamic CFD boundary condition system covering the flow field and thermal field. S4. Based on the monthly mining progress plan, the mine geometric model is updated synchronously to build a surface dust source parameter library covering the entire process of drilling, blasting, mining and loading, road transportation, crushing and disposal. The initial dust concentration field is generated by combining the dual dust source Gaussian diffusion model with specular reflection correction and gravity settlement correction. S5. Driven by the full-element dynamic CFD boundary conditions output in step S3, load the dynamic mine geometric model, surface dust source parameters and initial dust concentration field output in step S4. Couple the discrete phase model with the RNGk-ε turbulence model with Busineske thermal buoyancy correction to carry out transient gas-solid two-phase dust diffusion coupling simulation of fine particulate matter, inhalable particulate matter and total suspended particulate matter of multi-particle size classification, and output the hourly full mine particle size dust concentration field and turbulent kinetic energy field. S6. Use the in-situ monitoring data output in step S1 to update and iterate steps S1-S5.

[0005] Therefore, the beneficial effects of the above-mentioned open-pit coal mine dust diffusion simulation method are as follows: 1. Significantly improved meteorological field interpolation accuracy: By superimposing three correction factors of elevation difference, slope aspect angle and terrain shading, combined with wind power law vertical temperature and wind correction, it adapts to the complex steps and slope terrain of open-pit mines and outputs a three-dimensional layered meteorological field with a resolution of 1 hour, solving the problems of wind, temperature and air pressure interpolation distortion and inversion layer inability to be depicted in mountainous mining areas. 2. Achieve automatic updating of dynamic time-series CFD boundaries for all elements: The hourly inlet wind profile, net radiative heat flux of the wall, and stratified outlet static pressure are uniformly loaded through UDF custom functions, distinguishing the top surface boundary for sunny and cloudy days, without the need for manual interruption to modify the boundary, and fully coupling the dynamic meteorological drive of the flow field and thermal field. 3. Adapt to monthly dynamic mining conditions and synchronously iterate the mesh geometry: The mine pit DEM and solid model are updated with the mining progress. The structured mesh combined with the local densification scheme is adopted to meet the RNGk-ε near-wall turbulence calculation accuracy. Unlike the traditional fixed static mine pit model, it can truly reflect the flow field changes caused by the advancement of the steps and the rise of the spoil heap. 4. Construct a multi-process dynamic dust source database with full coverage and water spraying attenuation: uniformly quantify the source strength of the six major processes of drilling, mining, blasting, transportation, crushing and soil disposal, set three-stage dust suppression attenuation coefficients according to surface moisture content, distinguish between lignite and bituminous coal dust density, and adopt the Rossin-Lammle particle size distribution to fill the defects of traditional simulation that only has a single dust source, ignores water spraying dust reduction, and simplifies particle size. 5. Dual-source Gaussian initial field accelerates transient simulation convergence: It integrates the correction of physical dust gravity settling and the superposition of virtual dust sources of step mirror reflection to calculate the three-dimensional dust initial concentration field, replacing the zero-iteration mode without initial field, shortening the transient simulation calculation time, and restoring the dust reflection and accumulation effect of the mine wall.

[0006] 6. Complete particle force model for accurate dust diffusion tracking by particle size: The DPM discrete model incorporates resistance, gravity buoyancy, Saffman shear lift, and pressure gradient force to simulate the settling, diffusion, and wall capture of multi-sized dust particles such as PM2.5, PM5, PM10, and TSP in stages, achieving bidirectional coupling between gas and solid phases and improving the accuracy of dust concentration prediction at different scales.

[0007] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0008] Figure 1 This is a flowchart of a dust diffusion simulation method for open-pit coal mines according to the present invention. Detailed Implementation

[0009] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings and examples. It should be understood that the specific embodiments described herein are merely illustrative of the embodiments of the present invention and are not intended to limit the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of this application. Examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout.

[0010] It should be noted that the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion, such as a process, method, system, product, or server that includes a series of steps or units, not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such process, method, product, or device.

[0011] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0012] like Figure 1 As shown, a method for simulating dust diffusion in open-pit coal mines includes the following steps: S1. Deploy a multi-scale meteorological dust collaborative monitoring network in the mining area and obtain a standardized time-series monitoring dataset: Set up a 10-meter-high high-level benchmark meteorological station at the top of the highest spoil heap upstream of the prevailing wind direction, and deploy ground meteorological monitoring points and dust monitoring points within the mining area and within a set range of the mining area. Collect meteorological parameters of all elements and dust concentration data of different particle sizes. After outlier removal, consistency verification and standardization processing, generate a standardized time-series monitoring dataset.

[0013] S2. Based on the standardized time-series monitoring dataset output in step S1, and combined with the three types of terrain feature parameters extracted from the mine digital elevation model, namely elevation difference, slope aspect, and terrain shading, a three-dimensional inverse distance weighted interpolation model with terrain correction factor is constructed. Vertical constraints are then applied using the atmospheric boundary layer wind power law profile to generate a horizontally and vertically layered hourly three-dimensional meteorological parameter field for the entire mining area.

[0014] S3. Decompose the hourly three-dimensional meteorological parameter field of the entire mining area output in step S2 into three types of boundary parameters: inlet dynamic wind profile, wall grid net radiative heat flux, and outlet static pressure. Automatically load the boundary parameters hourly through user-defined functions to construct a full-element dynamic CFD boundary condition system covering the flow field and thermal field.

[0015] S4. Based on the monthly mining progress plan, the mine geometric model is updated synchronously to build a surface dust source parameter library covering the entire process of drilling, blasting, mining and loading, road transportation, crushing and disposal. The initial dust concentration field is generated by combining the dual dust source Gaussian diffusion model with specular reflection correction and gravity settlement correction.

[0016] S5. Driven by the full-element dynamic CFD boundary conditions output in step S3, load the dynamic mine geometric model, surface dust source parameters and initial dust concentration field output in step S4. Couple the discrete phase model with the RNGk-ε turbulence model with Busineske thermal buoyancy correction to carry out transient gas-solid two-phase dust diffusion coupling simulation of fine particulate matter, inhalable particulate matter and total suspended particulate matter of multi-particle size classification. Output the hourly full mine particle size dust concentration field and turbulent kinetic energy field.

[0017] S6. Use the in-situ monitoring data output in step S1 to update and iterate steps S1-S5.

[0018] The full range of meteorological parameters mentioned in step S1 include wind speed, wind direction, temperature, humidity, and air pressure.

[0019] Step S2 specifically includes the following steps: S21. Extraction of topographic feature parameters in the mining area: Based on the digital elevation model of the mine pit and the coordinates of each monitoring point output in step S1, and the standardized time-series monitoring data, extract the multi-dimensional topographic correlation parameters between the interpolation grid and the monitoring points.

[0020] S211. Traverse each interpolation grid point within the interpolation grid, calculate its horizontal distance, elevation difference, and slope parameters relative to all monitoring points, and form a point-by-point terrain parameter matrix. : ; in, ; ; ; In the formula, For the first Horizontal distance between each monitoring point and the interpolation grid point; For interpolation grid points and the first The absolute value of the elevation difference of each monitoring point; The maximum vertical height difference of the interpolation grid within the monitoring range of the monitoring point; , These represent the planar spans of the slope along the X and Y axes, respectively. The coordinates of the interpolation grid points; For the first Coordinates of each monitoring point; The total number of monitoring points participating in the interpolation; For the first The angle between the slope aspect of the monitoring point and the slope aspect of the interpolation grid point.

[0021] S212. Calculate the slope aspect of the slope where each monitoring point is located based on the digital elevation model of the mine pit: ; in, ; ; In the formula, For the first The slope direction of each monitoring point; , These are the topographic elevation gradients along the X and Y axes, respectively. , , , These are the elevations of the interpolated grids adjacent to the east, west, north, and south of the monitoring point grid, respectively. The edge length of the interpolation grid is the planar side length.

[0022] S213. Calculate the terrain obstruction coefficient: ; in, ; ; In the formula, For the first Terrain occlusion correction factor for each monitoring point; For the first The maximum occlusion height difference between each monitoring point and the interpolation grid point; For the first The horizontal straight line segment between each monitoring point and the interpolation grid point The actual terrain elevation; For the first The horizontal distance between each monitoring point and the interpolation grid point is [number] monitoring points. Theoretical spatial elevation at the location; The distance from the monitoring point to the location to be determined on the horizontal straight line segment between the monitoring point and the interpolation grid point is the horizontal length.

[0023] S214. Statistical dominant wind direction: Based on the full-element meteorological standardized time series data output in step S1, the dominant wind direction at a height of 10 meters within the interpolation period is statistically analyzed and used as the airflow direction benchmark for slope aspect weight correction.

[0024] S22. Based on the horizontal distance parameters extracted in step S211, calculate the basic interpolation weights of each monitoring point to the interpolation grid points: ; In the formula, For the first The basic interpolation weights of each monitoring point to the interpolation grid points; The distance power exponent is set to a fixed value of 2.

[0025] S23. Introduce elevation difference correction factor, slope aspect angle correction factor, and terrain shading correction factor to perform point-by-point correction on the basic interpolation weights obtained in step S22, and obtain the comprehensive interpolation weights: ; in, ; ; In the formula, For the first The comprehensive correction weight of each monitoring point to the interpolation grid points; For the first Elevation difference correction factor corresponding to each monitoring point; For the first Slope aspect angle correction factor corresponding to each monitoring point; For reference attenuation height, a fixed value of 100m is used; The dominant wind direction during the interpolation period.

[0026] S24. The full-element meteorological parameter values ​​of all monitoring points are weighted and summed according to the comprehensive correction weight to obtain the horizontal interpolated meteorological parameters of the interpolation grid points: ; In the formula, For the horizontal interpolated meteorological parameter values ​​of the interpolation grid points; For the first Measured meteorological parameter values ​​at each monitoring point.

[0027] S25. Based on the horizontal interpolation results obtained in step S24, the vertical wind speed distribution is corrected using the atmospheric boundary layer wind power law profile, and the vertical temperature distribution is corrected by combining seasonal and diurnal inversion characteristics, thereby improving the interpolation accuracy of meteorological parameters in the vertical direction.

[0028] S251. Vertical wind speed profile correction: The wind speed at different altitudes is constrained using the power law of the atmospheric boundary layer to ensure that the vertical wind speed distribution conforms to the laws of the near-surface atmospheric boundary layer. ; In the formula, Vertical height of interpolation grid points Corrected wind speed at the location; For reference height Horizontal interpolated wind speed at the location; For reference height, a monitoring height of 10m from the high-level benchmark meteorological station was used; The wind profile power exponent, and .

[0029] S252. Vertical Temperature Gradient Correction: The vertical temperature gradient is set according to the season and day / night scene: a positive decreasing gradient is used during the day in summer, with the temperature decreasing by 0.6℃ for every 100m increase; a negative inverse temperature gradient is used at night in winter, with the temperature increasing by 0.3℃ to 0.5℃ for every 100m increase, correcting the temperature distribution across the entire altitude layer.

[0030] S26. Generation of spatiotemporal dynamic three-dimensional meteorological parameter field: Based on the horizontal interpolation results of step S24 and the vertical correction results of step S25, a time-by-time three-dimensional meteorological parameter field covering the entire CFD computational domain is generated at a fixed time step.

[0031] The full CFD computational domain is a three-dimensional closed solution space defined by the gas-solid two-phase coupled simulation of open-pit coal mine dust. It is based on the geometric center of the mine pit, with a length and width of 5 times the east-west length of the mine pit. The entrance boundary is 1.5 times the mine length from the west edge of the mine pit, and the exit boundary is 2.5 times the mine length from the east edge of the mine pit. The total width of the computational domain is 5 times the north-south width of the mine pit, with 2 times the mine width reserved on each side. The total height of the computational domain is 4 times the maximum height difference of the mine pit. The full CFD computational domain completely includes the entire production area of ​​the mining area, spoil heap, transportation roads, and crushing station, and buffer spaces are reserved upstream and downstream, on both sides, and above the mine pit to avoid simulation errors caused by boundary cut-off of airflow, inversion layer, and dust diffusion flow field.

[0032] In step S26, hourly interpolation is first performed with a time step of 1 hour to generate a three-dimensional meteorological parameter sequence for 24 times throughout the day.

[0033] Then, a three-dimensional regular grid covering the entire CFD computation domain is constructed, with a horizontal resolution of 50m and a vertical resolution of 10m. Each grid point stores five types of meteorological parameters: wind speed, wind direction, temperature, humidity, and air pressure.

[0034] Finally, the spatial coordinate reference and timestamp are unified to form a complete hourly three-dimensional meteorological parameter field.

[0035] Step S3 specifically includes the following steps: S31. Based on the hourly three-dimensional meteorological parameter field output in step S2, generate a three-dimensional dynamic entrance wind profile with turbulence parameters. The entrance cross-section data is a time-series meteorological dataset extracted from all grid points of the vertical entrance elevation on the windward side of the full CFD computational domain from the hourly three-dimensional meteorological parameter field generated in step S2. The entrance cross-section elevation completely covers the north-south range of the mining area horizontally and runs vertically from the bottom to the top of the full CFD computational domain. Each grid point of the entrance cross-section stores the wind speed, wind direction, temperature, humidity, and static pressure parameters at the corresponding height.

[0036] S311. Extract meteorological parameters of the inlet section: Along the inlet section of the full CFD computation domain, extract the wind speed and wind direction data at different heights in the hourly three-dimensional meteorological parameter field output in step S2, and fit the inlet wind profile characteristic parameters at each time.

[0037] S312, Wind Direction Vector Decomposition: Converts the measured wind direction angle into three-dimensional velocity components in the CFD coordinate system, ensuring that the inlet wind direction matches the actual meteorological conditions, and the velocity component calculation formula is as follows: ; In the formula, The velocity component at the inlet section in the x-direction (east-west direction); The velocity component at the inlet section in the y-direction (north-south direction); The average wind speed at the height of the inlet section; It is the clockwise angle between the wind direction and due north.

[0038] S313. Calculation of inlet turbulence parameters: Calculate the inlet turbulent kinetic energy based on the inlet average wind speed and the power exponent of the wind profile. and turbulent dissipation rate : ; ; in, ; In the formula, The corrected inlet turbulence intensity is set to a fixed value of 0.1; This is a constant for the turbulence model, taken as a fixed value of 0.09; The integral scale for turbulence is taken as 0.07 times the inlet calculation height; This represents the original standard turbulence intensity.

[0039] S314. Generate hourly inlet boundary sequence: Generate 24 sets of inlet boundary parameter sequences containing velocity components, turbulent kinetic energy, and turbulent dissipation rate at a time step of 1 hour, forming a dynamic inlet boundary dataset.

[0040] S32. Differentiate between different underlying surface types in the mine pit and calculate the hourly net radiative heat flux of the wall grid.

[0041] S321. Divide the mine pit wall grid into three categories: mine pit bottom, step slope, and spoil heap slope, and then measure the surface albedo according to the zone. and surface emissivity Perform differentiated assignment: ; ; S322. Calculate the net radiative heat flux of the wall mesh: ; In the formula, for Net radiative heat flux of the wall grid at any given time; for Shortwave solar radiation at all times; for Atmospheric longwave radiation at all times; It is the Stefan Boltzmann constant; The surface thermodynamic temperature of the wall is denoted as .

[0042] S33. Based on the three-dimensional meteorological parameter field output in step S2, set the outlet pressure boundary and set the top surface boundary of the full CFD computational domain for each scenario to ensure that the boundary of the full CFD computational domain does not interfere with the internal flow field.

[0043] S331. Set the CFD outlet boundary as a pressure outlet, and adjust the relative static pressure of the outlet boundary according to the outlet-side gas pressure data output in step S2 to match the actual gas pressure gradient in the mining area: ; in, ; In the formula, The height corresponding to the CFD exit boundary Corrected relative static pressure; The reference static pressure is set to 0 Pa. For height Pressure correction at a relative reference altitude; The reference air density is taken as a fixed value of 1.225 kg / m³ under normal operating conditions. 3 ; This is the acceleration due to gravity.

[0044] S332. Define the top boundary of the entire CFD computational domain: In scenarios where cloud coverage is less than 40%, a symmetrical boundary is adopted on the top surface of the entire CFD computing domain to avoid artificial interference of the top surface boundary with the internal airflow. In scenarios where cloud coverage is no less than 40%, the top surface of the computational domain is set as a heat flux wall grid at the height of the cloud condensation nucleus, and long-wave reflected heat flux corresponding to the cloud amount is loaded to simulate the heating effect of the cloud layer on the lower air.

[0045] S34. Boundary temporal discretization and CFD interface integration: Write user-defined functions to automatically load three types of boundary parameters—inlet air profile, wall mesh heat flux, and outlet static pressure—time-by-time, construct a full-element dynamic CFD boundary condition system covering the flow field and thermal field, and integrate with the transient simulation solver in step S5.

[0046] Step S34 specifically includes the following steps: S341, Boundary temporal discretization: Discretize the three types of boundary parameters—inlet wind profile, wall grid heat flux, and outlet static pressure—into 24 transient boundary files with a 1-hour time step, unify the data format and timestamp, and align them with the meteorological field time dimension in step S2.

[0047] S342. Write a user-defined function: Use C language to write the DEFINE_PROFILE macro function to define the boundary update trigger time step, so that the boundary parameters at the corresponding time are read and the solver boundary conditions are updated every time an integer multiple of 1 hour is reached during the simulation.

[0048] S343. Boundary Condition Compilation and Loading: Compile the completed user-defined function into a dynamic link library and load it into the CFD solver to complete the interface between the dynamic boundary and the solver.

[0049] S344, Output a full-element dynamic boundary system: Integrate three types of boundaries: inlet dynamic wind profile, wall grid dynamic heat flux, and outlet static pressure, to form a full-element dynamic CFD boundary condition system covering the flow field and thermal field.

[0050] Step S4 specifically includes the following steps: S41. Based on the monthly mining progress plan, the mine topography and functional area layout are updated synchronously to generate a three-dimensional solid geometric model that matches the actual mining status.

[0051] S411 Monthly Topographic Data Update: Based on the monthly mining progress plan, the mining area excavation range, bench height and slope shape, transportation road direction, spoil heap height and range, and crushing station location are updated synchronously to generate a digital elevation model of the mine pit with a resolution of no less than 5m for the current month.

[0052] S412, 3D Solid Geometry Modeling: Import the updated digital elevation model of the mine pit into the 3D modeling software to construct the solid geometry model of the mine pit. Simultaneously adjust the range and terrain fit of the entire CFD calculation domain to ensure that the solid geometry model of the mine pit is consistent with the actual mining conditions.

[0053] S413, CFD Computation Mesh Adaptation and Update: Based on the updated mine entity geometry model, the CFD computation mesh is updated using a structured base mesh combined with local refinement (using a regular hexahedral structured mesh as the base mesh framework for the entire CFD computation domain, maintaining the basic coarse mesh size in the outer far-field region; subdividing and refining the mesh size for all dust-generating areas such as mine bench walls, mining, blasting, transportation, crushing, and waste disposal, and deploying a gradient boundary layer mesh on all mine walls), ensuring that the mesh accuracy of the near-wall and dust-generating regions meets the requirements of the RNG k-ε turbulence model, and outputting a dynamic mine mesh model.

[0054] S42. Set up the drilling face, blasting face, mining and loading face, transportation road, crushing station and spoil heap face of the open mine as independent face dust sources, and integrate the dust source location, face area, source strength, particle size distribution and dust particle density parameters of each face dust source to form a face dust source parameter library that can be dynamically updated with the mining progress.

[0055] The dust particle size distribution is defined as follows: the dust particle size in all processes is characterized using the Rosing-Lammler distribution, matched with the actual measurement results from a laser particle size analyzer on site; and the dust particle size distribution formula is: ; In the formula, For particle size greater than The proportion of dust mass; Dust particle size; The median particle size of the dust is taken as a fixed value of 10. The particle size distribution index is set to a fixed value of 1.12.

[0056] The formula for calculating the total dust source intensity on transportation routes is as follows: ; ; In the formula, The total dust source intensity (total dust release rate) of the entire transportation route. For linear surface source intensity per unit length of transportation road; This refers to the total length of the dust-generating transportation road.

[0057] Total dust source intensity at the mining and loading face, crushing station and spoil heap The calculation formula is as follows: ; In the formula, Provides a continuous source of power for the unit area of ​​the mining and loading face, crushing station, and spoil heap; This represents the area of ​​the exposed dust-generating area currently under operation. This serves as the baseline coefficient for work intensity. This is the dust conversion factor for the process.

[0058] Total dust intensity at blasting work site The calculation formula is as follows: ; in, ; In the formula, This refers to the total dust mass generated in a single blast. The duration of continuous release of blasting dust; The dust production rate coefficient for explosion-proof bankruptcy; This refers to the total amount of explosives loaded in a single blast.

[0059] The formula for calculating the total dust source intensity at the perforation work surface is as follows: ; in, ; In the formula, For perforation operations, the surface intensity per unit area is increased; This refers to the area of ​​the exposed dust-generating area during drilling operations. This is a dust generation conversion factor specific to the piercing process; This represents the strength coefficient for the perforation operation.

[0060] Dust particle density The expression is as follows: ; S43. Based on the surface dust source parameters output in step S42 and the hourly average wind speed output in step S2, the mirror reflection correction of the step wall and the particle gravity settling correction are introduced to construct the initial concentration field coupled with the dual dust sources, thereby accelerating the convergence speed of the transient simulation.

[0061] S431. Calculation of dust concentration at solid dust sources: Gravity settling is introduced to correct the vertical offset of the concentration center, and a Gaussian smoke flow model is used to calculate the dust concentration at various points in space at the solid dust source. ; in, ; ; ; In the formula, For physical dust sources The dust mass concentration at the point; For the dust to be abundant and strong, its extraction or ; The average wind speed in the mining area; The standard deviation of the lateral concentration distribution; The standard deviation of the vertical concentration distribution; These are the horizontal and vertical coordinates of the physical dust source, respectively. The final settling velocity of dust particles; The horizontal distance from the dust source along the direction of airflow; air density; Aerodynamic viscosity; The transverse turbulent diffusion coefficient; is the vertical turbulent diffusion coefficient.

[0062] S432. Calculation of Virtual Dust Source Concentration Based on Specular Reflection: Based on the principle of specular reflection, the reflection effect of the step wall on dust is equivalent to a virtual dust source at a symmetrical position on the wall. The diffusion concentration of the virtual dust source is calculated using the same diffusion formula as that for physical dust sources. ; In the formula, For virtual dust sources in The dust mass concentration at the point; These are the horizontal and vertical coordinates of the virtual dust source, respectively, which are obtained by symmetry of the physical dust source coordinates about the step wall.

[0063] S433. Calculation of terminal settling velocity of particles: The modified Stokes formula is used to calculate the terminal settling velocity of dust particles of different sizes. ; In the formula, The final settling velocity of the particles; The particle drag coefficient is taken as a fixed value of 0.5 under turbulent conditions.

[0064] S434. Superposition of initial concentration from multiple dust sources: The concentrations of physical dust sources from all processes are linearly superimposed with the corresponding virtual dust source concentrations to obtain the three-dimensional initial dust concentration field of the entire mining area. ; In the formula, For spatial points The initial dust mass concentration; The total number of dust sources is 6; For the first The physical dust source diffusion concentration of an individual dust source; For the first Virtual dust source diffusion concentration for each dust source; For the first The physical dust source diffusion concentration of each dust source.

[0065] S44. Based on the surface dust source parameter library output in step S42, set the discrete phase jet source and particle force parameters.

[0066] S441. Discrete phase jet source setting: Set the surface dust source of each process as a surface jet source. The mass flow rate of the dust source is consistent with the source strength parameter output in step S42, and the particle size follows the Rosing-Rammel distribution set in step S42.

[0067] S442, Two-way Coupling and Force Setting: Considering the reaction force of dust particles on the airflow, the two-way coupling between the discrete phase and the continuous phase is enabled to simulate the motion characteristics of particles in the turbulent field.

[0068] S443, Discrete Phase Boundary Settings: Ground, stepped slope, and spoil heap slope are set as capture boundaries, and particle contact with the wall is determined as settlement; the computational domain exit is set as an escape boundary, and tracking stops when a particle leaves the computational domain.

[0069] S444 Output Discrete Phase Parameter Set: Integrates three types of parameters—jet source, force, and boundary—to form the DPM discrete phase parameter set.

[0070] Step S5 specifically includes the following steps: S51. Driven by the full-element dynamic CFD boundary conditions output in step S3, load the dynamic mine geometric model, surface dust source parameters and initial dust concentration field output in step S4, enable the Bucineske approximation to introduce the thermal buoyancy effect, and construct a turbulent flow field solution system adapted to complex terrain and inversion scenarios.

[0071] S511. The RNG k-ε turbulence model is used to adapt to strong shear flows such as mine pit bench separation flow and slope flow to achieve thermal buoyancy correction. The equation for the transport is expressed as follows: ; in, ; ; ; In the formula, for The air velocity component in the direction of time; It is the reciprocal of the effective Prandtl number corresponding to the turbulent kinetic energy; For effective dynamic viscosity, and , Turbulent viscosity; The turbulent kinetic energy generated by the average velocity gradient; The turbulent kinetic energy generated by buoyancy; For custom turbulent energy items; For general coordinate component subscripts in tensor form, , respectively representing the X, Y, and Z directions of a three-dimensional Cartesian coordinate system; For general coordinate component subscripts in tensor form, , respectively representing the X, Y, and Z directions of a three-dimensional Cartesian coordinate system; The coefficient of thermal expansion of air; The Prandtl number represents the turbulent flow. This is the water spraying wetting attenuation coefficient, and .

[0072] Turbulent dissipation rate The transport equation is expressed as follows: ; in, ; ; In the formula, The effective Prandtl number corresponding to the turbulent dissipation rate is the reciprocal of the turbulent dissipation rate. , , All are constants of the RNG k-ε turbulence model, taken as 1.42, 1.68, and 0.0845 respectively; This is the near-wall correction term for the RNG k-ε turbulence model; For custom turbulent dissipation rate source terms; For strain rate, and ; The critical strain constant is fixed at 4.38. The strain correction factor is fixed at 0.012. Empirical coefficients are generated for the dissipation of operational disturbances; This represents the total dust source strength within the current grid cell. Calculate the mesh cell volume for the current CFD.

[0073] S512. The difference in air density caused by the temperature gradient is quantified by the Businsk approximation, and the thermal buoyancy is equivalent to the volume force source term and loaded into the Reynolds-averaged Navier-Stokes momentum equation to simulate the flow field difference between summer thermal convection and winter temperature inversion.

[0074] S513. Initial flow field loading: Using the initial three-dimensional meteorological field output in step S2 as the initial condition of the flow field, the dynamic boundary output in step S3 at the first moment is loaded to complete the steady-state flow field initialization and provide the initial flow field for transient simulation.

[0075] S52. Load the surface dust source parameter library and DPM discrete phase parameter configuration output from step S4, and carry out multi-particle-size graded gas-solid bidirectional coupling solution for fine particulate matter, inhalable particulate matter, and total suspended particulate matter to accurately restore the diffusion and settling characteristics of dust with different particle sizes.

[0076] S521. Taking into full account the four types of forces—drag, gravity, buoyancy, Saffman lift, and pressure gradient force—the following particle motion control equations are constructed: ; In the formula, This is the particle velocity vector; This is the gas phase velocity vector; This refers to the particle relaxation time. It is the vector of gravitational acceleration; The viscosity of air in motion; and These are the first and second tensors of the flow field strain rate tensor. OK Column components and the first OK Column components.

[0077] S522, Multi-size particle classification tracking: Based on the Rosing-Rammel particle size distribution in step S42, particle groups in four size ranges of PM2.5, PM5, PM10 and TSP are generated, and their motion trajectories and spatial concentrations are tracked respectively to realize the simulation of dust diffusion by particle size.

[0078] S53. Load the initial dust concentration field output in step S43, call the user-defined function in step S34 to update the boundary hourly, and complete the transient gas-solid two-phase coupling simulation of the whole day and night cycle according to the settings in step S44, and output the hourly particle size dust concentration field and turbulent kinetic energy field of the whole mine.

[0079] Step S512 specifically includes the following steps: S5121. Construct the following Bussinesk approximate linear correction formula for air density: ; In the formula, The coefficient of thermal expansion of air is a fixed value. ; As the global reference temperature, take (20℃).

[0080] S5122. The equivalent buoyancy source term in the vertical direction is derived from the air density difference: ; In the formula, It is a volumetric force of thermal buoyancy per unit volume, which acts only on the vertical momentum component along the Z-axis; We can see the physical laws: Summer daytime , This generates upward buoyancy, enhancing vertical convection turbulence in the mine; during winter temperature inversion periods... , It creates a downward pressure, weakening the vertical diffusion ability of dust.

[0081] S5123. Construct the three-dimensional Reynolds-mean momentum equations containing the Businesk buoyancy source term: ; In the formula, The static pressure of the fluid at any location within the entire flow field; for The direction of gravitational acceleration components, and the horizontal components in the X and Y directions. Vertical Z-axis ; The introduced Businesk thermal buoyancy volumetric force source term is used to achieve the regulation of flow field velocity by temperature gradient.

[0082] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. An open cut coal mine dust dispersion simulation method characterised by: Includes the following steps: S1. Deploy a multi-scale meteorological dust collaborative monitoring network in the mining area and obtain a standardized time-series monitoring dataset: Set up a 10-meter-high high-level benchmark meteorological station at the top of the highest spoil heap upstream of the prevailing wind direction, and deploy ground meteorological monitoring points and dust monitoring points within the mining area and within a set range of the mining area to collect meteorological parameters of all elements and dust concentration data of different particle sizes. After outlier removal, consistency verification and standardization processing, a standardized time-series monitoring dataset is generated. S2. Based on the standardized time-series monitoring dataset output in step S1, and combined with the three types of terrain feature parameters extracted from the mine digital elevation model, namely elevation difference, slope aspect, and terrain shading, a three-dimensional inverse distance weighted interpolation model with terrain correction factor is constructed. Vertical constraints are then applied using the atmospheric boundary layer wind power law profile to generate a three-dimensional meteorological parameter field of the entire mining area with horizontal and vertical layers for each hour. S3. Decompose the hourly three-dimensional meteorological parameter field of the entire mining area output in step S2 into three types of boundary parameters: inlet dynamic wind profile, wall grid net radiative heat flux, and outlet static pressure. Automatically load the boundary parameters hourly through user-defined functions to construct a full-element dynamic CFD boundary condition system covering the flow field and thermal field. S4. Based on the monthly mining progress plan, the mine geometric model is updated synchronously to build a surface dust source parameter library covering the entire process of drilling, blasting, mining and loading, road transportation, crushing and disposal. The initial dust concentration field is generated by combining the dual dust source Gaussian diffusion model with specular reflection correction and gravity settlement correction. S5. Driven by the full-element dynamic CFD boundary conditions output in step S3, load the dynamic mine geometric model, surface dust source parameters and initial dust concentration field output in step S4. Couple the discrete phase model with the RNGk-ε turbulence model with Busineske thermal buoyancy correction to carry out transient gas-solid two-phase dust diffusion coupling simulation of fine particulate matter, inhalable particulate matter and total suspended particulate matter of multi-particle size classification, and output the hourly full mine particle size dust concentration field and turbulent kinetic energy field. S6. Use the in-situ monitoring data output in step S1 to update and iterate steps S1-S5.

2. A method of simulating dust dispersion in an open cut coal mine according to claim 1, characterised in that: The full range of meteorological parameters mentioned in step S1 include wind speed, wind direction, temperature, humidity, and air pressure.

3. A method of simulating dust dispersion in an open cut coal mine according to claim 1, characterised in that: Step S2 specifically includes the following steps: S21. Extraction of terrain feature parameters in the mining area: Based on the digital elevation model of the mine pit and the coordinates of each monitoring point and the standardized time-series monitoring data output in step S1, extract the multi-dimensional terrain correlation parameters between the interpolation grid and the monitoring points. S211, traverse each interpolation grid point in the interpolation grid, calculate its horizontal distance, elevation difference and slope parameter with all monitoring points, form point-by-point terrain parameter matrix : ； in, ； ； ； In the formula, For the first Horizontal distance between each monitoring point and the interpolation grid point; For interpolation grid points and the first The absolute value of the elevation difference of each monitoring point; The maximum vertical height difference of the interpolation grid within the monitoring range of the monitoring point; , These represent the planar spans of the slope along the X and Y axes, respectively. The coordinates of the interpolation grid points; For the first Coordinates of each monitoring point; The total number of monitoring points participating in the interpolation; For the first The angle between the slope aspect of the slope where each monitoring point is located and the slope aspect of the slope where the interpolation grid point is located; S212. Calculate the slope aspect of the slope where each monitoring point is located based on the digital elevation model of the mine pit: ； in, ； ； In the formula, For the first The slope direction of each monitoring point; , These are the topographic elevation gradients along the X and Y axes, respectively. , , , These are the elevations of the interpolated grids adjacent to the east, west, north, and south of the monitoring point grid, respectively. The side length of the interpolation grid in the plane; S213. Calculate the terrain obstruction coefficient: ； in, ； ； In the formula, For the first Terrain occlusion correction factor for each monitoring point; For the first The maximum occlusion height difference between each monitoring point and the interpolation grid point; For the first The horizontal straight line segment between each monitoring point and the interpolation grid point The actual terrain elevation; For the first The horizontal distance between each monitoring point and the interpolation grid point is [number] monitoring points. Theoretical spatial elevation at the location; The horizontal length of the distance between the location to be determined and the monitoring point on the horizontal straight line segment between the monitoring point and the interpolation grid point; S214. Statistical dominant wind direction: Based on the full-element meteorological standardized time series data output in step S1, the dominant wind direction at a height of 10 meters within the interpolation period is statistically analyzed and used as the airflow direction benchmark for slope aspect weight correction. S22. Based on the horizontal distance parameters extracted in step S211, calculate the basic interpolation weights of each monitoring point to the interpolation grid points: ； In the formula, For the first The basic interpolation weights of each monitoring point to the interpolation grid points; The distance power exponent; S23. Introduce elevation difference correction factor, slope aspect angle correction factor, and terrain shading correction factor to perform point-by-point correction on the basic interpolation weights obtained in step S22, and obtain the comprehensive interpolation weights: ； in, ； ； In the formula, For the first The comprehensive correction weight of each monitoring point to the interpolation grid points; For the first Elevation difference correction factor corresponding to each monitoring point; For the first Slope aspect angle correction factor corresponding to each monitoring point; For reference attenuation height; The prevailing wind direction during the interpolation period; S24. The full-element meteorological parameter values ​​of all monitoring points are weighted and summed according to the comprehensive correction weight to obtain the horizontal interpolated meteorological parameters of the interpolation grid points: ； In the formula, For the horizontal interpolated meteorological parameter values ​​of the interpolation grid points; For the first Measured meteorological parameter values ​​at each monitoring point; S25. Based on the horizontal interpolation results obtained in step S24, the vertical wind speed distribution is corrected using the atmospheric boundary layer wind power law profile, and the vertical temperature distribution is corrected by combining seasonal and diurnal inversion characteristics. S251. Vertical profile correction of wind speed: The wind speed at different altitudes is constrained using the power law of atmospheric boundary layer wind. ； In the formula, Vertical height of interpolation grid points Corrected wind speed at the location; For reference height Horizontal interpolated wind speed at the location; For reference height, a monitoring height of 10m from the high-level benchmark meteorological station was used; The wind profile power exponent, and ; S252. Vertical Temperature Gradient Correction: The vertical temperature gradient is set according to the season and day / night scene: a positive decreasing gradient is used during the day in summer, with the temperature decreasing by 0.6℃ for every 100m increase; a negative inversion gradient is used at night in winter, with the temperature increasing by 0.3℃ to 0.5℃ for every 100m increase, correcting the temperature distribution across the entire altitude layer. S26. Generation of spatiotemporal dynamic three-dimensional meteorological parameter field: Based on the horizontal interpolation results of step S24 and the vertical correction results of step S25, a time-by-time three-dimensional meteorological parameter field covering the entire CFD computational domain is generated at a fixed time step. The full CFD computational domain is a three-dimensional closed solution space defined by the gas-solid two-phase coupled simulation of open-pit coal mine dust. It is based on the geometric center of the mine pit, with a length and width of 5 times the east-west length of the mine pit. The entrance boundary is 1.5 times the mine length from the west edge of the mine pit, and the exit boundary is 2.5 times the mine length from the east edge of the mine pit. The total width of the computational domain is 5 times the north-south width of the mine pit, with 2 times the mine width reserved on each side. The total height of the computational domain is 4 times the maximum height difference of the mine pit. The full CFD computational domain completely includes the entire production area of ​​the mining area, spoil heap, transportation roads, and crushing station, and buffer spaces are reserved upstream and downstream, on both sides, and above the mine pit.

4. The method for simulating dust diffusion in open-pit coal mines according to claim 3, characterized in that: In step S26, hourly interpolation is first performed with a time step of 1 hour to generate a three-dimensional meteorological parameter sequence for 24 times throughout the day; Then, a three-dimensional regular grid covering the entire CFD computation domain is constructed, with a horizontal resolution of 50m and a vertical resolution of 10m. Each grid point stores five types of meteorological parameters: wind speed, wind direction, temperature, humidity, and air pressure. Finally, the spatial coordinate reference and timestamp are unified to form a complete hourly three-dimensional meteorological parameter field.

5. The method for simulating dust diffusion in open-pit coal mines according to claim 4, characterized in that: Step S3 specifically includes the following steps: S31. Based on the hourly three-dimensional meteorological parameter field output in step S2, generate a three-dimensional dynamic entrance wind profile with turbulence parameters. The entrance cross-section data is a time-series meteorological dataset extracted from all grid points of the vertical entrance elevation on the windward side of the full CFD computational domain from the hourly three-dimensional meteorological parameter field generated in step S2. The entrance cross-section elevation completely covers the north-south range of the mining area horizontally and runs vertically from the bottom to the top of the full CFD computational domain. Each grid point of the entrance cross-section stores the wind speed, wind direction, temperature, humidity, and static pressure parameters at the corresponding height. S311. Extract meteorological parameters of the inlet section: Along the inlet section of the full CFD computation domain, extract the wind speed and wind direction data at different heights in the hourly three-dimensional meteorological parameter field output in step S2, and fit the inlet wind profile characteristic parameters at each time. S312, Wind Direction Vector Decomposition: Converts the measured wind direction angle into three-dimensional velocity components in the CFD coordinate system, ensuring that the inlet wind direction matches the actual meteorological conditions, and the velocity component calculation formula is as follows: ； In the formula, The velocity component at the inlet section in the x-direction; The velocity component at the inlet section in the y-direction; The average wind speed at the height of the inlet section; The angle between the wind direction and due north, clockwise. S313. Calculation of inlet turbulence parameters: Calculate the inlet turbulent kinetic energy based on the inlet average wind speed and the power exponent of the wind profile. and turbulent dissipation rate : ； ； in, ； In the formula, The corrected inlet turbulence intensity; These are constants for the turbulence model; The integral scale is for turbulence. The original standard turbulence intensity; S314. Generate hourly inlet boundary sequence: Generate 24 sets of inlet boundary parameter sequences containing velocity components, turbulent kinetic energy, and turbulent dissipation rate according to a 1-hour time step, forming a dynamic inlet boundary dataset. S32. Differentiate between different underlying surface types of mine pits and calculate the hourly net radiative heat flux of the wall grid. S321. Divide the mine pit wall grid into three categories: mine pit bottom, step slope, and spoil heap slope, and then measure the surface albedo according to the zone. and surface emissivity Perform differentiated assignment: ； ； S322. Calculate the net radiative heat flux of the wall mesh: ； In the formula, for Net radiative heat flux of the wall grid at any given time; for Shortwave solar radiation at all times; for Atmospheric longwave radiation at all times; It is the Stefan Boltzmann constant; The surface thermodynamic temperature of the wall; S33. Based on the three-dimensional meteorological parameter field output in step S2, set the outlet pressure boundary and set the top surface boundary of the full CFD calculation domain for each scenario to ensure that the boundary of the full CFD calculation domain does not interfere with the internal flow field. S331. Set the CFD outlet boundary as a pressure outlet, and adjust the relative static pressure of the outlet boundary according to the outlet-side gas pressure data output in step S2 to match the actual gas pressure gradient in the mining area: ； in, ； In the formula, The height corresponding to the CFD exit boundary Corrected relative static pressure; The reference static pressure is set to 0 Pa. For height Pressure correction at a relative reference altitude; The base air density; It is the acceleration due to gravity; S332. Define the top boundary of the entire CFD computational domain: In scenarios where cloud coverage is less than 40%, the top surface of the entire CFD computing domain adopts a symmetrical boundary. In scenarios where cloud coverage is no less than 40%, the top surface of the computational domain is set as a heat flux wall grid at the height of cloud condensation nuclei, and long-wave reflected heat flux corresponding to the cloud amount is loaded to simulate the heating effect of the cloud layer on the lower air. S34. Boundary temporal discretization and CFD interface integration: Write user-defined functions to automatically load three types of boundary parameters—inlet air profile, wall mesh heat flux, and outlet static pressure—time-by-time, construct a full-element dynamic CFD boundary condition system covering the flow field and thermal field, and integrate with the transient simulation solver in step S5.

6. The method for simulating dust diffusion in open-pit coal mines according to claim 5, characterized in that: Step S34 specifically includes the following steps: S341, Boundary temporal discretization: Discretize the three types of boundary parameters, namely inlet wind profile, wall grid heat flux, and outlet static pressure, into 24 transient boundary files with a time step of 1 hour, unify the data format and timestamp, and align them with the time dimension of the meteorological field in step S2. S342. Write a user-defined function: Use C language to write the DEFINE_PROFILE macro function to define the boundary update trigger time step, so that the boundary parameters at the corresponding time are read and the solver boundary conditions are updated every time an integer multiple of 1 hour is reached during the simulation. S343. Boundary Condition Compilation and Loading: Compile the completed user-defined function into a dynamic link library and load it into the CFD solver to complete the interface between the dynamic boundary and the solver. S344, Output a full-element dynamic boundary system: Integrate three types of boundaries: inlet dynamic wind profile, wall grid dynamic heat flux, and outlet static pressure, to form a full-element dynamic CFD boundary condition system covering the flow field and thermal field.

7. The method for simulating dust diffusion in open-pit coal mines according to claim 6, characterized in that: Step S4 Specifically, the following steps are included: S41. Based on the monthly mining progress plan, the mine pit topography and functional area layout are updated synchronously to generate a three-dimensional solid geometric model that matches the actual mining status. S411, Monthly Topographic Data Update: Based on the monthly mining progress plan, the mining area excavation range, bench height and slope shape, transportation road direction, spoil heap height and range, and crushing station location are updated synchronously to generate a digital elevation model of the mine pit with a resolution of no less than 5m for the current month. S412, 3D solid geometry modeling: Import the updated digital elevation model of the mine pit into the 3D modeling software, construct the solid geometry model of the mine pit, and simultaneously adjust the range and terrain fit of the full CFD calculation domain to ensure that the solid geometry model of the mine pit is consistent with the actual mining conditions. S413, CFD Computational Mesh Adaptation and Update: Based on the updated mine entity geometry model, the CFD computational mesh is updated using a structured basic mesh combined with local refinement to ensure that the mesh accuracy of the near-wall and dust-generating areas meets the requirements of the RNGk-ε turbulence model, and outputs a dynamic mine mesh model. S42. Set the drilling face, blasting face, mining and loading face, transportation road, crushing station and spoil heap face of open mine as independent face dust sources, and integrate the dust source location, face area, source strength, particle size distribution and dust particle density parameters of each face dust source to form a face dust source parameter library that can be dynamically updated with the mining progress. The dust particle size distribution is defined as follows: the dust particle size in all processes is characterized using the Rosing-Lammler distribution, matched with the actual measurement results from a laser particle size analyzer on site; and the dust particle size distribution formula is: ； In the formula, For particle size greater than The proportion of dust by mass; Dust particle size; The median particle size of the dust. It is the particle size distribution index; The formula for calculating the total dust source intensity on transportation routes is as follows: ； ； In the formula, The total dust source intensity of the entire transportation route; For linear surface source intensity per unit length of transportation road; This refers to the total length of the dust-generating transportation road; Total dust source intensity at the mining and loading face, crushing station and spoil heap The calculation formula is as follows: ； In the formula, Provides a continuous source of power for the unit area of ​​the mining and loading face, crushing station, and spoil heap; This represents the area of ​​the exposed dust-generating area currently under operation; This serves as the baseline coefficient for work intensity. This is the dust generation conversion factor for the process; Total dust intensity at blasting work site The calculation formula is as follows: ； in, ； In the formula, The total dust mass generated in a single blast; The duration of continuous release of blasting dust; The dust production rate coefficient for explosion-proof bankruptcy; This refers to the total amount of explosives loaded in a single blast. The formula for calculating the total dust source intensity at the perforation work surface is as follows: ； in, ； In the formula, For perforation operations, the surface intensity per unit area is increased; This refers to the area of ​​the exposed dust-generating area during drilling operations. This is a dust generation conversion factor specific to the piercing process; The strength coefficient for the perforation operation; Dust particle density The expression is as follows: ； S43. Based on the surface dust source parameters output in step S42 and the hourly average wind speed output in step S2, the step wall mirror reflection correction and particle gravity settling correction are introduced to construct the initial concentration field coupled with the dual dust sources, which accelerates the convergence speed of transient simulation. S431. Calculation of dust concentration at solid dust sources: Gravity settling is introduced to correct the vertical offset of the concentration center, and a Gaussian smoke flow model is used to calculate the dust concentration at various points in space at the solid dust source. ； in, ； ； ； In the formula, For physical dust sources The dust mass concentration at the point; For the dust to be abundant and strong, its extraction or ; The average wind speed in the mining area; The standard deviation of the lateral concentration distribution; The standard deviation of the vertical concentration distribution; These are the horizontal and vertical coordinates of the physical dust source, respectively. The final settling velocity of dust particles; The horizontal distance from the dust source along the direction of airflow; air density; Aerodynamic viscosity; The transverse turbulent diffusion coefficient; The vertical turbulent diffusion coefficient; S432. Calculation of Virtual Dust Source Concentration Based on Specular Reflection: Based on the principle of specular reflection, the reflection effect of the step wall on dust is equivalent to a virtual dust source at a symmetrical position on the wall. The diffusion concentration of the virtual dust source is calculated using the same diffusion formula as that for physical dust sources. ； In the formula, For virtual dust sources in The dust mass concentration at the point; These are the horizontal and vertical coordinates of the virtual dust source, respectively. S433. Calculation of terminal settling velocity of particles: The modified Stokes formula is used to calculate the terminal settling velocity of dust particles of different sizes. ； In the formula, The final settling velocity of the particles; This is the particle drag coefficient; S434. Superposition of initial concentration from multiple dust sources: The concentrations of physical dust sources from all processes are linearly superimposed with the corresponding virtual dust source concentrations to obtain the three-dimensional initial dust concentration field of the entire mining area. ； In the formula, For spatial points The initial dust mass concentration; Total number of dust sources; For the first The physical dust source diffusion concentration of an individual dust source; For the first Virtual dust source diffusion concentration for each dust source; For the first The physical dust source diffusion concentration of an individual dust source; S44. Based on the surface dust source parameter library output in step S42, set the discrete phase jet source and particle force parameters; S441. Discrete phase jet source setting: Set the surface dust source of each process as a surface jet source. The mass flow rate of the dust source is consistent with the source strength parameter output in step S42, and the particle size follows the Rosing-Rammel distribution set in step S42. S442, Two-way coupling and force setting: Considering the reaction force of dust particles on the airflow, the two-way coupling between the discrete phase and the continuous phase is enabled to simulate the motion characteristics of particles in the turbulent field. S443, Discrete Phase Boundary Settings: Ground surface, stepped slope, and spoil heap slope are set as capture boundaries, and particle contact with the wall surface is determined as settlement; the computational domain exit is set as an escape boundary, and tracking stops when a particle leaves the computational domain; S444 Output Discrete Phase Parameter Set: Integrates three types of parameters—jet source, force, and boundary—to form the DPM discrete phase parameter set.

8. The method for simulating dust diffusion in open-pit coal mines according to claim 7, characterized in that: Step S5 specifically includes the following steps: S51. Driven by the full-element dynamic CFD boundary conditions output in step S3, load the dynamic mine geometric model, surface dust source parameters and initial dust concentration field output in step S4, enable the Businesk approximation to introduce the thermal buoyancy effect, and construct a turbulent flow field solution system that is adapted to complex terrain and inversion scenarios. S511. The RNG k-ε turbulence model is used to adapt the strong shear flow in the mine pit steps to achieve thermal buoyancy correction. The equation for the transport is expressed as follows: ； in, ； ； ； In the formula, for The air velocity component in the direction of time; It is the reciprocal of the effective Prandtl number corresponding to the turbulent kinetic energy; For effective dynamic viscosity, and , Turbulent viscosity; The turbulent kinetic energy generated by the average velocity gradient; The turbulent kinetic energy generated by buoyancy; For custom turbulent energy items; For general coordinate component subscripts in tensor form, , respectively representing the X, Y, and Z directions of a three-dimensional Cartesian coordinate system; For general coordinate component subscripts in tensor form, , respectively representing the X, Y, and Z directions of a three-dimensional Cartesian coordinate system; The coefficient of thermal expansion of air; The Prandtl number for turbulence; This is the water wetting attenuation coefficient, and ; Turbulent dissipation rate The transport equation is expressed as follows: ； in, ； ； In the formula, The effective Prandtl number corresponding to the turbulent dissipation rate is the reciprocal of the turbulent dissipation rate. , , All are constants of the RNG k-ε turbulence model; This is the near-wall correction term for the RNG k-ε turbulence model; For custom turbulent dissipation rate source terms; For strain rate, and ; It is the critical strain constant; This is the strain correction factor; Empirical coefficients are generated for the dissipation of operational disturbances; This represents the total dust source strength within the current grid cell. Calculate the mesh cell volume for the current CFD. S512. The difference in air density caused by the temperature gradient is quantified by the Businesk approximation, and the thermal buoyancy is equivalent to the volume force source term and loaded into the Reynolds-averaged Navier-Stokes momentum equation to simulate the flow field difference between summer thermal convection and winter temperature inversion. S513, Initial flow field loading: Using the initial three-dimensional meteorological field output in step S2 as the initial condition of the flow field, the dynamic boundary output in step S3 at the first moment is loaded to complete the steady-state flow field initialization and provide the initial flow field for transient simulation; S52. Load the surface dust source parameter library and DPM discrete phase parameter configuration output from step S4, and carry out multi-size graded gas-solid bidirectional coupling solution for fine particulate matter, inhalable particulate matter, and total suspended particulate matter. S521. Taking into full account the four types of forces—drag, gravity, buoyancy, Saffman lift, and pressure gradient force—the following particle motion control equations are constructed: ； In the formula, This is the particle velocity vector; This is the gas phase velocity vector; This refers to the particle relaxation time. It is the vector of gravitational acceleration; The viscosity of air in motion; and These are the first and second tensors of the flow field strain rate tensor. OK Column components and the first OK Column component; Gas phase velocity vector The velocity gradient tensor; S522, Multi-size particle classification tracking: Based on the Rosing-Rammel particle size distribution in step S42, particle groups in four size ranges of PM2.5, PM5, PM10 and TSP are generated, and their motion trajectories and spatial concentrations are tracked respectively to realize the simulation of dust diffusion by particle size. S53. Load the initial dust concentration field output in step S43, call the user-defined function in step S34 to update the boundary hourly, and complete the transient gas-solid two-phase coupling simulation of the whole day and night cycle according to the settings in step S44, and output the hourly particle size dust concentration field and turbulent kinetic energy field of the whole mine.

9. The method for simulating dust diffusion in open-pit coal mines according to claim 8, characterized in that: Step S512 details Includes the following steps: S5121. Construct the following Bussinesk approximate linear correction formula for air density: ； In the formula, The coefficient of thermal expansion of air is a fixed value. ; This serves as the global reference temperature. S5122. The equivalent buoyancy source term in the vertical direction is derived from the air density difference: ； In the formula, The thermal buoyancy force per unit volume is the volumetric force. S5123. Construct the three-dimensional Reynolds-mean momentum equations containing the Businesk buoyancy source term: ； In the formula, The static pressure of the fluid at any location within the entire flow field; for The components of gravitational acceleration in the X and Y directions, and the horizontal components in the X and Y directions. Vertical Z-axis ; This is the introduced Businesk thermal buoyancy volumetric force source term.

Images