An adaptive simulation of nonlinear heat and moisture transfer as a boundary value problem

TL;DR

This paper introduces an innovative numerical simulation approach for nonlinear heat and moisture transfer in porous materials, transforming the problem into a boundary value problem solved efficiently with adaptive collocation methods, validated through case studies.

Contribution

It proposes a novel time discretization method that converts the problem into a boundary value problem, enabling more efficient and accurate simulations of diffusion in porous materials.

Findings

Effective treatment of nonlinearities and interfaces

High-order adaptive methods improve accuracy

Good agreement with experimental data

Abstract

This work presents an alternative view on the numerical simulation of diffusion processes applied to the heat and moisture transfer through porous building materials. Traditionally, by using the finite-difference approach, the discretization follows the Method Of Lines (MOL), when the problem is first discretized in space to obtain a large system of coupled Ordinary Differential Equations (ODEs). Thus, this paper proposes to change this viewpoint. First, we discretize in time to obtain a small system of coupled ODEs, which means instead of having a Cauchy (Initial Value) Problem (IVP), we have a Boundary Value Problem (BVP). Fortunately, BVPs can be solved efficiently today using adaptive collocation methods of high order. To demonstrate the benefits of this new approach, three case studies are presented, in which one of them is compared with experimental data. The first one considers nonlinear heat and moisture transfer through one material layer while the second one considers two material layers. Results show how the nonlinearities and the interface between materials are easily treated, by reasonably using a fourth-order adaptive method. Finally, the last case study compares numerical results with experimental measurements, showing a good agreement.

Tables

Table 1. Table 1. Maximum error ε ∞ subscript 𝜀 \varepsilon_{\,\infty} computed with different tolerances with the solver bvp4c .

$\mathrm{𝗍𝗈𝗅}$	${\epsilon }_{\infty ,u}$	mesh	${\epsilon }_{\infty ,v}$	mesh	CPU time
${10}^{-6}$	$3.5\cdot {10}^{-5}$	$19$	$2.1\cdot {10}^{-4}$	$21$	$21.9𝗌$
${10}^{-5}$	$1.5\cdot {10}^{-4}$	$10$	$4.5\cdot {10}^{-4}$	$15$	$15.9𝗌$
${10}^{-4}$	$1.5\cdot {10}^{-4}$	$10$	$8.9\cdot {10}^{-4}$	$11$	$15.3𝗌$
${10}^{-3}$	$1.5\cdot {10}^{-4}$	$10$	$5.4\cdot {10}^{-4}$	$10$	$13.9𝗌$
${10}^{-2}$	$1.5\cdot {10}^{-4}$	$10$	$2.3\cdot {10}^{-4}$	$10$	$13.3𝗌$
${10}^{-1}$	$1.5\cdot {10}^{-4}$	$10$	$7.2\cdot {10}^{-3}$	$10$	$13.1𝗌$

Table 2. Table 2. Features and CPU time of the methods implemented in the single layer case study.

	BVP	Euler implicit	Euler explicit
$\Delta t^{\star }$	${10}^{-1}$	${10}^{-2}$	${10}^{-5}$
$N_{x,u}$	$10$	$101$	$201$
$N_{x,v}$	$15$	$101$	$201$
${\epsilon }_{\infty ,u}$	$𝒪\left({10}^{-4}\right)$	$𝒪\left({10}^{-4}\right)$	$𝒪\left({10}^{-4}\right)$
${\epsilon }_{\infty ,v}$	$𝒪\left({10}^{-4}\right)$	$𝒪\left({10}^{-4}\right)$	$𝒪\left({10}^{-4}\right)$
CPU time $\left[𝗌\right]$	$18.5$	$57.5$	$2572$
CPU time $[-]$	$0.49$	$1$	$44$

Table 3. Table 3. Hygrothermal properties of the load bearing material [ 19 ] .

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

Table 4. Table 4. Hygrothermal properties of the finishing material [ 19 ] .

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

Table 5. Table 5. Hygrothermal properties of the wood fibre [ 40 ] .

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

Equations

$$\frac{\partial {\rho }_w}{\partial t}$$

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

$$P_c$$

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

$$\frac{\partial P_c}{\partial x}=\frac{\partial P_c}{\partial P_v}\frac{\partial P_v}{\partial x}+\frac{\partial P_c}{\partial T}\frac{\partial T}{\partial x}\simeq \frac{{\rho }_l{R}_{v}T}{P_v}\frac{\partial P_v}{\partial x}.$$

$$\frac{\partial {\rho }_w}{\partial t}=\frac{\partial {\rho }_w}{\partial \varphi }\frac{\partial \varphi }{\partial P_v}\frac{\partial P_v}{\partial t}+\frac{\partial {\rho }_w}{\partial T}\frac{\partial T}{\partial t}\simeq \frac{\partial {\rho }_w}{\partial \varphi }\frac{\partial \varphi }{\partial P_v}\frac{\partial P_v}{\partial t}.$$

$$\frac{\partial {\rho }_w}{\partial t}=\frac{w^{\prime }(\varphi )}{P_s(T)}\frac{\partial P_v}{\partial t}.$$

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

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

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

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

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

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

$$\displaystyle c_{\,TM}\$$

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

$$c_T(T,P_v)\frac{\partial T}{\partial t}+c_{TM}(T,P_v)\frac{\partial P_v}{\partial t}$$

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

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

$$\displaystyle\ +\ L_{\,v}\,h_{\,M}\,\Bigl{(}\,P_{\,v}-P_{\,v,\,\infty}\,(t)\,\Bigr{)}\,,$$

$$q_{\infty }=g_{\infty }\cdot H_l,$$

$$\displaystyle q_{\,\infty}\ =\ \sum_{j=1}^{m}\,\mathrm{s}\,\xi\,\sigma\,\biggl{[}\,\Bigl{(}\,T_{\,j}\,(x\,=\,0)\,\Bigr{)}^{\,4}\ -\ \Bigl{(}\,T\,(x\,=\,0)\,\Bigr{)}^{\,4}\,\biggr{]}\,,$$

$$P_v(x,t=0)=P_{v,i}$$

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

$$g(t)\stackrel{\mathrm{def}}{:=}-k_M(T,P_v){\frac{\partial P_v}{\partial x}|}_{x_0},$$

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

$$c_T^{\star }(v)\frac{\partial u}{\partial t^{\star }}+{\gamma }_1{c}_{TM}^{\star }(v)\frac{\partial v}{\partial t^{\star }}$$

$$\displaystyle\left[\begin{array}[]{cc}\dfrac{c_{\,M}^{\,\star}}{\mathrm{Fo}_{\,M}}&0\\[10.0pt] \dfrac{\gamma_{\,1}\,c_{\,TM}^{\,\star}}{\mathrm{Fo}_{\,T}}&\dfrac{c_{\,T}^{\,\star}}{\mathrm{Fo}_{\,T}}\end{array}\right]\cdot\left[\begin{array}[]{c}\dfrac{\partial v}{\partial t}\\[10.0pt] \dfrac{\partial u}{\partial t}\end{array}\right]\ =\ \dfrac{\partial}{\partial x}\left(\left[\begin{array}[]{cc}k_{\,M}^{\,\star}&0\\[10.0pt] \gamma_{\,2}\,k_{\,TM}^{\,\star}&k_{\,T}^{\,\star}\end{array}\right]\cdot\left[\begin{array}[]{c}\dfrac{\partial v}{\partial x}\\[10.0pt] \dfrac{\partial u}{\partial x}\end{array}\right]\right)\,.$$

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

$$\displaystyle\mathbf{n}\cdot\biggl{(}\,k_{\,T}^{\,\star}\,(v)\,\dfrac{\partial u}{\partial x^{\,\star}}\ +\ \gamma_{\,2}\,k_{\,TM}^{\,\star}\,(v)\,\dfrac{\partial v}{\partial x^{\,\star}}\,\biggr{)}\ =\$$

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

$$\displaystyle u\$$

$$\displaystyle u_{\,\infty}\$$

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

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

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

$$\displaystyle\mathrm{Fo}_{\,M}\$$

$$\displaystyle\gamma_{\,1}\$$

$$\frac{\partial u}{\partial t}=\nu \frac{{\partial }^{2}u}{\partial x^2},x\in \mathcal{I}\equiv [0,\ell ],t⩾0,$$

$$\mathcal{I}_{\,\Delta x}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \bigl{\{}x_{\,j}\ \equiv\ j\cdot\Delta x\ |\ j\ =\ 0,\,1,\,2,\,\ldots\,N,\,N+1\bigr{\}}\,,\qquad\Delta x\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \frac{\ell}{N+1}\,.$$

$$\frac{\mathrm{d}u_{\,j}}{\mathrm{d}t}\ =\ \frac{\nu}{\Delta x^{\,2}}\,\bigl{(}\,u_{\,j-1}\ -\ 2\,u_{\,j}\ +\ u_{\,j+1}\,\bigr{)}\,,\qquad j\ =\ 1,\,2,\,\ldots,\,N\,,$$

Full text

Key words and phrases: coupled heat and moisture diffusive transfer; boundary value problems; semi-discretization in time; semi-discretization in space; bvp4c

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Physical model

3. 2.1 Dimensionless representation

4. 3 The Method of Horizontal Lines

5. 3.1 Scalar nonlinear equations

6. 3.2 Problem reformulation as a BVP

7. 3.3 Numerical methods

8. 3.4 Extension to the coupled transfer

9. 4 Numerical benchmark

10. 4.1 Single layer case

11. 4.2 Multilayered domain

12. 4.2.1 Application

13. 5 Comparison of the numerical predictions with experimental data

14. 6 Conclusions

15. A Material properties

16. B Dimensionless parameters

17. B.1 Case from Section 4.1

18. B.2 Case from Section 4.2

19. B.3 Case from Section 5

20. C Analytical solution for the interface assumptions

Introduction

The hygrothermal transfer through porous structures is a matter of concern in many areas such as building physics, geophysics, environmental engineering and energy systems whose the transient evolution of heat and moisture migration plays an important role. Particularly, in the area of building physics, the heat and moisture transfer process through the porous envelope, roofing systems and the ground can strongly affect the energy efficiency, the thermal comfort of the occupants and the durability of the components [47, 44, 4]. Therefore, reliable assessment of hygrothermal transfer in building materials is a major issue, requiring efficient numerical tools for heat and moisture transfer in building materials [7].

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 discretization and incremental techniques are applied, such as the Euler implicit scheme in [29, 22] to solve large systems of equations. 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 [15]. The difficulties to compute the solution increase, particularly when using implicit schemes that require sub-iterations to treat those issues. In the literature [12, 9, 1], the important numerical costs of simulation tools are also mentioned and it is a matter of concern due to the substantial scale of buildings, where heat and moisture transfer phenomena have to be simulated.

In addition, in the models proposed in literature, the problem previously described is generally solved by traditional approaches such as the finite-difference method [17], the finite-volume method [13, 29, 28] and the finite-element method [45, 23, 39], which are well established in the fields of thermal sciences and building physics. In these classical approaches, the higher accuracy obtained for the space discretization of the numerical schemes is to order $\mathcal{O}\,\,(\Delta x^{\,2})\,$. For a space standard discretization $\Delta x\,=\,10^{\,-2}\,$, it implies that the error $\varepsilon$ on the solution of the equations cannot be lower than $\mathcal{O}\,(10^{\,-4})\,$. Within the issue of comparing the model numerical predictions with experimental observations as carried out for instance in [26, 34, 32], it is of major importance to control the accuracy of the solution.

For sure, the accuracy of the computed solution can be increased by reducing the space and time discretization parameters. However, as mentioned before, the standard approaches proposed in literature has a high degree of freedom. Therefore, increasing the number of spatial and temporal grid points will inevitably increase the computational time of the numerical model. With high-order numerical schemes, it is possible to have the same precision of low-order numerical schemes but with a lower computational cost, as shown in [30]. For this reason, these traditional methods have to be improved or even replaced by innovative and efficient ways of numerical simulation [33], particularly with the issue of comparing the model predictions with experimental observations.

Therefore, this article aims at contributing to the numerical development of hygrothermal transfer, by proposing a new approach to simulate the one-dimensional heat and moisture diffusive transfer trough single and multilayered building porous materials. The Method of Horizontal Lines is here proposed to solve the nonlinear heat and moisture transfer to increase significantly the accuracy in space with a low computational time of the numerical model. Usually, when using finite-differences, the discretization follows the Method of Lines (MOL). It means that the problem is first discretized in space to obtain a large system of coupled ODEs. Here, a different point of view is proposed based on discretizing first in the time domain to obtain a Boundary Value Problem (BVP). Such problems can be easily solved using adaptive collocation methods of high order. This approach is investigated in this paper to compute with high accuracy combined heat and mass transfer problems in porous materials.

The manuscript is organized as follows. Section 2 details the physical model of heat and moisture transfer while fundamentals of the proposed method are shown in Section 3. Numerical results are discussed in Section 4 and simulation are compared with experimental data in Section 5. Finally, in Section 6, the main conclusions are outlined with future perspectives.

Physical model

The physical problem considers one-dimensional heat and moisture transfer through a porous material defined by the spatial domain $\Omega_{\,x}\,=\,[\,0,\,L\,]$ and time domain $\Omega_{\,t}\,=\,[\,0,\,\tau\,]\,$. The following convention is adopted: $x\,=\,0$ corresponds to the surface in contact with the inside room and, $x\ =\ L\,$, corresponds to the outside surface. The moisture transfer occurs due to capillary migration and vapour diffusion. The heat transfer is governed by diffusion and latent mechanisms. The physical problem can be formulated as [27, 36]:

[TABLE]

where $\rho_{\,w}$ is the volumetric moisture content of the material, $\delta_{\,v}$ and $k_{\,l}\,$, the vapour and liquid permeabilities, $P_{\,v}\,$, the vapour pressure, $T\,$, the temperature, $R_{\,v}\,$, the water vapour gas constant, $P_{\,c}\,$, the capillary pressure, $c_{\,0}\,$, the material heat capacity, $\rho_{\,0}\,$, the material density, $c_{\,w}\,$, the water heat capacity, $\lambda\,$, the thermal conductivity, and, $L_{\,v}\,$, the latent heat of evaporation. Equation (2.1a) 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}$, and the Clausius–Clapeyron equation:

[TABLE]

where $P_{\,v}$ comes from the relation $\phi\,=\,\dfrac{P_{\,v}}{P_{\,s}\,(\,T\,)}\,$, in which $\phi$ is the relative humidity. Thus, by neglecting the variation of the capillary pressure and the mass content with temperature [39], the partial derivative of $P_{\,c}$ can be written as:

[TABLE]

In addition, we have:

[TABLE]

Considering the relation $\rho_{\,w}\,=\,w\,(\phi)\,$, obtained from the sorption isotherm, and from the relation between the vapour pressure $P_{\,v}$ and the relative humidity $\phi\,$, we get:

[TABLE]

We denote by

[TABLE]

Considering the previous notation, Equation (2.1) can be rewritten as:

[TABLE]

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

[TABLE]

where $P_{\,v,\,\infty}$ and $T_{\,\infty}$ stand for the vapour pressure and temperature of the air, $h_{\,M}$ and $h_{\,T}$ are the convective transfer coefficients and $\mathbf{n}$ is the normal that assumes $+\,1$ or $-\,1$ at the left or right boundary sides. If the bounding surface is in contact with the outside air, $g_{\,\infty}$ is the liquid flow from wind driven rain and $q_{\,\infty}$ is the sensible heat from the rain:

[TABLE]

where, $H_{\,l}$ is the water enthalpy. If the bounding surface is in contact with the inside building air then $g_{\,\infty}\,=\,0$ and the variable $q_{\,\infty}$ is the enclosure and long-wave radiative heat exchanged among the room surfaces:

[TABLE]

where $\mathrm{s}$ is the view factor between two surfaces, $\sigma$ is the Stefan–Boltzmann constant, $\xi$ is the emissivity of the wall surface, $j$ represents the $m$ bounding walls.

We consider a uniform vapour pressure and temperature distributions as initial conditions:

[TABLE]

One of the interesting outputs in the building physics framework are the sensible $q_{\,s}$ and the latent $q_{\,l}$ heat fluxes [46], which are respectively defined as:

[TABLE]

The moisture flow $g$ is similarly computed:

[TABLE]

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

Dimensionless representation

Before solving this problem directly, it is of great importance to get a dimensionless formulation. It enables us to determine important scaling parameters such as Biot and Fourier numbers. It allows also to estimate the relative magnitude of various terms in governing equations, and thus, eventually to simplify the problem using asymptotic methods [31]. In addition, the dimensionless form enables to manipulate numerically the quantities at the order of $\mathcal{O}\,\,(1)$ where the floating point arithmetics are designed to have minimal rounding errors [24]. Considering the temperature range of interest in building applications, temperature dependence was neglected when compared to their dependence with moisture content, with all transport coefficients considered as moisture content dependent. Therefore, the governing Equation (2.2) can be written in a dimensionless form as:

[TABLE]

which can be expressed by the following equation in a matrix form:

[TABLE]

The dimensionless formulation of the boundary conditions is:

[TABLE]

where the dimensionless quantities are defined as:

[TABLE]

where the subscript [math] represents a reference value, chosen according to the application problem and the superscript $\star$ represents a dimensionless quantity of the same variable. The dimensionless values of each numerical application can be found in the Appendix B.

The Method of Horizontal Lines

The perception of objects might change depending on the viewpoint. Keeping that in mind, in this study we apply the same idea to Partial Differential Equations (PDEs). Since our main applications include the heat and/or moisture transfer, we explain the idea on the simplest parabolic (or Fourier) equation, which is traditionally written as [16, 14]:

[TABLE]

where $\nu\,>\,0$ is the diffusion coefficient. The last equation has to be supplemented with right and compatible initial and boundary conditions to obtain a well-posed problem [18].

The way how Equation (3.1) is written conditions the tools and, thus, the way how it is solved numerically (or study it theoretically [14]). Indeed, the most popular approach nowadays is the so-called Method Of Lines (MOL) [42, 35, 25, 41, 21]. It consists in semi-discretization in space as the first step. With this goal in mind, let us consider a uniform (for the sake of simplicity) discretization in space:

[TABLE]

The classical second-order central finite difference scheme yields the following system of coupled ordinary differential equations (ODEs) [37]:

[TABLE]

with $u_{\,j}\,(t)\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,u\,(x_{\,j},\,t)$ is the solution nodal value. For simplicity of the discussion, we write here evolution equations only for interior nodes. There are still two remaining equations dependent on the boundary conditions which will be discussed in a further Section.

Hence, the Partial Differential Equation (3.1) solved with MOL yields approximatively to a large system of coupled differential equations. In typical accurate numerical simulations $N\,\gtrsim\,10^{\,2}\,\ldots\,10^{\,3}\,$. Then, the last system of equations can be discretized using various explicit or implicit time-marching schemes [6]. This philosophy is shown schematically in Figure 1(a). To make the language more accurate, we can refer to this approach as the Method of Vertical Lines (MOVL). From a brief literature survey we can conclude that the dominant majority of numerical simulations in the field of HAM transfer follow this philosophy. However, in this study we would like to propose an alternative view on diffusion equations. Namely, we claim that Equation (3.1) can be solved by a totally different numerical technique, if it is rewritten in Cauchy–Kovalevskaya form as:

[TABLE]

By defining a uniform time discretization:

[TABLE]

Equation (3.2) yields to the following system:

[TABLE]

The main term now is the second-order elliptic differential operator in space, while the time derivative can be seen (and treated) as a source term of secondary importance. In other words, instead of having an evolution problem, we have a Boundary Value Problem (BVP), which can be tackled by appropriate methods. This approach is schematically depicted in Figure 1(b), which will be called in this work as the Method of Horizontal Lines (MOHL). In the following Section 3.1, we propose a reformulation for scalar nonlinear equations.

Scalar nonlinear equations

Consider the following dimensionless nonlinear diffusion equation, which can describe the heat or the moisture transfer in a 1D section111In higher dimensional problems we would require the domain $\mathcal{I}\ \subset\ \mathds{R}^{\,m}\,$ \bigl{(}m\ \in\ \bigl{\{}2,\,3\bigr{\}}\bigr{)} to be bounded with a Lipschitz-continuous boundary. of a porous material:

[TABLE]

where the right-hand side \mathscr{F}\,(x,\,t)\,\in\,L^{\,2}\bigl{(}(0,\,T),\,L^{\,2}\,(\mathcal{I})\bigr{)} and describes eventual sources/sinks inside the material. The time horizon $T\ >\ 0$ is fixed. We assume that coefficients $c\,(u)$ and $\mathscr{D}\,(u)$ are positive definite:

[TABLE]

For the moisture diffusion, the coefficient $c\,(u)$ is the moisture storage coefficient $c_{\,M}^{\,\star}\,(u)$ and the coefficient $\mathscr{D}\,(u)$ is the diffusion coefficient $\mathrm{Fo}_{\,M}\cdot k_{\,M}^{\,\star}\,(u)\,$. For the heat transfer, $\mathscr{F}\,(x,\,t)$ incorporates the other terms. The subscripts denote partial derivatives, i.e., $u_{\,t}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{\partial\/u}{\partial\/t}\,$, $u_{\,x}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{\partial\/u}{\partial\/x}\,$, *etc.*Equation (3.3) has to be completed by initial and boundary conditions. The boundary conditions are described by Equation (2.13), which are of (mixed) Robin-type and the initial condition as $u_{\,0}\,(x)\,\equiv\,u\,(x,\,0)\,\in\,H^{\,1}(\mathcal{I})\,$. Then, the Boundary-Value Problem (BVP) (3.3) has a unique solution u\,(x,\,t)\,\in\,L^{\,2}\bigl{(}(0,\,T),\,H^{\,1}\,(\mathcal{I})\bigr{)}\,.

Problem reformulation as a BVP

For our purposes it will be more advantageous to recast Equation (3.3) in a non-conservative form. By assuming that function $\mathscr{D}\,(u)$ is continuously differentiable, we introduce its derivative:

[TABLE]

After taking one derivative, Equation (3.3) becomes:

[TABLE]

where we introduce the effective coefficients:

[TABLE]

and also the effective source term:

[TABLE]

Equation (3.5) is equivalent to Equation (3.3) for smooth solutions. Under assumptions (3.4) the solutions are indeed smooth.

Equation (3.5) is an evolution equation written in time–space: first we write the derivatives with respect to time, then the spatial ones. According to the philosophy of the approach described in Section 3, let us rewrite it in space–time:

[TABLE]

In other words, we see the time derivative $u_{\,t}$ as a source term for our convenience. Finally, we rewrite the last equation as a system of first-order differential equations in the spatial variable $x\,$:

[TABLE]

The last system can be rewritten with original coefficients as well:

[TABLE]

In order to have a true BVP, we have to get rid of the solution dependence on time. To do it, we (semi-)discretize the solution $u$ in time, i.e., we replace the surface $u\,(x,\,t)$ by a sequence of lines:

[TABLE]

where $\Delta t\,>\,0$ is a chosen time step. See Figure 1(b) for an illustration. This method is called the Method Of Horizontal Lines (MOHL). At each time layer $t^{\,n}$ we have to solve a true BVP:

[TABLE]

where \bigl{(}u_{\,t}^{\,\star}\bigr{)}^{\,n} is an approximation of the time derivative. Depending on the accuracy we would like to achieve, we can take the following backward finite-difference formulas:

[TABLE]

In this case, the second-order accuracy in time is implemented in the calculations in order to have a more precise solution, so that we can use $\Delta t^{\,\star}\,\geqslant\,10^{\,-1}\,$. The final scheme is unconditionally stable according to the construction of the method.

Numerical methods

All numerical results in this paper were computed using Matlab™2017b.The solver bvp4c [43] is employed to solve the BVP problem at every time step, other solvers also can be employed such as the bvp5c [43] and even bvp6c [20]. The codes for bvp4c & bvp5c are available within any standard Matlab™ distribution, while the code bvp6c was developed by Dr. Nick Hale (Stellenbosh University, South Africa) and is freely available. All these methods are based on finite-difference approximations that implement various orders of Lobatto IIIA formula. This is a collocation method and the corresponding collocation polynomial provides a $C^{\,1}\,(\mathcal{I})$ approximation of the uniform fourth, fifth or sixth orders of accuracy in $\Delta x$, respectively. Special attention should be given to the fact that other implementation details are quite different among the solvers, the order is not the only difference between them. In our numerical simulations we use adaptive methods of the fourth order which in most practical applications are enough.

Let us consider a system of ordinary differential equations of the form $u^{\,\prime}\,(x)\,=\,f\,(x,\,u)\,$, within the interval $[\,a,\,b\,]$ subjected to two boundary conditions $\psi\,(u\,(a),\,u\,(b))\,=\,0$. To use any of the BVP solvers, three inputs are required: the initial guess, a function with the boundary conditions and another function with the system of ordinary equations. They will return essentially three outputs: the spatial grid mesh, the solution approximation of $u\,(x)$ at the mesh points and the approximation to $u^{\,\prime}\,(x)$. Note that the solvers produce a solution that is continuous in the considered interval $[\,a,\,b\,]$ and with a continuous first derivative. For more details on the methods and their implementations refer to [43, 20].

The essential feature of these algorithms is the adaptive distribution of collocation nodes. The grid adaptation and error control are based on the residuals (there are two: one for the equation and one for the boundary conditions) of the $C^{\,1}\,(\mathcal{I})$ continuous solution. The convergence speed towards the BVP solution depends essentially on the quality of the initial guess. Here we deal with an IVP-BVP problem. Thus, the BVP-solution value on the previous time step can be a good approximation. However, some physics-based or vector-based extrapolation techniques might also be used to further reduce internal iterations on every time step. The boundary conditions are directly provide as one of the inputs of the solver. No special treatment need to be carried out before supplying it.

Extension to the coupled transfer

There are two ways of using the BVP approach in the heat and mass transfer System (2.3). The first one is to take advantage of the solver bvp4c that provides directly the derivative of the field as output. First, it is computed the solution for $v$ and then, it is computed the solution for $u\,$, which is the approach used in the application Section. This approach works because System (2.3) is weakly coupled. The procedure to implement the MOHL is described in Algorithm 1.

Therefore, using this approach, the system of first-order differential equations for the mass balance equation Equation (2.3a) is written as:

[TABLE]

while the system of the energy equation Equation (2.3b) as:

[TABLE]

where,

[TABLE]

In addition, the boundary conditions that complement Equation (3.8) are:

[TABLE]

and, the boundary conditions that complement Equation (3.11) are:

[TABLE]

where, $\mathrm{R}$ represents the residual and the subscripts $L$ and $R$ represent values of the left and right boundary conditions.

When dealing with highly coupled equations, the issue is to rewrite System (2.3) in a vectorized way, such as Equation (2.12), and solve it at as only one equation. In the next section, the MOHL is tested in two cases of numerical application.

Numerical benchmark

Simulations of one-dimensional coupled heat and moisture transport are carried out with the MOHL method, the Euler implicit method used in a numbers of works (see for instance [29, 22]) and a reference solution to verify the applicability and the advantages of the proposed method. The Euler implicit scheme approximates the continuous operator to the order $\mathcal{O}\,\,(\Delta x^{\,2}\,+\,\Delta t)\,$. The reference solutions are computed using the Matlab™ open source toolbox Chebfun [22].

To compare and validate the proposed method, the error between solution components $u^{\,\mathrm{num}}$ and $v^{\,\mathrm{num}}$, obtained by the MOHL method or the Euler implicit, and the reference solutions $u^{\,\mathrm{ref}}$ and $v^{\,\mathrm{ref}}$, are computed as functions of $x$ using the following Euclidean norms:

[TABLE]

where $N_{\,t}$ is the number of temporal steps. The uniform norm errors $\varepsilon_{\,\infty,\,u}$ and $\varepsilon_{\,\infty,\,v}$ are given by the maximal values of $\varepsilon_{\,2,\,u}\,(x)$ and $\varepsilon_{\,2,\,v}\,(x)\,$:

[TABLE]

In the aftermath, the numerical case studies are presented.

Single layer case

In order to contemplate the applicability of the presented method, the first case study considers the combined effects of the heat and moisture transfer. A physical application is performed in a $10$-$\mathsf{cm}$ monolithic bearing material. To test the robustness of the scheme with strong nonlinearities, the properties of the material are gathered from [19] and reminded in Table 3 of the Appendix A. Initial conditions are considered uniform over the spatial domain, with an initial temperature of $T_{\,i}\,=\,293.15\ \mathsf{K}$ and an initial vapour pressure of $P_{\,v,\,i}\,=\,1636.53\ \mathsf{Pa}\,$, referent to the relative humidity of $70\,\%\,$.

The boundary conditions, represented by the relative humidity $\phi$ and temperature $T\,$, are shown in Figure 2. They oscillate during $72$ hours of simulation. The convective mass and heat transfer coefficients are set to $h_{\,M}\,=\,2\cdot 10^{\,-7}\ {\mathsf{s}}/{\mathsf{m}}$ and $h_{\,T}\,=\,25\ {\mathsf{W}}/{\mathsf{(m^{\,2}\cdot K)}}$ at the left boundary, and to $h_{\,M}\,=\,3\cdot 10^{\,-8}\ {\mathsf{s}}/{\mathsf{m}}$ and $h_{\,T}\,=\,8\ {\mathsf{W}}/{\mathsf{(m^{\,2}\cdot K)}}\,$, at the right boundary. No additional source term is considered for the moment.

Simulations with the MOHL method were performed using the solver bvp4c with the relative and absolute tolerances set to $\mathsf{tol}\,=\,10^{\,-5}$ for moisture and $\mathsf{tol}\,=\,10^{\,-6}$ for temperature. The time step discretization is of $\Delta t^{\,\star}\,=\,10^{\,-1}$ and, for the space discretization an adaptive technique is employed, as explained in Section 3.3. The initial guess is composed by $10$ spatial nodes, for vapour pressure and temperature fields. For the Euler Implicit solution, simulations have been performed with $\Delta x^{\,\star}\,=\,10^{\,-2}$ and $\Delta t^{\,\star}\,=\,10^{\,-2}$.

Figure 3 presents the selected mesh determined by the solver for the last profile of temperature and vapour pressure fields, with the reference solution and with the refined solution. This last one is an evaluation of the solution that covers all the interval $[\,0,\,L\,]$ in more points than the selected mesh. Note that at each time layer the number of points of the selected mesh varies, with a maximum of $15$ points for the moisture field and $19$ points for the temperature field. Over the simulation time the number of nodes remained almost the same, which means it was sufficient for the set tolerance. If we decrease the tolerance, the number of nodes of the spatial grid would rise accordingly.

Simulation results of the temperature, vapour pressure and relative humidity evolution on the boundaries of the building component ($x\,=\,0\ \mathsf{m}$ and $x\,=\,0.1\ \mathsf{m}\,$) are presented in Figures 4(a), 4(b) and 4(c), respectively. The proposed method has been able to follow successfully the solicitations from the boundary conditions and represent the solution. The step on the relative humidity at the right boundary can be observed on the vapour pressure and relative humidity evolution. The vapour pressure diffuses over the material almost instantaneously with the step imposed. The level of moisture is high but the material is not saturated. It can be noted that for the simulation of such cases with really high levels of moisture, it is better to change the physical model using, for instance, the moisture content or the capillary pressure as potentials [29].

The computed temperature and relative humidity distribution profiles for $t\,=\,\{\,30,\,36,\,40\,\}\ \mathsf{h}$ are presented in Figures 5(a) and 5(b). As this material has a fast liquid transfer, the profiles of relative humidity are almost straight lines. The material is highly permeable but with a low hygroscopicity, as it does not retain the moisture. For the temperature, it can also be observed that it diffuses quickly and it has the opposite effect of the relative humidity. Furthermore, it can be noted a good agreement between the MOHL method, the Euler implicit and the reference solution.

The error $\varepsilon_{\,2}$ between the MOHL and the Euler implicit method with the reference solution are given in Figure 6. For the temperature distribution, the error is to the order of \mathcal{O}\,\,\bigl{(}10^{\,-5}\bigr{)} and, for the vapour pressure component is to the order of \mathcal{O}\,\,\bigl{(}10^{\,-3}\bigr{)}\,. The error is not of the same order for each solutions perhaps due to the following reasons: (i) the properties are more nonlinear for the moisture balance equations; (ii) or the residual effect of initial guess for the vapour pressure profile; (iii) or moisture diffusive is much slower than heat, for the present case; (iv) the order of the dimensionless fields $u$ and $v$ are not the same. Even with this difference, these results are consistent to the tolerance, that was set to $\mathsf{tol}\,=\,10^{\,-5}$ and to the discretization on time, that is second-order accurate $\mathcal{O}\,\,(\Delta t^{\,2})\,$. It is admissible that the actual error is slight larger than the prescribed tolerance. The solver estimates the residual and not the actual error strictu sensu (since the exact solution is unavailable). It is required that the error decrease when the tolerance decrease as shown in Table 1.

The accuracy of the solution depends on several factors, one of them is the absolute and relative tolerances of the solver. To compare different values of the tolerance with the solver bvp4c, Table 1 presents the error $\varepsilon_{\,\infty}$ for both solutions, the dimensionless temperature ($u$) and the dimensionless vapour pressure ($v$), as well as the number of mesh points used by the solver. As the tolerances get more restrictive, the solver adapts the solution making the number of mesh points to increase to provide a consistent solution with the specified accuracy. The initial number of points was $10\,$. For the tolerance $\mathsf{tol}\,\leqslant\,10^{\,-3}$, more than $10$ nodes are needed to compute the solution. For this case study, a tolerance of $\mathsf{tol}\,=\,10^{\,-3}$ would be enough to provide an accuracy to order of $\mathcal{O}\,(10^{\,-4})$.

The CPU time required to compute the solution of this case study is provided in Table 2. The compuational time required for the Euler explicit approach is also included. It has been evaluated using Matlab™ platform on a computer with a processor Intel® Core i$5$ CPU $2.80$ GHz and $15.6$ GB of RAM memory. The Euler explicit requires more time due to the CFL stability condition, since it is not possible to increase the time discretization. For the Euler implicit solution, it takes an average of $11$ iterations at each time step to converge to the solution which increases the global computational time of the method. The MOHL method is faster than the two Euler approaches implemented here. It can be remarked that for the same order of accuracy of the solution, the proposed method is two times faster. Moreover, as it can be seen in Table 1, the accuracy of the solution can be increased with a relatively low increase of the spatial mesh and of the CPU time.

One of the advantages of using the solver bvp4c is that it has as output the respective derivatives of the fields, making the computation of the fluxes straightforward. Thus, Figures 7(a) and 7(b) present the sensible and latent heat fluxes at the left and right boundaries. As illustrated in those figures, the sensible heat flux is the main component of the surface energy balance. The latent heat flux at the right boundary slightly changes with the step on the relative humidity at $t\,=\,36\ \mathsf{h}\,$, showing its little influence on the heat diffusion.

The total moisture flow (liquid plus vapour) at both boundaries are shown in Figure 8. As it can be observed, the flow follows variations of the boundary conditions. Furthermore, the left boundary has higher variations because its convective mass transfer coefficient is higher than the one at the other boundary. The step of relative humidity at the right boundary can also be observed on the moisture flow which suddenly changes.

Multilayered domain

In building constructions, multiple layers are commonly found. The configuration assumed at the interface between materials follows the hydraulic continuity ($P_{\,v,\,1}\,=\,P_{\,v,\,2}$ and $T_{\,v,\,1}\,=\,T_{\,v,\,2}\,$), which considers interpenetration of both porous structure layers [10]. Both materials are homogeneous and isotropic, and only heat and moisture transfer are simulated, through a perfectly airtight structure.

Consider Equation (2.3b) over a multidomain in which $x^{\,\star}\,\in\,[\,0,\,x_{\,\text{int}}^{\,\star}\,)\,\cup\,[\,x_{\,\text{int}}^{\,\star},\,1\,]\,$. The coefficients are written in a general form as illustrated for the moisture diffusion:

[TABLE]

where $x_{\,\text{int}}^{\,\star}$ is the location of the interface between materials and the indices $1$ and $2$ represent each material layer.

Thus, by using the MOHL approach, it is assumed that at the interface, the solution and the solution derivative of the problem are continuous. In this way, the method will search for the solution that can satisfy both conditions automatically. These interface assumptions are verified with an analytical solution in the Appendix C.

4.2.1  Application

This case study considers a porous wall formed by $2$ layers: $10$-$\mathsf{cm}$ of a load bearing material and $2$-$\mathsf{cm}$ of finishing material. Figure 9 shows schematically this physical situation. The selected materials further complicate the case, with the first layer having a faster liquid transfer while the second layer acts as a hygroscopic finish. The properties used for these materials were obtained from [19] which are given in Tables 3 and 4.

Initial conditions are considered uniform over the spatial domain, with an initial temperature of $T_{\,i}\,=\,293.15\ \mathsf{K}$ and an initial vapour pressure of $P_{\,v,\,i}\,=\,1.16\cdot 10^{\,3}\ \mathsf{Pa}\,$, referent to a relative humidity of $50\,\%\,$. The boundary conditions oscillate sinusoidally during $120$ hours of simulation, which are represented in Figure 10. The convective mass and heat transfer coefficients are set with the same values as in the previous case study.

The rain is present in the simulation between hours $[\,40,\,65\,]\ \mathsf{h}\,$, reaching a maximum value of $2\cdot 10^{\,-4}\ {\mathsf{kg}}/{\mathsf{(m^{2}\cdot s)}}$ at $52\ \mathsf{h}\,$, as presented in Figure 11(a), which generates a sensible heat flux of $12\ {\mathsf{W}}/{\mathsf{m^{2}}}$ at this boundary as shown in Figure 11(b). The rain flux is included at the left boundary, which causes a rapid increase of moisture within the material. The liquid flow from the rain can vary according to the wind speed [5]. In the literature [19, 39, 17, 22], it is possible to find values within the range $[\,2\cdot 10^{\,-5},\,8\cdot 10^{\,-4}\,]\ {\mathsf{kg}}/{\mathsf{(m^{2}\cdot s)}}\,$.

Simulations were performed using the bvp4c, with relative and absolute tolerances of the solver set to $10^{\,-5}\,$. The time is incremented with a discretization of $\Delta t^{\,\star}\,=\,10^{\,-1}\,$, with an initial guess composed by $20$ spatial nodes, for both the vapour pressure and the temperature. The main advantage of using the approach of MOHL is that there is no need to perform sub-iterations at each time step to consider the nonlinearities of the material properties. Besides, no special treatment is needed at the interface. Everything is handled automatically. For the Euler implicit solution, simulations have been performed with $\Delta x^{\,\star}\,=\,10^{\,-2}$ and $\Delta t^{\,\star}\,=\,10^{\,-1}\,$.

Figures 12(a) and 12(b) present the selected mesh for the last profile of temperature and vapour pressure fields, respectively, and also, the refined and the reference solutions. At each time iteration, the number of mesh points varies with an average of $498$ mesh points for the moisture field and $235$ mesh points for the temperature field. Due to high nonlinearities, a concentration of nodes is noted at the interface between materials and where the gradients are higher. For the domain of the second layer, with the finishing material, only $33$ nodes were used to compute the vapour pressure field, where $20$ are concentrated closer to the interface between the two materials. It indicates that the method proposed can be adapted according to the physical phenomenon.

The MOHL has demonstrated a good agreement with the reference solution given by Chebfun to reproduce the physical phenomenon. Distributions of the error $\varepsilon_{\,2}$ on function of $x$ are shown in Figure 13. The error of the MOHL is $\varepsilon_{\infty,\,v}\,=\,8.3\cdot 10^{\,-4}$ for the vapour pressure and $\varepsilon_{\infty,\,u}\,=\,10^{\,-4}$ for the temperature. These values depend on the order of the discretization in time and on the chosen tolerance of the solver.

Figure 14(b) presents the relative humidity profiles when the rain flow occurs. It is possible to observe a liquid concentration on the left edge, and, a stable situation on the right edge. This behaviour is observed due to different material properties. As the second layer is less hygroscopic than the first layer, the rain flow crosses both materials and arrives at the right boundary with little intensity. Additionally, Figure 14(a) presents the temperature profiles for the same time instants. The temperature increases in the middle of the material are due to the moisture content gradients on the storage coefficient.

The evolution of temperature, vapour pressure and relative humidity at the boundary surfaces ($x\,=\,0\ \mathsf{m}$ and $x\,=\,0.12\ \mathsf{m}\,$) are shown in Figures 15(a), 15(b) and 15(c), respectively. The vapour pressure varies according to the sinusoidal fluctuations of the boundary conditions affected by the rain flow. At $x\,=\,0\ \mathsf{m}\,$, the vapour pressure suddenly increases due to the rain flow imposed at the surface. The moisture due to the rain flow diffuses through both layers. Although, as the second layer is composed with a less hygroscopic material, the vapour pressure completely reaches this surface by $70\ \mathsf{h}\,$. In Figure 15(a), it is possible to observe a peak of temperature on the right boundary, which is caused because the heat transfer coefficient under vapour pressure gradient $k_{\,TM}^{\,\star}$ is ten times higher in the second layer than in the first one.

Figure 16(a) indicates the sensible, the latent and the total heat fluxes over both materials, at the instant $t\,=\,60\ \mathsf{h}\,$, a while after the peak of rain flow. Furthermore, Figure 16(b) presents the total moisture flow for the same time. In Material $1\,$, the sensible and latent heat fluxes have high values but with opposite signs. As the latent heat flux is slight higher they do not cancel each other. Material $2$ has a higher diffusive coefficient than Material $1\,$, which explains the horizontal line for the fluxes and flow profiles. Another important point is that the total flux and flow across the interface between the two materials are continuous, which is consistent with the assumption of hydraulic continuity [10, 11].

The moisture storage coefficient and the moisture diffusive coefficients are shown in Figures 17(a) and 17(b), respectively. They are presented on function of the space distribution for different time instant $t\,=\,\{\,40,\,50,\,60,\,70\,\}\ \mathsf{h}$. The differences are high between the material properties and it is possible to observe the non-linearities due to the moisture content gradients.

Comparison of the numerical predictions with experimental data

In this section, the physical model with the MOHL numerical model is compared with experimental data gathered from the French project HYGRO-BAT [2]. The measurements were performed at the French laboratory LOCIE (Laboratory of Optimisation of the Conception and Engineering of the Environment). One dimensional coupled heat and moisture transfer through a single-layered wall is monitored by sensors placed within the material and at its surfaces. The liquid transfer $k_{\,l}$ in the moisture balance equation — Equation (2.1a) — is neglected in this study. Surface sensors provide boundary conditions for the coupled simulation, while the other sensors provide reference measurements for the model confirmation.

The wall considered in their study is composed of a $16$-$\mathsf{cm}$ layer of wood fibre material, which is subjected to variations on the temperature and relative humidity for a $2$-week period. The material properties are given in Table 5 of the Appendix A. The reference data for the evaluation is provided by temperature and humidity sensors SHT75 Sensirion, located at $x\,=\,\{\,4,\,8,\,12\,\}\ \mathsf{cm}$ inside the wall. They are located sparsely within the wall to avoid interferences among them.

The sensors have an uncertainty of measurement of $\sigma_{\,T}^{\,\text{meas}}\,=\,0.3^{\circ}\ \mathsf{C}$ for temperature and of $\sigma_{\,\phi}^{\,\text{meas}}\,=\,0.018$ for relative humidity. In addition, the uncertainty regarding the position of the sensors is of $\sigma^{\,\text{pos}}\,=\,1\ \mathsf{cm}$ for the sensors located at $x\,=\,\{\,4,\,12\,\}\ \mathsf{cm}$ and of $\sigma^{\,\text{pos}}\,=\,0.5\ \mathsf{cm}$ for the other sensors. The uncertainty on position is higher for the sensor at $x\,=\,\{\,4,\,12\,\}\ \mathsf{cm}\,$, since according to the experiment design, they have been settled by perforating a whole in the wood fibre layer.

The temperature and relative humidity at the interior and exterior boundaries are measured and given in Figures 18(a) and 18(b), respectively. In these figures, the uncertainties related to the sensors measurements are shown in gray color. At the interior boundary, the temperature is set to approximately $24^{\circ}\ \mathsf{C}$ and the relative humidity set to $40\,\%\,$, at the first week, and to $70\,\%$ at the second week. The exterior temperature and relative humidity values are given by their measurement at the boundary, which is filtered by a $2$-$\mathsf{cm}$ coating layer, excluding solar radiation and driving rain phenomena. Thus, both boundaries are expressed as Dirichlet conditions for the model confirmation. For more details about the experiment one may refer to [38] and [40]. In this wok, the temporal derivative of $\rho_{\,w}$ in Equation (2.1b) was also taken into account, differently from [40].

For confirming simulations, the relative error - $\epsilon$ - between the solutions and the data are computed as:

[TABLE]

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

Simulations are performed on a single-layered wall of a material that separates two ambiances. The initial condition for the temperature and vapour pressure profiles are given by an interpolation of the measurements at $x\,=\,\{\,0,\,4,\,8,\,12,\,16\,\}\ \mathsf{cm}\,$, at the first time instant. The solver chosen to perform the MOHL method is the bvp4c, with relative and absolute tolerances set to $10^{\,-4}\,$. The solution has been computed with a time step of $\Delta t^{\,\star}\,=\,0.1\,$, the equivalent of $6\ \mathsf{min}\,$.

Results of the simulation are compared in terms of temperature and vapour pressure. Figure 19 presents the measured data and the simulation results in each one of the locations of the sensors, at $x\,=\,\{\,4,\,8,\,12\,\}\ \mathsf{cm}\,$. The error between the predicted solution and the experimental observations is between the uncertainties of the sensors during almost all simulation period. The uncertainties are displayed on gray in Figure 19 and they are calculated according to the position of the sensor ($\sigma^{\,\text{posi}}\,$) and according to the sensor measurement error ($\sigma^{\,\text{meas}}\,$):

[TABLE]

where $\sigma^{\,\text{total}}$ are the total uncertainties regarding the temperature and the vapour pressure fields. The detailed computation of the measurement uncertainties can be verified in [3]. For the temperature evolution, the highest difference is when the relative humidity at the left boundary change from $40\,\%$ to $70\,\%\,$. From the $7$th day, the measured value of temperature rises more than on the simulated one as can be observed in Figure 19(a). With the variations at the boundaries, the simulated temperature (v. Figures 19(c) and 19(e)) varies more than the measured temperature. It seems that the total diffusion coefficient of the temperature used for the simulation is higher than the real one. The difference between the simulated and measured values of temperature reaches a maximum of $2^{\circ}\ \mathsf{C}$ affected by the moisture flow.

A high influence of the step on the relative humidity can be observed at $x\,=\,4\ \mathsf{cm}$ in Figure 19(b), which decreases along the thickness. Simulations do not completely follow the variations of vapour pressure, but not because of the numerical method, but because of the physical model, which does not consider liquid transport neither hysteresis. In Figures  19(d) and 19(f), the differences between simulations and measurements get higher, with the incoming moisture flow reaching a maximum of $77\ \mathsf{Pa}$ and $80\ \mathsf{Pa}\,$, respectively. The absolute difference is higher at $x\ =\ 4\ \mathsf{cm}\,$, as it is closer to the left boundary and consequently to the incoming flow.

The relative error $\epsilon$ of the computed solutions are presented in Figure 20. Simulations with the MOHL method showed a good agreement with the reference data. The maximum relative error for temperature solution is $0.7\,\%$ and for vapour pressure is $5.7\,\%\,$. The average error on the dynamic profiles is $0.2\,\%$ for temperature and $1.7\,\%$ for vapour pressure, values which are close to the ones obtained by Rouchier et al.[40].

The MOHL method can provide a very accurate solution of the physical model. The difference that it is observed in Figure 19 is attributed to the model that describes the phenomenon, as previously mentioned, and to the values of the estimated properties. Nevertheless, despite these discrepancies, simulated results are considered satisfactory for the prediction of heat and moisture transport.

Conclusions

Mesh refinement is a great issue in the field of numerical simulations of thermal sciences problems. Therefore, the present article shows that the unsteady problem of heat and moisture transfer can be solved from a different perspective. Namely, we interchange the order of discretization by proposing a numerical model accurate to the order $\mathcal{O}\,\,(\Delta x^{\,4}\thicksim\Delta x^{\,6})$ with a moderate decrease of computational efforts. Traditionally, the numerical analysts discretize the differential governing equations first in space and, then, the coupled system of ODEs is solved in time. This approach is usually referred to as the Method Of Lines (MOL) [21].

Here, we proposed another point of view on the same problem, the so-called Method Of Horizontal Lines (MOHL). It consists of a semi-discretization; first, in time and only then in space using appropriate methods. So, the Initial-Boundary Value Problem (IVP–BVP) is replaced by a sequence of BVPs. Moreover, we took advantage of the fact that today there are well-tested robust adaptive numerical methods to tackle BVPs in one space dimension [43, 20]. The resulting discretization is fully implicit, thus, unconditionally stable. In other words, the semi-discretization in time is not subject to any kind of restriction regarding the time step [8]. Thus, the time step can be chosen based on the accuracy considerations solely. The distribution of spatial nodes is adaptive and can change at every time step. The grid is further refined (or unrefined) to meet a prescribed error tolerance. The numerical solution error is estimated by computing the continuous residual (inside the domain and at the boundaries).

This new method has been evaluated on two numerical case studies of heat and moisture transfer in porous media. Each case aimed at exciting the non-linear properties of the material to induce sharp profiles of temperature and vapour pressure. The first case considered a single material layer with sinusoidal variations of boundary conditions, taking into account a sudden rain effect at one boundary. The second case studied the heat and moisture diffusive transfer through a multi-layered material. For both cases, the numerical method has shown a high accuracy and perfect agreement with the respective reference solutions. The error was of the order $\mathcal{O}\,\,(10^{\,-4})$ or less. The advantage of the proposed method is the adaptive spatial mesh grid according to the solicitations of the physical phenomena. Moreover, it has been demonstrated that for the same order of accuracy of the solution, the MOHL numerical model is twice faster than the classical implicit Euler with central finite-difference approach. In the multilayered domain case study, for example, the nodes were concentrated at the interface between materials and at higher gradients, which allowed a solution with high order of accuracy.

The last case study aimed at highlighting the use of the numerical model to compare the numerical predictions with experimental data. The configuration represents a one layer wall exposed to climatic boundary conditions. On the other side, the boundary conditions are controlled so as to impose a step in relative humidity from $40\,\%$ to $70\,\%$ with a constant temperature. The comparison revealed a satisfactory agreement between the numerical results and the experimental data. The MOHL numerical model can provide a very accurate solution of the physical model to predict the heat and moisture transfer in porous materials.

As a technical improvement of the present method, there is another possible direction. Namely, one could notice that we have used a higher-order (at least, fourth order) $C^{\,1}$ continuous interpolation in space and only second-order discretization in time. It might appear as an imbalanced discretization. Hence, this point might be improved in future works by choosing higher-order differentiating in time as well. Another improvement can be achieved if we generalize the method to non-uniform time steps as well.

Acknowledgements

The authors acknowledge the Brazilian Agencies CAPES of the Ministry of Education for the financial support of this work, which was conducted during a scholarship supported by the International Cooperation Program CAPES/COFECUB (Grant $\#\,774/13$). 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]

[TABLE]

[TABLE]

Appendix A Material properties

Appendix B Dimensionless parameters

Case from Section 4.1

Problem (2.3) is considered with $\mathrm{Fo}_{\,T}\,=\,1.6\cdot 10^{\,-1}\,$, $\mathrm{Fo}_{\,M}\,=\,3.2\cdot 10^{-2}\,$, $\gamma_{\,1}\,=\,2.3\cdot 10^{-2}$ and $\gamma_{\,2}\,=\,1.58\cdot 10^{-1}\,$. The reference time is $t_{\,0}\,=\,1\ \mathsf{h}$, thus, the final simulation time is fixed to $\tau^{\,\star}\,=\,72\,$. The reference temperature and vapour pressure were taken the same as the initial condition. At the boundaries, the Biot numbers assume the following values: $\mathrm{Bi}_{\,M,\,L}\,=\,3.65$, $\mathrm{Bi}_{\,M,\,R}\,=\,0.55\,$, $\mathrm{Bi}_{\,T,\,L}\,=\,6.45\,$, $\mathrm{Bi}_{\,T,\,R}\,=\,2.06\,$, $\mathrm{Bi}_{\,TM,\,L}\,=\,0.13$ and $\mathrm{Bi}_{\,TM,\,R}\,=\,0.02\,$. The temperature and vapour pressure vary sinusoidally over the time as the following expressions:

[TABLE]

The properties are dimensionalised with the following reference values:

[TABLE]

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

[TABLE]

Case from Section 4.2

Problem (2.3) is considered with $\mathrm{Fo}_{\,T}\,=\,0.02\,$, $\mathrm{Fo}_{\,M}\,=\,0.07\,$, $\gamma_{\,1}\,=\,0.01$ and $\gamma_{\,2}\,=\,0.13\,$. The reference time is $t_{\,0}\,=\,1\ \mathsf{h}\,$, and, the final simulation time is fixed to $\tau^{\,\star}\,=\,120\,$. The reference temperature and vapour pressure were taken the same as the initial condition. At the boundaries, the Biot numbers assume the following values: $\mathrm{Bi}_{\,M,\,L}\,=\,4.4\,$, $\mathrm{Bi}_{\,M,\,R}\,=\,0.6\,$, $\mathrm{Bi}_{\,T,\,L}\,=\,6\,$, $\mathrm{Bi}_{\,T,\,R}\,=\,2\,$, $\mathrm{Bi}_{\,TM,\,L}\,=\,0.1$ and $\mathrm{Bi}_{\,TM,\,R}\,=\,0.02\,$. The ambient temperature, vapour pressure and rain flow are written as follows:

[TABLE]

The properties are scaled with the following reference values:

[TABLE]

Thus, the dimensionless properties of the Material $1$ are written as:

[TABLE]

and, the dimensionless properties of the Material $2$ are written as:

[TABLE]

Case from Section 5

Problem (2.3) is considered with $\mathrm{Fo}_{\,T}\,=\,0.06\,$, $\mathrm{Fo}_{\,M}\,=\,0.02\,$, $\gamma_{\,1}\,=\,5.17\cdot 10^{\,-3}$ and $\gamma_{\,2}\,=\,7.72\cdot 10^{\,-3}\,$. The reference time is $t_{\,0}\,=\,1\ \mathsf{h}\,$, and, the final simulation time is fixed to $\tau^{\,\star}\,=\,336\,$. The reference temperature is $T_{\,0}\,=\,293.15\ \mathsf{K}$ and the reference vapour pressure is $P_{\,v,\,0}\,=\,1166.91\ \mathsf{Pa}\,$.

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

[TABLE]

The dimensionless properties of the material can be written as:

[TABLE]

with the following reference values:

[TABLE]

Appendix C Analytical solution for the interface assumptions

Consider the simplified boundary value problem:

[TABLE]

with

[TABLE]

The problem at the left can be written as:

[TABLE]

And, the problem at right as:

[TABLE]

The analytical solution of Problems (C.3) and (C.6) are, respectively:

[TABLE]

whit the following derivatives:

[TABLE]

By assuming that the solution is continuous, $u_{\,1}(0)\ \equiv\ u_{\,2}(0)\,$, we can have the relation between the constants: $C_{\,1}\ =\ C_{\,2}\ =\ C$.

If we assumed that the flux is continuous at the interface $x\,=\,0$ we have:

[TABLE]

which leads to the following relation:

[TABLE]

Although, if we assumed that the derivative solution is continuous at the interface $x\,=\,0$ we have:

[TABLE]

which leads to the following relation:

[TABLE]

Now, we have two analytical solutions, one that considers the flux continuous, called here as Analytical-A, and the other analytical solution that consider the derivative solution continuous, called here as Analytical-B. Figure 21 presents both analytical solutions, with the Chebfun solution and the one computed with MOHL. It is shown that the Chebfun and MOHL solutions are in perfect agreement with the solution Analytical-B. This means that both numerical solutions compute the more regular solution and the flux is not continuous at the interface.