Automatic modeling method and system for complex terrains and buildings / structures considering surface roughness

Abstract

The invention discloses an automatic modeling method and system for complex terrains and buildings / structures in consideration of surface roughness, and the method is realized based on an open source tool, and comprises the steps: collecting the types of terrains, buildings and roughness, and the actual measurement data of a wind field; automatically processing the topographic data through a script and generating an STL file with a transition section; simplifying the building contour and removing abnormal points by using an RDP algorithm and a standard deviation method; filtering the small building and calculating equivalent roughness; generating a building STL in combination with an open source modeling tool; a roughness surface patch is generated on the basis of the DaviportWieringa classification and the equivalent roughness; fusing the terrain, the building and the roughness surface patch; and in OpenFOAM, accurate binding of roughness and computational grids is realized through grid patch association, and wind field simulation and verification are carried out. According to the method, full-process automation is achieved, the modeling efficiency can be improved, the precision deviation of near-surface wind field simulation is reduced through a double-path roughness processing mechanism, and the problems that a traditional method is time-consuming in modeling, insufficient in precision and dependent on commercial software are effectively solved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of wind field automatic modeling and CFD preprocessing, and particularly relates to an automatic modeling method and system for complex terrain and buildings / structures considering ground roughness. BACKGROUND

[0002] In the field of wind engineering computational fluid dynamics (CFD), high-precision modeling of complex terrain and building / structure groups is the core prerequisite for wind field simulation, pollutant dispersion prediction and wind resource assessment, but the existing technical process has long been subject to low-efficiency manual operation and tool chain fault. Terrain data (such as NASA SRTM elevation data) is generally stored in HGT binary format, and researchers need to manually write parsing scripts to convert it into CFD compatible format (such as STL), and rely on empirical rules to add boundary transition sections (usually need to be expanded by 5-10 times the characteristic height) to avoid numerical instability caused by backflow at the boundary of the calculation domain. According to statistics, the terrain data processing link alone accounts for 40% of the total research time, becoming a key bottleneck restricting industrial-level large-scale simulation.

[0003] Urban building / structure modeling faces more severe precision problems: open geographic data (such as GeoJSON files derived from OpenStreetMap) often carry noise vertices left over from LiDAR scanning, which will cause local cell distortion if directly used for mesh generation, leading to solver divergence; and if small buildings / structures are blindly removed to control the grid size, the overall roughness length ( ) of the building / structure group will be underestimated, significantly weakening the reliability of the near-surface wind speed profile prediction.

[0004] Current commercial solutions, such as using ArcGIS + ANSYS SpaceClaim, cannot automatically perform terrain-building / structure penetration Boolean fusion (need to manually repair millions of face intersections), and ignore the dynamic correlation of building / structure geometry simplification and compensation; secondly, the closed-source architecture blocks deep integration with open-source CFD tool chains (such as OpenFOAM), forcing engineers to rely on temporary scripts to bridge data, forming a "modeling-simulation" workflow breakpoint. SUMMARY

[0005] Therefore, the purpose of the present application is to provide an automatic modeling method and system for complex terrain and buildings / structures considering ground roughness, which not only overcomes the time-consuming problem of existing terrain data and building / structure data processing, but also accurately considers ground roughness, improves the accuracy of near-surface wind prediction, and provides technical support for automatic modeling of complex terrain and buildings and CFD simulation preprocessing.

[0006] To achieve the above purpose, the present application provides the following technical solutions: This invention first proposes an automated modeling method for complex terrain and buildings / structures that considers surface roughness, comprising the following steps: Step 1: Data Collection The system collects topographic elevation data, urban building data, ground roughness category data, and measured wind field data. The topographic elevation data is in HGT format; the urban building data is in GeoJSON format; the ground roughness category data is in GeoTiff format; and the measured wind field data is collected through a wind measurement tower. Step Two: Data Processing The collected terrain elevation data, urban building data, and ground roughness category data were processed separately, including: The method for processing terrain elevation data is as follows: use Python or Matlab scripts to convert binary terrain elevation data into point cloud data in ASCII format, use an improved Witoszynski wind tunnel contraction curve to generate transition sections around the terrain, and convert the point cloud data into triangular mesh data through Delaunay triangulation to generate STL files of the terrain. The processing method for urban building data is as follows: the Ramer-Douglas-Peucker algorithm is used to simplify the outline coordinates of buildings / structures, remove redundant points, and combine the standard deviation method to delete outlier points; small buildings / structures are filtered, the windward area and average height of the filtered buildings are recorded, and the equivalent roughness is calculated based on the Lettau formula; the shape coordinates of buildings / structures are extracted using Python scripts, OpenSCAD code is generated using the solidPython2 library, OpenSCAD is imported and negatively stretched to generate the STL file of the building / structure; The method for processing ground roughness category data is as follows: call Python's gdal library to read ground roughness category data in GeoTiff format, establish a mapping between land use type and roughness length based on Davenport classification and Wieringa correction, assign roughness to the terrain, calculate the equivalent roughness based on Lettau formula, generate roughness patches through Delaunay triangulation, and output a roughness attribute file with patch ID and roughness value association. Step 3: Data Fusion Use OpenSCAD's union() function to merge terrain STL files, building / structure STL files, and roughness patches, enabling Boolean operation verification and watertightness checks; Step 4: CFD Simulation Computational fluid dynamics simulations were performed using OpenFOAM, including setting up the computational domain, mesh generation, and computational settings. Roughness patches were defined in snappyHexMeshDict, and a Python script was written to batch assign roughness values ​​to the boundary condition file. The pimpleFoam solver and the SST k-ω model were used with wall roughness correction enabled. Measured wind profile information was applied at the inlet. The roughness length was achieved by mapping the roughness patches to the boundary conditions and associating them with the mesh. Binding to the computational grid and performing wind field calculations; Step 5: Result Output and Verification Output simulation results to verify the deviation between near-surface wind field information and measured wind field information, and ensure that the deviation does not exceed the preset threshold.

[0007] Furthermore, in step two, the formula for the transition curve generated based on the improved Witoszynski wind tunnel contraction curve is as follows: in: and These represent the horizontal and vertical coordinates of the transition curve, respectively. Indicates the height of the transition curve; Indicates the length of the transition curve; This represents the constant coefficient.

[0008] Furthermore, in step two, the epsilon value of the Ramer-Douglas-Peucker algorithm is set to 0.01-0.05 times the perimeter of the building / structure outline to control the degree of simplification; the standard deviation method uses a z-core value greater than ±3 as the criterion for outlier judgment and replaces it with the average height of buildings within a 300-meter radius.

[0009] Furthermore, in step two, the criteria for filtering small buildings / structures are that the roof area is less than a first set threshold or the height is less than a second set threshold; the equivalent roughness is calculated using the Lettau formula: in: This is the ratio of the projected area of ​​the windward side of a building / structure to the horizontal area of ​​the target region. This represents the average height of all filtered buildings / structures within the corresponding building area; and: in: This indicates the projected area of ​​the windward side of a building / structure.

[0010] Furthermore, in step two, the negative stretching extends downwards to z=-5 meters to ensure that the building / structure can penetrate the terrain and achieve a seamless connection with the terrain; the OpenSCAD code is automatically generated by the solidPython2 library and multiple buildings / structures are processed in batches.

[0011] Furthermore, in step two, the roughness length mapping of the Davenport classification corrected by Wieringa includes: a roughness length of 0.0002 meters for water bodies, a roughness length of 0.005 meters for bare land, and a roughness length of 0.03 meters for grassland; the triangular mesh area of ​​the roughness patch is no more than 100 square meters, and the normal direction is uniform.

[0012] Furthermore, in step three, self-intersection is avoided by checking Boolean operations in OpenSCAD, and water tightness is checked to ensure no leakage surfaces; the fused model is used for preprocessing in OpenFOAM simulation.

[0013] Furthermore, in step four, the CFD simulation includes: Computational domain settings: The blockage rate is less than 3%. The planar size of the computational domain is extended along the downwind direction by no less than 5 times the longest side of the target area along the downwind direction, and along the crosswind direction by no less than 3 times the longest side of the target area along the crosswind direction. The minimum height is 50 meters below the lowest point of the terrain, and the maximum height is no less than 1500 meters or 10 times the height of the tallest building. Mesh generation: A background mesh and local densification strategy is adopted. A structured background mesh is constructed using blockMesh, and adaptive densification around buildings and terrain is achieved by combining snappyHexMesh. The resolution of the background mesh is 1 / 300 of the size of the computational domain plane. The vertical direction adopts a non-uniform distribution. The mesh size is 5 meters in the 0-50 meter range near the ground, the mesh size increases by 1.1 times in the 50-500 meter range, and remains at 15 meters above 500 meters. Calculation settings: The pimpleFoam solver was selected, the governing equations were the three-dimensional Reynolds-averaged Navier-Stokes equations, the turbulence model was the SST k-ω model, roughness correction was enabled, and the boundary conditions included inlet wind speed profile, outlet free outflow, top slip wall and crosswind periodic boundary.

[0014] Furthermore, in step four, the roughness patch is mapped using boundary conditions and associated with the mesh to achieve the roughness length. The steps for binding with the computational grid are as follows: 41) Surface Data Format Conversion and Import Standardization: Convert triangular mesh surfaces containing roughness categories to OpenFOAM-compatible .stl format, and unify the coordinate system using the surfaceTransformPoints tool; write a Python script to read the geometric data of the triangular mesh surfaces and their corresponding coordinates. Value, generate "patch ID - roughness length" "Associated files; 42) Mesh and Patch Association Configuration: By configuring snappyHexMeshDict to work collaboratively with tools, ensure that the mesh fits the patch outline and generates a dedicated boundary patch, specifically: 421) In snappyHexMeshDict, define the roughness patch as a feature surface, specifying its surface type, file path, corresponding boundary patch information, and mesh refinement parameters; 422) Extract the feature edges of the roughness patch using the surfaceFeatureExtract tool and generate a feature edge file in .eMesh format; 423) Run blockMesh to build the background mesh, then run snappyHexMesh to generate associated meshes and boundary patches. In the splitting stage, split the mesh near the patch at the specified level. In the bonding stage, ensure that the surface mesh bonding error is ≤1 meter through iteration. 424) After the mesh is generated, use the checkMesh tool to verify whether the type of the roughness-specific patch is "wall" and to verify the integrity of the patch association; 43) Boundary condition mapping for roughness parameters: Modify the turbulenceProperties file to enable the wall roughness correction function of the SST k-ω model; in the boundary condition dictionary, set the velocity field on the roughnessPatch to no-slip condition, set the boundary condition type for turbulent kinetic energy and specific dissipation rate to wallFunction, and specify the use of kqRWallFunction and omegaWallFunction. 44) Batch assignment of roughness parameters based on patch ID: Develop an OpenFOAM auxiliary tool or write a Python script to read "patch ID - roughness length". "Associate files, establish topological associations between triangular face IDs and mesh faces in the roughnessPatch, and..." Values ​​are written in batches to the 0 / k and 0 / omega dictionaries to assign values ​​to the roughnessLength parameter, and the correctness of the assignment results is verified by the paraView tool; 45) Roughness Activation Mechanism During Simulation: During CFD simulation, the SST k-ω turbulence model calls the kqRWallFunction and omegaWallFunction wall functions to read the roughnessLength parameter of the roughnessPatch in real time, and introduces... The formula for calculating the friction velocity of the correction term corrects the turbulence field and wind speed field in the near-wall region, so that the simulated near-surface wind speed profile and turbulence intensity distribution match the actual surface roughness characteristics. The formula for calculating friction speed is as follows: in: The friction speed; Kármán's constant; Wind speed near the wall; To calculate the height of the point.

[0015] This invention also proposes an automated modeling system for complex terrain and buildings / structures that considers surface roughness, for performing the method described above, characterized by comprising: The data acquisition module is used to collect terrain elevation data, urban building data, ground roughness category data, and measured wind field data from open-source databases. Specifically, the terrain elevation data is in HGT format; the urban building data is in GeoJSON format; the ground roughness category data is in GeoTiff format; and the measured wind field data is collected through a wind measurement tower. The terrain processing module uses Python or Matlab scripts to convert terrain elevation data into point cloud data, and uses an improved Witoszynski wind tunnel contraction curve to generate transition sections, and generates STL files of terrain through Delaunay triangulation. The building processing module uses the Ramer-Douglas-Peucker algorithm and standard deviation method to process building / structure data, remove redundant and outlier points, and generate STL files of buildings / structures using the solidPython2 library and OpenSCAD. A filtering module is used to filter small buildings / structures and calculate the equivalent roughness based on the Lettau formula; The roughness processing module uses Python's gdal library to process ground roughness category data, assigns terrain roughness based on Davenport-Wieringa classification, generates roughness patches by combining equivalent roughness, and outputs a roughness attribute file. The fusion module uses OpenSCAD's union() function to merge terrain STL, building STL, and roughness patches, and performs Boolean operations for verification and watertightness checks. The CFD simulation module uses OpenFOAM for wind field simulation, including setting the computational domain, mesh generation, and computational parameters. Roughness patches are defined using snappyHexMeshDict, and a Python script is written to assign roughness values ​​to the boundary condition file. The roughness length is then achieved by mapping the roughness patches to boundary conditions and associating them with the mesh. Binding to the computational grid and performing wind field calculations; The verification module is used to output simulation results, verify the deviation between the near-surface wind field information and the measured wind field information, and ensure that the deviation does not exceed the preset threshold.

[0016] The beneficial effects of this invention are as follows: The present invention provides an automated modeling method for complex terrain and buildings / structures that takes into account surface roughness, and has the following technical effects.

[0017] (1) Full-process automation and efficiency leap have been achieved. This invention uses scripting tools such as Python to connect the traditionally commercial software-dependent and highly manual processes such as terrain analysis, building modeling, and model fusion into an automated pipeline, fundamentally solving the workflow gap problem of "modeling-simulation". It can effectively compress the preprocessing cycle, improve efficiency, and clear the key bottleneck for industrial-scale simulation.

[0018] (2) A precise processing paradigm for surface roughness was established, significantly improving simulation accuracy. This invention innovatively adopts a dual-path fusion mechanism of "inherent terrain roughness + equivalent roughness of small buildings": on the one hand, roughness length is automatically assigned to the terrain based on Davenport-Wieringa classification; on the other hand, the geometric influence of filtered small buildings is scientifically transformed into equivalent roughness using the Lettau formula. Most importantly, through a complete set of "mesh-patch" association and boundary condition mapping technology, lossless transfer and precise binding of roughness data to the OpenFOAM solver are achieved, making the roughness effect truly effective in the simulation. Ultimately, this effectively improves simulation accuracy and reduces the deviation between the simulation results and measured data of near-surface wind speed profiles and turbulence intensity.

[0019] (3) Achieving a synergistic optimization of modeling complexity and computational cost. By filtering out small buildings and replacing their influence with equivalent roughness, the number of building STL files in dense areas can be effectively reduced while ensuring physical realism, and the total amount of computational domain mesh can be reduced. This effectively avoids Boolean operation failures and mesh distortion, and shortens the mesh generation time.

[0020] (4) It breaks down technical barriers and lowers the application threshold. Based on open-source toolchains and public databases, this invention does not rely on any commercial software, which can effectively save procurement costs. This invention provides a standardized and reproducible technical path, enabling small and medium-sized teams and research institutions to conduct high-precision wind field simulations without a strong CFD background, which strongly promotes the popularization and application of this technology. Attached Figure Description

[0021] To make the objectives, technical solutions, and beneficial effects of this invention clearer, the following figures are provided for illustration: Figure 1 A flowchart illustrating the automated modeling method for complex terrain and buildings / structures that takes into account surface roughness, as described in this invention. Figure 2 This is a field diagram obtained from actual measurement data; Figure 3 A chart for Davenport roughness classification (corrected by Wieringa); Figure 4 This is a generated map of the terrain transition section; Figure 5 Modeling diagrams of terrain and buildings; Figure 6 A diagram illustrating a filtration method for small buildings and structures. Figure 7 A flowchart for determining the characteristics of the incoming wind field by integrating measured data and RANS simulations; Figure 8 Error comparison between RANS simulation inlet profile fitting results and measured results. Detailed Implementation

[0022] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, so that those skilled in the art can better understand and implement the present invention. However, the embodiments described are not intended to limit the present invention.

[0023] like Figure 1 As shown, this embodiment describes an automated modeling method for complex terrain and buildings / structures that considers surface roughness, including the following steps: Step 1: Data collection.

[0024] Collect topographic elevation data, urban building data, ground roughness category data, and measured wind field data; the topographic elevation data is in HGT format.

[0025] The terrain elevation data was obtained as HGT format binary data from the NASA:SRTM database (publicly available) for terrain modeling. The NASA:SRTM database can be found at: https: / / earthexplorer.usgs.gov / .

[0026] The city building data was obtained by downloading GeoJSON format files from OpenStreetMap (publicly available), and extracting the building outlines (x, y) and heights (z). The OpenStreetMap database can be found at: https: / / www.openstreetmap.org / .

[0027] The ground roughness category data is GeoTIFF format land use data obtained through Tsinghua University's Nebula Data Service Platform (publicly available), including categories such as cultivated land, grassland, water bodies, and densely built-up areas. Ground roughness category data: https: / / data-starcloud.pcl.ac.cn / .

[0028] The measured wind field data includes two locations within the target wind field: one at the inlet, used to fit the inlet profile; and the other within the target area, used to verify the simulation results. Wind field information was sampled at the target wind field using a wind measurement tower. The tower simultaneously acquired 10-minute average wind speed, wind direction data, and corresponding turbulence intensity data at five vertical heights. The measured data in this example is shown in the figure below. Figure 2 As shown.

[0029] Step Two: Data Processing The collected terrain elevation data, urban building data, and ground roughness category data were processed separately.

[0030] (1) Methods for processing terrain elevation data.

[0031] Convert binary terrain elevation data into ASCII point cloud data using Python or Matlab scripts. Specifically, convert DEM binary elevation data downloaded from public websites into point cloud data in Cartesian coordinates using a script (Python / Matlab), and output an ASCII format point cloud file (each line stores the x, y, z coordinates of one point).

[0032] An improved Witoszynski wind tunnel contraction curve (W50 curve) was used to add a transition segment around the terrain in the terrain simulation to avoid numerical instability caused by abrupt changes in terrain boundaries. The point cloud data was converted into triangular mesh data through Delaunay triangulation, and then the STL file of the terrain was obtained.

[0033] Specifically, the formula for the transition curve generated based on the improved Witoszynski wind tunnel contraction curve is as follows: in: and These represent the horizontal and vertical coordinates of the transition curve, respectively. Indicates the height of the transition curve; Indicates the length of the transition curve; This represents a constant coefficient, which in this embodiment is equal to 1 / 3.

[0034] In this embodiment, the transition curve formula used is: The specific operating steps are as follows: Figure 3 As shown.

[0035] (2) Methods for processing urban building data.

[0036] ① Data Processing. The Ramer-Douglas-Peucker (RDP) algorithm is used to simplify the coordinate points of the building / structure outline, removing redundant points to reduce LiDAR measurement noise. Anomalous outliers are removed using the standard deviation method to avoid sharp angles that could cause mesh distortion. In this embodiment, the RDP algorithm is used to simplify redundant points, specifically by setting the epsilon(ε) value (a threshold controlling the degree of simplification). In this embodiment, a proportional value is set. Specifically, in this embodiment, the epsilon value of the Ramer-Douglas-Peucker algorithm is set to 0.01-0.05 times the perimeter of the building / structure outline to control the degree of simplification. The standard deviation method uses a z-core value greater than ±3 as the criterion for outlier detection, and replaces it with the average height of buildings within a 300-meter radius.

[0037] ② Filtering Buildings / Structures. Small buildings / structures are filtered, and the windward area and average height of the filtered buildings are recorded to calculate the equivalent roughness, thus expressing the impact of the filtered small buildings / structures through equivalent roughness. In this embodiment, small buildings / structures are filtered, and the windward area and average height of the filtered buildings are recorded. The equivalent roughness is calculated based on the Lettau formula. Specifically, the criteria for filtering small buildings / structures are that the roof area is less than a first set threshold or the height is less than a second set threshold. In this embodiment, the first set threshold is 25 m. 2 The second threshold is set at 5 μm. The equivalent roughness is calculated using the Lettau formula: in: This is the ratio of the projected area of ​​the windward side of a building / structure to the horizontal area of ​​the target region. This represents the average height of all filtered buildings / structures within the corresponding building area; and: in: This indicates the projected area of ​​the windward side of a building / structure. For example... Figure 6 As shown, the projected area of ​​the windward side of the building / structure is: ③ Generate STL file. Extract the building / structure's shape coordinates using a Python script, generate OpenSCAD code using the solidPython2 library, import the generated code into OpenSCAD, and perform negative stretching (to z=-5 m) on the building / structure in OpenSCAD to ensure it can penetrate the terrain. Export the building / structure's STL file. Specifically, in this embodiment, negative stretching extends downwards to z=-5 meters to ensure the building / structure can penetrate the terrain and achieve a seamless connection with it. The OpenSCAD code is automatically generated using the solidPython2 library, and multiple buildings / structures are processed in batches.

[0038] (3) Methods for processing ground roughness category data.

[0039] The roughness processing is divided into two parts: the first is to call Python's gdal library to read GeoTiff data and establish a "land use type-" model based on Davenport classification (corrected by Wieringa). "Mapping, automatically assigning to each pixel" The first method uses a value to process the terrain roughness; the second method uses a formula. The small structures being filtered are subjected to equivalent roughness treatment. The resulting... The values ​​are used for subsequent mesh boundary condition settings. Then, the point cloud is converted into continuous triangular mesh patches using the Delaunay triangulation algorithm. At the same time, it is necessary to ensure that the normal directions of the triangular mesh patches are consistent, and generate a roughness property file in .stl format with patch ID and roughness value associated with each patch, i.e., the "patch ID-roughness value" association file (roughnessProperties.txt).

[0040] Specifically, in this embodiment, the roughness length mapping of the Davenport classification modified by Wieringa includes: a roughness length of 0.0002 meters for water bodies, a roughness length of 0.005 meters for bare land, and a roughness length of 0.03 meters for grassland; the area of ​​the triangular mesh of the roughness patch is no more than 100 square meters, and the normal direction is uniform.

[0041] Step 3: Data Fusion The OpenSCAD `union()` function is used to merge the terrain's STL file, the building / structure's STL file, and roughness patches. Boolean operations are enabled to prevent self-intersections, and a watertightness check is performed again to ensure no leaks. The merged model is then used for preprocessing in OpenFOAM simulations. Triangular mesh patches with roughness information are applied to the OpenFOAM simulation, and finally, the roughness length is achieved through boundary condition mapping and mesh-patch association. Binding to the computational grid.

[0042] Step 4: CFD simulation.

[0043] Computational fluid dynamics simulations were performed using OpenFOAM, including wind field calculations. This included setting up the computational domain, mesh generation, and computational settings. Roughness patches were defined in snappyHexMeshDict, and a Python script was written to batch assign roughness values ​​to the boundary condition file. The pimpleFoam solver and the SST k-ω model were selected to enable wall roughness correction. Measured wind profile information was applied at the inlet.

[0044] Specifically, in this embodiment, the roughness length is achieved by mapping the roughness patch to the boundary condition and associating it with the mesh. The steps for binding with the computational grid are as follows.

[0045] 41) Surface Data Format Conversion and Import Standardization: Convert triangular mesh surfaces containing roughness categories to OpenFOAM-compatible .stl format, and unify the coordinate system (to match the computational domain coordinate system) using the surfaceTransformPoints tool. Then, write a Python script to read the geometric data of the triangular mesh surfaces and their corresponding coordinates. Values ​​(terrain roughness from the Davenport-Wieringa classification mapping, equivalent roughness of small buildings from formula calculation) are used to generate "patch ID - roughness length". "The associated file (named roughnessProperties.txt) must clearly specify the unique ID of each facet and its corresponding..." The value (unit: m) provides a data basis for subsequent batch assignments.

[0046] 42) Configuration of mesh and patch association: By configuring snappyHexMeshDict to work in conjunction with tools, ensure that the mesh fits the outline of the patch and generates a dedicated boundary patch, as follows.

[0047] 421) In snappyHexMeshDict, define the roughness patch as a feature surface, specify its surface type, file path, corresponding boundary patch information, and mesh refinement parameters.

[0048] Specifically, open the system / snappyHexMeshDict configuration file in the computation domain, add a dedicated configuration for roughness patches under the surfaces dictionary, and specify their "feature surface" attribute, as follows: surfaces { roughnessPatchSurface / / The name of the patch surface, which needs to be associated with the name of subsequent patches. { type triSurfaceMesh; / / Surface type is triangular mesh file "constant / triSurface / roughnessPatch_aligned.stl"; / / Points to the patch file processed in 8.1 / / Define the boundary patch information corresponding to the patch regions { roughnessPatch / / The final name of the generated boundary patch (must match the subsequent boundary condition dictionary) { `type patch;` / / Sets the base patch type to "patch", which will need to be changed to "wall" later. } } / / Mesh refinement parameters: Ensure finer mesh near the facet to improve bonding accuracy refineLevel 4; / / Mesh refinement level for the area (consistent with the terrain refinement level, e.g., level 4) minLevel 2; / / The lowest mesh refinement level for the patch area to avoid excessively coarse local meshes. distance 10; / / Automatically activates encryption within 10m of the patch, covering key areas around the patch. } / / At the same time, define the STL surfaces of terrain and buildings here to achieve multi-entity collaborative encryption. terrainSurface { ...} buildingSurface { ...} } 422) Extract the feature edges of the roughness patch using the surfaceFeatureExtract tool (connect the surfaceFeatureExtract tool) and generate a feature edge file in .eMesh format.

[0049] Specifically, to avoid smoothing out the "boundaries between regions of different roughness and the lines connecting the surface to the terrain" during mesh generation, the feature edges of the surface need to be extracted separately: Add edge extraction rules for roughness surfaces in the system / surfaceFeatureExtractDict configuration file. roughnessPatchSurface { file "constant / triSurface / roughnessPatch_aligned.stl"; edges { boundaryEdges { type boundary;} / / Extracts the boundary edges of regions of a facet. featureEdges { type angle; angle 15;} / / Extract the steep edges where the patch meets the terrain. } writeObj true; } The terminal runs the surfaceFeatureExtract command to generate a feature edge file (.eMesh format) for the surface. This file will be recognized by snappyHexMesh as an "unsimplifiable feature" to ensure that the surface outline is accurately preserved.

[0050] 423) Run blockMesh to build the background mesh, then run snappyHexMesh to generate associated meshes and boundary patches. In the splitting stage, split the mesh near the facets at the specified level. In the bonding stage, ensure that the surface mesh bonding error is ≤1 meter through iteration.

[0051] Specifically, in this embodiment, the steps for generating associated meshes and boundary patches using snappyHexMesh are as follows.

[0052] First, construct a structured background mesh using blockMesh; The terminal runs `snappyHexMesh -overwrite`, and the tool performs the "mesh splitting → meshing → boundary layer generation" process according to its configuration. The core operations for roughness patches are as follows: Splitting phase: Split the mesh near the facet by refineLevel 4 (in conjunction with terrain refinement). Adhesion stage: 20 iterations are used to ensure that the mesh adhesion error on the surface of the patch is ≤1m; Boundary layer stage: Reuse the "3-layer boundary layer on terrain surface" setting (first layer thickness 0.5m, growth factor 1.2), no need to redefine it.

[0053] 424) After the mesh is generated, use the checkMesh tool to verify whether the type of the roughness-specific patch is "wall" and to verify the integrity of the patch association.

[0054] Specifically, after the mesh is generated, run checkMesh in the terminal to verify two key aspects: Boundary patch type: In the "Boundary Patches" section, confirm that the roughnessPatch patch type is "wall". If it is not "wall", force correction using the boundaryTagging dictionary of snappyHexMeshDict. Patch integrity: Verify that the number and area of ​​the roughnessPatch are consistent with the geometric parameters of the roughness patch (ensuring no patches are missing or redundant).

[0055] 43) Boundary condition mapping for roughness parameters: Modify the turbulenceProperties file to enable the wall roughness correction function of the SST k-ω model; in the boundary condition dictionary, set the velocity field on the roughnessPatch to no-slip condition, set the boundary condition type for turbulent kinetic energy and specific dissipation rate to wallFunction, and specify the use of kqRWallFunction and omegaWallFunction.

[0056] Specifically, in this embodiment, the method steps for mapping the boundary conditions of the roughness parameters are as follows.

[0057] 431) Modify the turbulenceProperties file: Open the constant / turbulenceProperties file, add a "roughness influence switch" in the turbulence model configuration section, and specify that the wall roughness correction function must be enabled for the SST k-ω model to ensure that the model can read subsequent assignments. The value is then applied to turbulent field calculations.

[0058] 432) Configure boundary condition dictionaries (U, k, ω): Enter the 0 / folder (the directory where initial boundary conditions are stored), and modify the boundary condition dictionaries for velocity field (U), turbulent kinetic energy (k), and specific dissipation rate (ω) respectively: For the roughnessPatch, the velocity field (U) is set to "noSlip" to reflect the actual obstruction of airflow by the ground surface. The boundary condition type for turbulent kinetic energy (k) and specific dissipation rate (ω) is set to "wallFunction", and kqRWallFunction (for k) and omegaWallFunction (for ω) are specified. These two types of wall functions support reading roughness parameters. For the follow-up The assignment is reserved for the interface.

[0059] 44) Batch assignment of roughness parameters based on patch ID: Develop an OpenFOAM auxiliary tool or write a Python script to read "patch ID - roughness length". "Associate files, establish topological associations between triangular face IDs and mesh faces in the roughnessPatch, and..." Values ​​are written in batches to the 0 / k and 0 / omega dictionaries to assign values ​​to the roughnessLength parameter, and the correctness of the assignment results is verified using the paraView tool.

[0060] This step implements the functionality described in "roughnessProperties.txt". The precise matching of "value" and "mesh boundary patch" avoids the inefficiency and error of manual assignment. The method and steps are as follows.

[0061] 441) Develop a coarse length mapping tool: Develop an OpenFOAM auxiliary tool or write a Python script with the core function of reading the roughnessProperties.txt file generated in section 8.1 and extracting the "patch" data. Establish the mapping relationship; through OpenFOAM's mesh interface (such as the libmeshTools library), establish a topological association between the "triangle patch ID" and the "mesh face in the roughnessPatch"; based on the association relationship, [the following steps are taken]. Values ​​are written in batches to the 0 / k and 0 / omega dictionaries to assign values ​​to the roughnessLength parameter of the roughnessPatch patch.

[0062] 442) Verification of assignment results: Launch paraView (the visualization tool accompanying OpenFOAM), load the computational mesh and 0 / k, 0 / omega boundary condition files, and view the roughnessLength distribution of the roughnessPatch using the "Boundary Condition Visualization" function: It is necessary to confirm the equivalent roughness regions of water bodies, grasslands, and small buildings. The values ​​are consistent with the design values, with no errors or omissions, ensuring that the roughness distribution conforms to physical reality.

[0063] 45) Roughness activation mechanism during simulation: During CFD simulation, the SST k-ω turbulence model calls the kqRWallFunction and omegaWallFunction wall functions to read the roughnessLength parameter of the roughnessPatch in real time (i.e., (value), by introducing The formula for calculating the friction velocity of the correction term corrects the turbulence field and wind speed field in the near-wall region, so that the simulated near-surface wind speed profile and turbulence intensity distribution match the actual surface roughness characteristics. For rough walls, friction speed The calculation will introduce Correction term: The formula for calculating friction speed is as follows: in: The friction speed; Kármán's constant; Wind speed near the wall; To calculate the height of the point.

[0064] The corrected friction velocity will further affect the generation of turbulent kinetic energy k and the dissipation process of specific dissipation rate ω, ultimately making the simulated near-surface wind speed profile and turbulence intensity distribution match the actual surface roughness characteristics, ensuring that the simulation deviation is ≤5%.

[0065] In this embodiment, the CFD simulation includes setting the computational domain, mesh generation, and computational settings. Specifically, a "roughnessPatch" is defined in snappyHexMeshDict, a Python script is written to read roughnessProperties.txt, and the roughness values ​​are batch-assigned to the 0 / k and 0 / omega files (using kqRWallFunction and omegaWallFunction); the pimpleFoam solver and the SST k-ω model (with wall roughness correction enabled) are selected, and the measured fitted wind profile information is applied at the inlet.

[0066] (1) Calculation domain settings: The blockage rate is less than 3%. The plane size of the calculation domain is extended along the downwind direction by no less than 5 times the longest side of the target area along the downwind direction and along the crosswind direction by no less than 3 times the longest side of the target area along the crosswind direction. The minimum height is 50 meters below the lowest point of the terrain and the maximum height is no less than 1500 meters or 10 times the height of the tallest building.

[0067] Specifically, the computational domain is set to meet the requirement of a blockage rate of <3% and to meet the requirement of a fully developed atmospheric boundary layer, including the following:

[0068] ① Setting the dimensions of the computational domain plane.

[0069] Based on the bounding rectangle of the target study area, the extension length along the downwind direction (ABL wind incidence direction) and the crosswind direction (perpendicular to the ABL wind incidence direction) shall not be less than 6 times the longest side of the target area along the downwind direction and 4 times the longest side of the crosswind direction, respectively. The target area refers to the core area of ​​complex terrain and building / structure distribution that needs to be simulated. The transition section has been integrated into the boundary of the computational domain in step three to ensure that the incoming wind has no boundary disturbance before entering the core area.

[0070] ② Set the height of the computing domain.

[0071] The height direction must completely cover the atmospheric boundary layer (ABL) and the tallest building / structure in the target area. The specific values ​​are: minimum height: set to 50m below the lowest terrain elevation in the target area (to avoid missing grids in terrain depressions); maximum height: 2000m.

[0072] ③ The total projected area of ​​all buildings and terrain within the computational domain shall be ≤3% of the total planar area of ​​the computational domain (blocking rate), and the ratio of the projected width of a single building in the downwind direction to the downwind length of the computational domain shall be ≤5%.

[0073] (2) Grid division.

[0074] A background mesh and local encryption strategy is adopted. A structured background mesh is constructed using blockMesh, and adaptive encryption around buildings and terrain is achieved by combining snappyHexMesh. The resolution of the background mesh is 1 / 300 of the size of the computational domain plane. It adopts a non-uniform distribution in the vertical direction. The mesh size is 5 meters in the 0-50 meter range near the ground, the mesh size increases by 1.1 times in the 50-500 meter range, and remains at 15 meters above 500 meters.

[0075] The OpenFOAM simulation process uses a layered strategy of "background mesh + local refinement" for mesh generation. A structured background mesh is constructed using blockMesh, and adaptive refinement of the area around buildings and terrain is achieved by combining snappyHexMesh. The specific operation is as follows.

[0076] ① Background mesh (blockMesh) parameter configuration.

[0077] Base resolution: A uniform base grid is set according to the planar dimensions of the computational domain. The grid size in the windward direction (x-axis) and crosswind direction (y-axis) is 1 / 300 of the length of the corresponding direction of the computational domain. A non-uniform distribution is used in the vertical direction (z-axis): near the ground (0-50m): grid size 5m (to capture the near-wall wind field affected by roughness); 50-500m: grid size increases linearly at a ratio of 1.1 (from 5m to 15m); above 500m: a constant size of 15m is maintained (to adapt to the characteristics of the upper atmospheric boundary layer). Grid partitioning: The computational domain is divided into 3×3×3 sub-blocks using blockMesh. The grid lines are uniformly distributed within each sub-block, and a continuous transition mode is used between sub-blocks to avoid abrupt grid changes. Simultaneously, it is ensured that the background grid coordinate system is completely aligned with the terrain, building, and roughness patch coordinate systems, with an error controlled within 0.1m.

[0078] ② Feature edge extraction (surfaceFeatureExtract).

[0079] For the terrain STL file and building STL file merged in step eight, extract two types of key geometric features: Building outline: Set the edge angle threshold to 15°, and extract the intersection of the building facade and the roof, and the adjacent edges with a spacing of ≤2m between buildings; Topographic steep slope lines: Extract topographic abrupt change lines with a slope greater than 30°, generate feature line files (.eMesh), set the "immutable" attribute to the extracted feature lines (to ensure that they are not simplified or deleted during mesh generation), and store them in the constant / triSurface / directory for snappyHexMesh to recognize in order to preserve geometric accuracy.

[0080] ③ Adaptive encryption (snappyHexMesh) parameter settings The encryption parameters are configured according to the principle of "building priority, terrain coordination", while reserving space for encryption level matching for roughness patch association: 1) Castellation (mesh splitting): Maximum splitting level: 6-level splitting is set near the building surface (base grid → 167m / 2) 6 ≈2.6m), the terrain surface is divided into 4 levels (167m / 2). 4 ≈10.4m), areas far from terrain and buildings remain at level 0 (no splitting); Split trigger conditions: When the distance between the grid and the building surface is less than 5m and the distance between the grid and the terrain surface is less than 10m, the corresponding level of splitting will be automatically initiated.

[0081] 2) Snapping (mesh bonding): Fitting accuracy: Fitting error on building surfaces ≤ 0.5m, fitting error on terrain surfaces ≤ 1m, achieved through 20 iterations; Non-orthogonality control: After each iteration, check the non-orthogonal angle of the mesh. If it is greater than 65°, automatically adjust the node position to ensure that the non-orthogonal angle of the entire domain is ≤70°.

[0082] 3) Layering (boundary layer mesh): Building wall: Set 5 boundary layers, the first layer is 0.2m thick, the growth factor is 1.2, the total thickness is 1.19m, and ensure that the y⁺ value is in the range of 30-100 (adapting to the wall function of the SST k-ω model). Terrain surface: Set 3 boundary layers, the first layer is 0.5m thick, the growth factor is 1.2, and the total thickness is 1.82m (the roughness patch coverage area reuses this boundary layer setting and does not need to be defined separately).

[0083] ④ Reservation of interfaces for roughness patch association.

[0084] In the "surfaces" dictionary of snappyHexMeshDict, an entry point for associating roughness patches is reserved: Ensure that the density level (level 4) and boundary layer settings (3 layers) of the terrain surface are consistent with the parameters of the roughness patch; The "faceZone" function preserves the topological correspondence between terrain grid surfaces and roughness patches, providing a data interface for batch association of "patch ID - grid surface".

[0085] ⑤ Mesh quality verification standards.

[0086] Quantitative metrics: Verified using the checkMesh tool to ensure that the following is met: minimum mesh volume > 1×10⁻ 6 m 3 (No negative volume); Volumetric distortion rate < 0.9; Aspect Ratio < 5 (near ground grid), < 10 (high-altitude grid); Visual verification: Use paraView to check the grid continuity at building corners and steep terrain edges to confirm that there is no overlap or gap; at the same time, check the smooth transition between the grid and the terrain grid in the area covered by the roughness patch.

[0087] (5) Calculation settings: The pimpleFoam solver was selected, the governing equation was the three-dimensional Reynolds-averaged Navier-Stokes equation, the turbulence model was the SST k-ω model, roughness correction was enabled, and the boundary conditions included inlet wind speed profile, outlet free outflow, top slip wall and crosswind periodic boundary.

[0088] After completing the computational domain setting and mesh generation, the computational settings are achieved through solver configuration, turbulence model selection, boundary condition definition, and convergence control to ensure that the roughness effect is accurately transferred to the flow field calculation, as detailed below.

[0089] ① Selection of solver and governing equations.

[0090] The pimpleFoam solver is used, based on the pressure-velocity coupled PIMPLE algorithm (PISO+SIMPLE hybrid algorithm). The governing equations are the three-dimensional Reynolds-averaged Navier-Stokes (RANS) equations, combined with the continuity equation. The momentum equation incorporates a gravity term and Coriolis force corrections, and the density is approximated using the Boussinesq approximation. ).

[0091] ② Turbulence model configuration.

[0092] The SST k-ω turbulence model with wall roughness correction is selected: The model's basic coefficients are set to the default settings. , , ), ensuring compatibility with the OpenFOAM solver; Roughness correction enabled: Add a "roughnessCorrection on" switch to the constant / turbulenceProperties file to explicitly require the model to read the roughnessLength parameter (i.e., ...). (value), providing a basis for correction in near-wall turbulence calculations; Initial turbulence parameters: turbulent kinetic energy k and specific dissipation rate ω obtained by fitting measured data. For the specific fitting method of the measured data in this example, please refer to [reference needed]. Figure 7 The fitting results are referenced. Figure 8 The error between the fitted data and the measured data is less than 3%, which meets the accuracy requirements.

[0093] ③ Boundary condition setting Define boundary conditions according to "functional zones" to ensure that the physical characteristics of each zone match: a) Inlet boundary (ablInlet, corresponding to the alongwind extension of the computational domain): Velocity (U): Wind speed profile fitted using measured data Pressure (p): Set to zero gradient; Turbulent kinetic energy (k) and specific dissipation rate (ω): both are vertical profiles fitted from measured data.

[0094] b) Outlet (corresponding to the downwind end of the computational domain): Velocity (U): Free outflow condition (inletOutlet); Pressure (p): Fixed value 0 Pa (relative pressure).

[0095] c) Roughness patch boundary: Velocity (U): No slip condition (noSlip); Pressure (p): zero gradient; Turbulence parameters (k, ω): The roughnessLength parameter is automatically read using two types of wall functions, kqRWallFunction (for k) and omegaWallFunction (for ω), and is calculated using the formula. By modifying the friction velocity, the influence of roughness on the turbulent field can be realized.

[0096] c) Top boundary (corresponding to the highest height of the computational domain): set as a slip wall with zero gradients for both velocity and pressure.

[0097] d) Crosswind boundary (side1 / side2): A cyclic boundary is used.

[0098] ④ Time and Iteration Control The system is configured in two phases: "initialization + effective simulation," to ensure that the flow field develops fully and the data is reliable. Calculation type: Unsteady-state simulation is used, with an initial time step of deltaT=0.5s, which is automatically adjusted by adjustTimeStep to ensure that the number of CFLs is ≤1.0 (maximum deltaT=2s). The total computation time for this simulation is 3600s (1 hour), of which the first 1800s is the flow field initialization stage (not included in the statistical results), and the last 1800s is the effective simulation stage, used for wind field characteristic analysis. The iteration settings for this example are as follows: the PIMPLE algorithm iterates 5 times in each time step, and the momentum equation and pressure equation are coupled and iterated 3 times.

[0099] ⑤ Numerical Algorithm Configuration Pressure equation: The GAMG (Geometric Multigrid) solver is used, with a residual convergence threshold of 1×10⁻. 6 Smoothing iteration count 2; Momentum equation (U): Solved using PBiCGStab, residual threshold 1×10⁻ 5 The preconditioner is DILU; Turbulence equations (k, ω): solved using the PBiCGStab solver, residual threshold 1×10⁻ 5 ; Discrete scheme: Spatial discretization uses the Gaussian linear scheme (central difference), with a limitedLinearV limiter (limiting coefficient 1.0) added to the convection term, and the time discretization uses a second-order implicit scheme (backward).

[0100] ⑥ Convergence verification and result output The reliability of the results is controlled using a dual standard of "quantitative convergence + physical verification": Convergence criterion: The root mean square error (RMS) of global velocity and pressure changes ≤ 1 × 10⁻ over 500 consecutive time steps. 4 And the difference in mass flow rate between the inlet and outlet is ≤1%.

[0101] Step 5: Result Output and Verification Output simulation results to verify the deviation between near-surface wind field information and measured wind field information, and ensure that the deviation does not exceed the preset threshold.

[0102] In this embodiment, the wind field data (U, p, k, ω) is saved every 60 seconds, and the wind field data of the last 600 seconds is saved and output to verify the comparison between the near-surface wind field information and the measured data, with a deviation of ≤5%.

[0103] This embodiment also proposes an automated modeling system for complex terrain and buildings / structures that considers surface roughness, used to execute the method described above in this embodiment, including a data acquisition module, a terrain processing module, a building processing module, a filtering module, a roughness processing module, a fusion module, a CFD simulation module, and a verification module.

[0104] The data acquisition module is used to collect topographic elevation data, urban building data, ground roughness category data, and measured wind field data from open-source databases. Specifically, the topographic elevation data is in HGT format; the urban building data is in GeoJSON format; the ground roughness category data is in GeoTiff format; and the measured wind field data is collected through a wind measurement tower.

[0105] The terrain processing module uses Python or Matlab scripts to convert terrain elevation data into point cloud data, and uses an improved Witoszynski wind tunnel contraction curve to generate transition sections, and generates STL files of the terrain through Delaunay triangulation.

[0106] The building processing module uses the Ramer-Douglas-Peucker algorithm and standard deviation method to process building / structure data, remove redundant and outlier points, and generate STL files of buildings / structures using the solidPython2 library and OpenSCAD.

[0107] The filtering module is used to filter small buildings / structures and calculates the equivalent roughness based on the Lettau formula.

[0108] The roughness processing module uses Python's gdal library to process ground roughness category data, assigns terrain roughness based on Davenport-Wieringa classification, and generates roughness patches by combining equivalent roughness, outputting a roughness attribute file.

[0109] The fusion module uses OpenSCAD's union() function to merge terrain STL, building STL, and roughness patches, and performs Boolean operations for verification and watertightness checks.

[0110] The CFD simulation module uses OpenFOAM for wind field simulation, including setting up the computational domain, mesh generation, and computational settings. Roughness patches are defined using snappyHexMeshDict, and a Python script is written to assign roughness values ​​to the boundary condition file. The roughness length is achieved by mapping the roughness patches to boundary conditions and associating them with the mesh. Bind to the computational grid and perform wind field calculations.

[0111] The verification module is used to output simulation results, verify the deviation between the near-surface wind field information and the measured wind field information, and ensure that the deviation does not exceed the preset threshold.

[0112] In summary, this embodiment relies on open-source toolchains such as Python and OpenSCAD to build a fully automated workflow from multi-source data processing to model generation: it automatically converts NASA SRTM terrain data (HGT format) to ASCII point clouds and extracts and simplifies OpenStreetMap building data (GeoJSON format) through scripts, generating terrain STL and building STL in batches without manual intervention; for multi-region parallel modeling scenarios, it supports the simultaneous processing of modeling tasks for ≥3 target regions, improving modeling efficiency by more than 60% compared to traditional manual operations (such as manual drawing in ArcGIS + manual repair in ANSYS SpaceClaim), and significantly shortening the preprocessing cycle for wind field simulation (from 15 days to within 3 days).

[0113] This embodiment innovatively adopts a dual-path automated processing mechanism of "inherent terrain roughness + equivalent roughness of small buildings": based on the Davenport-Wieringa classification system, it automatically assigns roughness lengths to GeoTiff land use data using Python's gdal library. (e.g., 0.0002m for water areas, 0.03m for grassland areas); small buildings (area < 25m²) will be filtered using the Lettau formula. 2 (or height < 5m) converted to equivalent Generate a patch "Associated roughness patches; subsequent implementation using OpenFOAM's wall functions (kqRWallFunction, omegaWallFunction)." Precise binding with the computational grid ensures that the simulation results of near-surface wind speed profile and turbulence intensity deviate from the measured values ​​by ≤5%, improving the simulation accuracy by 40% compared to the traditional "uniform roughness assignment" method.

[0114] In this embodiment, addressing the pain point of traditional modeling where "modeling small buildings one by one leads to a surge in mesh size and excessive computational cost," this invention replaces geometric modeling of individual small buildings by filtering them and calculating equivalent roughness. On the one hand, it reduces the number of STL files that need to be generated (more than 60% in dense areas), avoiding Boolean fusion failures and mesh distortion caused by an excessive number of buildings. On the other hand, the introduction of equivalent roughness eliminates the need for high-precision mesh refinement around small buildings, reducing the total amount of mesh in the computational domain by 30% to 50%. This shortens mesh generation time (from 24 hours to 8 hours) and reduces the hardware computing power requirements for CFD simulation (complex area simulation can be completed with ordinary workstations).

[0115] In this embodiment, addressing the common issues of measurement noise, outliers, and missing data in satellite images and LiDAR data, this invention employs a dual data cleaning mechanism to ensure modeling accuracy: The Ramer-Douglas-Peucker (RDP) algorithm is used to simplify building outline coordinates, with ε set to 0.01~0.05 times the building's outer perimeter, effectively removing redundant points from LiDAR scans (noise removal rate ≥85%); outliers in building height are identified and removed using the standard deviation method (z-core value > ±3), replacing them with the average height of buildings within a 300m radius, avoiding STL geometric distortion caused by outliers; even in scenarios with a missing rate of ≤10% in the original data, terrain elevation data can still be completed through Kriging interpolation, ensuring the reliability of the modeling results and overcoming the limitation of traditional methods where "data quality determines modeling accuracy."

[0116] In this embodiment, the entire development process is based on open-source tools (Python, Matlab, OpenSCAD, OpenFOAM) and public databases (NASA SRTM, OpenStreetMap, Tsinghua University Nebula Data Service Platform), without relying on commercial software such as ArcGIS and ANSYS, which can save software procurement costs (reducing costs by 50,000 to 100,000 yuan per project). At the same time, complete script code and parameter configuration files are provided, so that technicians do not need to have complex programming or CFD professional backgrounds to reproduce the modeling and simulation process by following the steps in the instruction manual. This significantly reduces the technical application threshold for small teams and research institutions in the field of wind engineering, and helps to popularize and promote the simulation technology of wind fields in complex terrain.

[0117] Therefore, this embodiment has the following advantages: (1) It can automatically and in batches model terrain and buildings / structures; (2) Automation takes into account the surface roughness, making the surface wind field simulation more accurate.

[0118] (3) The simplified small buildings are represented by equivalent roughness to reflect their impact on the surface wind field. Compared with traditional modeling methods, this can shorten the modeling time and reduce the modeling complexity. Previous modeling methods heavily relied on the quality of satellite imagery and LiDAR data. However, real-world datasets contain a large amount of measurement errors, faulty data, and even missing data, severely impacting modeling accuracy. This method uses the RDP algorithm combined with standard deviation to remove redundant points (reducing LiDAR measurement noise) and delete outliers.

[0119] In this article, "buildings / structures" includes buildings or structures.

[0120] The above-described embodiments are merely preferred embodiments provided to fully illustrate the present invention, and the scope of protection of the present invention is not limited thereto. Equivalent substitutions or modifications made by those skilled in the art based on the present invention are all within the scope of protection of the present invention. The scope of protection of the present invention is defined by the claims.

Claims

1. An automated modeling method for complex terrain and buildings / structures considering surface roughness, characterized in that: Includes the following steps: Step 1: Data Collection Collect topographic elevation data, urban building data, ground roughness category data, and measured wind field data; the topographic elevation data is in HGT format. Urban building data is in GeoJSON format; ground roughness category data is in GeoTIFF format; wind field measurement data is collected through a wind measurement tower. Step Two: Data Processing The collected terrain elevation data, urban building data, and ground roughness category data were processed separately, including: The method for processing terrain elevation data is as follows: use Python or Matlab scripts to convert binary terrain elevation data into point cloud data in ASCII format, use an improved Witoszynski wind tunnel contraction curve to generate transition sections around the terrain, and convert the point cloud data into triangular mesh data through Delaunay triangulation to generate STL files of the terrain. The processing method for urban building data is as follows: the Ramer-Douglas-Peucker algorithm is used to simplify the outline coordinates of buildings / structures, remove redundant points, and combine the standard deviation method to delete outlier points; small buildings / structures are filtered, the windward area and average height of the filtered buildings are recorded, and the equivalent roughness is calculated based on the Lettau formula; the shape coordinates of buildings / structures are extracted using Python scripts, OpenSCAD code is generated using the solidPython2 library, OpenSCAD is imported and negatively stretched to generate the STL file of the building / structure; The method for processing ground roughness category data is as follows: call Python's gdal library to read ground roughness category data in GeoTiff format, establish a mapping between land use type and roughness length based on Davenport classification and Wieringa correction, assign roughness to the terrain, calculate the equivalent roughness based on Lettau formula, generate roughness patches through Delaunay triangulation, and output a roughness attribute file with patch ID and roughness value association. Step 3: Data Fusion Use OpenSCAD's union() function to merge terrain STL files, building / structure STL files, and roughness patches, enabling Boolean operation verification and watertightness checks; Step 4: CFD Simulation Computational fluid dynamics simulations were performed using OpenFOAM, including wind field calculations. This included setting up the computational domain, mesh generation, and computational settings. Roughness patches were defined in snappyHexMeshDict, and a Python script was written to batch assign roughness values ​​to the boundary condition file. The pimpleFoam solver and the SST k-ω model were used with wall roughness correction enabled. Measured wind profile information was applied at the inlet. The roughness length was achieved by mapping the roughness patches to boundary conditions and associating them with the mesh. Binding to the computational grid and performing wind field calculations; Step 5: Result Output and Verification Output simulation results to verify the deviation between near-surface wind field information and measured wind field information, and ensure that the deviation does not exceed the preset threshold.

2. The automated modeling method for complex terrain and buildings / structures considering surface roughness according to claim 1, characterized in that: In step two, the formula for the transition curve generated based on the improved Witoszynski wind tunnel contraction curve is as follows: in: and These represent the horizontal and vertical coordinates of the transition curve, respectively. Indicates the height of the transition curve; Indicates the length of the transition curve; This represents the constant coefficient.

3. The automated modeling method for complex terrain and buildings / structures considering surface roughness according to claim 1, characterized in that: In step two, the epsilon value of the Ramer-Douglas-Peucker algorithm is set to 0.01-0.05 times the perimeter of the building / structure outline to control the degree of simplification; the standard deviation method uses a z-core value greater than ±3 as the criterion for outlier detection and replaces it with the average height of buildings within a 300-meter radius.

4. The automated modeling method for complex terrain and buildings / structures considering surface roughness according to claim 1, characterized in that: In step two, the criteria for filtering small buildings / structures are that the roof area is less than a first set threshold or the height is less than a second set threshold; the equivalent roughness is calculated using the Lettau formula: in: This is the ratio of the projected area of ​​the windward side of a building / structure to the horizontal area of ​​the target region. This represents the average height of all filtered buildings / structures within the corresponding building area; and: in: This indicates the projected area of ​​the windward side of a building / structure.

5. The automated modeling method for complex terrain and buildings / structures considering surface roughness according to claim 1, characterized in that: In step two, the negative stretching extends downwards to z=-5 meters to ensure that the building / structure can penetrate the terrain and achieve a seamless connection with the terrain; the OpenSCAD code is automatically generated by the solidPython2 library and multiple buildings / structures are processed in batches.

6. The automated modeling method for complex terrain and buildings / structures considering surface roughness according to claim 1, characterized in that: In step two, the roughness length mapping of the Davenport classification corrected by Wieringa includes: a roughness length of 0.0002 meters for water bodies, a roughness length of 0.005 meters for bare land, and a roughness length of 0.03 meters for grassland; the area of ​​the triangular mesh of the roughness patch is no more than 100 square meters, and the normal direction is uniform.

7. The automated modeling method for complex terrain and buildings / structures considering surface roughness according to claim 1, characterized in that: In step three, self-intersection is avoided by using Boolean operations in OpenSCAD, and water tightness is checked to ensure no leakage surfaces; the fused model is used for preprocessing in OpenFOAM simulation.

8. The automated modeling method for complex terrain and buildings / structures considering surface roughness according to claim 1, characterized in that: In step four, the CFD simulation includes: Computational domain settings: The blockage rate is less than 3%. The planar size of the computational domain is extended along the downwind direction by no less than 5 times the longest side of the target area along the downwind direction, and along the crosswind direction by no less than 3 times the longest side of the target area along the crosswind direction. The minimum height is 50 meters below the lowest point of the terrain, and the maximum height is no less than 1500 meters or 10 times the height of the tallest building. Mesh generation: A background mesh and local densification strategy is adopted. A structured background mesh is constructed using blockMesh, and adaptive densification around buildings and terrain is achieved by combining snappyHexMesh. The resolution of the background mesh is 1 / 300 of the size of the computational domain plane. The vertical direction adopts a non-uniform distribution. The mesh size is 5 meters in the 0-50 meter range near the ground, the mesh size increases by 1.1 times in the 50-500 meter range, and remains at 15 meters above 500 meters. Calculation settings: The pimpleFoam solver was selected, the governing equations were the three-dimensional Reynolds-averaged Navier-Stokes equations, the turbulence model was the SST k-ω model, roughness correction was enabled, and the boundary conditions included inlet wind speed profile, outlet free outflow, top slip wall and crosswind periodic boundary.

9. The automated modeling method for complex terrain and buildings / structures considering surface roughness according to claim 1, characterized in that: In step four, the roughness patch is mapped by boundary conditions and associated with the mesh to achieve the roughness length. The steps for binding with the computational grid are as follows: 41) Surface Data Format Conversion and Import Standardization: Convert triangular mesh surfaces containing roughness categories to OpenFOAM-compatible .stl format, and unify the coordinate system using the surfaceTransformPoints tool; write a Python script to read the geometric data of the triangular mesh surfaces and their corresponding coordinates. Value, generate "patch ID - roughness length" "Associated files; 42) Mesh and Patch Association Configuration: By configuring snappyHexMeshDict to work collaboratively with tools, ensure that the mesh fits the patch outline and generates a dedicated boundary patch, specifically: 421) In snappyHexMeshDict, define the roughness patch as a feature surface, specifying its surface type, file path, corresponding boundary patch information, and mesh refinement parameters; 422) Extract the feature edges of the roughness patch using the surfaceFeatureExtract tool and generate a feature edge file in .eMesh format; 423) Run blockMesh to build the background mesh, then run snappyHexMesh to generate associated meshes and boundary patches. In the splitting stage, split the mesh near the patch at the specified level. In the bonding stage, ensure that the surface mesh bonding error is ≤1 meter through iteration. 424) After the mesh is generated, use the checkMesh tool to verify whether the type of the roughness-specific patch is "wall" and to verify the integrity of the patch association; 43) Boundary condition mapping for roughness parameters: Modify the turbulenceProperties file to enable the wall roughness correction function of the SST k-ω model; in the boundary condition dictionary, set the velocity field on the roughnessPatch to no-slip condition, set the boundary condition type for turbulent kinetic energy and specific dissipation rate to wallFunction, and specify the use of kqRWallFunction and omegaWallFunction. 44) Batch assignment of roughness parameters based on patch ID: Develop an OpenFOAM auxiliary tool or write a Python script to read "patch ID - roughness length". "Associate files, establish topological associations between triangle face IDs and mesh faces in the roughnessPatch, and..." Values ​​are written in batches to the 0 / k and 0 / omega dictionaries to assign values ​​to the roughnessLength parameter, and the correctness of the assignment results is verified by the paraView tool; 45) Roughness Activation Mechanism During Simulation: During CFD simulation, the SST k-ω turbulence model calls the kqRWallFunction and omegaWallFunction wall functions to read the roughnessLength parameter of the roughnessPatch in real time, and introduces... The formula for calculating the friction velocity of the correction term corrects the turbulence field and wind speed field in the near-wall region, so that the simulated near-surface wind speed profile and turbulence intensity distribution match the actual surface roughness characteristics. The formula for calculating friction speed is as follows: in: The friction speed; Kármán's constant; Wind speed near the wall; To calculate the height of the point.

10. An automated modeling system for complex terrain and buildings / structures that takes into account surface roughness, for performing the method as described in any one of claims 1 to 11, characterized in that: include: The data acquisition module is used to collect terrain elevation data, urban building data, ground roughness category data, and measured wind field data from open-source databases. Specifically, the terrain elevation data is in HGT format; the urban building data is in GeoJSON format; the ground roughness category data is in GeoTiff format; and the measured wind field data is collected through a wind measurement tower. The terrain processing module uses Python or Matlab scripts to convert terrain elevation data into point cloud data, and uses an improved Witoszynski wind tunnel contraction curve to generate transition sections, and generates STL files of terrain through Delaunay triangulation. The building processing module uses the Ramer-Douglas-Peucker algorithm and standard deviation method to process building / structure data, remove redundant and outlier points, and generate STL files of buildings / structures using the solidPython2 library and OpenSCAD. A filtering module is used to filter small buildings / structures and calculate the equivalent roughness based on the Lettau formula; The roughness processing module uses Python's gdal library to process ground roughness category data, assigns terrain roughness based on Davenport-Wieringa classification, generates roughness patches by combining equivalent roughness, and outputs a roughness attribute file. The fusion module uses OpenSCAD's union() function to merge terrain STL, building STL, and roughness patches, and performs Boolean operations for verification and watertightness checks. The CFD simulation module uses OpenFOAM for wind field simulation, including setting the computational domain, mesh generation, and computational parameters. Roughness patches are defined using snappyHexMeshDict, and a Python script is written to assign roughness values ​​to the boundary condition file. The roughness length is then achieved by mapping the roughness patches to boundary conditions and associating them with the mesh. Binding to the computational grid and performing wind field calculations; The verification module is used to output simulation results, verify the deviation between the near-surface wind field information and the measured wind field information, and ensure that the deviation does not exceed the preset threshold.

Images