A System Coupled GIS and CFD for Atmospheric Pollution Dispersion Simulation in Urban Blocks

Abstract

Atmospheric pollution is a critical issue in public health systems. The simulation of atmospheric pollution dispersion in urban blocks, using CFD, faces several challenges, including the complexity and inefficiency of existing CFD software, time-consuming construction of CFD urban block geometry, and limited visualization and analysis capabilities of simulation outputs. To address these challenges, we have developed a prototype system that couples 3DGIS and CFD for simulating, visualizing, and analyzing atmospheric pollution dispersion. Specifically, a parallel algorithm for coordinate transformation was designed, and the relevant commands were encapsulated to automate the construction of geometry and meshing required for CFD simulations of urban blocks. Additionally, the Fluent-based command flow was parameterized and encapsulated, enabling the automatic generation of model calculation command flow files to simulate atmospheric pollution dispersion. Moreover, multi-angle spatial partitioning and spatiotemporal multidimensional visualization analysis were introduced to achieve an intuitive expression and analysis of CFD simulation results. The result shows that the constructed geometry is correct, and the mesh quality meets requirements with all values above 0.45. CPU and GPU parallel algorithms are 13.3× and 25× faster than serial. Furthermore, our case study demonstrates the developed system’s effectiveness in simulating, visualizing, and analyzing atmospheric pollution dispersion in urban blocks.

1. Introduction

Atmospheric pollution has long been a paramount concern that it can cause serious harm to human health and lead to substantial economic losses. The simulation of atmospheric pollution dispersion is necessary for the loss reduction. Computational Fluid Dynamics (CFD) models are powerful computational and modeling tools that have been widely used in simulation environments in many different fields, with the most common application being the simulation of atmospheric dispersion of hazardous emissions at an urban scale. For example, some studies have conducted simulation experiments on gas leakage using CFD and have shown that an increase in ambient temperature leads to an increase in the diffusion rate of the leaked gas in both vertical and horizontal directions [1,2]. The study by Cormier et al. [3] demonstrates how CFD can predict the consequences of liquefied natural gas (LNG) releases and indicates that CFD is an effective tool that can enhance the safety of LNG facilities. Tauseef et al. [4] and Yoshie et al. [5] simulated gas dispersion after natural gas leakage incidents using CFD models and studied the effects of obstacles, such as vegetation, buildings, and roads, on gas dispersion in the environment. Siddiqui et al. [6] used a CFD model to simulate gas leakage in indoor pipelines, and they delineated the hazardous zone indoors based on the simulation results. To study gas diffusion in soil caused by leakage accidents in urban buried natural gas pipelines, some researchers proposed suggestions for the proper utilization of soil near the pipelines by using fluid mechanics theory, porous media theory, and CFD models [7,8]. Using a CFD model, Feißel et al. [9] examined how vehicle emissions affect air quality and suggested a simplified prediction method. Schalau et al. [10] developed and validated a modified k-ε-turbulence model in OpenFOAM v5.0 for heavy gas dispersion in built-up areas under atmospheric conditions, and the result demonstrates that the model is reliable in prediction. Comparing CFD simulations of pollutant dispersion with Gaussian-type models, some researchers have found that CFD is better suited for simulating pollutant dispersion in areas with high building density [11,12]. Artificial intelligence and data-driven models are gaining attention in various applications due to the availability of large data. Hybrid models, such as LSTM-ALO, ANFIS-GBO, ELM-PSOGWO, SVR-SAMOA, and ANN-EMPA, have been proposed and applied to agriculture and environmental prediction [13,14,15,16,17]. Nowadays, deep learning models, such as the Short-Term Memory network (LSTM), the Convolutional Neural Network (CNN), and the Recurrent Neural Network (RNN), have found wide application in predicting atmospheric pollutants across different spatial and temporal scales [18,19,20,21,22,23,24,25,26]. Some researchers have combined LSTM and CNN models, considering the spatial and temporal correlations in data, to improve the performance of the prediction model [27,28].

Fluent software, as the most widely used CFD package, has accumulated various physical models and developed more efficient numerical techniques. However, due to the inclusion of numerous models and parameters to meet the demands of different fields, Fluent software is complex and has a high learning curve. To address this issue, researchers have explored custom development based on Fluent software and achieved certain results in various fields. For example, by studying the command flow log files of Gambit, a pre-processing software for Fluent, some researchers have conducted secondary development that improved modeling efficiency and made it easier for researchers to use Fluent software, enabling them to focus on model solving and analysis of computational results [29,30,31,32]. To reduce simulation workload for users, Xiao et al. [33] implemented batch processing for specific computational problems using parameters extracted from Fluent’s journal files. Li et al. solved the derived equations of motion electromagnetic field control through Fluent’s User-Defined Scalars (UDS) [34].

However, CFD calculation outputs often lack geospatial coordinate information, making it difficult to analyze spatiotemporal characteristics. Fortunately, Three-dimensional Geographic Information Systems (3DGIS) offer advantages in spatial data management and visual analysis, allowing for real-world objects to be abstracted into three-dimensional data for immersive display and powerful spatial analysis. Combining GIS and CFD models allows for exploitation of existing GIS data and the visualization and mapping capabilities of GIS. In recent years, there has been growing interest in the coupling of CFD and GIS. Murakami [35] concluded that GIS data structures are well suited for constructing complex CFD geometric models. Chu et al. [36] extracted the coordinates and heights of building polygons from GIS software and used them as input to the CFD pre-processing software for constructing a CFD geometric model. Coirier et al. [37] proposed a method to construct urban block models required for CFD simulations based on GIS data and conducted atmospheric pollution dispersion simulations using the model. Wong et al. [38] investigated the efficiency of using GIS data in CFD models and the sensitivity of the CFD results to different GIS data formats and discovered that raster format was more efficient, and contour data demanded considerable effort. Hooff [39] constructed a high-quality CFD geometric model of an urban residential area by importing GIS data into the CFD pre-processing software and simulated wind and thermal environments based on the geometric model. Maohui Zheng et al. [40] combined GIS and CFD software to simulate the dispersion of toxic gases in urban blocks and concluded that coupling GIS software with CFD software provides better visualization of the pollutant dispersion in complex buildings. Yunkai Xu [41] proposed a method to build complex urban building complexes for CFD simulation based on ArcGIS software and completed the loose coupling between Fluent software and ArcGIS software through intermediate data files. Xiaoyu Jiao [42] proposed a research scheme for coupling GIS and CFD for traffic pollutant dispersion simulation. Shouzhi Chang et al. [43] integrated CFD and GIS to identify and construct urban ventilation corridors (UVCs) in Changchun City, China. Wenzhong Wang et al. [44] investigated the import of Digital Elevation Model (DEM) data from GIS into the CFD pre-processing software ICEM CFD for three-dimensional model construction and meshing. They then simulated the target wind farm area based on the mesh file and analyzed the obtained wind velocity flow field vector map.

Although the coupling of GIS and CFD has been widely applied in various fields, such as atmospheric pollution dispersion simulation, and some achievements have been made, the coupling is mostly achieved through loose coupling, which involves complex interoperation and data conversion between independent software. Therefore, this study aims to present a tight coupling model of GIS and CFD for atmospheric pollution dispersion simulation. To achieve this, we utilized the following approaches. Firstly, to overcome the time-consuming pre-processing involved in CFD-based simulations of atmospheric pollution dispersion in urban blocks, we designed a high-performance CPU/GPU coordinate parallel transformation algorithm and encapsulated geometry construction and meshing commands using a parametric design method. This approach eliminated the data barrier between GIS and ICEM CFD, enabling rapid and automatic geometry construction and meshing. Secondly, to address the issue of CFD software requiring complex parameter inputs to cater to the demands of different fields, we examined CFD models related to atmospheric dispersion and proposed a Fluent GUI (Graphical User Interface) command flow parameterization design method and a TUI (Text-based User Interface) command flow encapsulation method. This made it possible to solve CFD atmospheric pollution dispersion models. Thirdly, to address the issue of CFD model outputs lacking geographic spatial coordinate information, which makes it difficult to display and analyze the spatiotemporal characteristics of pollution dispersion, we propose the use of multi-angle spatial partitioning and spatiotemporal multidimensional visualization analysis. This enables intuitive expression and analysis of simulation results, facilitating the understanding of pollution dispersion in urban blocks. Finally, leveraging the proposed methods and the tight coupling concept of GIS and CFD, we developed a user-friendly prototype system for simulating, visualizing, and analyzing atmospheric pollution dispersion in urban blocks, which significantly reduces the research threshold of CFD in pollution source dispersion. The system makes it easier for researchers or decision-makers to study and analyze air pollution in urban blocks and is expected to be an efficient tool to reduce property and human casualties caused by sudden pollution source leaks. For illustrative purposes, the three-dimensional characteristics of pollutant dispersion will be analyzed for a mixed residential and commercial area in Fuzhou.

2. Methodology

2.1. Introduction of ICEM CFD and Fluent Software

ICEM CFD and Fluent are two CFD software tools developed by ANSYS, Inc. Our coupling strategy relied on these two software tools to generate high-quality meshes and simulate fluid flow and pollutant dispersion in urban blocks. ICEM CFD is a pre-processing tool and primarily designed for mesh generation and geometry preparation, and it offers a range of mesh types, including structured and unstructured meshes. Fluent, on the other hand, is a powerful CFD solver that can simulate complex fluid flow and heat transfer problems, including atmospheric pollution dispersion. It features a variety of turbulence models, such as the k-ε turbulence model, to accurately capture the behavior of turbulent flows. Both ICEM CFD and Fluent are widely used in academic research and industrial applications for simulating fluid flow and heat transfer problems.

2.2. Automatic Geometry Construction and Meshing

Construction and meshing of urban block geometry is the foundation for atmospheric pollution dispersion simulations. Traditional CFD pre-processing software utilizes manual geometric modeling or a large number of coordinate inputs for geometry construction. For urban blocks containing considerable building complexes, this process is complex, time-consuming, and cannot meet the demands of near-real-time simulations. In this paper, we propose a parametric design method to generate a command flow file for geometry construction and meshing, and we ultimately achieve automatic construction and meshing of urban block geometry. The technical roadmap is illustrated in Figure 1.

2.2.1. Extraction of Geometry Coordinates

The research objects of urban block geometry are the building complexes within the entire simulation area. As the shapes of individual buildings can be complex, they are simplified as regular polyhedrons composed of vertices, edges, and faces, with the number of vertices varying depending on the shape of the building (In GIS, some features that appear to be a straight line at first glance are actually composed of multiple vertices when zoomed in). In this study, we assume that the buildings are solid with flat roofs. The constituent elements of a building with a rectangular parallelepiped shape are shown in Figure 2.

The geometry coordinates of the building complex in the urban block are obtained in the following three steps: (1) Based on the spatial topological containment relationship, using the simulation area boundary vector data as query object and the building vector data as queried object, all building objects within the simulation area are obtained, denoted as $B_i=\left\{B_i\right|i=0,1,2\dots m\}$;. (2) Each obtained building object is composed of faces made up of vertices, with each vertex represented by its longitude and latitude coordinates $\left(\lambda ,\phi \right)$. All vertex coordinates that make up each building object are saved in clockwise order, and the height attribute Z of each building object is read. The coordinate of each vertex is denoted as $P_j=\{{\lambda }_j,{\phi }_j,Z|j=0,1,2\dots n\}$;. (3) Repeat step 2 for all building objects in the simulation area to obtain all geometry coordinate data $B_i{P}_j=\{{\lambda }_{ij},{\phi }_{ij},Z_i|i=0,1,2\dots m,j=1,2\dots n\}$.

2.2.2. Parallel Approach-Based Coordinate Conversion of Building Complex

The process of computing the corresponding plane coordinates $\left(X,Y\right)$ from the given geographic coordinates $\left(\lambda ,\phi \right)$ in a spherical coordinate system using the Gaussian projection is called the forward Gaussian projection. The formula for the forward Gaussian projection is as follows:

$$\mathsf{\theta }=\{\begin{array}{c}\arctan \left(\frac{b^2}{a^2}\times \frac{\tan \phi }{\cos \left(\lambda -3zoneID+1.5\right)}\right)\\ \arctan \left(\frac{b^2}{a^2}\times \frac{\tan \phi }{\cos \left(\lambda -6zoneID+3\right)}\right)\end{array}$$

(1)

$$X={{\displaystyle \int }}_0^{\mathsf{\theta }}\mathrm{a}/\sqrt{\cos {\mathsf{\theta }}^2+\left(a^2/b^2\times \sin {\mathsf{\theta }}^2\right)}\mathrm{d}\mathsf{\theta }$$

(2)

$$Y=\mathrm{r}\cos \mathsf{\theta }\tan \left(\lambda -6\times zoneID+3\right)+500000$$

(3)

where $a$ and $b$ are the semi-major and semi-minor axes of the Earth ellipsoid, respectively, $\phi$ is the latitude (radian value), $\lambda$ is the longitude (radian value), $zoneID$ is the zone number of the current location, and $X,Y$ represent the corresponding plane coordinates to be computed.

As the forward Gaussian projection involves a definite integral solution process, and the simulation area of urban block typically contains a large number of building vertices, using conventional CPU-based serial algorithms based on Gaussian projection for the conversion of all acquired $B_i{P}_j$ data into plane coordinates under the Cartesian coordinate system can be time-consuming. To improve the efficiency of coordinate conversion, two approaches of coordinate conversion, namely, CPU parallel and GPU parallel, are investigated in this section.

The CPU serial algorithm converts the coordinate data of building vertices one by one. Only after the conversion of one vertex is completed, can the conversion of the next vertex proceed. As a result, each vertex conversion occupies only a small amount of CPU resources, leaving a significant portion of CPU resources idle. To address this issue, we propose a CPU-based parallel coordinate conversion algorithm that processes the conversion of coordinates for multiple buildings in parallel, with each building being treated as a single unit. This approach enables the simultaneous conversion of coordinates for multiple buildings at the same time. Using buildings as parallel units also effectively avoids the problem of reduced parallel efficiency caused by frequent thread switching due to excessively small parallel units.

As can be seen in Figure 3a, the building object is used as a parallel unit, and a cached thread pool technique is applied to call all available threads from the thread pool to perform forward Gaussian projection on the read geographic coordinates of multiple buildings simultaneously. A lock mechanism is then used to ensure the thread-safe storage of the converted coordinates until the conversion of all the read building geographic coordinates is completed.

Although CPUs serve as the control core of an operating system, they are not well-suited for large-scale computationally intensive parallel operations. In contrast, GPUs are specifically designed for concurrent computation of large-scale data with highly uniform and non-dependent types. Therefore, based on CUDA, a NVIDIA’s general-purpose parallel computing architecture, a GPU parallel coordinate conversion algorithm is proposed.

The CUDA framework employs a grid-based and block-based organization of GPU threads, with each grid and block organized in a two- or three-dimensional matrix. Nevertheless, as the number of vertices in each building varies, organizing GPU threads into a two-dimensional matrix may result in unused threads during coordinate conversion of buildings with a low number of vertices, as demonstrated in Figure 4.

To address this issue, as can be seen in Figure 3b, the two-dimensional set $B_i{P}_j$ needs to be converted into a one-dimensional set for storage. This involves reformatting the data of each row, which contain all vertex data for a single building, to a format where each row only contains data for a single vertex. Then, GPU threads are organized using a one-dimensional grid and one-dimensional block. In addition, the CUDA also provides four built-in variables, $gridDim$, $blockDim$, $blockIdx$, and $threadIdx$, to determine the unique identifier for each thread. Specifically, each thread with a unique ID executes the conversion of vertex coordinate at the corresponding ID in the one-dimensional set. The detailed algorithm is demonstrated in Algorithm 1.

Algorithm 1. Algorithm steps of the parallel GPU coordinate conversion of building complex vertices in CUDA-based framework. GPU parallel calculation of building vertex geographic coordinates to plane coordinates
Input: Coordinates of all building vertices λφZ_in ← {BiPj}, Number of all building vertices n
Output: Space Cartesian coordinates of all building vertices after conversion XYZ_out ← {XYZi}
1	function __global__ Kernel(XYZ_out, λφZ_in, n)//GPU-side functions
2		id ← (blockIdx.x * blockDim.x) + threadIdx.x//Get thread id
3		if id < n do//Perform coordinate conversion for threads that meet the condition
4			//Perform forward Gaussian projection
5			λφZ ← λφZ_in[id]
6			X, Y, Z = GaussianForward(λφZ)
7			XYZ_out[id] ← X, Y, Z//Store the converted coordinates
8		end if
9	end function
10	function __host__ buildsProjection(XYZ_out, λφZ_in, n)//CPU-side functions
11		//Apply for space on the device end
12		cudaMalloc(dev_λφZ_in, …), cudaMalloc(dev_ XYZ_out, …)
13		//Copy data from host side to device side
14		cudaMemcpy(dev_λφZ_in, λφZ_in, …)
15		cudaMemcpy(dev_XYZ_out, XYZ_out, …)
16		block ← blockMax/2, grid = (n − 0.5)/block + 1//Design CUDA thread organization
17		Kernel<<<grid,block>>>(dev_XYZ_out, dev_λφZ_in, n)//Call core functions
18		//Copy the calculation results from the device side back to the host side
19		cudaMemcpy(XYZ_out, dev_XYZ_out, …)
20		//Release space requested on the device side
21		cudaFree(dev_λφZ_in), cudaFree(dev_XYZ_out)
22		end function

2.2.3. Parametric Design Method for Geometry Construction and Meshing

The ICEM CFD software does not provide any secondary development APIs, but it does offer a command flow file for batch processing. In order to automate geometry construction and meshing, we have encapsulated the relevant commands using a parametric design method to generate command flow files for the relevant operations.

As the urban block buildings studied in this paper are treated as regular cubes, the construction of the simplified building in ICEM CFD software requires three steps: (1) creating the top vertices and bottom vertices, (2) creating the top, bottom, and side edges by connecting adjacent vertices, and (3) creating the top, bottom, and side faces of the building by connecting adjacent edges. To illustrate the script command for creating a point, the script command flow is as follows:

$$\mathrm{ic}\_\mathrm{point}\left\{\right\}\mathrm{GEOM}\mathrm{pnt}.000,0,0$$

(4)

where the command “ic_point {}” creates a point in the GEOM group with the name pnt.00 at the location (x, y, z) of (0, 0, 0). All variables, except for the “ic_point {}” command, can be changed. Therefore, all variables in the command flow can be parameterized and encapsulated in a “CreatePoint” method. Similarly, we can encapsulate methods for “edge”, “face”, “group”, and other command flows. The detailed flowchart of the automatic construction method of the geometry is shown in Figure 5.

ICEM CFD offers two meshing methods for geometry: structured and unstructured meshing. Structured meshing requires a predefined grid topology with a fixed number of nodes and elements, and it has strict geometric rules, such as the edges of adjacent cells needing to be aligned. This means that a significant amount of manual adjustments and corrections are required during structured meshing, increasing the difficulty and complexity of automated operations and making it unsuitable for large-scale simulations. On the other hand, unstructured meshing does not require predefined grids and does not need to follow strict geometric rules, making it easier to generate meshes for a variety of different geometries. It has a fixed sequence of steps and is easier to parameterize. Therefore, despite the fact that the building vector data in this study is composed of regular polyhedrons, unstructured meshing is more suitable due to the complexity and scale of the urban block as a whole and the need for automation. The unstructured meshing process typically involves the following steps: (1) creating a geometry, (2) setting meshing parameters, (3) saving a TETIN file representing the geometry (*.tin), (4) running the tetrahedral mesh generator, (5) smoothing the mesh elements, (6) saving the final unstructured mesh (*.uns), and (7) generating the mesh (*.msh). The necessary commands are shown in

Table 1. Unstructured mesh related operations corresponding to the commands.

Operations	Commands
setting the global meshing parameters	ic_set_meshing_params global
setting geometric parameters for a specific geometry family	ic_geo_set_family_params
generating a tetrahedral mesh	ic_run_tetra
smoothing the mesh elements	ic_smooth_elements
saving the generated unstructured mesh	ic_save_unstruct
executing a user-defined script or external program	ic_exec

.

For instance, to generate a TETIN file that represents the geometry, the command $\mathrm{ic}\_\mathrm{save}\_\mathrm{tetin}\mathrm{file}\mathrm{only}\_\mathrm{visible}$ is used. Here, “file” is a placeholder for the name of the TETIN file that will be generated (*.tin). The parameter “only_visible” indicates whether only the visible parts of the geometry will be included in the TETIN file. The “only_visible” option indicates that, during meshing or output, only the visible parts of the geometry will be saved according to the current view, while the parts that are invisible or obstructed by other objects in the geometry will not be saved or included in the output file. The default value for this parameter is 0, which saves all parts of the geometry. Similar to the above-mentioned strategy, the command can be encapsulated in the “SaveModel” method, and other commands can also be encapsulated in their corresponding methods.

2.3. Fluent-Based Solution for Atmospheric Dispersion Models

The CFD software ANSYS Fluent has complex parameter input designed to meet the needs of various fields, and it does not provide secondary development APIs. In this section, we examine CFD models related to atmospheric dispersion and propose a Fluent GUI (Graphical User Interface) command flow parameterization design method and a TUI (Text-based User Interface) command flow encapsulation method to indirectly provide secondary development APIs for Fluent. Finally, we implement automated generation of model solver command flow files for the calculation of atmospheric pollution dispersion models in Fluent. The technical roadmap is illustrated in Figure 6.

2.3.1. GUI Command Flow Parameterization Design Method

By manually operating the Fluent software, a Journal File (*.jou) can be generated to record a sequence of actions performed in the GUI. For example, a command flow generated by specifying a boundary condition for a particular zone is:

$$(\mathrm{cx}-\mathrm{gui}-\mathrm{do}\mathrm{cx}-\mathrm{set}-\mathrm{list}-\mathrm{selections}“\mathrm{Boundary}\mathrm{ConditionsTable}1\mathrm{List}2(\mathrm{Zone})”‘(1\left)\right)$$

(5)

This command selects the second item (index 1) in the list named “Boundary ConditionsTable1List2 (Zone)” in the Fluent GUI. While this index number is variable and not available to external programs, the name of the corresponding bounding area (e.g., inlet1) can be determined. Thus, the index can be replaced with the name, which is the key to parametric design of the GUI and allows for the parameterized encapsulation of Fluent-related operational command flows.

2.3.2. TUI Command Flow Encapsulation Method

The Fluent GUI command flow is not well organized and requires a call window display, making its execution relatively inefficient. In contrast, Fluent TUI commands can be directly entered and executed in Fluent’s Console or saved as a log file (*.jou) in a text editor for later execution by Fluent. A TUI command is composed of a directory and a command, where the directory is selected first, followed by executing the corresponding command under that directory. Fluent TUI provides 14 categories of first-level directories, each corresponding to a functional description listed in

Table 2. Fluent TUI first-level directory commands.

Directory Name	Functional Descriptions
adapt/	contains commands related to mesh adaptation
adjoint/	contains commands related to adjoint solver
define/	contains commands to define materials, boundary conditions, and other simulation parameters
display/	contains commands to control the display of the Fluent GUI
file/	contains commands to import or export files
mesh/	contains commands related to mesh generation and manipulation
parallel/	contains commands related to parallel
plot/	contains commands to create plots and animations
report/	contains commands to generate reports, such as force or mass reports
server/	contains commands to control the Fluent server
solve/	contains commands related to the solution of the fluid problem
surface/	contains commands related to surface modeling
turbo/	contains commands related to turbomachinery simulation
views/	contains commands to create, save, and restore views of the Fluent GUI

.

For example, the TUI command used to define a boundary condition is:

$$/\mathrm{define}/\mathrm{boundary}-\mathrm{conditions}/\mathrm{velocity}-\mathrm{inlet}\mathrm{inlet}1\mathrm{n}\mathrm{n}\mathrm{y}\mathrm{y}\mathrm{n}10\mathrm{n}101325$$

(6)

where “/define/boundary-conditions/” refers to the directory in which the command is located, “velocity-inlet” specifies the type of boundary condition, and “inlet1” is the name of the boundary. The input options for the command vary depending on the model being used. For instance, if the viscous model is set to Laminar and all other models are turned off, the necessary options would be the boundary name, velocity method (which has three options: “Magnitude and Direction”, “Components”, and “Magnitude, Normal to boundary”), reference frame (which has two options: “Absolute” and “Relative to Adjacent Cell Zone”), velocity magnitude, and initial gauge pressure. In the command, “n” and “y” represent “no” and “yes,” respectively. For example, “n n y” indicates the third one of velocity methods, and “y” indicates the first one of reference frames. The velocity magnitude is set to 10 m/s, and “n” preceding the velocity magnitude indicates no use of a velocity profile. Similarly, “n” preceding the initial gauge pressure indicates no use of a profile for the initial gauge pressure, which is set to 101,325 Pa. Note that, when selecting the first reference frame option, the command is simply “y”, not “y n”, and when selecting the second reference frame option, the command is “n y.” The other commands follow the same pattern.

As we can see, there are two major challenges in using TUI commands: (1) the parameters are not fixed and may vary depending on the specific models being used, and(2) the parameters are often represented in the form of “n” and “y”, which can be difficult to interpret. In this study, an object-oriented approach is used to encapsulate the TUI command flow of each function into easy-to-use functional objects.

As an example of the encapsulation process for the TUI command of setting boundary conditions mentioned above, the process is illustrated in Figure 7. Simulating urban air pollution dispersion involves four models: energy equation, turbulence model, radiation model, and species model. Based on our research on Fluent TUI commands, ten interfaces need to be designed. The boundary conditions involve five types: velocity-inlet, pressure-inlet, mass-flow-inlet, pressure-outlet, and mass-flow-outlet, which correspond to five methods under each of the ten interfaces. After designing the interfaces, an implementation class is provided to implement all the interfaces. However, due to the complexity of Fluent software operations, each method in the implementation class involves numerous parameters, so the parameters of each method are encapsulated into a parameter class. Finally, a boundary condition setting class (BoundaryConditionSet) is designed, which is accessible for developers.

2.4. Visualization and Analysis of Simulation Outputs

After the simulation is finished, Fluent stores the relevant calculation values of the entire study area in a binary text file, including node numbers, XYZ coordinates, pressure values, XYZ velocity components, velocity magnitude, velocity angle, and mass fractions of different materials at the corresponding coordinate positions. As a result of these simulation outputs being in text format, they are typically used for visualization and analysis outside of GIS, and they have not been effectively integrated within GIS for impact assessment. However, by coupling GIS and CFD, it is possible to perform spatiotemporal multi-dimensional visualization and spatial analysis on simulation outputs.

2.4.1. Conversion of Simulation Outputs to 3DGIS Data

The simulation outputs are typically saved as text data and cannot be directly imported into 3DGIS. To integrate the simulation outputs into 3DGIS, Inverse Gaussian projection is needed to convert the outputs into 3DGIS-compatible data.

Due to the large number of nodes contained in each time step of the simulation outputs, a considerable amount of conversion time is required. Therefore, the CPU parallel and GPU parallel algorithms mentioned in Section 2.2.2 are required to reduce the conversion time. The resulting converted CSV data is then transformed into a three-dimensional point Shapefile format using the geoprocessing tool interface (IGeoProcessor) in the ArcObjects library. The specific workflow is shown in Figure 8.

2.4.2. Multi-Angle Spatial Partitioning of Three-dimensional Spatial Simulation Outputs

After converting each node of the simulation output to a three-dimensional vector point, a large number of discrete three-dimensional points are spread across the entire study area, as shown in Figure 9. Visually, only the three-dimensional vector points distributed on the surface of the study area can be observed (blue-green points in Figure 9), while the internal three-dimensional points are difficult to observe, as they are obscured by the external points. To address this issue, our developed system implements the horizontal and arbitrary direction partitioning of three-dimensional spatial simulation outputs, which is especially useful for managing large-scale three-dimensional spatial datasets generated from atmospheric dispersion simulations. This approach enables more efficient processing of complex three-dimensional spatial data and facilitates more accurate and effective spatial analysis. The flowchart is shown in Figure 10.

The attribute query function in GIS can be used to obtain point objects located on the horizontal cutting plane. However, since three-dimensional points have volume, simply querying objects that satisfy the height attribute at the level of the horizontal cutting plane is not enough. Instead, we need to query objects that fall within a certain height range around the horizontal cutting plane, with the default range being the diameter size of the three-dimensional point. The default diameter size of three-dimensional points in GIS is typically 2 or 4 pixels, which is automatically adjusted based on the geographic extent and scale of the data. Red three-dimensional vector points represent the horizontal section data of simulation data in Figure 9 at a height of 1.5 m, as shown in Figure 11a.

However, using attribute query to perform arbitrary directional partitioning can be very complex or even impossible. Fortunately, GIS offers spatial topology query, which allows us to obtain point objects in cutting planes oriented at any angle with respect to the dataset. The specific steps are as follows: (1) create a baseline for the cutting plane; (2) set a symmetric buffer on both sides of the cutting plane baseline to form the cutting plane base, with the buffer distance defaulting to the diameter size of the discrete points; (3) stretch the cutting plane base to form the partitioning volume; (4) conduct a spatial containment query of the simulated vector points with the partitioning volume; (5) save the query results as a section layer. As shown in Figure 11b, the brownish-red three-dimensional vector points represent the cross-sectional data of the simulation data of Figure 9 at the yellow cutting plane.

2.4.3. Spatiotemporal Multidimensional Visualization and Analysis

By tightly coupling CFD with 3DGIS, it becomes possible to directly access and visualize the converted simulation results data in 3DGIS, enabling spatial analysis of the simulated data. However, since atmospheric pollution dispersion is a process-based phenomenon, the temporal expression of the simulation results is equally important as their spatial position. Therefore, this paper investigates two mixed forms of layer visualization and Image animation visualization to add a temporal dimension to the three-dimensional spatial simulation result data and achieve a spatiotemporal multidimensional visualization of the simulation result.

The steps for layer visualization are as follows: (1) partition the three-dimensional simulation outputs of all moments based on a spatial partitioning method to obtain the section data of interest; (2) render the section data at each moment using symbol methods, such as the gradient color method; (3) load and display the section data of interest at each moment in a time sequence. However, the layer-based loading approach does not provide a dynamic view of the spatial distribution of pollutants over time. Thus, image animation visualization is proposed based on layer visualization to solve this problem. The detailed steps are as follows: (1) set the visibility of the rendered section layers at each moment to ensure that only the layers at the same moment are visible and save the layer data at this moment in the current view as images named with the time of the current moment; (2) play the saved images in a time sequence.

2.5. Implementation of the CFD Coupled with 3DGIS

In this study, a prototype system is developed for simulating, visualizing, and analyzing atmospheric pollution dispersion based on the tight coupling model of GIS and CFD. The user interface is shown in Figure 12. The system is built using the .NET open-source development platform and leverages the ArcObjects component object set to create three-dimensional virtual geographic scenes.

Figure 13 presents the entire architecture of the prototype system, which is designed with a three-layer structure of C/S architecture (C/S architecture stands for Client/Server architecture. It refers to a computing model where the client-side and server-side of a system are separated into distinct layers, with each layer responsible for different functions. The client-side is typically responsible for user interface and presentation, while the server-side is responsible for data storage and processing). The data storage layer mainly uses local file system databases and commonly used geospatial databases in GIS, including building vector data with height attributes, ANSYS ICEM CFD command flow files (*.rpl) for automating geometry construction and meshing, ANSYS Fluent command flow files (*.jou) for atmospheric pollution dispersion simulation calculations, simulation result data from calculation, related data for spatial analysis of simulation results, and other data obtained from spatial analysis. The business logic layer mainly realizes the construction of three-dimensional virtual geographic scenes, the automatic construction of urban block geometry and mesh, the simulation calculation of atmospheric pollution dispersion in urban blocks, and the spatiotemporal multidimensional visualization and analysis of the simulation result data. The representation layer is a desktop client application built on the .NET platform.

The system is comprised of four major modules, namely: the creation of three-dimensional virtual geographic scenes, the automatic geometry construction and meshing for urban blocks (indicated by the red box in Figure 12), simulation of pollution dispersion (represented by the blue box in Figure 12), and spatiotemporal multi-dimensional visualization and analysis of simulation results. A more detailed breakdown of system functions is shown in the business logic layer of Figure 13.

Based on the approach investigated in Section 2.3, we designed two parameter-setting modes for atmospheric pollution dispersion simulation, namely, the brief mode and the expert mode. The former has pre-configured settings for specialized parameters, making it suitable for users who do not have specific requirements for the assumptions of fluid mechanics. For example, in the brief mode, CFD models only enable the energy equation, the standard k-epsilon model, and species transport model. After setting up the pollutant source and boundary of the study area (Figure 14), users can perform atmospheric pollution dispersion simulation by simply using the brief mode to specify parameters, such as pollutant source material, velocity-inlet boundary with its velocity magnitude, mass flux of pollutant source material, and the duration of the simulation. This simplifies the model parameter setting process and significantly reduces the learning curve for users. The latter provides more flexibility in parameter settings, allowing the user to select and configure parameters according to specific situations, which is ideal for users with adequate knowledge of fluid mechanics.

By tightly coupling GIS and CFD, users or decision makers can perform a variety of visualization and analysis operations on simulated results, including: (1) gradual color rendering, (2) conversion of simulated vector data to raster data using raster interpolation algorithms such as the inverse distance weight method (IDW) and Kriging method, (3) attribute query and spatial relationship query, (4) overlay analysis of the simulated data and residential block polygon data, (5) visualization of wind field vector maps at a specific moment, (6) multi-angle spatial partitioning, and (7) spatiotemporal multidimensional visualization and analysis. For instance, users can use the spatial partitioning method proposed in Section 2.3.2 to partition the simulation data for all moments, and then they can create animations for specified sections using the Animation Demo Window (Figure 15).

3. Experiments and Results

3.1. Study Area and Data

A rectangular area bordered by “Sangao Road”, “Gangtou Road”, “Daping Road”, and “Lianjiang Road” was selected as the target region for atmospheric pollution dispersion simulation in Fuzhou, China. The area consists of typical components found in an urban block, including seven residential areas, two schools, one shopping mall, and one integrated service department (which is a large building that houses various services). The satellite image of the area and road network of the area can be seen in Figure 16a. The experimental data comprises vector data of buildings in Fuzhou city and each residential area in the simulation area. The building vector data attribute contains the number of floors of each building, and the residential areas vector data attribute includes the name of each residential area and its population, as shown in Figure 16b.

3.2. Comparison of The Efficiency of Coordinate Conversion Algorithm for Building Complex

This study focuses on a simulation area containing 116 buildings and 4237 bottom vertices, where the conversion of building vertices from geographic to plane coordinates was implemented using CPU serial, CPU parallel, and GPU parallel Gaussian projection algorithms. The CPU and GPU hardware-related parameters are shown in

Table 3. Hardware-related parameter information of CPU and GPU.

Parameter	CPU	GPU
Model	Intel(R) Xeon(R) Gold 5228	NVIDIA GeForce RTX 2080 Ti
Other properties	Number of cores: 24; Number of logical processors: 48; L1 cache: 1.5 MB; L2 cache: 24.0 MB; L3 cache: 33.0 MB	Maximum number of blocks in each dimension of the grid: 2,147,483,647, 65,535, 65,535; Maximum number of threads in each dimension of a block: 1024, 1024, 64

.

As shown in Figure 17, the CPU serial algorithm required 4162.576 s (69.4 min) to perform the forward projection, whereas the CPU parallel algorithm took only 314.871 s (5.2 min), indicating an efficiency improvement of almost 13.3 times. On the other hand, the GPU parallel algorithm took 170.277 s (2.3 min) to perform the forward projection, resulting in an efficiency improvement of almost 25 times compared to the CPU serial algorithm and twice as much compared to the CPU parallel algorithm. Overall, these results demonstrate a significant enhancement in efficiency.

3.3. Validation of Geometry Construction and Meshing

Based on the geometry construction and meshing approach proposed in this study, the generated command flow file can be executed to automatically construct the geometry of buildings within the simulation area and perform meshing on the generated geometry. The constructed buildings are shown in Figure 18, where the left figure is the building vector data from the experimental data, and the right figure is the automatically generated building geometry suitable for CFD. The comparison shows the capability of the geometry construction function implemented in our system.

The results of meshing and the bar chart of mesh quality metrics are shown in Figure 19 and Figure 20, respectively. In this study, the mesh quality was evaluated using the Determinant metric, which is a measure of the distortion of each mesh element. The horizontal axis represents the quality of the mesh, which is calculated based on the determinant value of the Jacobian matrix, and the vertical axis represents the number of mesh elements with the corresponding mesh quality. Higher determinant values indicate better mesh quality. The bars with upward arrows indicate a higher number of elements. By clicking on a specific bar in the ICEM CFD software, the exact number of elements can be obtained. The bar chart shows that all mesh qualities are above 0.45, which meets the mesh quality requirement (0.3) of CFD simulation. Therefore, the proposed automatic meshing approach in this paper is reliable and effective.

3.4. Case Study of Chlorine Dispersion Simulation

The correctness and usability of the automated geometry construction and meshing method proposed in this paper have been verified in the previous section. Therefore, in this section, we utilize the automatically constructed geometry of buildings in the simulation area to simulate the dispersion of chlorine. Assuming that, in the simulated area of Fuzhou City, at the intersection of “Sangao Road” and “Lianjiang South Road” near the Metro Cash & Carry (Fuzhou Cangshan Store), an accident occurred with a truck carrying a large amount of chlorine, causing the leakage and dispersion. The wind speed at the time of the accident was about 8 m/s, with a southeast wind direction, and the impact of chemical reactions is not considered. The simulation area and boundary surfaces are shown in Figure 21. The dimensions of the study area are as follows: the length of the inlet-S boundary is 555 m, the length of the inlet-E boundary is 596 m, the length of the outlet-W boundary is 684 m, and the length of the outlet-N boundary is 546 m. The height of the study area is set to the height of the tallest building within the area, which is 20 m in this case.

3.4.1. Simulation Parameter Setting and Model Solving

To validate the effectiveness of the prototype system implementing the aforementioned methods, we focused on the system’s performance, rather than the complexity of the CFD simulation. Therefore, we opted for simple parameter settings using the brief mode of the system. Specifically, the model is set up with the energy equation, standard k-ε turbulence model, and species transport model. Gravity is considered in the -Z direction with a magnitude of −9.81 m/s². The solver type is pressure-based with the solver velocity attribute set to absolute velocity, and the time option is set to transient. Chlorine was added to the materials as a fluid, and the mixture consists of two gases, air and chlorine. The inlet-S and inlet-E surfaces were set to the velocity-inlet boundary condition, while the pollution truck was set to mass-flow-inlet boundary condition. To simplify the flow field, we chose to specify constant values for important parameters, such as mass flux and velocity magnitude, instead of using profiles as boundary conditions. Specifically, the velocity magnitude of the inlet-S was set to 8 m/s with air as the only species, while the velocity magnitude of the inlet-E was set to 0.6 m/s with air as the only species. The mass flux of the pollution truck was set to 6 kg/(m²-s) with chlorine gas as the only species. The turbulent intensity was set to 5% and the turbulent viscosity ratio was set to 10 for both the velocity-inlet and mass-flow-inlet. The outlet-W and outlet-N surface are set as outflow. The top and bottom surfaces of the study area and the building are set up as no-slip walls. The detailed parameters for the boundary conditions are shown in

Table 4. Boundary condition parameters for study area.

Zona Name	Type	Boundary Conditions
Pollution truck	mass-flow-inlet	Mass Flow Method: Mass Flux; Reference Frame: Absolute; Mass Flux: 6 kg/(m2-s); Initial Gauge Pressure: 0; Direction method: Normal to Boundary; Temperature: 300 k
inlet-S, inlet-E	velocity-inlet	Velocity Method: Magnitude, Normal to boundary; Reference Frame: Absolute; Initial gauge Pressure: 0; Temperature: 300 k
outlet-W, outlet-N	outflow	Flow Rate Weighting: 1
top and bottom surface, building	wall	Wall Motion: Stationary; Shear Condition: No Slip; Roughness Models: Standard; Temperature: 300 k

. The SIMPLEC pressure–velocity coupling scheme is selected as the solution method. The simulation results are saved as binary text files for the entire simulation area at each time step. The simulation runs for 120 s with a fixed time step of 1 s, and a maximum of 20 iterations are performed at each step.

Based on the approach proposed in Section 2.2, our system can automatically generate the command flow file for parameter setting described above and other related computing operations. The ANSYS Fluent software can be invoked in the background by our developed system to execute the command flow file and complete the simulation calculation. The residual plot after the simulation calculation is displayed in Figure 22, which shows relatively small iterative calculation errors for the continuity equation, energy equation, turbulence equation, velocity equation, and species transport model. Hence, the correctness of the simulation calculation in this case can be verified.

3.4.2. Analysis of Simulation Results

In this case, the state of chlorine dispersion was simulated for each second of the 120-s time period. Figure 23 shows the chlorine gas dispersion states at every 10-s interval, starting from the first second up to 100th second, at a height of 1.5 m on the horizontal section. The figure reveals that, as time progresses, the overall range of chlorine dispersion becomes increasingly larger. The high-concentration chlorine range (represented by the red area) rapidly expands in the first 50 s, but it reaches an equilibrium state after that under the prevailing conditions in the simulated area. This is due to the balance between the release and consumption of chlorine within this range, which is caused by the wind dispersing the chlorine and reducing its concentration. In contrast, the low and medium concentration chlorine ranges continue to expand, and their trend towards an equilibrium state is not evident. Therefore, timely actions should be taken during the early stages of chlorine leakage to prevent continuous leakage. Shortening the dispersion time can achieve a reduced range of low and medium concentrations of chlorine dispersion. Additionally, since the range of high-concentration chlorine will reach an equilibrium state quickly, its area should be quickly locked down, and rescue efforts should be focused on those within this range.

The spatiotemporal multidimensional visualization presented in Figure 23 was obtained through the following steps: (1) horizontal partitioning was performed on the simulation data at a height of 1.5 m using the horizontal partitioning method; (2) the chlorinous mass fraction field of the section data was converted into raster data using the raster interpolation method; (3) image animation visualization was used for visualization. In addition to applying the aforementioned methods to the simulation data, users can also perform other spatial analysis methods on the simulation results in the system, such as reclassifying the pollutant mass fraction based on pollution standards to classify the study area into regions of different pollution levels.

4. Discussion and Conclusions

Our team has been devoted to researching the application of GIS in atmospheric pollution accidents. For example, Huang et al. [45] integrated the atmospheric dispersion model CALPUFF with GIS, effectively integrating the numerical simulation results of the CALPUFF model into a three-dimensional terrain scene. They also realized the dynamic diffusion of model simulation results in the time dimension, as well as hierarchical classification in the spatial dimension. The disaster simulation information is data calculated by the model at a certain moment and a certain area, such as the concentration of pollutants and the probability of casualties. Tang et al. [46] organized the disaster information in GeoJSON format and designed a visualization method for sudden atmospheric pollution simulation disasters based on GeoJSON, achieving the visualization of atmospheric pollution disasters on mobile devices. However, we opted for a CFD model that is more suitable for simulating atmospheric pollution dispersion within urban blocks [11,12], rather than coupling GIS with a Gaussian model, and we successfully integrated CFD model outputs into GIS for spatiotemporal multidimensional visualization analysis, which effectively solved the issue of the lack of geospatial information in CFD model outputs. Furthermore, we introduced, for the first time, the tight coupling of GIS with CFD software ICEM CFD and Fluent, not only integrating the CFD model outputs into GIS, but also allowing users to rapidly construct the building complex geometry and meshing required for simulating pollution dispersion within urban blocks, which resolved the complexity and time-consuming nature of constructing the necessary geometry for CFD simulation of urban blocks. We also encapsulated the Fluent-based command flow, indirectly providing Fluent with secondary development APIs, and we achieved the automatic generation of model calculation command flow files to solve atmospheric pollution dispersion model, addressing the issue of Fluent’s high entry barrier due to its complex parameters tailored to various fields. The specific research results are as follows:

(1)

We propose an approach for automating the construction of geometry and meshing required for CFD simulations of urban blocks. Specifically, we design both CPU and CUDA-based GPU parallel algorithms to convert building vertex coordinates in spherical coordinate systems to Cartesian coordinates in the plane quickly. We also propose using a parametric design method to encapsulate geometry construction and meshing commands to achieve automatic and rapid construction of geometry and unstructured meshing. Through experiments on the study area, we validate that the constructed geometry is correct, and the mesh quality meets the requirements with all values above 0.45. Additionally, the CPU and GPU parallel algorithms are 13.3× and 25× faster than serial, respectively.

(2)

We investigated CFD models related to atmospheric dispersion and proposed a parameterization design method for Fluent GUI command flow and an encapsulation method for TUI command flow, providing secondary development APIs for Fluent customization. By designing simple parameter interaction interfaces, users enabled solving CFD atmospheric dispersion models in urban blocks.

(3)

We propose a spatio-temporal multidimensional visualization and analysis method based on three-dimensional GIS for simulation results. Specifically, we developed a method for converting CFD simulation results into three-dimensional GIS data and achieve the coupling of CFD simulation results with three-dimensional GIS. By using three-dimensional GIS attribute and spatial topological relationship queries, we enable multi-angle spatial partitioning of simulation results. We also propose two methods for achieving spatiotemporal multidimensional visualization and animation of simulation results: layer visualization and image animation visualization.

(4)

We integrated the above-mentioned methods to develop a system coupled GIS and CFD for atmospheric pollution dispersion simulation in urban blocks. This system provides a user-friendly and easy-to-use tool for relevant departments and researchers to simulate the dispersion of atmospheric pollutants in urban areas, as well as to easily explore deep-level information.

Due to the limitation of research conditions, the research in this paper has the following shortcomings:

(1)

The geometry is constructed only considering the buildings within the urban blocks, but not the topography, public facilities and green vegetation of the urban blocks.

(2)

The coupling of 3DGIS and CFD in this study is achieved through a tight integration of ArcGlobe and ANSYS software based on the .NET Framework technology framework. This involves converting GIS data into the required format for geometric construction in ICEM CFD and integrating Fluent simulation results back into GIS. All of thiese data are stored locally on the user’s computer. Consequently, our approach requires local access and conversion of shared data, which reduces efficiency and necessitates the installation of third-party software, thereby making the environment configuration complex.. In the future, open-source GIS libraries and CFD libraries can be chosen to fully integrate 3DGIS and CFD based on the method proposed in this paper.

(3)

For large simulation areas, the huge amount of vector point data generated by the simulation results can significantly decrease the efficiency of three-dimensional visualization, resulting in a poor user experience. In the future, it may be possible to explore the use of popular front-end two-dimensional/three-dimensional map engines, such as Cesium and Three.js, to enable three-dimensional visualization and related analysis of simulation results on the web, which could improve efficiency and user experience.