Using CFD to Evaluate Natural Ventilation through a 3D Parametric Modeling Approach

: Predicting building air change rates is a challenge for designers seeking to deal with natural ventilation, a more and more popular passive strategy. Among the methods available for this task, computational ﬂuid dynamics (CFD) appears the most compelling, in ascending use. However, CFD simulations require a range of settings and skills that inhibit its wide application. With the primary goal of providing a pragmatic CFD application to promote wind-driven ventilation assessments at the design phase, this paper presents a study that investigates natural ventilation integrating 3D parametric modeling and CFD. From pre- to post-processing, the workﬂow addresses all simulation steps: geometry and weather deﬁnition, including incident wind directions, a model set up, control, results’ edition, and visualization. Both indoor air velocities and air change rates (ACH) were calculated within the procedure, which used a test house and air measurements as a reference. The study explores alternatives in the 3D design platform’s frame to display and compute ACH and parametrically generate surfaces where air velocities are computed. The paper also discusses the effectiveness of the reference building’s natural ventilation by analyzing the CFD outputs. The proposed approach assists the practical use of CFD by designers, providing detailed information about the numerical model, as well as enabling the means to generate the cases, visualize, and post-process the results.


Introduction
Wind-driven natural ventilation is an attractive passive alternative in light of the challenges imposed by climate change and sustainable goals. To be effective, the strategy must deliver the required airflow quantities to satisfy both indoor air quality (IAQ) and thermal comfort criteria using outside air [1][2][3][4]. Nevertheless, natural ventilation is highly variable as it is driven by the climatic forces of wind (wind effect) and temperature (stack effect), which makes it challenging to estimate ventilation rates, an essential parameter for natural ventilation assessment. As the wind fluctuates and varies along the year, the investigations might verify if a specific place and time can provide the necessary quantity of air in a building/room. Different methods can be employed to address this task, including standards [5,6], guidelines [7,8], charts based on parametric analysis [9], empirical calculations [10], direct and indirect measurements [11], and building simulations [12][13][14][15]. An overview of these methods, with their respective outputs, can be found in [16,17].
With the advance of computing technology, the use of computational fluid dynamics (CFD) simulations has been widely increased to predict natural ventilation performance [18][19][20]. Wind-driven studies applying CFD have shown that the passive strategy can meet thermal comfort and reduce indoor temperatures and energy consumption by maximizing air velocity [21][22][23]. CFD is preferable to other methods due to its ability to address the natural ventilation phenomenon's complexity, sufficient validation with experimental data, and relatively low cost compared to the high-resolution generated data [24][25][26][27][28]. On the other hand, El Ahmar et al. [29] highlight that CFD simulation requires an extensive range of skills, being used primarily for research purposes and not much in practical applications as a widespread tool among architects and designers, who usually struggle to master it. Although accessible airflow simulation resources, as well as userfriendly approaches, have been developed [30,31], there are not many proposals regarding a framework that guides the process and combines parametric tools. When predicting natural ventilation within CFD, many parameters need to be considered, starting with the definition of the boundary conditions, followed by model set up and calculations and finalized by post-processing the simulation results. If designers wish to evaluate winddriven ventilation's potential and performance in either initial design or built environment through CFD, then many questions quickly arise.
In that light, this paper presents a study that assessed natural ventilation by integrating a multifaceted workflow and CFD building simulation. This research used existing programs and developed customized functions by benefiting from the growing interaction between building design and simulation through 3D parametric modeling platforms. The main objective of this paper is to provide practical use of CFD in wind-driven ventilation studies to promote CFD applications at the design phase through a structured procedure, guiding pre-to post-processing steps. The secondary aim is to provide different ways to treat and display the predicted air change rates (ACH) to exploit the available potentialities and alternatives for inspecting CFD outputs within the 3D environment.
All research steps are presented as follows: Section 1.1 reviews studies on wind-driven natural ventilation, and Section 1.2 provides an overview of how CFD models can be employed in these studies to predict airflows. Section 2 covers the workflow developed, presenting both the building and experimental protocol used as reference (2.2), as well as the taken steps to define the boundary conditions (2.3). A comprehensive description of the CFD model is detailed in Section 3, which addresses the domain, solver, mesh settings, and model verification. Section 4 compares different methodologies for post-processing, while Section 5 summarizes and discusses the study results. Finally, Section 6 compiles the research limitations, providing some remarks besides the pros and cons of running CFD within 3D parametric modeling, followed by a conclusion in Section 7.
This study's content should pave the way to promote CFD's application in winddriven natural ventilation studies through tool integration in parametric design platforms. The approach can assist CFD's use among non-experts, facilitating case generation and providing new ways to visualize and compute the results.

State of the Art in Wind-Driven Natural Ventilation Investigations
Natural ventilation occurs due to pressure differences in the openings caused by wind force, buoyance, or combining these two. When investigations are not on high building, and the wind forces can be considered sufficiently strong, the buoyance effect becomes negligible, and wind-driven natural ventilation can be the focus. Consequently, several studies are concentrated on measuring, evaluating, and predicting this passive strategy, especially in recent years because of sustainable and energy demand goals [16,32]. Some emphasize either cross-or single-sided ventilation, investigating specific topics such as the impact of the building dimensions [33], internal divisions [34], multiple windows [35], sheltering [36,37], adjacent obstacles [38,39], wind exchangers [40], large openings [41], and a greenhouse with a wind tower [42]. Moreover, Dai et al. [43] study air quality and virus propagation in interunit dispersion in multistory buildings. Experimentally, campaigns to assess wind-driven natural ventilation have been conducted on a real scale and with wind tunnels. They comprise pressure coefficients measured over facades for either isolated configurations [44] or in an urban context [45]. Lo and Novoselac [46] provide measurements in a multi-zone test building that, besides façade pressures, also includes wind properties, airflow through small window openings, and tracer-gas concentrations. A church typology is considered by Hayati et al. [47], while Tecle et al. [48] evaluate a low-rise building with different openings sizes and locations, room compartmentalization, and wind directions, measuring both ventilation rate and discharge coefficients. Extensive measurements for cross-ventilation in a sheltered building are presented by Shirzadi et al. [49], which can be used for the validation of CFD numerical models.

Assessing Natural Ventilation in CFD
Airflow prediction in CFD depends upon boundary/initial conditions as well on the turbulence model and the calculation method employed in the simulations. The three main turbulence models from lower to higher accuracy and computational cost are as follows: (i) Reynolds-averaged Navier-Stokes (RANS); (ii) large eddy simulation (LES); and (iii) direct numerical simulation (DNS) [50]. Due to the currently available computer hardware, DNS is still an unfeasible method for simulating natural ventilation in large-scale studies such as buildings; therefore, most of the studies adopt either RANS or LES models. Studies focused on evaluating turbulence models' performance, potentially suitable for indoor airflow calculations in terms of precision and computing demand, were investigated in References [51,52].
Of the calculation methods utilized, the two most popular are the integration of the velocity profile in the opening and the tracer-gas decay methods. The former is the simplest and widely-used method, [53,54], suitable only for incompressible (stable) flows, which is the case in most natural ventilation studies. It integrates the average normal velocities coming through an orifice (opening) with its area, and the number of velocities considered in the calculation corresponds to the number of small regions that constitute the opening mesh. Although the method applies to both RANS and LES simulations, an alternative, only valid for LES models, computes the instantaneous airflow rate at each time step by integrating the normal velocities at the openings and then time-averages the instantaneous airflow rates [55,56].
As for the tracer-gas decay method, the CFD calculates ventilation rates reproducing the tracer-gas decay technique used experimentally [57,58], and the dispersion of the gas tracer is modeled by introducing a further transport equation. While the method can overestimate airflow exchange rate in a numerical context since it considers all ventilation mechanisms [59], the tracer-gas technic is considered to be more reliable than the integration of the opening velocities method. Still, its high computational expense might not justify its relatively modest improvement of accuracy [59]. Moreover, the integration method is most precise when the incident wind flow is perpendicular to the opening. Otherwise, when the wind comes at an angle, the tracer-gas approach might be preferable [60].
On the other hand, airflow rates relative to the volume of the space may be preferable to the absolute ventilation flow, given in m 3 /s or l/s. In that sense, air change rate (ACH) is commonly employed as a target criterion to assess natural ventilation regarding acceptable indoor air quality and comfort. The measure is often used as a "rule of thumb" in ventilation design for building engineers [61] and, therefore, configures a desirable output.
Wind-driven CFD models can estimate ACH through the calculation of the airflow rate (Q r ) in m 3 /s by Equation (1) and later convert it to ACH by Equation (2). where V is the space volume (m 3 ); the subscript i represents the number of cells that make up the vertical plane across the opening; → v i is the mean velocity vector; → n i is the normal direction to the opening; A i is the corresponding opening area of cell i; and j is the opening j of the building.
Furthermore, indoor air velocity (m/s) is another considered measure when assessing natural ventilation [62][63][64]. Increasing air movement has the potential to shift thermal acceptability to higher operative temperature values, as demonstrated in research in hot and humid climates [65][66][67]. The adaptative comfort standard ASHRAE 55 [68] also addresses this wind cooling effect and foresees an increase in the average airspeed depending on temperature, limited to 1.2 m/s. Similarly, the European version (EN 16798-1) [69] states that values up to 1.5 m/s could extend the comfort limits in the spaces by 4 • C. Meanwhile, transient studies that investigated the wind effect in summer periods have also proven its effectiveness, using averaged air velocities calculated in horizontal planes across the spaces as an assessment method [70].
Thus, many measures are available to evaluate the performance of natural ventilation through CFD simulation. One method's choice over another depends on the investigation's approach and objectives, aside from the available initial data to feed the necessary conditions to build up an adequate model.

Overview
The method procedure consists of two main activities. First, geometry and location under study are defined (pre-processing), meaning the number of incident winds (angles and velocities) to be simulated is settled based on weather data analyses. Section 2.3 explains this process in greater detail. Second, calculation and post-processing steps are repeated for each of the investigated wind angles (procedure detailed in Sections 3 and 4).
The schematic diagram in Figure 1 compiles all steps proposed by this structured workflow. As illustrated, this study uses multiple software programs for pre-processing, simulations, results post-processing, and visualization. Most of the work was done within the commercial 3D modeler Rhinoceros 5 or Rhino [71] through one of its widely adopted graphical algorithm editor, Grasshopper [72], and its many plugins. Although numerous architectural design and simulation studies utilized these tools' power function [73][74][75][76], there is still a lack of connectivity between the parametric design platform and the building performance simulation, especially for natural ventilation assessment. Hence, a customized workflow was created to enable the automatization of the CFD simulations, together with a new alternative for post-processing and visualizing the airflow simulation results using Grasshopper native components. Rhino 5 was used for the geometric model because Grasshopper is settled for parametrization, and Butterfly, a Ladybug Tool [77], is used as an object-oriented python library that creates and runs CFD simulations using OpenFOAM. Different wind directions (building rotation angles) and other geometric configurations such as wind tunnel dimensions and object functions can be parameterized to generate their respective CFD cases and post-processing surfaces automatically. Additionally, the Supplementary Material used through these research steps is available for download in the Mendeley repository, and the link is available in the Supplementary Materials section in Back Matter of this article. Steps to assess wind-driven natural ventilation through computational fluid dynamics (CF)D within parametric 3D modeling.

Reference Building and Wind Velocities Measurements
The case-study building concerns one of the full-scale passive test houses from the INCAS experimental platform at the French National Institute for Solar Energy-INES facility, located near Chambéry, France (45°38′38.5″ N 5°52′27.4″ E). The edifice is used as a reference in this research because one of its goals is to support numerical simulations. Therefore, it is equipped with various sensors and has a simple design to ensure a more straightforward numerical verification process. The investigated building, also known as I-MA house, is a two-story rectangular cavity brick construction (7.5m × 8.5m), with the most massive openings facing south (34% glazed). Table 1 provides an overview of the investigated spaces within this study, while Figures 5b,c and 7 show the geometry of the house. The building floor plan is provided in the Appendix A ( Figure A1).  Steps to assess wind-driven natural ventilation through computational fluid dynamics (CF)D within parametric 3D modeling.

Reference Building and Wind Velocities Measurements
The case-study building concerns one of the full-scale passive test houses from the INCAS experimental platform at the French National Institute for Solar Energy-INES facility, located near Chambéry, France (45 • 38 38.5" N 5 • 52 27.4" E). The edifice is used as a reference in this research because one of its goals is to support numerical simulations. Therefore, it is equipped with various sensors and has a simple design to ensure a more straightforward numerical verification process. The investigated building, also known as I-MA house, is a two-story rectangular cavity brick construction (7.5 m × 8.5 m), with the most massive openings facing south (34% glazed). Table 1 provides an overview of the investigated spaces within this study, while Figures 5b,c and 7 show the geometry of the house. The building floor plan is provided in the Appendix A ( Figure A1). During a one-week experimental campaign (19-25.08.14) that assessed the building's thermal inertia, indoor air velocity data (m/s) at a 1.2 m height were recorded with an anemometer (DeltaOhm-accuracy ± 0.1 m/s) in the living room (sensor position is highlighted in Figure 9). The sampling rate was one recording per minute ( Figure 2). Although desirable, the experiment was unable to measure air velocities in other rooms or the building's ventilation rate due to resource and time limitations. During a one-week experimental campaign (19-25.08.14) that assessed the building's thermal inertia, indoor air velocity data (m/s) at a 1.2 m height were recorded with an anemometer (DeltaOhm-accuracy ± 0.1 m/s) in the living room (sensor position is highlighted in Figure 9). The sampling rate was one recording per minute ( Figure 2). Although desirable, the experiment was unable to measure air velocities in other rooms or the building's ventilation rate due to resource and time limitations. The measurement protocol involved fully opened internal doors, with one window in each orientation per floor tilted and opened during the night (9:00 p.m.-6:59 a.m.). Furthermore, the underfloor heating system was off, and the building ventilation system (noload supply airflow rate) was on for the whole period with approximately 135 m³/h flow rate, app. 0.5 h −1 . This building geometry, as well as its monitored data, was set as a reference for analyzing the outputs of CFD simulations, which used the climate file recorded on-site to perform the wind analyses and thus determine the boundary conditions. Alt- The measurement protocol involved fully opened internal doors, with one window in each orientation per floor tilted and opened during the night (9:00 p.m.-6:59 a.m.). Furthermore, the underfloor heating system was off, and the building ventilation system (no-load supply airflow rate) was on for the whole period with approximately 135 m 3 /h flow rate, app. 0.5 h −1 . This building geometry, as well as its monitored data, was set as a reference for analyzing the outputs of CFD simulations, which used the climate file recorded on-site to perform the wind analyses and thus determine the boundary conditions. Although the measurements are insufficient to carry out a rigorous numerical model validation, they provide proper orientation to guide the analyses concerning natural ventilation performance. A detailed description of the construction, as well as the climatic data used in the study, the experimental protocol, and the measurement equipment, can be found in Reference [78].

Wind Data Input for CFD Simulations
As wind direction fluctuates and significantly influences ventilation performance, experiments and simulations must consider its variation [79]. Thus, when defining the CFD boundary conditions for specific investigations/locations, it is essential to analyze the available wind data to outline the velocities and incident winds that are addressed at the numerical assessments. Considering this, the frequency distribution and wind speeds were initially extracted from the considered meteorological data. This analysis runs through the wind rose component from the parametric environmental plugin Ladybug [80], together with native Rhino/Grasshopper elements.
The component takes an EnergyPlus Weather file (.epw) as input, and how many cardinal points dividing the hourly average wind speed data must be defined. One needs at least four cardinal directions, and 36 angles (10 • interval) would represent a high resolution. Data distribution between 30 • to 40 • is recommended in [81] to establish a direct relationship between the available wind angles and the best building orientation that could benefit from the winds. Thus, 12 cardinal directions were considered in this investigation. As output, the component gives a list containing the time-frequency that the wind is coming from for each direction and the average wind velocity. However, a logarithmic normal distribution is considered to better describe the probability distribution of the wind speed frequency curves [82,83]. Therefore, both normal and lognormal distributions were used ( Figure 3) when analyzing the wind speed statistical distribution, and the best fit was selected. Each of these dominant incident wind directions corresponds to a building rotation angle that, together with its respective wind speed, configured a unique CFD model. in the study, the experimental protocol, and the measurement equipment, can be found in Reference [78].

Wind Data Input for CFD Simulations
As wind direction fluctuates and significantly influences ventilation performance, experiments and simulations must consider its variation [79]. Thus, when defining the CFD boundary conditions for specific investigations/locations, it is essential to analyze the available wind data to outline the velocities and incident winds that are addressed at the numerical assessments. Considering this, the frequency distribution and wind speeds were initially extracted from the considered meteorological data. This analysis runs through the wind rose component from the parametric environmental plugin Ladybug [80], together with native Rhino/Grasshopper elements.
The component takes an EnergyPlus Weather file (.epw) as input, and how many cardinal points dividing the hourly average wind speed data must be defined. One needs at least four cardinal directions, and 36 angles (10° interval) would represent a high resolution. Data distribution between 30° to 40° is recommended in [81] to establish a direct relationship between the available wind angles and the best building orientation that could benefit from the winds. Thus, 12 cardinal directions were considered in this investigation. As output, the component gives a list containing the time-frequency that the wind is coming from for each direction and the average wind velocity. However, a logarithmic normal distribution is considered to better describe the probability distribution of the wind speed frequency curves [82,83]. Therefore, both normal and lognormal distributions were used ( Figure 3) when analyzing the wind speed statistical distribution, and the best fit was selected. Each of these dominant incident wind directions corresponds to a building rotation angle that, together with its respective wind speed, configured a unique CFD model.  Besides the wind data from the experimental campaign (19-25.08.14), the investigations also considered the cooling season, identified from April to September [78], as a boundary condition. While the former compares simulated to measured data, the latter aims to take into account the period in which natural ventilation could be used as a cooling strategy, thus configuring a more representative analysis. In this sense, a statement condition function was applied to the Ladybug wind rose component for the cooling season, restricting the wind data assortment to the period when the dry-bulb temperature (DBT) was higher than 24 °C. Besides the wind data from the experimental campaign (19-25.08.14), the investigations also considered the cooling season, identified from April to September [78], as a boundary condition. While the former compares simulated to measured data, the latter aims to take into account the period in which natural ventilation could be used as a cooling strategy, thus configuring a more representative analysis. In this sense, a statement condition function was applied to the Ladybug wind rose component for the cooling season, restricting the wind data assortment to the period when the dry-bulb temperature (DBT) was higher than 24 • C. Figure 4 gathers the wind roses generated for both investigated scenarios, showing their respective wind speed and frequency values for each considered direction. Most of the wind speed distribution follows the logarithmic normal distribution, except for the southwest wind (θ = 240 • ) in the measurement week (Scenario 1) and θ = 60 • , θ = 90 • , and θ = 120 • in the cooling season (Scenario 2). Figure 3b illustrates this exception, showing that the normal distribution best describes the east wind (θ = 90 • ) of the cooling season. In most cases, however, the lognormal distribution provides a better fit, as shown in Figure 3a, which illustrates the data of the south wind (θ = 180 • ) recorded in the measurement week.
Energies 2021, 14, x FOR PEER REVIEW 8 of 27 Figure 4 gathers the wind roses generated for both investigated scenarios, showing their respective wind speed and frequency values for each considered direction. Most of the wind speed distribution follows the logarithmic normal distribution, except for the southwest wind (θ = 240°) in the measurement week (Scenario 1) and θ = 60°, θ = 90°, and θ = 120° in the cooling season (Scenario 2). Figure 3b illustrates this exception, showing that the normal distribution best describes the east wind (θ = 90°) of the cooling season. In most cases, however, the lognormal distribution provides a better fit, as shown in Figure 3a, which illustrates the data of the south wind (θ = 180°) recorded in the measurement week.  As 12 directions are considered for each scenario, 24 CFD models are generated. Nevertheless, as this represents a high computational cost, the study was restricted to the most occurring wind directions. Therefore, the wind angles with frequencies less than 1.5% in  As 12 directions are considered for each scenario, 24 CFD models are generated. Nevertheless, as this represents a high computational cost, the study was restricted to the most occurring wind directions. Therefore, the wind angles with frequencies less than 1.5% in the cooling period and 5% in the measurement week were not simulated, totaling 18 CFD models. The investigated angles (incident wind directions) in the study are highlighted in Figure 4.

Model Configuration
The CFD open-source code Field Operation and Manipulation OpenFOAM software [84], version 2006, was used to run a steady airflow through a three-dimensional model. As the study focused on wind-driven ventilation, the effects of thermal-driven ventilation (buoyancy) were not considered, meaning the simulations run under isothermal conditions. OpenFOAM has an extensive range of features to solve complex fluid flows and has been validated with both on-site measurements and wind-tunnel experimental data [85,86].
Conventionally, the simulation process takes place at the Grasshopper/Rhino interface, and the CFD engine runs in the background. This procedure, despite systematizing the calculation, has a downside: the program remains in processing until all simulation steps are finished. Therefore, considering the number of cases in this study, it was preferred to run them directly on the OpenFOAM terminal, using the Linux operating system. The use of Butterfly components was restricted to writing CFD model files, and instead of using command lines to execute the simulation steps, a script was developed to automatize the process.

Domain
Two geometries were developed for the study ( Figure 5). One uses the window configurations from the experimental campaign to enable comparison between simulated and measured data (CFD 1- Figure 5b). As opposed to the cooling season model, all openings were left completely open, without frames (CFD 2- Figure 5c). This approach provides the maximum airflow rates in the spaces to assess if natural ventilation may be employed as a cooling strategy.
Energies 2021, 14, x FOR PEER REVIEW 10 of 27 algorithm) [95]. The under-relaxation factors (α, 0 < α ≤ 1) were respectively set to 0.3 to pressure p, and 0.7 to velocity U, turbulent kinetic energy k, and dissipation rate ε. The under-relaxation technique is employed to improve computational stability, limiting the amount that the variable changes from one iteration to the next [96]. Second-order discretization schemes were used for both convection and diffusion terms of the mathematical governing equations. The selected numerical schemes set for these terms were correspondingly bounded Gauss linearUpwind [97] and the Gauss linear limited corrected, where the explicit non-orthogonal correction was set to Ψ = 0.333.
Both residuals values and monitoring results at specific locations (probes) were set as convergence criteria. For the latter, points located at the window surface area were selected, and the velocity and pressure outputs in those points were checked. The solutions converted approximately at the maximum residual level of 10 −5 with solution imbalances of less than 0.5%. All CFD simulations were run using a desktop with 20 Cores and 288 GB of RAM memory. The computational time for each incident wind was approximately 7 min.

Mesh processing and boundary conditions
The accuracy of the CFD results directly depends on the mesh quality, which was generated by the snappyHexMesh and the blockMesh utilities of OpenFOAM. First, block-Mesh creates a background mesh that defines the domain extensions and a base level mesh. Afterward, snappyHexMesh constructs a three-dimensional polyhedral mesh around the investigated geometry described by a tri-surface mesh made of hexahedra cells. The computational grid was fully structured, being more refined near the area of interest, the immediate surroundings of the building model. Minimum and maximum layer thickness were 0.01 and 0.3, respectively, with an expansion ratio of 1.1, guaranteeing a smooth transition from fine mesh near the wall surfaces and avoiding significant aspect ratios that can compromise convergence. The resulting grid and mesh for both CFD models are shown in Figure 5b,c, which had four refinement levels and three cells between levels. Although the snappyHexMeshDict already contains numerous quality control settings, the mesh quality was also visually verified in ParaView [98] version 5.8.0, an opensource, multi-platform data analysis for visualization applications; and statistically assessed with the checkMesh utility of OpenFOAM. Other quality controls, such as maximum non-orthogonality and feature angles were also changed, set respectively to 50 and varying form 70-90, considering each of the incident wind direction simulated.
(a) Moreover, grid sensitivity studies were conducted to ensure that the selected grid results did not vary considerably with finer grid resolutions. Therefore, coarse, medium, and fine meshes were created for CFD model 2 considering the north wind direction (θ = 0°) by changing the number of refinement levels. Consequently, the ratio between the smallest and largest mesh is beyond a conventional range, although the differences in results among them differ to the third decimal place. Table 2 presents the total number of cells, model residuals, and the calculated mean scalar velocity values for the I-MA first floor. As the variance between the coarse and fine grid is close to 1.5% and between the medium and fine mesh is below 1%, all simulations utilized the medium mesh properties. Regarding the boundary conditions, vertical profiles for the mean velocity U, turbulent kinetic energy k, and dissipation rate ε were assigned at the inlet through the atmospheric boundary layer (ABL) function [99]. As for the wind velocity, applied normal to the inlet, the values were set according to the wind data selection described previously in Section 2.3. With the ABL function, the roughness class and length (Z0) of the area around the reference building is taken into account, where Z0 =0.25 corresponds to an area with scattered obstacles, e.g., low buildings [100], which is the case for the experimental building too.
Furthermore, standard wall functions were assigned for the ground and building surfaces, and the near-wall velocity profiles were modeled so that a transition between the fully turbulent area and the region near the wall is considered due to the turbulence models requiring complementary information in the near-wall region [101]. The non-dimensional parameter y+ that describes the treatment near the walls was measured in 58, which is within the best practice interval from 30 to 300 for k-ε turbulence model [102,103]. A zero-static pressure (fixedValue) was applied at the outlet plane, also assigned with a velocity inlet/outlet type. Zero normal velocities and zero-gradient pressures were set to the ground and building surfaces, while the top and lateral sides of the domain were set as slip, meaning the viscous effects at the wall are negligible on those regions. The model setup encompasses a coupled indoor-outdoor CFD simulation, where the same computational domain captures the dynamic interaction between external and internal environments around openings [61,87,88]. Based on the best practices guidelines [89], the computational domain size was determined: windward, height, and lateral sides = 5H, and leeward = 15 H to allow flow redevelopment; where H is the building height (8.5 m). As a result, the computational domain ( Figure 5a) had the following measurements: length × width × height, of 183 × 98 × 42 m 3 , resulting in a blockage ratio of 1.8%, which meets the recommended values by previous studies of less than 3%, to avoid the effect of compressed flow [89,90].

Solver Settings
We considered the flow as turbulent; therefore, a turbulence model, RANS equations [91], was employed for determining the wind pressure variations on the I-MA house, solved in combination with the standard k-ε turbulence model. Numerous studies investigated different CFD turbulence models to predict airflow in naturally ventilated buildings and revealed that the two-equation k-ε turbulence model is optimal from a performance/cost perspective compared to other models [52,[92][93][94]. The k-ε turbulence model is one of the most commonly used two-equations turbulence model in CFD for natural ventilation problems, and it is described as stable, giving reasonably accurate results for most indoor airflows. However, it may have difficulties dealing with special room situations (e.g., high-buoyancy effect and large temperature gradient [51]), which is not the subject of this research.
For the velocity-pressure coupling, a steady-state solver for incompressible and turbulent flows was used-the semi-implicit method for pressure-linked equations (SIMPLE algorithm) [95]. The under-relaxation factors (α, 0 < α ≤ 1) were respectively set to 0.3 to pressure p, and 0.7 to velocity U, turbulent kinetic energy k, and dissipation rate ε. The under-relaxation technique is employed to improve computational stability, limiting the amount that the variable changes from one iteration to the next [96]. Second-order discretization schemes were used for both convection and diffusion terms of the mathematical governing equations. The selected numerical schemes set for these terms were correspondingly bounded Gauss linearUpwind [97] and the Gauss linear limited corrected, where the explicit non-orthogonal correction was set to Ψ = 0.333.
Both residuals values and monitoring results at specific locations (probes) were set as convergence criteria. For the latter, points located at the window surface area were selected, and the velocity and pressure outputs in those points were checked. The solutions converted approximately at the maximum residual level of 10 −5 with solution imbalances of less than 0.5%. All CFD simulations were run using a desktop with 20 Cores and 288 GB of RAM memory. The computational time for each incident wind was approximately 7 min.

Mesh Processing and Boundary Conditions
The accuracy of the CFD results directly depends on the mesh quality, which was generated by the snappyHexMesh and the blockMesh utilities of OpenFOAM. First, blockMesh creates a background mesh that defines the domain extensions and a base level mesh. Afterward, snappyHexMesh constructs a three-dimensional polyhedral mesh around the investigated geometry described by a tri-surface mesh made of hexahedra cells. The computational grid was fully structured, being more refined near the area of interest, the immediate surroundings of the building model. Minimum and maximum layer thickness were 0.01 and 0.3, respectively, with an expansion ratio of 1.1, guaranteeing a smooth transition from fine mesh near the wall surfaces and avoiding significant aspect ratios that can compromise convergence. The resulting grid and mesh for both CFD models are shown in Figure 5b,c, which had four refinement levels and three cells between levels. Although the snappy-HexMeshDict already contains numerous quality control settings, the mesh quality was also visually verified in ParaView [98] version 5.8.0, an open-source, multi-platform data analysis for visualization applications; and statistically assessed with the checkMesh utility of OpenFOAM. Other quality controls, such as maximum non-orthogonality and feature angles were also changed, set respectively to 50 and varying form 70-90, considering each of the incident wind direction simulated.
Moreover, grid sensitivity studies were conducted to ensure that the selected grid results did not vary considerably with finer grid resolutions. Therefore, coarse, medium, and fine meshes were created for CFD model 2 considering the north wind direction (θ = 0 • ) by changing the number of refinement levels. Consequently, the ratio between the smallest and largest mesh is beyond a conventional range, although the differences in results among them differ to the third decimal place. Table 2 presents the total number of cells, model residuals, and the calculated mean scalar velocity values for the I-MA first floor. As the variance between the coarse and fine grid is close to 1.5% and between the medium and fine mesh is below 1%, all simulations utilized the medium mesh properties. Regarding the boundary conditions, vertical profiles for the mean velocity U, turbulent kinetic energy k, and dissipation rate ε were assigned at the inlet through the atmospheric boundary layer (ABL) function [99]. As for the wind velocity, applied normal to the inlet, the values were set according to the wind data selection described previously in Section 2.3. With the ABL function, the roughness class and length (Z 0 ) of the area around the reference building is taken into account, where Z 0 = 0.25 corresponds to an area with scattered obstacles, e.g., low buildings [100], which is the case for the experimental building too.
Furthermore, standard wall functions were assigned for the ground and building surfaces, and the near-wall velocity profiles were modeled so that a transition between the fully turbulent area and the region near the wall is considered due to the turbulence models requiring complementary information in the near-wall region [101]. The nondimensional parameter y+ that describes the treatment near the walls was measured in 58, which is within the best practice interval from 30 to 300 for k-ε turbulence model [102,103]. A zero-static pressure (fixedValue) was applied at the outlet plane, also assigned with a velocity inlet/outlet type. Zero normal velocities and zero-gradient pressures were set to the ground and building surfaces, while the top and lateral sides of the domain were set as slip, meaning the viscous effects at the wall are negligible on those regions.

Model Verification
To demonstrate the accuracy of the numerical simulations, velocity components predicted within the RANS model were compared against experimental results. The flow condition and building dimensions followed the wind tunnel experiment from Karava et al. [104], which used particle image velocimetry (PIV) to measure the velocity field inside a single zone. The full-scale model width, length, and height equal 20 × 20 × 16 m, respectively, with one opening (9.2 m × 3.6 m) at the center of the inlet and outlet facades (Configuration E1 [104]), modeled in 1:2 scale. All other settings of the CFD model were maintained as described in the previous sections. Figure 6 shows the simulation results and the experimental data of the normalized velocity U/U ref distribution on a horizontal measurement line (L). The square symbols represent the PIV measurements of Karava et al. [104], and the x symbols are the simulation results of Ramponi and Blocken [87] with shear-stress transport (SST) k-ω model (their reference case). The mean absolute deviation between simulation and experimental data is 0.090 for the whole section but drops to 0.047 if restricted to the readings inside the investigated geometry. The comparison shows that Ux results agree with data from the reference study, and CFD can reproduce the velocity tendency.

Model Verification
To demonstrate the accuracy of the numerical simulations, velocity components predicted within the RANS model were compared against experimental results. The flow condition and building dimensions followed the wind tunnel experiment from Karava et al. [104], which used particle image velocimetry (PIV) to measure the velocity field inside a single zone. The full-scale model width, length, and height equal 20 x 20 x 16 m, respectively, with one opening (9.2 m x 3.6 m) at the center of the inlet and outlet facades (Configuration E1 [104]), modeled in 1:2 scale. All other settings of the CFD model were maintained as described in the previous sections. Figure 6 shows the simulation results and the experimental data of the normalized velocity U/Uref distribution on a horizontal measurement line (L). The square symbols represent the PIV measurements of Karava et al. [104], and the x symbols are the simulation results of Ramponi and Blocken [87] with shear-stress transport (SST) kω model (their reference case). The mean absolute deviation between simulation and experimental data is 0.090 for the whole section but drops to 0.047 if restricted to the readings inside the investigated geometry. The comparison shows that Ux results agree with data from the reference study, and CFD can reproduce the velocity tendency.

Post-Processing
Natural ventilation is evaluated in three ways in this study: (i) average air velocities (m/s) in horizontal planes across the spaces; (ii) ACH, using either all velocity vectors that pass through the openings; or (iii) filtering recirculation of return air with the assistance of the parametric design platform, which also displays the velocity vectors.

Air Velocities
After the case has run and converged, OpenFOAM surface field value object functions are added to the system file, and the case runs once again so that average indoor air velocity (m/s) across all spaces are calculated at 1.2 m height horizontal plans (.stl files).

Post-Processing
Natural ventilation is evaluated in three ways in this study: (i) average air velocities (m/s) in horizontal planes across the spaces; (ii) ACH, using either all velocity vectors that pass through the openings; or (iii) filtering recirculation of return air with the assistance of the parametric design platform, which also displays the velocity vectors.

Air Velocities
After the case has run and converged, OpenFOAM surface field value object functions are added to the system file, and the case runs once again so that average indoor air velocity (m/s) across all spaces are calculated at 1.2 m height horizontal plans (.stl files). These surfaces, shown in yellow in Figure 7a, are parametrically exported from the Rhino geometry file using Grasshopper plugins, which transformed the drawn plans into OpenFOAM object functions. The average wind speed calculated on each of the plans/rooms is saved as a sub-folder among the post-processing files. A MatLab script accesses these folders and collects the recorded results, grouping them according to their respective simulated incident wind direction.

ACH-Integration Method (All Velocity Vectors)
The multi-platform ParaView was used to fulfill two purposes: qualitative and quantitative analyses of results. While one allows a visual examination, observing airflows/velocity vectors, the other aims to calculate the ACH through the integration method. For this, the airflow rates (m 3 /s) are determined by applying the integrate variables filter in the vertical planes forming the openings crosssection of the investigated rooms within the three-dimensional CFD model. However, in ParaView, all vectors are computed, not just the positive ones, which can lead to calculation errors. Figure 7b shows these vertical planes (solid fill), and the red rectangle represents the slice plane that cuts one of the openings. The airflow is calculated in each of the room openings that, when added together, give the room ventilation rate (Equation (1)), converted to ACH using Equation (2). These surfaces, shown in yellow in Figure 7a, are parametrically exported from the Rhino geometry file using Grasshopper plugins, which transformed the drawn plans into Open-FOAM object functions. The average wind speed calculated on each of the plans/rooms is saved as a sub-folder among the post-processing files. A MatLab script accesses these folders and collects the recorded results, grouping them according to their respective simulated incident wind direction.

ACH-Integration Method (All Velocity Vectors)
The multi-platform ParaView was used to fulfill two purposes: qualitative and quantitative analyses of results. While one allows a visual examination, observing airflows/velocity vectors, the other aims to calculate the ACH through the integration method. For this, the airflow rates (m³/s) are determined by applying the integrate variables filter in the vertical planes forming the openings crosssection of the investigated rooms within the three-dimensional CFD model. However, in ParaView, all vectors are computed, not just the positive ones, which can lead to calculation errors. Figure 7b shows these vertical planes (solid fill), and the red rectangle represents the slice plane that cuts one of the openings. The airflow is calculated in each of the room openings that, when added together, give the room ventilation rate (Equation (1)), converted to ACH using Equation (2).

ACH-Integration Method (Without Recirculation)
The primary goal was to visualize the velocity vectors coming through the openings, allowing for a comprehensive airflow assessment. For this, the following protocol was adopted: first, vector velocities (x, y, and z axes) that cross each point of the vertical openings section are exported as .csv files; second, a Grasshopper/Rhino definition imports both vectors and points, so one can visualize and analyze it. Grasshopper reads the Para-View .csv files and uses predefined geometry information in Rhino to plot the vector data in its respective vertical sections.
Additionally, the 3D parametric modeling tool allows filtering of the velocity vectors, accounting for just the positive velocities and so removing air recirculation. By using native Grasshopper components, the openings normal planes are identified, and only the velocity vectors entering the rooms are computed, discarding the recirculation flows, represented as red arrows in Figure 8.With this approach, the scalar velocity (m/s) in each opening area (m²) is estimated, determining volume flow (m³/s) and, consequently, the ACH in each of the investigated rooms. Recirculation flows are accounted as 0 m/s, and the total opening area is considered.
One limitation of this approach is that the surface opening must be completely free, which prevents the analysis of different window frames/configurations. Due to this limitation, the procedure was only applied to CFD 2 (cooling season scenario), where all openings were modeled without frames.

ACH-Integration Method (Without Recirculation)
The primary goal was to visualize the velocity vectors coming through the openings, allowing for a comprehensive airflow assessment. For this, the following protocol was adopted: first, vector velocities (x, y, and z axes) that cross each point of the vertical openings section are exported as .csv files; second, a Grasshopper/Rhino definition imports both vectors and points, so one can visualize and analyze it. Grasshopper reads the ParaView .csv files and uses predefined geometry information in Rhino to plot the vector data in its respective vertical sections.
Additionally, the 3D parametric modeling tool allows filtering of the velocity vectors, accounting for just the positive velocities and so removing air recirculation. By using native Grasshopper components, the openings normal planes are identified, and only the velocity vectors entering the rooms are computed, discarding the recirculation flows, represented as red arrows in Figure 8.With this approach, the scalar velocity (m/s) in each opening area (m 2 ) is estimated, determining volume flow (m 3 /s) and, consequently, the ACH in each of the investigated rooms. Recirculation flows are accounted as 0 m/s, and the total opening area is considered.
Moreover, both ParaView and Grasshopper working states allow for settings adjusted for one case to be replicated in others. Thus, only one post-processing definition was created to perform qualitative analyses in ParaView and import the .csv files in Grasshopper/Rhino. The most intensive task was generating the vector velocities (.csv files) for the ACH calculations, where the rotation angles needed to be updated for each of the different incident directions simulated at the CFD cases.  Figure 9 illustrates the flow distribution for the most occurring wind direction (θ = 180°) at the I-MA house in a qualitative perspective, modeled after the window settings employed during measurements [78]. The dashed lines represent the measuring position over which the predicted indoor air velocities for each of the investigated wind directions were plotted (Figures 10 and 11), and the rectangle defines the region covered by the anemometer. First, the indoor air velocities (m/s) recorded in the living room were screened so that the experimental campaign measurements could be related to the CFD 1 outputs. In this sense, data recorded with closed windows (7:00 a.m.-8:59 p.m.) were excluded, and the remaining readings were grouped and averaged given the same angle range used in CFD simulations, resulting in a single measured value for each wind direction. Therefore, the One limitation of this approach is that the surface opening must be completely free, which prevents the analysis of different window frames/configurations. Due to this limitation, the procedure was only applied to CFD 2 (cooling season scenario), where all openings were modeled without frames.

CFD 1-Measurement Week
Moreover, both ParaView and Grasshopper working states allow for settings adjusted for one case to be replicated in others. Thus, only one post-processing definition was created to perform qualitative analyses in ParaView and import the .csv files in Grasshopper/Rhino. The most intensive task was generating the vector velocities (.csv files) for the ACH calculations, where the rotation angles needed to be updated for each of the different incident directions simulated at the CFD cases. Figure 9 illustrates the flow distribution for the most occurring wind direction (θ = 180 • ) at the I-MA house in a qualitative perspective, modeled after the window settings employed during measurements [78]. The dashed lines represent the measuring position over which the predicted indoor air velocities for each of the investigated wind directions were plotted (Figures 10 and 11), and the rectangle defines the region covered by the anemometer. Moreover, both ParaView and Grasshopper working states allow for settings adjusted for one case to be replicated in others. Thus, only one post-processing definition was created to perform qualitative analyses in ParaView and import the .csv files in Grasshopper/Rhino. The most intensive task was generating the vector velocities (.csv files) for the ACH calculations, where the rotation angles needed to be updated for each of the different incident directions simulated at the CFD cases.  Figure 9 illustrates the flow distribution for the most occurring wind direction (θ = 180°) at the I-MA house in a qualitative perspective, modeled after the window settings employed during measurements [78]. The dashed lines represent the measuring position over which the predicted indoor air velocities for each of the investigated wind directions were plotted (Figures 10 and 11), and the rectangle defines the region covered by the anemometer. First, the indoor air velocities (m/s) recorded in the living room were screened so that the experimental campaign measurements could be related to the CFD 1 outputs. In this sense, data recorded with closed windows (7:00 a.m.-8:59 p.m.) were excluded, and the remaining readings were grouped and averaged given the same angle range used in CFD simulations, resulting in a single measured value for each wind direction. Therefore, the First, the indoor air velocities (m/s) recorded in the living room were screened so that the experimental campaign measurements could be related to the CFD 1 outputs. In this sense, data recorded with closed windows (7:00 a.m.-8:59 p.m.) were excluded, and the remaining readings were grouped and averaged given the same angle range used in CFD simulations, resulting in a single measured value for each wind direction. Therefore        As the building ventilation system (no-load supply) was in operation throughout the experiment, even in the occurrence of external calm winds, air velocity measurements remained around 0.065 m/s, which does not necessarily indicate natural ventilation but rather an effect of the mechanized system. This is the case of incident winds at θ = 150°, θ = 180°, θ = 210°, and θ = 240° that have wind speed boundary condition assigned as less than 1.5 m, and the mechanized ventilation prevails over natural. As the artificial flow was not modeled, CFD simulations output indoor air velocities much slower than those measured, differing by about 30%. Therefore, these wind directions are not used as a reference comparison to the numerical model results and are not presented in Figure 12.

CFD 1-Measurement Week
Nevertheless, while the air velocities across the spaces were non-uniform, varying significantly along the longitudinal axis, the average values calculated in the horizontal planes are a reasonable representation of indoor airflows. The representation efficiently captures the room ventilation performance, providing a more comprehensive analysis than that coming from a single sensor, restricted to a small accuracy. Hence, calculation planes configure a convenient and reliable method when assessing natural ventilation performance, in both building assessment or even in the project development phase.

CFD 2-Cooling Season
• Air velocities In Figure 13, the air velocities calculated in the horizontal planes (1.2 m above the floor) for the cellar (C), living room (LV), and the three bedrooms (BR 1-3) are shown as bars for all ten incident wind directions investigated, and the inlet wind speed used in each angle is presented as a dotted line. The internal airflow varies significantly due to wind direction fluctuations being more sensitive to the direction than the outdoor wind speed, as demonstrated in [62]. For example, the cellar presents both the highest and lowest predicted indoor air velocities, 0.88 m/s and 0.08 m/s for incident wind directions θ = 0° and θ = 270°, respectively.
Considering all the investigated rooms, air velocities vary from 0.1 to 0.4 m/s, with an average value around 0.34 m/s on the ground floor (cellar (C) and living room (LR)) and 0.23 m/s on the first floor (bedroom 1-3 (BR1-BR3)), which, besides smaller openings, also has a balcony. Within this air velocity range, internal operative temperatures up to 28 °C would be acceptable according to comfort standards. However, other passive or active strategies must be utilized in higher temperatures, such as buoyancy-driven airflows or fans. As the building ventilation system (no-load supply) was in operation throughout the experiment, even in the occurrence of external calm winds, air velocity measurements remained around 0.065 m/s, which does not necessarily indicate natural ventilation but rather an effect of the mechanized system. This is the case of incident winds at θ = 150 • , θ = 180 • , θ = 210 • , and θ = 240 • that have wind speed boundary condition assigned as less than 1.5 m, and the mechanized ventilation prevails over natural. As the artificial flow was not modeled, CFD simulations output indoor air velocities much slower than those measured, differing by about 30%. Therefore, these wind directions are not used as a reference comparison to the numerical model results and are not presented in Figure 12.
Nevertheless, while the air velocities across the spaces were non-uniform, varying significantly along the longitudinal axis, the average values calculated in the horizontal planes are a reasonable representation of indoor airflows. The representation efficiently captures the room ventilation performance, providing a more comprehensive analysis than that coming from a single sensor, restricted to a small accuracy. Hence, calculation planes configure a convenient and reliable method when assessing natural ventilation performance, in both building assessment or even in the project development phase.

•
Air velocities In Figure 13, the air velocities calculated in the horizontal planes (1.2 m above the floor) for the cellar (C), living room (LV), and the three bedrooms (BR 1-3) are shown as bars for all ten incident wind directions investigated, and the inlet wind speed used in each angle is presented as a dotted line. The internal airflow varies significantly due to wind direction fluctuations being more sensitive to the direction than the outdoor wind speed, as demonstrated in [62]. For example, the cellar presents both the highest and lowest predicted indoor air velocities, 0.88 m/s and 0.08 m/s for incident wind directions θ = 0 • and θ = 270 • , respectively.
Considering all the investigated rooms, air velocities vary from 0.1 to 0.4 m/s, with an average value around 0.34 m/s on the ground floor (cellar (C) and living room (LR)) and 0.23 m/s on the first floor (bedroom 1-3 (BR1-BR3)), which, besides smaller openings, also has a balcony. Within this air velocity range, internal operative temperatures up to 28 • C would be acceptable according to comfort standards. However, other passive or active strategies must be utilized in higher temperatures, such as buoyancy-driven airflows or fans. The wind direction's influence on the air flows and the predicted air velocities can be observed in Figure 14. The image combines the distribution of normalized mean wind speed Unorm (the mean wind speed values were rescaled to have mean = 0 and standard deviation = 1) and streamlines patterns around the ground and first floor. For the sake of space, only the three most frequent incident wind directions (θ = 150°, 180°, 210°) and the highest speed (θ = 0°) are portrayed. Both predicted Unorm and streamlines vary significantly with incident wind directions in terms of eddies' existence, size, and position in the rooms. Depending on the incident wind angle, the flows in space may present either a pattern from a cross-ventilated room, with a more uniform distribution (Figure 14c,e, or of a single-sided one, showing eddy formations in the room corners (Figure 14a,g).
Moreover, higher indoor air velocities are observed when the building orientation is oblique to an incident wind angle. The case with the south wind (θ = 180°), for example, perpendicular to the building's openings (Figure 14e The wind direction's influence on the air flows and the predicted air velocities can be observed in Figure 14. The image combines the distribution of normalized mean wind speed U norm (the mean wind speed values were rescaled to have mean = 0 and standard deviation = 1) and streamlines patterns around the ground and first floor. For the sake of space, only the three most frequent incident wind directions (θ = 150 • , 180 • , 210 • ) and the highest speed (θ = 0 • ) are portrayed. Both predicted U norm and streamlines vary significantly with incident wind directions in terms of eddies' existence, size, and position in the rooms. Depending on the incident wind angle, the flows in space may present either a pattern from a cross-ventilated room, with a more uniform distribution (Figure 14c  The wind direction's influence on the air flows and the predicted air velocities can be observed in Figure 14. The image combines the distribution of normalized mean wind speed Unorm (the mean wind speed values were rescaled to have mean = 0 and standard deviation = 1) and streamlines patterns around the ground and first floor. For the sake o space, only the three most frequent incident wind directions (θ = 150°, 180°, 210°) and the highest speed (θ = 0°) are portrayed. Both predicted Unorm and streamlines vary signifi cantly with incident wind directions in terms of eddies' existence, size, and position in the rooms. Depending on the incident wind angle, the flows in space may present either a pattern from a cross-ventilated room, with a more uniform distribution (Figure 14c,e, o of a single-sided one, showing eddy formations in the room corners (Figure 14a,g).
Moreover, higher indoor air velocities are observed when the building orientation i oblique to an incident wind angle. The case with the south wind (θ = 180°), for example perpendicular to the building's openings (Figure 14e,f)), presents lower air velocity value than the oblique one, with the south-southwest direction (θ = 210°) (Figure 14g,h), alt hough both have the same inlet wind speed of 1.55 m/s. While the bedrooms 1 and 2 show air velocities of 0.11 and 0.25 m/s for θ = 180°, when θ = 210°, the values are 0.26 and 0.34 m/s, respectively. This influence of the building orientation concerning the wind inciden angle is also observed in [18].

Ground Floor
First Floor  •

ACH-Integration Method
A comparison between the ACH calculated with the two integration methods described in Sections 4.2 and 4.3 is presented in Figure 15. for the living room, cellar, and bedrooms 1 and 2. Bedroom 3 is omitted for brevity and because it shows a similar airflow behavior as bedroom 2. Moreover, higher indoor air velocities are observed when the building orientation is oblique to an incident wind angle. The case with the south wind (θ = 180 • ), for example, perpendicular to the building's openings (Figure 14e,f)), presents lower air velocity values than the oblique one, with the south-southwest direction (θ = 210 • ) (Figure 14g respectively. This influence of the building orientation concerning the wind incident angle is also observed in [18].

• ACH-Integration Method
A comparison between the ACH calculated with the two integration methods described in Sections 4.2 and 4.3 is presented in Figure 15. for the living room, cellar, and bedrooms 1 and 2. Bedroom 3 is omitted for brevity and because it shows a similar airflow behavior as bedroom 2.
Energies 2021, 14, x FOR PEER REVIEW All velocity vectors Without recirculation Figure 15. Predicted ACH at living room, cellar, and bedrooms 1 and 2-comparison between two approaches of integration method.
A direct relationship can be observed between the predicted velocities calcu the horizontal plane ( Figure 14) and the ACH estimated with the integration me the same time, the ventilation rates restricted to the velocity vectors entering t (dotted line with square markers) are in general lower in the ground floor and b the first floor than those using all velocities that pass through the openings (con lines with cross-markers). As the ACH values depend on the room volume, th room (65.7 m³) presents the lowest air rates, varying from 1.66 to 2.8 (ACH −1 ), u method described in the Section 4.2, against 1.6 to 3.3 (ACH −1 ) when disregarding lation flows (Section 4.3). Despite the different ventilation rates found in each wind direction, filtered ACH values calculated by combining CFD results to Gras components vary, on average, between −11% to +6%, which can impact perform sessments. In Appendix A, Figure A1 shows the visualization in Rhino of the normal velocities coming through the openings for the south wind direction (θ =

Research Constraints and Remarks
This study employs CFD to explore wind-driven natural ventilation through ametric modeling platforms. Accuracy of the assessment process is ensured by c defining the boundary conditions and numerical settings, besides comparing CFD (velocity components and indoor air velocities) with data from a wind tunnel exp and a large-scale measurement campaign. Nevertheless, future studies and desig cations must consider some limitations regarding this investigation, including: A direct relationship can be observed between the predicted velocities calculated in the horizontal plane ( Figure 14) and the ACH estimated with the integration method. At the same time, the ventilation rates restricted to the velocity vectors entering the room (dotted line with square markers) are in general lower in the ground floor and bigger in the first floor than those using all velocities that pass through the openings (continuous lines with cross-markers). As the ACH values depend on the room volume, the living room (65.7 m 3 ) presents the lowest air rates, varying from 1.66 to 2.8 (ACH −1 ), using the method described in the Section 4.2, against 1.6 to 3.3 (ACH −1 ) when disregarding recirculation flows (Section 4.3). Despite the different ventilation rates found in each incident wind direction, filtered ACH values calculated by combining CFD results to Grasshopper components vary, on average, between −11% to +6%, which can impact performance assessments. In Appendix A, Figure A1 shows the visualization in Rhino of the average normal velocities coming through the openings for the south wind direction (θ = 180 • ).

Research Constraints and Remarks
This study employs CFD to explore wind-driven natural ventilation through 3D parametric modeling platforms. Accuracy of the assessment process is ensured by carefully defining the boundary conditions and numerical settings, besides comparing CFD results (velocity components and indoor air velocities) with data from a wind tunnel experi- • Assessment of natural ventilation occurred in wind-driven conditions. Different solver settings and boundary conditions are required to investigate buoyancy-driven ventilation.

•
The CFD domain did not include immediate surroundings. For analyses in an urban context, modeling the neighborhood should be considered.

•
The influence of window opening degree was not investigated, since all openings were modeled without frames opened. This approach delivers the maximum ventilation potential of the building. However, as window operation influences natural ventilation performance [105], window opening schedules could be modeled for a more precise assessment, when known.

•
The current study measured air velocity (m/s) and used RANS and k-ε turbulence models in the investigations. Despite its effectiveness, a different approach would be required to properly verify the ACH calculated by the adapted integration method described in Section 4.3. Some alternatives include using LES or other turbulence models such as shear-stress transport (SST) k-ω associated with physical experimentation involving airflow prediction through the tracer-gas decay method. Testing different turbulence models comprise the scope of future work.

•
Displaying velocity vectors within the 3D parametric design software is an alternative in light of qualitative airflow evaluations. However, further experimental campaigns should be programmed for tracer-gas measurements besides deeper investigations to allow for checking the accuracy of the approach when employed to solve the integration method excluding recirculation flows.
Finally, some advantages and disadvantages of integrating parametric 3D modeling with CFD are highlighted, compared to conventional CFD models: Pros: • The CFD workflow and tasks are organized with visual programming resources rather than code lines, allowing a better understanding of their steps, especially for users unfamiliar with conventional programming languages.

•
The integration between the project's geometry and the CFD model's automatic generation facilitates performing CFD simulations during either project development or building performance assessment. As the tools are embedded in a single platform, parameterized geometric configurations can be easily changed and tested. • Several configurations are pre-established within the CFD plugin in the 3D software and have default values based on natural ventilation studies' best practices. Therefore, designers can perform CFD simulations with more reliability. Nonetheless, it is possible to change or add new functions using all OpenFOAM capabilities. As a result, advanced users can adapt their models as desired and profit from the 3D parametric resources.

Cons:
• Although OpenFOAM is a free, open-source platform, the proposed CFD workflow needs a commercial 3D modeling software. • Compared to other commercial CFD programs, such as Ansys-Fluent and CFX, the number of predetermined features to set up a model with the 3D modeling platform's CFD plugin is smaller. Hence, the user might have to make a more significant effort when trying to model a problem outside the available files and settings.

Conclusions
This paper presented a framework to systematize the analyses of a building's winddriven natural ventilation potential through CFD simulation. Combining available programs and developing customized functions, the study provided a pragmatic use of CFD for either building design or assessment. A full-scale test house was used as a base case, its geometry and the climatic characteristics of its site were settled as a reference, and air velocity measurements performed within the experimental house served as validation data. Insights regarding the procedure that covers from pre-to post-processing steps include:

•
The proposed workflow assists practical use of CFD, combining more accessible winddriven airflow simulation resources and an oriented method, which can promote CFD applications at the design phase and encourage comprehensive building evaluations.

•
The detailed information regarding wind data definition and numerical model settings can help future studies and applications willing to use CFD in natural ventilation assessments.

•
Combining 3D modeling platform utilities with CFD objective functions to generate parameterized calculation surfaces automates post-processing and accelerates the analysis. It provides an effective quantitative way to investigate naturally ventilated buildings with CFD.
Lastly, the CFD simulations provided information on the ventilative effectiveness of the test house, used as a reference building in this study. The results lead to the following conclusions: 1.
Since air velocities across the room are non-uniform, the horizontal calculation plans provide a quick but rather general space reference value. For a more detailed analysis, one must check velocity contours and streamlines.

2.
For the cooling season, 10 of the 12 wind directions considered at the climactic analysis were simulated, preferring those most frequent. Although the number of wind directions to be considered is adjustable, it is noted that wind direction fluctuations considerably influence the internal airflow. Therefore, the greater the number of incident wind directions, the better the building characterization. However, if one must limit cases, prioritizing the oblique winds based on occurrence would be recommended as they are more effective for open windows.
Supplementary Materials: The following are available online at https://data.mendeley.com/datasets/ hy6mjv8f6t/draft?a=54555508-b14e-4c6a-88a7-379e836d394d, 3D model (Rhino file), Grasshopper, and OpenFoam files (θ = 0 • ) of the natural ventilation investigations (scenario 2) performed with the I-MA experimental house. The documents used to post-process ACH, including the ParaView description, Grasshopper, Rhino, and .csv files, are also accessible within the link. Funding: This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.