Determination of a Methodology to Derive Correlations Between Window Opening Mass Flow Rate and Wind Conditions Based on CFD Results

: This paper presents a methodology for the development of an empirical equation which can provide the air mass ﬂow rate imposed by single-sided wind-driven ventilation of a room, as a function of external wind speed and direction, using the results from Computational Fluid Dynamics (CFD) simulations. The proposed methodology is useful for a wide spectrum of applications, in which no access to experimental data or conduction of several CFD runs is possible, deriving a simple expression of natural ventilation rate, which can be further used for energy analysis of complicated building geometries in 0-D models or in object-oriented software codes. The developed computational model simulates a building, which belongs to Rheinisch-Westfälische Technische Hochschule (RWTH, Aachen University, Aachen, Germany) and its surrounding environment. A tilted window represents the opening that allows the ventilation of the adjacent room with fresh air. The derived data from the CFD simulations for the air mass ﬂow were ﬁtted with a Gaussian function in order to achieve the development of an empirical equation. The numerical simulations have been conducted using the Ansys Fluent v15.0 ® software package. In this work, the k-w Shear Stress Transport (SST) model was implemented for the simulation of turbulence, while the Boussinesq approximation was used for the simulation of the buoyancy forces. The coe ﬃ cient of determination R 2 of the curve is in the range of 0.84–0.95, depending on the wind speed. This function can provide the mass ﬂow rate through the open window of the investigated building and subsequently the ventilation rate of the adjacent room in air speed range from 2.5 m / s to 16 m / s without the necessity of further numerical simulations. along the y-axis for 6, 10 and 15 m / s wind speed, respectively. These ﬁgures show both a general view of the building (a) and a magniﬁed view of the area of interest in front of the window (b). The negative values of velocity z-component represent the direction of the ﬂow towards the window.


Introduction
Nowadays, natural ventilation of a building is critical in order to reduce energy consumption for space conditioning (cooling and ventilation). Natural ventilation is the process by which clean air, normally outdoor air, is intentionally provided to a space and stale air is removed, without using mechanical systems [1]. The purpose of natural ventilation is to achieve maximum human comfort in indoor spaces by ensuring maximum use of natural resources [2].
There are two types of natural ventilation in buildings: wind-driven and buoyancy-driven ventilation [3]. Wind-driven ventilation arises from the different pressure created by external wind conditions, while openings in the perimeter of the building permit the flow infiltration [4]. Buoyancy-driven ventilation occurs because of temperature difference between the interior and exterior air. Normally, there is a combination of these two phenomena of ventilation. In order for the designers to improve the energy efficiency of the buildings, they use two mechanisms for improving the natural ventilation: (a) the single-sided ventilation (SS) and (b) the cross-flow ventilation (CR) [5,6]. In both 2 of 21 cases, ventilation is either wind or buoyancy-driven or both and occurs through the openings of the building. Cross-flow ventilation (CR) is achieved using windows on both sides of a building. SS ventilation describes a situation in which wind enters and leaves the building through one or two openings located on the same side of the building or the room. Ventilation within the building is mainly affected by the geophysical morphology and the surrounding buildings [7], the used building side for the ventilation mechanism, the type of the openings and the external wind conditions [8]. Therefore, the numerical study of the natural ventilation process is complicated due to the fact that the airflow is affected from multiple factors simultaneously.
As understood, SS ventilation uses only one side of the building, so it is less efficient compared to CR ventilation. Although single-sided ventilation is less efficient, it is more suitable for cellular room environments, such as offices, because they have not openings on both wall sides and cannot implement the cross-flow ventilation mechanism [9].
In general, there are two main methods to define the mass flow rate and the flow parameters through a building's opening; the first one is the calculation by experiments either in full-scale buildings or in wind tunnel [10] and the second one is the numerical prediction [11]. The numerical approach gives the flexibility to study several cases of building structures and environmental conditions. In this paper, our efforts are focused on the development of a function that can predict the ventilation air mass flow rate through a building's opening as a function of external wind direction and speed.
Several investigations have been carried out in order to develop empirical equations that are capable of predicting the ventilation rate of buildings. In 1980, Phaff and De Gids et al. [12] proposed an empirical expression to calculate the airflow rate in a single-sided ventilation zone, based on opening area, wind speed and air temperature. The empirical expression is based on measurements that have been performed on the first floor of buildings, which are surrounded by other buildings up to four floors high.
Warren and Parkins [13,14] also proposed an empirical expression for buoyancy-driven and one for wind-driven single-sided natural ventilation. These expressions are function of the opening dimensions, such as window's area and height, the gravitational acceleration, the average temperature difference between indoor spaces and outdoor environment and wind speed. A combination of these two expressions by quadrature function yields the final equation for the calculation of ventilation flow rate. In 2005, the American Society of Heat, Refrigerating and Air-Conditioning Engineers [15] proposed the first mathematical expression to calculate airflow in single-sided ventilation that takes into account the wind angle of incidence.
Many investigations used the empirical expressions that have already been mentioned for validation purposes. In this framework, Alloca et al. [16] investigated the single-sided natural ventilation of a building by using a CFD model together with an analytical/empirical model. The analytical method involves the equations from Bernoulli theory for the buoyancy-driven flow with an empirical discharge coefficient using the empirical model of Phaff and de Grids [12], for single-sided ventilation. This computational study investigates two cases. The first case is a typical student dormitory with two openings; an upper and a lower window on the same wall side. The second case is a three-story apartment composed of three identical dormitory rooms stacked vertically above one another. The CFD model follows the same trends as the empirical model, but underestimates the ventilation rates by 35%.
Asfour et al. [17] presented a comparison between the airflow rate calculated with a network mathematical model and the airflow rate calculated with CFD simulations. In this work, three different rooms with nearly the same volume, but different aspect ratios, were studied. Each case examines two wind speeds, 1 and 5 m/s, and two wind angles, namely 0 • and 45 • . The discrepancy percentage between estimated and calculated airflow rate is in the range of −11.5% and 5.3%. Due to the good agreement between the results of the two models, the network mathematical model can be used as a validation tool of CFD studies that have no access to experimental data. Caciolo et al. [18] examined the accuracy of the air change rate predictions by the already existing empirical correlations of Warren, Phaff and De Gibs [12], Larsen and Heiselberg [19] and Dascalakis [20]. They conclude that in case of leeward opening, all correlations overestimate the airflow rate. This is attributed to the fact that these correlations do not consider the reduction of the stack effect due to the existence of turbulent diffusion at the opening. On the contrary, in the case of windward opening, the Warren's correlation shows the best agreement with experiments. In a later investigation, Caciolo et al. [21] presented a new correlation for the case of the leeward side by using CFD simulation results. The new correlation shows a better agreement with experiment results, compared to existing correlations.
Tang et al. [22] proposed a new hypothetical correlation based on an experimental study for prediction of airflow rates in a primary school in Beijing in the case of low air speed values and insignificant buoyancy effects, which implement the development of unorganized airflow. Conducting 168 h of experiments, they compare the measured airflow rates against the values derived from the already existing correlations [12,14,19,20]. In a second step, they propose a new hypothetical correlation capable of predicting more accurately the airflow rates in the case of unorganized airflow. The new correlation shows a good prediction of airflow rate. The average deviation is reduced to 17.37%, which is 7% less than the lowest deviation attained from existing correlations.
Wang et al. [23] presented an empirical model capable of predicting the mass flow rate induced by single-sided wind-driven ventilation due to the pressure difference along an opening height. For validation purposes, CFD simulations are performed. The maximum difference between the empirical model prediction and the CFD results is less than 25%. The largest difference is found in the case of leeward side, in which the flow field near the opening is much more complicated. In a later investigation, Wang [24] studied the impact of three types of window, i.e., hopper, awning and casement, in the case of SS natural ventilation. The expression of airflow rate is a function of window and building geometry, opening angle, wind incident angle and speed. These semi-empirical models are based on the previously analytical model and on pressure coefficients. The validation of these expressions has been achieved by using experimental measurements with the tracer-gas method combined with CFD numerical simulations. The new semi-empirical model for predicting the aeration rate for the three types of window presents a good agreement between the measured data and the CFD results.
Pan et al. [25] presented a model for calculating the ventilation rate in SS natural ventilation of an apartment due to wind-and buoyancy-driven effects, based on their previous empirical model [23]. The model is validated by using measured data and is able to predict natural ventilation rate with an average error of 12.7%. The air temperature difference between the indoor and outdoor space is ranged from −2.3 K to 13.2 K. Compared to other six empirical correlations [12,13,[18][19][20]22], this model performs well due to the fact that the other models calculate the ventilation rates with average errors ranging from 12.9% to 46.1%. Moreover, this model takes into account the impact of both positive and negative buoyancy forces along with outside air pressure on natural ventilation through a single opening in contrast with the other models available.
As it was expected, the preceding literature survey shows that none of the models available in the literature is ideal. Although the existing models are sophisticated and have functional forms, the majority of them requires the knowledge of pressure coefficients, discharge coefficients or correction factors for each type of opening. These coefficients are usually obtained experimentally or from standard pressure coefficients. Therefore, it seems that there is an obvious gap in the literature regarding the existence of a simple, but versatile and credible, methodology for the derivation of an expression, which can provide the aeration rate through a building's opening with no dependence on sophisticated experimental or numerical data. Furthermore, this work has examined a wider range of wind speeds compared to existing studies in order to derive the mathematical function of the airflow rate. The selected wind speeds are equal to 6, 10 and 15 m/s, which correspond to 2, 5 and 7 bft wind speeds on the Beaufort scale, respectively. Furthermore, there is no symmetry in the model, in contrast to other studies, since this is a real building and not a theoretical one. Finally, for each case of wind speed, five different wind directions have been applied.
The implemented model has already been developed and validated against experimental data in a well-controlled environment inside a wind tunnel by some members of the present research group in past studies [26][27][28]. However, the direct comparison of the CFD results against the experimental data regarding the aeration rate of the examined room is not possible, since the wind speed and direction that are used as time-averaged input values in the CFD model are measured by a weather station that is located on the top of the building inside the developed boundary layer. Therefore, due to recirculations inside the boundary layer, these values cannot be considered as representative to time-averaged values of the free airflow conditions and the exerted aeration rates of the investigated room cannot be supposed that have been resulted by these wind conditions. More information regarding the problem with the location of the weather station and the experimental data will be further provided in the Results section. Moreover, an extra effort of validity has been made using previous empirical correlations. However, this effort was not successful, because the existing empirical correlations include pressure, discharge and correction coefficients. These parameters are case-dependent, significantly affected by the opening type and usually defined by experiments. In this specific work, there is no available experimental data for these specific coefficients, so the implementation of the respective equations is not possible either. Additionally, the semi-empirical model of Wang includes the term of neutral plane (plane with zero air velocity), which cannot be defined in this work, due to the type and the opening of the window. Therefore, it was also impossible to make a direct comparison of the results provided by this work with this specific model.
In order to achieve accurate simulation of air-flow, a three-dimensional model has been chosen, based on a building that is located at RWTH Aachen University. An appropriate expression for aeration rate as a function of wind speed and angle of incidence is fitted to the normalized data. In this study, the best fit was achieved using a type of Gaussian function. In order to ensure the verification of the derived equation, three additional wind speeds have been selected; one inside the range of 6 to 15 m/s and two outside the limits of this range, to compare the calculated airflow rates against the estimated ones by the mathematical function. The agreement is good since the maximum relative difference is below 10%, except the cases with wind flow parallel to the building, where the maximum relative difference can be as high as 38%. This high relative difference is attributed to recirculations that are developed in front of the window opening and the empirical correlation cannot take into account. More information regarding this issue can be found in the Results section. To sum up, this methodology is useful for a wide spectrum of applications, in which no access to experimental data or conduction of several CFD runs is possible. Moreover, with this methodology a simple expression of natural ventilation rate can be exported. This expression can be used for further energy analysis of complicated building geometries in 0-D models or in object-oriented software codes. Finally, even if the obtained correlation is not general and can only be used for this specific window type and this specific building envelope, the methodology is generic and can be followed in all cases.

Mathematical Model
This work simulates: (a) the developed flow field affected by the wind conditions and the natural convection, (b) turbulence and (c) the energy transfer due to convection and diffusion. The natural convection mechanism is attributed to the implemented temperature difference between ambient air and the room wall temperature. In this work, the applied buoyancy forces are calculated by using the Boussinesq approximation. In combined radiation and convection heat transfer problems, the Boltzmann number represents the ratio of convection to radiation heat transfer [29]. In this work, this specific ratio is very high, i.e., equal to Bo = 334, and thus radiation effects can be neglected. All numerical simulations are solved in steady-state conditions, assuming that the implemented free airflow conditions represent time-averaged values. Since the free airflow conditions (i.e., wind speed and wind direction) are considered to be steady and representative to time-averaged values, it is expected that the transient analysis would eventually have provided the same results with the steady-state analysis regarding the mass flow rate through the window opening, when the calculation time has sufficiently flowed. Therefore, in order to save significant calculation time and avoid convergence issues that might arise in the transient calculation, the steady-state approximation has been followed. Furthermore, since all CFD simulations were steady-state cases and the mesh in tilted window area has high skewness value, the pressure-velocity coupling is achieved by using the Semi-Implicit Method (SIMPLE) algorithm [30,31]. The momentum, energy, turbulent kinetic energy and specific dissipation rate for the first 800 iterations are spatially discretized by using the first-order upwind scheme. After 800 iterations, second-order accuracy schemes are used for momentum and energy equations. The convergence of the steady-state simulation is controlled by monitoring the mass flow rate through the opening using a User-Defined Function (UDF). Figure 1a presents the convergence of the calculated air mass flow rate that gets into the room during the simulation process, assuming 6m/s wind speed and all five cases for the angle of incidence. The simulation of each case is considered converged when the mass flow rate through the opening tend to oscillate around a constant value, see Figure 1b. The residuals of continuity, velocity and energy are below 10 −5 , 10 −5 and 10 −8 , respectively. The mass flow rate for each case is defined by its mean value during the last 500 iterations. state analysis regarding the mass flow rate through the window opening, when the calculation time has sufficiently flowed. Therefore, in order to save significant calculation time and avoid convergence issues that might arise in the transient calculation, the steady-state approximation has been followed. Furthermore, since all CFD simulations were steady-state cases and the mesh in tilted window area has high skewness value, the pressure-velocity coupling is achieved by using the Semi-Implicit Method (SIMPLE) algorithm [30,31]. The momentum, energy, turbulent kinetic energy and specific dissipation rate for the first 800 iterations are spatially discretized by using the first-order upwind scheme. After 800 iterations, second-order accuracy schemes are used for momentum and energy equations. The convergence of the steady-state simulation is controlled by monitoring the mass flow rate through the opening using a User-Defined Function (UDF). Figure 1a presents the convergence of the calculated air mass flow rate that gets into the room during the simulation process, assuming 6m/s wind speed and all five cases for the angle of incidence. The simulation of each case is considered converged when the mass flow rate through the opening tend to oscillate around a constant value, see Figure 1b. The residuals of continuity, velocity and energy are below 10 -5 , 10 -5 and 10 -8 , respectively. The mass flow rate for each case is defined by its mean value during the last 500 iterations.
where p is static pressure, τ is the stress tensor and ρg ⃗ is the gravitational body force. In this study, the term of gravitational body force, ρg ⃗, contains the Boussinesq approximation (Equation (6)). In Equation (3), T denotes the temperature and keff is the effective conductivity.

Turbulence Model
The advantage of the k-w model over the k-epsilon model [32] is the improved performance for the approximation of the boundary layers under adverse pressure gradients and the more accurate predictions regarding: (a) internal flows, (b) flows that exhibit strong curvature, (c) separated flows, and d) jets. However, k-w model also presents a major disadvantage. More specifically, the boundary layer computations are very sensitive to the values of ω in the free stream. In order to overcome this restriction, it is necessary to use the Shear-Stress-Transport (SST) model. Furthermore, Hooff et al. [33] and Ramponi and Blocken reported in [34] that Shear-Stress-Transport (SST) k-ω model provides high accuracy in predicting wind profiles. Thus, the turbulence of the flow was modelled using in all where p is static pressure, = τ is the stress tensor and ρ → g is the gravitational body force. In this study, the term of gravitational body force, ρ → g , contains the Boussinesq approximation (Equation (6)). In Equation (3), T denotes the temperature and k eff is the effective conductivity.

Turbulence Model
The advantage of the k-w model over the k-epsilon model [32] is the improved performance for the approximation of the boundary layers under adverse pressure gradients and the more accurate predictions regarding: (a) internal flows, (b) flows that exhibit strong curvature, (c) separated flows, and d) jets. However, k-w model also presents a major disadvantage. More specifically, the boundary layer computations are very sensitive to the values of ω in the free stream. In order to overcome this Energies 2019, 12, 1600 6 of 21 restriction, it is necessary to use the Shear-Stress-Transport (SST) model. Furthermore, Hooff et al. [33] and Ramponi and Blocken reported in [34] that Shear-Stress-Transport (SST) k-ω model provides high accuracy in predicting wind profiles. Thus, the turbulence of the flow was modelled using in all cases the two-equation Shear-Stress-Transport (SST) k-ω based model, developed by Menter [35]. The turbulence kinetic energy, k, and the specific dissipation rate, ω, are obtained from Equations (4) and (5), respectively:

Boussinesq Model
In natural-convection flow problems, Boussinesq approximation provides a faster convergence of the solution procedure compared to the default case where fluid density is considered as a function of temperature. Therefore, this model has also been implemented in this specific case. This model treats density as a constant value in all solved equations, except for the buoyancy term in the momentum equation, Equation (6), which provides the volume forces due to buoyance: where ρ ref is the constant density of the flow equal to 1.16 kg/m 3 , T ref represents the buoyancy reference temperature, i.e., 286.88 K, and β = 1/T ref is the thermal expansion coefficient equal to 0.00343 K −1 . The Boussinesq approximation is valid when the product β(T − T ref ) is lower than unity [36]. Thus, as the room temperature T is equal to 292.38 K, the condition is fulfilled and the Boussinesq approximation can be used.

Normalization
In order to derive the necessary correlations, it is also necessary to define the appropriate dimensionless parameters. The first parameter refers to the normalized mass flow rate and the second one to the dimensionless direction (angle relative to North-to-South direction). The mathematical formulas for the definition of these dimensionless quantities are given by Equations (7) and (8), where . m max is the maximum mass flow rate (kg/s) among the simulations belonging to the group of the same wind speed and different wind directions (as already explained), . m is the actual mass flow rate (kg/s) numerically calculated at the window opening for each specific case and θ is the angle (o) defining the free airflow direction relative to North-to-South one, resulting in the normalized angles of −0.5, −0.25, 0, 0.25 and 0.45 for 45 • , 90 • , 135 • , 180 • and 215 • wind direction, respectively.

Geometric Model
The geometric model for the computational simulations has been developed in ANSYS DesignModeler ® . The dimensions of the building are 17.98 m × 70.73 m × 11.64 m (L × W × H). According to the existing best practice guidelines of Franke et al. [37], the domain in the flow direction must be extruded by at least eight times the height of the building, when the flow profiles are not available and the flow is blocked to a large extent (e.g., 10%). In this work, the domain has been extruded by approximately 15.5 times the height of the building (H) along z-axis, since the building blockage is quite high (16%). Furthermore, based on the best practices proposed for the case of a single building, the distance between the top of the computational domain and the roof of the building must be 4-10 times the height of the building. In this work, the extrusion of the domain was 2.72 times the height of the building. Even if the distance is lower than the proposed range, it has been ensured, based on the results, that this distance does not insert an artificial acceleration of the flow over the building. Finally, the distance between the building's sidewalls and the lateral boundaries of the computational domain has been selected to be equal to approximately three times the height of the building envelope, since the blockage is quite high (16%). Therefore, the total length of the domain is 197.87 m (z-axis), the total width 142.69 m (x-axis) and the total height 31.64 m (y-axis). The extrusion of the domain far from the building envelop is necessary to simulate the fully developed flow field and to ensure that this is only dependent on the imposed boundary conditions (BCs) and not affected by the building. The extrusion of the domain far from the building envelop is necessary to simulate the fully developed flow field and to ensure that this is only dependent on the imposed boundary conditions and not affected by the building.  2.72 times the height of the building. Even if the distance is lower than the proposed range, it has been ensured, based on the results, that this distance does not insert an artificial acceleration of the flow over the building. Finally, the distance between the building's sidewalls and the lateral boundaries of the computational domain has been selected to be equal to approximately three times the height of the building envelope, since the blockage is quite high (16%). Therefore, the total length of the domain is 197.87 m (z-axis), the total width 142.69 m (x-axis) and the total height 31.64 m (yaxis). The extrusion of the domain far from the building envelop is necessary to simulate the fully developed flow field and to ensure that this is only dependent on the imposed boundary conditions (BCs) and not affected by the building. The extrusion of the domain far from the building envelop is necessary to simulate the fully developed flow field and to ensure that this is only dependent on the imposed boundary conditions and not affected by the building. Figure 2 and Figure    In addition, the numerical domain includes the investigated room with the specific window opening. The dimensions of the room and the location of the window are shown in Figure 4a. The window itself includes a solid boundary that represents the glass, while its frame is included in the wall boundaries of the building envelop. Furthermore, the opening angle is equal to 5.8 o and the effective flow area is equal to 0.485 m 2 (red region in Figure 4a,b). Figure 4b presents the replication of the window opening. The black color represents the window glass, the grey the building walls and the red the window opening. Figure 4c presents the dimensions of the window opening. This type of window was selected as the actual type of the window opening under investigation. In addition, sloping windows are usually found in European homes allowing the efficient ventilation of the building regardless of the weather conditions. The methodology presented in this paper is replicable though, since it can be followed also for other types of windows: awning windows, horizontal and vertical pivot windows or turn windows. 2.72 times the height of the building. Even if the distance is lower than the proposed range, it has been ensured, based on the results, that this distance does not insert an artificial acceleration of the flow over the building. Finally, the distance between the building's sidewalls and the lateral boundaries of the computational domain has been selected to be equal to approximately three times the height of the building envelope, since the blockage is quite high (16%). Therefore, the total length of the domain is 197.87 m (z-axis), the total width 142.69 m (x-axis) and the total height 31.64 m (yaxis). The extrusion of the domain far from the building envelop is necessary to simulate the fully developed flow field and to ensure that this is only dependent on the imposed boundary conditions (BCs) and not affected by the building. The extrusion of the domain far from the building envelop is necessary to simulate the fully developed flow field and to ensure that this is only dependent on the imposed boundary conditions and not affected by the building.   In addition, the numerical domain includes the investigated room with the specific window opening. The dimensions of the room and the location of the window are shown in Figure 4a. The window itself includes a solid boundary that represents the glass, while its frame is included in the wall boundaries of the building envelop. Furthermore, the opening angle is equal to 5.8 o and the effective flow area is equal to 0.485 m 2 (red region in Figure 4a,b). Figure 4b presents the replication of the window opening. The black color represents the window glass, the grey the building walls and the red the window opening. Figure 4c presents the dimensions of the window opening. This type of window was selected as the actual type of the window opening under investigation. In addition, sloping windows are usually found in European homes allowing the efficient ventilation of the building regardless of the weather conditions. The methodology presented in this paper is replicable though, since it can be followed also for other types of windows: awning windows, horizontal and vertical pivot windows or turn windows. In addition, the numerical domain includes the investigated room with the specific window opening. The dimensions of the room and the location of the window are shown in Figure 4a. The window itself includes a solid boundary that represents the glass, while its frame is included in the wall boundaries of the building envelop. Furthermore, the opening angle is equal to 5.8 • and the effective flow area is equal to 0.485 m 2 (red region in Figure 4a,b). Figure 4b presents the replication of the window opening. The black color represents the window glass, the grey the building walls and the red the window opening. Figure 4c presents the dimensions of the window opening. This type of window was selected as the actual type of the window opening under investigation. In addition, sloping windows are usually found in European homes allowing the efficient ventilation of the building regardless of the weather conditions. The methodology presented in this paper is replicable though, since it can be followed also for other types of windows: awning windows, horizontal and vertical pivot windows or turn windows.

Boundary Conditions
The bottom surface is defined as a stationary wall with no-slip condition. On the upper surface a symmetry boundary condition is applied. Moreover, a symmetry boundary assumes that the normal velocity and the normal gradients of all variables are equal to zero. The bottom surface (ground) and the interior and exterior walls of the building are modeled as adiabatic walls with noslip condition. Ambient temperature is defined as equal to 286.88 K and the room walls have a constant temperature equal to 292.38 K. The number and the location of the inlet velocity surfaces of the computational domain are dependent on the wind direction. In general, for every wind speed studied, five different incidence angles are considered, i.e., 45 o , 90 o , 135 o , 180 o and 215 o ( Figure 5). In this study, the northwestern wind directions are not examined since transient phenomena of recirculations appear in front of the window area. The angle of incidence of the air flow is relative to the angle defined by North-to-South direction. Therefore, the north wind presents 0 o angle of incidence, while the south wind 180 o . The outlet surfaces of the computational domain are considered as pressure outlet. The rest surfaces in each case have symmetry boundary conditions. The implemented BCs regarding these crucial operating parameters are given in detail in Table 1.

Boundary Conditions
The bottom surface is defined as a stationary wall with no-slip condition. On the upper surface a symmetry boundary condition is applied. Moreover, a symmetry boundary assumes that the normal velocity and the normal gradients of all variables are equal to zero. The bottom surface (ground) and the interior and exterior walls of the building are modeled as adiabatic walls with no-slip condition. Ambient temperature is defined as equal to 286.88 K and the room walls have a constant temperature equal to 292.38 K. The number and the location of the inlet velocity surfaces of the computational domain are dependent on the wind direction. In general, for every wind speed studied, five different incidence angles are considered, i.e., 45 • , 90 • , 135 • , 180 • and 215 • ( Figure 5). In this study, the northwestern wind directions are not examined since transient phenomena of recirculations appear in front of the window area. The angle of incidence of the air flow is relative to the angle defined by North-to-South direction. Therefore, the north wind presents 0 • angle of incidence, while the south wind 180 • . The outlet surfaces of the computational domain are considered as pressure outlet. The rest surfaces in each case have symmetry boundary conditions. The implemented BCs regarding these crucial operating parameters are given in detail in Table 1.

Boundary Conditions
The bottom surface is defined as a stationary wall with no-slip condition. On the upper surface a symmetry boundary condition is applied. Moreover, a symmetry boundary assumes that the normal velocity and the normal gradients of all variables are equal to zero. The bottom surface (ground) and the interior and exterior walls of the building are modeled as adiabatic walls with noslip condition. Ambient temperature is defined as equal to 286.88 K and the room walls have a constant temperature equal to 292.38 K. The number and the location of the inlet velocity surfaces of the computational domain are dependent on the wind direction. In general, for every wind speed studied, five different incidence angles are considered, i.e., 45 o , 90 o , 135 o , 180 o and 215 o ( Figure 5). In this study, the northwestern wind directions are not examined since transient phenomena of recirculations appear in front of the window area. The angle of incidence of the air flow is relative to the angle defined by North-to-South direction. Therefore, the north wind presents 0 o angle of incidence, while the south wind 180 o . The outlet surfaces of the computational domain are considered as pressure outlet. The rest surfaces in each case have symmetry boundary conditions. The implemented BCs regarding these crucial operating parameters are given in detail in Table 1.    According to Aachen meteorological data, the average wind speed over the course of the year varies in a range of 2 to 16 m/s [38,39]. The mathematical correlation has been derived by using the CFD results in the cases of three different wind speeds, i.e., 6, 10 and 15 m/s. These three values have been commonly agreed with RWTH Aachen University, since for these specific cases the experimental values of aeration rates were available. Nevertheless, as explained in the results section below, the experimental data cannot be used for validation of the CFD model. As a further step, it was necessary to test how accurately the mathematical correlation can predict the aeration rate for any additional wind speed. Therefore, three more wind speeds were used to check the agreement of the CFD results with the results of the empirical expression as regards the mass flow rate through the window opening. For this purpose, three additional wind speeds have been selected; the first is close to the lower bound of the wind range that is typical for Aachen Region, the second one is an intermediate value and the third one is the upper bound of the provided wind range. The selection of the wind directions is arbitrary, except 215 • , which is selected due to the fact that for this wind direction the RWTH Aachen University has some experimental values of aeration rates. Nevertheless, as already mentioned, the experimental data cannot be used for validation of the CFD model.

Numerical Grid
The numerical grid is developed in ANSYS Meshing ® . The mesh consists of approximately 8.7 million cells, all of which are hexahedrons. This type of mesh elements can provide smooth solution convergence and validity of the derived results, as compared to the experimental values or the real operating conditions. Mesh shows high quality, since the skewness factor does not exceed in any case the upper acceptable limit of 0.94. Figure 6 presents the developed numerical grid. More specifically, Figure 6a is a general view of the whole domain, while Figure 6b is an enlarged view of the room, showing the numerical grid among the external environment, the room and the tilted window. According to Aachen meteorological data, the average wind speed over the course of the year varies in a range of 2 to 16 m/s [38,39]. The mathematical correlation has been derived by using the CFD results in the cases of three different wind speeds, i.e., 6, 10 and 15 m/s. These three values have been commonly agreed with RWTH Aachen University, since for these specific cases the experimental values of aeration rates were available. Nevertheless, as explained in the results section below, the experimental data cannot be used for validation of the CFD model. As a further step, it was necessary to test how accurately the mathematical correlation can predict the aeration rate for any additional wind speed. Therefore, three more wind speeds were used to check the agreement of the CFD results with the results of the empirical expression as regards the mass flow rate through the window opening. For this purpose, three additional wind speeds have been selected; the first is close to the lower bound of the wind range that is typical for Aachen Region, the second one is an intermediate value and the third one is the upper bound of the provided wind range. The selection of the wind directions is arbitrary, except 215 o , which is selected due to the fact that for this wind direction the RWTH Aachen University has some experimental values of aeration rates. Nevertheless, as already mentioned, the experimental data cannot be used for validation of the CFD model.

Numerical Grid
The numerical grid is developed in ANSYS Meshing ® . The mesh consists of approximately 8.7 million cells, all of which are hexahedrons. This type of mesh elements can provide smooth solution convergence and validity of the derived results, as compared to the experimental values or the real operating conditions. Mesh shows high quality, since the skewness factor does not exceed in any case the upper acceptable limit of 0.94. Figure 6 presents the developed numerical grid. More specifically, Figure 6(a) is a general view of the whole domain, while Figure 6(b) is an enlarged view of the room, showing the numerical grid among the external environment, the room and the tilted window.

Correlation
In the current study is presented the optimal number of simulations in order to develop the mass flow rate function, since after an assiduous investigation there is no significant effect on the type or coefficients of the fitted functions. Figure 7(a) presents the CFD results' absolute values of real mass flow rate and wind direction for each case, while Figure 7(b) presents the respective normalized values. Due to the above explained, mass flow rate normalization formula used for Figure 7(b), the mass flow rate normalized values cannot be used to compare the flow in a specific angle among the various wind speeds. For instance, at the angle of 45 o (-0.5 of normalized angle), the mass flow rate corresponding to the velocity of 6 m/s is greater than the other two cases, while in the CFD results of Figure 7(a), it has the lowest value of the three. In a general view, it can be seen that the mass flow rate follows the same trend for the range between 90 o and 180 o wind direction. More specifically, the maximum mass flow rate in each case of wind direction is observed in the case with the maximum wind speed, while the minimum one in the case with the minimum wind speed. This is normal, since

Correlation
In the current study is presented the optimal number of simulations in order to develop the mass flow rate function, since after an assiduous investigation there is no significant effect on the type or coefficients of the fitted functions. Figure 7a presents the CFD results' absolute values of real mass flow rate and wind direction for each case, while Figure 7b presents the respective normalized values. Due to the above explained, mass flow rate normalization formula used for Figure 7b, the mass flow rate normalized values cannot be used to compare the flow in a specific angle among the various wind speeds. For instance, at the angle of 45 • (−0.5 of normalized angle), the mass flow rate corresponding to the velocity of 6 m/s is greater than the other two cases, while in the CFD results of Figure 7a, it has the lowest value of the three. In a general view, it can be seen that the mass flow rate follows the same trend for the range between 90 • and 180 • wind direction. More specifically, the maximum mass flow rate in each case of wind direction is observed in the case with the maximum wind speed, while the minimum one in the case with the minimum wind speed. This is normal, since the mass flow rate is proportional to velocity. However, the case that correspondingly represents the lower bound of the wind direction's range does not follow the same trend. More specifically, in the case with the wind direction of 45 • , it can be observed that the mass flow rate in the case with 15 m/s wind speed is located between the respective mass flow values in the cases with 10 m/s and 6 m/s wind speed. In addition, in wind direction of 215 • , the cases with 10 m/s and 15 m/s wind speed have almost the same mass flow rate, in contrast to the differences that are detected between these two cases in other wind directions. Finally, the maximum mass flow rate is presented in the case with wind direction of 135 • , since in this case the air/mass flow is perpendicular to the window opening, so the normal component of the velocity vector takes the maximum value. The wind conditions were measured by a weather station that is located on the building's roof. However, in Figure 8, case of 135 o with 10 m/s wind speed, it can be observed (view of the building from southwestern) that the location is inside the developed boundary layer, so the measurements cannot be considered as representative for time-averaged wind conditions and the measurements of aeration rate as accurate enough for validation purposes.  Based on the general view of flow patterns around the building, a large recirculation area (kidney vortex) at the leeward side of the building can be observed in all cases examined. As the velocity is increased, the vortices become more intense. Inside the room a clockwise vortex is developed as a result of the closed door. More specifically, due to the closed door, mass flow rates The wind conditions were measured by a weather station that is located on the building's roof. However, in Figure 8, case of 135 • with 10 m/s wind speed, it can be observed (view of the building from southwestern) that the location is inside the developed boundary layer, so the measurements cannot be considered as representative for time-averaged wind conditions and the measurements of aeration rate as accurate enough for validation purposes. The wind conditions were measured by a weather station that is located on the building's roof. However, in Figure 8, case of 135 o with 10 m/s wind speed, it can be observed (view of the building from southwestern) that the location is inside the developed boundary layer, so the measurements cannot be considered as representative for time-averaged wind conditions and the measurements of aeration rate as accurate enough for validation purposes.  Based on the general view of flow patterns around the building, a large recirculation area (kidney vortex) at the leeward side of the building can be observed in all cases examined. As the velocity is increased, the vortices become more intense. Inside the room a clockwise vortex is developed as a result of the closed door. More specifically, due to the closed door, mass flow rates   Based on figures 9-11, the flow in all three investigated cases is parallel to the building envelope, due to the wind direction. In the case of the lowest wind speed, i.e., 6 m/s, the air flow is decelerated Based on figures 9-11, the flow in all three investigated cases is parallel to the building envelope, due to the wind direction. In the case of the lowest wind speed, i.e., 6 m/s, the air flow is decelerated Based on figures 9-11, the flow in all three investigated cases is parallel to the building envelope, due to the wind direction. In the case of the lowest wind speed, i.e., 6 m/s, the air flow is decelerated Based on the general view of flow patterns around the building, a large recirculation area (kidney vortex) at the leeward side of the building can be observed in all cases examined. As the velocity is increased, the vortices become more intense. Inside the room a clockwise vortex is developed as a result of the closed door. More specifically, due to the closed door, mass flow rates through the window opening are equal. This mutual and simultaneous mass exchange between the external and internal environment implies the development of a recirculation inside the room.
Based on Figures 9-11, the flow in all three investigated cases is parallel to the building envelope, due to the wind direction. In the case of the lowest wind speed, i.e., 6 m/s, the air flow is decelerated by the building's friction forces and the z-component of wind speed is very low; due to wind direction and the low wind speed. These are the main reasons why the case of 6 m/s wind speed presents the minimum mass flow rate among the three cases. By comparing the rest two cases (10 and 15 m/s wind speed), it can be concluded that the air flow with the maximum velocity is significantly deflected from the building envelope, since the thickness and the required length of fully-developed boundary layer increases, as the velocity increases too. Therefore, the values of the z-component of the velocity (towards the building's window) in the area of interest are higher in the case with 10 m/s wind speed compared to the one with 15 m/s. This results in higher mass flow rate in the case with 10 m/s wind speed compared to the case with 15 m/s. Figure 12 presents the velocity magnitude and the direction of the airflow through the open area of the tilted window for the case of wind direction of 45 • and all the wind speeds investigated. In all three cases, the main mass of air enters the room from the right and the top side of the opening, while it exits through the left one. Furthermore, the two cases with 10 and 15 m/s wind speed present almost identical vectors. However, a difference between these cases can be observed on the upper-left corner, where the flow from the room to the environment in the case of 10 m/s wind speed presents higher velocity compared to the case of 15 m/s, as a result of the increased mass flow rate and the moving of the natural plane of the flow towards the window top side.      As expected, a recirculation zone is developed close to the windward side of the building. This recirculation structure is transferred towards the window, as the wind speed increases. As already explained, when the wind speed increases, so does the thickness and the necessary length of fullydeveloped boundary layer. The development of a recirculation zone exactly in front of the window induces uncertainties in the developed flow field and significantly affects the induced mass flow rate through the respective opening. This is the main reason why the difference in the mass flow rate between the cases of 10 and 15 m/s wind speed is not so significant as in the other direction cases. As expected, a recirculation zone is developed close to the windward side of the building. This recirculation structure is transferred towards the window, as the wind speed increases. As already explained, when the wind speed increases, so does the thickness and the necessary length of fullydeveloped boundary layer. The development of a recirculation zone exactly in front of the window induces uncertainties in the developed flow field and significantly affects the induced mass flow rate through the respective opening. This is the main reason why the difference in the mass flow rate between the cases of 10 and 15 m/s wind speed is not so significant as in the other direction cases. As expected, a recirculation zone is developed close to the windward side of the building. This recirculation structure is transferred towards the window, as the wind speed increases. As already explained, when the wind speed increases, so does the thickness and the necessary length of fully-developed boundary layer. The development of a recirculation zone exactly in front of the window induces uncertainties in the developed flow field and significantly affects the induced mass flow rate through the respective opening. This is the main reason why the difference in the mass flow rate between the cases of 10 and 15 m/s wind speed is not so significant as in the other direction cases. Finally, the external wind direction affects the flow development inside the room, creating counterclockwise vortices. Figure  The fitted function depends on the type of window opening, building dimensions, and wind conditions. Thus, in the current study the best fit to the normalized data is achieved by using the Gaussian function provided by Equation (9): The constants a and b are selected equal to 1 and 0, respectively, in order for the Gaussian curve to pass through the point (0,1). Depending on the incident velocity that characterizes each group of computational runs, the value of c parameter that ensures the best agreement between the derived dimensionless curves and the Gaussian one is c = 0.4055, c = 0.4292 and c = 0.3361 for the wind speeds of v = 6 m/sec, v = 10 m/sec and v = 15 m/sec, respectively. The curve used for fitting for each specific velocity case along with the coefficient of determination, R 2 , is given in Figure 17  The fitted function depends on the type of window opening, building dimensions, and wind conditions. Thus, in the current study the best fit to the normalized data is achieved by using the Gaussian function provided by Equation (9): The constants a and b are selected equal to 1 and 0, respectively, in order for the Gaussian curve to pass through the point (0,1). Depending on the incident velocity that characterizes each group of computational runs, the value of c parameter that ensures the best agreement between the derived dimensionless curves and the Gaussian one is c = 0.4055, c = 0.4292 and c = 0.3361 for the wind speeds of v = 6 m/s, v = 10 m/s and v = 15 m/s, respectively. The curve used for fitting for each specific velocity case along with the coefficient of determination, R 2 , is given in Figure 17.
The constants a and b are selected equal to 1 and 0, respectively, in order for the Gaussian curve to pass through the point (0,1). Depending on the incident velocity that characterizes each group of computational runs, the value of c parameter that ensures the best agreement between the derived dimensionless curves and the Gaussian one is c = 0.4055, c = 0.4292 and c = 0.3361 for the wind speeds of v = 6 m/sec, v = 10 m/sec and v = 15 m/sec, respectively. The curve used for fitting for each specific velocity case along with the coefficient of determination, R 2 , is given in Figure 17. Because of the fact that the c parameter changes among the investigated cases, depending on the wind speed, it is necessary to formulate a mathematical expression that connects the c parameter with the velocity magnitude. This mathematical expression is given in Equation (10): where v is the wind speed in m/s. The final expression of the mass flow rate prediction as a function of the normalized wind direction θ and wind speed is described by the Equation (11): Because of the fact that the c parameter changes among the investigated cases, depending on the wind speed, it is necessary to formulate a mathematical expression that connects the c parameter with the velocity magnitude. This mathematical expression is given in Equation (10): where v is the wind speed in m/s. The final expression of the mass flow rate prediction as a function of the normalized wind direction θ and wind speed v is described by the Equation (11): Finally, Equation (12) provides the maximum mass flow rate for each group of the examined cases as a function of the incident velocity, with R 2 = 0.995: . m max (v) = 0.04599v (12) The coefficient of determination of Equations (10) and (12) is equal to 1 and 0.995, respectively.

Verification of Empirical Function
In order to validate the derived mathematical correlation, it is necessary to conduct additional simulations using different wind speeds to assess the agreement of the provided values by the mathematical expression against the CFD results. In this framework, three values of velocity magnitude have been selected; the first one smaller than 6 m/s, the second one in-between the range of 6 and 15 m/s and the third one above the maximum selected value of 15 m/s. Table 2 shows the percentage relative error between the mass flow rate numerically calculated and the one estimated by the derived function in the case of 2.5 m/s wind speed. The relative error has been calculated by Equation (13). The maximum numerical errors can be seen in the two limit values of wind direction, i.e., 45 • and 215 • , and are equal to 38.2% and 8.8%, respectively. Contrary to the limit values of wind speed, the interval ones present very good agreement. This fact can be also observed in Figure 18, which presents both the derived graph of the function for this specific wind speed (blue line) and the CFD results (red dots): values of wind direction, i.e., 45 o and 215 o , and are equal to 38.2% and 8.8%, respectively. Contrary to the limit values of wind speed, the interval ones present very good agreement. This fact can be also observed in Figure 18, which presents both the derived graph of the function for this specific wind speed (blue line) and the CFD results (red dots):  The second wind speed that has been selected is equal to 12 m/s and is an interval value of the range between 6 and 15 m/s. It can be observed that in this case a better agreement between the CFD and empirical results can be achieved for all the wind directions, since the maximum relative difference is approximately equal to 3.5% (Table 3). This agreement is also noticeable in Figure 19, The second wind speed that has been selected is equal to 12 m/s and is an interval value of the range between 6 and 15 m/s. It can be observed that in this case a better agreement between the CFD and empirical results can be achieved for all the wind directions, since the maximum relative difference is approximately equal to 3.5% (Table 3). This agreement is also noticeable in Figure 19, which presents both the derived graph of the function for this specific wind speed (blue line) and the CFD results (red dots). which presents both the derived graph of the function for this specific wind speed (blue line) and the CFD results (red dots). For the last case, the wind speed is 16m/s. In Table 4 the relative error between the estimated and CFD results can be seen. When the wind direction is parallel to the window i.e., 45 o and 215 o the relative error is significant, i.e., 35.5% and 36%, respectively. Moreover, an underestimation of mass flow rate can be observed for 90 o and 180 o as the velocity of 16 m/s does not belong to the range of For the last case, the wind speed is 16 m/s. In Table 4 the relative error between the estimated and CFD results can be seen. When the wind direction is parallel to the window i.e., 45 • and 215 • the relative error is significant, i.e., 35.5% and 36%, respectively. Moreover, an underestimation of mass flow rate can be observed for 90 • and 180 • as the velocity of 16 m/s does not belong to the range of selected values for Gaussian fitting. This fact can be also observed in Figure 20, which presents both the derived graph of the function for this specific wind speed (blue line) and the CFD results (red dots).  Based on the results, it can be concluded that the mathematical correlation cannot accurately predict the mass flow rate through the window opening when both of the following two arguments are valid: a) the wind speeds are out of the range of the values that have been selected for the function development and b) the wind direction is parallel to the window's surface. This can be attributed to the fact that the Gaussian function is symmetrical around the central value, while the developed flow field presents recirculations and flow deflection in front of the window that induce uncertainties and dissimilarities, so the differences between the Gaussian function and the CFD results become more significant. Thus, it is difficult to predict the mass flow rate with reasonable accuracy in the cases of wind direction parallel to the window and wind speed out of the range of the selected values for the derivation of the function. These results are consistent with the findings of Wang [23], who observed that only when the mass flow is perpendicular to the tilted window the proposed semi-empirical model agreed with his CFD simulations. Furthermore, this model refers to opening angles 30 o 45 o , while in the current study an angle of 5.8 o is examined. This is important because due to the complicated geometry it is not clear the location of the neutral plane which is a term in the semiempirical model for the calculation of the ventilation rate. Thus, the results of this study cannot be compered by the semi-empirical model.

Conclusions
This paper presents a simple, versatile methodology for the development of an empirical equation, which can provide the air mass flow rate imposed by single-sided wind-driven ventilation of a room, as a function of external wind speed and direction, using the results from CFD simulations. k-w SST turbulence model and Boussinesq approximation have been used for the simulation of turbulence and buoyancy forces, respectively.
In order to achieve the derivation of a function from CFD simulations for prediction of the mass Based on the results, it can be concluded that the mathematical correlation cannot accurately predict the mass flow rate through the window opening when both of the following two arguments are valid: a) the wind speeds are out of the range of the values that have been selected for the function development and b) the wind direction is parallel to the window's surface. This can be attributed to the fact that the Gaussian function is symmetrical around the central value, while the developed flow field presents recirculations and flow deflection in front of the window that induce uncertainties and dissimilarities, so the differences between the Gaussian function and the CFD results become more significant. Thus, it is difficult to predict the mass flow rate with reasonable accuracy in the cases of wind direction parallel to the window and wind speed out of the range of the selected values for the derivation of the function. These results are consistent with the findings of Wang [23], who observed that only when the mass flow is perpendicular to the tilted window the proposed semi-empirical model agreed with his CFD simulations. Furthermore, this model refers to opening angles 30 • 45 • , while in the current study an angle of 5.8 • is examined. This is important because due to the complicated geometry it is not clear the location of the neutral plane which is a term in the semi-empirical model for the calculation of the ventilation rate. Thus, the results of this study cannot be compered by the semi-empirical model.

Conclusions
This paper presents a simple, versatile methodology for the development of an empirical equation, which can provide the air mass flow rate imposed by single-sided wind-driven ventilation of a room, as a function of external wind speed and direction, using the results from CFD simulations. k-w SST turbulence model and Boussinesq approximation have been used for the simulation of turbulence and buoyancy forces, respectively.
In order to achieve the derivation of a function from CFD simulations for prediction of the mass flow rate, it was necessary to use three wind speeds, namely 6, 10 and 15 m/s. In each case of wind speed five different wind directions were simulated. The normalized mass flow rates were fitted using a type of Gaussian function. The validation of the empirical function has been performed by conducting additional simulations with wind speed equal to 2.5, 12 and 16 m/s. In contrast to the case of the velocity of 12 m/s, whose predictions have a very good agreement with the simulation results, the other two cases present significant relative error when the airflow is parallel to the window. With these wind directions, the CFD results have showed the development of recirculations near the window and the deflection of the flow from the building. Since these phenomena are complicated and the function cannot accurately take them into account, the relative error between the simulation and the prediction in these cases is increased. Moreover, the selected velocities of 2.5 and 16 m/s do not belong in the range of the values which the correlation is based on. The broad range of wind speeds that have been examined and the non-symmetrical building formation (contrary to the symmetrical conditions of pilot simulations) distinguish the present work from previous publications, that use CFD simulations for an empirical correlation for a shorter range of wind speeds and symmetrical conditions. However, an interesting follow-up work of this study could be to use experimental data for further accuracy of the correlation, because in the present case the experimental values are not valid, since the weather station is located inside the developed boundary layer.
Author Contributions: In this work, the collaboration of many contributors were necessary. P.S. was responsible for the investigation, P.D. was responsible for the methodology, D.R. was responsible for the writing (review and editing) and N.N. was responsible for the project administration.
Funding: The analysis has been performed in the framework of PLUG-N-HARVEST research project, grant number 768735, funded by EU's Horizon2020 program.