An open-source solver for coupled physical processes modeling urban microclimate based on OpenFOAM

Dr. Clément Nevers, Email: [email protected]
Dr. Aytaç Kubilay, Email: [email protected]

Contributors:
Dr. Andrea Ferrari
Dr. Lento Manickathan

Table of Contents:

• A. Introduction
• B. Models and governing equations
• C. Implementation in urbanMicroclimateFoam
• D. Simulation settings
• E. Running urbanMicroclimateFoam simulations
• F. Additional notes
• References

urbanMicroclimateFoam is a multi-region solver consisting of an "air" region, representing the subdomain for turbulent transport, and multiple "solid" regions, representing subdomains for built environment, e.g. facades, roofs, street, etc. It is originally based on chtMultiRegionFoam, which is a multi-region solver included in the standard distribution of OpenFOAM, with additional functionality.

urbanMicroclimateFoam includes:

• CFD: A computational fluid dynamics model solving the turbulent, convective air flow and heat and moisture transport in the air subdomain.

• HAM: A coupled heat and moisture transport model solving the absorption, transport and storage of heat and moisture in the porous building materials.

• RAD: A radiation model solving the net longwave and shortwave radiative heat fluxes using a view factor approach.

• VEG: A vegetation model solving heat balance for urban trees and green surfaces.

The source code for the solver urbanMicroclimateFoam, and several examples cases, can be accessed at the following address: https://github.com/OpenFOAM-BuildingPhysics/urbanMicroclimateFoam

The following documentation assumes that the user has a working OpenFOAM installation and a basic understanding of how it works. For a list of alternative ways to install OpenFOAM, as well as the official user guide, please refer to: https://www.openfoam.org/

In order to use urbanMicroclimateFoam, one must go into the downloaded folder that includes the source code:

cd /path/to/urbanMicroclimateFoam

and compile it together with the necessary libraries by typing:

./Allwmake

In the air region, the solver uses steady Reynolds-averaged Navier-Stokes (RANS) to calculate wind-flow field. The general structures of the governing equations for continuity, momentum and energy are similar to the standard OpenFOAM solvers buoyantSimpleFoam and chtMultiRegionFoam. Additionally, the turbulent transport of humidity is solved.

Continuity and momentum equations are as follows, respectively:

$$ \nabla \cdot (\rho u) = 0~~~~~~~~~~(1) $$

$$ \rho (u \cdot \nabla) u = - \nabla p + \nabla \cdot (\tau + \tau_t) + \rho g + s_u~~~~~~~~~~(2) $$

where $\rho$ is air density, $u$ velocity vector, $p$ pressure, $\tau$ and $\tau_t$ viscous and turbulent stress tensors, respectively, $g$ gravitational acceleration vector and $s_u$ momentum sources.

The transport of energy and humidity in the air are as follows, respectively:

$$ \rho (u \cdot \nabla) h + \rho (u \cdot \nabla) K = - \nabla \cdot (q + q_t) + s_h~~~~~~~~~~(3) $$

$$ \rho (u \cdot \nabla) w = - \nabla \cdot (g + g_t) + s_w~~~~~~~~~~(4) $$

where $h$ is enthalpy, $K$ specific kinetic energy, $q$ and $q_t$ laminar and turbulent heat fluxes, respectively, $s_h$ energy sources, $w$ humidity ratio, $g$ and $g_t$ laminar and turbulent humidity fluxes, respectively, and $s_w$ moisture sources.

In solid regions, the solver uses a continuum modeling approach for the coupled transport of heat and moisture (Whitaker, 1977). The coupled heat and moisture transport equations are given in equations (1) and (2), respectively (Janssen et al., 2007):

$$ (c_0 \rho_0 + c_l w_s) \frac{\partial T}{\partial t} = - \nabla (q_c + q_a)~~~~~~~~~~(5) $$

$$ \frac{\partial w_s}{\partial p_c} \frac{\partial p_c}{\partial t} = - \nabla (g_l + g_v)~~~~~~~~~~(6) $$

where $c_0$ denotes the specific heat of dry material, $\rho_0$ the density of dry material, $c_l$ the specific heat of liquid water, $w_s$ the moisture content, $T$ absolute temperature and $p_c$ capillary pressure. The derivative $\frac{\partial w_s}{\partial p_c}$ represents the moisture capacity of the porous material. Here, the capillary pressure is defined as the pressure difference between the liquid and the gaseous phase and, therefore, negative for unsaturated conditions.

$q_c$ and $q_a$ denote the conductive and advective heat fluxes, respectively. $g_l$ and $g_v$ denote liquid and vapor moisture fluxes, respectively.

It is possible to model part of the urban surfaces as impermeable, or to turn off moisture transport in porous urban media completely. In such case, equation (1) for the heat transport simplifies to pure heat conduction:

$$ (c_0 \rho_0) \frac{\partial T}{\partial t} = - \nabla (q_c)~~~~~~~~~~(7) $$

Implementation of the HAM model in solid regions and the material library for moisture storage and transport properties are identical to the solver hamFoam. Further information is available in hamFoam wiki.

Long-wave (thermal) and short-wave (solar) radiative exchanges are calculated with view factors using a radiosity approach. The implementation considers multiple reflections of both long-wave and short-wave radiation.

Direct solar irradiation is modeled with ray casting to capture the variation of shading and incidence angle over the day. Diffuse solar radiation is modeled based on the surface orientation.

The model assumes that surfaces are gray and diffuse, i.e. radiative properties such as emissivity and albedo are independent of wavelength and direction. All surfaces are opaque to both long-wave and short-wave radiation, i.e. transmissivity is zero and, therefore, windows on building facades are not considered. Further, due to the nature of view factor approach, absorption, scattering and emission of radiation by air are neglected.

Trees are modeled as porous fluid zones by setting non-zero values for a field variable defining leaf area density (LAD).

Heat balance on a leaf is defined based on sensible, radiative and latent heat fluxes. Transpiration rate is a function of stomatal resistance, which takes into account the physiological response of trees based on species, solar radiation, vapor pressure deficit (relative humidity) in the air and available soil moisture content. Once leaf temperatures are obtained with an iterative energy balance calculation, the resulting impact of cooling can be estimated.

Source/sink terms are defined within these porous zones for momentum, heat, moisture.

Coupling between the air and solid regions happens at the coupled boundaries of mappedWall type. Required data are sampled from the samplePatch of the sampleRegion. These inputs are defined in the boundary files of the mesh for the air and solid regions, following the available sampleMode's in OpenFOAM, e.g. as also in chtMultiRegionFoam.

• Data exchange between coupled patches occurs at "exchange timesteps", i.e. once every deltaT as defined in file controlDict. Between CFD and HAM, converged data from one model is used for the next time step of the other model.

• CFD uses Dirichlet conditions at the coupled boundaries, i.e. fixed value temperature and humidity, calculated by HAM. HAM uses Neumann conditions at the coupled boundaries, i.e. heat and moisture fluxes, based on the transfer coefficients calculated by CFD.

• CFD solves for the steady air flow, while HAM solves for the unsteady transport in order to take into account dynamic storage. This is based on the assumption that timescales for air transport in CFD are much shorter, i.e. air flow adapts to changes relatively faster compared to the response rate in the HAM domain.

The following part summarizes the directory structure and input files in a typical case folder for urbanMicroclimateFoam:

• 0/air/: Initial time directory with initial conditions and boundary conditions for the field variables defined for the air region.

also cloudCover, Tambient, LAD.

• 0/"solidRegion"/: Initial time directory with initial conditions and boundary conditions for the solid region named "solidRegion". There can be multiple of these.

• constant/: Includes the mesh files for each region and input files for simulation settings.

• constant/air/: Input files for simulation options specific to the air region.

• constant/"solidRegion"/: Input files for simulation options specific to the solid region named "solidRegion".

• system/: Input files for solver settings.

• system/air/: Input files for solver settings specific to the air region.

• system/"solidRegion"/: Input files for solver settings specific to the solid region named "solidRegion".

The solver is compatible with the standard boundary conditions available in OpenFOAM, and also other third party boundary conditions, e.g. groovyBC. This part describes the boundary conditions, specific to urbanMicroclimateFoam, which provide additional functionality.

• Coupled boundary conditions:

  • CFDHAMfluidMoistureCoupledMixed

    Humidity ratio boundary condition for fluid side of a coupled boundary. Obtains fixed value humidity ratio from solid side.

    Example usage:

    patchName
    {
        type            compressible::CFDHAMfluidMoistureCoupledMixed;
        value           uniform 0.008;
    }

    Value entry is only an initial value and will be overwritten by the boundary condition.

  • CFDHAMfluidMoistureCoupledImpermeable

    Humidity ratio boundary condition for fluid side of a coupled boundary. Imposes zero gradient for an impermeable wall. Used together with CFDHAMsolidMoistureCoupledImpermeable.

    Example usage:

    patchName
    {
        type            compressible::CFDHAMfluidMoistureCoupledImpermeable;
        value           uniform 0.008;
    }

    Value entry is only an initial value and will be overwritten by the boundary condition.

  • CFDHAMfluidTemperatureCoupledMixed

    Temperature boundary condition for fluid side of a coupled boundary. Obtains fixed value temperature from solid side.

    Example usage:

     patchName
     {
         type            compressible::CFDHAMfluidTemperatureCoupledMixed;
         value           uniform 300.0;
     }

    Value entry is only an initial value and will be overwritten by the boundary condition.

  • CFDHAMsolidMoistureCoupledMixed

    Moisture content boundary condition for solid side of a coupled boundary. Obtains moisture flux from fluid side.

    Example usage:

     patchName
     {
         type            compressible::CFDHAMsolidMoistureCoupledMixed;
         value           uniform -1e7;
     }

    Value entry is only an initial value and will be overwritten by the boundary condition.

  • CFDHAMsolidMoistureCoupledImpermeable

    Moisture content boundary condition for solid side of a coupled boundary. Imposes zero gradient for an impermeable wall. Used together with CFDHAMfluidMoistureCoupledImpermeable.

    Example usage:

     patchName
     {
         type            compressible::CFDHAMsolidMoistureCoupledImpermeable;
         value           uniform -1e7;
     }

    Value entry is only an initial value and will be overwritten by the boundary condition.

  • CFDHAMsolidTemperatureCoupledMixed

    Temperature boundary condition for solid side of a coupled boundary. Obtains heat flux from fluid side.

    Example usage:

     patchName
     {
         type            compressible::CFDHAMsolidTemperatureCoupledMixed;
         qrNbr           qr;
         qsNbr           qs;
         value           uniform -1e7;
     }

    Value entry is only an initial value and will be overwritten by the boundary condition.

  • CFDHAMsolidMoistureTransferCoeff

    Moisture content boundary condition for solid side of a coupled boundary. Calculates moisture flux from the provided transfer coefficient, betacoeff, and exterior vapor pressure, pv_o.

    Intended for use without coupling with CFD, e.g. during initialization of moisture content field before the actual coupled simulation.

    Example usage:

     patchName
     {
         type            compressible::CFDHAMsolidMoistureTransferCoeff;
         betacoeff       1e-8;
         pv_o            "$FOAM_CASE/0/pv_o.dat";
         value           uniform -1e7;
     }

    betacoeff is a constant. pv_o is a list, varying over time.

    Value entry is only an initial value and will be overwritten by the boundary condition.

  • CFDHAMsolidTemperatureHeatCoeff

    Temperature boundary condition for solid side of a coupled boundary. Calculates heat flux from the provided transfer coefficients, hcoeff and betacoeff, and exterior temperature and vapor pressure, Tamb and pv_o.

    Intended for use without coupling with CFD, e.g. during initialization of temperature field before the actual coupled simulation.

    Example usage:

     patchName
     myInterfacePatchName
     {
         type            compressible::CFDHAMsolidTemperatureTransferCoeff;
         qrNbr           qr;
         qsNbr           qs;
         hcoeff          8.0;
         Tamb            "$FOAM_CASE/0/Ta.dat";
         betacoeff       1e-8;
         pv_o            "$FOAM_CASE/0/pv_o.dat";            
         value           uniform 293.15;
     }

    hcoeff abd betacoeff are constants. Tamb and pv_o are lists, varying over time.

    Value entry is only an initial value and will be overwritten by the boundary condition.

• Boundary conditions related to blending layer:

  • blendingLayerVel

    Fixed value boundary condition for velocity at the lateral boundaries of CFD domain. Intended for use when imposing boundary values from available results of a mesoscale simulation.

    Velocity is read from input files that contain lists of boundary values. Input files are located in folder constant/air/Utarget_"patchName"/. Names of input files are Utarget_"patchName"_"timeValue", e.g. Utarget_east_10800.

    Example usage:

     patchName
     {
         type            blendingLayerVel;
         inputTimeStep   3600;
         value           uniform (1 1 0);
     }

    inputTimeStep is the interval of input files, e.g. 3600 for hourly input from mesoscale simulations. In case simulation deltaT is smaller, boundary values will be linearly interpolated between two input files.

    Value entry is only an initial value and will be overwritten by the boundary condition.

  • blendingLayerTemp

    Fixed value boundary condition for temperature at the lateral boundaries of CFD domain. Intended for use when imposing boundary values from available results of a mesoscale simulation.

    Temperature is read from input files that contain lists of boundary values. Input files are located in folder constant/air/Ttarget_"patchName"/. Names of input files are Ttarget_"patchName"_"timeValue", e.g. Ttarget_north_36000.

    Example usage:

     patchName
     {
         type            blendingLayerTemp;
         inputTimeStep   3600;            
         value           uniform 300.0;
     }

    inputTimeStep is the interval of input files, e.g. 3600 for hourly input from mesoscale simulations. In case simulation deltaT is smaller, boundary values will be linearly interpolated between two input files.

    Value entry is only an initial value and will be overwritten by the boundary condition.

Solver settings specific to urbanMicroclimateFoam are defined in file system/controlDict:

• time step and output frequency:

  • deltaT: "exchange timestep" in seconds, as described in part C.1. CFD and HAM will be solved one after the other for the duration of deltaT. Typically, deltaT is between 600-3600 s. Shorter deltaT will increase the frequency of data exchange between CFD and HAM solutions, increasing accuracy and, also, computational cost.

  • writeInterval: This should be kept as 1, so that the solver writes output at each deltaT.

• settings for fluid solution:

  • minFluidIteration: minimum number of internal iterations during steady CFD solution for each "exchange timestep". The solver will continue air flow calculation until it reaches the number of iterations given here even if the SIMPLE convergence criteria in file *system/air/fvSolution" is reached.

  • maxFluidIteration: maximum number of internal iterations during steady CFD solution for each "exchange timestep". The solver will stop air flow calculation once it reaches the number of iterations here and continue with HAM calculation,

• settings for solid solution:

  • initialSolidTimestepFactor: this factor, multiplied with the value given for solver deltaT, defines the first timestep in the unsteady HAM calculation. The solver uses adaptive timesteps based on the Picard iteration scheme (Celia et al., 1990). For further information on the Picard iteration, refer to hamFoam wiki.

  • minDeltaT: the minimum timestep before the solver fails. In case the Picard iteration fails after reaching this timestep, the timestep is not decreased further. Instead, the solver stops and the solver writes the field variables at that time step for inspection.

  • maxDeltaT: the upper limit for the timesteps.

• settings for blending layer (in case blending is enabled):

  • dampingThickness: the distance from the lateral boundaries of CFD domain, where the blending layer is activated.

  • alphaCoeffU: damping strength for velocity

  • alphaCoeffT: damping strength for temperature

Material properties for the porous solid urban materials are provided in constant/"solidRegionName"/transportProperties. urbanMicroclimateFoam allows multiple materials to be defined within a solid region.

An example setup with two materials defined for building facade is as follows:

buildingMaterials
(
	{
	    name       loadBearing;
	    buildingMaterialModel    HamstadBrick;
	    rho        2005;
	    cap        840;
	    lambda1    0.5;
	    lambda2    4.5e-3;
	}
	{
	    name        finishing;
	    buildingMaterialModel    HamstadPlaster;
	    rho        790;
	    cap        870;
	    lambda1    0.2;
	    lambda2    4.5e-3;
	}
);

Different materials are handled by different cellZones. The name entry is the name of the cellZone, which includes the cells that are assigned for the given material. The buildingMaterialModel entry is the material name from the material library, which defines moisture transport properties such as moisture retention, liquid permeability and vapor permeability. The entries rho and cap are the density (rho [kg/m³]) and the specific heat (cap [J/kgK]) of the dry material, respectively. Thermal conductivity [W/mK] is defined as:

$$ \lambda = \lambda_1 + \lambda_2 \cdot w_s~~~~~~~~~~(8) $$

where $\lambda_1$ is the thermal conductivity of the dry material and $\lambda_2$ is the coefficient for the thermal conductivity in moist state.

For a complete list of materials, refer to the folder _LIB/buildingMaterialModel in the source code.

Settings for long-wave and short-wave radiation are located in case files constant/air/radiationProperties and constant/air/solarLoadProperties, respectively.

Due to the current implementation of radiation modeling, these settings are less flexible compared to the other submodels in the solver and may not work other available models in OpenFOAM.

It is recommended that solverFreq is set to 1, so that the values are updated at each "exchange timestep".

Radiation models can be turned off by setting the model name to none, e.g.:

radiationModel  none;

or by removing the file for model properties from case folder.

To begin the simulation, first, go into the case folder:

cd /path/to/caseFolder

For mesh generation, please refer to various examples provided in the tutorials. Steps for mesh generation can be followed in files named as Allprepare.

There are several pre-processing utilities that are required to run before the simulation starts:

faceAgglomerate -region air
viewFactorsGen -region air
solarRayTracingGen -region air

Once these the simulation can be started (in serial calculation mode) by executing:

urbanMicroclimateFoam

Note that, for simplicity, this manual is written for serial calculation. In parallel run, all commands must use the following way, including the pre-processing tools:

mpirun -np N faceAgglomerate -region air -parallel >log.faceAgglomerate
mpirun -np N viewFactorsGen -region air -parallel >log.viewFactorsGen
mpirun -np N solarRayTracingGen -region air -parallel >log.solarRayTracingGen
mpirun -np N urbanMicroclimateFoam -parallel >log

where N is the number of processors.

For further information on parallel computing with OpenFOAM, please refer to the official OpenFOAM documentation.

OpenFOAM provides several approaches to create cellZones. Some examples are given below.

• topoSet

topoSet utility creates or modifies cellSets/faceSets/pointSets by using a variety of sources for selecting cells/faces/points. These sources include boxes, defined with two opposite corners, and surfaces, defined by an .stl file.

An input file topoSetDict is available within OpenFOAM source code as an example and can be copied to your simulation case folder by this command:

foamGet topoSetDict

which may offer multiple files as option. The "annotated" one, located in OpenFOAM folder OpenFOAM/OpenFOAM-6/etc/caseDicts/annotated/topoSetDict provides an extensive list of usage examples and available sources for selection.

You can find below an example of using boxToCell to select cells within a box of rectangular prism. First, a cellSet named facadeCells is created. Then, this cellSet is converted to a cellZone with the same name. This name should be identical to the name provided in transportProperties as shown in part D.2.

actions
(
    {
        name    facadeCells;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box   (70.0 0 100.0)(70.1 10.0 150.0);
        }
    }
    {
        name    facadeCells;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set    facadeCells;
        }
    }           
);

Another example below shows the usage of surfaceToCell to select cells within a distance of an .stl file. This example creates two cellZones named pavement and soil.

actions
(
    // create a cellSet named pavement within a distance of 10 cm to ground surface
    {
        name    pavement;
        type    cellSet;
        action  new;
        source surfaceToCell;
        sourceInfo
        {
            file            "constant/ground/triSurface/ground.stl";
            outsidePoints   ((505.1 345.5 619));
            includeCut      false;
            includeInside   false;
            includeOutside  false;
            nearDistance    0.1;
            curvature       -100;
        }
    }
    // create a cellZone named pavement from the cellSet pavement
    {
        name    pavement;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set    pavement;
        }
    }  
    
    // create a new cellSet named soil, identical to the cellSet pavement
    {
        name    soil;
        type    cellSet;
        action  new;
        source  cellToCell;
        sourceInfo
        {
            set pavement;
        }
    }
    // invert the selection of cells in the cellSet soil. 
    // this will select all the remaining cells in the solid region mesh
    {
        name    soil;
        type    cellSet;
        action  invert;
    } 
    // create a cellZone named soil from the cellSet soil
    {
        name    soil;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set    soil;
        }
    }              
           
);

After placing the input file topoSetDict under the relevant region folder system/"regionName"/, topoSet can be run with the following command:

topoSet -region regionName

• setSet

setSet is another utility that can create or modifies cellSets. It can run interactively or by the use of batch files.

An example input in batch files can be as follows, which crates a cellZone named facadeCells:

cellSet facadeCells new boxToCell (70.0 0 100.0)(70.1 10.0 150.0)
cellZoneSet facadeCells new setToCellZone facadeCells

After placing the batch file with the content above under the relevant region folder system/"regionName"/, setSet can be run with the following command:

setSet -region regionName -batch system/regionName/setset.batch

urbanMicroclimateFoam will solve each steady CFD calculation for an "exchange timestep" until the convergence criteria defined by residualControl settings in case file system/air/fvSolution are reached or the maximum number of internal iterations defined by maxFluidIteration in case file system/controlDict is reached.

The solver writes only the final output of each steady CFD calculation and not the intermediate states of flow fields between each "exchange timestep", for example as in a calculation of solver buoyantSimpleFoam. In order to monitor the convergence of field variables at a number of monitoring points, function objects available in OpenFOAM can be used. The following entry in system/controlDict will write the values at the given list of probe locations under case folder postProcessing/probes/air/.

functions
{
    probes
    {   
        // Where to load it from
        libs            ("libsampling.so");

        type            probes;

        // Name of the directory for probe data
        name            probes;

        // Write at same frequency as fields
        writeControl    timeStep;
        writeInterval   10;  

        // Fields to be probed
        fields          (T U); 
        region air;

        probeLocations
        (
            (100.595 86.3716 2.0)
            (119.398 86.1579 2.0)
        );
    }    
}

The example above writes the values of air temperature and velocity at two probe locations every 10th internal interation.

Celia, M.A., Bouloutas, E.T., Zarba, R.L., 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26, 1483–1496.

Janssen, H., Blocken, B., Carmeliet, J., 2007. Conservative modelling of the moisture and heat transfer in building components under atmospheric excitation. Int. J. Heat Mass Transf. 50, 1128–1140.

Whitaker, S., 1977. Simultaneous Heat, Mass, and Momentum Transfer in Porous Media: A Theory of Drying. Adv. Heat Transf. 13, 119–203.