Modelling urban airflow and natural ventilation using a GPU-based lattice-Boltzmann method

Simulation of urban air ﬂ ow and ventilation potential is desirable for building design, however the complex and transient nature of ﬂ ows in urban environments makes this a challenging task. This study aims to evaluate the capability of a lattice-Boltzmann method (LBM) code deployed on a graphical processing unit (GPU) using a large eddy sub-grid turbulence model for cross-ﬂ ow ventilation of an idealised cubical building at wind-tunnel scale. ANSYS Fluent is used as a numerical comparison. Façade pressure and ventilation of the cube are investigated for parallel and perpendicular wind directions with the building in isolation and regular array format. Pressures, velocities and ventilation rates are compared to experimental data from wind tunnel and full-scale experiments of the Silsoe cube. Simulations compare favourably with experimental values and between each other. When the cube was surrounded by other cubes, simulations suggest that vortex shedding from up-wind buildings provides pulsating ventilation, improving air ﬂ ow ingress in the parallel wind cases. A parametric study showed that doubling surrounding building height had a small negative effect on ventilation but was mitigated by high levels of downdraft and ﬂ ow ﬂ uctuations in the vertical plane. Comparatively, doubling the central building height had a net positive effect but caused high internal airspeeds for both angles. The LBM code running on one GPU was several orders of magnitude faster than Fluent with similar accuracy. Simulation time using the LBM approach was several orders of magnitude lower than Fluent. © 2017 The Authors. Published by Elsevier Ltd.


Introduction
Naturally ventilated buildings are common worldwide and are advocated as part of sustainable and resilient infrastructure development [1], however the relationship between external airflow and indoor air quality is still an area of much debate and challenging research [2].Even for simple building geometries, naturally induced airflow patterns can be highly complex [3].Since modelling outdoor and indoor air is a problem of scale/time, where large eddies dominate external flow [4] but smaller eddies dominate inside buildings [5], care needs to be taken when modelling these phenomena together [6].The traditional simulation approach using finite-volume computational fluid dynamics (CFD) to capture the detail of urban flows and transient behaviour requires increasingly substantial computing resources [7].However, graphical processing units (GPUs) are becoming increasingly powerful with massively parallel capabilities, and therefore lend themselves to the airflow simulation process using a novel lattice-Boltzmann method (LBM) [8,9].
This study explores the applicability of a non-traditional lattice-Boltzmann method [10] on the graphical processing unit [11] using an LES Smagorinsky sub-grid model for simulating natural ventilation evaluated against experimental wind-tunnel and full-scale experimental data.This is compared against transient finitevolume simulations carried out using ANSYS Fluent 16 (ANSYS Inc., Canonsburg, PA, USA).A parametric study considers airflow in and around a cubical hollow building with central façade openings, which is investigated in isolation and in an array format for parallel and perpendicular approach flows.Results are used to explore the effect of wind and neighbouring buildings on façade pressures, velocities and ventilation rates when the building is smaller or larger than the surrounding structures.

CFD modelling in urban environments
Computational simulations are becoming increasingly popular in the area of urban physics such as in Refs.[7,12e14].In particular, urban pedestrian wind environments [6,15,16] often predominate early building design stage, with modelling wind-driven rain coming later on [17].More recently, the interest in pollutant dispersion [18] within the urban area from car exhausts or toxic releases is becoming increasingly prevalent [12,13,19,20].The ability to be able to model air ingress into buildings through openable façade elements, such as windows, is critical to informing ventilation potential as well as pollutant exposure [6,21].However, one of the major drawbacks of current CFD approaches is the tradeoff between computational time and accuracy, with results of airflow patterns often taking days or weeks to calculate.Emergency responders, risk analysis and building estates management require an almost instantaneous response to changes in the environment such as chemical leaks, infection transmission or pollutant dispersion, and so real-time (or close to real-time) simulations are of particular interest.There is also a need for simulation to be able to capture transient changes as well as provide mean parameters.
Vortex shedding, recirculation and reattachment [22] are commonplace in the urban environment and are well known to pose difficulties to traditional CFD turbulence modelling practices employing Reynolds' Averaged Navier-Stokes (RANS) [12].Despite known inaccuracies [23e25], indoor air modelling relies heavily on eddy viscosity k-ε turbulence closure models because of their computationally inexpensive two-equation approach [26].External flow warrants a time-dependent approach and highly computationally expensive large eddy simulation (LES) models have the ability to capture key flow structures [5] but have a slow turnaround time.The necessity for coupling indoor and outdoor airflow simulations for naturally ventilated buildings is becoming increasingly clear, and a combined approach to turbulence is necessary.
Over the last decade, lattice-Boltzmann methods (LBM) have emerged as a new technique for the computation of fluid flow phenomena [8,10].This approach has proven to be an efficient tool, which is especially well suited to a parallel implementation [10].Contrary to the conventional method of solving the Navier-Stokes (NS) equation, the lattice-Boltzmann methods (LBM) represent a bottom-up approach by starting at a discrete microscopic model.By construction, this conserves the desired macroscopic quantities such as the hydrodynamic variables of mass and momentum, which behave according to the dynamics described by the NS equation [27].
The LBM originates from the lattice gas automata method and can be regarded as an explicit discretisation of the Boltzmann equation [10].In principle, an unfeasible number of particles would be required to represent the true nature of fluid flow so statistical distributions are employed on regular lattices instead.The LBM has several advantages over the NS equations, such as its numerical stability [11] and accuracy [8] and the capacity to efficiently handle complex geometries.Thus the LBM is an explicit numerical scheme with only local operations.It has the advantage of being easy to implement and is especially well suited for massively parallel machines like graphics processing units (GPUs) [9].As a consequence of the GPU's architecture, computational speed-up of many orders of magnitude can be achieved in comparison to traditional CPU simulation [10].
Recent LBM modelling of indoor air quality [18,28,29] has shown the validity and applicability of the method, particularly for modelling indoor contaminant and particle dispersion [18].The LBM has been shown to capture important indoor flow features in real-time and compares well to simulations carried out using finite-volume methods [8,10].Due to the large scale of external airflows, this methodology also lends itself well to computation on GPU.In particular, Obrecht et al. [9] show the applicability to flow around solid arrays of cubical structures with encouraging results.However, a gap exists for the lattice-Boltzmann method in coupling indoor and outdoor airflows, including validating the approaches against experiments.

Methodology
The first stage of the study applied the two simulation approaches to a series of scale model scenarios, and compared the results to published data from wind-tunnel experiments.

Mesh sensitivity analysis: Meinders et al. [30]
A mesh sensitivity analysis was conducted initially using experimental data from an isolated solid surface-mounted cube in a boundary layer wind-tunnel published in Meinders' et al. [30] (Case 0).A wind-tunnel of height h ¼ 50 mm and width 600 mm with cube size H ¼ 15 mm (h/H ¼ 3.3) was placed on one of the walls along the channel centre-line.The domain can be seen in and shown in Fig. 1.In-flow conditions were controlled by tripping the boundary layer 75 cm upstream of the cube leading face (see Fig. 1) and vertical velocity profiles (Fig. 5) are given in section 3.3.2.A uniform profile of 4.4 m/s was used at the inlet with zero pressure gradients at the outlet.

Study scenarios
Subsequently, the scale model computational set-ups were then expanded to consider a cube with facade openings representing an idealised building.Airflow patterns were investigated in cases 1e4 for a cube, width ¼ 0.2 m, length ¼ 0.2 m, height ¼ 0.16 m with openable façade elements (width ¼ 0.092 m and height ¼ 0.036 m).Simulations were conducted for cross-flow ventilation with the wind approaching perpendicularly to (0 ) and parallel to (90 ) the openings.In both cases the cube was considered in isolation (case 1a/b, see Fig. 2a) and in a nine-cube array format (cases 2e4 see Fig. 2bed).Cases 1a and 2a are similar to Karava et al. [31] and Tominaga et al. [32] and were therefore compared quantitatively against their experimental data.The effect of surrounding building height was investigated by doubling the array height in case 3a and b (see Fig. 2c).The central cube was then placed on a plinth, whereby doubling its height but maintaining its volume; this was investigated in cases 4a and 4b (see Fig. 2d).
The computational domains for all cases were similar to case 0, based on dimensional guidelines set out by Blocken et al. [33]: 4H upwind, 3H either side, 7H downwind of the cube and 4H high.Table 1 describes the case names and description of all numerical experiments.
Fig. 1.Solid isolated cube set-up in boundary layer flow for Case 0, similar to Meinders et al. [30].
Following the scale model investigations, both computational approaches were applied to a full-scale ventilated cube scenario.This study was carried out using the "Silsoe cube", a 6 m Â 6 m x 6 m metallic cube situated in an open field which is an idealised experimental scenario that has been used for wind engineering research since the early 1990's [8].Experiments were carried out to characterise external and internal flow when the cube had openings (1 m Â 0.4 m) located on opposite façades (see Fig. 3).Façade pressures, P, were measured using 32 differential pressure sensors (Honeywell, USA) located as in Fig. 3a.A reference velocity, U ref and reference pressure P ref were measured at a mast 41 m upstream of the cube in the prevailing WSW wind direction, at a height of 6 m (roof height).Pressure coefficients, C p , were determined using equation (1): Ventilation rates, Q, within the cube were calculated from the windward pressure coefficients based on the opening areas, A, using equation ( 2) with an experimentally derived discharge coefficient specific to this opening: C d ¼ 0.658.
Measurements of wind speeds and façade pressures on the Silsoe cube were made constantly from 30 th May 2015 until 7th July 2015.Wind angles were averaged into ±2.5 segments.All variables were averaged over 30 min periods.Further details can be found in King et al. [34].

Lattice-Boltzmann method
The lattice-Boltzmann method is a discrete version of the Boltzmann transport equation: where f (x, e, t) describes the evolution in time of the distribution of one particle in phase space, F is the external force field, m the mass of the particle, and U the collision operator.It is well known, using ChapmaneEnskog expansion, that the incompressible NaviereStokes equations can be recovered from the Boltzmann transport equation for small Knudsen numbers.Discretisation occurs both in time, with constant time steps Dt, and phase space, generally using a regular orthogonal grid of mesh size Dx and a finite set of N þ 1 particle velocities x i with x 0 ¼ 0 (equation ( 4)).The latter is commonly a subset of the velocities linking any node of the grid to its nearest neighbours as the D3Q19 stencil that was used for current simulations (see Fig. 4).
The evolution of the distribution functions using the Bhatna-gareGrosseKrook [35] collision is described by: where viscosity (n) is related to particle relaxation time (t), Density and velocity were recovered as follows: where c ¼ Dx Dt is the lattice speed.The local equilibrium functions are calculated by Taylor expansion: Pressure was recovered through P ¼ c 2 s r, where c s ¼ c=√3 is the speed of sound.

Turbulence modelling
The smallest Eddies were modelled through the Smagorinsky sub-grid model, previously implemented into the LBM in Delbosc et al. [10] This model uses a positive turbulent eddy viscosity, n t , to represent small scale energy damping.This viscosity n t is computed from the local stress tensor S ab as follows: where C > 0 is the Smagorinsky constant, which typically takes on a value between 0.01 and 0.2.D represents the filter width and is the magnitude of the local stress tensor: The total viscosity of the fluid equals the sum of the physical viscosity and the eddy viscosity n total : In the LBM, the effect of the eddy viscosity is incorporated into a local relaxation time t s given by: This modified relaxation time is then used in the relaxation process of the LBGK equations; so each node of the lattice relaxes at different rates.The local stress tensor is relatively easy to compute within the LBM, compared to traditional schemes (requiring finite difference computations) and can be computed locally from the non-equilibrium stress tensor: where a and b runs over three-dimensional space.The intensity of the local stress tensor S ab is calculated thus:

Domain and boundary conditions
A uniform velocity profile was imposed at the inlet.In the scalemodel cases, an atmospheric boundary layer was established by using surface-mounted ribs upwind of the cubes to induce turbulence in the flow (see Fig. 1).The mean streamwise velocity of the approach flow was compared in Fig. 5 against experimental values from Meinders et al. [30].The streamwise velocity (u/U ref ) and velocity fluctuations (u' 2 /U ref ) were quantitatively comparable to those found by Meinders et al. (1999) (see Fig. 5).Quantitative comparison for each point shows variation to be less than 10% throughout.
In the full-scale Silsoe case, the approach flow was based on the logarithmic boundary layer profile (equation ( 16)).In accordance  with Richards and Hoxey [36] we used a reference velocity (U ref ) of 10 m/s at z ¼ 6 m and ground roughness length z 0 of 0.01 m to define the inlet velocity condition.

UðzÞ u
here u * is the turbulent friction velocity (0.656 m/s) and k is von Karman's constant (0.41).
The outlet condition was obtained by imposing a constant massflux outflow.No near wall treatment was implemented for the current LBM based LES simulation and so sensitivity close to the wall was investigated as part of a sensitivity analysis.The no slip conditions of the macroscopic velocities are obtained through the bounce back conditions applied on the distribution functions at the walls.This implicitly imposes the linear profile near the wall.To ensure temporal convergence to a statistically steady state, the variables were averaged over time from 50T 0 to 650T 0 , where T 0 ¼ H/u 0 is the turn-over time, H is the cube height and u 0 is the maximum inlet velocity.
The LBM code was run on the graphical processing unit (GPU) as detailed in Delbosc et al. [10], in this case an Nvidia GTX 780Ti, using the CUDA language.Visualisation was made in OpenGL and Paraview (Kitware, USA) for post-processing.

Fluent model
ANSYS Fluent 16.2 was used to perform the finite-volume based CFD analysis, using a transient isothermal approach with the LES Smagorinsky sub-grid scale turbulence model.The value of the Smagorinsky constant was chosen to be 0.02 based on previous studies [9].Domains for all simulations were identical in dimension to the LBM models.
The sides and top of the computational domain were modelled as no-slip walls and were far enough way not to impact the flow around the cube, which was tested during a sensitivity analysis [17].The outlet was set at 0Pa gauge pressure in all cases with zero gradient for all other variables.
Horizontal inhomogeneity in the vertical plane was first investigated by performing a simulation on an empty domain to establish appropriate location of the cube.It was found that three building heights in from the inlet was sufficient to reduce decay of turbulent properties in the approach flow.
Central differencing schemes were used for both convective terms and diffusion terms.For time integration, a central bounded second-order accurate scheme was utilised with a time step of 0.001s.This is equivalent to a maximum Courant-Friedrichs-Levy (CFL) value of 2 based on Dx ¼ H/32 and U ref ¼ 10 m/s.This CFL value was used in order to maintain equality with the LBM simulation.
Initially, a steady state simulation was created initially using the k-u SST model (SIMPLE algorithm) and then after 1000 iterations the model was changed to the LES model.Additionally, the PISO (Pressure Implicit with Splitting of Operators) algorithm was employed for pressure-velocity coupling in Fluent.30 iterations per time-step were utilised in Fluent to ensure statistical convergence of flow variables to less than 10% variation as the CFL number is above the recommended maximum value of 1 [37].
To ensure temporal convergence to a statistically steady state, the variables were averaged over time from 50T 0 to 650T 0 .Computation time for each Fluent simulation was approximately two days on 16 AMD cores and 32 GB RAM at the ARC2 high performance computing facility at the University of Leeds.

Meshing
Meshing for the Fluent simulations was conducted in ANSYS Meshing using unstructured hexahedral cells.Solution sensitivity to mesh size was explored previously [34] leading to selection of a mesh containing ~1.8 million cells.
Meshing for the LBM simulation was created by dividing the flow domain into cubical cells of uniform size.The cells were tagged according to whether it represented a solid obstacle or a fluid domain.We used Lua scripting language to build a custom toolkit to generate the structured hexahedral cells along with all the boundary conditions.Mesh sensitivity analysis was carried out using the Case 0 scenario, by comparing simulations to data presented in Meinders et al. [30] for a solid cube (H ¼ 0.015 m) in a turbulent channel.Three subsequent mesh refinements were tested (see Table 2), where the cell size is a ratio of the cube height.

External airflow
Fig. 6a and b shows a sensitivity analysis for the two Fig. 5. Averaged approach flow characteristics from Meinders et al. [30].
computational methods compared to wind-tunnel experimental data [30] for approach flow of the isolated Case 0 cube at z/H ¼ 0.5 and z/H ¼ 0.3 respectively.Results show that both the Fluent and LBM code capture the decaying speed as the flow approaches the cube wall, and equally are able to predict the velocity in the downstream areas.
Pressure is recovered using the relationship P ¼ c 2 r, where c 2 is the speed of sound and r is density.Pressure coefficients (equation (11)) are compared between experimental and numerical simulations plotted along the centreline of the cube in Fig. 6.These are calculated using P ref and U ref are reference pressure and reference velocity at a point upwind of the cube at eve height.CIBSE guide A [38] C p values are superimposed as indicative values used in conventional engineering ventilation calculations (see Fig. 7).
C p variations between mesh resolutions for the LBM code, particularly on the front face of the cube arise because no wallfunction is used in the code.Consequently, flow variables are sensitive to the mesh right up to the wall.Far away from the cube (Fig. 6) little improvement was obtained for velocity data above a mesh ratio of H/16.However, pressure predictions on the cube's façade benefitted from using a smaller cell size ratio of H/32 and this was used henceforth.
Fig. 8 shows average and standard deviation of velocity magnitude plotted parallel (line x in Fig. 8) and perpendicular (line y in Fig. 8) to the prevailing wind within the array for case 2a and case 3a at z ¼ 0.08 m (half way up the cube).The averaging period is 600H/U ref (s).Experimental data from the Architectural Institute Japan (AIJ) [39] is included as a comparison, however it there are slight differences in wind tunnel configuration, so this is used as a guide for trends and behaviours rather than quantitative validation.
Results show that both the Fluent and LBM models predict velocities that are similar to those measured experimentally.Along line y, flow accelerates as it enters the canyon for both models and concurs in trend with the experiment.The LBM code shows slightly slower speeds past the central building, when the buildings are the same size but matches experiments well when it is twice as high (Fig. 8d).A 3% difference exists between LBM and experiment at point of closest comparison, which increases up to a maximum up 23% at end of the canyon.Fluent represents the horizontal trends of the flow (line x), on average, better than the LBM code (6% vs 17% error with respect to experiment).Flow speed is noticeably higher behind the central building in the experimental data but this is not picked up in the Fluent or the LBM results, although the speeds are within the variation shown by both codes (Fig. 8aeb).

Internal airflow
Fig. 9aec shows time-averaged normalised velocity plots inside the ventilated scale-model cube for both LBM and Fluent codes, and compared to experimental data from Tominaga et al. [32].Velocity is plotted on three vertical lines x/D ¼ 0.125, 0.5 and 0.875, where D is the horizontal length of the cube for the isolated cube case1a (Fig. 9d).The angle of the flow dips as it enters the window and discharges slightly towards the ground, which is represented by the decay in velocity along the streamwise direction, and is apparent in the LBM contour plot in Fig. 9d.The LBM code tends to overestimate the velocity in the peak regions, whereas Fluent is more conservative.Results from both codes show the same trends as the experiments, although there are differences with respect to experimental values [32] in both codes close to the wall.
Experimental data for normalised internal velocity was not available for the array and so comparison is made against turbulent kinetic energy (k/U 2 ref ) [32], as shown in Fig. 10.According to CIBSE guide A [38], kinetic energy should be preserved after passing through large openings (as in this case).Both experimental data and numerical predictions show that this is not the case, and results are more in-line with expected behaviour for a small opening [38].Both models show a good comparison with experimental data, although the LBM code tends to slightly over predict turbulent kinetic energy values close to the inlet and under predict those close to the floor.
Fig. 11 shows contours of normalised velocity magnitude ( U =U ref ) on a central vertical plane for 0 wind cases 1a, 2a, 3a and 4a and a horizontal plane at mid-window height for 90 wind cases 1b, 2b, 3b and 4b.All images show the results from the LBM simulations.
Velocity inside the cube decreases as the wind shifts from perpendicular (case 1a) to parallel (case 1b) to the opening (Fig. 11a and b).When the wind is perpendicular, flow enters the cube where the jet subsequently dips towards the ground.This behaviour has been seen previously [31] where the jet is only quasi-impinging.Similar behaviour exists where the cube is twice as high as surrounding buildings in case 4a (Fig. 11c).Experimental anemometry data [32] support this feature as shown in Fig. 9.This jet is not as evident when the cube is smaller and surrounded by other cubes (case 2a) (Fig. 11c and d) where the other structures shield the cube from direct wind; this has also been found experimentally [14].As a consequence, even at 0 wind (perpendicular), only a moderate short jet is seen.Fig. 11c shows how vertical downdrafts from the taller outer cubes create reversed airflow patterns within the central cube leading to ventilation occurring from the rear window (Fig. 11c).This reduces the internal jet's speed, however it creates a ventilation regime that is not intuitive.
A small lateral pulsing phenomenon in the wind is found at 90 (see Fig. 12a for a snapshot of instantaneous flow) when the cube is in isolation, and this becomes magnified in the array case (Fig. 11d).In isolation, the flow detaches from the leading side edges of be cube and re-attaches towards the back half of the window (Fig. 11b and c), which is known to be a typical phenomenon [40].Fig. 12bec shows flow along the right and left window, split into u, v and w components.In the array case at 90 , the same phenomenon is repeated (see Fig. 11d).This is probably due to strong vortex shedding from up-wind cubes impacting on façade pressures of down-wind cubes [36].Fig. 12a shows velocity contours for horizontal wind movement (v), and These results highlight the coupled motion of the air blowing in one opening and being sucked out the other.The vertical component (w) of wind can be seen to be the same order of magnitude as u and v for all cases at the window opening (see Fig. 12b and c).This is found to exert a significant effect on flow and has been reported previously [41].The phenomena seems to be strongest for the cases when the outer cubes are twice as high as the central cube (case 3a and case 3b).

Ventilation potential
Comparison of the numerical approaches with full-scale experiments using the Silsoe cube was carried out by determining normalised ventilation rates within the cube, For the experimental data, Q is as defined by equation ( 2) using the average of four pressure taps around the windward window.A is the opening area, and U ref the streamwise wind speed collected experimentally for all wind directions at eave height measured 41 m upstream of the cube.Equivalent values from the numerical results were calculated using the area averaged air speed at the windward opening in the perpendicular cases and an average (ignoring negative spanwise velocities) from both left and right windows in the parallel cases.Fig. 13 shows the distribution of Q* with wind angle determined from the experiments, together with values calculated from the Fluent and LBM simulations at the two wind angles, 0 and 90 .Error bars on the CFD results represent root mean squared errors either side of the mean.In both cases the simulated values compare well to the experimental data.Similar results are found for both LBM and Fluent simulations, and both models predict ventilation rates which are of a comparable magnitude to the experimental results and show similar trends with angle.Differences for case 1a of less than 5.6% and 10% for Fluent and LBM respectively, were found between experimental median and computational prediction.On the other hand, relatively few experimental points were available for the case 1b case and therefore comparison was less accurate with up to a 60% difference.
Both Fluent and LBM results were compared against CIBSE guide A suggested values for each of the cases.Fig. 13 shows these Q* values compared for all simulation cases.As the incident wind moves from parallel (90 ) to perpendicular (0 ), the CFD model air change rate increases for all cases.This is most noticeable in the isolated scenario (case 1a), and concurs with findings in Chu et al. [28].Such patterns are less noticeable in the array case of equal height buildings (case 2a), where experimentally, a lower ventilation rate and range is found, but computationally a high level of variation is predicted.Parallel wind cases show low Q* values, with the isolated case (case 1b) exhibiting the lowest LBM prediction value.When in an array format (case 2b), ventilation rate is predicted to be almost double the isolated case (case 1b), but variation was higher.In general, upwind structures block the prominent effects of oncoming wind causing ventilation rates to drop.Under these conditions vortex shedding from upwind cubes are seen to form a pulsating ventilation pattern.Therefore, variations increase in air change rates which are then driven predominantly by local fluctuations and turbulent structures (see Fig. 12a).Fig. 6.LBM mesh sensitivity analysis of U/U ref at (a) z/H ¼ 0.5 cube heights upwind, (b) z/H ¼ 0.3 downwind.Simulations compared to experimental data from Ref. [30] and CBISE guide A values [38].Fig. 7. C p values calculated from Fluent and LBM models along the centreline of the case 0 cube, compared to experimental [30] and CBISE guide A values [38].
Q* decrease between 32% and 69% compared to the evenly spaced array when the outer buildings are taller than the central one for the perpendicular and parallel cases respectively.On the other hand, results suggest that Q* may increase up to 187% when the central building is taller than the surrounding ones under perpendicular flow.Less than 6% difference is seen between these cases under parallel flow.In addition, the parallel wind case (case 4b) appears to behave in a similar manner to the isolated case (case 1b), which is in line with Chu et al.'s findings [28] (see Fig. 14).

Computational time comparison
Computational time for calculating one second of flow was compared between the two codes for both the isolated (case 1a) and the array case (2a).Fig. 15 shows a direct comparison between both codes.Times are based on calculation time recorded by Fluent and the LBM internally and does not include rendering times.

Discussion
The experimental and CFD results show that mean flow structure and turbulence statistics depend significantly on the presence of additional structures as well as the wind angle.Unsteady transient effects are important, especially in the lower urban canopy layer where turbulent fluctuations dominate over the mean flow (Lim et al., 2009).Measurements in a full-scale set-up of this kind  are challenging, not least due to the variation of wind direction during averaging periods.
Both numerical codes are able to capture the pressure coefficients of the small-scale solid cube (Fig. 4).However, the roof has always presented a complicated area to model due to turbulent shear layers and re-attachment zones [42e44], which is apparent in results for both CFD codes.However, this is likely to be the least important face for ventilation calculation in most buildings.Adding a dynamic sub-grid model and implementing wall functions in LBM, which will be done in the future, should capture the wall boundary layers more accurately.Inclusion of low-Re wall functions may also improve Fluent results, an inclusion of more realistic surface roughness and dynamic mesh refinement may allow for increasing approximation between both codes and experimental results, especially for exterior surface pressures.
CIBSE guide A [38] provides some insight into the dependence of C p on the building location but the orifice equation for calculating ventilation does not account for varying façade pressures with height.When the canyon gap is equal to the building width in the cross-wind direction the building standards suggest using isolated building C p values for ventilation calculations.It must be highlighted that Fig. 13 suggests that CIBSE [38] provides an overestimation of the C p for arrays.
For all configurations tested, two distinct ventilation flow regimes were found: quasi-impinging jet and pulsating flow.When the openings face the wind, an impinging jet dominates the indoor Normalised air change rate, Q*, appears slightly higher in the array case with flow parallel to windows (case 2b) in comparison to the perpendicular flow case of the same set-up (case 2a).This is largely due to the vortex shedding from upstream buildings causing turbulent eddies to enter the façade openings, as seen by high levels of variation in spanwise velocity.Upwind buildings block the prominent effects of oncoming wind, which, in general cause ventilation rates to drop.Variations increase in air change rates and these are highly dependent on local turbulent structures, which can be dominated by v' and w'.
In general, Fluent tends to predict higher Q* values than the LBM code for perpendicular cases and predict lower values for parallel cases.CIBSE guide A comparison values are meant as guidance as Q* are calculated based on an isolated building or an urban environment (array cases).Root mean square errors are generally higher for Fluent than LBM throughout.One reason for this may be that this particular case exhibits typically higher v' and w' values than all others.As a consequence, Fluent values may need to be used conservatively when predicting ventilation potential.
Higher internal velocities are displayed during perpendicular wind conditions for the isolated case.This in turn may reduce internal mixing as the air short-circuits between the windows [25].Buildings in central array locations could benefit from cross-flow ventilation to improve air change effectiveness and rely on  upwind surrounding buildings to provide pulsating flow on days when the wind is parallel to the openings.As a disadvantage, higher indoor air speeds may be found.Feasibility will also depend on the building internal design and other considerations such as noise and external air pollution.Buildings that are smaller than their surrounding buildings may be classed with a relatively low ventilation potential but may benefit from downdrafts to increase local turbulence and hence ventilation rates.This may also result in unexpected reversals of ventilation airflow patterns within the building.

Future considerations
The authors acknowledge that the cubical structures tested with large openings are not fully representative of buildings, which may also include complexities such as eaves and pitched roofs.It does however allow a simplification of the urban environment which allows for individual parameter sensitivity analysis.Only isothermal cases were considered, as the focus of the work was on wind driven ventilation and flow effects, but it is acknowledged that heat in the urban environment may have significant effects on airflow at low wind speeds [7].As a consequence, high wind speeds of 10 m/s were modelled in the simulations, and in the full-scale experiments care was taken to ensure comparisons were made under conditions where thermal effects were minimal.In addition, uneven building height has been shown to exert a significant effect on urban airflow patterns [45], but the effect on ventilation potential has yet to be ascertained.
Mean flow and turbulence structure are found to depend significantly on surrounding cubes.Unsteady transient effects are important, especially in the lower canopy layer of urban environments where turbulent fluctuations dominate over the mean flow [2,14].Skimming flow dominates in the array scenario when the incident flow is at a right angle to the street canyons and all buildings are of the same size.This seemed to restrict vertical air movement within the canyons and therefore is known to have implications for pollutant dispersal [46,47].This will be investigated in future research with the use of a passive scalar.
While it is an idealised study, that only considers wind-driven effects, the coupling of symmetrical arrays, with openable façade elements at full-scale provides a clear opportunity to investigate ventilation under realistic wind conditions.Speed-up by using additional GPU devices in the future will have a marked impact on increased performance [10] as at the moment the code is restricted to a single GPU.In comparison to the cost of high performance computing hardware such as additional CPUs, the cost of adding an extra GPUs is comparatively small.Full parallelisation of the LBM code could create real-time airflow simulations comparable in accuracy to traditional CFD.Incorporating grid refinement techniques may help speed of calculation.In addition, wall functions may be able to help capture complex flow patterns such as reattachment zones on roofs, and will be incorporated in future developments.Accurate real time full-building simulation could be a useful tool for future building designers, indoor air quality assessment, estates management and risk analysts.

Conclusions
This paper presents a CFD analysis using ANSYS Fluent and GPUbased LBM code compared against experimental measurements of airflow in and around a cubical building in isolation and within an array of similar buildings.The CFD modelling gives insight into the external flow patterns and effect of differing building height and wind direction on ventilation potentials.
The study shows that both Fluent and LBM approaches using the LES model, are capable of providing comparable results for pressures, velocities and ventilation rates in all cases.Both computational approaches are also shown to compare well to wind-tunnel and full-scale experimental data.The LBM code also reproduces the inflow jet found by Tominaga et al. [48], and therefore appears to be a good approach for modelling the outdoor-indoor coupling.In addition, the LBM code can significantly reduce simulation times compared to traditional CFD methods, while maintaining comparable accuracy.
Air change rates increase as the incident wind becomes perpendicular to the window for cross-flow cases.This is reflected in both experiment and computation simulation.The isolated cube showed comparatively similar ventilation rates to the array cube at this angle.Increasing surrounding building height may result in reversed interior flow patterns but vertical air fluctuations may help mitigate lower ventilation potential by pushing in fresh air from outside.CIBSE guide A calculations may overestimate ventilation rates in urban environments, and there is a need for further research to understand whether the findings from the idealised geometries in this study can be translated to real environments.

Fig. 8 .
Fig. 8. U/U ref plotted along a horizontal and vertical line at z ¼ 0.08 m for the array of cubes case 2a and case 4a and compared with experimental data from the AIJ [36].Error bars are one standard deviation.

Fig. 10 .Fig. 11 .
Fig. 10.Average turbulent kinetic energy (k/U 2 ref ) for the array case 2a (case 2a) for the LBM and Fluent models compared to published experimental data.

Fig. 12 .Fig. 13 .
Fig. 12. Case 1b.a) Contours of instantaneous velocity magnitude in the v direction on a horizontal plane at z ¼ 0.22m b/c) Velocity magnitude plots in three directions for the left and right window at z ¼ 0.22 m.

Fig. 14 .
Fig. 14.Q* predictions for the LBM and Fluent codes compared to experimental values determined using methods in CIBSE guide A [38].Error bars show root mean squared error.

Table 1
Case set-up and description.

Table 2
Cell count for the LBM mesh sensitivity tests.