Application developed for PhD thesis entitled: "Software development for the optimization of the influence of wind flows within energy applications and sustainable town planning"

The thesis is available by clicking on the corresponding badge.

The open-source tool is developed and integrated into the Kratos Multiphysics software package (the GitHub page). Kratos is designed as an Open-Source framework for the implementation of numerical methods to solve engineering problems. It is written in C++ and designed to allow collaborative development by researchers focusing on modularity and performance.

Six separate modules are realized to make the tool as robust as possible, each capable of performing specific operations. In detail, the modules are:

• geo_data: it allows obtaining the terrain and building geometry data providing only the coordinates of the center and the radius of the domain to be analyzed.
• geo_preprocessor: it performs preliminary operations on the data used for realizing the numerical model.
• geo_importer: it imports information from different formats such as STL, OBJ, XYZ etc. However, the list can be extended to make the process more versatile.
• geo_mesher: it allows computing the volumetric mesh. For this modulus, two different approaches can be followed. The first method is based on calculating the distance to the terrain surface for performing a metric-based re-mesh procedure. This step can be repeated iteratively to increase the mesh quality. The open-source software MMG via API in Kratos is used for this purpose. On the other hand, the second approach sets the minimum and maximum dimensions of the contour surface elements.
• geo_building: it performs operations on the building geometry, such as, for example, positioning within the domain and removing buildings outside a given boundary or below a selected elevation.
• geo_model: it allows both to split the lateral surface of the domain to study different incoming wind directions and create boundary conditions following the conventions in Kratos.

The information for the realization of the numerical model is obtained through two different datasets: one for the terrain and one for the buildings. Both databases are managed by the geo_data module.

GDEM (Global Digital Elevation Model) images obtained through the ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) file are used. The ASTER instrumentation is onboard NASA's Terra spacecraft.

A digital elevation model (DEM) is a digital raster representation of ground surface topography or terrain. Each raster cell (or pixel) has a value corresponding to its altitude above sea level.

ASTER-GDEM (Version 3) files are divided into 1-degree by 1-degree data area tiles. The names of the individual data tiles refer to the latitude and longitude at the geometric center of the lower left (southwest) corner pixel. GDEM files can be retrieved through the Earthdata website or API.

Figure 1 shows the tile named ASTGTMV003_N42E014, related to the area under evaluation (Pescara).

Figure 1 - ASTER-GDEM tile of interest area

The geometric model concerning the buildings is obtained through OpenStreetMap (OSM) website or the Overpass API. The OSM provides maps of the world with free content. Geographic data in the OSM is under a free license (the Open Database License). This license makes possible the following operations:

• Copy, distribute, use the database,
• Create works from the database,
• Modify, transform, and develop the database.

Each element within the OSM can be represented through four different types:

• node: it indicates both point objects such as trees or traffic lights and, in a simplified way, commercial activities, places of interest etc.,
• way: it is a set of unclosed points and can describe a road, river, path, etc.,
• polygon: it is a closed line and can define a building, lake, park, etc.,
• relation: it is a set of the previous elements. It can be defined as a container used to describe complex objects. For buildings with complex geometry, the relationships are of the Multipolygon type. An example is shown in Figure 2 where the insides polygon represents the interior patio of the building.
  

Figure 2 - Multipolygon in OSM. way_1: outer perimeter,
way_2: inner perimeter

Considering the open-source nature of OSM, some areas may be not sufficiently mapped, and, therefore, some buildings are not present. In this case, the buildings can be inserted manually into the OSM through the dedicated section.

Figure 3 shows the area analyzed in this work and the zones in which some buildings are included are highlighted in red.

Figure 3 - Area analyzed in this work: pre-existing buildings (in yellow)
and the areas where buildings are added (in red)

The numerical model, generated for performing CFD analyses, is created using information referred to both terrain and building datasets. The operations necessary for the realization of the final numerical model are described in the following subsections.

Data obtained through the Earthdata portal can be converted into 3D files using the DEMto3D plugin within the QGIS (Quantum Geographic Information System) software. This tool is designed for the conversion of DEM files into 3D printing files.

Tests conducted on small models (with a resolution of a few hundred meters) produce satisfactory results. The tool, however, is strongly limited because it does not allow the conversion of large areas of a few square kilometers useful for realizing numerical models of urban quarters.

For this reason, an automated procedure is created (within the geo_data module) for downloading and converting the DEM file into the OBJ format.

After setting the bounding box of the area of interest, the ASTER-GDEM tiles in which the domain falls are downloaded (Figure 4 left) via Earthdata API and the tiles are cut according to the coordinates of the bounding box (Figure 4 right). Subsequently, the DEM file is converted into an OBJ file, thus obtaining a 3D terrain model (Figure 5).

Figure 4 - ASTER-GDEM tile in which the domain falls (left),
and the bounding box of the domain (right)

Figure 5 - 3D model of the terrain

A specific condition can appear when the area of interest is situated on several tiles. The worst situation is when it is at the intersection of 4 different tiles (Figure 6). In this case, therefore, it is necessary to merge the tiles and, afterwards, crop the domain of interest.

Figure 6 - Domain positioned at the intersection of 4 tiles

The obtained OBJ file is imported using the geo_importer module.

With the geo_preprocessor, however, a circular portion of the terrain surface is cropped and divided into three different areas:

• A central circle with the real orography, where the buildings are placed,
• An intermediate annular portion with only the real orography (without the buildings),
• An outer annular portion, subject to a smoothing procedure, is modeled with a simplified orography with extreme nodes located at zero elevation.

This configuration avoids a discontinuity in the incoming wind speed profile.

The smoothing procedure consists of changing the z-coordinate of the nodes according to the following equation:

$$ z_i = {(z_i - z_{min}) \beta + z_{min}} $$

where $z_i$ is the coordinate of the i-th node, $z_{min}$ is the minimum elevation of the domain, and $\beta$ is a reductive coefficient that can take values between 0.0 and 1.0. It is calculated as:

$$ \beta = 1 - \frac{d_i - r_{ground}}{r_{bondary} - r_{ground}} $$

where $d_i$ is the distance of the i-th node from the center of the domain, $r_{ground}$ is the inner radius of the annular portion with simplified orography, and $r_{boundary}$ is the radius of the domain.

The result of the operations just described is a cloud of points that is converted into a 2D mesh using the Triangle library (available in Python).

In Figure 7 it is shown the scheme of the points mentioned above.

Figure 7 - Portions of the terrain surface

The geo_data module allows, in addition, to download (via the Overpass API) and convert the data provided by the OSM into three-dimensional models.

Data via the Overpass API are downloaded by providing the bounding box coordinates of the area of interest and the type of data to be extracted. In the specific case, only the information about the buildings is downloaded and stored in a GeoJSON file.

GeoJSON is an open format that is used to store spatial geometries whose attributes are described through the JavaScript Object Notation. Possible geometries are points, lines, polygons, and multiple collections of these types. The typical structure of a GeoJSON file is like the following:

{
  "type": "way",
  "id": 333362603,
  "bounds": {
    "minlat": 42.4424610,
    "minlon": 14.2067793,
    "maxlat": 42.4429295,
    "maxlon": 14.2073653
  },
  "nodes": [
    3404740665,
    3404740666,
    3404740667,
    3404740668,
    3404740665
  ],
  "geometry": [
    {"lat": 42.4429295, "lon": 14.2072623},
    {"lat": 42.4428541, "lon": 14.2073653},
    {"lat": 42.4424610, "lon": 14.2068823},
    {"lat": 42.4425326, "lon": 14.2067793},
    {"lat": 42.4429295, "lon": 14.2072623}
  ],
  "tags": {
    "building": "apartments",
    "building:levels": "6"
  }
}

Information regarding the height of buildings can be provided either as number of floors or as total height. In the first case, the full height is calculated as the product of the number of floors and the storey height, which, by default, is set at 3.0 m.

Within the geo_data module, the GeoJSON file is converted into an OBJ format. To do this, first, each building is processed iteratively, and the ground footprint geometry is created; then, it is extruded upwards, obtaining the 3D model. However, this procedure has some limitations:

• tapered buildings are not correctly represented.
• In the case of Multipolygon, only the external perimeter is considered (buildings with internal courtyards are not well modelled).
• Adjacent buildings are combined into a single building, and the joint surface is removed.
• Overlapping buildings are combined into a single building.
• Since the community uploads the information, it may happen that some buildings do not have height information. In this case, a default height is set to create the 3D model of the building.

The result of the procedure just described is shown in Figure 8 and Figure 9.

Figure 8 - 3D model of the buildings obtained through the developed procedure
(Top view - city of Pescara)

Figure 9 - 3D model of the buildings obtained through the developed procedure
(Axonometric view - city of Pescara)

As mentioned in the previous section, the buildings cover only the central portion of the final numerical model. For this reason, a function (in geo_building module) is implemented to remove structures outside the given perimeter.

Before merging the geometry of the buildings with one of the terrain, an additional step is necessary to place the buildings at the correct elevation. The OSM, in fact, does not provide the elevation information and, therefore, the buildings are all located on a zero level.

Thanks to the potential provided by Kratos, a function that can calculate the distance (on the z-axis) that each building has with the geometry of the ground is developed within the geo_building module. This function allows to shift the building at the exact altitude. The outline of what has just been described is shown below:

Figure 10 - Shifting of buildings on the terrain surface

The operations described in the previous subsection allow the generation of the 3D model of a portion of a city located anywhere in the world (fully automatically). However, the geometry of the buildings and the terrain are still “separated”, as reported in Figure 11 by using different colours.

Figure 11 - Axonometric view of the terrain surface and building geometries

Two different procedures are implemented in the geo_mesher module to obtain the final numerical model.

The first method involves the creation of an initial coarse volumetric mesh in which the geometries of the buildings are inserted. Then, the distance field is calculated with the level-zero on the building's surface through an iterative process. Finally, the refinement process is executed with the MMG library (via API in Kratos).

Figure 12 - Refinement mesh procedure: (a) distance field step 1,
(b) distance field step 2, (c) distance field step 3, and (d) result

Figure 12 shows a small test carried out on a group of buildings to evaluate the effectiveness of the meshing process. The iterative refinement is performed to improve the quality of the mesh as shown from Figure 12a to Figure 12c. After some iterations, elements with negative distance are removed. The result is a body-fitted mesh (Figure 12d) where smaller elements are placed near the buildings while coarser elements are located away from them.

Good results are also obtained in a complex test case. Figure 13 depicts the result achieved on many buildings with complex shapes. The section on a building (Figure 14) underlines the different sizes of the elements (smaller near the surface of the building and larger away from the building).

Figure 13 - Large scale test with refinement process: axonometric view

Figure 14 - Large scale test with refinement process: section on a building

The procedure just described allows to obtain a good quality mesh but with very high computational efforts. This problem can be solved by adopting the parallel version of the used re-mesh tool (ParMMG). At the present stage, however, the parallel version is under development. For this reason, a different procedure is created and followed to generate a numerical model by combining terrain and building geometries.

The first step performed with this new procedure involves drilling the terrain surface with the building footprints (Figure 15). This operation allows the creation of new nodes on the terrain surface to insert the geometry of the buildings (Figure 16). Figure 17 shows the terrain surface before and after the "drilling" step focusing on the inserted nodes.

Figure 15 - Terrain surface with the building footprints

Figure 16 - Terrain surface with the building geometries

Figure 17 - (a) terrain surface before the "drilling" procedure with the new building
nodes highlighted, (b) terrain surface after the “drilling” procedure

The next step is the creation of the lateral and upper surface of the domain. Then, the minimum and maximum dimensions that the elements can reach are set for each surface. This step is schematized in Figure 18.

Figure 18 - Size of elements for each mesh. In parentheses the default values

Subsequently, the tetrahedral mesh is created through TetGen (available with the MeshPy library in Python). Figure 19 reports the final mesh of a test conducted on a few buildings. The domain has a diameter of 2.0 km and a height of 300.0 m.

Figure 19 - Tetrahedral mesh on a test case: axonometric view

Figure 20 - Tetrahedral mesh on a test case: detail view of buildings

As noted in Figure 19 and Figure 20, the element sizes change in different domain regions. It is due to the dimensions previously set.

To obtain the result reported below, it is necessary to provide only the coordinates of the center of the domain (42.442679 latitudes and 14.207077 longitudes), its diameter and height, and the diameter of the portion of the buildings.

Figure 21 - Numerical model generated with the developed application:
terrain and building surfaces (top view)

Figure 22 - Numerical model generated with the developed application:
boundary surfaces with the top surface in transparency (axonometric view)

Figure 23 - Numerical model generated with the developed application:
detail on building surfaces (axonometric view)

To create the volumetric mesh, the minimum and maximum dimensions of the elements for each region are set. This information is indicated in the table reported below:

Table 1 - Size of regions and meshes

where the Terrain inner is the central circle where the buildings are located, the Terrain middle is the intermediate annular portion without buildings, and the Terrain outer is the outer annular portion subject to the smoothing procedure.

A tetrahedral mesh is generated on the base of the parameters reported in Table 1. The total number of nodes is 1,379,229, the total number of triangles is 635,948, and one of tetrahedra is 7,764,938.

Triangular elements are the outer faces of the tetrahedra positioned on the external surfaces of the domain. On the triangular elements, the boundary conditions are appropriately set.

Figure 24, Figure 25 and Figure 26 depict, respectively, the axonometric view of a portion of the generated model, a zoom on the buildings and a detail on the building mesh. Lastly, a section is shown in Figure 27 where different elements' sizes relative to the buildings' surface can be seen.

Figure 24 - Tetrahedral mesh of the area under study generated with the
developed application: axonometric view

Figure 25 - Tetrahedral mesh of the area under study generated with the
developed application: axonometric view of the buildings

Figure 26 - Tetrahedral mesh of the area under study generated with the
developed application: detailed mesh on the buildings

Figure 27 - Tetrahedral mesh of the area under study generated with the
developed application: section on the buildings

Before assigning BCs, the lateral surface is divided into several sectors. The default value is 12, but different values can be set.

The BCs are assigned according to the settings reported in Table 2.

Table 2 - BCs assigned

The definition of the BCs is handled automatically by the geo_model module. Once the analysis to be performed is set (wind from North, East, South or West), the lateral surfaces are fixed as Velocity Inlet or Pressure Outlet, respectively.

The building and terrain surfaces are set as Wall No-Slip, while the top surface is defined as Wall Slip.

All the information about BCs is enclosed in a JSON file. An example is reported below:

"boundary_conditions_process_list": [
	{
		"Parameters": {
			"model_part_name": "FluidModelPart.Buildings"
		},
		"kratos_module": "KratosMultiphysics.FluidDynamicsApplication",
		"process_name": "ApplyNoSlipProcess",
		"python_module": "apply_noslip_process"
	},
	{
		"Parameters": {
			"model_part_name": "FluidModelPart.TopModelPart"
		},
		"kratos_module": "KratosMultiphysics.FluidDynamicsApplication",
		"process_name": "ApplySlipProcess",
		"python_module": "apply_slip_process"
	},
	{
		"Parameters": {
			"model_part_name": "FluidModelPart.BottomModelPart"
		},
		"kratos_module": "KratosMultiphysics.FluidDynamicsApplication",
		"process_name": " ApplyNoSlipProcess",
		"python_module": "apply_noslip_process"
	},
	
	...
]

The numerical model, obtained with the above-described application, is tested by performing a CFD simulation. In the present analysis, an incoming wind from the North direction is set, and the obtained results are depicted from Figure 28 to Figure 31.

The results show a very similar trend to those obtained with the microscale commercial software. Specifically, in Figure 28, the velocity magnitude result is shown; in Figure 29, the pressure field with streamlines is reported, in Figure 30, the velocity vectors on a longitudinal section on the reference building are depicted, and in Figure 31, the streamlines on the neighborhood under investigation are plotted.

Figure 28 - Top view of velocity magnitudes: incoming wind from North (-y)

Figure 29 - Top view of pressure levels and streamlines:
incoming wind from North (-y)

Figure 30 - Section view of velocity magnitudes on the reference building:
incoming wind from North (-y)

Figure 31 - Axonometric view of streamlines in the urban fabric:
incoming wind from North (-y)