Physics-aware wind surrogate model for UAV aerodynamic response assessment

In the coming years, the use of unmanned aerial vehicles (UAVs) to transport goods and people will become increasingly relevant in urban environments (Jordi et al., 2022). This rise in air mobility technologies has highlighted the potential of the urban air mobility (UAM) sector to revolutionize urban transportation, improve efficiency, reduce traffic and reduce carbon emissions (Garrow et al., 2021; Zhao and Feng, 2024).

However, implementing this new mode of mobility requires addressing safety issues, as UAVs can pose risks to pedestrians (Tepylo et al., 2023). Complex air phenomena, such as turbulent buildings wakes and vortex shedding, can compromise UAV stability (Reiche et al., 2021). Addressing these challenges requires advances in modelling, simulation and navigation systems, enabling UAVs to detect and adapt to environmental conditions in real time (Oo et al., 2023; Giersch et al., 2022; Tonti et al., 2024).

Computational fluid dynamics (CFD) simulations are widely used to analyse airflow dynamics across different locations. Although high-resolution models such as large eddy simulation (LES) are effective in optimizing urban flight paths (Rienecker et al., 2023), their high computational cost limits their use in real-time prediction tools. Thus, Reynolds-averaged Navier−Stokes (RANS) models (Sousa et al., 2018) are preferred, as they provide essential flight design data, including flow variables (e.g. velocity and pressure) and turbulence properties (e.g. kinetic energy and dissipation). These models simulate phenomena like boundary layer separations and flow instabilities (e.g. shear layers), which may compromise the UAV’s stability.

Even using simplified RANS models, their computational time makes them unfeasible for real-time predictions. To overcome this problem, machine learning (ML) techniques, particularly data-driven reduced-order models (ROM), efficiently extract information from CFD databases. Proper orthogonal decomposition (POD), an unsupervised ML technique, simplifies complex data sets by identifying patterns and structures (Brunton et al., 2020), and has been applied in fields such as non-linear structural and aerodynamic flow analysis (Ayoub et al., 2022; Xu et al., 2023). In the context of fluid dynamics, these dimensionality reduction models, in combination with other mathematical techniques, allow real-time predictions of aerodynamic problems (Moreno et al., 2016; Qin et al., 2019). Other recent techniques, like, for example, physics-informed neural networks (Cai et al., 2021), which provide a powerful framework for modelling complex, nonlinear dynamics (Gao et al., 2024) by embedding physical laws into the learning process, were discarded due to their training complexity and high computational needs, especially for large-scale urban simulations. On the contrary, POD offers a straightforward and robust approach for analysing and reconstructing urban wind flow fields, particularly when data quality is high, and the focus is on dominant flow structures.

Among the different algorithms used for dimensionality reduction, high order singular value decomposition (HOSVD) is a powerful and robust algorithm, particularly suited for breaking down high-dimensional tensors into simpler components (Brunton et al., 2020). By representing key patterns through orthogonal basis matrices and a core tensor that weights the influence of each basis, HOSVD enables effective reconstruction of the original tensor. If numerical interpolation techniques are combined with HOSVD, real-time predictions can be made for new cases not contained in the initial data set. Examples of use in different industrial applications can be found in (Viéitez et al., 2019; Martín et al., 2012), demonstrating the effectiveness of HOSVD.

Another limitation of CFD wind flow simulations (and the ROM developed from them) is that they are typically designed for a specific site, requiring manual construction of the computational domain and mesh generation, as well as user-defined properties (e.g. domain extension, terrain classification, building geometries). Most works focus on terrains with simple topography near domain boundaries, often flat ground, making it easier to implement atmospheric boundary layer (ABL) conditions (Hågbo et al., 2020; Toparlar et al., 2015). These constraints make the numerical tool impractical for flight path design in different urban environments. LiDAR point clouds can be used to reconstruct the terrain, but some studies omit buildings or generate them manually (Mirzaei, 2021; Kastner and Dogan, 2020). Another challenge is the treatment of vegetation, often modelled as porous media by adding resistance to airflow. Consequently, for simulations involving complex urban landscapes, it is essential to automate the classification of cells that are considered porous media. However, some studies that use LiDAR data face difficulties in accurately characterizing vegetation (Huo and Chen, 2023).

This study presents a ROM numerical tool (based on HOSVD) for real-time wind flows prediction that, using meteorological data as inputs, delivers precise air flow forecasts with a resolution of less than 1 m, achieving predictions 7200 times faster than traditional CFD simulations. This tool, combined with a state-of-the-art UAV simulator, is used to assess the flight stability of UAVs under different wind conditions. All steps needed in surrogate model generation, from the CFD geometrical domain and mesh generation, cells classification, numerical database generation and ROM model generation are automated, allowing the direct application of the tool to different locations.

The manuscript is organized as follows. Section 2 presents an overview of the study site (selected as an example of use of the tool in complex locations), the domain and mesh automatic generation tool, and the RANS model used to build the numerical database. Section 2 also presents the validation of the CFD model with experimental measurements extracted from four anemometers located at the study site. Section 3 explains the process for generating the reduced-order model in detail, while Section 4 shows examples of use of the generated tool to analyse the impact of meteorological conditions on the UAV’s flight path stability along specific flight trajectories in the study site. Finally, Section 5 discusses the relevant conclusions and proposes future directions of research.

The selected study site is located around the University of Vigo “Campus del agua” building^[1], at Ourense city, Spain. In this case, a fictional vertiport is proposed on this building. Figure 1(a) shows part of the vertiport area (left) around the study building (right), where the complexity of the terrain can be appreciated.

The computational domain was generated following the methodology introduced in previous work (Aldao et al., 2024), which presented an automatic 3D modelling approach using open-source data from the Spanish Geographic Institute (IGN) (Instituto Geográfico Nacional, 2024). Firstly, the terrain surface was generated using IGN’s MDT05, a digital elevation model (DEM) that maps the elevation of Spanish territory with a 5-m resolution. As shown in Figure 1(b), the city of Ourense has significant elevation changes that impact local wind patterns. To account for these effects in the simulations, a relatively large area was selected, covering a circle of radius R₁ = 1 km centred at the vertiport building location. To prevent numerical instabilities during CFD calculations and avoid flow recirculation near the boundaries, the computational domain was artificially extended, as depicted in Figure 1(c). The domain extensions consist of two parts: a cosenoidal transition surface, which ensures smooth, gradual altitude variation up to a constant level, and a completely flat external surface, which promotes orderly inflow and outflow and allows the development of atmospheric boundary layer (ABL) conditions. In addition, to account for the roughness properties of different ground surfaces, the terrain mesh was semantically classified using SIOSE AR, a vector-based land-use model provided by IGN, which categorizes the Spanish territory based on diverse data sources such as cadastral records, LiDAR point clouds and aerial imagery.

The building geometries were generated by extracting their footprint polygons from the SIOSE AR model and calculating volumes using base heights from the MDT05 model and rooftop heights (calculated from LiDAR data), modelling each building as a simple block. Only buildings within 500 m of the vertiport were modelled in detail, as they have the greatest impact on wind flow. This process led to the final 3D model. This process resulted in the final 3D model shown in Figure 1(d). The entire process of domain generation and surface faces characterization was done using open-software Python.

Finally, the mesh was generated in Ansys Fluent Meshing using a polyhexcore mesh, combining polyhedrons and hexahedrons for good orthogonality. After conducting a sensitivity analysis (Section 2.2.4), the final mesh comprised 34 millions of cells.

The finite volume method (FVM), which has been extensively used in the simulation of urban wind flows, was selected to solve the steady RANS equations (1) and (2) for incompressible flow inside the computational domain:

where 〈u〉 represents the Reynolds-averaged air velocity, µ_∞ corresponds to a constant ^[2] air dynamic viscosity, equal to 1.802 × 10⁻⁵kg/(m s), 〈p〉 denotes the averaged pressure generated due to the flow motion, ρ_∞ is the air density, considered constant (compressibility effects are negligible in this study), equal to 1.22 kg/m³ and S represents a momentum source term to apply inside cells classified as porous media.

To model the Reynolds stress tensor in equation (2), the standard κ − ε was chosen:

where $Cμ,Cε1,Cε2$⁠, σ_k and σ_ε are constants. Therefore, the Reynolds stress tensor $τ¯¯Re$ is modelled using the turbulent viscosity μ_t as follows:

where $I¯¯$ is the identity matrix.

This model was selected because of its effectiveness in predicting turbulence in steady-state scenarios. This makes it suitable for performing quick and accurate analyses of the Reynolds-averaged flow around the structure. The κ − ε model is widely used for these types of problem, either in its traditional formulation or in more recent adaptations. Furthermore, many studies highlight its robustness and the potential to improve its performance through a careful selection of coefficients (Shirzadi et al., 2017; Xiong et al., 2022).

To achieve a more realistic CFD analysis, as IGN has a vegetation LiDAR point cloud classification, a voxelization of the LiDAR point cloud was performed in a structured grid to classify the cells of the CFD mesh within each voxel, obtaining the mesh cells that correspond to the vegetation (only large trees, such as conifers, oaks and similar types, were modelled). A new cellZone was created to group all these cells for porous medium modelling.

Furthermore, this methodology was also applied to cells inside buildings that are not geometrically modelled (e.g. R > 500 m). As mentioned previously, only the buildings near the study area were geometrically represented (Liu et al., 2018); the remaining buildings were excluded to reduce computational costs (modelling more buildings generates smaller cells, resulting in more cells). Once the mesh cells corresponding to these buildings were filtered, a porous model was assigned to represent them as impermeable media. This approach ensured that the wind profile entering the campus area (the study area) was more realistic. As an example, the left plot of Figure 2 shows the terrain, while the right plot shows the addition of the cells classified as porous media. The purple cells represent buildings (red buildings are those that were geometrically modelled), while the dark green cells indicate vegetation.

All porous cells were modelled as simple homogeneous porous media, and a source term S is incorporated into the momentum conservation equation using the Darcy−Forchheimer model:

In this equation, D and F represent the laminar (Darcy term) and inertial loss (Forchheimer term) parameters. Non-geometrically modelled buildings are treated as impermeable media, with high coefficients of D = 2000 1/m² and F = 2000 1/m. For vegetation cells, only inertial losses are considered, with F calculated using F = 2, LAD, C_D, an expression obtained by combining the porous model with (SIMSCALE, 2024), where Leaf Area Density (LAD) is 1.25 [m²/m³] from the ORNL DAAC database (ORNL DAAC, 2014), and C_D is the drag coefficient set at 0.83 based on previous studies (Bekkers et al., 2022; Gonçalves et al., 2023). Thus, F = 2.075 1/m was assigned to all vegetation cells.

The boundary conditions for the surface boundaries indicated in Figure 1(d) are the following:

• Top: a symmetry condition was applied to the top boundary of the domain.

• Terrain and buildings: All geometrically modelled terrain and building surfaces were treated as walls. In addition, the transition region and the flat surface were treated in the same way. The atmospheric wall law for rough surfaces (Hargreaves and Wright, 2007) was applied to the surface-tangent average velocity at the centroid of every cell in contact with the boundary wall:

where d_c is the distance from the centroid of the cell to the boundary surface and z₀ is the roughness of the surface. Since the terrain was categorized in nine different types [see Figure 1(d), Figure 2], different z₀ were assigned based on Silva et al. (2007). Turbulence properties k_c and ε_c were established using the wall law, with a zero-gradient for turbulent kinetic energy and a dissipation rate defined as $εc=(u*)3κdc$⁠:

• Lateral: An inlet-outlet boundary condition was applied to the lateral surface of the computational domain, where a neutral ABL was imposed to define a horizontal velocity (Hargreaves and Wright, 2007; Yang et al., 2009), characterized by equations (9)–(11):

with C₁ = 0 and C₂ = 1 (Richards and Norris, 2019). In equation (9), u* represents the frictional velocity of the ABL in terms of the wind intensity u_ref at the reference height z_ref. For these simulations, a reference height of 10 m was chosen. Vector t_wind is the unit vector of the wind of the incoming air flow on the Lateral boundary, which can be expressed in terms of the orientation angle α: t_wind = (cos α, sin α, 0). The angle may vary from 0° (indicating the wind coming from the north) to 359°.

The model equations were solved using FVM with OpenFOAM (version 10) and the SIMPLE algorithm through the simpleFoam solver. Spatial discretization schemes were second-order accurate. The pressure was solved using the GAMG iterative method, while the remaining variables were solved with the smoothed Gauss−Seidel method. In addition, two non-orthogonal corrector loops were set, a uniform tolerance of 10⁻⁴ and relaxation factors were defined (0.8 for velocity and 0.3 for the other variables).

The computational domain and mesh generation for the 34 million-cell mesh took 1 h using 24 cores of an AMD EPYC 9454 processor, while the simulation took 2 h with 30 cores of the same processor.

A sensitivity analysis was performed to evaluate the impact of mesh resolution on the results using boundary conditions u_ref = 10 m/s and α = 0°. Three mesh sizes were tested: Mesh #1 (20 million cells), Mesh #2 (34 million cells) and Mesh #3 (41 million cells). Velocity at a probe located in a recirculation region past the building of interest was compared.

Mesh #1 and Mesh #2 showed a deviation from Mesh #3 of 6.1% and 1.3%, respectively. Mesh #1, Mesh #2 and Mesh #3 required 1.3, 2 and 7.5 h of computation, respectively. Mesh #2 achieved a good balance between accuracy and computational cost and was selected to build the numerical database.

A validation of the numerical model against the experimental measurements was carried out. Four different anemometers, one fixed at the top of the vertiport building [named A0, see Figure 3(a)], and three portable [named A1, A2 and A3, see Figure 3(b)] distributed throughout the campus site, as indicated in Figure 3(c), were used. The fixed anemometer is a 4-blade propeller probe ^[3], while the portable anemometers are two-dimensional ultrasonic probes ^[4].

The data collected by the anemometers were segmented into 1-h intervals to calculate the average wind speeds and directions for comparison with the CFD model. These measurements, from an experimental campaign on 11 June 2024, are summarized in Table 1(a). A ROM (explained in Section 3) was coupled to an optimization algorithm to identify the wind inlet conditions (u_ref and α) that minimized the error between the experimental and the numerical results at the A0-A3 anemometer locations. The optimization was guided by the following cost function that accounts for differences, both in magnitude and direction, of the measured and numerical wind speed vectors:

In equation (12), 〈u_i〉 and 〈v_i〉 are the numerical and measured mean velocities at the ith probe location, and h_i is the weight coefficient for the ith term. After a sensitivity analysis, the same weight h = 0.5 was assigned to all anemometers, ensuring balanced contributions from each, as well as equal weighting for wind intensity and direction measurements.

The method used to minimize the cost function f was the Nelder−Mead algorithm, which optimizes by iteratively refining a simplex of points in the search space. This method has also been applied effectively in the validation and optimization of CFD problems (Thompson et al., 2017). As a result, specific CFD wind flow boundary condition values (u_ref and α) were identified for each hour measurement, as shown in Table 1(a). Relative deviations between the CFD predictions at A0-A3 locations and the anemometers measurements are also indicated in Table 1(a). The averaged deviations for the wind velocity magnitude and direction are 8% and 5%, respectively.

Once the appropriate values for the inlet wind flow magnitude (u_ref) and direction (α) were identified, CFD simulations were performed using the identified inlet boundary conditions to predict wind velocities at two key locations [M0 and M1, see Figure 3(d)]. These points, at 126.8 m above ground level, correspond to locations where the 1-km resolution Weather Research and Forecasting (WRF) model from MeteoGalicia provides wind data (MeteoGalicia, 2024), allowing direct comparison of CFD and meteorological predictions for the day of the experiment. Table 1, shows the resulting absolute errors in wind speed magnitude and direction, respectively.

The mean deviation in magnitude of the CFD results [obtained from Table 1(b)] is 1.25 m/s, while the mean deviation in direction [Table 1(c)] is 20.78°. These values are consistent with the typical error margins observed in numerical weather prediction models. In particular, the WRF model has an estimated error of 1.6 m/s in wind intensity and 35° in wind direction (MeteoGalicia, 2024). Considering these benchmark values, our results fall within the expected range, confirming that the CFD simulations provide reliable and comparable estimates for both wind speed magnitude and direction.

A total of 360 simulations were performed to build the CFD database. These simulations varied the reference wind velocity u_ref and wind direction t_wind, with different values of u_ref ranging from 2 to 20 m/s and increments of 2 m/s. For each velocity, 36 different wind incidence angles α were tested, with increments of 10°.

After each simulation reached convergence, the solution was mapped to a coarser mesh with 2.2 million cells, as the original mesh contained very small cells near the walls (to ensure y⁺ < 400). However, since the UAV does not fit within these small boundary layers, the solution was transferred to a coarse mesh without boundary layer cells, consisting of 2.2 million cells, thereby decreasing the computation time of the predictive model without sacrificing essential information needed for UAV path planning.

To enhance computational efficiency and minimize calculation times, the simulations were parallelized into 30 subdomains and run on a high-performance computing station with 756 GB of RAM and two AMD EPYC 9454 processors (each with 48 cores running at 3.8 GHz). Three simulations were run simultaneously, and each simulation took approximately 2 h. The entire process, including mapping, took around 240 h to complete.

Before creating the reduced-order model, it was essential to process the numerical results for the 360 CFD simulations. For each simulation, each flow variable u_x, u_y, u_z, p, k and ε is stored in a list. The arrangement of the elements inside these lists reflects the values of the variables at the centroid of each cell within the mapped mesh. Since the same grid is used for all cases in the numerical database, the order of the cells remains consistent for all files.

A three-dimensional tensor is then build for each flow variable (namely, $Tux,Tuy,Tuz,Tp,Tk$ and T^ε) by piling the list according to the input velocity of each case and the angle of incidence. Therefore, the first index of the tensor is the spatial location of the cells, the second index corresponds to the inlet velocity parameter u_ref and the third index denotes the wind incidence angle parameter α. Consequently, the tensor dimensions for each variable is N_cells × 10 × 36, where N_cells is the number of cells of the mapped mesh (approx. 2.2 M), 10 is the number of simulated velocities and 36 is the number of simulated wind incidence angles. As an example, element $Tijkux$ indicates the velocity in the x direction at ith cell, for the jth inlet velocity u_ref and the kth angle of incidence α.

After obtaining the tensors for each variable, HOSVD, which extends the capabilities of the singular value decomposition (SVD) to multi-dimensional tensors, was applied to each tensor.

Decomposition of a tensor T of size I₁ × I₂ ×…× I_N is as follows:

where S is the tensor core, U⁽ⁱ⁾ are the mode matrices and ×_i represents the n-mode product (tensor product per matrix). As in the SVD, this expression can be expressed in terms of components:

where r_i denotes the rank of each matrix U⁽ⁱ⁾.

To compute U⁽ⁱ⁾ matrices, each mode’s unfolding matrix (denoted as Bⁱ) must be calculated. Computation of the unfolding matrices can be seen in (De Lathauwer et al., 2000). Then, SVD has to be performed for each of these Bⁱ matrices. In this way, the U⁽ⁱ⁾ matrices are obtained. Finally, to obtain the core tensor, the following calculation must be performed [see equation (15)]:

However, for this problem, it is not necessary to compute the full tensor core, S, to make predictions. Instead, it is sufficient to obtain the unfolded matrix of the tensor core associated with the first parameter. The matrix S₁ can be directly calculated using equation (16):

where B¹ is the first unfolding of the data tensor (T) and ⊗ denotes the Kronecker product.

To compress the tensor, a reduced set of modes, $r˜i$⁠, can be selected:

To evaluate approximation accuracy following HOSVD truncation, the relative root mean square error (RRMSE) is estimated by:

Here, N represents the number of parameters (in this case, N = 3) and $σ1n,…,σrnn$ are the singular values of each unfolding matrix, Bⁿ. In this study, the following accumulated energy AEn is used for each parameter n to determine the number of singular values (⁠$r˜n$⁠) that must be retained to preserve a specified percentage of information:

Once the HOSVD of the tensor has been calculated and an adequate truncation has been selected, it is necessary to perform interpolations on the columns of the modal matrices to reconstruct (e.g. to predict) the solution for a specific case not contained in the database.

In this study, the interpolations are performed on the columns of the truncated matrices. These truncated matrices are denoted as $U˜(i)$⁠. The dimensions of $U˜(1)$ are $Ncells×r˜1$⁠, the dimensions of $U˜(2)$ are $10×r˜2$⁠, while for $U˜(3)$ are $36×r˜3$⁠, where $r˜1<r1,r˜2<r2$⁠, and $r˜3<r3$⁠. In addition, the truncated core tensor, $S˜$⁠, will also be used for the predictions.

To reconstruct the solution for each cell of the mesh, the vector solution (named q_out), with N_cells entries, is computed using equation (20):

where S₁ represents the truncated matrix of the unfold matrix of S associated with the first parameter, calculated similarly to the matrix Bⁱ, while $iU˜(j)$ are the interpolated vectors. For this problem, $iU˜(2)$ and $iU˜(3)$ have dimensions $1×r˜2$ and $1×r˜3$⁠, respectively.

These interpolated vectors are computed using cubic 1D splines in each of the columns of $U˜(2)$ and $U˜(3)$⁠.

After developing the predictive model using HOSVD in combination with truncation and interpolation techniques, local real-time meteorological data are integrated as input of the predictive model to forecast air flow behaviour in the urban environment.

As the predictive model inputs are the wind speed and direction (at 10 m from the ground) at the artificial lateral boundary condition, while the available MeteoGalicia WRF forecast wind data ^[5] are located at the points M0 and M1 (at ground altitude of 126.8 m), it is then necessary to use the optimization tool described in Section 2.3 to select the appropriate ROM inputs that provide the smallest wind deviations at M0 and M1. Thus, a functional similar to (12) is defined for M0 and M1 locations, and its minimum value is found using the ROM tool coupled to the optimization Nelder−Mead algorithm.

In this study, an in-house Fortran-95 code is used to execute the HOSVD, truncation and interpolation processes, while the integration of the weather forecast predictions is done using a Python script.

The effect of HOSVD truncation on reconstruction accuracy is shown in Figure 4, which depicts the behaviour of the singular values and accumulated energy for each parameter. From Figure 4(a), achieving 99% of accumulated energy for the spatial parameter requires retaining only 79 modes out of 360, or 1/5 of the modes. For the other two parameters, wind intensity and wind direction, retaining 99% of information requires 3 (out of 10) and 28 (out of 36) modes, respectively, equivalent to 1/3 and 7/9 of the modes [see Figure 4(b) and 4(c)].

Table 2(a) summarizes the effect of retaining a given number of modes on accumulated energy, computation time for the predictions and the RRMSE error [given by equation (18)] introduced by HOSVD truncation for each parameter. It highlights the importance of truncating the number of modes, as this can significantly reduce the computation time without substantially increasing error.

In addition, a case not included in the original database was selected for prediction, allowing comparison with a CFD simulation. The selected case has an input reference velocity of u_ref = 4.5 m/s and α = 25°. Figure 5(a) shows 3D streamlines predicted by the surrogate model, while Figure 5(b) and 5(c), show a cross-sectional velocity view predicted by the CFD and the absolute error between the prediction of the ROM and the exact CFD solution. The absolute error is generally very small, with slightly higher errors in building wake regions. However, the prediction is quite accurate.

To evaluate the precision of the prediction of the surrogate mode, the mean absolute error (MAE) for the velocity magnitude ‖u‖ is presented in Table 2(b). The MAE, calculated across all grid cells, quantifies the mean deviation between the ROM predictions and the CFD results. It is defined as:

In addition, the 90% confidence interval for the absolute error is provided to estimate the range of error variation. This interval, defined by the minimum and maximum absolute error values (percentile limits), offers insight into the stability and accuracy of predictions for different numbers of retained modes in the HOSVD. When retaining 99% of the modes, the average error in velocity magnitude is 2.5 × 10^–2 m/s.

Therefore, for the final implementation of the ROM, 99% accumulated energy was selected, retaining the first [79, 3, 28] modes for space, wind velocity and wind direction, respectively. Across tested cases, the mean relative errors for ROM predictions of wind velocity in the vertiport region were 2%, compared to CFD results.

As a practical application in a UAM scenario, the impact of wind conditions on UAV behaviour was assessed. The purpose of these analyses is to quantify the impact of wind intensity and turbulence on the UAV stability and identify potential areas that may compromise the safety of operations. To this end, the predictions of the surrogate model were integrated into RotorPy, a state-of-the-art UAV simulator (Folk et al., 2023) as indicated in Figure 6.

Wind gust profiles are incorporated using a Dryden turbulence model, a widely accepted mathematical model of continuous gusts used in aircraft design and simulation. The model can be adjusted based on the variables 〈u〉, k and ε from the ROM model, replicating wind speed fluctuations through random sampling of power spectral density functions. As described in Figure 6, these wind fluctuations are integrated into the RotorPy UAV simulation module, which computes the aircraft dynamics using unsteady aerodynamic models. The simulator also replicates the flight controller response, which, based on state estimation derived from sensor measurements, computes the control actions necessary to maintain a given groundspeed and direction.

For this analysis, different test trajectories were generated, defined by a set of relevant turning points. Between these points, it is assumed that the aircraft follows a constant reference groundspeed v_g,ref, except during take-off and landing, where a constant acceleration or deceleration is applied to prevent abrupt manoeuvres, respectively. To analyse the aircraft’s behaviour at each point along the trajectory, a statistical motion assessment is conducted following this methodology.

As an example, the 650 m length trajectory displayed in Figure 7(a) was analysed following the previous methodology. The flight path was discretized at 1-m-length intervals, and the wind flow variables were computed at each trajectory point using the ROM prediction (the weather conditions for analysis were u_ref = 4 m/s and α = 240°). Figure 7(b) shows the turbulence kinetic energy and turbulence integral length scale at each point on the trajectory path. Then, wind speed fluctuations were simulated using the Dryden model, and RotorPy is used to calculate the response of the aircraft, considering the reference flight conditions (|v_g| = 6 m/s except for take-off and landing) at each discretization point. For each of these points, a statistical analysis was performed and the standard deviations in groundspeed and angular rates were computed for the given flight conditions [see Figure 7(c) and (d),]. For the calculations, AsTec Hummingbird is considered, a UAV model integrated within RotorPy. For further information about the aircraft, RotorPy’s documentation can be consulted (Folk et al., 2023).

For this example, high values of angular rate deviations are obtained during take-off and landing operations due to the relatively small integral turbulence length, which, in turn, increases the oscillation frequency, generating oscillations in the drone that may hinder the instability of the UAV flight. During the rest of the flight, these deviations remain smaller, but a region of relatively high angular rate deviations is encountered at half of the trajectory, generated again by the change in the wind turbulence values within the trajectory region.

An automatic numerical tool is presented for fast (1 s) and precise (under 1 m) wind flow predictions in complex urban environments. The generated tool, which efficiently exploits a (previously computed) CFD database using HOSVD techniques, uses the forecast wind predictions from the local meteorological agency as inputs. The tool has been validated against different results obtained from a CFD model (showing deviations of 2%), which, in turn, has been validated against experimental campaigns using anemometers located at a complex topographically vertiport site.

The surrogate model can run on standard computers, making costly high-performance workstations unnecessary (except maybe for the preprocessing of the CFD database) and enabling the incorporation of real-time meteorological predictions. This integration of real-time data, combined with the high precision of the surrogate model, makes the tool particularly suitable for planning drone routes in complex environments, such as cities.

The numerical tool presented here integrates an automated module for generating simulation domains directly from geospatial data, significantly enhancing its adaptability to any environment. This feature not only allows for the simulation of various terrains, but also facilitates a more versatile and general-use tool. The module can automatically categorize terrain to assign specific roughness values and classify volumetric cells to model far-located buildings and large vegetation areas as porous media. This flexibility extends to potential future applications, such as the generation of fast tools for the prediction of wildfire evolution in forested areas, e.g. by assigning different combustion properties depending on the characteristics of the terrain vegetation.

Predicting wind flows in urban settings is a challenging task with ample scope for future development. Although this work focused on neutral ABL conditions with no analysis of thermal effects (e.g. buoyancy plumes generated by isolation effects on the urban site), future work should include unstable and stable ABL conditions and analysis of thermal effects occurring, for example during day and night, respectively. As the inclusion of additional effects will also increase the number of parameters on which the problem depends, incorporating techniques such as Gappy Proper Orthogonal Decomposition (Gappy-POD) could further reduce the required CFD database size.

Funding: The authors would like to acknowledge the financial support for this work provided by grants PID2021-125060OB-100 and TED2021-129757B-C3, funded by MICIU/AEI/10.13039/501100011033 and the “European Union NextGenerationEU/PR TR”.