Advanced reduced-order models for moisture diffusion in porous media

TL;DR

This paper compares two model reduction techniques, Spectral-ROM and PGD, demonstrating their effectiveness in accurately and efficiently modeling moisture diffusion in porous media with significantly reduced computational cost.

Contribution

The article provides a detailed methodology for applying Spectral-ROM and PGD to moisture diffusion problems, serving as a numerical benchmark for future reduced-order modeling in porous materials.

Findings

Both methods achieve around tenfold reduction in model order.

The approaches provide solutions with less than 1% error compared to reference models.

They effectively handle linear, non-linear, and parametric moisture diffusion problems.

Abstract

It is of great concern to produce numerically efficient methods for moisture diffusion through porous media, capable of accurately calculate moisture distribution with a reduced computational effort. In this way, model reduction methods are promising approaches to bring a solution to this issue since they do not degrade the physical model and provide a significant reduction of computational cost. Therefore, this article explores in details the capabilities of two model-reduction techniques - the Spectral Reduced-Order Model (Spectral-ROM) and the Proper Generalised Decomposition (PGD) - to numerically solve moisture diffusive transfer through porous materials. Both approaches are applied to three different problems to provide clear examples of the construction and use of these reduced-order models. The methodology of both approaches is explained extensively so that the article can be used as a numerical benchmark by anyone interested in building a reduced-order model for diffusion problems in porous materials. Linear and non-linear unsteady behaviors of unidimensional moisture diffusion are investigated. The last case focuses on solving a parametric problem in which the solution depends on space, time and the diffusivity properties. Results have highlighted that both methods provide accurate solutions and enable to reduce significantly the order of the model around ten times lower than the large original model. It also allows an efficient computation of the physical phenomena with an error lower than 10^{-2} when compared to a reference solution.

Tables

Table 1. Table 1. Computational cost of the methods for the nonlinear case.

Method	CPU time ($𝗌$)	CPU time ($\%$)
Euler implicit	36	100
PGD using SVD	5.29	15
PGD using DEIM	1.9	5
Spectral–ROM	0.35	0.9

Equations

$$\frac{\partial U}{\partial t}$$

$$c_b{\rho }_0\frac{\partial T}{\partial t}$$

$$\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}=\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 P_c}{\partial x}=\frac{R_{v}T{\rho }_l}{P_v}\cdot \frac{\partial P_v}{\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}\simeq \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}=f^{\prime }(\varphi )\frac{1}{P_s}\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{)}\,\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$$

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

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

$${\mathrm{Bi}}_{v,\mathrm{R}}=\frac{h_{v,\mathrm{R}}\cdot L}{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$$

$$\displaystyle\frac{\partial u}{\partial t}\ =\ \nu\,\frac{\partial}{\partial x}\biggl{(}\,\frac{\partial u}{\partial x}\,\biggr{)}\,,$$

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

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

$$u(x,t)\approx u_n(x,t)=\sum _{i=0}^n{a}_i(t){\varphi }_i(x).$$

$${\varphi }_i(x)\equiv T_i(x).$$

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

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

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

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

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

$$⟨f,g⟩={\int }_{-1}^1\frac{f(x)g(x)}{\sqrt{1-x^2}}\mathrm{d}x.$$

$${\dot{a}}_i(t)-\nu {\widetilde{\stackrel{~}{a}}}_i(t)=0,$$

$$\sum _{i=0}^n{\tilde{a}}_i(t)T_i(-1)-{\mathrm{Bi}}_{v,\mathrm{L}}\sum _{i=0}^n{a}_i(t)T_i(-1)+{\mathrm{Bi}}_{v,\mathrm{L}}{u}_{\mathrm{L}}(t)$$

$$-\sum _{i=0}^n{\tilde{a}}_i(t)T_i(1)-{\mathrm{Bi}}_{v,\mathrm{R}}\sum _{i=0}^n{a}_i(t)T_i(1)+{\mathrm{Bi}}_{v,\mathrm{R}}{u}_{\mathrm{R}}(t)$$

$${\dot{a}}_i(t)=\mathcal{A}{a}_i(t)+\mathrm{b}(t),$$

$$a_{0,i}\equiv a_i(0)=\frac{2}{\pi c_i}{\int }_{-1}^1\frac{u_0(x)T_i(x)}{\sqrt{1-x^2}}\mathrm{d}x,$$

Full text

Key words and phrases: reduced-order modelling; moisture diffusion; numerical methods; spectral methods; Proper Generalised Decomposition (PGD)

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Moisture transfer in porous materials

3. 3 Methodology

4. 3.1 Reduced Spectral method

5. 3.1.1 Application of the Spectral–ROM

6. 3.2 Proper Generalised Decomposition

7. 3.2.1 Iterative resolution

8. 3.2.2 Computing R(t) and S(x)

9. 3.2.3 Convergence of the global enrichment

10. 3.3 Comparison of the numerical solution

11. 4 Linear transfer in porous material

12. 5 Computing parametric solution using reduced-order models

13. 5.1 Extension of model reduction techniques to parametric problems

14. 5.1.1 Spectral reduced-order model

15. 5.1.2 PGD reduced-order model

16. 5.2 Case study

17. 6 Nonlinear transfer in porous material

18. 6.1 Extension of model reduction techniques to nonlinear problems

19. 6.1.1 Spectral reduced-order model

20. 6.1.2 PGD reduced-order model

21. 6.2 Case study

22. 7 Conclusion

23. A Dimensionless values

24. A.1 Linear case

25. A.2 Parametric case

26. A.3 Nonlinear case

27. B Solving the Spectral ODE reduced system

Introduction

Mathematical models representing the physical phenomena of heat and moisture transfer in porous media have been developed since the works of Philip and De Vries [58] and Luikov [45]. According to Luikov [45, Chapter 6], the following system of differential equations represents the physical phenomenon of heat and mass transfer through capillary porous materials:

[TABLE]

where $U$ is the specific mass content in the porous body, corresponding to the ratio of the mass of the water (in the vapor and the liquid phase) to the mass of dry-basis body. The other quantities are $T\,$, the temperature, $a_{\,m}\,$, the mass transfer coefficient for vapor (denoted with the subscript $1$) and liquid inside the body, $\delta\,$, the thermal-gradient coefficient, $\rho_{\,0}\,$, the specific mass of the dry body, $c_{\,b}\,$, the specific heat of the body and, $r_{\,12}\,$, the latent heat of vaporization.

Particularly, in the context of building components, ones aims at estimating accurately the moisture transfer since it is related to the occupants’ health, to the durability of building porous elements and to the building energy consumption and demand. However, the elaboration of numerical methods for the mathematical problem (1.1) faces several issues. First, the physical phenomena of moisture and heat transfer do not evolve on the same time scale. The characteristic time may scales, for many materials, from minutes for the heat diffusion and to a month for the moisture diffusion. Moreover, the time horizon of the studies is of the order of $1\ \mathsf{year}$ when evaluating the building energy consumption. The characteristic times of the physical phenomena have differences of several orders of magnitude. Another complexity of the problem arises from the material properties, which strongly vary with moisture and temperature, with space when considering multi-layers porous walls and sometimes with times when looking at the durability of the materials. Therefore, these complexities challenge the elaboration of numerical models.

Several works have been presented in the literature to propose numerical methods for some mathematical model (1.1). A list of the main tools can be obtained in [47, 68]. Most of the hygrothermal simulation tools are built using numerical approaches and discrete representations of the continuous equations by means of standard and incremental techniques (Finite-Difference, Finite-Volume and Finite-Element methods) to compute the solution. Due to stability conditions, most of the approaches are based on implicit schemes as described for instance in [48, 39]. Therefore, these approaches do not circumvent the complexities of the mathematical problem. These schemes require the solution of large systems of equations (an order of $10^{\,6}$ for three-dimensional problems). Moreover, when considering nonlinear building material properties, sub-iterations at each time step are induced, increasing significantly the computational cost as mentioned for examples in [29, 28]. It is important to mention that the computational resources available in a computer are not increasing anymore [66]. Thus, it is worth investigations to develop numerical models based on an optimal usage of the available computational resources.

Therefore, the challenge relies on proposing efficient numerical models, capable of computing accurately the solution with a reduced computational effort. These models are of great concern particularly for investigations requiring important computations, as for instance sensitivity analysis or parameter estimation problems. In this way, model reduction techniques appear as a very interesting alternative to substantially reduce the number of operations, saving computational resources (CPU time and memory) with no loss of accuracy of the solution. It is important to note that the fidelity of the mathematical model is conserved. In recent years, reduced-order modeling techniques have proven to be powerful tools, providing accurate predictions while dramatically reducing computational time, for a wide range of applications, covering from fluid mechanics, heat transfer, structural dynamics and other fields.

Therefore, the scope of this work is limited to the application of two advanced a priori reduction-order model techniques, the PGD [13] and the Spectral–ROM [32] techniques, which have been successfully implemented for coupled heat and moisture diffusion problems. The comparison is carried out for moisture transfer phenomena providing an important primary evaluation for anyone intending to build a reduced-order model for diffusion problem. Several features of the methods are analyzed in terms of both model order reduction and accuracy of the computed solution. The investigation is carried for three case studies: (i) linear transfer; (ii) parametric problems, aiming at computing a model whose solution depends on the material properties and (iii) nonlinear transfer with moisture dependent material properties.

Finally, the manuscript is organized as follows: first, Section 2 presents the description of the physical phenomena; then Section 3 gives an explanation of the PGD and Spectral reduced-order model techniques; and, further sections present the three different case studies: (i) linear transfer in Section 4; (ii) parametric problem whose depends on time, space and on the diffusion parameter in Section 5, and (iii) nonlinear transfer in Section 6.

Moisture transfer in porous materials

The physical problem involves unidimensional moisture diffusion through a porous material defined by the spatial domain $\Omega_{\,x}\,=\,[\,0\,,L\,]$. The moisture transfer only occurs according to the liquid and vapor diffusion. The physical problem can be formulated as [38, 2, 63]:

[TABLE]

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

Eq. (2.1) can be written using the vapor 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]

Thus, we have:

[TABLE]

As we consider the mass transfer under isothermal conditions, the second term vanishes and we obtain:

[TABLE]

In addition, we have:

[TABLE]

Considering the relation $\rho_{\,l+v}\ =\ f\,(\phi)\ =\ f\,(\,P_{\,v}\,,T\,)\,$, obtained from material properties and from the relation 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 vapor 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}}$ the moisture storage coefficient.

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

[TABLE]

where $P_{\,v,\,\mathrm{L}}$ and $P_{\,v,\,\mathrm{R}}$ are the vapor pressure of the ambient air, $g_{\,l,\,\mathrm{L}}$ and $g_{\,l,\,\mathrm{R}}$ are the liquid flow (driving rain) at the two bounding surfaces. The initial condition is consider with a uniform vapor pressure distribution:

[TABLE]

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 [33]. Therefore, we define the following dimensionless parameters:

[TABLE]

In this way, the dimensionless governing equations are then written as:

[TABLE]

Finally, this is the problem of interest considered in this work. The procedure of the methods used for the problem solution is described in the next section.

Methodology

For the sake of simplicity, without losing the generality, the methods are first explained considering $d_{\,m}^{\,\star}$ and $c_{\,m}^{\,\star}$ as constants, noting $\nu\mathop{\ \stackrel{{\scriptstyle\mathrm{def}}}{{:=}}\ }\dfrac{d_{\,m}^{\,\star}}{c_{\,m}^{\,\star}}\,$, and we drop $\star$ for the sake of notation compactness. Thus, considering a linear diffusion equation:

[TABLE]

for x\,\in\,\big{[}\,0,\,1\,\big{]}. The boundary conditions are:

[TABLE]

Using a standard discretization method, such as Euler or Crank–Nicolson, to compute the solution of Eq. (3.1) yields in computing a solution $u\,(\,x\,,t\,)$ for each point of the discretised spatial and time domains. The following is adopted: $N_{\,x}$ and $N_{\,t}\,$, which stand for the number of elements according to the discretization of the space and time domains. Thus, the order of the so-called large original model is $p\ =\ N_{\,x}\cdot N_{\,t}\,$.

Model reduction aims at decreasing the degrees of freedom present in a numerical model. It aims at approximating the solution by a lower order model $N\,\ll\,p$ without reducing drastically the fidelity of the physical model. One of the features is to significantly decrease the computational resources space (CPU time and memory).

The scope of this work is limited to the two a priori reduction-order model techniques, the PGD [13] and the Spectral–ROM [32] techniques, which have been successfully implemented in building physics to reduce computational cost while maintaining high fidelity solutions. Both techniques assume separated tensorial representation of the solution by a finite sum of function products. The Spectral–ROM 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, while the PGD aims at computing directly the basis of functions by minimizing the residual. Further details on both methodologies are given in the next sections.

Reduced Spectral method

Spectral methods are successfully applied in studies of wave propagation, meteorology, computational fluid dynamics, quantum mechanics and several other fields [19]. Some works on the transport phenomena can be found in literature involving diffusive [36, 67], convective [22, 61] and radiative [44, 21, 46] heat transfer. Spectral techniques applied in these works are varied, adopted according to the geometry, boundary conditions, field of application, and other necessities. In recent works, researchers have implemented spectral methods for solving heat and moisture transfer in food engineering [56] and on fluid flows [50]. According to the authors’s knowledge, there are no results in the literature regarding the application of spectral methods for solving diffusive moisture transfer in building physics applications other than [32].

Spectral methods consider a global representation of the solution, which means the derivative at a certain spatial point depends on the solution of the entire domain and not only on its neighbors [18]. Besides, Spectral methods consider a sum of polynomials that suit for the whole domain, almost like an analytical solution, providing a high approximation of the solution. Therefore, as its error decreases exponentially, it is possible to have the same accuracy of other methods but with a lower number of Galerkin modes, which makes this method memory minimizing, allowing to store and operate a lower number of degrees of freedom [65].

3.1.1  Application of the Spectral–ROM

The idea of the Spectral method is to assume that the unknown $u\,(\,x\,,t\,)$ from Eq. (3.1) can be approximatively represented as a finite sum [47]:

[TABLE]

Here, $\{\phi_{\,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 [32]. The Chebyshev polynomials are chosen as the basis functions as they are optimal in $\mathcal{L}_{\,\infty}$ approximation norm [34]. Therefore, we have:

[TABLE]

For more details on Chebyshev polynomials the readers may refer to [57, 18].

As we have chosen the basis functions, we can compute 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’s 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 [57].

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

[TABLE]

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

[TABLE]

which is solved via the Tau–Galerkin method that requires Eq. (3.5) to be orthogonal to the Chebyshev basis functions $\langle\,R\,,T_{\,i}\,\rangle\,=\,0$. The scalar product is defined by:

[TABLE]

Thus, it leads to the following relation between the spectral coefficients:

[TABLE]

Finally, after the projection and expansion of the residual, the result is a system of Ordinary Differential Equations (ODE), 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). 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. Therefore, two more equations are added to the system:

[TABLE]

with $T_{\,i}\,(-1)\,=\,(-1)^{\,i}$ and $T_{\,i}\,(\,1\,)\,\equiv\,1$ (see [57]). Eqs. (3.6a) and (3.6b) are written in an explicit way, with coefficients $a_{\,n}$ and $a_{\,n-1}$ expressed in terms of 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 can be explicitly built. Otherwise, we have a system of Differential-Algebraic Equations. Moreover, the reduced system of ordinary differential equations has the following form:

[TABLE]

where, $\mathrm{b\,}(\,t\,)\in\mathds{R}^{\,N}$ is a vector coming usually from boundary conditions and $\mathcal{A}\in\mathrm{Mat}\,_{N\times N}(\mathds{R})\,$ is a matrix with constant coefficients and with $\mathcal{O}\,(\,N\,)\simeq 10\,$. The main advantage of a Spectral–ROM is that $N\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).

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

[TABLE]

where, $u_{\,0}\,(\,x\,)$, is the dimensionless initial condition. After solving the reduced system of ODEs (Eqs. (3.7) and (3.8)), 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]

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.9), depending on the cases considered. The most straightforward way to solve it is to apply a numerical integration scheme, e.g., an adaptive Runge–Kutta with moderate accuracy. So, with an embedded error control and not so stringent tolerances, it can be done very efficiently. A detailed presentation on how to solve this system can be found in Appendix B.

Proper Generalised Decomposition

The Proper Generalised Decomposition (PGD) is also a model reduction method. It originates in the radial space-time separated representation proposed by Ladevèze in $1985$ [42]. In $2006\,$, the separated representations were extended to the multidimensional case by Chinesta and co-workers [8]. Interested readers may see [25, 26] for additional details on the method as well as [24] for an introduction. This strategy has been successfully applied and validated for various industrial applications. For instance, the PGD method was applied to quantum mechanics (Schrödinger equation) [7], kinetics theory (Fokker–Planck equation non-Newtonian fluids) [59, 9], phase separation in heterogeneous mixtures (Langer equation) [43], virtual surgery (forces, vibrations, etc.) [52], nonlinear stochastic problems (Burgers equation, 2D nonlinear diffusion problems) [53], multi-scale and multiphysics problems (visco plasticity, damage, etc.) [51, 16], computational fluid dynamics (anisotherm Navier–Stokes problems) [31], and, more recently to, heat and moisture transfer in building materials [15].

The PGD solution to problem Eq. (3.1) is sought as a separated representation of functions of time $t$ and space $x\,$:

[TABLE]

The order of PGD ROM scales with \mathcal{O}\,(\,p\,)\ =\ M\cdot\bigl{(}\,N_{\,x}\ +\ N_{\,t}\,\bigr{)}.

3.2.1  Iterative resolution

Solving problem (3.1) numerically using the PGD method consists in calculating modes $(F^{\,i}\,,G^{\,i})$ iteratively from $i\ =\ 1$ to $M\,$. The first mode $(F^{\,1}\,,G^{\,1})$ is initialised in order to satisfy the initial and boundary conditions in all zones. At enrichment step $m\,<\,M$, we assume that a former approximation of $u\,(\,x\,,t\,)$ is known and the new couple $F^{\,m+1}\,(\,x\,)\ =\ R\,(\,x\,)$ and $G^{\,m+1}\,(\,t\,)\ =\ S\,(\,t\,)$ has to be calculated according to:

[TABLE]

Eq. (3.11) is introduced into Eq. (3.1). Thanks to the separated representation of the solution $u\,(\,x\,,t\,)$ for dimensions $t$ and $x$ (Eq. (3.10)), we get:

[TABLE]

where $\mathrm{Res}^{\,m+1}$ is the residual of Eq. (3.1) because Eq. (3.11) is an approximation of the solution.

3.2.2  Computing R(t) and S(x)

We are at step $m$ searching for the new couple $R$ and $S$ by solving Eq. (3.12). To compute them, Eq. (3.12) will be successively projected on $R$ and $S\,$. For this, we note the scalar product $\langle\,\bullet,\,\bullet\rangle_{\,y}$ in the domain $\Omega_{\,y}\,$, defined by:

[TABLE]

evaluated numerically using the discrete values of functions $f$ and $g$ and a trapezoidal approximation for the integral. Here, the scalar product is defined for both time and space domains.

Eq. (3.12) is projected on $R\,$, assuming $\langle\,\mathrm{Res}^{\,m+1},R\,\rangle_{\,x}\ =\ 0$ to obtain:

[TABLE]

with:

[TABLE]

Eq. (3.12) is now projected on $S\,$, assuming $\langle\,\mathrm{Res}^{\,m+1}\,,S\rangle_{\,t}\ =\ 0$ and gives:

[TABLE]

with:

[TABLE]

After theses projections, to solve Eqs. (3.13) and (3.14), an alternating direction fixed-point algorithm is used. The stopping criterion, assuming the algorithm has converged, is:

[TABLE]

where $q$ is the index of iteration of the fixed-point algorithm and $\eta_{\,1}$ is a tolerance parameter chosen by the user.

3.2.3  Convergence of the global enrichment

Functions $R$ and $S$ have just been computed by a fixed-point algorithm. The PGD basis is enriched, noting $F^{\,m+1}\ =\ R$ and $G^{\,m+1}\ =\ S$ the new modes. The field of interest $u$ can be written as:

[TABLE]

The enrichment of the PGD solution stops when the norm of the residual \bigl{|}\bigl{|}\,\mathrm{Res}^{\,m+1}\,\bigr{|}\bigr{|} is assumed negligible with respect to $\eta_{\,2}\,$, another tolerance parameter chosen by the user:

[TABLE]

Comparison of the numerical solution

To compare and validate the proposed methods, the error between the solution $u^{\,\mathrm{num}}\,(\,x\,,t\,)$, obtained by one of numerical the methods, 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 uniform 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 further Sections.

Linear transfer in porous material

The first case of linear moisture transfer is considered from [62, 12] to analyze the moisture buffer effects in a $500$-$\mathsf{mm}$ aerated concrete under isothermal condition, at a temperature of $23\,\,\mathsf{{}^{\circ}C}\,$. The vapor permeability is $d_{\,m}\,=\,3\cdot 10^{\,-11}\,\,\mathsf{s}\,$ and its moisture storage is $c_{\,m}\,=\,1.85\cdot 10^{\,-4}\,\,\mathsf{kg/m^{3}/Pa}\,$ [12]. The uniform initial vapor pressure in the material is $P_{\,v}^{\,i}\,=\,842\,\,\mathsf{Pa\,}\,$, corresponding to a relative humidity of $30\,\%\,$. The total time of simulation corresponds to $120\,\,\mathsf{h}\,$. The left boundary is set to a constant vapor pressure, identical to the initial condition. At the right boundary, the relative humidity varies sinusoidally between $33\,\%$ and $75\,\%\,$, with a period of $24\,\,\mathsf{h}\,$. The convective vapor transfer coefficient is set to $2\cdot 10^{\,-8}\,\,\mathsf{s/m}\,$.

The solution of the problem has been computed for a discretization of $\Delta x^{\,\star}\ =\ 5\cdot 10^{\,-3}$ and $\Delta t^{\,\star}\ =\ 10^{\,-1}\,$. The physical phenomena are well represented, as illustrated in Figure 1(a) with the time evolution of the relative humidity at $x\,=\,0.47\,\mathsf{m}$. The variations follow the ones of the right boundary conditions and, the diffusion process goes towards the periodic regime. It can be noted a good agreement between the two reduced-order models. Furthermore, the vapor pressure profile is shown in Figure 1(b) for $t\,=\,\{20\,,\,80\,,\,120\}\,\mathsf{h}$, enhancing the good accuracy of the solution to represent the physical phenomenon.

The absolute error of the reduced-order model is of the order of $\mathcal{O}\,\simeq(\,10^{\,-4}\,)$, as illustrated in Figure 2(a). The methods are compared with a reference solution, which has been computed using the Matlab open source toolbox Chebfun [30]. To give a solution with this order of accuracy, the PGD needed $22$ modes, while the Spectral only required $9$ modes. It corresponds to $21$ and $7$ degrees of freedom, respectively. Figure 2(b) presents the error $\varepsilon_{\,\infty}$ as a function of the number of modes. As we increase the number of modes, the solution of the Spectral–ROM converges faster than the PGD solution because the convergence of the Spectral method is exponential. To illustrate the convergence of the solution, the profile of the vapor pressure for the last time instant of simulation is represented with different numbers of modes in Figures 3(a) and 3(b) for the PGD and the Spectral–ROM, respectively.

The reduced system of ODEs, of size $7$ is implemented in Matlab, and the Spectral coefficients $\{a_{\,n}\,(\,t\,)\}$ are calculated for any intermediate time instants by the solver ODE45 [64]. These Spectral coefficients are shown in Figure 4(a). It can be seen that the first coefficients have the highest magnitude making the solution to converge with a few modes (an order of $10$), which happens thanks to the fact that the Chebyshev polynomials have excellent approximation properties for smooth function. In addition, the last coefficient determines the magnitude of the residual, implying that the error will not be lower than the magnitude of the last coefficient $a_{\,n}$, as explained in [32]. A brief comparison with an analytical solution, built on Fourier decomposition [55], reveals that the eigenvalues of the Spectral method decrease faster, as shown in Figure 4(b). Note that eigenvalues of the analytical solution do not have to coincide with the ones of the Spectral method since the basis functions are not the same for the Chebyshev polynomials and the trigonometric ones. Figure 5 shows the Fourier power spectrum function of the signal frequency per unit of time, generated by the fast Fourier transform. In this Figure, one peak is observed in the signal frequency which is consistent with the step of the relative humidity on the right boundary.

Regarding the PGD, Figures 6(a) and 6(b) present the first modes, depending on time and space, respectively. They do not have a physical meaning and constitute a separated representation of the solution. Their tensorial product enables to compute the solution of the problem. It is a similar approach for the Spectral–ROM, where the coefficients $a_{\,n}\,(\,t\,)$ are multiplied by the Chebyshev polynomials.

Remarks on a posteriori POD method

Here, the purpose is to compare the PGD and Spectral model order reduction to the well-established POD approach. For this, the reference solution $u$ obtained with the Matlab open source toolbox Chebfun was used to compute the correlation matrix R , with elements $\{r_{\,ik}\}$ given by:

[TABLE]

Then, the singular values $\lambda$ of the correlation matrix are computed. A truncation is operated in the eigenvectors basis to define $\phi_{\,i}$ the spatial basis function composed of the $Q$ eigenvectors of the correlation matrix $r\,$. Thus, the solution of the POD reduced-order model (POD–ROM) approximates the solution of the problem by:

[TABLE]

where $Q$ is the number of modes corresponding to the model order reduction. For this case study, Figure 7(a) shows the convergence of the solution obtained with a POD reduced-order model. Moreover, Figure 7(b) presents the error with the reference solution as a function of the modes $Q\,$. The solution of the POD–ROM converges faster than the PGD and Spectral approaches. Only $Q\ =\ 5$ modes are sufficient to compute a solution accurate to the order $\mathcal{O}\,\,(10^{\,-3})\,$, which is lower than the one for the PGD or Spectral–ROM for the same accuracy. However, as underlined in the Introduction section, the POD approach is a posteriori. To build the POD reduced-order model, a preliminary computation of the solution was required, which is a non-negligible restriction.

Computing parametric solution using reduced-order models

ding performance assessment. Consider, for instance, the analysis of the wall behavior in terms of heat and mass transfer, as a function of different parameters such as thermal inertia, vapor permeability, insulation thickness, among others. In the context of environmental issues and thermal regulations, the wall behavior may be optimized as a function of those parameters. Several studies of parametric simulations are presented in the literature. In [11, 54, 6, 37, 17, 5, 69], numerical methods are used to determine the optimum insulation thickness of different wall configurations. In [10], the influence of wall thermal inertia on the energy consumption was investigated by using the EnergyPlus program for $24$ construction types. In [41], the MBV of five hemp concrete materials were assessed using a numerical method. Those parametric simulations are based on models using numerical methods due to almost no restriction in terms of boundary conditions, geometry, material properties, among other considerations. Nevertheless, for parametric studies, they require large numbers of simulations. Indeed, the numerical model is not dependent on the parameters of interest. Thus, a computation of the numerical model is required for each value of the parameters within their domain of variation, demanding a high calculation cost, even after the dramatic evolution of computer hardware since the $1970$’s. Therefore, reduced-order model can be used to perform efficient parametric studies with a limited computational costs.

Extension of model reduction techniques to parametric problems

The issue of solving a parametric problem is to compute the solution $u$ of Eq. (3.1) depending on the usual time and space coordinates $x$ and $t$ as well as the parameter $\nu$ of the problem. Thus, the solution is seek as $u\,(\,x\,,t\,;\nu\,)\,$, where $\nu$ is defined as a coordinate of the problem within a given interval \nu\ \in\ \bigl{[}\,\nu_{\,\mathrm{min}}\,,\,\nu_{\,\mathrm{max}}\,\bigr{]}. Here we note $N_{\nu}$ the number of elements (cardinal) of the domain \bigl{[}\,\nu_{\,\mathrm{min}}\,,\,\nu_{\,\mathrm{max}}\,\bigr{]}\,.

5.1.1  Spectral reduced-order model

For the parametric study, we want to compute the solution as a function of $(\,x\,,t\,;\nu\,)\,$:

[TABLE]

The basis function will always depend only on $x$, that is why the parameter $\nu$ is calculated with Spectral coefficients $\{a_{\,n}\,(\,t\,;\nu\,)\}\,$. As it is not straightforward to compute the coefficients $a_{\,n}\,(\,t\,;\,\nu\,)$ depending on both parameters, we compute the solution for each value of the parameter $\nu\,$, using a loop. The latter can easily be parallelised on high-performance computer systems. It would be possible to vectorize the computation of the parametric Spectral–ROM solution although the method would significantly loose its speed calculation.

5.1.2  PGD reduced-order model

To compute a parametric solution of Eq. (3.1) $u\,(\,x\,,t\,;\nu\,)\,$, the PGD approach assumes a separated tensorial representation of the solution:

[TABLE]

Functions \bigl{(}\,F\,,\,G\,,\,H\,\bigr{)} are computed following the methodology described in Section 3.2. Interested readers may refer to [25, 23, 14] for complementary details on the methodology. The important point, is that the solution $u$ is computed at once as a function of the coordinates $x\,$, $t$ and $\nu\,$.

Case study

For this case, we seek for a parametric solution of problem Eq. (2.4). The vapor pressure is computed as a function of time $t$, space $x$ and moisture storage capacity $c_{\,m}$. As in the previous case, simulations are preformed in order to reproduce experiments that estimate the moisture buffer value of the materials. Thus, the right boundary condition is exposed to cyclic changes of relative humidity between $33\,\,\mathsf{\%}\,$ and $75\,\,\mathsf{\%}\,$, with a $24\;\,\mathsf{h}\,$ period. The total time of simulation is still $120\;\,\mathsf{h}\,$. The convective vapor coefficient is $h_{\,v}\ =\ 2\cdot 10^{\,-8}\;\,\mathsf{s/m}\,$. The left boundary is set to a constant vapor pressure, identical to the initial condition $P_{\,v}^{\,i}\ =\ 842\;\,\mathsf{Pa\,}\,$. Simulations undergo at a constant $23\,\mathsf{{}^{\circ}C}\,$ temperature. All $500\text{-}\,\mathsf{mm}\,$ materials have the same vapor permeability, $d_{\,m}\ =\ 2.4\cdot 10^{\,-11}\;\,\mathsf{s}\,$, while the moisture storage capacity varies in the interval $\Omega_{\,c}\ =\ 1.2\cdot 10^{\,-3}$ and $6\cdot 10^{\,-3}\;\,\mathsf{kg/m^{3}/Pa}\,$ [62].

First, we perform a simulation for $10$ different values of moisture storage capacity in the interval $\Omega_{\,c}$, representing different kinds of materials. Two different techniques of reduction order models were employed, the PGD and the Spectral–ROM. To validate these methods, results were compared to the reference solution, constructed with the Matlab open source toolbox Chebfun.

Figure 8(a) shows the mass content of each material among the simulation time. Even with a low difference between the highest and the lowest values of storage capacity, it is possible to observe significant variations of the weight as the the storage capacity increases, retaining more moisture. Furthermore, Figure 8(b) presents the weight for the two highest values of moisture storage capacity, corresponding to $6\cdot 10^{\,-3}$ and $4.15\cdot 10^{\,-3}$ $\mathsf{kg/m^{3}/Pa}$. The last profile of the pressure vapor for all values of moisture storage capacity is represented in Figure 9(a). In these figures, the PGD and the Spectral–ROM are in a good agreement with the reference solution.

To compute the parametric study, with the same order of accuracy, the PGD needed around $M\ =\ 100$ modes, while the Spectral–ROM used only $N\ =\ 11$ modes. This difference comes from the nature of the methods, that are constructed by different ways. The error $\varepsilon_{\,\infty}$ is shown as a function of the storage capacity values in Figure 9(b). The methods were constructed in order to give the same order of accuracy, around $\mathcal{O}\,\,(\,10^{\,-3}\,)\,$. It should be noted that the Spectral–ROM can give more accurate results, with the same degrees of freedom, by increasing the tolerance in the ODE45 Matlab function, when the reduced system is being computed. The degrees of freedom of the Spectral–ROM were predetermined based on the previous case and by the order of the parameters values. Meanwhile, the PDG computes its solution if the residual is lower than a given tolerance. The rate of convergence of the PGD approach is illustrated in Figure 10, presenting the error as a function of the number of modes.

The numbers of operation for each approach can be estimated, remembering $N_{\,x}$ and $N_{\,t}$ stand for the number of elements according to the discretization of the space and time domains, respectively. The quantity $N_{\,\nu}$ represents the number of elements of parameter $c_{\,m}$ considered for the parametric study. A standard approach based on implicit Euler schemes requires $N_{\,x}\cdot N_{\,t}\cdot N_{\,\nu}$ operations. Considering the discretization parameters to reach the given accuracy $N_{\,t}\ =\ 1.2\cdot 10^{\,5}$ and $N_{\,x}\,=\,4\cdot 10^{\,2}\,$, the number of operations scales with [33]:

[TABLE]

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

[TABLE]

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

[TABLE]

For this parametric case, knowing that the tolerance was set to $10^{\,-3}\,$, with $N\ \simeq\ 11$ modes, the number of operations performed by the Spectral–ROM is expressed as:

[TABLE]

For the PGD, the number of operations depends on the number of modes $M\,$, the number of iterations $S$ required for the fixed point algorithm to converge and the sum of the number of spatial, temporal and parameters elements $N_{\,x}\,$, $N_{\,t}$ and $N_{\,\nu}\,$, respectively [15]. For this application, the PGD fixed point algorithm requires around $S\ =\ 20$ iterations at each mode and $M\ =\ 100$ modes. The discretization parameter used for the PGD are $N_{\,t}\ =\ 1200$ and $N_{\,x}\ =\ 100\,$. Therefore, the number of operations for the PGD approaches scales with:

[TABLE]

Therefore, the number of operations increases faster for the Spectral–ROM and for the standard Euler methods than for the PGD approach. After some kind of initial investment, the increase in the number of operations for the PGD is much slower than for other approaches. The advantage of the PGD, in this case, it is related to its ability to compute at once the solution depending on the three coordinates, whereas the Spectral–ROM computes the solution for each value of moisture storage capacity independently, by a loop. It should be noted that the Euler approach, based for instance on backward time centered space, is a low order approximation of the solution, providing a less accurate solution than the Spectral–ROM. An interesting advantage of the PGD is the low storage cost of the solution thanks to the tensorial representation of the solution. This feature may be particularly interesting for real time applications. These features also impact the CPU time of each algorithm, which has been evaluated using Matlab platform on a computer with Intel i$7$ CPU and $8$GB of RAM. Figure 11 shows the CPU time as a function of the number of elements of the parameter $c_{\,m}$. For comparison, the CPU time required using the Euler implicit scheme is also reported. Since the Spectral–ROM has a reduced system to solve, its computational time drops significantly when compared to traditional methods. However, the loop to simulate the parametric study increases the CPU time linearly with the number of elements. For few numbers of elements of parameter $c_{\,m}$, around $20\,$, the Spectral–ROM is faster than the PGD. Yet, if the number of elements increases, the PGD is a more attractive method. In addition, the large original model, based on implicit Euler scheme, requires an important extra CPU time to compute the parametric solution.

Nonlinear transfer in porous material

The last case considers nonlinear transfer with moisture-dependent material properties. Therefore, the diffusion coefficient $\nu$ depends on the field $u$ and Eq. (3.1) becomes:

[TABLE]

Extension of model reduction techniques to nonlinear problems

Some adaptations of the methods must be carried out in order to consider nonlinear diffusion transfer. For this, in what follows the reduced-order models are described.

6.1.1  Spectral reduced-order model

In order to apply better the Spectral method, Eq. (6.1) is written in the non-conservative form:

[TABLE]

where,

[TABLE]

By using Spectral methods the unknown $u\,(\,x\,,\,t\,)$ is approximated by the finite sum (3.3) and, the derivatives can be written so that the Chebyshev polynomials remain the same, as in the linear case of Eq. (3.4). Thus, Eq. (6.2) becomes:

[TABLE]

The nonlinear terms \nu\,\Bigl{(}\,\sum_{i\,=\,0}^{n}\,a_{i}\,(\,t\,)\,T_{i}\,(\,x\,)\,\Bigr{)} and \lambda\,\Bigl{(}\,\sum_{i\,=\,0}^{n}\,a_{\,i}\,(\,t\,)\,T_{\,i}\,(\,x\,)\,\Bigr{)} are treated by applying the Tau–Galerkin method and the Chebyshev–Gauß quadrature [32]. Contrary to 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 adding two 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. (6.1.1). The initial condition is given by Eq. (3.8) and the DAE system is solved by ODE15s or ODE23t from Matlab. For further details, interested readers may consult [32, 64].

6.1.2  PGD reduced-order model

To treat the nonlinearity of the problem, at the enrichment step $m\ <\ M\,$, the nonlinear term $\nu\,(\,u\,)$ is approximated using the solution from previous steps:

[TABLE]

Then, the matrix of the coefficient $\nu$ is separated into a tensorial product in the space and time directions, using a Singular Value Decomposition (SVD) [35] or a Discrete Empirical Interpolation Method (DEIM) [20, 1]:

[TABLE]

This decomposition enables us to separate the coefficient into a component depending on the coordinate of the problem. Therefore, Eq. (3.12) becomes:

[TABLE]

Functions $R$ and $S$ are then computed using a similar approach as the one described in Section 3.2.

Case study

The material investigated is the wood fiber, which properties have been presented in [60]. The moisture transport coefficient $d_{\,m}$ is assumed as a first-degree polynomial of the relative humidity, while the moisture capacity $c_{\,m}$ as a second-degree polynomial:

[TABLE]

In terms of boundary conditions, a Robyn-type is assumed for both sides of the material, as described in Eq. (2.3). The variation of the relative humidity of the ambient air is given in Figure 12. Variations were chosen in order to excite the material in the hygroscopic region of the properties. The vapor convective transfer coefficients are set to $h_{\,v,\,\mathrm{L}}\ =\ \,1\cdot 10^{\,-8}\mathsf{s/m}$ and $h_{\,v,\,\mathrm{R}}\ =\ 1.5\cdot 10^{\,-8}\,\mathsf{s/m}$. As for the previous case, the time simulation is fixed to $5$ days.

Results have been computed using discretization parameters $\Delta x^{\,\star}\,=\,10^{\,-2}$ and $\Delta t^{\,\star}\,=\,10^{\,-2}$ for both methods. The Spectral–ROM has been built for $N\ =\ 8$ modes while the PGD for $M\ =\ 30$ modes. both results have been compared to a reference solution computed with the Matlab open source toolbox Chebfun. Profiles of relative humidity in the material are shown in Figure 13(a). The time evolution of relative humidity at $x\ =\ 0.074\,\mathsf{m}$ is given in Figure 13(b). A very good agreement is highlighted between the solutions. The physical phenomena are accurately represented. The relative humidity at $x\ =\ 0.074\,\mathsf{m}$ increases according to the variation of the boundary conditions. The reduced-order models have been built to give the same order of accuracy $\mathcal{O}\,\,(\,10^{\,-3}\,)$ as illustrated in Figure 14(a). It can be noted that the PGD needs $M\ =\ 30$ modes to compute the solution of the nonlinear problem, while only $M\ =\ 22$ modes were required in the linear case. The Spectral–ROM only needs one extra mode compared to the linear case. The error with the reference solution is given as a function of the number of modes for both reduced-order models in Figure 14(b). Again, the Spectral–ROM converges faster to an accurate solution than the PGD approach. A limit is observed in the error of the Spectral–ROM, around $\mathcal{O}\,\,(\,10^{\,-4}\,)\,$, due to the tolerance of the Matlab solver that was set to this value.

The CPU time of each method has been evaluated for the same order of accuracy of the solution and are reported in Table 1. For the comparison, the CPU time required with the classical Euler implicit scheme is also indicated. Both methods enable significant computational savings, $95\,\,\mathsf{\%}\,$ and $99.1\,\,\mathsf{\%}\,$ for the PGD and Spectral, respectively. The PGD requires more time than the Spectral–ROM, mainly due to the treatment of the nonlinearity of the problem. It can be noted that using the DEIM for the treatment of the nonlinearity permits to reduced by more than two the CPU time of the PGD, compared to the SVD. Indeed, in the latter case, at each iteration, the solution has to be composed to evaluate the nonlinear coefficients and then separate them along each coordinate of the problem. The CPU time of each approach is related to the number of operations. For the PGD, it scales with:

[TABLE]

where $N_{\,\mathrm{nl}}$ represents the number of operations for the treatment of nonlinearities. Depending on the method used for the decomposition of the solution, the number of operations scales with:

[TABLE]

where $\mathcal{O}\,\,(\,K\,)\,\simeq\,3$ is the order of the decomposition of the coefficients. It can be understood why the CPU time of the PGD using the DEIM is lower.

For the Spectral–ROM, thanks to the analytical pre-treatment of the solution, there is almost no increase in the number of operations. According to Eq. (6.2), the number of operations is only multiplied by two:

[TABLE]

For the large original model, using the Euler implicit scheme, the number of operations equals [33]:

[TABLE]

For this case, the average number of sub-iterations required required for the implicit scheme to treat the nonlinearity was around $\mathcal{O}\,\,(N_{\,\mathrm{nl}})\,\simeq\,(\,12\,)$ for a tolerance fixed to $0.01\cdot\Delta t^{\,\star}\,$.

Conclusion

Due to moisture-dependent material properties and weather driven boundary conditions, numerical methods are used to compute the solution of moisture transfer problems. Usual approaches based for instance on Euler or Crank–Nicolson schemes require the solution of large systems of equations, which imposes important numerical costs. Model reduction techniques appear then as efficient alternatives, enabling to reduce the model order without deteriorating the representation of the physical phenomena. These methods aim at preserving the computation resources in terms of CPU time and memory. Among the a priori model reduction techniques applied for diffusion problems, this paper intended to compare the Spectral and the PGD methods. It provided extensive details in order to provide a numerical benchmark for the community intending to building reduced-order models of diffusion problems in porous material. The case studies provided clear example of the construction and the use of these reduced-order models. Since the dimensionless properties of each case were provided in the Appendix A, these numerical benchmarks can be used for the whole community of transport in porous media.

The two reduced-order methods assume that the solution is approximated by a finite sum of functions products. The Spectral method fixes a set of basis function for the space domain. Here the Chebyshev polynomials have been chosen. Analytical preliminary treatment of the Spectral solution has been operated to set an ordinary system of equations to compute the temporal coefficients of the solution. On the other hand, the PGD has no assumptions and attempt to compute directly the basis functions by minimizing the equation residual. The comparison was carried out for three cases, commonly found in building physics. The first one deals with linear moisture transfer. The second one aimed at computing a parametric solution, whose model outputs depend not only on the space and time coordinates, but also on the moisture capacity of the material. The last case dealt with a nonlinear transfer problem, with moisture dependent material properties.

Results have demonstrated that both reduced-order models, Spectral and PGD, accurately represent the moisture transfer and both approaches provide an important reduction of the model order. While the order of the large original model scales with several hundred, the one of the ROMs is proportional to a few tens or even less. For the nonlinear case, thanks to this order reduction, the model reduction techniques enable to save more than $95\,\,\mathsf{\%}\,$ of the CPU time, compared to a large original model based on a Euler implicit scheme. If both methods are efficient, some distinctions between the two ROMs have been highlighted. For the linear and nonlinear cases, the Spectral–ROM has a lower order than the PGD, $\mathcal{O}\,\,(\,N\,)\,\simeq\,8$ against $\mathcal{O}\,\,(\,M\,)\,\simeq\,20$ to $30\,$, for the same accuracy. Moreover, the error of the Spectral–ROM decreases faster with the number of modes than the PGD. For these reasons, the Spectral–ROM is faster than the PGD. For the nonlinear case, the CPU time of the Spectral–ROM is divided by $5$ compared to the PGD. For the parametric case study, the two approaches compute the solution by different processes. The Spectral–ROM computes a solution for each numerical value of the material properties within a defined interval. Then, a loop is operated to compute the solution for each value of the material properties. The PGD approach considers directly the material properties as a coordinate of the problem within a defined interval of values. The solution is approximated by a tensorial representation and a basis of functions of each of the three coordinates is computed. Thus, the parametric solution is obtained at once. If the PGD ROM needs more modes than the Spectral method, the number of operations to compute the solution is smaller. Moreover, the increase of operations for the PGD is much slower, in this case, thanks to the tensorial representation of the solution. For a parametric solution depending on $50$ values of material properties, the CPU time of the PGD is six times faster.

To conclude the comparison of the two model reduction techniques, results have highlighted that the Spectral approach is more efficient in terms of order reduction, preserving computational resources for linear and nonlinear moisture diffusion problems. For the computation of parametric solutions, the PGD appears to be more efficient. These promising results encourage further investigation for two- or three-dimensional problems including combined heat and moisture transfer phenomena in porous media, where the order of the large original model becomes even higher.

Acknowledgements

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.

Nomenclature

[TABLE]

Appendix A Dimensionless values

Linear case

Problem (2.4) is taken into account with $g_{\,l,\,\mathrm{L}}^{\,\star}\ =\ g_{\,l,\,\mathrm{R}}^{\,\star}\ =\ 0\,$ and a Dirichlet condition on the left side:

[TABLE]

The dimensionless properties of the material are $d_{\,m}^{\,\star}\ =\ 1$ and $c_{\,m}^{\,\star}\ =\ 430\,$. The reference time is $t^{\,0}\ =\ 1$ $\mathsf{h}\,$, thus the final simulation time is fixed to $t^{\,\star}\ =\ 120\,$. The Biot number is $\mathrm{Bi}_{\,v,\,\mathrm{R}}\ =\ 333\,$. The boundary conditions are expressed as:

[TABLE]

Parametric case

Problem (2.4) is taken into account with $g_{\,l,\,\mathrm{L}}^{\,\star}\ =\ g_{\,l,\,\mathrm{R}}^{\,\star}\ =\ 0\,$ and a Dirichlet condition on the left side, the same as in the previous case. The reference time is $t^{\,0}\ =\ 1$ $\mathsf{h}\,$, thus the final simulation time is fixed to $t^{\,\star}\ =\ 120\,$. The Biot number is $\mathrm{Bi}_{\,v,\,\mathrm{R}}\ =\ 100\,$. The boundary conditions are expressed as:

[TABLE]

The dimensionless properties of the materials are $d_{\,m}^{\,\star}\ =\ 1$ and $c_{\,m}^{\,\star}$ assume the following values:

[TABLE]

Nonlinear case

Problem (2.4) is taken into account with $g_{\,l,\,\mathrm{L}}^{\,\star}\ =\ g_{\,l,\,\mathrm{R}}^{\,\star}\ =\ 0\,$ and Robin condition on both boundaries. The Biot number are $\mathrm{Bi}_{\,v,\,\mathrm{L}}\ =\ 10$ and $\mathrm{Bi}_{\,v,\,\mathrm{R}}\ =\ 15\,$. The boundary conditions are expressed as:

[TABLE]

The reference time is $t^{\,0}\ =\ 1$ $\mathsf{h}$, thus the final simulation time is fixed to $t^{\,\star}\ =\ 120\,$. The dimensionless properties of the materials are:

[TABLE]

Appendix B Solving the Spectral ODE reduced system

Consider a more general situation of the reduced-order system (3.9):

[TABLE]

where $\mathrm{b\,}\,(\,t\,,a\,(\,t\,)\,;\,\nu\,)$ depends on the solution $a\,(\,t\,)$ via nonlinear boundary conditions, or it contains problem’s nonlinearities, if there are some. The dependence on parameters is the most accurate within the chosen Spectral framework. The general analytical solution to problem (B.1) can be written as:

[TABLE]

The exponential matrix is defined as the limit:

[TABLE]

in which $\mathrm{Id}\in\mathrm{Mat}\,_{(n-2)\times(n-2)}(\mathds{R})\,$ is the identity matrix. However, this method is not the best way to compute the exponential matrix. In some particular cases, the solution of Eq. (B.1) can be simplified and thus better exploited [49].

Case I:

If we have homogeneous boundary conditions, problem (B.1) becomes:

[TABLE]

and it can be analytically solved as:

[TABLE]

Using modern methods, the exponential matrix can be computed using $\sim 48\,n^{\,3}$ floating point operations per second (FLOPS) [3]. As an information to the reader, the previous result was $\sim 538\,n^{\,3}$ FLOPS [40].

However, one can notice that we do not really need to build the exponential matrix, but we want to compute its action on the initial state vector $a_{\,0}$. Nowadays, it can be directly done, without forming $\mathrm{e}^{\,(\,t\ -\ t_{\,0}\,)\,\mathcal{A}\,(\,\nu\,)}\,$, explicitly to a prescribed accuracy that can be set significantly lower than the standard machine precision $\sim 10^{\,-16}$ [4]. If computing Eq. (B.1) by a Matlab solver, for example ODE45, the standard tolerance is of order of $\sim 10^{-6}\,$.

Case II:

If we have inhomogeneous boundary conditions constant in time, the problem from Eq. (B.1) becomes:

[TABLE]

which can also be analytically solved:

[TABLE]

Case III:

If we have inhomogeneous boundary conditions are linear in time, problem Eq. (B.1) becomes:

[TABLE]

the solution is given by:

[TABLE]

Case IV:

Boundary condition is polynomial in time:

[TABLE]

Then, the solution is given by the following analytical formula:

[TABLE]

Above we introduced the so-called matrix $\varphi-$functions:

[TABLE]

A few first functions are given below explicitly:

[TABLE]

The general power series representation of $\varphi-$functions is

[TABLE]

The exponential definitions of $\varphi_{\,m}\,(\,z\,)$ should not be used for practical simulations, because of severe cancellation errors for $z\,\ll\,1\,$. Efficient methods for computation of $\varphi-$functions have been developed based on Padé-type expansions, to give an example, Matlab’s function expm() is based on such approximations [27].

Case V:

For a general case of linear boundary conditions, the solution of problem Eq. (B.1) is:

[TABLE]

To exploit the last formula, one might employ a quadrature formula to discretize the integral (I):

[TABLE]

where we employed rectangle formula for simplicity. We note that the sequence $\{\mathrm{e}^{\,\Delta t\,\mathcal{A}\,(\,\nu\,)}\}_{\,j\,=\,1}^{\,m}$ can be entirely computed in an efficient manner [4].

Case VI:

Considering a general nonlinear case of boundary conditions from problem Eq. (B.1) and the general solution Eq. (B.2). To exploit a better solution, we can develop the function $\tau\,\mapsto\,\mathrm{b\,}\,(\,\tau\,,a(\tau)\,;\nu\,)$ in Taylor expansion series and integrate it exactly:

[TABLE]

where,

[TABLE]

Finally, the series solution can be exploited by truncating it at some finite order:

[TABLE]