A spectral method for solving heat and moisture transfer through consolidated porous media

TL;DR

This paper introduces a spectral reduced-order numerical method for efficiently simulating heat and moisture transfer in multilayered porous media, demonstrating high accuracy and computational speed through case studies and experimental validation.

Contribution

It presents a novel spectral reduced-order approach that simplifies the complex coupled heat and moisture transfer equations, improving efficiency and accuracy over traditional methods.

Findings

Method accurately predicts nonlinear heat and moisture transfer.

Spectral approach reduces computational cost significantly.

Validated against experimental data with high reliability.

Abstract

This work presents an efficient numerical method based on spectral expansions for simulation of heat and moisture diffusive transfers through multilayered porous materials. Traditionally, by using the finite-difference approach, the problem is discretized in time and space domains (Method of lines) to obtain a large system of coupled Ordinary Differential Equations (ODEs), which is computationally expensive. To avoid such a cost, this paper proposes a reduced-order method that is faster and accurate, using a much smaller system of ODEs. To demonstrate the benefits of this approach, tree case studies are presented. The first one considers nonlinear heat and moisture transfer through one material layer. The second case - highly nonlinear - imposes a high moisture content gradient - simulating a rain like condition - over a two-layered domain, while the last one compares the numerical prediction against experimental data for validation purposes. Results show how the nonlinearities and the interface between materials are easily and naturally treated with the spectral reduced-order method. Concerning the reliability part, predictions show a good agreement with experimental results, which confirm robustness, calculation efficiency and high accuracy of the proposed approach for predicting the coupled heat and moisture transfer through porous materials.

Tables

Table 1. Table 1. Hygrothermal properties of water.

Property	Value	Unit
Heat capacity, $c_w$	$4180$	$\left[𝖩/\left(\mathrm{𝗄𝗀}\cdot 𝖪\right)\right]$
Latent heat of evaporation, $L_v$	$2.5\cdot {10}^{\mathrm{\hspace{0.17em}6}}$	$\left[𝖩/\mathrm{𝗄𝗀}\right]$
Water gas constant, $R_v$	$462$	$\left[𝖩/\left(\mathrm{𝗄𝗀}\cdot 𝖪\right)\right]$
Density, ${\rho }_l$	$1000$	$\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{𝟥}}\right]$
Saturation pressure, $P_s\left(T\right)$	$997.3\cdot {\left(\frac{T-159.5}{120.6}\right)}^{8.275}$	$\left[\mathrm{𝖯𝖺}\right]$

Table 2. Table 2. Some features of the one-layer heat and moisture transfer case.

	IMEX	Spectral $N=10$
$\Delta t^{\star }$	$1.00\cdot {10}^{-2}$	$1.00\cdot {10}^{-1}$
${\epsilon }_{\infty }^u$	$1.48\cdot {10}^{-5}$	$3.30\cdot {10}^{-5}$
${\epsilon }_{\infty }^v$	$1.43\cdot {10}^{-4}$	$2.31\cdot {10}^{-4}$
DOF	$200$	$20$
CPU time $\left(𝗌\right)$	$6.46$	$0.94$
CPU time $(\%)$	$\mathrm{𝟏𝟎𝟎}$	$\mathrm{𝟏𝟓}$

Table 3. Table 3. Computer run time required for the numerical schemes to perform simulations.

Case study	Spectral	IMEX	Chebfun	Ratio (IMEX/Spectral)
$1$ layer	$0.9𝗌$	$6.4𝗌$	$46.6𝗌$	$7.1$
$2$ layers	$1.8𝗌$	$12.5𝗌$	$227.9𝗌$	$6.9$

Nomenclature
Latin letters
$c_{\mathrm{\hspace{0.17em}0}}$	material specific heat	$\left[𝖩/\left(\mathrm{𝗄𝗀}\cdot 𝖪\right)\right]$
$c_w$	liquid water specific heat	$\left[𝖩/\left(\mathrm{𝗄𝗀}\cdot 𝖪\right)\right]$
$c_M$	moisture storage coefficient	$\left[{𝗌}^{\mathrm{𝟤}}/{𝗆}^{\mathrm{𝟤}}\right]$
$c_T$	energy storage coefficient	$\left[𝖩/\left({𝗆}^{\mathrm{𝟥}}\cdot 𝖪\right)\right]$
$h_M$	convective vapour transfer coefficient	$\left[𝗌/𝗆\right]$
$h_T$	convective heat transfer coefficient	$\left[𝖶/\left({𝗆}^{\mathrm{𝟤}}\cdot 𝖪\right)\right]$
$H_l$	liquid water enthalpy	$\left[𝖩/\mathrm{𝗄𝗀}\right]$
$g_{\infty }$	liquid flow	$\left[\mathrm{𝗄𝗀}/\left({𝗆}^{\mathrm{𝟤}}\cdot 𝗌\right)\right]$
$g$	flow	$\left[\mathrm{𝗄𝗀}/\left({𝗆}^{\mathrm{𝟤}}\cdot 𝗌\right)\right]$
$k_l$	liquid permeability	$\left[𝗌\right]$
$k_M$	moisture transf. coeff. under vap. press. grad.	$\left[𝗌\right]$
$k_{TM}$	heat transf. coeff. under vap. press. grad.	$\left[{𝗆}^{\mathrm{𝟤}}/𝗌\right]$
$k_T$	heat transf. coeff. under temp. grad.	$\left[𝖶/\left(𝗆\cdot 𝖪\right)\right]$
$L$	length	$\left[𝗆\right]$
$L_v$	latent heat of vaporization	$\left[𝖩/\mathrm{𝗄𝗀}\right]$
$𝐧$	normal space, that assumes either $+1$ or $-1$	$[-]$
$P_s$	saturation pressure	$\left[\mathrm{𝖯𝖺}\right]$
$P_v$	vapour pressure	$\left[\mathrm{𝖯𝖺}\right]$
$q$	heat flux	$\left[𝖶/{𝗆}^{\mathrm{𝟤}}\right]$
$R_v$	water gas constant	$\left[𝖩/\left(\mathrm{𝗄𝗀}\cdot 𝖪\right)\right]$
$T$	temperature	$\left[𝖪\right]$
$w$	moisture content	$\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{𝟥}}\right]$
Greek letters
${\delta }_v$	vapour permeability	$\left[𝗌\right]$
$\varphi$	relative humidity	$[-]$
$\rho$	specific mass	$\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{𝟥}}\right]$
$\lambda$	thermal conductivity	$\left[𝖶/\left(𝗆\cdot 𝖪\right)\right]$
Dimensionless parameters
$\mathrm{Bi}$	Biot number	$[-]$
$u$	temperature	$[-]$
$v$	vapour pressure	$[-]$

Table 5. Table 4. Hygrothermal properties of the load bearing material [ 23 ] .

Property	Value	Unit
Volumetric heat capacity	${\rho }_{\mathrm{\hspace{0.17em}0}}{c}_{\mathrm{\hspace{0.17em}0}}=2005\cdot 840$	$\left[𝖩/{𝗆}^{\mathrm{𝟥}}\cdot 𝖪\right]$
Sorption isotherm	$w\left(\varphi \right)=47.1{\left[1+{\left(-1692.94\cdot \ln \left(\varphi \right)\right)}^{1.65}\right]}^{-0.39}+109.9{\left[1+{\left(-2437.83\cdot \ln \left(\varphi \right)\right)}^6\right]}^{-0.83}$	$\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{𝟥}}\right]$
Vapour permeability	${\delta }_v\left(\varphi \right)=6.413\cdot {10}^{-9}\cdot \frac{\left(1-\frac{w\left(\varphi \right)}{157}\right)}{0.503{\left(1-\frac{w\left(\varphi \right)}{157}\right)}^2+0.497}$	$\left[𝗌\right]$
Liquid permeability	$k_l\left(\varphi \right)=2.52\cdot {10}^{-4}\cdot \exp \left(-1.55\cdot {10}^{\mathrm{\hspace{0.17em}6}}\cdot \varphi \right)$	$\left[𝗌\right]$
Thermal conductivity	$\lambda \left(\varphi \right)=0.5+0.0045\cdot w\left(\varphi \right)$	$\left[𝖶/\left(𝗆\cdot 𝖪\right)\right]$

Table 6. Table 5. Hygrothermal properties of the finishing material [ 23 ] .

Property	Value	Unit
Volumetric heat capacity	${\rho }_{\mathrm{\hspace{0.17em}0}}{c}_{\mathrm{\hspace{0.17em}0}}=790\cdot 870$	$\left[𝖩/{𝗆}^{\mathrm{𝟥}}\cdot 𝖪\right]$
Sorption isotherm	$w\left(\varphi \right)=209{\left[1+{\left(-2.7\cdot {10}^{\mathrm{\hspace{0.17em}14}}\cdot \ln \left(\varphi \right)\right)}^{1.27}\right]}^{-0.21}$	$\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{𝟥}}\right]$
Vapour permeability	${\delta }_v\left(\varphi \right)=6.413\cdot {10}^{-9}\cdot \frac{\left(1-\frac{w\left(\varphi \right)}{209}\right)}{0.503{\left(1-\frac{w\left(\varphi \right)}{209}\right)}^{\mathrm{\hspace{0.17em}2}}+0.497}$	$\left[𝗌\right]$
Liquid permeability	$k_l\left(\varphi \right)=\exp \left[-33+0.0704\cdot \left(w-120\right)-1.742\cdot {10}^{-4}\cdot {\left(w-120\right)}^2-2.795\cdot {10}^{-6}\cdot {\left(w-120\right)}^3-1.157\cdot {10}^{-7}\cdot {\left(w-120\right)}^4+2.597\cdot {10}^{-9}\cdot {\left(w-120\right)}^5\right]$	$\left[𝗌\right]$
Thermal conductivity	$\lambda \left(\varphi \right)=0.2+0.0045\cdot w\left(\varphi \right)$	$\left[𝖶/\left(𝗆\cdot 𝖪\right)\right]$

Table 7. Table 6. Hygrothermal properties of the wood fibre [ 45 ] .

Property	Value	Unit
Volumetric heat capacity	$c_{\mathrm{\hspace{0.17em}0}}{\rho }_{\mathrm{\hspace{0.17em}0}}=161.1\cdot {10}^{\mathrm{\hspace{0.17em}3}}$	$\left[𝖩/\left({𝗆}^{\mathrm{𝟥}}\cdot 𝖪\right)\right]$
Sorption isotherm	$f\left(\varphi \right)=7.063\cdot {10}^{-5}\cdot {\varphi }^{\mathrm{\hspace{0.17em}3}}-0.00736\cdot {\varphi }^{\mathrm{\hspace{0.17em}2}}+0.4105\cdot \varphi +0.2688$	$\left[\mathrm{𝗄𝗀}/{𝗆}^{\mathrm{𝟥}}\right]$
Vapour permeability	${\delta }_v\left(\varphi \right)=6.36\cdot \varphi +2.16\cdot {10}^{-11}$	$\left[𝗌\right]$
Thermal conductivity	$\lambda (\varphi ,T)=0.038+0.192\cdot \frac{f\left(\varphi \right)}{{\rho }_l}+1.08\cdot {10}^{-4}\cdot T$	$\left[𝖶/\left(𝗆\cdot 𝖪\right)\right]$

Equations

$$\frac{1}{P_s}\frac{\partial w}{\partial \varphi }\frac{\partial P_v}{\partial t}$$

$$\displaystyle\bigl{(}\,\rho_{\,0}\ c_{\,0}\ +\ w\ c_{\,w}\,\bigr{)}\ \dfrac{\partial T}{\partial t}$$

$$\displaystyle k_{\,M}\$$

$$\text{under vapour pressure gradient},$$

$$\displaystyle k_{\,TM}\$$

$$\displaystyle k_{\,T}\$$

$$\displaystyle c_{\,M}\$$

$$\displaystyle c_{\,T}\$$

$$c_M(P_v)\frac{\partial P_v}{\partial t}$$

$$c_T(P_v)\frac{\partial T}{\partial t}$$

$$\displaystyle\mathbf{n}\cdot\Biggl{(}k_{\,M}\,\dfrac{\partial P_{\,v}}{\partial x}\Biggr{)}$$

$$\mathbf{n}\cdot\Biggl{(}k_{\,T}\,\dfrac{\partial T}{\partial x}\ +\ k_{\,TM}\,\dfrac{\partial P_{\,v}}{\partial x}\Biggr{)}\ =\ h_{\,T}\,\Bigl{(}\,T\ -\ T_{\,\infty}\,(t)\,\Bigr{)}\\ \ +\ L_{\,v}\,h_{\,M}\,\Bigl{(}\,P_{\,v}\ -\ P_{\,v,\,\infty}\,(t)\,\Bigr{)}\ -\ H_{\,l}\,g_{\,\infty}\,(t)\,,$$

$$P_v(x,t=0)$$

$$T(x,t=0)$$

$$P_{v,1}(x_{\text{int}},t)$$

$$T_1(x_{\text{int}},t)$$

$$\displaystyle\Biggl{(}\,k_{\,M,\,1}\ \dfrac{\partial P_{\,v,\,1}}{\partial x}\,\Biggr{)}\Bigg{|}_{x_{\,\text{int}}}$$

$$\displaystyle\Biggl{(}\,k_{\,T,\,1}\ \dfrac{\partial T_{\,1}}{\partial x}\ +\ k_{\,TM,\,1}\ \dfrac{\partial P_{\,v,\,1}}{\partial x}\,\Biggr{)}\Bigg{|}_{x_{\,\text{int}}}$$

$$q_s(t)\stackrel{\mathrm{def}}{:=}-k_T{\frac{\partial T}{\partial x}|}_{x_0}$$

$$\displaystyle g\,(t)\$$

$$\displaystyle u\$$

$$t^{\star }$$

$$\displaystyle k_{\,M}^{\,\star}\$$

$$\displaystyle\mathrm{Bi}_{\,T}\$$

$$\displaystyle g^{\,\star}_{\,\infty}\$$

$$c_M^{\star }\frac{\partial v}{\partial t^{\star }}$$

$$c_T^{\star }\frac{\partial u}{\partial t^{\star }}$$

$$\displaystyle\mathbf{n}\cdot\Biggl{(}k_{\,M}^{\,\star}\,\dfrac{\partial v}{\partial x^{\,\star}}\Biggr{)}\ =\$$

$$\displaystyle\mathbf{n}\cdot\Biggl{(}k_{\,T}^{\,\star}\,\dfrac{\partial u}{\partial x^{\,\star}}\ +\ k_{\,TM}^{\,\star}\,\dfrac{\partial v}{\partial x^{\,\star}}\Biggr{)}\ =\$$

$$\displaystyle\ +\ \mathrm{Bi}_{\,TM}\,\Bigl{(}\,v\ -\ v_{\,\infty}\,(t^{\,\star})\,\Bigr{)}\ -\ H^{\,\star}_{\,l}\,g^{\,\star}_{\,\infty}\,(t^{\,\star})\,,$$

$$u(x^{\star },t^{\star }=0)$$

$$v(x^{\star },t^{\star }=0)$$

$$v_1(x_{\text{int}}^{\star },t^{\star })$$

$$u_1(x_{\text{int}}^{\star },t^{\star })$$

$$\displaystyle k_{\,M,\,1}^{\,\star}\ \dfrac{\partial v_{\,1}}{\partial x^{\,\star}}\Bigg{|}_{x^{\,\star}_{\,\text{int}}}$$

$$\displaystyle\Biggl{(}\,k_{\,T,\,1}^{\,\star}\ \dfrac{\partial u_{\,1}}{\partial x^{\,\star}}\ +\ k_{\,TM,\,1}^{\,\star}\ \dfrac{\partial v_{\,1}}{\partial x^{\,\star}}\,\Biggr{)}\Bigg{|}_{x^{\,\star}_{\,\text{int}}}$$

$$c_M(v)\frac{\partial v}{\partial t}$$

$$c_T(v)\frac{\partial u}{\partial t}$$

$$\displaystyle\mathbf{n}\cdot\Biggl{(}k_{\,M}\,(v)\,\dfrac{\partial v}{\partial x}\Biggr{)}\ =\$$

$$\displaystyle\mathbf{n}\cdot\Biggl{(}k_{\,T}\,(v)\,\dfrac{\partial u}{\partial x}\ +\ k_{\,TM}\,(v)\,\dfrac{\partial v}{\partial x}\Biggr{)}\ =\$$

$$\displaystyle\ +\ \mathrm{Bi}_{\,TM,\,R}\,\Bigl{(}\,v\ -\ v_{\,\infty,\,R}\,(t)\,\Bigr{)}\,.$$

$$\frac{\partial v}{\partial t}$$

Full text

Key words and phrases: Spectral methods; Chebyshev polynomials; Tau–Galerkin method; numerical simulation; heat and moisture transfer; reduced-order modeling

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Mathematical model

3. 2.1 Dimensionless formulation

4. 3 Spectral reduced-order model

5. 3.1 Method description

6. 3.2 Extension to multilayered domains

7. 4 Numerical benchmark

8. 4.1 Single-layered domain

9. 4.2 Multi-layered domain

10. 5 Validation of the model

11. 5.1 Experimental setup

12. 5.2 Simulation

13. 6 Conclusions

14. Acknowledgements

15. A Material properties

16. B Coefficients of the different materials

17. C Dimensionless values

18. C.1 Case from Section 4.1

19. C.2 Case from Section 4.2

20. C.3 Case from Section 5.2

Introduction

Energy consumption of conditioned spaces is strongly related to heat transfer through porous structures, which are dependent of external (weather) and internal conditions. Moreover, moisture migration and accumulation within the porous structures considerably affect the transient hygrothermal performance of porous elements, especially in buildings. Thus, coupled hygrothermal simulations are important to estimate envelope energy performance and risks associated to the presence of moisture such as material degradation, mold growth and related health aspects of occupants.

Initially, the major hypothesis was of no coupling between governing equations of heat and moisture transfers. Although, in porous materials, moisture transfer and accumulation have a direct impact on the heat transfer, especially when there is phase change [37, 5, 43]. As presented by Deru [13], to precisely determine the heat losses, simultaneous calculations with moisture content are required, as they are closely interdependent. Temperature and moisture contents are highly interconnected phenomena and, hence, must be simulated together [4, 33]. For instance, moisture can affect the effective thermal conductivity by a factor of ten.

However, there are some challenges on how to characterize mathematically those physical processes, due to the complexity of the physical phenomena and to the porous structure. Studies of heat and moisture transfer have been published since $1950$’s by Philip and De Vries [40] and Luikov [29]. They represented the coupled processes of heat and moisture transfer by a system of two nonlinear second-order partial differential equations which uses as driving potentials the temperature and the moisture content gradients. The system is nonlinear mainly since of the phase change term in the energy conservation equation and also due to the fact the transport coefficients are highly moisture content and/or temperature dependent. Although, the paramount importance of accurately predict moisture content and temperature fields in several fields of science, detailed simulation of the heat and moisture transfer has only been performed with the improvement of the computer systems. The main reason for that is due to the complexity of the problem and to the high computational cost to solve the coupled heat and moisture transfer equations, which is at least two orders of magnitude greater than that needed to solve only the heat conduction equation [13].

The numerical models used to predict heat and moisture transfer have to deal with multi-scale effects in both space and time domains, with different characteristic times and lengths. For example, simulation of building behavior is generally analyzed on the scale of one year (or more). However, the phenomena and particularly the boundary conditions evolve on a time scale of minutes or even seconds. The geometric configurations of the buildings require three-dimensional modeling of lengthy elements such as facades and ground. Furthermore, when dealing with heat and moisture, the nonlinear behavior of the materials should be taken into account. The combination of all those factors implies the use of more robust and efficient numerical methods since no analytical solution does exist for those problems.

For modeling purposes, the problem previously described is generally solved by the finite-difference methods [19, 16], the finite-volume methods [14, 35, 32] and the finite-element methods [47, 27, 44], which are well-established methods in the field of heat and mass transfer. The basics of a numerical approach is the idea of making an approximation of the solution, which takes a solution for a finite number of degrees of freedom (DOFs). The greater the number of discrete points, the closer to the exact solution will be the approximated solution [31].

However, when the solution of one problem requires long simulation periods (years), considering an entire building, with a very fine time step and mesh refinement, computation becomes too time-consuming. To solve this problem, model reduction techniques can be used as an alternative to approach the solution of the problem and reduce the computational cost.

The intent to construct reduced-order models (ROMs) is to provide an accurate description of the physical phenomena by decreasing the number of degrees of freedom, while retaining the model’s fidelity, at a computational cost much lower than the large original model [42]. In recent years, reduced-order modeling techniques have proven to be powerful tools for solving various problems. Important efforts have been dedicated to developing reduced-order models that can provide accurate predictions while dramatically reducing computational time, for a wide range of applications, covering different fields such as fluid mechanics, heat transfer, structural dynamics, among others [28, 24, 2]. Reduced-order models – such as POD (Proper Orthogonal Decomposition), MBR (Modal Basis Reduction) and PGD (Proper Generalized Decomposition) – have shown a relevant reduction of the computational cost and have been successfully employed by the building physics community [3]. In those works, they have applied reduced-order models to build accurate solutions with less computational effort than the complete original model. Reduced-order models can be classified as a priori or a posteriori methods. The a posteriori approaches need a preliminary computed (or even experimental) solution data of the large original problem to build the reduced one. Whereas the a priori ones do not need preliminary information on the studied problem. The reduced-order model is unknown a priori and is directly built. A careful attention must be paid regarding the definition of ROMs since sometimes it is related to the degradation of the physical model [46], which is not the case of the present work.

Another promising approach to solve the coupled transfer problem is the spectral method, which is a robust and highly accurate method that has been applied to solve partial differential equations since the $70$’s, but lost its spot due to the difficulties to treat nonlinearities, complex geometries, irregular domains, and non-periodic boundary conditions. However, the Spectral methods have overcome some of the mentioned difficulties and now they are successfully applied in studies of wave propagation, meteorology, computational fluid dynamics, quantum mechanics and other fields [8]. The main attractiveness of this method is the superior rate of convergence and the low dissipation and dispersion errors, making its use also attractive to the industry. Nonetheless, spectral methods still have some constraints. For example, complex geometries are one of their main drawbacks as they work better when the geometry of the problem is fairly smooth and regular [6], while finite-element methods are particularly well suited to problems in very complex geometries. Moreover, spectral methods can offer higher accuracy in geometries like boxes and spheres, which can be combined into more complex shapes [17, 9]. In applications where geometry-related disadvantages are not present, the classic finite-element, finite-volume and finite-difference methods do not come close in terms of efficiency.

Some works related to transport phenomena can be found in literature involving diffusive [22, 50], convective [11, 41] and radiative [26, 10, 30] heat transfer. The spectral techniques applied in these works are diversified, adopted according to the geometry, boundary conditions and field of application. In recent works, researchers have implemented spectral methods for solving heat and moisture transfer in food engineering [38] and on fluid flow [36]. Recently, in [20], the authors have studied the moisture transfer in porous building materials considering layered domains, and in [18], they have compared the Spectral method to others ROMs, applied to parametric problems of the building physics field.

Therefore, the scope of this work is to continue the investigations presented in [20] and [18], extending it to the coupled heat and mass transfer. Here, the Spectral method is used to compute one-dimensional heat and moisture diffusion transfer in porous materials, which is validated against experimental data from the literature. The problems treated here involve irregular domains, high nonlinear transport coefficients and non-periodic boundary conditions. The objective is to significantly reduce the computational cost while maintaining high fidelity solutions. This technique assumes separated tensorial representation of the solution by a finite sum of function products. It fixes a set of spatial basis functions to be the Chebyshev polynomials and then, a system of ordinary differential equations is built to compute the temporal coefficients of the solution using the Tau–Galerkin method.

The efficiency of the Spectral approach will be analyzed for simple and multilayered domains with highly nonlinear properties and with sharp boundary conditions and profiles of solutions. For this purpose, the manuscript is organized as follows. First, the description of the physical phenomena is presented (Section 2). Then, the Spectral technique is described (Section 3). In the sequence, the proposed method is applied to two different cases: (i) considering heat and moisture transfer through a single layer (Section 4.1) and (ii) focusing on the heat and moisture transfer through a multilayered domain (Section 4.2). Finally, the method and the model are compared with experimental data from literature in Section 5, considering a single material with real boundary conditions. The main conclusions of the study are outlined in Section 6.

Mathematical model

The physical problem considers heat and moisture transfer through a porous material defined in the one-dimensional spatial domain $\Omega_{\,x}\ =\ [\,0,\,L\,]$ and in the time horizon $\Omega_{\,t}\ =\ [\,0,\,\tau\,]\,$. The moisture transfer occurs due to capillary migration and vapour diffusion. The heat transfer is governed by diffusion and latent mechanisms. The physical model of the problem can be formulated as [3]:

[TABLE]

where $w\ [\mathsf{kg/m^{3}}]$ is the material volumetric moisture content, $\phi\ [-]$, the relative humidity, $\delta_{\,v}\ [\mathsf{s}]$ and $k_{\,l}\ [\mathsf{s}]$, the vapour and liquid permeabilities, $P_{\,v}\ [\mathsf{Pa}]$, the vapour pressure, $P_{\,s}\ [\mathsf{Pa}]$, the saturation pressure, $T\ [\mathsf{K}]$, the temperature, $R_{\,v}\ [\mathsf{J/(kg\cdot K)}]\,$, the water vapour gas constant, $c_{\,0}\ [\mathsf{J/(kg\cdot K)}]$, the material specific heat, $\rho_{\,0}\ [\mathsf{kg/m^{3}}]$, the material specific mass, $\rho_{\,l}\ [\mathsf{kg/m^{3}}]\,$, the water specific mass, $c_{\,w}\ [\mathsf{J/(kg\cdot K)}]\,$, the water specific heat, $\lambda\ [\mathsf{W/(m\cdot K)}]$, the thermal conductivity, and, $L_{\,v}\ [\mathsf{J/kg}]\,$, the latent heat of vaporization. Table 1 presents the values of the water properties considered in this work.

The relation between the moisture content $w$ and the relative humidity $\phi$ is given by the sorption isotherm and the relation between the vapour pressure $P_{\,v}$ and the relative humidity $\phi$ is given by $\phi\,=\,\nicefrac{{P_{\,v}}}{{P_{\,s}(T)}}\,$. In addition, the following assumptions are adopted in this study: (i) no hysteresis effect; (ii) temperature within the range $[0,\,40]^{\circ}\mathsf{C}\,$; (iii) no temperature dependency on the mass balance equation and (iv) properties are dependent only on the vapour pressure field.

Thus, considering the following notation:

[TABLE]

system (2.1) can be rewritten in one-dimensional form as:

[TABLE]

Finally, the problem of interest is a coupled system of two nonlinear parabolic partial differential equations, with vapour pressure $P_{\,v}$ and temperature $T$ gradients as driving potentials.

Boundary conditions.

The moisture exchange between the environment and the surface is driven by vapour exchange, including evaporation and condensation, and by driving rain:

[TABLE]

where $P_{\,v,\,\infty}\ [\mathsf{Pa}]$ stands for the vapour pressure far from the surface and, $h_{\,M}\ [\mathsf{s/m}]$, is the convective moisture transfer coefficient. If the bounding surface is in contact with the outside air, $g_{\,\infty}\ [\mathsf{kg/(m^{2}\cdot s)}]$ is the liquid flow from wind driven rain. The normal $\mathbf{n}$ assumes $+\,1$ or $-\,1$ at the left or right boundary sides.

The heat balance at the boundary includes the convective exchange, the latent heat transfer due to vapour exchange, and the sensible heat transfer due to precipitation, which is expressed as:

[TABLE]

where $T_{\,\infty}\ [\mathsf{K}]$ is the temperature of the air that varies over time and $h_{\,T}\ [\mathsf{W/(m^{2}\cdot K)}]$ is the convective heat transfer coefficient. Regarding to the moisture part, $L_{\,v}\ [\mathsf{J/kg}]$ is the latent heat of vaporization of water and $H_{\,l}\,=\,c_{\,w}\,(T_{\,\infty}\ -\ T_{\,\text{ref}})\ [\mathsf{J/kg}]$ is the liquid water enthalpy, with $T_{\,\text{ref}}=273\,\mathsf{K}\,$.

Initial conditions.

The initial conditions can either have a uniform distribution or a profile more appropriated to the boundary conditions to reduce a warm-up simulation period, which can be very significant depending on the material hygrothermal properties and on the thickness of the building component.

[TABLE]

Interface.

The configuration assumed at the interface between materials follows the hydraulic continuity [12], which considers interpenetration of both porous structure layers. Consider two different materials, both of them are homogeneous and isotropic, and the coupled heat and moisture transfer are simulated, through a perfectly airtight structure. The hydraulic continuity assumes that there is a continuous moisture distribution of vapour content and temperature:

[TABLE]

a continuous moisture flow and a continuous heat flux across the interface verify:

[TABLE]

where $x_{\,\text{int}}\ \in\ \Omega_{\,x}$ represents the location of the interface between materials. Therefore, we can split the spatial domain in two parts, $\Omega_{\,x,\,1}\,=\,[\,0\,,x_{\,\text{int}}\,]$ and $\Omega_{\,x,\,2}\,=\,(\,x_{\,\text{int}}\,,\,L\,]\,$, which represent the spatial domain of material $1$ and material $2\,$.

Fluxes and flows.

One of the interesting outputs in the building physics framework is the heat flux, divided into sensible $q_{\,s}$ and latent $q_{\,l}$ heat fluxes $[\mathsf{W/m^{2}}]$, which are defined as:

[TABLE]

The moisture flow $g$ $[\mathsf{kg/(m^{2}\cdot s)}]$ is similarly computed:

[TABLE]

where $x_{\,0}\,\in\,[\,0\,,L\,]\,$.

Dimensionless formulation

Before solving directly the problem, it is of capital importance to get a dimensionless formulation of the problem under consideration [25]. In this way, we define the following dimensionless quantities:

[TABLE]

where the subscript ref represents a reference value, chosen according to the application problem and the superscript $\star$ represents a dimensionless quantity of the same variable. Therefore, the governing system (2.2) can be written in a dimensionless form as:

[TABLE]

The dimensionless formulation of the boundary conditions are:

[TABLE]

and of the initial conditions are:

[TABLE]

Interface conditions in the dimensionless form are written as:

[TABLE]

In the following, we drop $\star$ for the sake of clarity.

Spectral reduced-order model

Spectral methods consider a sum of polynomials that suit for the whole domain, providing a high approximation of the solution. The smoother a function is, the faster the convergence of its spectral series [6]. For considerably smooth problems, the error decreases exponentially, making the solution with the same order of accuracy of other methods but with a much lower number of degrees of freedom. As a result, this method has a low memory usage, allowing to store and operate a lower number of variables [48]. The Spectral methods used in this work are the Chebyshev polynomials on the basis function and the Tau–Galerkin method to compute the temporal coefficients.

Method description

Problem (2.7) has an important difficulty in dealing with the nonlinearities of the storage $c$ and diffusion $k$ coefficients, all of them depending on the moisture content field. These coefficients are usually given by empirical functions from experimental data. For this reason, Eq. (2.7a) and (2.7b) are recalled with a simplified notation:

[TABLE]

In addition, boundary conditions are simplified and also written with a shorter notation, just for the sake of explaining the method in a pedagogical way/manner:

[TABLE]

A special attention must be given to the spatial domain because the Chebyshev Spectral method is traditionally presented on the canonical interval \big{[}-1\,,1\,\big{]}\,. Thus, if the dimensionless interval is not within \big{[}-1\,,1\,\big{]}\,, a change of variables (domain transformation) must be performed for the computational domain.

In order to apply better the spectral method, Eqs. (3.1a) and (3.1b) are written in the non-conservative form as:

[TABLE]

where,

[TABLE]

The unknowns $u\,(\,x,\,t\,)$ and $v\,(\,x,\,t\,)$ from Eq. (3.3) are accurately represented as a finite sum [34, Chap. 6]:

[TABLE]

Here, $\{\mathsf{T}_{\,i}\,(x)\}_{\,i\,=\,0}^{\,n}$ are the Chebyshev polynomials, $\{a_{\,i}\,(t)\}_{\,i\,=\,0}^{\,n}$ are the corresponding time-dependent spectral coefficients and $n$ represents the number of degrees of freedom of the solution component. Eqs. (3.4a) and (3.4b) can be seen as a series truncation from $N\,=\,n\,+\,1$ modes. The Chebyshev polynomials are chosen as the basis functions since they are optimal in $\mathcal{L}_{\,\infty}$ approximation norm [21]. Therefore, the expression of the derivatives in the Chebyshev basis are:

[TABLE]

where the dot denotes $\dot{a}_{\,i}\,(t)\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\dfrac{\mathrm{d}\,a\,(t)}{\mathrm{d}t}$ according to Newton notation. Note that the derivatives are re-expanded in the same basis function. As a result, coefficients $\{\tilde{a}_{\,i}\,(t)\}$ and $\{\tilde{\tilde{a}}_{\,i}\,(t)\}$ must be re-expressed in terms of coefficients $\{a_{\,i}\,(t)\}\,$. The connection is given explicitly from the recurrence relation of the Chebyshev polynomial derivatives [39]:

[TABLE]

with,

[TABLE]

The derivatives of $u\,(x,\,t)$ are written in a similar way by just replacing $v$ by $u$ and $a$ by $b\,$. Thus, the derivatives defined by Eqs. (3.5a), (3.5b) and (3.5c) are replaced into Eqs. (3.3a) and (3.3b) to provide the residuals:

[TABLE]

which are considered a misfit of the approximate solution. The purpose is to minimize the residual so the solution satisfies the governing equations. To this end, the residual is minimized via the Tau–Galerkin method, which requires that Eqs. (3.6a) and (3.6b) be orthogonal to the Chebyshev basis functions $\langle\,R\,,\mathsf{T}_{\,j}\,\rangle\,=\,0\,$:

[TABLE]

As a result, the project residuals are:

[TABLE]

where, $\mathcal{M}$ is a diagonal and the singular matrix ($\mathrm{rank}\,(\,\mathcal{M}\,)\,=\,N\,-\,2$) which contains the coefficients of the Chebyshev weighted orthogonal system. The matrix $\mathcal{M}$ has the following form:

[TABLE]

and, the matrices with indices ($i\,,$ $j$) written as:

[TABLE]

with indices $i\,$, $j$ are the ones defined in Eqs. (3.4) and (3.7).

By using the Chebyshev–Gauß quadrature, the integrals are also approximated by a finite sum:

[TABLE]

where,

[TABLE]

and $x_{\,k}$ are the Chebyshev nodes:

[TABLE]

The value of $m$ is approximately the same as the number of modes, as discussed in [20].

To complete the problem, the boundary conditions are also written in the form of the residuals. First, the boundary conditions associated to the moisture transport — Eq. (3.2a):

[TABLE]

Then, the boundary conditions for Eq. (3.2c), regarding the heat transport:

[TABLE]

with $T_{\,i}\,(\,-\,1)\,=\,(\,-\,1)^{\,i}$ and $T_{\,i}\,(\,1\,)\,\equiv\,1$ (see [39] for more details).

In this way, it is possible to compose the system of ODEs to be solved, plus the four additional algebraic expressions regarding the boundary conditions. Finally, the system of differential–algebraic equations (DAEs) has the following form:

[TABLE]

where $\mathrm{b\,}_{\,1}\,(t)$ and $\mathrm{b\,}_{\,2}\,(t)$ are vectors containing the boundary conditions, previously defined by $\omega_{\,1}\,$, $\omega_{\,2}\,$, $\kappa_{\,1}$ and $\kappa_{\,2}\,$:

[TABLE]

Matrix $\mathcal{A}$ is written from the right member of Eq. (3.8a), and, matrices $\mathcal{B}$ and $\mathcal{C}$ are written from the right member of Eq. (3.8b).

Initial values of the coefficients $\{a_{\,i}\,(t\,=\,0)\}$ and $\{b_{\,i}\,(t\,=\,0)\}$ are calculated by the orthogonal projection of the initial condition [8]:

[TABLE]

where, $v_{\,0}\,(x)$ and $u_{\,0}\,(x)$, are the dimensionless initial condition.

Therefore, the reduced system of ODEs composed from Eqs. (3.19) and (3.20) can be solved. Different approaches can be used to solve the system of ODEs (3.19). The most straightforward solution is to apply a numerical integration scheme, with moderate accuracy. So, with an embedded error control and not so stringent tolerances, it can be done very efficiently. In this work, the Matlab™ environment was used to perform simulations, and the solvers ODE15s or ODE23t were used to solve the differential-algebraic system of equations (DAEs). The output are the vectors of spectral coefficients $\{a_{\,i}\,(t)\}_{\,i\,=\,0}^{\,n}$ and $\{b_{\,i}\,(t)\}_{\,i\,=\,0}^{\,n}\,$. Then, it enables to reconstruct the solution thanks to spectral representations.

Extension to multilayered domains

We present the idea for two materials, but it can be generalized to any finite number of subdomains straightforwardly. Consider that the original spatial domain $\Omega_{\,x}\,=\,[\,0,\,L\,]$ is decomposed into two subdomains $\Omega_{\,x,\,1}\,=\,[\,0,\,x_{\,\text{int}}\,]$ and $\Omega_{\,x,\,2}\,=\,(\,x_{\,\text{int}},\,L\,]\,$, which represent each material surface and $x_{\,\text{int}}$ represents the location of the interface, as previously defined. From this, the unknowns $v\,(x,\,t)$ and $u\,(x,\,t)$ are then written as:

[TABLE]

for $t\geqslant 0\,$. Thus, $v_{\,1}\,$, $v_{\,2}\,$, $u_{\,1}$ and $u_{\,2}$ are written respectively as:

[TABLE]

Note that the Chebyshev polynomials $\mathsf{T}\,(x)$ are always the same and $x$ must always be in the closed set $[\,-1,\,1\,]\,$. To assure this, the subdomains $\Omega_{\,x,\,1}$ and $\Omega_{\,x,\,2}$ are linearly transformed to the spectral domains $\bar{\Omega}_{\,x,\,1}$ and $\bar{\Omega}_{\,x,\,2\,}$, as illustrated in Figure 1, so they can fit within the interval of interest, which should be seen as a horizontal coordinate in each material. In addition, by separating the real domain in two spectral domains, the smoothness of the solution is assured.

The condition at the interface between the two materials states the continuity of the fields and the flows as in Eq. (2.8). It implies that the derivatives of the fields $u$ and $v$ are not continuous at the interface. This important remark has to be taken into account in the construction of the spectral reduced-order model. Indeed, the domain is decomposed in subdomains to maintain a smooth solution and particularly a continuous derivative on each sub-domain. In this way, the model order reduction is optimal and ensures the error of the Spectral-ROM to decrease exponentially. It is possible to build the reduced order model considering the whole domain (without decomposition). However, the convergence is undermined since the solution and its derivatives are not smooth at the interface between two materials. Many more modes would be necessary to reach the same accuracy, as detailed in Theorem $1$ of [49, Chap. 4].

By considering the two materials, the DAE system represented in Eq. (3.19) becomes:

[TABLE]

The interface conditions – Eq. (2.8) – are then written in the spectral form as:

[TABLE]

which are included in vectors $b_{\,1}$ and $b_{\,2}$ and set equal to zero. In the same way, the boundary conditions are written in the spectral form as:

[TABLE]

which are included in $b_{\,3}$ and $b_{\,4}$ and set equal to zero. Vectors $b_{\,1}$, $b_{\,2}$ , $b_{\,3}$ and $b_{\,4}$ are $N\times 1$ column vectors with the form:

[TABLE]

With all elements listed before, it is possible to set the system to be solved. The system of ODEs for solving the coupled heat and moisture transfer with two layers has the size of $4\cdot N$ with four additional algebraic expressions for the boundary and interface conditions. The initial condition is also given by Eq. (3.20) and the DAE system is solved by the Matlab™ solver ODE15s. In this work, the approach was presented for a two-layered wall for the sake of clarity, knowing it can be easily extended to any number of layers and any other boundary conditions at their interfaces.

Numerical benchmark

To analyze the accuracy of the proposed method, the error between the solution $Y^{\mathrm{num}}$, obtained by the Spectral or IMEX methods, and the reference solutions $Y^{\mathrm{ref}}$, are computed as functions of $x$ using the following Euclidean norms:

[TABLE]

where $N_{\,t}$ is the number of temporal steps and $Y$ can be the temperature $T$ or the vapour pressure $P_{\,v}\,$. The reference solution $Y^{\mathrm{ref}}\,(\,x\,,t\,)$ is computed using the Matlab™ open source toolbox Chebfun [15]. Moreover, the uniform norm error – $\varepsilon_{\,\infty}^{\,Y}$ – is given by the maximal values of $\varepsilon_{\,2}^{\,Y}\,(\,x\,)$ :

[TABLE]

Single-layered domain

Simulations of one-dimensional coupled heat and moisture transport are carried out with the spectral method, the IMplicit-EXplicit (IMEX) method and a reference solution to verify the applicability of the spectral method. The IMEX scheme approximates the continuous operator to order $\mathcal{O}\,\,(\Delta x^{\,2}\ +\ \Delta t)\,$. The advantage of this semi-implicit scheme over the fully implicit one is to avoid sub-iterations in the solution procedure and, at the same time, being stable and consistent.

This case study considers moisture-dependent coefficients $c_{\,M}\,$, $k_{\,M}\,$, $c_{\,T}\,$, $k_{\,T}$ and $k_{\,M}$ as illustrated in Figure 17 of Appendix B. Their variations are similar to the load bearing material from Hagentoft [23], which are presented in Table 4 of Appendix A. The material of $0.1\ \mathsf{m}$ of length has uniform vapour pressure and temperature initial conditions, in which $P_{\,v,\,0}\,=\,1.16\cdot 10^{\,\,3}\,\mathsf{Pa}$ and $T_{\,0}\,=\,293\,\mathsf{K}$. The ambient temperature and relative humidity at the boundaries vary sinusoidally as illustrated in Figures 2(a) and 2(b). The convective vapour transfer coefficients are $h_{\,M,\,L}\,=\,2\cdot 10^{\,-7}\,\mathsf{s/m}$ and $h_{\,M,\,R}\,=\,3\cdot 10^{\,-8}\,\mathsf{s/m}$, for the left and right boundaries, while the convective heat transfer coefficients are set to $h_{\,T,\,L}\,=\,25\,\mathsf{W/(m^{2}\cdot K)}$ and $h_{\,T,\,R}\,=\,8\,\mathsf{W/(m^{2}\cdot K)}$. At the boundaries, only the convective exchange is considered. The simulation is performed for $7\,\mathsf{days}\,$. The dimensionless values can be found in Appendix C.

The Spectral method is composed by $N\,=\,10$ modes with $m\,=\,15\,$. The ODE15s was used to solve System (3.19), with a tolerance set to $10^{\,-\,5}\,$. This solver is adaptive in time, although it can provide the integration in time for the instants required. For this case, the Spectral method was compared to the IMEX and to a reference solution computed using the Chebfun Matlab™ toolbox [15]. All solutions have been computed with $\Delta x^{\,\star}\,=\,10^{\,-\,2}\,$. For the time discretization, the IMEX solution was computed with $\Delta t^{\,\star}\,=\,10^{\,-\,2}$ and the Spectral with $\Delta t^{\,\star}\,=\,10^{\,-\,1}\,$, to have the same order of error of the solution.

Vapour pressure variations at the boundaries are shown in Figure 3(a). In the first four days, it is possible to observe the influence of the initial condition in which the vapour pressure rises significantly, meaning that when the simulation started the material was not in balance with the ambient environment. The vapour pressure at $x\,=\,0\,\mathsf{m}$ and $x\,=\,0.1\,\mathsf{m}$ oscillates according to the boundary conditions while also retaining the moisture. As the convective vapour transfer coefficient is higher at the left side, the material will exchange more with the external environment making variations more visible.

Figure 3(b) shows the temperature variations at the boundaries of the material. The temperature slowly oscillates according to the temperature ambient conditions and also with the vapour pressure variations. As one could expect, vapor pressure and temperature values inversely oscillate.

Results of the error $\varepsilon_{\,2}$ as a function of $x$ are shown in Figure 4(a). The errors associated to the IMEX scheme and to the Spectral method are of order of $\mathcal{O}\,\,(10^{-5}\,)$, for temperature, and of order of $\mathcal{O}\,(10^{-4}\,)$ for moisture. The computer run time used to perform this case and the error $\varepsilon_{\,\infty}$ are given in Table 2, which shows the spectral method is more efficient, even compared to the IMEX, which is an improved method based on the Euler scheme as it does not need sub-iterations. The number of degrees of freedom (DOF) of the Spectral solution is considerably lower if compared to the IMEX solution, which impacts directly on the computational time, making the Spectral approach $7$ times more efficient for this case.

Figure 4(b) displays the last coefficient $a_{\,n}$ of the Spectral-ROM solution for the temperature $u$ and vapour pressure $v$ solutions. The last spectral coefficient is the smallest one and it gives the order of approximation of the residual [6, Page 51]. It works as an upper bound on the error. For example, the error of the $u$ solution in Figure 4(a) is of order of $\mathcal{O}\,\,(10^{\,-5})$, which has the same magnitude of the last coefficient, as can be seen in Figure 4(b).

To illustrate the rate of convergence of the methods, Figure 5(a) presents the error $\varepsilon_{\,\infty}$ as a function of the number of modes $N$ for the Spectral solution, and, as a function of the number of spatial nodes $N_{\,x}$ for the IMEX solution. As the number of modes increases, the solution of the Spectral method converges exponentially and stabilizes with approximatively $13$ modes, reaching an accuracy of order of $\mathcal{O}\,\,(10^{-5})\,$, which is equivalent to the tolerance set on the solver ode15s. On the other hand, the IMEX method converges slower than the Spectral method. It needs at least $100$ spatial nodes to reach the same accuracy of the Spectral method: the CPU time of each simulation was measured and it is presented in Figure 5(b). The computational time increases faster for the Spectral method. However, the spectral solution does not need many modes to converge to an acceptable accuracy ($N\,\simeq\,10$).

The Spectral method has demonstrated good agreement to represent the physical model and the fidelity of the model is not deteriorated with the use of this approach.

Multi-layered domain

This case study considers a porous wall formed by $2$ layers: $8$-$\mathsf{cm}$ load bearing material and $2$-$\mathsf{cm}$ finishing material. Figure 6 shows its schematic representation. The first layer has a faster liquid transfer, while the second layer acts as a hygroscopic finish. The properties of these materials are given in Tables 4 and 5. Considering the temperature range of interest in building applications, temperature dependence was neglected when compared to their dependence on moisture content, with the transport coefficients calculated as a function of the moisture content.

Initial conditions are considered to be uniform over the spatial domain, with the initial temperature of $T_{\,0}\,=\,293.15\ \mathsf{K}$ and the initial vapour pressure of $P_{\,v,\,0}\,=\,1.16\cdot 10^{\,3}\ \mathsf{Pa}$, regarding to a relative humidity of $50\,\%\,$. The boundary conditions oscillate sinusoidally during $7$ days of simulation, which are represented in Figures 2(a) and 2(b). The convective mass and heat transfer coefficients are set to $h_{\,M}\,=\,2\cdot 10^{\,-7}\,\mathsf{s/m}$, $h_{\,M,\,R}\,=\,3\cdot 10^{\,-8}\,\mathsf{s/m}$, $h_{\,T,\,L}\,=\,25\,\mathsf{W/(m^{2}\cdot K)}$ and $h_{\,T,\,L}\,=\,8\,\mathsf{W/(m^{2}\cdot K)}$. The liquid water flow (rain) has two peaks, one at $42\,\mathsf{h}$ and the other one at $126\,\mathsf{h}$, reaching a maximum value of $1.53\cdot 10^{\,-4}\ {\mathsf{kg}}/{\mathsf{(m^{2}\cdot s)}}$ as shown in Figure 7(a), which generates a sensible heat flux of $11\,{\mathsf{W}}/{\mathsf{m^{2}}}$ and $14\,{\mathsf{W}}/{\mathsf{m^{2}}}$ as displayed in Figure 7(b).

Simulations with the Spectral-ROM were performed using the ode15s, with a tolerance set to $\mathsf{tol}\,=\,10^{\,-\,5}\,$, with $N\,=\,8$ modes and $m\,=\,13\,$. The time is incremented with a discretization of $\Delta t^{\,\star}\,=\,10^{\,-\,1}$ and the Spectral solution is project to $N_{\,x}\,=\,101$ spatial nodes. The IMEX solution is computed for a discretization parameter of $\Delta t^{\,\star}\,=\,10^{\,-2}$ and $\Delta x^{\,\star}\,=\,10^{\,-\,2}\,$.

Figure 8(a) presents three profiles of vapour pressure field, at the time instants $t\,=\,\{20,\,30,\,40\}\,\mathsf{h}\,$, which are included on the interval of the first incoming raining flow. Variations of vapour pressure are more significant on the first layer because this material is more hygroscopic than the second one. At the interface, one can notice a discontinuity on the derivative. However, the continuity of the field and the one of the flow are assured. Additionally, Figure 8(b) presents the temperature profiles for the same time instants. Variations on the boundaries occur mainly due to the variations of $P_{\,v,\,\infty}$ in both fields, and not only due to the rain flow, which has an opposite effect on the temperature.

The evolution of the temperature and the vapour pressure at the boundary surfaces ($x\,=\,0\ \mathsf{m}$ and $x\,=\,0.1\ \mathsf{m}$) is shown in Figures 9(b) and 9(a), respectively. The vapour pressure varies according to the sinusoidal fluctuations of the boundary conditions until the rain hits the surface. In the first peak, it is not possible to observe high changes but in the second peak, as the material accumulates moisture, the vapour pressure is suddenly augmented and it diffuses through both layers. In fact, only when the vapour pressure is higher inside both material that is possible to observe the impact of the rain flux because of the material properties, more precisely due to the sorption isotherm. As can be observed in Figure 9(b), the temperature at the boundaries of the composite wall varies according to the gradients of vapour pressure rather than the temperature gradients.

By examining the profiles and evolution of temperature and vapour pressure profiles, it is possible to notice that the solutions of the Spectral and IMEX methods are both in a good agreement with the reference solution given by Chebfun. The distribution of the error $\varepsilon_{\,2}$ as a function of $x$ are shown in Figure 10. The $\mathcal{L}_{\infty}$ error of the spectral solution is $\varepsilon_{\infty}^{\,v}\,=\,2.9\cdot 10^{\,-3}$ for the vapour pressure and $\varepsilon_{\infty}^{\,u}\,=\,1.29\cdot 10^{\,-4}$ for the temperature. These values depend on the order of the time discretization, on the chosen tolerance of the solver and on the number of modes. For the IMEX solution, the $\mathcal{L}_{\infty}$ error is $\varepsilon_{\infty}^{\,v}\,=\,2.6\cdot 10^{\,-3}$ for the vapour pressure and $\varepsilon_{\infty}^{\,u}\,=\,1.1\cdot 10^{\,-4}$ for the temperature. To obtain the same accuracy of the solutions, the Spectral-ROM was $7$ times more efficient than the IMEX approach, as it can be seen in Table 3. One should recall that the IMEX has the same accuracy of the Euler explicit and is more efficient, in terms of computer run time, than the largely used Crank–Nicolson scheme. Just for curiosity, the computational time of the Chebfun simulations have been added. As one can notice they are much higher than the values obtained by the other methods. In fact, Chebfun is made for numerical computing of a wide range of problems. As it does not know about the specific problem, it does not take advantage of the problem structure. Anyway, it is a good tool for comparing purposes.

Figure 11(a) presents the last spectral coefficients of $u$ and $v$ solutions for material $1$, while Figure 11(b) presents the fourth and fifth spectral coefficients for material $2$. The magnitude of the last spectral coefficient acts as an error estimator, determining the error upper limit, as previous mentioned. In material $1$, spectral coefficients needed more modes to have an acceptable solution, while in the other material, an accurate solution can be built with a lower number of modes. This difference occurs due to the properties of material $1$ that are more nonlinear than those of material $2$. One may notice in Figure 11(b) a peak after day $5$, which disappears shortly after. This may occur due to aliasing errors. However, the magnitude of this error is very small if compared with the magnitude of the solution and it does not occur on the other spectral coefficients. As the problem demands only a few modes (around $10$), this kind of error does not really affect the final solution.

In addition, Figure 12 displays the Fourier power spectrum function of the signal frequency per unit of time, on the left and right boundaries, for $u$ and $v$ solutions. Oscillations occurring for $u\,$, from $2\cdot 10^{\,-\,1}\,\mathsf{Hz}$, are attributed to aliasing errors. However, the power of this frequency is very low if compared with the highest peak, corresponding to a difference of $6$ orders of magnitude.

The results for the sensible, latent and total heat fluxes at the left boundary $x\,=\,0\ \mathsf{m}$ are given graphically in Figure 13(a), the boundary that receives the rain flow. The sensible and latent heat fluxes have high values but with opposite signs. Although, as they do not have the same value, they do not cancel each other, as shown for the total heat flux, which is the sum of the latent and sensible heat fluxes. Furthermore, Figure 13(b) presents the total moisture flow at the same boundary. The rain flow is observed in the two maximal peaks, one right after day $1$ and the other right after the day $5\,$. The other variations are caused by the boundary conditions.

Validation of the model

In this section, the physical model and the spectral solution are compared with experimental data gathered from the French project HYGRO-BAT [1]. The measurements were performed at the French laboratory LOCIE (Laboratory of Optimisation of the Conception and Engineering of the Environment) [45]. One-dimensional coupled heat and moisture transfer through a single-layered wall is monitored by sensors placed inside of the material and on its surfaces. In this case, they have not considered the liquid transfer in the moisture balance equation – Eq. (2.1a). Surface sensors provide boundary conditions for the coupled simulation, while the other sensors provide reference measurements for the model validation.

The relative error $\epsilon$ is computed to compare simulations with the experimental data and it is defined as:

[TABLE]

where $Y_{\,k}^{\,\mathrm{num}}$ is the computed solution and $Y_{\,k}^{\mathrm{\,meas}}$ is the measured data.

Experimental setup

The wall considered in this study is composed of a $16-\mathsf{cm}$ layer of wood fibre material, which is subjected to variations in terms of temperature and relative humidity for a $14-$days period. The material properties are given in Table 6. The reference data for the evaluation is provided by temperature and humidity sensors SHT75 Sensirion, located at $x\,=\,\{4,\,8,\,12\}\,\mathsf{cm}$ within the wall.

The sensors have a measurement uncertainty of $\pm\,0.3^{\circ}\mathsf{C}$ for temperature and of $\pm\,0.018$ for relative humidity. Furthermore, the uncertainty regarding the position of the sensors is of $\pm\,1\,\mathsf{cm}$ for the sensors located at $x\,=\,\{4,12\}\,\mathsf{cm}$ and of $\pm\,0.5\,\mathsf{cm}$ for the other sensor. The uncertainty on the positions are different because the sensor at $x\,=\,\{4,12\}\,\mathsf{cm}$ have been settled by perforating a whole in the material layer.

The measured temperature and relative humidity at the interior and exterior boundaries are given in Figures 14(a) and 14(b). The gray color around the curves represent the uncertainties related to the measurements. At the interior boundary, the temperature is set to approximately $24^{\circ}\mathsf{C}$ and the relative humidity set to $40\%$ in the first week and to $70\%$ in the second week. The exterior temperature and relative humidity values are given by their measurement at the boundary, which is filtered by a $2\,\mathsf{cm}$ of a coating layer, excluding solar radiation and driven rain phenomenon. Thus, both boundaries are expressed as Dirichlet-type conditions for the model validation.

Simulation

Simulations are performed for a single-layer wall of a material that separates two well defined environments. The initial condition for temperature and vapour pressure fields are given by an interpolation of the measurements at $x\,=\,\{0,\,4,\,8,\,12,\,16\}\,\mathsf{cm}$ at the first time instant.

The Spectral method is composed of $N\,=\,8$ modes with $m\,=\,13\,$. The ODE15s was used to solve numerically System (3.19), with a tolerance set to $10^{\,-\,5}\,$. The solution has been computed with a time step of $\Delta t^{\,\star}\,=\,0.1\,$, the equivalent of $6\,\mathsf{min}$.

Simulations are compared with the dynamic temperature and vapour pressure. Figure 15 presents the measured data and the simulation results in each one location of the sensors, at $x\,=\,\{4,\,8,\,12\}\,\mathsf{cm}$.

The error between the predicted solution and the experimental observations is between the uncertainties of the sensors during almost all simulation period. For the temperature evolution, the highest discrepancy is when the moisture at the left boundary changes from $40\%$ to $70\%$ of relative humidity. From the $7$th day, the measured temperature rises more than the calculated one as shown in Figure 15(a), which are higher than the uncertainties of the sensors. The difference between the simulated and measured temperature reaches a maximum of $2^{\circ}\mathsf{C}$ in this period. The temperature simulated in at $x\,=\,\{8,\,12\}\,\mathsf{cm}$ are presented in Figures 15(c) and 15(e). The predicted values are closer to the measurements than in $x\,=\,4\,\mathsf{cm}$, and they follow considerably well variations of the outside boundary. It seems that the total diffusion coefficient of the temperature used for the simulation is higher than the real one.

The influence of the step on the relative humidity can strongly be observed at $x\,=\,4\,\mathsf{cm}$ in Figure 15(b) and less far from this boundary. Simulations are able to follow the variations of the measured vapour pressure better than the temperature. In Figures 15(d) and 15(f), the difference between simulations and measurements become more important, with the incoming moisture flow, reaching a maximum difference of $70\,\mathsf{Pa}$ and $80\,\mathsf{Pa}$, respectively. The absolute difference is higher for $x\ =\ 4\,\mathsf{cm}$ as it is closer to the left boundary and consequently to the incoming flow. The discrepancies, come from the physical model, which does not consider liquid transport neither hysteresis. Another explanation is the fact that the properties also were estimated with and admissible error, which influences the predicted solution. This difference on the vapour pressure simulations impacts on the predictions of temperature.

The relative error $\epsilon$ of the computed solutions are present in Figure 16. Solutions simulated with the Spectral method showed a good agreement with the reference data. The maximum relative error for the temperature solution is of $0.7\%$ and for the vapour pressure solution is of $5.6\%$. The average error on the dynamic profiles is $0.2\%$ for temperature and $1.75\%$ for vapour pressure, which are close to the values obtained by Rouchier et al.[45].

The spectral method can provide a very accurate solution of the physical model. In this case, the last spectral coefficient for the temperature is of order of $b_{\,n}\,\simeq\,\mathcal{O}\,(10^{\,-\,5})$ and for the vapour pressure of order of $a_{\,n}\,\simeq\,\mathcal{O}\,(10^{\,-\,3})\,$, for $N\,=\,8\,$. The difference observed in Figure 15(a) can be caused by differences in the mathematical model and by the uncertainties of the estimated properties that are not taken in consideration [7]. Despite of this disparity, simulated results can be considered satisfactory to predict the heat and moisture transport.

As the boundary conditions are of the Dirichlet-type, the IMEX approach is more efficient in this case, taking only $4$ times more CPU time than the Spectral reduced-order model approach. For a single simulation, it may not appear important, however, when the direct model has to be simulated thousands of times as in [45], the gains become much more expressive.

Conclusions

In the present work, we showed that the unsteady heat and moisture transfers problem is solved efficiently with the spectral reduced-order model. This approach has been evaluated on two numerical unidimensional case studies of heat and moisture transfer in porous media. Each case aimed at exciting the nonlinear properties of the material to induce sharp profiles of temperature and vapour pressure. The first case considered a single material layer with sinusoidal boundary conditions. The second case took into account the rain effect at one boundary through a multi-layered material. For all cases, the spectral method has shown a high accuracy and perfect agreement with the respective reference solutions. The maximum global error was of the order of $\mathcal{O}\,(10^{\,-\,3})\,$. The advantage of the proposed method is the low computational burden $7$ times more efficient than the IMEX approach. In sensitivity analysis and inverse problems, when one has to simulate the system of differential equation many times, the spectral approach becomes very attractive.

Besides the numerical benchmark, one case with experimental data was performed to confirm the physical model and the spectral-based solution method. Results of the simulation showed a good agreement with the experimental data, which confirms all the process of simulation. Results can be improved by adding the hysteresis and the liquid transport effects. Although, this information was not available for this material and this experiment.

Further research should be dedicated to multiple space dimensions and to the application considering both the diffusive and advective transfer, which complexity increases considerably due to the nonlinearities of the problem.

Acknowledgements

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – in the framework of the International Cooperation Program CAPES/COFECUB (Grant $\#774/13$). The Authors acknowledge also the support from CNRS/INSIS and Cellule Énergie under the grant MN4BAT-2017. Finally, Professor Mendes thanks the Laboratory LAMA UMR 5127 for the warm hospitality during his visits in 2018, which were supported by the project MN4BAT under the AAP Recherche 2018 programme of the University Savoie Mont Blanc.

Nomenclature

Appendix A Material properties

Appendix B Coefficients of the different materials

See Figure 17.

Appendix C Dimensionless values

Some reference values were used on all the simulation cases. The reference time is $t_{\,\text{ref}}\,=\,1\,\mathsf{h}\,$, the equivalent to $3600\,\mathsf{s}$. The reference temperature was $T_{\,\text{ref}}\,=\,293.15\,\mathsf{K}$ and the reference of the vapour presure was $P_{\,v,\,\text{ref}}\,=\,1166.9\,\mathsf{Pa}\,$. The reference length is the total length of the spatial domain $L_{\,\text{ref}}\,=\,L\,\mathsf{m}\,$, so it is possible to have a dimensionless domain between $x^{\,\star}\,=\,[\,0,1\,]\,$.

Case from Section 4.1

The temperature boundary conditions are expressed as:

[TABLE]

and, of the vapour pressure boundary conditions as:

[TABLE]

For the dimensionless properties of the material, they can be written as:

[TABLE]

with $k_{\,M,\,\text{ref}}\,=\,5.4712\cdot 10^{\,-9}\,\mathsf{s}$ and $k_{\,T,\,\text{ref}}\,=\,0.5021\,\mathsf{W/(m\cdot K)}\,$.

Case from Section 4.2

The dimensionless temperature and vapour pressure at the boundaries are written as in the previous case and also the materials properties of material 1. Thus, for the second material, properties are written as:

[TABLE]

with $k_{\,M,\,\text{ref}}$ and $k_{\,T,\,\text{ref}}$ equal to the previous case. In addition, the rain flow is expressed as:

[TABLE]

Case from Section 5.2

The boundary conditions are gathered from the experimental data and just admensionalized. The dimensionless form of the initial condition are:

[TABLE]

The dimensionless properties of the material can be written as:

[TABLE]

with $k_{\,M,\,\text{ref}}\,=\,3.34\cdot 10^{\,\,-\,11}\,\mathsf{s}$ and $k_{\,T,\,\text{ref}}\,=\,6.98\cdot 10^{\,\,-\,2}\,\mathsf{W/(m\cdot K)}\,$.