An open-source solver for coupled heat and moisture transport in porous media based on OpenFOAM

Dr. Aytaç Kubilay, Email: [email protected]

Contributors:
Dr. Andrea Ferrari, Email: [email protected]

Table of Contents:

• A. Overview
• B. Governing equations and boundary conditions
• C. Implementation in hamFoam
• D. Solver settings
• E. Material properties
• F. Running hamFoam simulation
• References

hamFoam is mainly intended for the absorption, transport and storage simulations of coupled heat and moisture within porous building materials. The source code for the solver hamFoam can be accessed at the following address:
https://github.com/OpenFOAM-BuildingPhysics/hamFoam

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

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

cd /path/to/hamFoam

and compile it together with the necessary libraries by typing:

./Allwmake

The website mentioned above also includes a tutorial case available for download, which models a 1D facade, composed of brick and plaster, that is submitted to subsequent relative humidity, heat and moisture loads at the inner and outer surfaces. This is a case included within the HAMSTAD benchmark cases (case 4: response analysis) (Hagentoft, 2002).

The following directories and files exist in a typical OpenFOAM case folder. Depending on the solvers and utilities used, the file structure can be slightly different and/or additional files can be required.

• 0/: This is the initial time directory, which includes the files for the field variables that are used by the solver. Field variables can have scalar, vector, and tensor values. They are named depending on where they are defined in the mesh, e.g. volVectorField for a vector field defined at cell centers, surfaceScalarField for a scalar field defined at face centers of a given grid cell. The files in the 0/ directory include the initial conditions and boundary conditions for the field variables. New time directories are written as the solver runs depending on the output frequency set by the user.

• constant/polyMesh/: Includes the mesh files.

• system/: The files located here define the solver settings.

  • controlDict: Settings such as timesteps, output frequency, precision, etc.
  • fvSchemes: Discretization schemes for each term in the equations are given here.
  • fvSolution: Linear solver settings, convergence criteria, under-relaxation, etc.

hamFoam uses the continuum modeling approach (Whitaker, 1977), where the different phases are not distinguished separately at a certain point in the material but, instead, the macroscopic behavior of the porous material is modeled. 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) \frac{\partial T}{\partial t} = - \nabla (q_c + q_a)~~~~~~~~~~(1) $$

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

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$ the moisture content, $T$ absolute temperature and $p_c$ capillary pressure. The derivative $\frac{\partial w}{\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.

On the right hand side in Eq. (2), $g_l$ and $g_v$ denote liquid and vapor moisture fluxes [kg/m²s], respectively, as given below:

$$ g_l = - K_l (\nabla p_c - \rho_l g)~~~~~~~~~~(3) $$

$$ g_v = - K_{v,p_c} \nabla p_c - K_{v,T} \nabla T~~~~~~~~~~(4) $$

where $K_l$ denotes the liquid permeability, $g$ the gravitational acceleration $\rho_l$ the liquid density. Vapor flux includes transport due to capillary pressure gradient and due to temperature gradient. $K_{v,p_c}$ and $K_{v,T}$ denote the vapor permeabilities for transport due to capillary pressure gradient and due to thermal gradient, respectively.

In Eq. (1), $q_c$ and $q_a$ denote the conductive and advective heat fluxes [W/m²], respectively:

$$ q_c = - \lambda \nabla T~~~~~~~~~~(5) $$

$$ q_a = (c_lT)g_l + (c_vT + L_v)g_v~~~~~~~~~~(6) $$

where $\lambda$ denotes the thermal conductivity, $c_v$ the specific heat of water vapor and $L_v$ the heat of vaporization. Here, the advective heat transfer represents the sensible and latent heat flow due to vapor and liquid flow.

The governing Eqs. (1) and (2) have been derived under several assumptions:

• no convective air transport occurs
• no liquid transfer due to thermal gradient occurs
• the effect of gravity is negligible
• radiative transfer does not occur
• moisture storage is independent of temperature
• the gaseous phase does not contribute to the storage of moisture and heat
• the temperatures remain below the boiling temperature of water
• thermal equilibrium between all phases
• material properties are independent of temperature

For more detailed information, the reader is referred to Janssen et al. (2007).

For building physics applications, typically, the moisture exchange at the boundaries comprises the convective vapor exchange, $q_{conv}$ [kg/m²s]:

$$ g_{conv} = \beta (p_{v,e} - p_{v,s})~~~~~~~~~~(7) $$

where $\beta$ denotes the convective mass transfer coefficient, $p_{v,e}$ and $p_{v,s}$ the exterior and surface vapor pressure, respectively. In case of exterior boundaries, in addition to the vapor exchange, there is also the wind-driven rain flux, $g_{wdr}$ [kg/m²s].

The heat exchange at the boundaries can be defined by the convective heat transfer, $q_{conv}$ [W/m²], the radiative heat transfer, $q_{rad}$ [W/m²], and the latent and sensible heat transfer due to vapor exchange, $q_{vap}$ [W/m²]. $q_{conv}$ and $q_{vap}$ are defined as below:

$$ q_{conv} = \alpha (T_e - T_s)~~~~~~~~~~(8) $$

$$ q_{vap} = (c_vT_s + L_v) g_{conv}~~~~~~~~~~(9) $$

where $\alpha$ denotes the convective heat transfer coefficient, $T_e$ and $T_s$ the exterior and surface temperature, respectively. In case of wind-driven rain flux, the sensible heat transfer due to rain, $q_{wdr}$, is given as:

$$ q_{wdr} = c_l T_{wdr} g_{wdr}~~~~~~~~~~(10) $$

where $T_{wdr}$ denotes the rainwater temperature.

Transient heat and moisture transport are solved using adaptive timesteps. At each timestep, Eqs. (1) and (2) are solved iteratively based on the Picard iteration scheme (Celia et al., 1990). The Picard iteration is performed until temperature and moisture content values converge (when the changes in the absolute values are smaller than a certain threshold between consecutive Picard iterations) for that timestep. If the solver diverges or does not converge within a given amount of maximum Picard iterations, then the same timestep is tried again with a smaller timestep. If the solver successfully converges within the given maximum Picard iterations, then the solution moves on with a larger timestep for the next iteration. The workflow diagram for the (outer) timestep iteration and the (inner) Picard iteration is given in Fig. 1. This workflow can be found in code file hamFoam.C.

Fig. 1. Workflow diagram for hamFoam solver.

The governing equations for heat and moisture transport are defined in code files TsEqn.H and pcEqn.H, respectively. Note that the formulation given in Eq. (2) is "pc-based", meaning that the dependent variable is capillary pressure for all terms. This formulation is the default option in hamFoam. However, past research shows that in some cases large mass conservation errors can be caused by the pc based formulation (Celia et al., 1990; Janssen et al., 2007; Liu, 2012).The mass conservation errors are mainly due to the nonlinear nature of Eq. (2), the moisture capacity $\frac{\partial w}{\partial p_c}$ strongly depends on $p_c$. Taking a constant moisture capacity for a timestep where it changes rapidly, can lead to a deviation from its exact value. Therefore, in code file pcEqn.H, there is also a "mixed" formulation as given in Eq. (11).

$$ \frac{\partial w}{\partial t} = - \nabla (g_l + g_v)~~~~~~~~~~(11) $$

The use of the "mixed" formulation is explained within the solver settings in part D. Corresponding modification is also made for the heat transport equation in code file TsEqn.H in case the mixed formulation is used.

As for the boundary conditions, the fluxes given in Eqs. (7-10) are defined under folder _LIB/externalLoadModel/derivedFvPatchFields/ using the boundary condition types "HAMexternalHeatFlux" and "HAMexternalMoistureFlux". The boundary conditions require variables, such as $\alpha$, $\beta$, $p_{v,e}$, $T_e$, $q_{rad}$, $R_{wdr}$, $T_{wdr}$, which are read from external text files. For more information, please refer to part E.

An important clarification needs to be done in terms of runoff modeling. In typical building physics applications, the moisture content in the porous building materials is limited with the capillary moisture content. In cases where the surface reaches capillary saturation and the rain flux is larger than the possible absorption by the material, the excess water flows over the surface as runoff. In such a case, the boundary condition switches from a flux condition (Neumann) to a fixed capillary pressure value (Dirichlet), as coded at the end of the code file HAMexternalMoistureFluxFvPatchScalarField.C where the boundary moisture fluxes are defined. A similar modification is done for the heat flux as in the code file HAMexternalHeatFluxFvPatchScalarField.C. Currently, there is no film runoff modeling included in hamFoam and runoff water is assumed to disappear from the system.

Here, several solver settings that are specific to hamFoam are introduced. Some of the settings are specific to the Picard iteration, while others are related to the way hamFoam works. The usage of these settings are listed below, as well as their default values in case they are not entered by the user.

As defined in code file readSolverControls.H, the following Picard controls are defined:

• PicardTolerancews: the convergence criteria for Picard iterative scheme, defined as the maximum difference of absolute moisture content values between consecutive Picard iterations. Default value: 0.01 kg/m³.
• PicardToleranceTs: the convergence criteria for Picard iterative scheme, defined as the maximum difference of absolute temperature values between consecutive Picard iterations. Default value: 0.01 K.
• nIterPicardMax: maximum number of Picard iterations allowed. If the Picard iteration does not converge within this limit, the solver will decrease the timestep and try again. Default value: 10.
• increase_factor: the factor by which to increase the timestep after a successful Picard iteration. Default value: 1.2.
• decrease_factor: the factor by which to decrease the timestep after an unsuccessful Picard iteration. Default value: 0.5.
• 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. Default value: 0.1 s.
• maxDeltaT: the upper limit for the timesteps. Default value: 600 s.

In addition to the Picard controls, there are also the following settings:

• pcEqnForm: defines whether to use the pc-based or mixed formulations for the governing equations. Default is pc-based.
• minCrel: defines whether to use a lower limit for moisture capacity. This is to avoid very small timesteps at very dry states of materials, where the moisture capacity decreases a lot. Default value is VSMALL (1E-37).

These settings can be entered under case folder in file system/controlDict. For the tutorial case HamstadBenchmarkCase4, the following settings are used:

Picard
{
    minDeltaT    1E-6;
    maxDeltaT    30;

    PicardToleranceTs    0.1;
    PicardTolerancews    0.1;
}
pcEqnForm    pc-based; 

hamFoam allows multiple materials to be defined within a given computational mesh. Material properties (both thermal and moisture) for each of these materials are defined under case folder in file constant/transportProperties.

For the tutorial case HamstadBenchmarkCase4, the settings are 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 separate cellZones. The building facade in the benchmark case consists of two layers: 1) brick with a thickness of 10 cm as load bearing material at the exterior side and 2) plaster with a thickness of 2 cm as finishing at the interior side. For each material, the density (rho [kg/m³]) and the specific heat (cap [J/kgK]) of the dry state are defined. Thermal conductivity [W/mK] is defined as:

$$ \lambda = \lambda_1 + \lambda_2 w~~~~~~~~~~(12) $$

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

The entry for "buildingMaterialModel" defines the moisture retention, moisture capacity, liquid permeability and vapor permeability relations. For a complete list of materials, refer to the folder _LIB/buildingMaterialModel. The moisture transport properties are typically results of measurements on each specific material. The material library also includes two generic definitions based on van Genuchten relationships: VanGenuchten and VanGenuchtenVapDiff.

In case file constant/transportProperties, the "name" variable simply acts as the name of the cellZone that selects the cells where a certain material is defined. The cellZones can be generated using the standard OpenFOAM pre-processing tool "setSet". Necessary input is included in case file system/setset.batch in the form as follows:

cellSet loadBearing new boxToCell (0.00 0.00 0)(0.1 1 1)

which defines a cellSet with the name of "loadBearing" within a box defined by the two diagonal points. Then, this cellSet is converted to a cellZone by the following line:

cellZoneSet loadBearing new setToCellZone loadbearing

The name of the cellZone here should be identical to the one defined in constant/transportProperties.

To begin the simulation, one must go into the folder of HamstadBenchmarkCase4:

cd /path/to/HamstadBenchmarkCase4

In the tutorial case, mesh is generated using "blockMesh" utility, which generates a structured hexahedral mesh based on the instructions in case file "system/blockMeshDict". In order to generate the mesh, type:

blockMesh

which will create the computational grid for the 1D building facade. Mesh files are written under "constant/polyMesh/". Once the mesh is generated, the cellZones can be generated by running:

setSet –batch system/setset.batch

The necessary text files for the boundary conditions are given under time folder 0/ for each boundary under a subfolder with the boundary name, e.g. for exterior boundary "0/outside/" and for interior boundary "0/inside/". Fig. 2 shows the variation of exterior equivalent temperature, $T_e$, and the wind-driven rain flux, $R_{wdr}$ during a duration of 5 days. The remaining variables are constant for the whole duration with the following values: convective heat transfer coefficient $\alpha$ = 25 W/m²K, convective moisture transfer coefficient $\beta$ = 2E-7 s/m, exterior vapor pressure $p_{v,e}$ = 1150 Pa, rainwater temperature $T_{wdr}$ = 10 °C. Note that the exterior equivalent temperature already accounts for ambient air temperature and short and long wave radiation. Therefore, the radiative heat flux is set to 0.

Fig. 2. Climatic loads at the exterior boundary.

Fig. 3. shows the variation of interior vapor pressure, $p_{v,i}$. The remaining variables are constant for the whole duration with the following values: convective heat transfer coefficient $\alpha$ = 8 W/m²K, convective moisture transfer coefficient $\beta$ = 3E-8 s/m and interior air temperature $T_i$ = 20 $^{\circ}$C. Interior radiative heat flux is set to 0.

Fig. 3. Climatic loads at the interior boundary.

After setting the boundary conditions, the simulation can be started (in serial calculation mode) by executing:

hamFoam

Note that, for simplicity, this manual is written for serial calculation. For parallel computation, please refer to OpenFOAM documentation.

The variation of surface temperature and moisture content values over time can be viewed in ParaView by typing:

paraFoam

Fig. 4. Variations of surface temperature and moisture content at the exterior boundary.

The resulting temperature and moisture content at the exterior surface are given in Fig. 4. At the beginning, as the exterior temperature decreases, surface temperature decreases and condensation occurs at the surface. This is visible by the increase in surface moisture content. In the following days, there is a cycle of wetting by wind-driven rain and drying by solar radiation. At certain times, the exterior surface reaches the capillary saturation value of brick, where the excess water runs off. Similarly, the variations of temperature and moisture content at the interior surface are given in Fig. 5.

Fig. 5. Variations of surface temperature and moisture content at the interior boundary.

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.

Hagentoft, C-E., 2002. HAMSTAD – Final report: Methodology of HAM-​modeling, Report R-​02:8. Gothenburg, Department of Building Physics, Chalmers University of Technology.

Liu, X., 2012. suGWFoam: an open source saturated-unsaturated ground water flow solver based on OpenFOAM, Civil and Environmental Engineering Studies Technical Report No. 01-01, Department of Civil and Environmental Engineering, University of Texas at San Antonio.

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