Efficient hygrothermal wall transport models inside simulations of conditioned buildings

With regard to room comfort and potential moisture damage hygrothermal conditions inside a building are of high interest, especially in highly humid or dry climatic regions. In order to achieve satisfying climatic conditions not only thermal storage, but also moisture storage of walls needs to be taken into consideration. Dynamic effects of walls and the impact on thermal room climate require a detailed thermal wall simulation. In the simulation program NANDRAD a Finite Volume approach is used to discretized each building construction. The balance equation system is now extended with coupled moisture mass balance, both for zone air volumes and all construction elements. This leads to a duplication of system unknowns with highly nonlinear interactions between solution variables, resulting in strong increase of simulation time. In order to improve performance for larger simulation models, a mixed grid approach is evaluated. Hereby, a detailed grid is used for energy balances and coarser grid for moisture balances. In the article we discuss the procedure, evaluate the impact on accuracy of solution results, and discuss numerical problems. We show that the procedure indeed helps reducing the overall number of unknowns, which, however, appears to pay off only for large building models.


Introduction
Whole building energy simulation primarily involves analysis of energy balances and thermal comfort. However, there are applications where humidity conditions in building zones are also of interest. This is particularly true for hot and humid climatic zones, where perfect air conditioning may not be possible or desired in order to reduce energy consumption. Further, evaluation of moisture buffering effects to improve indoor air quality and/or prevention of moisture related damage may be motivations to include moisture balances in whole building models.
Modelling coupled heat and moisture transport in porous building construction is state of the art since many years and simulation models such as DELPHIN [9], WUFI [8], etc. can be used to analyse individual wall constructions in a very detailed manner. With respect to whole building simulation programs, very few use discretized construction models with detailed calculation of thermal storage capacity. NANDRAD is such a whole-building energy simulation engine and uses a traditional multi-zone approach with balance equations for all air spaces (zones). All constructions (walls, ceilings, floors, roof) are modelled as one-dimensional multi-layered detailed wall constructions (see Fig. 1). WUFI Plus [3] uses a similar approach.
After application of the Finite Volume Method for the partial differential energy balance equations of the construction, a typically large system of ordinary differential equations is obtained. These stiffly coupled equations are integrated with an error controlled, variable time step integrator CVODE [7]. The embedded nonlinear equation system is solved with a modified Newton method, hereby exploiting the sparse nature of the equation system through use of either sparse direct solvers [5] or preconditioned Krylov subspace methods [10].
The NANDRAD model was now extended such, that each energy balance equation is complemented by a moisture balance equation. Conserved quantity is the moisture mass density, corresponding to absolute humidity in zonal mass balances. The new moisture mass balance equations are tightly coupled to their respective energy balance equations.
Formulation of the moisture storage and transport equations is based on DELPHIN model, yet somewhat simplified because application in building simulation typically excludes excessive wetting and drying Note, that in NANDRAD the entire system of equations is solved in a fully coupled manner, thus ensuring full conservation of energy and moisture masses and time-accurate integration. Other simulation models, for example WUFI Plus [3] use a co-simulation type calculation algorithm to interface zonal balances with wall construction sub-models.

Grid variation study
It can be expected that reducing the number of unknowns will improve performance again. Therefore, in a first attempt a grid variation study was performed. The test case is based on a residential building model. A zone with a single outside wall (brick with plaster coating), and several indoor wall constructions and a floor construction is simulated. Ambient climate with typical boundary conditions is used, and a constant ventilation rate of 0.5 1/h is applied. On weekday from 8 a.m. till 5 p.m. the room is considered as occupied. The zone is heated ideally with a constant setpoint of 20°C. We apply heat sources to room air (maximum 20 person, heat load 80 W per person) and additional moisture sources (150 g/h). For non-occupied day time heating setpoint is reduced to 16°C and person heat sources as well as moisture sources are deactivated. This test case allows analysis of grid related effects as studied in this article.
In the grid sensitivity study, the overall number of Finite Volume cells in all constructions was reduced (through larger stretch factors in the grid generation algorithms). For the coarse grid, this resulted in an inacceptable damping of zone relative humidity compared to the detailed/fine grid variant (Fig. 2). This effect is caused by wall moisture storage, which is clearly overestimated in this case. Much better results are achieved by using a medium discretization.   2 shows three grid sensitivity variants and lists the grid cell count used for the outside wall construction. The very coarse grid algorithm uses each wall layer as single grid cell and simply splits the boundary elements into two (in order to allow surface value extrapolation). Thus, a 3layered construction has 5 cells. For the other grid variants, a tanh double-sided stretch algorithm is used, with different grid density parameter. The coarse grid variant can be roughly compared to the surface and deep moisture buffer layer approach in TRNSYS [6]. Clearly, the accuracy of this approach is limited.
As illustration, the profiles computed inside the outside wall construction are shown in Fig. 3. At the first view, application of a medium size discretization can be considered as sufficient. Indeed, even for this non-critical test case we observe significant deviations from detailed solution. Note, that there exist application cases with demand on a fine thermal wall discretization. For instance, NANDRAD allows modelling of internal thermal wall sources [11]. This feature supports advanced simulation cases such as simulation of underfloor heating systems with regard to heat delay and phase change materials.
For this reason, the authors do not advice application of a medium grid discretization without a preceding sensitivity study per default. Instead, we introduce a new approach: We focus on hygric response of room air climate and the damping effect of moisture storage inside the wall as shown in the presented example. Assuming the related moisture wetting and drying processes to be rather slow compared to temperature changes in wall constructions in these applications, it could be useful to reduce the number of Finite Volume cells for the moisture balance only, and keep the detailed grid for the energy balance. We call this a mixed-grid Finite Volume method. There are several research questions associated with this approach: 1. How can a mixed grid calculation be implemented with regard to a good accuracy? 2. Is it worth the extra effort, i.e. do we see a significant performance gain? 3. Do we work out a general rule both guaranteeing accuracy and high simulation performance for any application case?
We find it useful to use an exemplary approach for mixed grid simulation, including a simple coarse grid generation rule for moisture balance. In a first step we illustrate some implementation details and check validity of simulation results for our test case. The next step includes performance analysis. If our procedure shows significant impact on simulation performance, we can analyze the other questions in order to obtain a generalpurpose grid selection and flux calculation algorithm.

Mixed grid approach
As already mentioned, we use a fine discretization for energy balance equations in order to keep an accurate representation of the temperature distribution and storage capacity of the construction. The grid density required for this is determined in discretization study using thermal balance equations, only.
Previous studies related to suitable wall grids when using energy balance alone showed that the surface elements in exchange with the room air need to be fairly thin (1 to max. 5 mm). Otherwise, the room air balance has access to too much thermal capacity and calculated room temperature oscillations are dampened quite substantially. With respect to moisture buffering behavior of wall constructions, the surface elements must be equally thin, in order to avoid excessive dampening of moisture oscillations in the room air. Consequently, we can only remove/enlarge grid cells inside the construction, but should keep a fine grid near the surface.
For the moisture balance inside the wall, we use a rather coarse numerical grid (see Fig. 4) and solve for moisture mass density. The coarse and fine grids are aligned such, that the grid cells share common faces for heat and moisture transport. Basically, we just merge the fine grid cells from the thermal balances inside the construction. To distinguish grid cell numbering, we use cell indexes i for the energy balance grid, and j for the moisture balance grid. Applying the standard Finite Volume method for either grid and balance equation, we need to formulate heat flows across cell boundaries for the fine grid, and moisture flows (specifically vapor diffusion fluxes) across the coarse grid cell boundaries. We use a standard formulation of heat and vapor diffusion flux quantities [1,6,9].
When evaluating these fluxes across Finite Volume cell boundaries, we need to average transport coefficients and determine gradients. With respect to thermal transport, the resistance-based flux computation has proven to be best (corresponding to harmonic averaging of thermal conductivities). Equation (1) In this study, we restrict ourselves to consider only vapor transport inside the construction. For the flux calculation we use vapor pressure gradient as driving potential. There are different ways to approximate a pressure gradient across the cell interface. In our study, we chose a classical FV approach, by evaluating the vapor pressures in each cell the coarse grid and use finite-difference quotients as approximation for the respective gradient (Equation (2)). The vapor conductivity is a distance-weighted arithmetic average of the cell-specific vapor conductivities (harmonic averaging would also be possible).
The vapor pressures in the coarse grid cells are computed with an area-weighted average of the temperatures in the corresponding set of fine grid cells. The enthalpy fluxes across the cell interfaces are computed using upwinding,

Implementation
In the first, traditional Finite Volume implementation, the solution variables for the balance equations of energy and moisture were stored one after another in the global solution vector. The composition of the Jacobian matrix was easily obtained, since instead of scalar variables in the thermal-only mode, we could simply use 2x2 block matrices in the same locations. However, with the mixedgrid approach, the composition of the global equation system is much more complex. Indeed, knowledge of the Jacobian pattern is essential for usage of sparse-matrix solver techniques -either sparse direct solvers such as the KLU, or preconditioned Krylov subspace methods. The preconditioners used in the latter require detailed dependency information of variables and equations. In our opinion, exploiting the sparse nature of the problem with such techniques is mandatory to achieve sufficient simulation performance.
A major task in implementing the mixed grid approach was hence the composition of the new global solution vector and the corresponding Jacobian pattern. We found it practical to sort the solution variables such, that first the conserved quantities of the energy balance are stored and afterwards all moisture states. With the knowledge of the grid mapping (e.g. which fine grid cells were merged into coarse grid cells), the Jacobian pattern could be defined. The Jacobian matrix itself is computed using finite difference quotient approximation.

Simulation results
We apply our mixed grid approach to the simulation example of the previous section. For the mixed grid, we need to generate a grid refinement as seen in Fig. 4 avoiding overlapping volumes of fine and coarse grid. For this purpose, we split up a thin layer of 5 mm from the boundary material layer inside the construction of each wall.
The resulting temperature and relative humidity profiles are shown in Fig. 6. Note, that relative humidity is calculated from moisture mass density with the help of sorbtion isotherm. Therefore, relative humidity of the mixed grid is directly correlated to moisture content and kept constant over the thickness of each coarse grid layer and the small boundary layers. Relative humidity value at wall boundaries is very well approximated by this approach. We assume, that in the current test case dissipation of vapor towards outside/the inside room is approximated very well. Further, temperature profile of mixed and fine grid approach are almost identical.
As a result, zone climate in Fig. 5 is approximated very well by mixed grid calculation. Note, that in the current simulation moisture storage concentrates at wall boundaries. Obviously, our grid refinement is sufficient for reconstructing moisture and vapor pressure profile. In contrast, thermal storage effects demand on a high discretization density of energy balance grid. This fact explains the deviations of medium detailed grid as observed in Fig. 2.

Performance tests
As a next step we study effect of mixed approach on simulation performance. For this purpose, we reproduce current room in two dimensions, see Fig. 7. This modelling approach guarantees each zone to contact minimum two outside walls. We start with a moderate building size: -Medium: Building with outside, inside and floor construction, bar of rooms 16 x 3 in length and height The building geometry is shown in Fig. 7. The test case generates 10400 unknown states including zone and wall energy and moisture balance over a fine wall discretization grid. Mixed discretization reduces number of states to 6731. We simulate over the period of one month (January).  The governing equations form a system of nonlinear differential equations. We solve the equations with the help of time step controlled implicit multi-order integration method CVODE [7]. The CVODE solver transforms the differential equations into a system of nonlinear equations, that are solved by a modified Newton method. The method includes the solution of large linear equation systems. It is well known, that solution of these system is a performance critical step.
In our approach, discretization of wall elements generates large and sparse equation systems. We use the KLU direct linear equation solver [5]. This solution method has proved to be suitable for large scale problems [4]. Solver steps and their percentage of simulation time are monitored in Table 1. KLU performs a linear equation setup (KLU setup) and a solution step (KLU solve). These computations dominate simulation time.
KLU setup includes Jacobian assembly and Jacobian factorization. Further effort is needed for Jacobian condition number calculation, Jacobian copy and reordering operations, that are not shown in this table. Jacobian assembly 50% 57% Jacobian factorization 10% 7% Fig. 9. Simulation performance for medium test size Simulation results are represented in Fig. 9. We monitored simulation speed (simulation time/ real time) every hour of the simulation period. The fine grid approach reaches even a higher simulation speed in the test case. This effect can be explained with the help of detailed performance analysis in Table 1. Jacobian assembly is the most critical simulation step for fine grid approach and even higher for the mixed grid approach.
This observation is clearly related to implementation details: We generate Jacobian with the help of Finite difference quotient. This calculation determines a large number of system function evaluations. Each system function evaluation recalculates heat and moisture fluxes inside zones and walls. The corresponding implementation was optimized for detailed grid, while the mixed grid flux calculation was implemented without regard to efficiency. Latter is both completely performed at fine and coarse approximation grid. For this reason, Jacobian assembly is more expensive than for the fine grid. This effect is responsible for the reduced simulation speed of the mixed grid case.
Indeed, Jacobian factorization is implemented highly effectively and does not exceed 10% of complete simulation time. This indicates an advantageous choice of sparse matrix solver and a precise evaluation of sparse pattern. For the highly optimized NANDRAD solver framework, these solver steps are no longer time critical. Consequently, the evaluation of balance equations obviously remains as optimization potential. An efficient implementation of flux calculation including parallelization is the key for further improvement of simulation performance. These demands are not fulfilled by mixed grid approach. It is common knowledge, that computational effort for Jacobian factorization of direct linear equation solvers may increase drastically for large system size. In order to check significance of this effect, we increased the number of zones until we observed performance losses of the fine grid. This study leads to the second test case: -Large: Building with outside, inside and floor construction, bar of rooms 40 x 8 in length and height The simulation generates 66958 unknown states for detailed and 43352 states for mixed discretization. Performance of the fine grid breaks down, see Fig. 10. It can be stated, that KLU direct linear equation solver performs well, until a large scale system size is reached. Hygrothermal simulations of very complex buildings could therefore benefit from the of reduced number states as generated by the mixed grid approach.
At the current state of art, large scale hygrothermal problems are only generated for testing and have no relevance in practical application. Note further, that the presented result is restricted to a specific solver configuration (CVODE + KLU). More detailed studies including alternative linear equation solvers, i.e. iterative methods should be performed in order to characterize large scale simulation case. This is not the focus of current article.

Conclusions
Consideration of moisture storage is an important aspect for different application cases, i.e. simulation of wall drying under specific inside climate conditions or the control of relative humidity inside rooms with high moisture sources (museums, libraries). Moisture and thermal storage of wall constructions damp temperature and relative humidity peaks inside the room in times of high occupancy and usage. Consequently, detailed and dynamic consideration of wall construction is an essential task of building simulation. The nonlinear hygrothermal equations cannot be solved analytically and need a spatial discretization grid. The simulation program NANDRAD implements these detailed equations. Indeed, number of unknowns drastically increases compared to nodal models.
Typically, a single room simulation will be considered as sufficient and simulation performance is in an acceptable range. However, taking into consideration central ventilation with air conditioning generates a more complex simulation case. Humidity control of whole building is part of air conditioning unit, and multi-zone simulation of the building may be necessary. With regard to future applications, simulation performance is an important aspect of modelling.
We studied different approaches for reducing the number of unknowns in building simulation by applying different simulation grids. First, we reduced discretization density. Indeed, this approach always demands a case specific grid sensitivity analysis. In a next step, we generated a mixed grid as combination of a fine discretization grid for energy and a coarse grid for moisture balance. This idea is motivated by the high dynamics of temperature inside a typical wall construction compared to the slow moisture transport. For moisture balance, we investigated a grid refinement at wall boundaries in order to represent moisture infiltration into the wall with a sufficient accuracy. For our single test case, the mixed grid simulation generates conform results.
Indeed, performance tests show no improvement in the case of moderate size buildings. Using high performance linear solvers, the solution of internal linear equation systems is highly efficient and has only small impact on simulation speed. Contrarily, evaluating the hygrothermal states and fluxes dominates computational effort. Due to an improved performance, an efficient implementation of model equations is of higher priority than the reduction of system unknowns by a mixed simulation grid.
The performance advantage of mixed approach can only be proven in a theoretical large simulation case that will barely appear in planning practice. Note, that the performance improvements are only shown for a single solver configuration: we used the direct sparse solver KLU, that performs expensive matrix operations such as Jacobian factorization. Following a completely different numerical concept, iterative sparse solvers, i.e. preconditioned Krylov subspace methods may perform better for large scale buildings. The performance results for the fine and mixed grid approach may completely vary in another solver configuration.
As a result, the authors clearly prefer to apply a fine discretization grid for coupled hygrothermal equations. The risk of obtaining a bad accuracy with the mixed grid approximation cannot be excluded, while performance is not sufficiently improved. With regard to large problems, we advice to investigate in large scale numerical solution methods. Such methods typically include error tests and therefore, they guarantee simulation speed up with a controlled solution quality.