A hybrid analytical–numerical method for computing coupled temperature and moisture content fields in porous soils

Abstract

This work is devoted to proposing a hybrid numerical–analytical method to address the problem of heat and moisture transfer in porous soils. Several numerical and analytical models have been used to study heat and moisture transfer. The complexity of the coupled transfer in soils is such that analytical solutions exist only for limited problems, while numerical solutions can deal with more realistic ones but at a higher computational cost. Therefore, we propose to implement analytical solutions where variations of temperature and moisture content are known to be almost nonvarying, while the numerical solution is implemented in the remaining region, near the boundaries. The coupling between solutions is performed assuming the continuity of both fields and fluxes at each interface. This strategy allows assuring the physical phenomenon occurring at the interface. Numerical experiments are performed, showing the accuracy, the efficiency, and the great potential of the method regarding applications in nonlinear soil problems.

Keywords

Coupled heat and mass transfer hybrid method soil simulation analytical solution

Introduction

The ground may play an important role in estimating the energy consumption of low-rise buildings. Significant heat losses, from inside through the ground floor slab into soil, may imply a rise in the energy consumption (Rees et al., 2000). The heat loss via the ground is mainly calculated with pure heat transfer, as presented by Cichota et al. (2004), Claesson and Hagentoft (1991), Dos Santos and Mendes (2004), Nowamooz et al. (2015), and Rees et al. (2000). Most of these researches are based on the major hypothesis of none coupling effects between soil heat and moisture transfer. In particular, some cases consider the properties of soil as constant.

Although the soil is a porous material, the moisture transfer and its accumulation have a direct impact on the heat transfer. Moisture migration keeps the ground slab moist through the year (Rantala and Leivo, 2009). Janssen et al. (2004) showed that to precisely determine the heat losses, it requires simultaneous calculations with moisture content as they are closely related.

Other works (Bouddour et al., 1998; Janssen et al., 2002; Mulay and Worek, 1990; Rees et al., 2001) also proved that moisture content is a relevant parameter to be considered due to its mechanism of transport: as liquid water evaporates or condensates, it can absorb or release latent heat of vaporization, according to the process.

The heat and moisture transfer problem can be solved by numerical or analytical approaches, depending on the specified problem (Qin et al., 2006). The most used solution method is the numerical one due to almost none restrictions of boundary conditions (BCs), geometry, material properties, among other considerations. Examples can be found in the literature such as the finite-volume method, the finite-element method, and the finite-difference method (Berger et al., 2014; Mendes et al., 2002; Mikhailov and Özisik, 1994; Thomas et al., 1980; Yamazawa, 2001).

Concerning problems of heat and moisture transfer in porous soils, Liu et al. (2005) presented a 1D mathematical model based on the volume-averaging method. The authors validated their numerical method with an experiment conducted under natural environmental conditions. Dos Santos and Mendes (2006) presented a complete analysis of soil simulation with building interactions. They applied a numerical method to compute the 3D heat and moisture transfer in soils, with the finite-volume method. Their work included the effects of rain on the heat losses from the ground to the building. Coelho et al. (2009) developed a succeeded approach based on neural network and clustering method to compute the temperature and moisture fields in unsaturated soils. Mendes et al. (2012) compared in situ soil measurements with two different numerical models, the coupled Philip and De Vries model prescribed in Domus (Mendes et al., 2008) and a decoupled one, in order to validate them. Both numerical models succeed in predicting the temperature and the moisture content profiles.

Analytical methods as the ones proposed in Abahri et al. (2011), Chang and Weng (2000), Qin et al. (2006), and Ribeiro et al. (1993) have some disadvantages such as simplifications of the real physical phenomenon and the fact that variable properties and transport coefficients have to be considered linear. Nevertheless, the analytical solution has important advantages (Cichota et al., 2004) such as the reduction of the computational cost, especially when the variations of temperature and moisture content profiles must be investigated over significantly long periods (Lee et al., 2010).

Even with the development of these methods, modeling the coupled phenomenon of heat and moisture transfer in soil faces important issues. Indeed, the computational cost rises importantly when computing the solution due to the following problems among others:

The definition of the BCs at the bottom surface of the soil domain is a difficult question. Most of the works (Coelho et al., 2009; Dos Santos and Mendes, 2006; Janssen et al., 2004; Thomas et al., 1980) consider an adiabatic condition, assuming that after a certain depth, there are no temperature and moisture content gradients. This hypothesis imposes the modeling of a large spatial domain (around 10/20 m), whereas the interesting dynamics region lies at the first decimeters (Rantala, 2005).

The consideration of moisture transfer imposes strong nonlinearities (Dos Santos and Mendes, 2006; Thomas et al., 1980). Contrarily to the analytical approach, numerical methods can consider the moisture dependency of physical properties. Depending on the type of numerical approach, one requires very short time steps to compute the solution with sufficient accuracy. Thus, for long-time simulation (for one or several years), considering large spatial domain, the computational cost increases.

Thus, to deal with all these issues, traditional numerical methods have to be improved or even replaced by other models (Pallin and Kehrer, 2013). The objective of this work is to propose an innovative mathematical method to solve nonlinear heat and moisture transfer in soils, with a special focus on the application to large domains. A hybrid method is proposed based on the combination of analytical and numerical approaches. This type of approach was already applied by Godinho et al. (2011) for solving the problem of underwater sound scattering by an elastic shell. Lee et al. (2010) proposed a similar approach, dealing only with analytical solutions, for solving heat and moisture transfer through a two-layered wall.

As a matter of fact, it can be noted that temperature and moisture content gradients are greater near the surface, and deep into the ground, temperature and moisture distributions are almost uniform (Coelho et al., 2009). This behavior provides a good opportunity for the hybrid method. The strength of the numerical models can be used near the surface, to perform accurate computations where the effects of the nonlinearities are more intense. However, the analytical solution can be used for the deeper layers where the diffusion process is relatively smooth.

In this work, section “Coupled heat and moisture transfer in soils” details the mathematical model of heat and moisture transfer in unsaturated porous soils. Fundamentals of the hybrid method are presented in section “Hybrid analytical–numerical method.” Numerical results are performed and discussed in section “Computational simulation.” Then, in section “Final remarks,” conclusions are addressed and future perspectives are provided.

Coupled heat and moisture transfer in soils

Mathematical formulation

Dos Santos and Mendes (2006) adopted the Philip and De Vries equations for the heat and moisture transfer in soils, with $θ$ and T as driving potentials. Based on them, the energy conservation equation is written as

ρ_{0} c_{m} \frac{\partial T}{\partial t} = \nabla [(λ + L_{v} ρ_{ℓ} D_{Tv}) \nabla T + L_{v} ρ_{ℓ} D_{θ v} \nabla θ]

(1)

and the mass conservation as

\frac{\partial θ}{\partial t} = \nabla [D_{T} \nabla T + D_{θ} \nabla θ]

(2)

where $ρ_{0}$ is the solid matrix density (kg/m³), $c_{m}$ the specific heat (J/kg/K), T the temperature (K), $λ$ the thermal conductivity (W/m/K), $L_{v}$ the latent heat of vaporization (J/kg), $θ$ the volumetric moisture content (m³/m³), $ρ_{ℓ}$ the liquid water density (kg/m³), $D_{Tv}$ the vapor-phase transport coefficient associated with a temperature gradient (m²/s/K), and $D_{θ v}$ the vapor-phase transport coefficient associated with a moisture content gradient (m²/s). Coefficients $D_{T}$ and $D_{θ}$ are the mass transport coefficients associated with the temperature gradient (m²/s/K) and the moisture content gradient (m²/s), respectively, which are related only with the vapor-phase transport, $D_{Tv}$ and $D_{θ v}$ . The properties $λ$ , $D_{Tv}$ , $D_{θ v}$ , $D_{T}$ , and $D_{θ}$ are functions of temperature and moisture content $(T, θ)$ (See notation in Appendix 1).

The BCs at the surface are considered to be in contact with the internal air of the building and can be mathematically described by the following equations (Dos Santos and Mendes, 2006; Mendes et al., 2012)

{((λ + L_{v} ρ_{ℓ} D_{Tv}) \frac{\partial T}{\partial y} + L_{v} ρ_{ℓ} D_{θ v} \frac{\partial θ}{\partial y})}_{y = 0} = h_{T} (T_{y = 0} - T_{\infty}) + L_{v} h_{m} (ρ_{v, y = 0} - ρ_{v, \infty})

(3a)

{(D_{T} \frac{\partial T}{\partial y} + D_{θ} \frac{\partial θ}{\partial y})}_{y = 0} = \frac{h_{m}}{ρ_{ℓ}} (ρ_{v, y = 0} - ρ_{v, \infty})

(3b)

with $ρ_{v}$ as the water vapor density (kg/m³), $ρ_{v, \infty}$ the water vapor of the surrounding air, $h_{T}$ the convective heat transfer coefficient (W/m²/K), and $h_{m}$ the convective mass transfer coefficient (m/s). The term $(ρ_{v} - ρ_{v, \infty})$ can be written as a linear combination of temperature and moisture content (Dos Santos and Mendes, 2006; Mendes et al., 2002)

(ρ_{v} - ρ_{v, \infty}) = M_{1} (T - T_{\infty}) + M_{2} (θ - θ_{\infty}) + M_{3}

(4)

where

M_{1} = A \frac{M}{ℜ} ϕ_{\infty}

(5a)

M_{2} = \frac{M}{ℜ} {(\frac{P_{s}}{T})}_{y = 0}^{prev} {(\frac{\partial ϕ}{\partial θ})}_{y = 0}^{prev}

(5b)

M_{3} = \frac{M}{ℜ} [{(\frac{P_{s}}{T})}_{y = 0}^{prev} \cdot R (θ_{y = 0}^{prev}) + ϕ_{\infty} \cdot (R (T_{\infty}) - R (T_{y = 0}^{prev}))]

(5c)

with

\frac{P_{s}}{T} = - 3813.8 + 41.7635 \cdot T - 0.1531 \cdot T^{2} + (1.87 \cdot 10^{- 4}) \cdot T^{3}

(6a)

R_{1} (T) = \frac{P_{s}}{T} - (A_{1} T + B_{1})

(6b)

R_{2} (θ) = f (θ) - (A_{2} θ + B_{2})

(6c)

In these equations, the subscript $\infty$ represents the surrounding air, R is a residual function, $P_{s}$ is the saturated pressure (Pa), ℜ is the universal gas constant (J/mol/K), M is the molecular mass (kg/mol), $ϕ$ is the relative humidity, and ‘prev’ is the previous iteration. $A_{1}$ and $B_{1}$ are the straight-line coefficients from the approximation $(P_{s} / T = A_{1} T + B_{1})$ , their values correspond to $A_{1} = 0.5874 (Pa)$ and $B_{1} = - 1619708 (Pa / K)$ . Similarly, $R (θ_{y = 0}^{prev})$ is the residue calculated by the sorption isotherm, evaluated at the previous iteration. Its approximation coefficients are calculated according to the material. The other constants admit the respective values: $ρ_{ℓ} = 1000 (kg / m^{3})$ , $M = 0.018 (kg / mol)$ , $ℜ = 8.314 (J / K / mol)$ , and $L_{v} = 2.5 \cdot 10^{6} (J / kg)$ . $T_{\infty}$ , $ϕ_{\infty}$ , $h_{T}$ , and $h_{m}$ are given by the problem.

On the other side, an adiabatic condition is assumed

{((λ + L_{v} ρ_{ℓ} D_{Tv}) \frac{\partial T}{\partial y} + L_{v} ρ_{ℓ} D_{θ v} \frac{\partial θ}{\partial y})}_{y = L} = 0

(7a)

{(D_{T} \frac{\partial T}{\partial y} + D_{θ} \frac{\partial θ}{\partial y})}_{y = L} = 0

(7b)

At the initial condition, the temperature and the moisture content values are

T (y, 0) = T_{0} (y)

(8a)

θ (y, 0) = θ_{0} (y)

(8b)

These values are chosen according to the problem. They can be constant values or functions of y. Other conditions will be necessary to the use of the hybrid method and need to be presented. The Neumann condition gives

{((λ + L_{v} ρ_{ℓ} D_{Tv}) \frac{\partial T}{\partial y} + L_{v} ρ_{ℓ} D_{θ v} \frac{\partial θ}{\partial y})}_{y = i} = - q_{T}

(9a)

{(D_{T} \frac{\partial T}{\partial y} + D_{θ} \frac{\partial θ}{\partial y})}_{y = i} = - q_{m}

(9b)

with $q_{T}$ as the heat flux density (W/m²) and $q_{m}$ as the moisture flux (m/s). And the Dirichlet condition

T (y = i, 0) = T_{i}

(10a)

θ (y = i, 0) = θ_{i}

(10b)

The dimensionless representation

This subsection presents the dimensionless form of the linear mathematical model. To adimensionalize the problem, the following dimensionless quantities are defined:

$χ = \frac{y}{L}$	Dimensionless coordinate
$τ = \frac{t}{t^{*}}$	Dimensionless time
$u (χ, τ) = \frac{T - T_{b}}{T_{a} - T_{b}}$	Dimensionless temperature
$v (χ, τ) = \frac{θ - θ_{b}}{θ_{a} - θ_{b}}$	Dimensionless moisture content
$Fo = \frac{α t^{*}}{L^{2}}$	Fourier number
$Ko = r \frac{θ_{a} - θ_{b}}{T_{a} - T_{b}}$	Kossovich number
$Lu = \frac{D_{θ}}{α}$	Luikov number
$Pn = δ \frac{T_{a} - T_{b}}{θ_{a} - θ_{b}}$	Posnov number
$B i_{m} = \frac{h_{m} L}{D_{θ}}$	Dimensionless convective mass transfer coefficient (Biot)
$B i_{T} = \frac{h_{T} L}{λ}$	Dimensionless convective heat transfer coefficient (Biot)
$Q_{T} = \frac{q_{T} L}{λ (T_{a} - T_{b})}$	Dimensionless heat flux
$Q_{m} = \frac{q_{m} L}{D_{θ} (θ_{a} - θ_{b})}$	Dimensionless moisture flux

where $T_{a}$ , $T_{b}$ , $θ_{a}$ , and $θ_{b}$ represent constant temperature and moisture values chosen to define the dimensionless variables. $T_{a}$ and $θ_{a}$ are selected to be the highest value of the temperature and moisture content, respectively, of one cycle. Even so, $T_{b}$ and $θ_{b}$ are the mean values of temperature and moisture content of one cycle, respectively. The characteristic time and the characteristic length are defined as $t^{*}$ and $L$ . Coefficient $δ = D_{Tv} / D_{θ v}$ represents the thermogradient coefficient, $α = λ / ρ_{0} c_{m}$ represents the thermal diffusivity, and $r = L_{v} ρ_{ℓ} / ρ_{0} c_{m}$ .

The dimensionless problem of equations (1) and (2) is written as (Mikhailov and Özisik, 1994)

\frac{\partial u}{\partial τ} = (Fo + Fo \cdot Ko \cdot Lu \cdot Pn) \frac{\partial^{2} u}{\partial χ^{2}} + Fo \cdot Ko \cdot Lu \frac{\partial^{2} v}{\partial χ^{2}}

(11a)

\frac{\partial v}{\partial τ} = Fo \cdot Lu \cdot Pn \frac{\partial^{2} u}{\partial χ^{2}} + Fo \cdot Lu \frac{\partial^{2} v}{\partial χ^{2}}

(11b)

Considering the subsequent BCs at the surface

{\frac{\partial u}{\partial χ} |}_{χ = 0} = B i_{T} (u - u_{\infty})

(12a)

{\frac{\partial v}{\partial χ} |}_{χ = 0} = (B i_{m} {\bar{M}}_{1} - PnB i_{T}) (u - u_{\infty}) + {\bar{M}}_{2} B i_{m} (v - v_{\infty}) + B i_{m} {\bar{M}}_{3}

(12b)

where ${\bar{M}}_{1}$ , ${\bar{M}}_{2}$ , and ${\bar{M}}_{3}$ are new dimensionless numbers

{\bar{M}}_{1} = \frac{M_{1}}{ρ_{ℓ}} \frac{T_{a} - T_{b}}{θ_{a} - θ_{b}}

(13a)

{\bar{M}}_{2} = \frac{M_{2}}{ρ_{ℓ}}

(13b)

{\bar{M}}_{3} = \frac{M_{3}}{ρ_{ℓ} (θ_{a} - θ_{b})}

(13c)

At the bottom, the dimensionless forms are

{\frac{\partial u}{\partial χ} |}_{χ = ℓ} = 0

(14a)

{\frac{\partial v}{\partial χ} |}_{χ = ℓ} = 0

(14b)

Initial conditions, represented by equations (8a) and (8b), are now written as follows

u (χ, 0) = 0

(15a)

v (χ, 0) = 0

(15b)

The other BCs in the nondimensional form are

{\frac{\partial u}{\partial χ} |}_{χ = i} = - Q_{T} + Lu \cdot Ko \cdot Q_{m}

(16a)

{\frac{\partial u}{\partial χ} |}_{χ = i} = - Pn \cdot Q_{T} - Q_{m}

(16b)

u = u_{i}

(16c)

v = v_{i}

(16d)

Hybrid analytical–numerical method

The hybrid method proposed here consists in decomposing the domain of the differential governing equations in two subdomains, separated by an interface. The location of this interface is not randomly chosen, and its location is determined by a pre-simulation of the problem, with the numerical method. After this, each one of the subdomains is solved analytically or numerically, developing their solution procedures alternately. Thus, through an iterative process, the approximation of the original problem is achieved by a composition of both subdomain solutions (Quarteroni and Valli, 1999; Shanthikumar and Sargent, 1983).

Analytical method

There are several methods of analytical solutions that have been applied to solve the heat and mass transfer problem (Abahri et al., 2011; Chang and Weng, 2000; Lee et al., 2010), as the separation of variables, the integral transforms, the Laplace transform, and others techniques. The one applied in this work is the separation of variables (Özisik, 1993), which can be found in Appendix 2. This approach is used due to the use of the hybrid method, which requires the update of the initial condition of the solution at each increment of time. The analytical method is particularly interesting for its robustness, which allows reducing the number of spatial and temporal nodes without losing accuracy. There is almost no restriction of stability.

Numerical method

The numerical method applied in this study is based on the finite-difference approach, with an implicit second-order scheme, of the Crank–Nicolson type (Crank and Nicolson, 1947). It is used as part of the hybrid method and as a reference solution to validate the features of the hybrid approach.

Illustration

To clarify the hybrid method, consider the domain $Ω$ over equations (11a) and (11b) decomposed into two parts, $Ω_{1}$ and $Ω_{2}$ as represented in Figure 1. In this case, BCs 0 and 1 are of the Robin type and Neumann type, respectively. Subdomains $Ω_{1}$ and $Ω_{2}$ do not overlap one another. The purpose of the domain decomposition is to solve one of the subdomains ( $Ω_{1}$ ) numerically, whereas the second subdomain ( $Ω_{2}$ ) is solved analytically.

Figure 1.

Schematic diagram of the domain decomposition in two parts, where N corresponds to the Neumann condition and D to the Dirichlet condition at the interface (i). Besides, the BC corresponds to the boundary condition at outer surfaces with 0 and 1 corresponding to each side.

It is important to mention that material properties are considered constant in subdomain 2 due to the analytical solution but nonlinear in subdomain 1 due to the numerical solution. This separation will be defined according to the variation of the material properties of the soil (see “Introduction”).

When the domain is decomposed, an interface between subdomains is created (called hybrid interface henceforth), and (i) the continuity of the field and (ii) the conservation of the flux must be ensured (Özisik, 1993). The mathematical formulation of these assumptions for the temperature is given by (De Freitas et al., 1996)

{\begin{matrix} u_{1, χ = i^{-}} = u_{2, χ = i^{+}} \\ {λ_{1} \nabla_{n} u_{1} |}_{χ = i^{-}} = {λ_{2} \nabla_{n} u_{2} |}_{χ = i^{+}} \end{matrix}

(17)

Superscripts − and + represent the left and the right sides of the interface, respectively. The continuity of the field and the conservation of the flux are also assumed for the moisture content.

To satisfy equation (17), at one side of the interface, a Neumann BC is assumed. Then, on the other side, a Dirichlet BC is considered. It would be more common to assume a Robin condition at the interfaces. Although, it is not employed here because of the eigenvalue problem originated from the analytical solution. It yields a transcendental equation whose resolution cannot be systematically computed for any particular value of the BC without degrading seriously the computational advantage of using an analytical solution.

Thus, taking into account the hybrid interface, the boundaries of the $Ω_{1}$ domain are written as follows

B C_{0} for χ = 0

(18a)

{\frac{\partial u_{1}}{\partial χ} |}_{χ = i} = - ψ_{u, i}, for χ = i

(18b)

{\frac{\partial v_{1}}{\partial χ} |}_{χ = i} = - ψ_{v, i}, for χ = i

(18c)

where

ψ_{u, i} = {Q_{T, 2} |}_{i^{+}} - Lu \cdot Ko \cdot {Q_{m, 2} |}_{i^{+}}

(19a)

ψ_{v, i} = {Pn \cdot Q_{T, 2} |}_{i^{+}} + {Q_{m, 2} |}_{i^{+}}

(19b)

And the boundaries of $Ω_{2}$ are defined as

u_{2} (χ, τ) = u_{i}, for χ = i

(20a)

v_{2} (χ, τ) = v_{i}, for χ = i

(20b)

B C_{1} for χ = 1

(20c)

Subscripts 1 and 2 represent subdomains $Ω_{1}$ and $Ω_{2}$ . As, in this example, we are working with only one material, the continuity of the moisture field can be assured for equation (20b), without considering the discontinuity of the moisture content profile.

Since $ψ_{u, i}$ , $ψ_{v, i}$ , $u_{i}$ , and $v_{i}$ are unknowns, the algorithm executes iterations within every time step. These conditions must be initialized with a first guess. A natural choice for the first guess is given by the initial conditions; as we are working with incrementation on time, the choice is the previous value. If the first guess is initiated with a random value, then the number of iterations may increase severally, or even do not converge.

u_{i}^{k = 0} = u_{2} (χ = i, τ = m - 1)

(21a)

v_{i}^{k = 0} = v_{2} (χ = i, τ = m - 1)

(21b)

ψ_{u, i}^{k = 0} = {\frac{\partial u_{1}}{\partial χ} |}_{χ = i}^{m - 1}

(21c)

ψ_{u, i}^{k = 0} = {\frac{\partial u_{1}}{\partial χ} |}_{χ = i}^{m - 1}

(21d)

where $k = 1, 2, 3, \dots$ is the number of iterations and $k = 0$ is the first guess value; $m - 1$ represents the value at the previous time.

For the first time step, the algorithm runs the numerical solution considering the first guess of the flux at the hybrid interface (equation (21c)). Once the temperature and moisture content profiles are computed for $Ω_{1}$ , its values at the interface $u_{i}^{k}$ and $v_{i}^{k}$ are reinitialized. Then, the analytical solution is computed for $Ω_{2}$ , considering the value just computed of the temperature and moisture content at the interface. It returns the temperature and moisture content profile for $Ω_{2}$ and a new value for fluxes at the interface $ψ_{u, i}^{k}$ and $ψ_{v, i}^{k}$ . Then, the numerical method runs once again considering an updated BC at the hybrid interface. Both methods keep iterating until the convergence at the interface is reached

\max [| u_{i}^{k} - u_{i}^{k - 1} |; | v_{i}^{k} - v_{i}^{k - 1} |; | ψ_{u, i}^{k} - ψ_{u, i}^{k - 1} |; | ψ_{v, i}^{k} - ψ_{v, i}^{k - 1} |] < ε

(22)

where $ε$ is a fixed tolerance. This tolerance cannot be smaller than the order of error of the numerical solution; otherwise, it will not converge. When convergence is reached, the algorithm goes forward to the next physical time step $τ : = τ + d τ$ .

The hybrid method algorithm can be synthesized in Algorithm 1.

Algorithm 1

Solution of the hybrid method for the two subdomains.

whileτ < total timedo

Initializing of $u_{i}^{k = 0}$ ; $v_{i}^{k = 0}$ ; $ψ_{u, i}^{k = 0}$ and $ψ_{v, i}^{k = 0}$ ;

while $\max [u_{i}^{k} - u_{i}^{k - 1}; v_{i}^{k} - v_{i}^{k - 1}; ψ_{u, i}^{k} - ψ_{u, i}^{k - 1}; ψ_{v, i}^{k} - ψ_{v, i}^{k - 1}] < ε$ do

Step 1: Run numerical method in Ω₁ with $ψ_{u, i}^{k = 0}$ and $ψ_{v, i}^{k = 0}$ ;

Returns $u_{i}^{k}$ and $v_{i}^{k}$ ;

Step 2: Run analytical method in Ω₂ with $u_{i}^{k}$ and $v_{i}^{k}$ ;

Returns $ψ_{u, i}^{k}$ and $ψ_{v, i}^{k}$ ;

Interface update: $u_{i}^{k = 0} : = u_{i}^{k}$ , $v_{i}^{k = 0} : = v_{i}^{k}$ , $ψ_{u, i}^{k = 0} : = ψ_{u, i}^{k}$ , and $ψ_{v, i}^{k = 0} : = ψ_{v, i}^{k}$ ;

end

Time increment: $τ : = τ + d τ$ ;

Fields update: $u^{m - 1} : = u^{m}$ and $v^{m - 1} : = v^{m}$ ;

end

Computational simulation

Case study

To study the features of the hybrid method for heat and moisture transfer in soils, a unidimensional case study of 20 m length is considered. Figure 2 presents the geometric configuration of the problem to be simulated. The dry-basis properties of the sandy–silt soil are presented in Table 1, from Oliveira et al. (1993). The nonlinear properties can be found in Table 3.

Figure 2.

Schematic representation of the ground problem.

Table 1.

Dry-basis properties of the soil.

Soil	$ρ_{0} (m^{3} / kg)$	$c_{m} (J / kg / K)$	$λ_{m} (W / m / K)$
Sandy–silt	1280	880	0.3452

The surface, at y = 0, is subjected to Robin BCs, mathematically represented by equation (3). For the surface conditions, it was considered a yearly average temperature of 288.15 K with a daily and a yearly variation of 5K. For the relative humidity, it was considered a daily variation between 60% and 80% (Dos Santos and Mendes, 2006; Mendes et al., 2012). Summarizing it

T_{\infty} = 288.15 + 5 \sin (2 π t / 31, 536, 000) + 5 \sin (2 π t / 86, 400) (K)

(23a)

ϕ_{\infty} = 0.7 + 0.1 \sin (2 π t / 86, 400)

(23b)

h_{T} = 5 W / m^{2} / K

(23c)

h_{m} = 5 \cdot 10^{- 3} m / s

(23d)

The soil is initially at a constant temperature of $T_{0} = 288.15 K$ and with a moisture content of $θ_{0} = 0.061 m^{3} / m^{3}$ . At the bottom (y = 20 m), into the ground, a Neumann adiabatic condition is assumed, equation (9).

Adimensionalization

To write the problem in the nondimensional form, some characteristic values must be adopted. The characteristic time is $t^{*} = L^{2} / a$ ; in the linear case, this implies that $Fo = 1$ . The characteristic length is equivalent to half of the domain L = 10 m, and the temperature and moisture content references are $T_{a} = 298.15 K$ , $T_{b} = 288.15 K$ , $θ_{a} = 0.0713 m^{3} / m^{3}$ and $θ_{b} = 0.061 m^{3} / m^{3}$ . The following values of boundary and initial conditions are resulting from the adimensionalization process

u_{\infty} = 0.5 \sin (\frac{2 π τ}{τ_{year}}) + 5 \sin (\frac{2 π τ}{τ_{day}})

(24a)

v_{\infty} = \sin (\frac{2 π τ}{τ_{day}})

(24b)

u_{0} = 0

(24c)

v_{0} = 0

(24d)

$τ_{day}$ represents the daily period and $τ_{year}$ the yearly period of variation of the BC, respectively. These values will not change for every simulation, while the specific values of the dimensionless numbers are given when the results are presented.

Simulation procedure

The simulation procedure performed by the hybrid method is divided into two steps: (i) perform the pre-simulation to chose the hybrid interface and (ii) simulate the hybrid analytical–numerical method to compute the fields.

In the framework of the hybrid method, the whole domain is divided into two parts, subdomains $Ω_{1}$ and $Ω_{2}$ , separated by an interface, as illustrated in Figure 2. As previously mentioned, to determine the location of the hybrid interface a pre-simulation is performed considering linear properties. This was the way found to compute a location for the interface without spending much computational time. To generate the worst scenario, the constant properties considered are their maximum values inside the range of relative humidity cycle, where the process of diffusion will reach the most deeper as possible. The interface location is then selected according to the moisture penetration depth $δ_{v}$ , because the properties vary according to this potential. For the second step, subdomain $Ω_{1}$ is computed considering the nonlinearities of the properties and the Robin condition, while subdomain $Ω_{2}$ is computed with a simplification of the problem, using the analytical method.

The linear pre-simulation is performed with the numerical method presented in section “Numerical method.” The dimensionless numbers are Ko = 2.3, Lu = 0.01, Fo = 1, and Pn = 13. These values are obtained from the maximum property values (Table 3) referent to the range of the relative humidity. Biot numbers are assumed to be $B i_{T} = 1$ and $B i_{m} = 10$ . The dimensionless time $τ = 0.1$ is equivalent to 1 year of simulation, and, $τ = 2.6 \cdot 10^{- 4}$ and $τ = 10^{- 5}$ are equivalent to 1 day and to 1 h, respectively. Concerning the dimensionless discretization of the numerical method, the time step $d τ$ is equivalent to 1 h of the dimensionless time and the spatial discretization to $d χ = 5 \cdot 10^{- 4}$ .

Some additional simulations are executed, considering the linear case and the nonlinear, just to the necessity of considering the nonlinearities. With these, it is possible to understand how well the hybrid method is suitable in this case. The difference between the linear and the nonlinear pre-simulation is that the dimensional numbers Ko, Lu, Pn, and Fo are varying at each time step according to the moisture content.

Numerical results and discussion

The importance of considering nonlinear properties

In Figure 3, the differences for the case study between a pre-simulation with constant properties and nonlinear ones are highlighted. It displays the evolution of moisture content over the time, presented on the soil surface ( $χ = 0$ ), where it is known to present stronger nonlinearities (Dos Santos and Mendes, 2006). It is evident that this case cannot be simulated considering constant properties, and the root means square error is 0.8503, which means both models are very different.

Figure 3.

Evolution of the moisture content at the soil surface.

In addition, Figure 4 compares the linear with the nonlinear case but presenting the dimensionless moisture as a function of the space. The nonlinearities are stronger at the surface in contact with air, that is, $χ = 0$ . Nevertheless, it can be noted that after a certain depth, results of the linear case are similar to the nonlinear case.

Figure 4.

Moisture content profile for $τ = 0.5$ dimensionless time of simulation, comparing the linear and the nonlinear cases.

Pre-simulation: definition of the hybrid interface

This type of behavior is particularly suitable for the hybrid method. The numerical part will be implemented in the part where the nonlinearities are stronger. As noted in Figure 4, the nonlinearities are stronger near the surface. The problem is the location of the interface that will separate both subdomains. For this, a pre-simulation is performed with the numerical method, for a dimensionless time $τ = 0.5$ , considering constant properties. The values of the properties used are the highest possible inside the range of the BCs.

Given the results of the pre-simulation, the mass and temperature fields and the fluxes are analyzed. First, it seeks where the fluxes are near to zero (with a given tolerance) and then where the second derivatives of the fields are close to zero

\arg_{χ \in [0, ℓ]} ({\frac{\partial v}{\partial χ} |}_{χ = (0, ℓ)} \leq ε_{1}, {\frac{\partial^{2} v}{\partial χ^{2}} |}_{χ = (0, ℓ)} \leq ε_{2})

(25)

If equation (25) is satisfied, then the location $χ = χ_{i}$ is chosen to be the interface. However, if there exist more than one location where equation (25) is satisfied, then the location selected is the one farthest from the exterior surface. The same criterion is applied for field u. As the properties only depend on v, the interface for moisture and temperature fields is the same. Here, a tolerance $ε = 10^{- 2}$ is chosen. There is a restriction on these tolerances. Their values cannot be lower than the order of error of the numerical solution $O = (d χ^{2} + d τ^{2})$ . If this occurs, then the convergence will never be achieved.

Figure 5 presents the location where the condition given by equation (25) is satisfied as a function of time. There are two logical values to be selected, the highest value $δ_{v}^{\max} = 0.9273$ or the average $δ_{v}^{mean} = 0.7036$ . The temperature and moisture content have different penetration depths, depending on the Luikov number and the moisture content. In this case, the temperature penetrates more deeply than the moisture content because these parameters have low values, the Luikov number is on order of $10^{- 2}$ .

Figure 5.

Penetration depth of temperature and moisture content gradients.

Comparison of hybrid and numerical methods

Using previous results, the interface location is first set as $δ_{v}^{\max} = 0.9273$ . The simulation using the hybrid method is performed for a dimensionless time of $τ = 1$ , equivalent to 10 years. A perfect agreement can be seen, in Figures 6 and 7, of the hybrid and the reference finite-difference methods. As it can be observed, the temperature variations penetrate into the ground more deeply than the moisture content. Even with the hybrid model, simplified with an analytical solution in one part, it does not imply lower accuracy. It can be seen that after a certain depth, there is no flux arriving from the surface, justifying the use of the hybrid method. Although the pre-simulation does not ensure a good location of the hybrid interface, it can be noted that the temperature flux is very low after $χ = 0.5$ and the pre-simulation pointed to a penetration around $χ = 1.2$ . The same behavior is observed for the moisture content. One thing that can be improved is the criteria to find the location of the interface, by restricting the tolerance or by using a nonlinear model as pre-simulation.

Figure 6.

Temperature profiles for different time instants.

Figure 7.

Moisture content evolution at different depths, $χ = {0, 0.025, 0.1, 0.5}$ .

If the pre-simulation is considered on the computational cost of the hybrid method then, in this case, 20s must be added to the final simulation time. More details of the computational cost are provided in section “Computational cost.”

Figure 8 presents the root mean squared error (RMSE), calculated among the time for each node of the grid. This error is given by the following formulation (Cichota et al., 2004)

RMSE = \sqrt{\frac{1}{N} \sum_{j = 1}^{N} (Y_{j}^{hyb} (x, t) - Y_{j}^{ref} (x, t))}

(26)

where $Y_{j}^{hyb}$ is the predicted value, which is calculated by the hybrid method and $Y_{j}^{ref}$ is the reference value, determined by the numerical method.

Figure 8.

Error computed for temperature (u) and moisture content (v), considering $δ_{v}^{\max}$ .

For the hybrid simulation with the interface located at $δ_{v}^{\max}$ , it can be seen in Figure 8 that the highest differences in RMSE values are near the interface and near the upper surface. Although for the interface, these differences are acceptable knowing that tolerance at the interface (equation (22)) is in the order of $10^{- 3}$ and their values are within this range, for the hybrid simulation with the interface located at $δ_{v}^{mean}$ , the error increases at this location because there are higher gradients arriving around this point, as can be observed in Figure 9. Even with the increase in the error, the results are acceptable for both locations of the interface because they are in the order of $10^{- 4}$ and lower than the order of the numerical scheme.

Figure 9.

Error computed for temperature (u) and moisture content (v), considering $δ_{v}^{mean}$ .

In Figure 10, the error is presented as a function of time for the case where the interface is located at $δ_{v}^{mean}$ . It varies cyclically according to the incoming heat and moisture fluxes originated by the BCs. It is still increasing because temperature and moisture content profiles have not reached the steady periodical solution; for better results, the total simulation time must be extended. However, as previously mentioned, the errors are considerably low, in the order of $10^{- 4}$ .

Figure 10.

Error in function of time computed for temperature (u) and moisture content (v), considering $δ_{v}^{mean}$ .

A numerical model using adaptive mesh has a similar approach. The mesh discretization is more refined where the nonlinearities are higher at the surface, whereas the mesh size increases deep in the ground. Nevertheless, the analytical solution may not have restrictions regarding stability or accuracy.

Computational cost

As previously mentioned, the hybrid method is faster than the traditional Crank–Nicolson scheme, when they are directly compared as shown in Table 2 for the CPU time. Simulations were performed in a computer with a processor Intel^® Core^™ i5-4460, with 8 GB of RAM memory, and coded in MATLAB (2013). If we add the computational time of the pre-simulation, the total gain with the hybrid method drops from 16% to 12%.

Table 2.

Computational time required for the nonlinear case.

Numerical scheme	CPU time (s)	CPU time (%)
Hybrid	361.4	84
Crank–Nicolson	430.4	100

Table 3.

Properties of the sandy–silt soil.

$θ (m^{3} / m^{3})$	$ϕ (-)$	$λ (W / m / K)$	$D_{θ} (m^{2} / s)$	$D_{T} (m^{2} / s / K)$	$D_{θ v} (m^{2} / s)$	$D_{Tv} (m^{2} / s / K)$
0	4.00E−03	3.00E−01	0	0	1.63E−10	9.50E−14
2.59E−02	1.96E−01	3.10E−01	0	0	2.19E−09	3.71E−12
5.17E−02	6.10E−01	3.30E−01	1.10E−31	4.70E−36	2.10E−09	1.16E−11
7.76E−02	8.60E−01	3.50E−01	3.26E−11	1.38E−15	9.04E−10	1.70E−11
1.034E−01	8.68E−01	3.90E−01	2.43E−09	9.26E−14	3.34E−10	1.97E−11
1.293E−01	8.76E−01	4.60E−01	1.30E−08	7.46E−13	7.70E−11	2.13E−11
1.551E−01	8.85E−01	6.00E−01	4.32E−08	3.28E−12	2.65E−11	2.23E−11
1.810E−01	8.93E−01	8.10E−01	1.09E−07	1.03E−11	1.14E−11	2.33E−11
2.068E−01	9.01E−01	1.11E+00	2.35E−07	2.64E−11	5.67E−12	2.42E−11
2.327E−01	9.09E−01	1.40E+00	4.52E−07	5.83E−11	3.09E−12	2.52E−11
2.585E−01	9.18E−01	1.59E+00	8.03E−07	1.16E−10	1.80E−12	2.61E−11
2.844E−01	9.26E−01	1.66E+00	1.34E−06	2.12E−10	1.10E−12	2.70E−11
3.102E−01	9.34E−01	1.67E+00	2.17E−06	3.65E−10	6.98E−13	2.78E−11
3.361E−01	9.42E−01	1.67E+00	3.41E−06	5.96E−10	4.53E−13	2.87E−11
3.619E−01	9.51E−01	1.68E+00	5.27E−06	9.29E−10	2.98E−13	2.95E−11
3.878E−01	9.59E−01	1.68E+00	8.12E−06	1.39E−09	1.97E−13	3.03E−11
4.136E−01	9.67E−01	1.68E+00	1.26E−05	2.02E−09	1.29E−13	3.10E−11
4.395E−01	9.75E−01	1.68E+00	2.04E−05	2.83E−09	8.21E−14	3.15E−11
4.653E−01	9.84E−01	1.68E+00	3.54E−05	3.77E−09	4.87E−14	3.15E−11
4.912E−01	9.92E−01	1.68E+00	7.47E−05	4.55E−09	2.38E−14	3.02E−11
5.169E−01	1.00E+00	1.68E+00	1.93E−03	2.69E−10	2.14E−16	4.38E−13

In the aftermath, the number of operations will be analyzed, since it presents another way to compare the computational cost. Considering the hybrid method, the domain over equations (11a) and (11b) is separated into two parts. The part calculated with the analytical solution has some simplifications, which leads to the following governing equations for this subdomain

\frac{\partial u}{\partial τ} = (Fo + Fo \cdot Ko \cdot Lu \cdot Pn) \frac{\partial^{2} u}{\partial χ^{2}}

(27a)

\frac{\partial v}{\partial τ} = Fo \cdot Lu \cdot Pn \frac{\partial^{2} u}{\partial χ^{2}}

(27b)

Thus, this representation substitutes the original problem, given by equations (11a) and (11b), by two new differential linear equations. With these assumptions, there are no more updates of properties for this part in which the coupling effects of moisture over the temperature are neglected. It is important to note that these differential linear equations can be solved using analytic solution, as presented in the hybrid method. As a consequence, solving the two linear differential equations implies a lower computational cost than solving the original ones. At the same time, this is not the only part concerning to the hybrid method, and a deeper analysis must be carried out.

Notations $N_{χ}$ and $N_{τ}$ denote the spatial and time discretization parameter, and $N_{prop}$ denotes the number of material properties to update as a function of the field at each time step iteration. The computational cost of the hybrid method is estimated to be $C^{hyb} = ({N_{χ,}}_{Ω_{1}} \cdot N_{τ}) (N_{prop} \cdot {N_{χ,}}_{Ω_{1}}) + ({N_{χ,}}_{Ω_{2}} \cdot N_{τ}) + (N_{χ} \cdot N_{τ})$ . The first term represents the numerical part as ${N_{χ,}}_{Ω_{1}}$ discrete points are computed for the $N_{τ}$ time iterations with $N_{prop}$ material properties updated for the ${N_{χ,}}_{Ω_{1}}$ discrete points. The second term is due to the analytical part of the hybrid method and the last term due to the pre-simulation. For the finite-difference method, the computational cost is $C^{num} = (N_{χ} \cdot N_{τ}) (N_{prop} \cdot N_{χ})$ . The computational gain with the hybrid method is not so obvious. For this reason, the occurrences when $C^{hyb} \leq C^{num}$ are analyzed

\underset{numerical part}{\underset{︸}{\frac{N_{χ, Ω_{1}}^{2}}{N_{χ}^{2}}}} + \underset{analytical part}{\underset{︸}{\frac{{N_{χ,}}_{Ω_{2}}}{N_{prop} \cdot N_{χ}^{2}}}} + \underset{pre - simulation}{\underset{︸}{\frac{1}{N_{prop} \cdot N_{χ}}}} \leq 1

(28)

The first term corresponds to the numerical part of the hybrid method. Thus, to make the sentence truthful, it is necessary that ${N_{χ,}}_{Ω_{1}} << N_{χ}$ , and to make it true, the mass diffusion process must be very slow, leading to a very small Luikov number, $Lu << 1$ and, consequently, to a short penetration depth, which is likely to occur in soil heat and moisture problems, making an appropriate case for the application of the hybrid method.

The analytical part corresponds to the second term of equation (28). The higher the number of properties to update, the lower this term. In addition, as this part is simulated analytically, the number of nodes does not influence the accuracy of the solution, so that they can be decreased as much as necessary. However, as the analytical solution is running with the numerical part, its advantage is lost due to the way it is implemented, with time steps.

Furthermore, the third term of equation (28) concerns to the pre-simulation cost. Similar to the second term, the higher the number of properties to update, the lower this term. In addition, if the properties vary also with the temperature, the computational cost of the hybrid method will be even lower. In summary, the hybrid method offers an interesting approach when the diffusion process is very slow, and a great number of properties must be updated at each iteration.

Final remarks

Conclusion

This work proposed a hybrid method to compute the heat and moisture transfers in porous soils. It consists of finding an area where the second derivative of the moisture content is approximately null. In this area, a simplification of the problem is considered and calculated with an analytical solution, leading to a reduction of the computational cost.

The method has been illustrated in a 1D nonlinear case of heat and moisture transfer in unsaturated soils. A 20-m domain was considered with a 10-year time simulation length. The analytical solution is implemented deep into the ground, where variations of properties can be neglected. The numerical solution is implemented where the effect of the nonlinearities is more intense. Both solutions are connected through an iterative process which is performed until reaching the convergence.

Since no analytical solution exists for the complete problem, its verification was performed against a numerical model. It has shown an excellent agreement.

Furthermore, the domain separation may reduce calculations, considering the properties to be constant 5 m below the surface and thus ignoring the second derivative of moisture content in this part. The hybrid method offers an interesting reduction of the computational cost when the diffusion process is very slow, and a great number of properties varies with the field. These conditions are reached in most problems of soil heat and moisture transfer and, therefore, suit well to the hybrid method.

Further investigation is needed to improve the hybrid method. Particularly, the determination of the hybrid interface location can be enhanced with a process calculating alongside the hybrid method, updating the position of the interface among the time iterations. Another point is the way both domains are coupled; perhaps, if the domains are overlapped, the model can become more robust. In addition, modeling the problem in two or three dimensions with the hybrid method, the computational cost can be significantly reduced.

Outlooks

The methodology presented in this manuscript considered only unidimensional problems. Given the promising results, it can be applied in the context of whole-building energy simulation, aiming at estimating the building energy consumption, where the diffusion transfer in the spatial domains (wall and soils) is assumed as unidimensional. However, in some practical cases, 2D (or even 3D) transfer might be considered for a more accurate assessment of diffusive transfer phenomena in building foundations.

The hybrid approach can be used in such cases. The two subdomains are set as illustrated in Figure 11. The first domain is defined where the field has high magnitudes of variations. In this domain, a numerical model considering 2D transfer and nonlinear material properties can be used. In the second domain, the variations of the dependent variable are assumed to be negligible. The analytical solution as the ones proposed by Comini and Lewis (1976) and Guimarães et al. (2010) can be used in this second domain.

Figure 11.

Schematic representation of the ground problem considering 2D transfer.

Another possibility that deserves investigation is to consider an analytical solution for unidimensional diffusive transfer for the simplified model. Indeed, the magnitudes of the variations of the field have to be assumed small compared to the ones within the first domain. Moreover, due to adiabatic BCs, the fluxes are unidimensionally directed (over ψ).

An important issue when considering 2D transfer is the interface of the hybrid model. As for the unidimensional approach, the continuity of the fluxes and of the dependent variables fields is ensured at the interface by equation (17). It may be advantageous to have overlapped domains to guarantee this condition. The location of the interface can be determined with a pre-simulation, similar to what has been done for the unidimensional case presented here, or, for practical reasons, it could be determined by other means. The option is implementing the hybrid algorithm in a such way the interface moves along the simulation, supported by the computations of the fluxes to determine its location. As we believe that the results are promising, research for 2D transfer will be conducted in the near future.

Footnotes

Appendix 1

Appendix 2 Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Brazilian Agencies CAPES of the Ministry of Education and CNPQ of the Ministry of Science, Technology and Innovation.

References

Abahri

Belarbi

Trabelsi

(2011) Contribution to analytical and numerical study of combined heat and moisture transfers in porous building materials. Building and Environment 46(7): 1354–1360.

Berger

Chhay

Guernouti

et al . (2015) Proper generalized decomposition for solving coupled heat and moisture transfer. Journal of Building Performance Simulation 8: 295–311.

Bouddour

Auriault

J-L

Mhamdi-Alaoui

et al . (1998) Heat and mass transfer in wet porous media in presence of evaporation-condensation. International Journal of Heat and Mass Transfer 41(15): 2263–2277.

Chang

Weng

(2000) An analytical solution to coupled heat and moisture diffusion transfer in porous materials. International Journal of Heat and Mass Transfer 43(19): 3621–3632.

Cichota

Elias

de Jong van Lier

(2004) Testing a finite-difference model for soil heat transfer by comparing numerical and analytical solutions. Environmental Modelling & Software 19(5): 495–506.

Claesson

Hagentoft

C-E

(1991) Heat loss to the ground from a building—I. General theory. Building and Environment 26(2): 195–208.

Coelho

LDS

Freire

Dos Santos

et al . (2009) Identification of temperature and moisture content fields using a combined neural network and clustering method approach. International Communications in Heat and Mass Transfer 36(4): 304–313.

Comini

Lewis

(1976) A numerical solution of two-dimensional problems involving heat and mass transfer. International Journal of Heat and Mass Transfer 19(12): 1387–1392.

Crank

Nicolson

(1947) A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Mathematical Proceedings of the Cambridge Philosophical Society 43: 50–67.

10.

De Freitas

Abrantes

Crausse

(1996) Moisture migration in building walls—analysis of the interface phenomena. Building and Environment 31(2): 99–108.

11.

Dos Santos

Mendes

(2004) Multidimensional effects of ground heat transfer on the dynamics of building thermal performance. ASHRAE Transactions 110(2): 345–353.

12.

Dos Santos

Mendes

(2006) Simultaneous heat and moisture transfer in soils combined with building simulation. Energy and Buildings 38(4): 303–314.

13.

Godinho

Amado-Mendes

Pereira (2011) A hybrid analytical-numerical model based on the method of fundamental solutions for the analysis of sound scattering by buried shell structures. Mathematical Problems in Engineering 2011: e710623.

14.

Guimarães

Delgado

de Freitas

(2010) Mathematical analysis of the evaporative process of a new technological treatment of rising damp in historic buildings. Building and Environment 45(11): 2414–2420.

15.

Janssen

Carmeliet

Hens

(2002) The influence of soil moisture in the unsaturated zone on the heat loss from buildings via the ground. Journal of Building Physics 25(4): 275–298.

16.

Janssen

Carmeliet

Hens

(2004) The influence of soil moisture transfer on building heat loss via the ground. Building and Environment 39(7): 825–836.

17.

Lee

Kim

Chung

et al . (2010) Analytical solution for heat and moisture diffusion in layered materials. Canadian Geotechnical Journal 47(6): 595–608.

18.

Liu

Peng

(2005) Study of heat and moisture transfer in soil with a dry surface layer. International Journal of Heat and Mass Transfer 48(21–22): 4579–4589.

19.

MATLAB (2013) Version 8.1.0.604 (R2013a). Natick, MA: The MathWorks Inc.

20.

Mendes

Abadie

Akinyemi

(2012) Numerical and experimental assessment of ground temperature and water content under mild weather conditions. Special Topics & Reviews in Porous Media: An International Journal 3(4): 329–339.

21.

Mendes

Barbosa

Freire

et al . (2008) A simulation environment for performance analysis of HVAC systems. Building Simulation 1(2): 129–143.

22.

Mendes

Philippi

Lamberts

(2002) A new mathematical method to solve highly coupled equations of heat and mass transfer in porous media. International Journal of Heat and Mass Transfer 45(3): 509–518.

23.

Mikhailov

Özisik

(1994) Unified Analysis and Solutions of Heat and Mass Diffusion. New York: Dover.

24.

Mulay

Worek

(1990) Simultaneous transport of heat and moisture in a partially saturated porous media. International Communications in Heat and Mass Transfer 17(2): 135–145.

25.

Nowamooz

Nikoosokhan

Lin

et al . (2015) Finite difference modeling of heat distribution in multilayer soils with time-spatial hydrothermal properties. Renewable Energy 76: 7–15.

26.

Oliveira

Jr Freitas

Prata

(1993) Influence of medium properties in diffusivity of Phillip and de Vries model in unsaturated soils. Technical report, University of Santa Catarina, Florianópolis.

27.

Özisik

(1993) Heat Conduction, 2nd edn.New York: A Wiley-Interscience.

28.

Pallin

Kehrer

(2013) Hygrothermal simulations of foundations: part 1: soil material properties. Journal of Building Physics 37(2): 130–152.

29.

Qin

Belarbi

Ait-Moktar

et al . (2006) An analytical method to calculate the coupled heat and moisture transfer in building materials. International Communications in Heat and Mass Transfer 33: 39–48.

30.

Quarteroni

Valli

(1999) Domain Decomposition Methods for Partial Differential Equations. Oxford: Oxford University Press.

31.

Rantala

(2005) Estimation of the mean temperature distribution underneath a slab-on-ground structure. Journal of Building Physics 29(1): 51–68.

32.

Rantala

Leivo

(2009) Heat, air, and moisture control in slab-on-ground structures. Journal of Building Physics 32(4): 335–353.

33.

Rees

Adjali

Zhou

et al . (2000) Ground heat transfer effects on the thermal performance of earth-contact structures. Renewable and Sustainable Energy Reviews 4(3): 213–265.

34.

Rees

Zhou

Thomas

(2001) The influence of soil moisture content variations on heat losses from earth-contact structures: an initial assessment. Building and Environment 36(2): 157–165.

35.

Ribeiro

Cotta

Mikhailov

(1993) Integral transform solution of Luikov’s equations for heat and mass transfer in capillary porous media. International Journal of Heat and Mass Transfer 36(18): 4467–4475.

36.

Shanthikumar

Sargent

(1983) A unifying view of hybrid simulation/analytic models and modeling. Operations Research 31(6): 1030–1052.

37.

Thomas

Morgan

Lewis

(1980) A fully nonlinear analysis of heat and mass transfer problems in porous bodies. International Journal for Numerical Methods in Engineering 15(9): 1381–1393.

38.

Yamazawa

(2001) A one-dimensional dynamical soil–atmosphere tritiated water transport model. Environmental Modelling & Software 16(8): 739–751.