WO2021112756A1 - Method and apparatus for simulating a temperature distribution of a surface of a built environment - Google Patents

Abstract

A computer-implemented method for modelling a surface temperature distribution of a surface of a built environment comprises obtaining input data comprising latitude, longitude, geometrical data indicative of a subdivision of the surface into a plurality of surface patches, and weather data. The method comprises iteratively, over time and over the surface patches, initialising a temperature of each surface patch; obtaining a sky model that describes a shortwave sky radiance distribution; using the input data and the sky model, computing the temperature of the surface patches according to a time-dependent energy balance equation that comprises a sum of separately modelled contributions from a plurality of heat transfer processes; re-initialising the temperature of each surface patch for the next time step to be the computed surface temperature for that surface patch for the current time step; and outputting the surface temperature distribution by aggregating the computed patch surface temperatures.

Description

METHOD AND APPARATUS FOR SIMULATING A TEMPERATURE DISTRIBUTION OF A SURFACE OF A BUILT ENVIRONMENT

TECHNICAL FIELD

The present disclosure relates to a method and apparatus for modelling a surface temperature distribution of a surface of a built environment, for example for the purpose of optimising the design of such environments.

BACKGROUND

In many locales, increased heat absorption due to urbanisation has had a significant impact. For example, in Singapore, the local temperature has risen twice as fast as the global average over the past 60 years. In the tropical dense urban environment of Singapore, the surfaces of artificial constructions absorb the energy from intense solar radiation in the daytime and emit radiative heat during the night to induce an urban heat island (UHI) effect. The heat stress due to radiant heat gain of the surface has a significant impact on human health and outdoor activity during the day and dramatically increases energy consumption from use of air conditioning.

In view of the above, it has become increasingly important to consider thermal comfort when designing urban spaces. To this end, mean radiant temperature (MRT) is one of the most important indicators of outdoor thermal comfort. MRT quantifies radiant heat exchange between the human body and the surrounding environment.

A number of different ways of measuring MRT have been developed, including through the use of globe thermometers, integral radiation sensors, two-sphere radiometers, and constant air temperature sensors. However, using measurements to determine MRT as a function of space and time around an urban area is impractical and costly. Accordingly, numerical simulation methods have been developed as an alternative. Accurate simulation of MRT requires modelling of radiation exchange, i.e. the energy transfer among surfaces in the form of electromagnetic waves over the shortwave (0.2-4.0μm) and longwave (4.0-10.0μm) regions of the spectrum. The shortwave radiation originates from the Sun and travels through the atmosphere and reaches the terrestrial surfaces. The atmosphere and the terrestrial surfaces are heated up by the shortwave radiation and emit heat at thermal infrared frequencies, also called the longwave radiation.

Many techniques have been developed for modelling radiation exchange, such as the view factor method, the discrete ordinate method (DOM), the Radiosity algorithm, and ray tracing techniques. Each of these techniques has unique applicability and shortcomings. For instance, the view factor method is suitable for quick evaluation with compromised accuracy, the DOM technique has high accuracy but suffers computational inefficiency for non-participating media (for example, because of needing to solve the full radiative transfer equation in every cell of a 3D scene), the Radiosity algorithm offers reasonable accuracy but limited scalability due to computer memory constraints on matrix dimensions, and the ray tracing technique has high accuracy but may lack the ability to model in three dimensions due to platform constraints. Further, existing techniques need to extract surface temperature values from measurements using handheld IR instruments or thermocouples at selected sites, estimate based on a simple scaling of a visible radiation exchange simulation, or estimate based on view factors. Moreover, existing techniques (such as Fluent™ of ANSYS, Inc.) that generate surface temperature based on surface energy balance are often coupled with computationally expensive CFD computations for modelling of convective heat transfer.

It is desirable therefore to overcome or alleviate one or more of the above difficulties, or at least to provide a useful alternative.

SUMMARY Disclosed herein is a computer-implemented method for modelling a surface temperature distribution of a surface of a built environment, the method comprising, using at least one processor: obtaining input data comprising: a latitude and longitude of the surface; material parameters for one or more materials of the surface; geometrical data indicative of a subdivision of the surface into a plurality of surface patches, each surface patch having a surface normal; and weather data; initialising a temperature of each surface patch; and for each of a plurality of time steps: obtaining, using the weather data, the latitude, and the longitude, a sky model that describes a shortwave sky radiance distribution; for each surface patch, performing a surface temperature computation operation comprising: using the input data and the sky model, computing the temperature of the surface patch for the current time step according to a time-dependent energy balance equation that comprises a sum of separately modelled contributions from a plurality of heat transfer processes; after all patches have been processed, re-initialising the temperature of each surface patch for the next time step to be the computed surface temperature for that surface patch for the current time step; and outputting the surface temperature distribution by aggregating the computed surface temperature of each surface patch after the plurality of time steps.

In certain embodiments, the heat transfer processes are shortwave radiation exchange, longwave radiation exchange, convection, and conduction.

In certain embodiments, the contribution from shortwave radiation exchange and/or the contribution from longwave radiation exchange is, or are, determined based on Monte Carlo ray-tracing.

The method may comprise generating an octree structure based on the input geometrical data; wherein said Monte Carlo ray-tracing is performed in accordance with the octree structure, such that for each ray, only surface patches contained in octree cubes in the path of the ray are used for a ray-patch intersection test.

In certain embodiments, each surface patch is characterised by multiband material parameters. The multiband material parameters of at least some respective patches may be different than the multiband material parameters of the other respective patches.

In certain embodiments, the contribution from shortwave radiation exchange includes a contribution from multiple diffuse scattering.

The sky model may be the Perez all-weather sky model, and the clearness index of the Perez all-weather sky model may be used to infer cloud cover ratio for determining the contribution from longwave radiation exchange.

A convection coefficient for determination of the contribution from convection may be determined by optimized fitting from unshaded surface temperature measurements, or by derivation based on wind speed from the weather data. Also disclosed herein is an apparatus for modelling a surface temperature distribution of a surface of a built environment, the apparatus comprising: a kernel module that is configured to: obtain input data comprising: a latitude and longitude of the surface; material parameters for one or more materials of the surface; geometrical data indicative of a subdivision of the surface into a plurality of surface patches, each surface patch having a surface normal; and weather data; initialise a temperature of each surface patch; and for each of a plurality of time steps: obtain, by a sky model component, using the weather data, the latitude, and the longitude, a sky model that describes a shortwave sky radiance distribution; for each surface patch, perform a surface temperature computation operation comprising: using the input data and the sky model, computing the temperature of the surface patch for the current time step according to a time-dependent energy balance equation that comprises a sum of separately modelled contributions from a plurality of heat transfer processes; after all patches have been processed, re-initialising the temperature of each surface patch for the next time step to be the computed surface temperature for that surface patch for the current time step; and output the surface temperature distribution by aggregating the computed surface temperature of each surface patch after the plurality of time steps.

The heat transfer processes modelled by the apparatus may be shortwave radiation exchange, longwave radiation exchange, convection, and conduction.

The apparatus may comprise a shortwave component that is configured to determine the contribution from shortwave radiation exchange based on Monte Carlo ray-tracing; and/or a longwave component that is configured to determine the contribution from longwave radiation exchange based on Monte Carlo raytracing.

The apparatus may also comprise an octree generator that is configured to generate an octree structure based on the input geometrical data; wherein the shortwave component and/or the longwave component is configured to perform said Monte Carlo ray-tracing in accordance with the octree structure, such that for each ray, only surface patches contained in octree cubes in the path of the ray are used for a ray-patch intersection test. In certain embodiments, each surface patch is characterised by multiband material parameters. The multiband material parameters of at least some respective patches may be different than the multiband material parameters of the other respective patches. The contribution from shortwave radiation exchange may include a contribution from multiple diffuse scattering.

In certain embodiments of the apparatus, the sky model is the Perez all-weather sky model, and the sky model component is configured to use the clearness index of the Perez all-weather sky model to infer cloud cover ratio for determining the contribution from longwave radiation exchange.

The apparatus may comprise a convection component that is configured to determine a convection coefficient by optimized fitting from unshaded surface temperature measurements, or by derivation based on wind speed from the weather data.

BRIEF DESCRIPTION OF THE DRAWINGS

Some embodiments of a method and a system for modelling of surface temperature of a built environment, in accordance with present teachings will now be described, by way of non-limiting example only, with reference to the accompanying drawings in which:

Figure 1 is a flow diagram of an example method for modelling of surface temperature of a built environment;

Figure 2 is a schematic illustration of heat transfer processes modelled and key features in the method of Figure 1;

Figure 3 is a system architecture of an example apparatus for modelling of surface temperature of a built environment;

Figure 4 is a block diagram of modules of the apparatus of Figure 3;

Figure 5 shows geometry of a scene with four selected grid points on top, west, east and mid-east building surfaces, used to test the method of Figure 1;

Figure 6 is a graph showing absorbed shortwave heat flux computed by a method according to certain embodiments, plotted against absorbed shortwave heat flux for Radiance, on the top, west, east and mid-east building surfaces between 0700hrs and 1800hrs on 21 June; Figure 7 is a graph showing a comparison of fvDOM with calculated longwave heat flux from a method according to certain embodiments, on the top, west, east and mid-east building surfaces over the period of one generic hour;

Figure 8 shows a breakdown of different components of heat transfer and surface temperature on a) top, b) west, c) east and d) mid-east building surfaces;

Figure 9 shows a breakdown of multiband radiation exchange on a) top, b) west, c) east and d) mid-east building surfaces;

Figure 10 shows surface temperature distribution at (a) 1030hrs, (b) 1330hrs, (c) 1630hrs on 21 June; Figure 11 shows a scene of two industrial buildings (a) top view, (b) side view; Figure 12 shows hourly averaged meteorological weather input for 20 February 2018; and

Figure 13 shows a comparison of field measured surface temperature and calculated surface temperature according to certain embodiments, at four sites on building I south wall, building II roof, west and south wall between 0700hrs and 2000hrs on 20 February 2018.

DETAILED DESCRIPTION Advantageously, embodiments of the present invention are able to model shortwave and longwave radiation exchange in a physically accurate manner, and furthermore, are able to accurately generate surface temperatures of buildings and other structures (such as roads, pavements, vegetation, etc.) in a built environment at iterative time steps, without requiring measurements, and in more computationally efficient manner than prior art approaches.

Embodiments make use of Monte Carlo ray-tracing (MCRT) coupled with the Perez all-weather sky model (R. Perez et al., "All-weather model for sky luminance distribution - Preliminary configuration and validation", Solar Energy, 1993, 50(3): p. 235-245) for modelling of the shortwave radiation exchange.

Additionally, the clearness index from the Perez sky model may be used for the inference of cloud cover ratio when deriving the sky emissivity for the isotropic sky longwave radiation. In some embodiments, MCRT is applied to compute longwave radiation exchange for physical accuracy, and geometrical and material flexibility. Advantageously, the incorporation of Monte Carlo hemispherical sampling diminishes the aliasing effect and increases physical accuracy, as the actual photon scattering and emission process are random in nature.

Embodiments of the invention are able to output the surface temperature at each time step using the time-dependent energy balance equation based on accurately modelled spectral radiation exchange and carefully chosen conduction and convection models catered to scenarios where the radiation exchange dominates over other modes of heat transfer. It is important to note that the modules in embodiments of the present invention are designed to be independent and can be customised to accommodate the actual solar and wind weather conditions. In short, the present invention enables the physically accurate and computationally efficient modeling of multiband radiation exchange and surface temperature distribution in a complex urban environment.

Referring initially to Figure 1, a method 100 will be described. Some key features of model components of the method are illustrated schematically in Figure 2 and described in further detail below.

A method 100 for modelling a temperature distribution of a surface of a built environment comprises a first step 110 of obtaining input data. The input data may include the following:

• Geometry data relating to the geometry of the surface of a 3D scene corresponding to the built environment to be modelled (including one or more building structures, and surrounding features such as ground and urban features such as pavements, roads, grass, trees, etc.). For example, the geometry data may be provided in a 3D format such as STL format, and may comprise a plurality of triangular patches (or patches of other shapes), each having a surface normal. The patches can be of non-uniform lengths and arbitrary orientations. The 3D surface geometry is particularly suitable for ray tracing as the ray-surface intersection can be resolved for an arbitrary patch size and orientation. • Material properties for the surface to be simulated. This may include, for each material from which one or more buildings in the built environment is composed, and/or from which surrounding surfaces such as trees, ground surfaces and the like are composed, the following: multiband shortwave albedo, longwave emissivity, specific heat capacity, material density, thermal conductivity and material thickness, in SI units. Some commonly used urban materials may be pre-stored in a material library including concrete, asphalt, glass, grass, and waterbody, etc.

• The date and time for the simulation to start, the time over which the simulation is to be conducted, and the duration of each iteration of a temperature computation operation (to be described in more detail below).

• The geographical location of the built environment, including latitude and longitude, which are used to compute the solar incident angle in a shortwave radiation exchange model.

• Weather data (for example, in EPW format, the ASHRAE International Weather for Energy Calculations (IWEC) weather data representing a typical meteorological year of the chosen location). The weather data may include global horizontal radiation (Wh/m²), direct normal radiation (Wh/m²), diffuse horizontal radiation (Wh/m²), dry bulb temperature (T_db), dew point temperature (T_dp) (°C), relative humidity (%), horizontal infrared radiation intensity (Wh/m²), wind speed (m/s), etc.

The weather data may be used for several purposes as part of embodiments of the method 100. For example:

• Direct normal irradiance (DNI, Ia_\r) can be obtained from hourly averaged direct normal radiation to compute a direct contribution from the Sun.

• Diffuse horizontal irradiance (DHI) can be obtained from hourly averaged diffuse horizontal radiation.

• DNI, DHI and T_dp may serve as inputs for constructing the Perez all-weather sky model.

• T_db, also called air temperature, can be used for initializing the patch temperatures before sunrise. It is also used to represent the air temperature and indoor surface temperature in the convection and conduction models, respectively. • The longwave heat flux emitted by the sky I _IR can be obtained in two ways: 1) directly adopt the hourly averaged horizontal infrared radiation intensity from input weather data, 2) compute based on Stefan Boltzmann's law using T_dP from input weather data, and the sky clearness index from the Perez sky model.

• The convection coefficient H can be obtained in two ways: 1) load from user-definable input file defined by fitting with GHI and the hourly averaged horizontal infrared radiation intensity from input weather data, 2) derive based on wind speed from input weather data.

Next, the method may comprise performing a series of pre-processing operations. One or more of the pre-processing operations may make use of at least some of the input data.

For example, at step 120, the method 100 may comprise generating a plurality of grid points from the geometry data. Typically, each grid point will have coordinates of the centroid of one of the patches of the surface. Additionally, each grid point may have the same normal vector as the surface patch with which it is associated.

At step 124, method 100 may comprise generating an octree structure, based on the geometry data. The use of octrees is described in A. S. Glassner, "Space subdivision for fast ray tracing," IEEE Computer Graphics and Applications, vol. 4, no. 10, pp. 15-24, Oct. 1984, the content of which is incorporated herein by reference.

In one example, the octree structure may be generated according to the following: -The surface patches of the whole scene are input.

-The whole scene is enclosed with a cube, which is then recursively subdivided into 8 equal sized child cubes, until the number of surface patches in each child cube falls below a predefined limit or the size of the sub-cube is smaller than a certain size. -For each octree cube, its voxel coordinates, size, its parent cube, and child cubes are stored.

The octree structure may advantageously be used to increase the computational efficiency of ray tracing processes that are implemented as part of shortwave and longwave radiation exchange modelling as will be described in more detail below.

At step 126, an initialisation operation may be performed. Each surface patch, and its associated grid point generated at step 120, may be assigned an initial surface temperature. For example, if the simulation begins before sunrise, each patch and grid point may be assigned the air temperature as the initial temperature. In some embodiments, the initial surface temperature may be assigned based on one or more measurements of the surface temperature.

Next, an iterative computation is performed for each of a plurality of time steps, beginning at a simulation start time, and ending after a predetermined duration. For example, each of N time steps may have a fixed duration equal to (end time - start tim e)//V.

At step 128, method 100 may comprise obtaining a sky luminance prediction model. The sky model provides the distribution of shortwave sky radiance. For example, the sky model may use, as inputs: DNI, DHI, and T_dP from weather input; month, day, time; and geographical location (latitude and longitude). Typically, the distribution provided by the sky model will be anisotropic. Advantageously, the sky model obtained at step 128 may be the Perez all-weather sky model. The Perez model has been shown to provide the best overall performance for accurate prediction of global and diffuse solar irradiance and illuminance over a wide range of geographical locations.

Next, at steps 130 to 138, an iterative surface temperature computation operation is performed. The iterations are performed over all surface patches of the surface being modelled. The temperature computation operation comprises, for each patch :

• obtaining the parameters of the patch (step 130); • computing contributions from a plurality of heat transfer processes for the patch (step 132); and

• computing the change in patch surface temperature using the sum of the computed contributions from step 132 (step 134).

The parameters of each surface patch include its associated grid point, its normal vector, and material parameters (such as shortwave albedo, longwave albedo, thickness, emissivity, etc. as discussed above). The amount by which the patch surface temperature is to be updated in a given iteration can be determined according to the following time-dependent energy balance equation:

where c_p is the specific heat capacity (Jkg ^-1K ^-1), p is the material density (kgm^-3), and d is the material thickness (m) of the surface patch of interest.

The manner in which the four (in this case) respective contributions are computed in step 132 will now be described.

Shortwave radiation exchange Q_sw

The shortwave radiation exchange model may be based on Monte Carlo ray- tracing coupled with the Perez all-weather sky model obtained at step 122. Advantageously, the shortwave radiation model described herein provides coverage of the full shortwave spectrum, i.e. the UV, visible and near IR frequencies, compared to previous approaches such as implemented in the ray- tracing program suite Radiance , which mainly focuses on illuminance analysis in the visible spectrum.

In the present shortwave model, the shortwave radiation exchange Q_sw is computed as the summation of absorbed shortwave radiant flux on a surface coming from the Sun I_dir, the sky I_sw,sky and multiple scattering /_SW,scat in the urban scene, given by

where ωⁱ is the weighting coefficient of each frequency band extracted from the solar spectrum (AM 1.5 Global), is the shortwave albedo of the patch, and θ

_S is the angle between the solar incident direction and the patch normal. The superscript i denotes the spectral band, which ranges from 0 to 5, corresponding to the albedo values of the full shortwave spectrum, the UV, blue, green, red, near IR bands, respectively. It should be noted that the value of the weighting coefficient ranges between 0 and 1, where the weighting coefficient of the full shortwave spectrum is given by

. Correspondingly, the full spectrum shortwave radiation exchange always equals to the summed contributions from the five spectral bands is applied to Eq. (1)

for the derivation of patch surface temperature. The solar irradiance in each spectral band is processed in parallel. Different spectral bands share the same ray path, but adopt different weighting coefficients ωⁱ and apply different spectral albedo values

during scattering event on urban surfaces.

To sample diffuse sources, the hemisphere above the grid point is discretised into M-N number of sky patches, such that the projected area of each sky patch on the plane below is identical. This is more computationally efficient as the rays can be given equal weighting during summation. For each grid point associated with a surface patch, sampling rays are sent from the grid point towards the hemisphere, one ray per patch. Advantageously, randomness is introduced by the Monte Carlo Inversion Technique that allows the direction of the sampling ray to vary within the azimuthal θ and zenith Φ boundaries of a sky patch. The radiance value L is summed up to obtain the irradiance I on the grid point, given by

Compared to fixed sky patch sampling, the Monte Carlo hemispherical sampling avoids aliasing, and is more accurate as it more closely mimics photon scattering, which is random in nature. The ray-tracing procedure may make use of the octree structure generated at step 124. When using the octree structure for ray tracing, each ray emanating from a grid point is traced from the grid point to the light source, and the sequence of smallest octree cubes that the ray pierces through is identified. Starting from the octree cube enclosing the grid point, ray-patch interaction tests are carried out until a blockage is identified. Advantageously, the octree structure accelerates the ray tracing process as only the patches contained in the octree cubes in the path of the ray are used for ray-patch intersection tests. The computational load of the most time-consuming component of ray tracing, the ray-patch intersection test, is therefore greatly reduced.

Eq. (3) may be used to compute the contributions of the diffuse sky radiation I_SW,sky and the diffuse scattering in the scene Isw,_Scat, given by

The radiance from each sky patch L_SW,sky is obtained from the Perez all-weather sky model obtained at step 122. The radiance of the i^th frequency band from the scattering surface L_SW,scat,i is derived from the radiant flux on the surface in Eq. (2), · The factor of π is applied as the surfaces are treated as

Lambertian surfaces, i.e. perfectly diffusely reflecting, which is applicable for most urban features, i.e. concrete, pavement, grass etc. Multiple scattering in the scene is accounted for by treating the ray-surface interception point as a secondary grid point and hemispherical sampling can be recursively applied.

The contribution from the direct Sun can be deterministically obtained, as the Sun's direction is predictable from the date, time and geographical location of the scene as provided by the input data at step 110. Longwave radiation exchange Q_LW

In certain embodiments, Monte Carlo ray-tracing (MCRT) may be used to model the longwave radiation exchange, covering the thermal IR spectral band. In the presently proposed longwave model, randomness in both the azimuthal and zenith directions is introduced, making the sky sampling stochastic. The modelling of the longwave contribution is thus less prone to aliasing and more accurate than previous approaches such as implemented in the ray-tracing simulator QESRadiant, which adopts a uniform pattern for sky sampling rotated with a random angle only in the azimuthal direction.

In the presently proposed longwave model, the longwave radiation exchange QLW accounts for absorbed longwave radiant flux on a surface coming from radiant emission of the sky /_LW,sky, surrounding urban surfaces /_LW,surf, scattering in the scene /_LW,scat, and also longwave emitting flux from the surface Q_{LW, emit}, given by

where e is the longwave absorption coefficient, assumed to be the same as surface emissivity. The scattering of the longwave radiation in the scene /_LW,scat is assumed to be negligible.

Eq. (3) is again applied for the computation of /_LW,sky and /_LW,surf, given by

The incoming radiant emittance from the sky is given by

whereby the hourly averaged horizontal infrared radiation from the sky I_IR can be either directly extracted from the ASHRAE International Weather for Energy Calculations (IWEC) weather data representing a typical meteorological year of the chosen location, or computed based on the Stefan Boltzmann's law with the cloudy sky emissivity which combines clear sky emissivity modeled (for example, using the model of Berdahi and Martin, "Emissivity of clear skies," Solar Energy, Vol. 32, 5, 1984) using the dew point temperature near the ground and the cloud cover ratio from the sky clearness index provided by the Perez all-weather sky model. The incoming radiant emittance from surrounding surfaces is given by

The factor of p is applied as the surfaces are treated as Lambertian

surfaces, i.e. perfectly diffusely emitting.

In addition, the outgoing longwave emittance is defined according to the Stefan- Boltzmann law,

Convection Q_conv

Following Newton's law of cooling, the convection heat transfer is given by

, where H is the heat transfer coefficient. Based on an assumption

of well-spaced, low- to mid-rise buidings, the convection coefficient H is assumed to be constant everywhere.

Advantageously, H can be derived prior to starting the simulation method 100, by numerically solving Eq. (1) to obtain a calculated rooftop surface temperature, and fitting the calculated temperature to measured temperatures of the unshaded rooftop at the same times. The surface temperature of the unshaded rooftop is chosen because Eq. (1) can be simplified by simply determining the incoming shortwave and longwave irradiance (W/m²) from the hourly average of the input weather data, i.e. the global horizontal irradiation (Wh/m²) and the horizontal infrared radiation intensity (Wh/m²), respectively, without going through the computationally expensive ray-tracing process. This is because for the rooftop surface, ambient scattering and incoming longwave radiation from surrounding surfaces can be ignored.

Alternatively, H can be derived from the hourly varying wind speed weather input. Any one of a number of models may be used for this purpose, as outlined in Qin, Y. & Hiller, J.E., "Ways of formulating wind speed in heat convection significantly influencing pavement temperature prediction," Heat Mass Transfer, vol. 49, 745, 2013. For more accurate convection modeling, the current model can be replaced by computational fluid dynamics simulation of wind flow, which is highly computationally demanding.

It has been found, though, that the ID convection model proposed herein is capable of providing adequately accurate results in the context of well-spaced, low- to mid-rise buidings, without needing to resort to computationally expensive CFD calculations.

Conduction Q_cond

In certain embodiments, a simplified conduction model may be used to determine the contribution from conduction, as the heat transfer due to conduction is insignificant compared to radiation exchange. Following Fourier's law of heat transfer, the heat transfer by conduction for a building is proportional to the negative temperature gradient from the temperature of the external building wall, T_surf, to that of the internal building wall, T_in, given by , where

d is the thickness of the wall, and k is the thermal conductivity. The conduction model can be replaced with more comprehensive models for more accurate conduction modeling, though as noted, the contribution from conduction is relatively small in the context of a tropical climate such that a conduction model of increased accuracy is unlikely to have a large impact on the simulation results.

Once the respective contributions from the plurality of heat transfer processes have been computed at step 132, the temperature of the surface patch at the end of the time step can be computed according to Eq. (1), at step 134.

At step 136, a check is performed as to whether all patches have been processed.

If not, the method returns to step 130 and the above process (steps 130-134) is repeated for the next patch.

If so, an output operation is performed, at step 138. In particular, the four heat transfer components at the current time step and the surface temperature at the end of the current time step may be output for each patch. At step 139, the surface temperature of each surface patch, and its associated grid point generated at step 120, may be re-initialised to be the computed surface temperature at the end of the current time step.

Next, at step 140, a check is performed as to whether the simulation duration is complete.

If not, the method 100 returns to step 122, but starting again at the first surface patch (not shown) in step 130.

If so, the method 100 ends, and the surface temperature distribution is output by aggregating the computed surface temperature of each surface patch after the plurality of time steps.

The method and system according to certain embodiments has numerous advantages.

Embodiments provide physically accurate microscale modeling of radiation exchange based on ray tracing, accommodating realistic 3D building surface geometry with arbitrary size/orientation, and ability to accommodate different materials for different surface patches. This is especially valuable for architects and developers looking at the design details and orientations of a cluster of new buildings. It should be noted that embodiments of the present invention can be scaled up without losing the level of detail.

Further, the present inventors have realised for modelling of parameters that are important for comfort level in urban environments (such as mean radiant temperature), the contributions of convection and/or conduction may be ignored in many situations. For example, in low-wind conditions which represent a "worst case" scenario for outdoor comfort levels, convection is and radiation exchange is dominated by the shortwave and longwave components of Eq. (1). Accordingly, it is possible to use simplified conduction and convection models that are much less computationally expensive than (for example) CFD-based methods, whilst retaining physical accuracy. For example, a simplified convection model can be obtained by deriving a convection coefficient from fitting of measured and simulated surface temperatures for an unshaded rooftop. By making the ray directions stochastic within each solid angle during ray tracing, embodiments of the invention better mimic the random nature of photon traveling for better accuracy.

Importantly, the presently disclosed approach provides physically accurate radiation exchange modeling with little on-site measurement effort.

Further, the method 100 enables individual heat transfer components to be identified, thereby enabling the identification of the main source of thermal load and facilitating design of targeted heat control measures for enhanced building energy efficiency. Furthermore, it enables the evaluation of building fagade materials during initial design and retrofitting, i.e. ageing due to UV exposure and effectiveness of solar control film.

As will be appreciated, each of the components of Eq. (1) may be replaced with a number of different variants.

For example, instead of Monte Carlo ray tracing, the longwave component QLW may be computed by a method such as fvDOM (G.D. Raithby and E.H. Chui, "A Finite-Volume Method for Predicting a Radiant Heat Transfer in Enclosures With Participating Media", J. Heat Transfer, 112(2), pp. 415-423, 1990). fvDOM is a conservative method to solve the radiative transfer equation for a finite number of discrete solid angles in their corresponding directions.

In another example, the Perez all-weather model may be implemented in OpenFOAM of OpenCFD Ltd (openfoam.com). To fully account for the radiation emitted by the sky in the solver, the skydome is decomposed into 145 evenly- distributed subdivisions for the radiation boundary. The diffuse shortwave radiation and sky longwave radiation can be respectively evaluated based on Perez all-weather and isotropic transposition factor. Turning now to Figure 3, an example architecture of an apparatus 300 for modelling a surface temperature distribution of a built environment is shown. The components of the apparatus 300 can be configured in a variety of ways. The components can be implemented entirely by software to be executed on standard computer server hardware, which may comprise one hardware unit or different computer hardware units distributed over various locations, which may communicate over a network. Some of the components or parts thereof may also be implemented by application specific integrated circuits (ASICs) or field programmable gate arrays.

In the example shown in Figure 3, the apparatus 300 is computing device in the form of a commercially available server computer system based on a 32 bit or a 64 bit Intel architecture, and the processes and/or methods executed or performed by the computing device 300 are implemented in the form of programming instructions of one or more software components or modules 322 stored on non-volatile (e.g., hard disk) computer-readable storage 324 associated with the computing device 300. At least parts of the software modules 322 could alternatively be implemented as one or more dedicated hardware components, such as application-specific integrated circuits (ASICs) and/or field programmable gate arrays (FPGAs).

The computing device 300 includes at least one or more of the following standard, commercially available, computer components, all interconnected by a bus 335:

(a) random access memory (RAM) 326;

(b) at least one computer processor 328, and

(c) at least one network interface connector (NIC) 330 which connects the computer device 300 to a data communications network and/or to external devices.

The computing device 300 includes a plurality of standard software modules, including:

(a) an operating system (OS) 336 (e.g., Linux or Microsoft Windows); and (b) structured query language (SQL) modules 342 (e.g., MySQL, available from http://www.mysql.com), which allow data to be stored in and retrieved/accessed from an SQL database 344.

Advantageously, the database 344 forms part of the computer readable data storage 324. Alternatively, the database 316 is located remote from the computing device 300 shown in Figure 3. The database 344 may store input data 402 (Figure 4), sky model data generated or obtained as part of method 100, and so on.

The boundaries between the modules and components in the software modules 322 are exemplary, and alternative embodiments may merge modules or impose an alternative decomposition of functionality of modules. For example, the modules discussed herein may be decomposed into submodules to be executed as multiple computer processes, and, optionally, on multiple computers. Moreover, alternative embodiments may combine multiple instances of a particular module or submodule. Furthermore, the operations may be combined or the functionality of the operations may be distributed in additional operations in accordance with the invention. Alternatively, such actions may be embodied in the structure of circuitry that implements such functionality, such as the microcode of a complex instruction set computer (CISC), firmware programmed into programmable or erasable/programmable devices, the configuration of a field- programmable gate array (FPGA), the design of a gate array or full-custom application-specific integrated circuit (ASIC), or the like.

Each of the blocks of the flow diagrams of method 100 performed by the apparatus 300 may be executed by a module (of software modules 322) or a portion of a module. The processes may be embodied in a non-transient machine-readable and/or computer-readable medium for configuring a computer system to execute the method. The software modules may be stored within and/or transmitted to a computer system memory to configure the computer system to perform the functions of the module. For example, as shown in Figure 4, the modules 322 may receive input data 402 and may include:

• Grid generator 410, which generates the grid points based on the geometry data;

• Sky model generator 412;

• Octree generator 414, which generates an octree structure for ray tracing based on the geometry data;

• Kernel module 420, which includes a plurality of sub-modules including: o Flow control module 422 for controlling the iterative temperature update operation and providing outputs from said operation; o Convection component 424 for modelling the contribution of convection to the change in surface temperature; o Conduction component 426 for modelling the contribution of conduction to the change in surface temperature; o Shortwave component 428 for modelling the contribution of shortwave radiation exchange to the change in surface temperature, in accordance with the discussion above with reference to Eq. (2)- (5); and o Longwave component 430 for modelling the contribution of longwave radiation exchange to the change in surface temperature, in accordance with the discussion above with reference to Eq. (6)-(8).

The computing device 300 normally processes information according to a program (a list of internally stored instructions such as a particular application program and/or an operating system) and produces resultant output information via input/output (I/O) devices 330. A computer process typically includes an executing (running) program or portion of a program, current program values and state information, and the resources used by the operating system to manage the execution of the process. A parent process may spawn other, child processes to help perform the overall functionality of the parent process. Because the parent process specifically spawns the child processes to perform a portion of the overall functionality of the parent process, the functions performed by child processes (and grandchild processes, etc.) may sometimes be described as being performed by the parent process. Example results obtained by an embodiment of the method 100, as described above, will now be presented. The results were benchmarked against with existing well-established radiation exchange models. In particular, the shortwave component was benchmarked with Radiance software, and the longwave component was compared with fvDOM (Finite Volume Discrete Ordinates Method) from the open-source software OpenFOAM®.

Validation of shortwave heat flux

The case study involved a scene (.STL) consisting of a 3x3 array of building blocks (W3.lmxL12.0mxH2.6m), spaced by a 2.3m gap in the east direction and 2.0m in the north direction, and a ground plane, as shown in Figure 5. The local latitude was 1.37° and longitude was 103.98°, typical of Singapore. The building surface consisted of 142,448 triangular surface patches and the ground consisted of 39,394 triangular surface patches. One grid point was generated at the centroid of each patch. Four grid points were selected on the top, west, east and mid-east surfaces to represent the radiation exchanges in open space, on sparsely spaced building walls and in street canyons.

Firstly, we benchmarked the shortwave radiation module 428 with Radiance. The weather input of hourly averaged direct normal irradiance and diffuse horizontal irradiance was extracted from TMY weather data of Singapore (. EPW ) which is freely accessible from the EnergyPius website. For simplicity, it was assumed that concrete material is used for both buildings and ground with albedo value of 0.21. For Radiance and for the shortwave component 428, 1024 Monte Carlo sampling rays were used and up to two ambient bounces were considered for diffuse scattering of shortwave radiation.

The model performance was evaluated against known observations using three statistical metrics: the coefficient of determination ( R ²) and the refined index of agreement ( RIA ) and the root mean square error (RMSE). R² is a widely used relative measure of fit and it indicates how well the regression line explains the variation in model predictions. It varies between 0 and 1, where 1 indicates the model is in perfect numerical agreement with measurements and 0 indicates no correlation. Unlike R² that is highly sensitive to outliers, RIA focuses more on systematic bias between a model and the known observations. RIA varies between -1 and 1, where 1 indicates a perfect fit and -1 indicates poor prediction. Unlike R² and RIA that are relative measures of fit, RMSE is an absolute measure of fit. It is also highly sensitive to outliers and can be regarded as the standard deviation of the unexplained variance. Lower values of RMSE indicate better fit and 0 means perfect fit. It can be observed from Figure 6 that the absorbed shortwave heat flux results from RESim agree perfectly with Radiance, with R² and RIA values close to 1 and a low RMSE value of 5.4095. The slight difference can be attributed to the stochastic nature of the Monte Carlo sampling technique adopted in both Radiance and the shortwave component 428 of method 100.

Validation of longwave heat flux

Next, the longwave radiation model 430 was benchmarked against fvDOM in OpenFOAM®. The Discrete Ordinate Method (DOM) is regarded as a fairly accurate radiation exchange model. The same geometry was used as shown in Figure 5. To observe the cooling down effect of the buildings, the initial surface temperature of the buildings was set to a high value of 126.75°C (400K) and the ground temperature was set to 26.75°C (300K). The radiation exchange process was carried out over a generic period of one hour, with the iterative time step of 1 minute and the output time interval of 10 minutes. 154 sampling rays were used in both the longwave model 430 and fvDOM. It was assumed that the scene is placed in a black box such that there is no incoming radiation from extraterrestrial space and atmosphere and any outgoing radiation is absorbed by the boundary of the black box. In this case, the surface temperature variation is driven by two modes of radiation exchange, namely the incoming longwave radiation from surrounding surfaces of urban features and the emitted longwave radiation from the surface. The conduction and convection effects were deliberately disregarded.

The energy balance equation was therefore reduced to

We assumed that concrete material is used everywhere, whereby the material properties are listed in Table 1. In addition, we adopted F/=15Wm^-2K ^-1 for air. Table 1 Material properties of concrete.

From Figure 7, the modeled longwave heat flux is in excellent agreement between RESim and fvDOM with R² and RIA values close to 1 and a low RMSE value of 18.73.

Outputs from method 100

This section presents the output options from the method 100.

80 sampling rays were used for both shortwave and longwave components 428 and 430. Up to 1 ambient bounce was considered for shortwave heat flux and ambient scattering for longwave heat flux was ignored. ID conduction and convection models were included for completeness.

In Figure 8, different components of heat transfer on the top, west, east and middle east building surfaces are shown, including the absorbed shortwave heat flux, absorbed longwave heat flux, emitted longwave heat flux, as well as conduction and convection are shown.

In Figure 9, spectral components of radiation exchanges are presented on the top, west, east, and middle east building surfaces, including absorbed UV, visible, near IR, thermal IR, and emitted thermal IR. The breakdown of heat transfer components enables the identification of the main source of thermal load enabling design of targeted heat control measures for enhanced building energy efficiency. Furthermore, it enables the evaluation of building façade materials during initial design and retrofitting, i.e. ageing due to UV exposure and effectiveness of solar control film.

In Figure 10, the 3D surface temperature distribution at three selected timings 1030hrs, 1330hrs and 1630hrs on 21 June are displayed.

Case study

We applied the proposed method 100 to an actual scene in Singapore consisting of two industrial buildings, as shown in Figure 11. For simplicity, a LoD1 (Level of Detail 1) building geometrical model is adopted. Building I is oriented such that the south facing wall is rotated θ=6.3° to the north and the building has a central air-well. It is worth noting that the geometry of building I is fully compatible with the method 100 but could not be handled by the existing ray-tracing radiation transfer solver QESRadiant which is limited to axis aligned geometry. Four sites were selected on the south wall of building I and on the roof, west and south wall of building II for comparison with measured surface temperature. The material properties of the building and ground are tabulated in Table 2. The roof was assigned concrete, building walls were assigned metals, and the ground was assigned asphalt based on site inspection.

The hourly averaged weather inputs, including the direct normal irradiance (DNI), diffuse horizontal irradiance (DHI), and ambient air temperature T_air, were derived from the on-site measured meteorological data on 20 February 2018, as shown in Figure 12. For the purpose of illustration, a constant convection coefficient was used H = 9.84J/m²K to account for wind flow. The simulation was performed at 1- minute intervals and the results were output at 20-minute intervals. 80 Monte Carlo sampling rays were used, and up to 1 ambient bounce was computed for shortwave radiation and ambient bounce was ignored for longwave radiation. From Figure 13, the modeled and measured surface temperature at four sites are in excellent agreement with R² and RIA values close to 1 and a low RMSE value of 2.49.

Table 2 Material properties.

Many modifications will be apparent to those skilled in the art without departing from the scope of the present invention.

Throughout this specification, unless the context requires otherwise, the word "comprise", and variations such as "comprises" and "comprising", will be understood to imply the inclusion of a stated integer or step or group of integers or steps but not the exclusion of any other integer or step or group of integers or steps.

The reference in this specification to any prior publication (or information derived from it), or to any matter which is known, is not, and should not be taken as an acknowledgment or admission or any form of suggestion that that prior publication (or information derived from it) or known matter forms part of the common general knowledge in the field of endeavour to which this specification relates.

Claims

1. A computer-implemented method for modelling a surface temperature distribution of a surface of a built environment, the method comprising, using at least one processor: obtaining input data comprising: a latitude and longitude of the surface; material parameters for one or more materials of the surface; geometrical data indicative of a subdivision of the surface into a plurality of surface patches, each surface patch having a surface normal; and weather data; initialising a temperature of each surface patch; and for each of a plurality of time steps: obtaining, using the weather data, the latitude, and the longitude, a sky model that describes a shortwave sky radiance distribution; for each surface patch, performing a surface temperature computation operation comprising: using the input data and the sky model, computing the temperature of the surface patch for the current time step according to a time-dependent energy balance equation that comprises a sum of separately modelled contributions from a plurality of heat transfer processes; after all patches have been processed, re-initialising the temperature of each surface patch for the next time step to be the computed surface temperature for that surface patch for the current time step; and outputting the surface temperature distribution by aggregating the computed surface temperature of each surface patch after the plurality of time steps.

2. A computer-implemented method according to claim 1, wherein the heat transfer processes are shortwave radiation exchange, longwave radiation exchange, convection, and conduction.

3. A computer-implemented method according to claim 2, wherein the contribution from shortwave radiation exchange and/or the contribution from longwave radiation exchange is, or are, determined based on Monte Carlo ray-tracing.

4. A computer-implemented method according to claim 3, comprising generating an octree structure based on the input geometrical data; wherein said Monte Carlo ray-tracing is performed in accordance with the octree structure, such that for each ray, only surface patches contained in octree cubes in the path of the ray are used for a ray-patch intersection test.

5. A computer-implemented method according to claim 1, wherein each surface patch is characterised by multiband material parameters.

6. A computer-implemented method according to claim 5, wherein the multiband material parameters of at least some respective patches are different than the multiband material parameters of the other respective patches.

7. A computer-implemented method according to claim 2, wherein the contribution from shortwave radiation exchange includes a contribution from multiple diffuse scattering.

8. A computer-implemented method according to claim 2, wherein the sky model is the Perez all-weather sky model, and wherein the clearness index of the Perez all-weather sky model is used to infer cloud cover ratio for determining the contribution from longwave radiation exchange.

9. A computer-implemented method according to claim 2, wherein a convection coefficient for determining the contribution from convection is determined by optimized fitting from unshaded surface temperature measurements, or by derivation based on wind speed from the weather data.

10. An apparatus for modelling a surface temperature distribution of a surface of a built environment, the apparatus comprising: a kernel module that is configured to: obtain input data comprising: a latitude and longitude of the surface; material parameters for one or more materials of the surface; geometrical data indicative of a subdivision of the surface into a plurality of surface patches, each surface patch having a surface normal; and weather data; initialise a temperature of each surface patch; and for each of a plurality of time steps: obtain, by a sky model component, using the weather data, the latitude, and the longitude, a sky model that describes a shortwave sky radiance distribution; for each surface patch, perform a surface temperature computation operation comprising: using the input data and the sky model, computing the temperature of the surface patch for the current time step according to a time-dependent energy balance equation that comprises a sum of separately modelled contributions from a plurality of heat transfer processes; after all patches have been processed, re-initialising the temperature of each surface patch for the next time step to be the computed surface temperature for that surface patch for the current time step; and output the surface temperature distribution by aggregating the computed surface temperature of each surface patch after the plurality of time steps.

11. An apparatus according to claim 10, wherein the heat transfer processes are shortwave radiation exchange, longwave radiation exchange, convection, and conduction.

12. An apparatus according to claim 11, comprising a shortwave component that is configured to determine the contribution from shortwave radiation exchange based on Monte Carlo ray-tracing; and/or a longwave component that is configured to determine the contribution from longwave radiation exchange based on Monte Carlo ray-tracing.

13. An apparatus according to claim 12, comprising an octree generator that is configured to generate an octree structure based on the input geometrical data; wherein the shortwave component and/or the longwave component is configured to perform said Monte Carlo ray-tracing in accordance with the octree structure, such that for each ray, only surface patches contained in octree cubes in the path of the ray are used for a ray-patch intersection test.

14. An apparatus according to claim 10, wherein each surface patch is characterised by multiband material parameters.

15. An apparatus according to claim 14, wherein the multiband material parameters of at least some respective patches are different than the multiband material parameters of the other respective patches.

16. An apparatus according to claim 11, wherein the contribution from shortwave radiation exchange includes a contribution from multiple diffuse scattering.

17. An apparatus according to claim 11, wherein the sky model is the Perez all- weather sky model, and wherein the sky model component is configured to use the clearness index of the Perez all-weather sky model to infer cloud cover ratio for determining the contribution from longwave radiation exchange.

18. An apparatus according to claim 11, comprising a convection component that is configured to determine a convection coefficient by optimized fitting from unshaded surface temperature measurements, or by derivation based on wind speed from the weather data.

19. A system for modelling a surface temperature distribution of a surface of a built environment, comprising at least one processor in communication with computer-readable storage, the computer-readable storage having stored thereon instructions for causing the at least one processor to perform a method according to claim 1.