Solving nonlinear diffusive problems in buildings by means of a Spectral Reduced-Order Model

TL;DR

This paper introduces a Spectral Reduced-Order Model for simulating nonlinear moisture transfer in building materials, achieving high accuracy and significant computational efficiency compared to traditional methods.

Contribution

The paper presents a novel Spectral ROM approach for nonlinear diffusive problems in buildings, demonstrating superior accuracy and computational speed over classical schemes.

Findings

Spectral ROM accurately simulates moisture transfer in porous materials.

The method reduces CPU time by up to 100 times in nonlinear cases.

Performance is validated against classical Euler and Crank-Nicolson schemes.

Abstract

This paper proposes the use of a Spectral method to simulate diffusive moisture transfer through porous materials as a Reduced-Order Model (ROM). The Spectral approach is an a priori method assuming a separated representation of the solution. The method is compared with both classical Euler implicit and Crank-Nicolson schemes, considered as large original models. Their performance - in terms of accuracy, complexity reduction and CPU time reduction - are discussed for linear and nonlinear cases of moisture diffusive transfer through single and multi-layered one-dimensional domains, considering highly moisture-dependent properties. Results show that the Spectral reduced-order model approach enables to simulate accurately the field of interest. Furthermore, numerical gains become particularly interesting for nonlinear cases since the proposed method can drastically reduce the computer run time, by a factor of 100, when compared to the traditional Crank-Nicolson scheme for one-dimensional applications.

Tables

Table 1. Table 1. Absolute error ε ∞ subscript 𝜀 \varepsilon_{\,\infty} for different number of modes N 𝑁 N and different truncations m 𝑚 m\, .

	$N=5$	$N=6$	$N=7$	$N=8$	$N=9$
$m=2$	$7.0\cdot {10}^{-2}$	—	$2.90\cdot {10}^{-1}$	—	—
$m=3$	$2.8\cdot {10}^{-3}$	$2.02\cdot {10}^{-2}$	$1.94\cdot {10}^{-2}$	—	—
$m=4$	$2.7\cdot {10}^{-3}$	$1.43\cdot {10}^{-3}$	$3.06\cdot {10}^{-3}$	$2.46\cdot {10}^{-2}$	$2.89\cdot {10}^{-1}$
$m=5$	$2.6\cdot {10}^{-3}$	$1.54\cdot {10}^{-3}$	$9.30\cdot {10}^{-4}$	$1.05\cdot {10}^{-3}$	$4.13\cdot {10}^{-3}$
$m=6$	$2.6\cdot {10}^{-3}$	$1.39\cdot {10}^{-3}$	$7.07\cdot {10}^{-4}$	$3.40\cdot {10}^{-4}$	$4.41\cdot {10}^{-4}$
$m=7$	$2.6\cdot {10}^{-3}$	$1.40\cdot {10}^{-3}$	$6.59\cdot {10}^{-4}$	$3.20\cdot {10}^{-4}$	$2.60\cdot {10}^{-4}$
$m=8$	$2.6\cdot {10}^{-3}$	$1.39\cdot {10}^{-3}$	$6.70\cdot {10}^{-4}$	$3.20\cdot {10}^{-4}$	$1.90\cdot {10}^{-4}$
$m=9$	—	—	—	$3.09\cdot {10}^{-4}$	$2.40\cdot {10}^{-4}$
$m=10$	—	—	—	$3.09\cdot {10}^{-4}$	$2.40\cdot {10}^{-4}$

Table 2. Table 2. Computational time required for the numerical schemes perform the nonlinear case ( Δ ​ x ⋆ =  10 − 2 Δ superscript 𝑥 ⋆ superscript 10 2 \Delta x^{\,\star}\,=\,10^{-2} and Δ ​ t ⋆ =  10 − 2 Δ superscript 𝑡 ⋆ superscript 10 2 \Delta t^{\,\star}\,=\,10^{-2} ).

Numerical Scheme	CPU time ($𝗌$)	CPU time ($\%$)	Average number of iterations
Spectral $N=\mathrm{\hspace{0.17em}9}$	$3$	$1$	—
Chebfun (Reference)	$96$	$29$	—
Crank–Nicolson	$327$	$100$	1

Equations

$$\frac{\partial {\rho }_{l+v}}{\partial t}=\frac{\partial }{\partial x}\left(k_l\frac{\partial P_c}{\partial x}+k_v\frac{\partial P_v}{\partial x}\right),$$

$$P_c$$

$$\frac{\partial P_c}{\partial P_v}$$

$$\frac{\partial P_c}{\partial x}=\underset{=\frac{R_{v}T{\rho }_l}{P_v}}{\underbrace{\frac{\partial P_c}{\partial P_v}}}\cdot \frac{\partial P_v}{\partial x}+\frac{\partial P_c}{\partial T}\cdot \frac{\partial T}{\partial x}.$$

$$\frac{\partial {\rho }_{l+v}}{\partial t}=\frac{\partial {\rho }_{l+v}}{\partial \varphi }\cdot \frac{\partial \varphi }{\partial P_v}\cdot \frac{\partial P_v}{\partial t}+\frac{\partial {\rho }_{l+v}}{\partial T}\cdot \frac{\partial T}{\partial t}.$$

$$\frac{\partial P_c}{\partial x}$$

$$\frac{\partial {\rho }_{l+v}}{\partial t}$$

$$\frac{\partial {\rho }_{l+v}}{\partial t}=f^{\prime }(\varphi )\frac{1}{P_s}\cdot \frac{\partial P_v}{\partial t}.$$

$$\displaystyle f^{\,\prime}\,(\,\phi\,)\,\frac{1}{P_{\,s}}\cdot\frac{\partial P_{\,v}}{\partial t}\ =\ \frac{\partial}{\partial x}\biggl{[}\,\biggl{(}\,k_{\,l}\,\frac{R_{\,v}\,T\,\rho_{\,l}}{P_{\,v}}\ +\ k_{\,v}\,\biggr{)}\cdot\frac{\partial P_{\,v}}{\partial x}\,\biggr{]}\,.$$

$$\displaystyle c_{\,m}\,\frac{\partial P_{\,v}}{\partial t}\ =\ \frac{\partial}{\partial x}\biggl{[}\,d_{\,m}\,\frac{\partial P_{\,v}}{\partial x}\,\biggr{]}\,.$$

$$d_m\frac{\partial P_v}{\partial x}$$

$$-d_m\frac{\partial P_v}{\partial x}$$

$$P_v(x,t=0)=P_v^i(x).$$

$$u=\frac{P_v}{P_v^0},$$

$$d_m^{\star }=\frac{d_m}{d_m^0},$$

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

$$d_m^{\star }\frac{\partial u}{\partial x^{\star }}$$

$$-d_m^{\star }\frac{\partial u}{\partial x^{\star }}$$

$$u$$

$$\frac{\partial u}{\partial t}=\nu \frac{{\partial }^{2}u}{\partial x^2},$$

$$\frac{\partial u}{\partial x}$$

$$-\frac{\partial u}{\partial x}$$

$$u(x,t)=\sum _{i=0}^{\infty }{a}_i(t){\phi }_i(x)=\underset{=u_n(x,t)}{\underbrace{\sum _{i=0}^n{a}_i(t){\phi }_i(x)}}+\underset{\ll 1}{\underbrace{\sum _{i=n+1}^{\infty }{a}_i(t){\phi }_i(x)}}.$$

$${\phi }_i(x)\equiv {\mathsf{T}}_i(x).$$

$${\mathsf{T}}_0(x)=1,$$

$${\mathsf{T}}_{i+1}(x)=2x{\mathsf{T}}_i(x)-{\mathsf{T}}_{i-1}(x).$$

$$\frac{\partial u_n}{\partial x}$$

$$\frac{{\partial }^2{u}_n}{\partial x^2}$$

$$\frac{\partial u_n}{\partial t}$$

$${\tilde{a}}_i=\frac{2}{c_i}\sum _{\begin{array}{c}p=i+1\\ p+i\text{odd}\end{array}}^{n}p{a}_p,$$

$${\tilde{a}}_n\equiv 0,$$

$$\displaystyle\tilde{\tilde{a}}_{\,i}\ =\ \dfrac{1}{c_{\,i}}\sum_{\begin{subarray}{c}p\,=\,i\,+\,2\\ p\,+\,i\;\text{even}\end{subarray}}^{\,n}\,p\,\Bigl{(}\,p^{\,2}\ -\ i^{\,2}\,\Bigr{)}\,a_{\,p}\,,$$

$${\widetilde{\stackrel{~}{a}}}_{n-1}\equiv {\widetilde{\stackrel{~}{a}}}_n\equiv 0,$$

$$c_i=\{\begin{array}{ccc}2,& \text{if}& i=0,\\ 1,& \text{if}& i>0.\end{array}$$

$$\displaystyle R\,(\,x\,,\,t\,)\ =\ \sum_{i\,=\,0}^{n}\,\Bigl{[}\,\dot{a}_{\,i}\,(\,t\,)\ -\ \nu\ \tilde{\tilde{a}}_{\,i}\,(\,t\,)\,\Bigr{]}\,\mathsf{T}_{\,i}\,(\,x\,)\,,$$

$$\displaystyle\Bigl{\|}\;R\,(\,x\,,\,t\,)\;\Bigr{\|}_{\,2}\ \longrightarrow\ \min\,,$$

Full text

Key words and phrases: Spectral methods; Chebyshev polynomials; Tau-Galerkin method; numerical simulation; diffusive phenomena; reduced-order modelling

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Moisture transfer in porous materials

3. 3 Spectral reduced-order model for linear transfer

4. 3.1 Method description

5. 3.2 Validation of the numerical solution

6. 4 Numerical application

7. 4.1 Linear case

8. 4.2 Weakly nonlinear case

9. 4.2.1 Case study

10. 4.3 Numerical cost estimation

11. 4.3.1 Solving the system of ODEs

12. 5 Treating general nonlinearities

13. 5.1 A highly nonlinear case

14. 6 Multilayer domain

15. 6.1 Adaptation of the reduced Spectral Method

16. 6.2 A multilayer case

17. 7 Conclusions

18. Acknowledgments

19. Nomenclature

20. A Dimensionless values

21. A.1 Case from Section 4.1

22. A.2 Case from Section 4.2

23. A.3 Case from Section 5.1

24. A.4 Case from Section 6.2

Introduction

Moisture transfer through porous materials is a matter of concern in many areas, such as building physics, food engineering, hydrology, agriculture, geophysics, environmental engineering, energy systems, among others, where the transient evolution of moisture may play a role of paramount importance. Particularly, in the area of building physics, moisture transfer process through the porous envelope, roofing systems and ground may strongly affect energy and hygrothermal performance of those elements and, at the same time, it can influence the building occupants’ health, the material’s durability and the energy consumption and demand of the edifice.

The mechanisms that control the transport of moisture in porous building materials occurs simultaneously in its different phases. In the vapour phase, the moisture transfer is mostly governed by diffusive and convective transport while in the liquid phase it is governed mainly by capillarity, which is strongly influenced by weather conditions [13].

Over the last decades, several models were proposed in the literature to mathematically describe the moisture transport as described in [30] and for the assessment of moisture effects, numerical tools have been developed to accurately simulate the processes of moisture transfer in building materials [4]. Since the $1990$’s, many computer-based tools for the prediction of the hygrothermal performance were developed, such as DELPHIN [5], MATCH [41], MOIST [8], WUFI [16] and UMIDUS [38, 32]. Moisture models have also been implemented in whole-building simulation tools and tested in the frame of the International Energy Agency Annex $41$, which reported on most of the detailed models and their successful applications for accurate assessment of hygrothermal transfer in buildings [51].

As building material properties are temperature- and moisture-dependent and the boundary conditions are driven by weather variables, the models included in those tools are based on numerical approaches using discrete representations of the continuous equations. To compute the solution, standard discretisation and incremental techniques are applied such as the Euler implicit scheme in [31, 16, 5, 46, 42, 22, 23] to solve large systems of equations (of an order of $10^{\,6}$ for three-dimensional problems). Furthermore, when dealing with nonlinearities, hygrothermal properties of porous materials have to be updated as a function of the temperature and moisture content fields at each iteration. The difficulties to compute the solution increase, particularly when using implicit schemes that require sub-iterations to treat those issues. In the literature, the important numerical costs of simulation tools [14, 12, 33, 2] is also mentioned and it is a matter of concern due to the great scale of buildings, where heat and moisture transfer phenomena have to be simulated. For those reasons, innovative and efficient ways of numerical simulation are worthy of further investigation and model reduction techniques can be an interesting alternative approach to deal with this problem.

The intent in constructing reduced-order models (ROMs) is to provide accurate description of the physical phenomena by decreasing the degrees of freedom, while retaining the model’s fidelity, at a computational cost much lower than the large original model [40]. 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 [26, 21, 3]. Examples with finite-element and finite-volume applications can be found in [49] and [45], respectively. In their work, they apply reduced-order models to build accurate solutions with less computational effort than the large original model. A careful attention must be paid regarding the definition of ROMs since sometimes it is related to degradation of the physical model [43], which is not the case of the present work.

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. Since the $2000$th, aiming to reduce the computational cost, reduction model techniques started to take place in the context of heat and moisture transfer for building physics applications, as an alternative to traditional methods. Different kinds of approaches can be considered, such as the a posteriori Proper Orthogonal Decomposition (POD), the Modal Basis Reduction (MBR) and the a priori Proper Generalized Decomposition (PGD), which has shown a relevant reduction of the computational cost for successful applications in the building physics area [6].

Spectral methods are successfully applied in studies of wave propagation, meteorology, computational fluid dynamics, quantum mechanics and several other fields [9]. Some works on the transport phenomena can be found in literature involving diffusive [19, 50], convective [11, 39] and radiative [24, 10, 29] heat transfer. Spectral techniques applied in these works are varied, 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 [36] and on fluid flow [34]. According to the authors’ knowledge, there is no research in the literature so far regarding the application of spectral methods for solving diffusive moisture transfer in building physics.

Therefore, the scope of this work is to present an innovative approach, applied for the first time in the context of building physics, i.e., the a priori Spectral reduced-order model technique. In this work, the method is used to compute one-dimensional moisture diffusion in porous materials. 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. In this work, aiming at proposing the use of a Spectral method to simulate a complex and time-consuming phenomenon of diffusive moisture transfer, the temperature effect has been disregarded, but it must be the next step of the investigation on a highly efficient numerical method to simulate combined heat and moisture transport through porous building elements.

The efficiency of the Spectral approach will be analyzed/proven/demonstrated for simple and multilayered domains with highly nonlinear properties 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 four different cases in one dimension. The first one considers linear transfer (Section 4.1) to validate the method. The second one focuses on a weak nonlinear transfer (Section 4.2), in which some simplifications are considered, while the third one presents a strongly nonlinear transfer case with moisture-dependent material properties (Section 5). Finally, the last case study (Section 6.2) considers a multilayered wall with important interface conditions imposed.

It is important to mention that the case studies of this article are directly related to two recent papers, focusing on the establishment of efficient numerical models for nonlinear moisture transfer in terms of accuracy and reduced computational effort. The first work [22] discusses the choice of the physical potential for the formulation of the physical problem, while the second one [17] provides a discussion of the numerical methods enabling to build a large original model to solve a nonlinear moisture diffusion problem.

Moisture transfer in porous materials

The physical problem involves one-dimensional moisture diffusion through a porous material defined by the spatial domain $\Omega_{\,x}\,=\,[\,0,\,L\,]\,$. The moisture transfer occurs according to liquid and vapour diffusion processes. The physical problem can be formulated as [27, 1]:

[TABLE]

where $\rho_{\,l+v}$ is the volumetric moisture content of the material and $k_{\,v}$ and $k_{\,l}\,$, the vapour and liquid permeabilities.

Eq. (2.1) can be written using the vapour pressure $P_{\,v}$ as the driving potential. For this, we consider the physical relation, known as the Kelvin equation, between $P_{\,v}$ and $P_{\,c}\,$:

[TABLE]

where $\phi\,=\,\nicefrac{{P_{\,v}}}{{P_{\,s}\,(\,T\,)}}$ is the relative humidity. Thus, the derivative of $P_{\,c}$ regarding $P_{\,v}$ is expressed as:

[TABLE]

and the derivative regarding $x$ can be written as:

[TABLE]

In addition, the left-hand term of Eq. (2.1) can also be expressed in terms of $P_{\,v}$ and $T\,$:

[TABLE]

As the problem has been assumed isothermal, the temperature derivatives vanishes so that Eqs. (2.2) and (2.3) can be written as:

[TABLE]

Considering the relation $\rho_{\,l+v}\ =\ f\,(\,\phi\,)\ =\ f\,(\,P_{\,v}\,,T\,)$, obtained from material properties and from the relation $\phi\,=\,\nicefrac{{P_{\,v}}}{{P_{\,s}(T)}}$ between the vapour pressure $P_{\,v}$ and the relative humidity $\phi$, we get:

[TABLE]

Eq. (2.1) can be therefore rewritten as:

[TABLE]

The material properties $f\,$, $k_{\,l}$ and $k_{\,v}$ depend on the vapour pressure $P_{\,v}$. Therefore, we denote $d_{\,m}\mathop{\stackrel{{\scriptstyle\ \mathrm{def}}}{{:=}}\ }k_{\,l}\,\dfrac{R_{\,v}\,T\,\rho_{\,l}}{P_{\,v}}\ +\ k_{\,v}$ as the global moisture transport coefficient and $c_{\,m}\mathop{\stackrel{{\scriptstyle\ \mathrm{def}}}{{:=}}\ }f^{\,\prime}\,(\,\phi\,)\;\dfrac{1}{P_{\,s}}$ as the moisture storage coefficient. Thus, considering the previous notation, Eq. (2.4) can be written as:

[TABLE]

At the material bounding surfaces, Robin-type boundary conditions are considered:

[TABLE]

where $P_{\,v,\,\mathrm{L}}$ and $P_{\,v,\,\mathrm{R}}$ are the vapour pressure of the ambient air, $g_{\,l,\,\mathrm{L}}$ and $g_{\,l,\,\mathrm{R}}\,$, the liquid water flow (driving rain) and $R$ and $L$ stand for the right and left bounding surfaces. For the initial condition, the vapour pressure distribution is written in function of the space:

[TABLE]

The initial condition 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 (at the order of years) depending up on the material hygrothermal properties and on the thickness of the building component.

It is important to obtain a unitless formulation of governing equations while performing mathematical and numerical analysis of given practical problems, due to a certain number of reasons already discussed in [17]. Therefore, we define the following dimensionless parameters:

[TABLE]

where the superscript [math] represents a reference value, chosen according to the application problem and the superscript $\star$ represents a dimensionless quantity of the same variable. In this way, the dimensionless problem is written as:

[TABLE]

Finally, this is the problem of interest considered here for resolution. Now, the procedure of the Spectral method will be described to propose a Reduce Order Model for the solution of this problem.

Spectral reduced-order model for linear transfer

While finite-difference and finite-element methods are based on a local representation of functions, using low-order approximations, Spectral methods consider a global representation of the solution, which yields beyond all orders approximations [7]. In the global representation approach, the value of the derivative at a certain spatial location depends on the solution on the entire domain and not only on its neighbours. Spectral methods consider a sum of polynomials that suit for this whole domain, almost like an analytical solution, providing a high approximation of the solution. As its error decreases exponentially it is possible to have the same accuracy of other methods but with a lower number of modes, which makes this method memory usage minimized, allowing to store and operate a lower number of degrees of freedom [47]. 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

For the sake of simplicity and without loosing the generality, this method is first explained considering the dimensionless coefficients $d_{\,m}^{\,\star}$ and $c_{\,m}^{\,\star}$ as constants, noting $\nu\,\mathop{\stackrel{{\scriptstyle\ \mathrm{def}}}{{:=}}\ }\,\dfrac{d_{\,m}^{\,\star}}{c_{\,m}^{\,\star}}$ and thus, considering the linear diffusion equation:

[TABLE]

for $t\,>\,0\,$ and x\,\in\,\big{[}-1\,,1\,\big{]}\,; the $\star$ symbol was dropped for the purpose of conciseness to explain the method. A special attention must be given to the spatial domain, because the Chebyshev Spectral method we use is described between the interval \big{[}-1\,,1\,\big{]}\,. Thus, if the dimensionless interval is not in this interval, a change of variables (domain transformation) must be performed for the computational domain.

The boundary conditions are written as:

[TABLE]

The Spectral method assumes that the unknown $u\,(\,x,\,t\,)$ from Eq. (3.1) can be accurately represented as a finite sum [30]:

[TABLE]

Here, $\{\varphi_{\,i}\,(\,x\,)\}_{\,i\,=\,0}^{\,n}$ is a set of basis functions that remains constant in time, $\{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. Eq. (3.3) can be seen as a series truncation after $N\,=\,n\,+\,1$ modes. The Chebyshev polynomials are chosen as the basis functions since they are optimal in $\mathcal{L}_{\infty}$ approximation norm [18]. It should be observed that other bases can be used, such as the Fourier and Legendre polynomials. Therefore, we have:

[TABLE]

The first Chebyshev polynomials are:

[TABLE]

and, higher order polynomials can be constructed using a recursive relation [37]:

[TABLE]

As we have chosen the basis functions, now we can write the derivatives:

[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 Chebyshev 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 [37]:

[TABLE]

with,

[TABLE]

Using the expression of the derivatives provided by Eqs. (3.4b) and (3.4c), the residual of the diffusion equation (3.1) is:

[TABLE]

which is considered a misfit of the approximate solution. The purpose is to minimize the residual:

[TABLE]

which is realized via the Tau–Galerkin method, which requires Eq. (3.5) to be orthogonal to the Chebyshev basis functions $\langle\,R\,,\mathsf{T}_{\,i}\,\rangle\,=\,0\,$:

[TABLE]

namely,

[TABLE]

Then, taking the orthogonality property of the Chebyshev polynomials into account [37], it leads to the following relation among the spectral coefficients:

[TABLE]

Finally, after the projection and expansion of the residual, the result is a system of Ordinary Differential Equations (ODEs), with $n\,-\,2$ equations to be solved as a function of time. The two extra coefficients are obtained by substituting the derivative (3.4a) into the boundary conditions (3.2):

[TABLE]

with $\mathsf{T}_{\,i}\,(-1)\,=\,(-1)^{\,i}$ and $\mathsf{T}_{\,i}\,(\,1\,)\,\equiv\,1$ (see [37]). Eqs. (3.7a) and (3.7b) are written in an explicit way, with coefficients $a_{\,n}$ and $a_{\,n-1}$ expressed in terms of all the other coefficients.

Therefore, the original partial differential equation (3.1) is reduced to a system of ODEs plus two algebraic expressions. For linear problems, the system of ODEs is explicitly built. Moreover, the reduced system of ordinary differential equations has the following form:

[TABLE]

where, $\mathcal{A}\in\mathrm{Mat}\,_{(n-2)\times(n-2)}(\mathds{R})\,$, with constant coefficients and with $n\simeq\,\mathcal{O}\,(\,10\,)\,\,$. Besides, $\mathrm{b\,}(\,t\,)\in\mathds{R}^{\,(n-2)}$ is a vector coming usually from boundary conditions.

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

[TABLE]

where, $u_{\,0}\,(\,x\,)$, is the dimensionless initial condition. After solving the reduced system of ODEs (Eqs. (3.8) and (3.9)), it is possible to compose the solution along with the Chebyshev polynomial.

Thus, by using the Spectral–ROM approach to build the reduced-order model, the time-dependent coefficients $\{a_{\,i}\,(\,t\,)\}$ are computed by solving the following system:

[TABLE]

remembering that $\mathcal{A}\in\mathrm{Mat}\,_{s\times s}\,(\mathds{R})$ is a constant coefficient matrix, $\mathrm{b\,}(\,t\,)\in\mathds{R}^{s}$ is a vector coming from the boundary conditions and $a_{\,0}$ is the vector of initial spectral coefficients. The main advantage of a Spectral–ROM is that $s\ll p$, where $p$ is the number of degrees of freedom needed to solve problem (3.2) by means of conventional methods (finite-differences, finite-elements and finite-volumes). We note that the matrix $\mathcal{A}$ and the vector $\mathrm{b\,}(\,t\,)$ might depend on problem parameters, such as the diffusion coefficient $\nu\,$:

[TABLE]

Different approaches can be used to solve the system of ODEs (3.10), depending on the cases considered. The most straightforward way to use the Spectral–ROM from Eq. (3.10) is to apply a numerical integration scheme, e.g., an adaptive Runge–Kutta with moderate accuracy, since Eq. (3.10) is just a ROM. So, with an embedded error control and not so stringent tolerances, it can be done very efficiently. In this study, we shall employ ODE solvers for simplicity, since we are interested in the whole trajectory.

Validation of the numerical solution

To compare and validate the proposed method, the error between a solutions obtained by one of the numerical methods $u^{\,\mathrm{num}}\,(\,x,\,t\,)$, and the reference solution $u^{\,\mathrm{ref}}\,(\,x,\,t\,)$, is computed as a function of $x$ by the following formulation:

[TABLE]

where $N_{\,t}$ is the number of temporal steps. The global error $\varepsilon_{\,\infty}$ is given by the maximum value of $\varepsilon_{\,2}\,(\,x\,)$:

[TABLE]

The computation of the reference solution $u^{\,\mathrm{ref}}\,(\,x\,,t\,)$ is detailed in Sections 4.1, 4.2 and 5.1.

Numerical application

Linear case

The first case considers linear moisture transfer in a material with $0.1\,\mathsf{m}$ of length. The moisture transport coefficient has a value of $d_{\,m}\,=\,1.97\cdot 10^{-10}\,\mathsf{s}$ and the moisture storage a value of $c_{\,m}\,=\,7.09\cdot 10^{-3}\,\mathsf{kg/m^{3}/Pa}$ [17]. The initial vapour pressure across the material is considered to be uniform as $P_{\,v}^{\,i}\,=\,1.16\cdot 10^{\,3}\,\mathsf{Pa\,}$, corresponding to a relative humidity of $50\,\%$ and to a temperature of $20^{\,\circ}\mathsf{C}$. Simulations are performed for a total time of $120\,\mathsf{h}\,$. The boundary conditions, represented by the relative humidity $\phi$ are given in Figure 1. The sinusoidal variations oscillate between dry and moist states during the total simulation time. The convective vapour transfer coefficients are set to $h_{\,v,\,\mathrm{L}}\,=\,2\cdot 10^{-7}\,\mathsf{s/m}$ and $h_{\,v,\,\mathrm{R}}\,=\,3\cdot 10^{-8}\,\mathsf{s/m}$ for the left and right boundaries, respectively. As the readers may be interested in simulating the proposed case, dimensionless values are provided in Appendix A.

This case study is performed with the Spectral–ROM using $N\,=\,6$ modes and with two central finite-difference approximations schemes: (i) the Euler implicit and (ii) the Crank–Nicolson. The reference solution is computed using the Matlab open source toolbox Chebfun [15].

The reduced system of ODEs is implemented in Matlab and the spectral coefficients $\{a_{\,n}\,(\,t\,)\}$ are calculated for any intermediate time instant by the solver ODE45. The solver is set with an absolute and relative tolerance of $\mathsf{tol}\,=\,10^{\,-4}\,$. The integration in time is based on an explicit Runge–Kutta formula for ODE45. The inputs are the initial time, the final time and the time step (optional) and the solver supplies the integration at the given time. One should recall that computations of the Spectral solution are performed for the reference spatial domain of $[\,-1,\,1\,]$ and then transformed to the interested one.

It can be seen that the physical phenomena are well represented, as illustrated in Figure 2(a) with the evolution of the vapour pressure at $x\,=\,0.04\,\mathsf{m}$. The variations follow the conditions of the left boundary and with the diffusion process going towards the periodic regime. It can be noted a good agreement between the Spectral–ROM and the other methods. Furthermore, the vapour pressure profile is shown in Figure 2(b) for the instants $t\,=\,\{8,\,50,\,120\}\,\mathsf{h}$, enhancing the good accuracy of the solution to represent the physical phenomena.

The absolute error $\varepsilon_{\,2}$ of the numerical methods applied and the reference solution is of the order of $\mathcal{O}\,(\,10^{\,-4}\,)\,$, as illustrated in Figure 3. The solutions of the problem have been computed for discretisation parameters of $\Delta x^{\,\star}\,=\,1\cdot 10^{\,-2}$ and $\Delta t^{\,\star}\,=\,1\cdot 10^{\,-1}$ for the Spectral–ROM and the Crank–Nicolson methods. However, the Euler implicit scheme needed more refinement to reach the same order of accuracy, with $\Delta x^{\,\star}\,=\,1\cdot 10^{\,-2}\text{ and }\Delta t^{\,\star}\,=\,1\cdot 10^{\,-2}\,$.

Figure 4(b) presents the absolute error $\varepsilon_{\,2}$ for the Spectral–ROM using different number of modes. As we increase the number of modes, the solution of the Spectral–ROM gets more accurate with the solution converging within few modes (less than 10). To illustrate the convergence of the solution, the profile of the vapour pressure for the last time instant of simulation is presented for a different number of modes in Figure 4(a). In this case, if we compare the solution with $3$ modes to the solution with $5$ modes a significant difference can be noticed. With $5$ modes we already have a satisfactory solution to the problem, with the absolute error of the order of $\mathcal{O}\,(\,10^{-3}\,)\,$, while the solution with $3$ modes is still oscillating. The number of modes of the Spectral method is predetermined in order to build the system of ODEs. In this case, a number of six modes proved to be good enough.

Spectral coefficients $a_{\,n}\,(\,t\,)$ are shown in Figures 5(a) and 5(b). It can be seen the first coefficients have the most significant values. For this reason, the Spectral method needs few modes to converge to the solution (an order of $10$) because its first modes have the highest magnitudes. A brief comparison with an analytical solution, built on Fourier decomposition [35], reveals that the eigenvalues of the Spectral method decrease faster, as shown in Figure 6. Note that the eigenvalues of the analytical solution do not have to coincide with the ones of the Spectral method since the eigenfunctions are not the same for the Chebyshev polynomials and the trigonometric ones. Furthermore, the magnitude of the last spectral coefficient acts as an error estimator, determining the error upper limit.

The global absolute error $\varepsilon_{\,\infty}$ for the conventional numerical methods applied is calculated as a function of spatial discretisation $\Delta x^{\,\star}\,$. Fig. 7 shows that the Spectral-ROM has the same accuracy for all values of $\Delta x^{\,\star}$. It is due to the fact that the Spectral method is based on Chebyshev polynomials, which enables to calculate the solution in each spatial node, as an analytical solution. For this reason, the error of the Spectral solution is almost a straight line, not depending on the spatial discretisation. However, for the conventional methods, the solution gets inaccurate when the value of $\Delta x^{\,\star}$ increases. It should be noted that the Spectral–ROM can provide even more accurate results, by increasing the number of modes or by decreasing the tolerance in the ODE Matlab solver to certain limits.

Weakly nonlinear case

This case is called weakly nonlinear because the boundary conditions remain linear and, only the diffusion coefficient has a slight dependency on the moisture field. Thus, the diffusion equation is written as:

[TABLE]

where, $\nu\,(\,u\,)\ =\ \nu_{\,0}\ +\ \nu_{\,1}\cdot u\,$. Since we have the diffusion coefficient depending on the field $\nu\,(\,u\,)$, the diffusion equation — Eq. (4.1) — can be rewritten as:

[TABLE]

where the residual has the following form:

[TABLE]

Then, by applying the Tau–Galerkin method, the residual is minimized by assuming it orthogonal to the basis functions $\langle\,R\,,\mathsf{T}_{\,k}\,\rangle\,=\,0\,$, which is defined in Eq. (3.6), leading to the following equation:

[TABLE]

where,

[TABLE]

Eq. (4.2) is a closed system of ODEs. Coefficients $c_{\,i,\,j,\,k}$ are calculated at once, and coefficient $\tilde{\tilde{a}}_{\,i}$ are related to $a_{\,i}$ though a linear transformation $\tilde{\tilde{a}}\ =\ \textsc{D}_{\,2}\cdot a$, in which $\textsc{D}_{\,2}\in\mathrm{Mat}\,_{(n-2)\times(n-2)}(\mathds{R})$ is a second order derivative matrix.

4.2.1  Case study

This case considers that $c_{\,m}$ and $d_{\,m}$ have a slight dependency on the moisture. The material piece has a length of $0.1$ $\mathsf{m}$, with a relative humidity-dependent diffusion coefficient:

[TABLE]

The initial vapour pressure in the material is considered uniform $P_{\,v}^{\,i}\ =\ 1.16\cdot 10^{\,3}$ $\mathsf{Pa\,}$, corresponding to a relative humidity of $50\,\%$ and to temperature of $20^{\circ}\mathsf{C\,}$. Simulations are performed for a total time of $72\,\mathsf{h}$, the equivalent of three days. The boundary conditions, represented by the relative humidity $\phi$ are given in Figure 8. The relative humidity oscillates sinusoidally between $50\,\%$ and $75\,\%$ on the left boundary and between $50\,\%$ and $80\,\%$ on the right boundary. The convective vapour coefficients are set to $h_{\,v,\,\mathrm{L}}\,=\,3\cdot 10^{\,-8}\,\mathsf{s/m}$ and $h_{\,v,\,\mathrm{R}}\,=\,2\cdot 10^{\,-7}\,\mathsf{s/m}$ for the left and right boundaries, respectively. The dimensionless values of this case are also provided in Appendix A.

The Spectral reduced-order model is composed by $N\,=\,6$ modes and its coefficients $\{a_{\,n}(\,t\,)\}$ are obtained through the use of the solver ODE45, with a tolerance set to $\mathsf{tol}\,=\,10^{\,-4}$. The discretisations used to compute the Spectral solution are $\Delta x^{\,\star}\,=\,10^{\,-2}$ and $\Delta t^{\,\star}\,=\,10^{\,-1}\,$. The reference solution is computed with the open source toolbox Chebfun in Matlab.

The evolution of the vapour pressure in the middle of the material, at $x\,=\,0.05$ $\mathsf{m}$, is shown in Figure 9(a). The vapour pressure varies according to the sinusoidal fluctuations from both boundary conditions. The vapour pressure profiles at different times are illustrated in Figure 9(b) for $t\,=\,\{9,\,38,\,72\}\,\mathsf{h}$, highlighting the good agreement of the Spectral solution in representing the variations.

The absolute error $\varepsilon_{\,2}$ has been computed between the reference solution and the Spectral–ROM for a different number of modes, as illustrated in Figure 10. For $N\,=\,6$ and $N\,=\,5$ modes the absolute error is of the same order, $\mathcal{O}\,(\,10^{\,-3}\,)\,$, proving the accuracy of the solution and showing that $5$ modes are good enough.

Figures 11(a) and 11(b) present the first and the last three coefficients $a_{\,n}\,$, respectively. The magnitude of the coefficient, in the total contribution of the solution, decreases with the order of the coefficient. The last coefficient determines the magnitude of the error, implying that the error will be lower than $a_{\,n}\,$. It is due to the truncation in the number of terms in the spectral representation of the solution and the fact that the solution is smooth. Thus, the higher the number of modes, the higher the accuracy. For this case, we cannot have a more precise solution than \sup_{\,t\,\in\,\bigl{[}\,0\,,T\,\bigr{]}}|a_{\,6}\,(t)|\,=\,1.3\cdot 10^{\,-3}\,.

Numerical cost estimation

The number of operations for each approach can be estimated. We denote by $N_{\,x}$ and $N_{\,t}$ the number of nodes according to the discretisation in both space and time domains. For explicit methods, it can be related to the CFL type conditions. A standard approach based on the Euler implicit scheme requires $N_{\,x}\cdot N_{\,t}\,$, operations while the Crank–Nicolson scheme requires at least twice as many, as it is built on both implicit and explicit parts. Considering the same discretisation parameters $N_{\,t}\,=\,1200$ and $N_{\,x}\,=\,100$ for both methods, the number of operations for the linear case scales with:

[TABLE]

Considering these discretisation parameters, the order of accuracy is not the same for both methods, it is at the order of $\varepsilon_{\infty}\simeq\mathcal{O}\,(10^{-3})$ for the Euler implicit and of $\varepsilon_{\infty}\simeq\mathcal{O}\,(10^{-4})$ for the Crank–Nicolson solution.

For the Spectral-ROM, the number is related to the solution of the system of ODEs (Eq. (3.8)), computed in this case with the Matlab solver ODE45. It is based on the iterative Runge–Kutta method to approximate the solution. The number of operation depends on the tolerance (tol) of the solver, which has a maximum tolerance of $\sim 10^{\,-5}$ for ODE45. Thus, we have:

[TABLE]

where $\tau$ is the total time of simulation. At each time step, the Runge–Kutta needs to compute the vector product $\mathcal{A}_{\,s\times s}\,$, where $s$ depends on the degree of freedom $N$ of the solution ($s\,=\,N\,-\,2$). Thus, it leads to $6\cdot s^{\,2}$ operations to perform, knowing that $s$ is of the order of $10\,$. Consequently, the total number of operations for the Spectral-ROM scales with approximately:

[TABLE]

Considering the first case, knowing that the tolerance was set to $10^{\,-4}\,$, with $N\,=\,6$ modes the number of operations performed by the Spectral–ROM is expressed by:

[TABLE]

Comparing the number of operations of this case, we can already see that the Spectral–ROM is less costly than the other methods applied. Notice that the number of degrees of freedom necessary to solve the diffusion problem by means of the Spectral method is inferior to the ones necessary to solve the whole system of partial differential equations. Using Euler or Crank–Nicolson methods, the computational complexity scales with $p\,=\,N_{\,x}\,$, whereas the one Spectral–ROM is $s\,=\,(\,N\,-\,2\,)\,$. For this case, the numerical application gives $p\,\equiv\,10^{\,2}$ and $s\,\equiv\,4\,$. Moreover, we can note the reduction of the order of the solution, using the Spectral approach. According to the previous results, the fidelity of the model is not degraded but only the order of the solution. Besides, $N_{\,t}$ can not be reduced due to accuracy issues.

4.3.1  Solving the system of ODEs

The time spent on simulations is also related to the solver used to compute system of ordinary differential equations. For the weakly nonlinear case, different solvers were employed with different values of tolerances. Between all Matlab ODE solvers, the ODE15s was the most efficient, combining accuracy and speed. By decreasing the tolerance of the solver, more accurate results can be obtained. But sometimes, if we improve the precision of the solver, the error can be limited by the magnitude of the last spectral coefficient. Thus, another way to get more precise solutions is through increasing the number of modes.

Therefore, depending on the accuracy sought on the results, several options are available. The choice of the ODE solver is related to the nature of the problem. For example, if the problem has two components which vary drastically on different time scales, then the problem is stiff, or difficult to evaluate. The solvers are then classified according to the problem type. For non-stiff problems, ODE45, ODE23 and ODE113 are the most appropriate, but for stiff problems, the other ODE solvers are recommended (ODE15s, ODE23s, ODE23t and ODE23tb). Further information can be found in [44].

Treating general nonlinearities

Problem (2.5) has an important difficulty in dealing with the nonlinearities of the moisture storage coefficient $c_{\,m}$ and of the diffusion coefficient $d_{\,m}\,$, both depending on the moisture content field. These coefficients are usually given by empirical functions from experimental data. Due to those nonlinearities, some modifications in the way of using the Spectral method have to be taken into account. For this reason, Eq. (2.5a) is recalled with a simplified notation:

[TABLE]

In order to better apply the Spectral method, Eq. (5.1) is rearranged as follows:

[TABLE]

where,

[TABLE]

By using Spectral methods the unknown $u\,(\,x\,,t\,)$ is approximated by the finite sum (3.3) with Chebyshev polynomials as basis functions. The derivatives are written as in the linear case, by Eqs. (3.4a), (3.4b) and (3.4c). Thus, substituting them into Eq. (5.2) gives:

[TABLE]

By applying the Galerkin projection we have:

[TABLE]

where,

[TABLE]

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 determined according to numerical investigations and will be discussed for the next case study.

In addition, we have the expressions of the nonlinear boundary conditions:

[TABLE]

Different from the linear case, the boundary conditions cannot provide an explicit expression for the two last coefficients $a_{\,n}\,(\,t\,)$ and $a_{\,n-1}\,(\,t\,)\,$. Thus, it is not possible to compute the solution in the same way. Although, with all elements listed before, it is possible to set the system to be solved by composing a system of ODEs with two additional algebraic expressions for the boundary conditions. It results in a system of Differential–Algebraic Equations (DAEs) with the following form:

[TABLE]

where, $\mathcal{M}$ is a diagonal and singular matrix ($\mathrm{rank}\,(\,\mathcal{M}\,)\,=\,n\ -\ 2$) containing the coefficients of the Chebyshev weighted orthogonal system, $\mathrm{b\,}\,(\,t\,)$ is a vector containing the boundary conditions and, $\mathcal{A}\cdot a_{\,n}\,(\,t\,)$ is composed by the right member of Eq. (5.3). The initial condition is given by Eq. (3.9) and the DAE system is solved by ODE15s or ODE23t from Matlab.

A highly nonlinear case

This case study considers moisture dependent coefficients $c_{\,m}$ and $d_{\,m}$, illustrated in Figures 12(a) and 12(b). Their variations are similar to the load bearing material from [22]. The initial vapour pressure is uniform $P_{\,v}^{\,i}\,=\,1.16\times 10^{\,\,3}$ $\mathsf{Pa}$. No moisture flow is taken into account at the boundaries. The ambient vapour pressures at the boundaries are illustrated in Figure 13. At the left boundary, it has a fast drop until the saturation state and at the right boundary, it has a sinusoidal variation. The material is thus excited until the capillary state. The convective vapour transfer coefficients are set to $h_{\,v,\,\mathrm{L}}\,=\,2\cdot 10^{\,-7}$ $\mathsf{s/m}$ and $h_{\,v,\,\mathrm{R}}\,=\,3\cdot 10^{\,-8}$ $\mathsf{s/m}$ for the left and right boundary, respectively. The simulation is performed for $120\,\mathsf{h}\,$. As in the previous case study, the dimensionless values can be found in Appendix A.

The Spectral method is composed by $N\,=\,9$ modes with $m\,=\,7\,$. The ODE15s was used to solved System (5.5), with a tolerance set to $10^{\,-4}\,$. For this case, the Spectral method was compared to the Crank–Nicolson [17] and to a reference solution computed using the Chebfun Matlab toolbox [15]. All solutions have been computed with the following discretisation parameters: $\Delta t^{\,\star}\,=\,10^{\,-2}$ and $\Delta x^{\,\star}\,=\,10^{\,-2}\,$.

Vapour pressure variations in the boundaries are shown in Figure 14(a). The vapour pressure at $x\,=\,0.1\,\mathsf{m}$ slowly oscillates according to the right boundary condition. It also increases within the material according to the step imposed at the left boundary $x\,=\,0\,\mathsf{m}$. This increase can be also observed on three profiles of vapour pressure illustrated in Figure 14(b), in which the diffusion process is represented going from left to right.

All methods have demonstrated good agreement to represent the physical phenomenon. Again, the fidelity of the model does not deteriorate with the use of a Spectral approach. Results of the error $\varepsilon_{\,2}$ in function of $x$ are shown in Figure 15(a). The error of the Crank–Nicolson scheme is proportional to $\mathcal{O}\,(\,\Delta t^{\star\,2}\,)\,$. The Spectral method with $N\,=\,9$ modes is one order of magnitude more accurate and faster than the Crank–Nicolson method, even considering the same discretization parameters $\Delta t^{\,\star}$ and $\Delta x^{\,\star}\,$. Although, if we decrease the number of modes to $N\,=\,6$ and maintaining the same discretization parameters $\Delta t^{\,\star}$ and $\Delta x^{\,\star}\,$, we reach the same order of accuracy of the Crank–Nicolson method, as observed in Figure 15(b).

The solution of the Spectral methods becomes more accurate with the increase of the number of modes, as shown in Figure 15(b). With $6$ modes, we have satisfactory results, with the error of the order of $\mathcal{O}\,(\,10^{\,-3}\,)\,$. As we increase only the number of modes, without changing other parameters, the error begins to stabilize, and with $8$ and $9$ modes the error remains the same.

As already observed in the linear case, the Spectral method does not depend on the number of spatial points, but on the order of the ODE solver tolerance and also on the number of modes. For the nonlinear case, the error also depends on the truncation of the sum $\sum_{k\,=\,1}^{m}\,$. For this reason, the error $\varepsilon_{\,\infty}$ in function of $m$ is shown in Table 1. The optimal $m$ number is approximated by numerical experimentation, and as can be seen in the Table 1, the best value for $m$ is the one equivalent to the number of modes.

Figures 16(a) and 16(b) represent the first and last three coefficients $a_{\,n}$ of the Spectral–ROM solution. The step in the left boundary can be also seen in these figures for the first days, and after that, the values tend to stabilize. It is possible to see the reduction in the magnitude of the coefficient with the increase of the number of coefficients. As for the previous cases, the last coefficients are always the smallest ones.

A parametric study is performed in order to verify the computational cost of the proposed method. The discretisation parameters are set to $\Delta x^{\,\star}\,=\,10^{\,-2}$ and $\Delta t^{\,\star}\,=\,10^{\,-2}\,$, while the number of modes $N$ of the Spectral solution and the tolerance of the solver vary. Figure 17(a) presents the maximum absolute error $\varepsilon_{\,\infty}$ in function of the number of spectral modes. As we increase the number of modes, the solution gets more accurate. Although, after a certain number of modes, the solution converges to a minimum value, that is related to the tolerance of the ODE solver. The time to perform each spectral simulation is presented in Figure 17(b). For this numerical application, the CPU time has been evaluated using Matlab platform on a computer with Intel i7 CPU and 8GB of RAM. The computational effort to perform the simulation increases linearly with the number of modes. However, it remains extremely low. To better appreciate the computational cost of each approach, Table 2 provides the CPU time to compute the solution using the Crank–Nicolson scheme, the Chebfun toolbox for the same discretisation parameters. The Spectral solution has been computed with $N\,=\,9$ modes. It is preferable to focus on the ratio of computer run time rather than on absolute values, that is system-dependent. Even with an average number of sub-iterations is $\,N_{\,\mathrm{NL}}\,\,\simeq\,\mathcal{O}\,(1)$ of the Crank–Nicolson scheme, the Spectral method is substantially faster than the other methods. It represents only $1\,\%$ of the CPU time needed using the Crank–Nicolson approach.

Multilayer domain

In constructions, multiple layers are commonly found. The configuration assumed at the interface between materials follows the hydraulic continuity [13], which considers inter-penetration of both porous structure layers. Both materials are homogeneous and isotropic, and only moisture transfer is simulated, through a perfectly airtight structure. The hydraulic continuity establishes that there must be a continuous moisture flow through the interface and a continuous distributions of vapour content:

[TABLE]

where $x_{\,\text{int}}$ represents the location of the interface between materials and subscripts $1$ and $2$ stand for Material $1$ and Material $2\,$, respectively.

Adaptation of the reduced Spectral Method

The original spatial domain $\Omega_{\,x}\,=\,[\,0\,,\,L\,]$ is decomposed in two sub-domains $\Omega_{\,x,\,1}\,=\,[\,0\,,x_{\,\text{int}}\,]$ and $\Omega_{\,x,\,2}\,=\,[\,x_{\,\text{int}}\,,\,L\,]\,$, which represent each material surface. These sub-domains are linear transformed to the spectral domain $\bar{\Omega}_{\,x,\,1}\,=\,[\,-1\,,1\,]$ and $\bar{\Omega}_{\,x,\,2}\,=\,[\,-1\,,1\,]$ so they can fit within the interval of interest as illustrated in Figure 18. From this, the unknown $u\,(\,x,\,t\,)$ is then defined as:

[TABLE]

in which $u_{\,1}\,(x,\,t)$ is the solution defined over domain $\bar{\Omega}_{\,x,\,1}$ and $u_{\,2}\,(x,\,t)$ is the solution defined over domain $\bar{\Omega}_{\,x,\,2}\,$. Thus, $u_{\,1}$ and $u_{\,2}$ are written respectively as:

[TABLE]

which represent the solution for Material $1$ and Material $2\,$, respectively. Note that the Chebyshev polynomials are always the same and the transformations always occur with the temporal coefficients.

The condition at the interface between two materials states the continuity of the fields and the flows. It implies that the derivative of the field $u$ is not continuous at the interface between two materials. This important remark has to be taken into account in the construction of the Spectral reduced order model. Indeed, the domain is decomposed in sub-domains to maintain a smooth solution and particularly a continuous derivative on each sub-domain. In this way, the model order reduction is optimal and ensure the error of the Spectral-ROM to decrease exponentially. It is totally 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. More modes would be necessary to reach the same accuracy, as detailed in Theorem $1$ of [48, Chap. 4].

By considering the two materials, Eq. (5.3) becomes:

[TABLE]

The interface conditions — Eqs. (6.1a) and (6.1b) — are adimensionalized and written in the spectral form as:

[TABLE]

which are included in vector $b_{\,1}$ 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_{\,2}$ and set equal to zero. Vectors $b_{\,1}$ and $b_{\,2}$ are column vectors of size $N\times 1$ with the form:

[TABLE]

With all elements listed before, it is possible to set the system to be solved. Different from the previous case, here the system of ODEs has the double of the size $2\,N$ and it has four additional algebraic expressions for the boundary and interface conditions. The initial condition is also given by Eq. (3.9) and the DAE system is solved by ODE15s from Matlab. In this work, the approach was presented for a wall with two layers for the sake of clarity, knowing that it can be extended to any number of layers.

A multilayer case

This case study considers a porous wall formed by $2$ layers: $10\ \mathsf{cm}$ of a load bearing material and $2\ \mathsf{cm}$ of a finishing material, as illustrated in Figure 19. The selected materials complicate the case, with the first layer having a faster liquid transfer while the second layer acts as an hygroscopic finish. The properties used for these materials were obtained from [20] and are presented in Figures 20(a) and 20(b). Temperature dependence was neglected and transport coefficients modeled as a function of moisture content. Boundary and initial conditions are set with the same values as in the previous case study: initial vapour pressure of $P_{\,v}^{\,i}\,=\,1.16\cdot 10^{\,3}\ \mathsf{Pa}$ on both materials and boundary conditions represented in Figure 13.

By using the Spectral approach, it is assumed that at the interface, the solution and the flow of the problem are continuous. In this way, the method will search for the solution that can satisfy both conditions. Simulations were performed using the ODE15s, with a tolerance set to $\mathsf{tol}\,=\,10^{\,-4}$ and with $N\,=\,10$ modes. These values were chosen based on the previous numerical study. The time is incremented with a discretization of $\Delta t^{\,\star}\,=\,10^{\,-2}\,$, which is equivalent to $36\ \mathsf{s}\,$.

The evolution of vapour pressure at the boundary surfaces ($x\,=\,0\ \mathsf{m}$ and $x\,=\,0.12\ \mathsf{m}$) is shown in Figure 21(a). At $x\,=\,0\ \mathsf{m}$ , the vapour pressure suddenly increases due to the step imposed at the surface. The moisture from the vapour pressure step diffuses through both layers. Although, as the second layer is composed with a less hygroscopic material, the vapour pressure completely reaches this surface by $65\ \mathsf{h}\,$. At $x\,=\,0.12\ \mathsf{m}\,$, the vapour pressure varies according to the sinusoidal fluctuations of the boundary conditions until the flow arrives. This increasing can also be observed on three profiles of vapour pressure illustrated in Figure 21(b). Different from the previous case, the moisture flow takes more time to reach the right boundary due to the material properties of the second layer. Finally, at $t\,=\,120\ \mathsf{h}$, it is still possible to observe the influence of the step on the vapour pressure.

The Spectral–ROM has demonstrated a good agreement with the reference to represent the moisture diffusion trough composed walls. Distribution of the error $\varepsilon_{\,2}$ on function of $x$ is given in Figure 22. The order of the error is the same as in the previous case $\varepsilon_{\,2}\simeq\mathcal{O}(10^{\,-4})\,$, but the error $\varepsilon_{\,\infty}$ is higher here. This is explained since we keep the same numerical configurations of the other case but the nonlinearities increase compared to the previous configuration. Nonetheless, results provided by the Spectral reduced-order model are still acceptable.

Regarding the CPU time, this case has also presented competitive outputs. The way in which the spatial domain was split, makes the reduced system to have the double of the size, if compared with a single layer simulation. Now, the matrix $\mathcal{A}$ has the double of its size, making the computer run time twice as high $(7\ \mathsf{s}\,)$, as shown in Figure 17(b).

Conclusions

Most of the numerical methods applied to mathematical models used in building physics are commonly based on implicit schemes to compute the solution of diffusive problems. Its main advantage is due to the stability conditions for the choice of the time step $\Delta t^{\,\star}\,$. However, implicit schemes require important sub-iterations when treating nonlinear problems. This work was therefore devoted to exploring the use of an innovative reduced-order approach based on the Spectral method. Spectral methods are well-known in other applications, such as meteorology and wave propagation, although it was not used before as a reduced-order model. Thus, in this work, we showed that they can be applied in some one-dimensional building physics problems to compute a reduced-order model.

The first case study considered a linear diffusive moisture transfer through a porous material. The Spectral–ROM was compared to the classical Euler implicit scheme, to the Crank–Nicolson scheme and to a reference solution obtained using Chebfun toolbox for Matlab. Results have shown the dynamics and amplitude of hygrothermal fields are perfectly represented by the Spectral-ROM solution. The fidelity of the physical model is totally conserved by the Spectral-ROM. Only the order of the solution is highly reduced. Using standard approaches, the order of the solution rises with $10^{\,2}$ whereas with the Spectral method, the order of the solution scales with $6\,$. In the second case, a weak nonlinear problem was treated, which has a field dependent diffusion coefficient. To build the reduced system of ODEs, the same features of the linear case were used. Its reduced system was written with an explicit formulation and then implemented in Matlab. In the highly nonlinear case, the reduced system is numerically obtained as the system of ODEs cannot be explicitly expressed. The third case study focused on such general highly nonlinear transfer model, with material properties strongly dependent on the relative humidity field. To treat the nonlinearities, the Chebyshev–Gauß quadrature was employed to solve the integrals. Again, the accuracy of the approach has been demonstrated by representing accurately the physical phenomenon, with an absolute error of the order of $\varepsilon_{\,2}\simeq\mathcal{O}\,(10^{\,-4})$ comparing to the reference solution. A parametric study on the number of modes and the tolerance of the ODE solver has also been carried out. Moreover, when comparing the CPU time of the different approaches, the Crank–Nicolson is one hundred times longer than the Spectral method to compute the solution. To bring applications closer to building physics problems, a wall with two materials is used for the last case study. By using the Spectral reduced-order model the spatial domain is decomposed and the interface conditions can be easily imposed. As the complexity of the problems rises, the Spectral method needs more modes, with still a very low computational effort compared to standard approaches, and yet it does not mean that the Spectral method loses its efficiency.

The application of Spectral methods is not straightforward, neither intuitive as for example for the finite-difference method. Although, the efforts used in its implementation are compensated by the results, which showed to be very promising. In other domains, Spectral methods have also demonstrated their great potential for solving more complex problems [28, 36, 25], which instigates the development of further work in the building physics field on the solution of combined heat and moisture transfer and through multidimensional geometries.

Acknowledgments

The authors acknowledge the Brazilian Agencies CAPES of the Ministry of Education and CNPQ of the Ministry of Science, Technology, and Innovation, for the financial support. The authors also acknowledge the Junior Chair Research program “Building performance assessment, evaluation and enhancement” from the University of Savoie Mont Blanc in collaboration with The French Atomic and Alternative Energy Center (CEA) and Scientific and Technical Center for Buildings (CSTB).

Nomenclature

[TABLE]

Appendix A Dimensionless values

Case from Section 4.1

Problem (2.5) is considered with $g_{\,l,\,\mathrm{L}}^{\,\star}\,=\,g_{\,l,\,\mathrm{R}}^{\,\star}\ =\ 0$ and the dimensionless properties of the material are equal to $d_{\,m}^{\,\star}\,=\,1$ and $c_{\,m}^{\,\star}\,=\,8.6\,$. The reference time is $t^{\,0}\,=\,1$ $\mathsf{h}$, thus the final simulation time is fixed to $\tau^{\,\star}\,=\,120\,$. The Biot numbers are $\mathrm{Bi}_{\,v,\,\mathrm{L}}\,=\,101.5$ and $\mathrm{Bi}_{\,v,\,\mathrm{R}}\,=\,15.2\,$. The boundary conditions are expressed as:

[TABLE]

Case from Section 4.2

Simplification of the problem (2.5) are carried out, considering $g_{\,l,\,\mathrm{L}}^{\,\star}\,=\,g_{\,l,\,\mathrm{R}}^{\,\star}\,=\,0\,$ and $\nu\,(\,u\,)\,=\,d_{\,m}^{\,\star}\,(\,u\,)/c_{\,m}^{\,\star}\,(\,u\,)$. In this way, the problem is written as:

[TABLE]

The reference time is $t^{\,0}\,=\,1$ $\mathsf{h}\,$, thus the final simulation time is fixed to $\tau^{\,\star}\,=\,72\,$. The Biot numbers are $\mathrm{Bi}_{\,v,\,\mathrm{L}}\,=\,15.2$ and $\mathrm{Bi}_{\,v,\,\mathrm{R}}\,=\,101.5\,$. The boundary conditions are expressed as:

[TABLE]

and, the dimensionless property of the material is:

[TABLE]

Case from Section 5.1

Problem (2.5) is considered with $g_{\,l,\,\mathrm{L}}^{\,\star}\,=\,g_{\,l,\,\mathrm{R}}^{\,\star}\,=\,0\,$. In this way, the dimensionless governing equations are written as:

[TABLE]

in which, the dimensionless properties of the material are:

[TABLE]

Simulations are performed for a total time of $\tau^{\,\star}\,=\,120\,$. The ambient water vapour pressure at the boundaries are different from the previous case study. At the left boundary, $u_{\,\mathrm{L}}$ has a fast jump until the saturation state u_{\,\mathrm{L}}\,=\,2,\;\forall t\in\,\bigr{[}\,10\,,40\,\bigl{]} and at the right boundary, $u_{\,\mathrm{R}}\,(\,t^{\,\star}\,)\,=\,1\ +\ 0.8\;\sin\left(\,\dfrac{2\pi\,t^{\,\star}}{4}\,\right)\,$, with $\mathrm{Bi}_{\,v,\,\mathrm{L}}\,=\,101.5$ and $\mathrm{Bi}_{\,v,\,\mathrm{R}}\,=\,15.2\,$. Reference values are $d_{\,m}^{\,0}\,=\,1.98\cdot 10^{\,-10}\ [\mathsf{s}]\,,$ $t^{\,0}\,=\,3600\ [\mathsf{s}]$ and $L\,=\,0.1\ [\mathsf{m}]\,$.

Case from Section 6.2

The dimensionless properties of Material 1 are:

[TABLE]

and of Material 2 are:

[TABLE]

with $u\,=\,[0,\,2]$. The ambient water vapour pressure at the boundaries are the same from the previous case study, with $\mathrm{Bi}_{\,v,\,\mathrm{L}}\,=\,4.4$ and $\mathrm{Bi}_{\,v,\,\mathrm{R}}\,=\,0.65\,.$ Reference values are $d_{\,m}^{\,0}\,=\,5.5\cdot 10^{\,-9}\ [\mathsf{s}]\,,$ $t^{\,0}\,=\,3600\ [\mathsf{s}]$ and $L\,=\,0.12\ [\mathsf{m}]\,$.