math.math-mp
Docteur De 
  
Suelen Gasparin 
email: suelengasparin@hotmail.com
  
Thèse Dirigée Par 
  
Denys Dutykh 
  
M Jean-Jacques 
  
Roux Professeur 
  
M Aristeu Da Silveira 
  
Neto Professeur 
  
M Carsten 
  
Rode Professeur 
  
Mme Monika 
  
Woloszyn Professeur 
  
M Paulo César 
  
Philippi Professeur 
  
M Denys Dutykh 
  
M Nathan 
  
Mendes Professeur 
  
Dr Nathan Mendes 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
THÈSE Pour obtenir le grade de

Keywords: numerical simulation, heat and moisture transfer in porous materials, nonlinear diffusion problems, reduced-order model, spectral methods, parametric problems, multi-zone building simulation, two-dimensional transfer. vii numerical simulation, building physics, heat and moisture transfer, nonlinear diffusion problems, reduced-order model, spectral methods Programming Latex, Maple, Matlab

Building energy consumption is directly impacted by weather parameters such as temperature, solar radiation, atmospheric pressure, relative humidity and wind velocity. The knowledge of the building hygrothermal performance enables the design of energy efficient buildings and the prediction of overall durability and sustainability of envelopes. Therefore, designers and builders are interested in modeling the long-term performance of the envelopes by means of accurate, reliable and fast simulation tools.

Several numerical models have been proposed in the literature to study the heat and moisture transfer in building materials. In general, this problem is solved by traditional methods, such as finite-difference and finite-volume methods, using mainly implicit schemes. Nevertheless, these methods impose costly sub-iterations to treat the nonlinearities and very fine discretization, which increase substantially the simulation computational cost. Therefore, this research has been focused on the development and analyses of numerical methods for efficiently simulate the problem of heat and mass transfer through porous materials.

In the first part of this thesis, improved schemes of the traditional numerical methods have been developed to eliminate costly sub-iterations to treat nonlinearities, to improve the order of accuracy and to save computer run time. Despite the great progress with the new numerical schemes, the conclusion of the first part shows that we still have to deal with large systems of equations, particularly when treating multi-dimensional transfer problems. For this reason, to reduce even more the computational burden and the size of the system, a reduced-order model, based on spectral methods is proposed in the sequence to provide an accurate description of the physical phenomena. The degrees of freedom of the solution is strongly decreased while maintaining the model fidelity. It ensures a computational cost much lower than the complete original model.

All these methods are applied to problems related to building physics, such as single and multilayer nonlinear transfer, the determination of optimum insulation thickness, the process of moisture buffer effects and transfer in one-or two-zone building models. In conclusion, we show how to build efficient numerical models, in terms of computational cost and accuracy, to investigate the heat and mass transfer in porous materials.

Resumé

La consommation énergétique des bâtiments est directement impactée par les variables climatiques et en particulier par la température, la radiation solaire, la pression atmosphérique, l'humidité relative et la vitesse du vent. La maîtrise de la performance hygrothermique des bâtiments permet la conception d'ouvrage à faible consommation énergétique intégrant des enveloppes durables et dépourvues de risques de pathologies liées à l'humidité. Ainsi, les architectes et constructeurs nécessitent de disposer d'outils, pour simuler les performances de l'enveloppe sur le long terme, qui reposent sur des modèles précis, fiables et rapides.

Plusieurs modèles ont été proposés dans la littérature pour investiguer les transferts de chaleur et de masse dans les matériaux poreux composant les bâtiments. Généralement, les équations régissant les phénomènes sont résolues en employant des méthodes conventionnelles basées sur les différences finies ou les volumes finis combinés avec des schémas temporels implicites. Cependant, ces méthodes imposent des sous-itérations, coûteuses en temps de calculs, pour traiter les non linéarités du problème. De plus, elles exigent des discrétisations très fines qui augmentent significativement les temps de calcul des simulations. Par conséquent, ces travaux de recherches se concentrent sur le développement et l'analyse de méthodes numériques pour simuler de manière efficiente les problèmes de transferts de chaleur et de masse dans les matériaux poreux.

Dans la première partie de cette thèse, des schémas améliorés par rapport aux approches traditionnelles sont proposés pour éliminer les sous-itérations coûteuses pour traiter les non linéarités, pour améliorer l'ordre de précision et pour diminuer les temps de calculs. Nonobstant les progrès réalisés avec les nouveaux schémas numériques, les conclusions de cette première partie montrent que nous manipulons toujours des systèmes d'équations de grandes dimensions, en particulier lorsqu'il s'agit de problèmes de transferts multidimensionnels. Pour cette raison, afin de diminuer d'avantage les coûts de calcul et la taille du modèle, un modèle d'ordre réduit basé sur les méthodes spectrales est proposé subséquemment. Il garantit une prédiction précise des phénomènes physiques. De plus, le degré de liberté de la solution est fortement réduit tout en préservant la fidélité du modèle. Le coût de calcul est nettement inférieur à celui du modèle complet originel.

L'ensemble de ces méthodes sont appliquées à des problèmes relatifs aux phénomènes physiques intervenant dans les bâtiments, tels que les transferts mono ou multi-couches, la recherche de l'épaisseur optimale d'isolation, les mécanismes de tampons hygroscopiques ou encore les transferts à l'échelle de bâtiments mono-ou bi-zones. En conclusion, nous ix démontrons comment construire des modèles numériques efficaces, en termes de temps de calculs et précision, pour l'investigation des transferts de chaleur et de masse dans les matériaux poreux.

Mots clés: Simulation numérique, transferts de chaleur et de masse dans les matériaux poreux, problèmes de diffusion non-linéaires, modèles d'ordres réduits, méthodes spectrales, problèmes paramétriques, simulation de bâtiments multizones, transferts bidimensionnels x

Resumo

O consumo de energia em edificações é diretamente influenciado por parâmetros climáticos como temperatura, radiação solar, pressão atmosférica, umidade relativa e velocidade do vento. O conhecimento do desempenho higrotérmico das edificações permite a concepção de edifícios energeticamente eficientes e a previsão da sua durabilidade global. Portanto, designers e construtores estão interessados em modelar o desempenho higrotérmico a longo prazo dos edifícios por meio de ferramentas de simulação precisas, confiáveis e rápidas.

Diversos modelos numéricos fora propostos na literatura para estudar a transferência de calor e umidade em materiais de construção. Em geral, esse problema é resolvido por métodos tradicionais, como os métodos de diferenças finitas e volumes finitos, usando principalmente esquemas implícitos. No entanto, esses métodos impõem custosas subiterações para tratar as não-linearidades e discretizações muito refinadas, que aumentam substancialmente o custo computacional da simulação. Portanto, esta pesquisa foi focada no desenvolvimento e análise de métodos numéricos para simular eficientemente o problema de transferência de calor e massa através de materiais porosos.

Na primeira parte desta tese, esquemas aprimorados dos métodos numéricos tradicionais foram desenvolvidos para eliminar as sub-iterações causadas pelas não-linearidades, para melhorar a ordem de precisão e para economizar tempo de simulação do computador. Apesar do grande progresso com os novos esquemas numéricos, a conclusão da primeira parte mostra que ainda temos que lidar com grandes sistemas de equações, principalmente ao tratar problemas de transferência multidimensional. Por esta razão, para reduzir ainda mais o custo computacional e o tamanho do sistema de equações, um modelo de ordem reduzida, baseado em métodos espectrais, é proposto para fornecer uma descrição precisa dos fenômenos físicos. Garantindo assim, um custo computacional muito inferior ao modelo original completo.

Todos os métodos propostos são aplicados a problemas relacionados à física do edifício, tais como a transferência não-linear através de uma e de múltiplas camadas de material, a determinação da espessura ótima do isolamento, a retenção e liberação da umidade na parede e a transferência em construção de uma e duas zonas. Para concluir, mostra-se como construir modelos numéricos eficientes, em termos de custo computacional e precisão, para investigar a transferência de calor e massa em materiais porosos.

Palavras-chaves: simulação numérica, transferência de calor e umidade, materiais porosos, problemas de difusão não-linear, modelos de ordem reduzida, métodos espectrais, xi simulação multi-zonal em edifícios, transferência bidimensional, problemas paramétricos. xiv

Introduction

The energy efficiency in buildings has been a large theme of discussion in the past five decades since the energy crises in the seventies. According to the International Energy Agency (2019), there is a forecast of growth in the global energy demand of 30%, from 2017 until 2040. Therefore, it is important to manage the energy consumption efficiently by causing less impact to the planet.

To ensure reasonable and sustainable use of energy in buildings, it is necessary to understand the physical phenomena involved. Buildings are constantly exposed to variations of temperature, moisture and air pressure which impact directly on their energy consumption. Designers and builders are interested in knowing the long-term performance of building envelopes, using simulation tools. Simulating the hygrothermal behaviour allows to make a good design and to determine the overall durability and sustainability of the building envelope [START_REF] Tariku | Transient model for coupled heat, air and moisture transfer through multilayered porous media[END_REF][START_REF] Woloszyn | Tools for performance simulation of heat, air and moisture conditions of whole buildings[END_REF].

Initially, the major hypothesis was that energy consumption was related only to the heat transfer. Although, in porous materials, moisture transfer and its accumulation have a direct impact on the heat transfer [START_REF] Bouddour | Heat and mass transfer in wet porous media in presence of evaporation-condensation[END_REF][START_REF] Mulay | Simultaneous transport of heat and moisture in a partially saturated porous media[END_REF][START_REF] Rees | The influence of soil moisture content variations on heat losses from earth-contact structures: an initial assessment[END_REF]. As presented by [START_REF] Deru | A model for ground-coupled heat and moisture transfer from buildings[END_REF], to precisely determine the heat losses, simultaneous calculations with moisture content are required. Temperature and moisture contents are highly interconnected phenomena and, hence, must be simulated together (Berger et al., 2015b;[START_REF] Mendes | Models for prediction of heat and moisture transfer through porous building elements[END_REF].

Fundamental studies of heat and moisture transfer have been published since 1950's by [START_REF] Philip | Moisture movement in porous materials under temperature gradients[END_REF] and [START_REF] Luikov | Heat and mass transfer in capillary-porous bodies[END_REF]. They represented the coupled process of heat and moisture transfer by a system of two nonlinear partial differential equations which uses the temperature and the moisture content gradients as driving potentials. However, detailed simulation of the coupled heat and moisture transfer has only been performed with the improvement of the computer systems. Indeed, the computational cost to solve the coupled problem is at least the double than the computational cost to solve only the heat conduction equation [START_REF] Deru | A model for ground-coupled heat and moisture transfer from buildings[END_REF].

Computational software such as Domus [START_REF] Mendes | A simulation environment for performance analysis of hvac systems[END_REF], EnergyPlus [START_REF] Crawley | Energyplus: Energy simulation program[END_REF], DELPHIN (Bauklimatik Dresden, 2011), MOIST [START_REF] Burch | An analysis of moisture accumulation in walls subjected to hot and humid climates[END_REF], WUFI (Fraunhofer IBP, 2005), TRNSYS (A. [START_REF] Klein | TRNSYS, a transient simulation program[END_REF]) make use of hygrotermal simulations to estimate the energy efficiency of a building [START_REF] Woloszyn | Tools for performance simulation of heat, air and moisture conditions of whole buildings[END_REF]. The model in some of those tools considers the energy transfer governed by diffusion and latent mechanisms while the moisture transfer occurs due to capillary migration and vapour diffusion. The mathematical problem can be written by a system of two partial differential equations [START_REF] Luikov | Heat and mass transfer in capillary-porous bodies[END_REF]: (1.1b) where H [J/m 3 ] is the enthalpy, t [s] is the time, j q [W/m 2 ] is the total heat flux, w [kg/m 3 ] is the moisture content, j m [kg/(s • m 2 )] is the moisture flow and S [kg/(m 3 • s)] is a volumetric source term. Different potentials can be used to develop these Equations (1.1), such as temperature T , vapour pressure P v , capillary pressure P c , volumetric moisture content w, etc. To solve this System of Equations (1.1), some important challenges have to be faced:

∂H ∂t = -∇ • j q , (1.1a) ∂w ∂t = -∇ • j m + S ,
• The storage and diffusivity coefficients are highly nonlinear, varying relatively to both of water moisture and temperature fields, and sometimes also varying according to the space, in the case of multi-layered materials.

• The boundary conditions must be of Robin-type, which are function of time according to weather data. In addition, phenomena such as solar radiation, wind and rain vary also as function of wall orientation.

• There is a difference on the time scales in building physics: the horizon of simulation is of the order of the years, while the characteristic time of diffusion is the order of hours and minutes.

• The integration of complex geometries as corners, foundation, roofs and walls in the whole building simulation requires two-and three-dimensional simulations which increases considerably the computational difficulty.

Due to the complexities of the problem raised above, the use of analytical solutions is impossible. Therefore, the approach used by the software to solve these equations is based on numerical methods, which make use of discrete representations of the continuous equations. The basics of a numerical approach is the idea of making an approximation of the solution, which takes a solution for a finite number of Degrees Of Freedom (DOF). The greater the number of discrete points, the closer to the exact solution will be the approximated solution if the numerical method converges to the true solution [START_REF] Maliska | Heat Transfer and Computational Fluid Mechanics[END_REF].

Problem statement

In literature, the heat and moisture transfer are generally solved by traditional approaches, such as the finite-difference method [START_REF] Fan | An improved model of heat and moisture transfer with phase change and mobile condensates in fibrous insulation and comparison with experimental results[END_REF][START_REF] Qin | Coupled heat and moisture transfer in multi-layer building materials[END_REF], the finite-volume method [START_REF] Manz | Experimental and numerical study of a mechanically ventilated glass double façade with integrated shading device[END_REF][START_REF] Mendes | A method for predicting heat and moisture transfer through multilayered walls based on temperature and moisture content gradients[END_REF][START_REF] Santos | Combined heat, air and moisture (HAM) transfer model for porous building materials[END_REF] and the finite-element method [START_REF] Lu | Modelling of heat and moisture transfer in buildings: I. model program[END_REF][START_REF] Rouchier | Influence of concrete fracture on the rain infiltration and thermal performance of building facades[END_REF]Thomas et al., 1980). The discretization of the problem regarding the time normally uses implicit schemes, such as the Crank-Nicolson and the Euler implicit scheme [START_REF] Janssen | Simulation efficiency and accuracy of different moisture transfer potentials[END_REF][START_REF] Mendes | A method for predicting heat and moisture transfer through multilayered walls based on temperature and moisture content gradients[END_REF]. They are chosen due to the supposed unconditional stability property, which allows the use of larger time steps. Unlike explicit schemes that are limited by the CFL stability condition, which requires very refined time step. Nevertheless, when treating the nonlinearities of the problem, implicit schemes require important sub-iterations to reach the convergence of the solution, which makes the simulation highly time-consuming [START_REF] Abuku | Impact of wind-driven rain on historic brick wall buildings in a moderately cold and humid climate: Numerical analyses of mould growth risk, indoor climate and energy consumption[END_REF]. In addition, Dalgliesh et al. (2005); Dos [START_REF] Santos | Simultaneous heat and moisture transfer in soils combined with building simulation[END_REF] reported the high computational cost of simulation tools to compute the heat and moisture transfer phenomena on the scale of buildings. When the solution of one problem requires long simulation horizons (years), considering an entire building, with a mesh grid refinement, the computation becomes too time-consuming to be useful in practice. Therefore, how to simulate efficiently, i.e., accurately and rapidly, nonlinear problems of heat and moisture transfer in buildings?

First attempts of alternative approaches

Model reduction techniques can be used as an alternative approach to solve the coupled heat and moisture transfer problem. The intent is to construct reduced-order models (ROMs) to provide an accurate description of the physical phenomena by decreasing the number of degrees of freedom, while retaining the model fidelity, at a computational cost much lower than the large original model [START_REF] Reddy | Constrained reduced-order models based on proper orthogonal decomposition[END_REF]. In recent years, reduced-order modelling techniques have proven to be powerful tools for solving various problems. Important efforts have been dedicated to developing reduced-order models that can provide accurate predictions while dramatically reducing computational time, for a wide range of applications, covering different fields, including fluid mechanics, heat transfer, structural dynamics, among others [START_REF] Bai | Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems[END_REF][START_REF] Herzet | Model reduction from partial observations[END_REF][START_REF] Lucia | Reduced-order modeling: new approaches for computational physics[END_REF]. Reduced-order models, such as Proper Orthogonal Decomposition (POD) [START_REF] Tallet | Optimal flow control using a POD-based reduced-order model[END_REF], Modal Basis Reduction (MBR) [START_REF] Gao | Dynamical building simulation: A low order model for thermal bridges losses[END_REF][START_REF] Kim | Decomposition and coupling of soil domain for modeling vertical ground heat exchangers using the state model size reduction technique[END_REF] and Proper Generalized Decomposition (PGD) [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF], have shown a relevant reduction of the computational cost and have been successfully employed by the building physics community (Berger et al., 2015a). In those works, they have applied reducedorder models to build accurate solutions with less computational effort than the complete original model. Reduced-order models can be classified as a priori or as a posteriori. The a posteriori based approaches need a preliminary computed (or even experimental) solution data of the large original problem to build the reduced one. Whereas the a priori ones do not need preliminary information on the studied problem. The reduced-order model is unknown a priori and is built directly. A careful attention must be paid regarding the definition of ROMs since sometimes it is related to the degradation of the physical model [START_REF] Schilders | Introduction to Model Order Reduction[END_REF], which is not the case of the present work.

Another promising approach to solve the coupled transfer problem is the spectral method, which is a robust and highly accurate method that has been applied to solve partial differential equations since the 70's, but lost its spot due to the difficulties to treat nonlinearities, complex geometries, irregular domains, and non-periodic boundary conditions. However, Spectral methods have overcome some of the mentioned difficulties and now they are successfully applied in studies of wave propagation, meteorology, computational fluid dynamics, quantum mechanics and other fields [START_REF] Canuto | Spectral Methods -Fundamentals in Single Domains[END_REF]. The difference between the Spectral methods and the other reduced-order methods is the exponential convergence and the low dissipation and dispersion errors, making their use very attractive.

Some works related to transport phenomena can be found in literature involving diffusive [START_REF] Guo | The Application of the Chebyshev-Spectral Method in Transport Phenomena[END_REF][START_REF] Wang | Collocation methods for uncertain heat convection-diffusion problem with interval input parameters[END_REF], convective [START_REF] Chen | Spectral collocation method for natural convection in a square porous cavity with local thermal equilibrium and non-equilibrium models[END_REF][START_REF] Ramreddy | A spectral relaxation method for linear and non-linear stratification effects on mixed convection in a porous medium[END_REF] and radiative [START_REF] Chen | Chebyshev collocation spectral method for solving radiative transfer with the modified discrete ordinates formulations[END_REF][START_REF] Li | Iterative and direct Chebyshev collocation spectral methods for one-dimensional radiative heat transfer[END_REF][START_REF] Ma | Thermal radiation heat transfer in one-and two-dimensional enclosures using the spectral collocation method with full spectrum k-distribution model[END_REF] heat transfers. The spectral techniques applied in these works are diversified, adopted according to the geometry, boundary conditions and field of application. In recent works, researchers have implemented spectral methods for solving heat and moisture transfer in food engineering [START_REF] Pasban | Spectral method for simulating 3D heat and mass transfer during drying of apple slices[END_REF] and on fluid flow [START_REF] Motsa | On the new bivariate local linearisation method for solving coupled partial differential equations in some applications of unsteady fluid flows with heat and mass transfer[END_REF][START_REF] Nascimento | A comparison of fourier pseudospectral method and finite volume method used to solve the burgers equation[END_REF]. It appears as an interesting alternative to reduce the computational time and increase accuracy that has never been applied to the field of building physics.

Scope

Therefore, the scope of this thesis consists in the development of several efficient computational models, based on original mathematical methods, to be applied to the context of building simulations. The main objective is to extend improved standard techniques and reduced order models to uni-and bi-dimensional nonlinear heat and moisture trans-fer problems, in order to reduce computational cost of simulations without the loss of accuracy. This thesis analyses several numerical methods for simulating transfer problems in building physics. Some of these works have already been published along the doctoral period. For this reason, this document presents an overview of the work, and the extension of the problems here presented that can be found in the articles in the Appendix C. The main purpose is to present our contribution in a more compact document. The following list presents the published articles already published:

1. Gasparin et al. (2018) This document is followed by Chapter 2 that describes each numerical method. First, we take the one-dimensional nonlinear diffusion equation to present the DuFort-Frankel method, the Method of Horizontal Lines, the Quasi-Uniform Nonlinear Transformation and the Spectral method. Then, the spectral method is extended to solve other types of problem: one-dimensional parametric simulation and two-dimensional transient transport.

Chapter 3 applies the methodology of Chapter 2 in real case studies. The objective is to consider simple (but realistic) cases to illustrate the strength of the proposed numerical methods. The first test case considers one-dimensional moisture transfer through a porous material. Each solution method is separately analysed in terms of accuracy and computer run time and at the end, they are all compared against each other. The case study is elaborated to present the features of each method. More complex numerical experiments can be found in the enclosed articles.

Chapter 4 addresses the concluding remarks and the future perspectives. In this chapter, different innovative numerical methods are described, focused on the diffusion transport phenomena. They are divided into two groups: Group A, with the Improved Complete Original Models and Group B, with the Reduced-Order Models. Group A has two sub-groups, one with a fixed spatial grid, composed by the DuFort-Frankel method and another sub-group with the adaptive grid, which comprises both the Quasi-Uniform Nonlinear Transformation and the Method of Horizontal Lines. In the Group B, the Spectral method is the reduced-order model that has been chosen. The Spectral techniques here addressed are the Tau-Galerkin method, the Collocation method and the Chebyshev polynomials. The Spectral method is extended to two-dimensional applications and to parametric problems. Figure 2.1 shows graphically the classification of the numerical methods to be presented in this thesis manuscript. The first section of this chapter contains the definitions related to the general notation, to the mathematical model and to the discretization parameters. Then, each section contains one numerical method: Section 2.2 describes the DuFort-Frankel approach; Section 2.3 the Method of Horizontal Lines; Section 2.4 the method of Quasi-Uniform Nonlinear Transformation; and Section 2.5 presents the Spectral methods. After the description of each method, their strengths and drawbacks are highlighted in the section Important Features and how they were applied in the context of building physics specified in the section Extended Application.

Definitions

Notation

Consider x, x and t as being the independent variables, which domains are defined as:

x ∈ I x ≡ [ 0, 1 ] , x ∈ I x ≡ [ -1, 1 ] , and t ∈ I t ≡ [ 0, τ ] .
The field of interest u represents the unknown function, which can be defined as:

u : I x × I t -→ IR , (x, t) -→ u (x, t) ,
for one-dimensional problems. Some partial derivatives are written in a simplified way to have a short-hand notation:

u t def := ∂u ∂t , u x def := ∂u ∂x , and u xx def := ∂ 2 u ∂x 2 .

The mathematical problem

For explaining the numerical methods, consider the following nonlinear parabolic equation that can represent the dimensionless heat or the moisture transfer inside a porous material:

c • ∂u ∂t = Fo • ∂ ∂x k • ∂u ∂x , ( 2.1) 
where Fo = const > 0 is the Fourier number that characterizes the ratio between the reference time and time of the diffusion phenomenon, c is a dimensionless storage coefficient and k is the dimensionless diffusivity transfer coefficient. The transport coefficients c and k are assumed to be real positive functions:

c : IR -→ IR > 0 , k : IR -→ IR > 0 , u -→ c (u) , u -→ k (u) ,
with IR > 0 = { y ∈ IR | y > 0 } being the set of all positive real numbers. In addition, let us assume that the function k is continuously differentiable, so its derivative is given by:

k (u) def := ∂k (u) ∂u .
The linear form of Equation (2.1) can be written as the Fourier equation [START_REF] Evans | Partial Differential Equations[END_REF][START_REF] Fourier | Théorie analytique de la chaleur[END_REF]:

∂u ∂t = α • ∂ 2 u ∂x 2 , (2.2)
where α = (Fo • k)/c is a constant diffusivity coefficient. Equations (2.1) and (2.2) have to be supplemented with compatible initial and boundary conditions to obtain a well-posed problem [START_REF] Hadamard | Sur les problèmes aux derivées partielles et leur signification physique[END_REF]. Thus, boundary conditions of Robin-type are considered on both sides of the spatial domain:

k • ∂u ∂x = Bi L • u -u L ∞ , for x = 0 and t > 0 , (2.3) -k • ∂u ∂x = Bi R • u -u R ∞ , for x = 1 and t > 0 .
(2.4)

The dimensionless Biot numbers Bi L and Bi R represent the ratio of transfer resistances within the material and the one at its surface. They are considered as constant in time:

Bi L = const > 0 and Bi R = const > 0 .
Furthermore, the ambient fields u L ∞ and u R ∞ are time dependent:

u L ∞ : t -→ u L ∞ (t) , u R ∞ : t -→ u R ∞ (t)
where the superscript L represents the left side (at x = 0) and the superscript R represents the right side (at x = 1) . At the initial state, the field is given by the initial condition:

u (x, t = 0) = u 0 ,
where u 0 is a given function of space:

u 0 : x -→ u 0 (x) .
This is the dimensionless formulation of the physical problem. It is normally presented using the superscript star but it has been removed for the sake of clarity to present the numerical methods.

Discretization

Complete discretization

The continuous spatial domain I x is uniformly discretized, with x j being the spatial nodes:

x j ≡ j • ∆x, j = 0, 1, 2, . . . , N -1, N Nx elements , ∆x = const > 0 ,
where N x = N + 1 is the total number of spatial nodes and ∆x is the spatial step size:

∆x def := N x ,
that can also be related with two spatial nodes of the grid: ∆x = x j+1 -x j . In this case, ≡ 1 since the spatial domain is dimensionless. The temporal nodes t n are defined by a uniform discretization in the time domain I t so that:

t n ≡ n • ∆t , n = 0, 1, 2, . . . N t -1 Nt elements , ∆t = const > 0 ,
where ∆t is the time step that can be defined by:

∆t def := τ N t ,
with N t being the total number of temporal nodes. Thus, the values of the function u (x, t) in discrete nodes will be denoted by:

u n j def := u (x j , t n ) ,
also known as a point value on the grid. The Cartesian representation of the discretization is given in Figure 2.2.

Semi-discretization

To semi-discretize the solution u in space, the surface u (x, t) is replaced by a sequence of vertical lines (the directions are chosen as in Figure 2.3(a)): where ∆x > 0 is the grid spacing and u j is the evolution in time solution. Moreover, to semi-discretize the solution u in time, the surface u (x, t) is replaced by a sequence of horizontal lines:

u j (t) def := u (x j , t) , x j = j • ∆x , j ∈ N 0 ,
u n (x) def := u (x, t n ) , t n = n • ∆t , n ∈ N 0 ,
where ∆t > 0 is a chosen time step. 

Coefficients of the transport equation

In addition, the material properties at the discrete points can be projected on grid nodes as well:

c n j def := c (u n j ) , k n j def := k (u n j ) , k n j+ 1 2 def := k u n j+ 1 2 , k n j-1 2 def := k u n j-1 2 .
The field u evaluated at intermediate nodes j + 1 2 and j -1 2 are defined as:

u n j+ 1 2 = u n j + u n j+1 2 , u n j-1 2 = u n j + u n j-1 2 .
They represent the mean value of the field at the face between two cells, as depicted in Figure 2.4. 

The DuFort-Frankel scheme

Consider the discretization of Equation (2.1) by means of the finite-difference method, with Forward approximation in Time and Central approximation in Space (FTCS), also known as the explicit scheme:

c n j • u n+1 j -u n j ∆t = Fo ∆x 2 • k n j+ 1 2 • u n j+1 -k n j+ 1 2 + k n j-1 2 • u n j + k n j-1 2 • u n j-1 .
(2.5) By using the so-called DuFort-Frankel method [START_REF] Dufort | Stability conditions in the numerical treatment of parabolic differential equations[END_REF], the discrete point u n j from the explicit scheme is replaced by its mean value on two successive time layers:

u n j = u n-1 j + u n+1 j 2 .
Thus, Equation (2.5) is modified to obtain the DuFort-Frankel scheme:

u n+1 j = λ 1 λ 0 + λ 3 • u n j+1 + λ 2 λ 0 + λ 3 • u n j-1 + λ 0 -λ 3 λ 0 + λ 3 • u n-1 j , n 1 , (2.6) with λ 0 def := 1 , λ 1 def := 2 • ∆t • Fo ∆x 2 • c n j • k n j+ 1 2 , λ 2 def := 2 • ∆t • Fo ∆x 2 • c n j • k n j-1 2 , λ 3 def := ∆t • Fo ∆x 2 • c n j • k n j+ 1 2 + k n j-1 2 ,
which is also an explicit scheme but with a time discretization on three levels instead of two. The scheme (2.6) has the stencil depicted in Figure 2.5. When dealing with the nonlinearities of the material properties, the greatest advantage of explicit schemes is that they do not require sub-iterations. At the time layer n + 1, the material properties

k n j+ 1 2 , k n j-1 2
are explicitly calculated at t n . Moreover, the property c n j is also explicitly computed, being know as a frozen coefficient.

Stability analysis.

The standard von Neumann stability analysis on the linear equation shows that the DuFort-Frankel scheme is unconditionally stable [START_REF] Richtmyer | Difference Methods for Initial-Value Problems[END_REF][START_REF] Taylor | The stability of the Du Fort-Frankel method for the diffusion equation with boundary conditions involving space derivatives[END_REF]. The consistency error analysis of the scheme (2.6) for the linear Equation (2.2) shows the following result:

L n j = α • ∆t 2 ∆x 2 τ • ∂ 2 u ∂t 2 + ∂u ∂t -α • ∂ 2 u ∂x 2 + 1 6 • ∆t 2 • ∂ 3 u ∂t 3 - 1 12 • α • ∆x 2 • ∂ 4 u ∂x 4 - 1 12 • α • ∆t 2 • ∆x • ∂ 5 u ∂x 3 ∂t 2 + O ∆t 4 ∆x 2 , (2.7)
where:

L n j def := u n+1 j -u n-1 j 2 • ∆t -α • u n j-1 -u n-1 j + u n+1 j + u n j+1 ∆x 2 .
So, from the asymptotic expansion for L n j one can see the DuFort-Frankel scheme is second-order accurate in time and:

• first-order accurate in space if ∆t ∝ ∆x 3/2 ;

• second-order accurate in space if ∆t ∝ ∆x 2 .

It is important to note that the spatial derivative of the boundary condition has to be discretized to the second order of accuracy O(∆x 2 ) to maintain the stability of the method [START_REF] Taylor | The stability of the Du Fort-Frankel method for the diffusion equation with boundary conditions involving space derivatives[END_REF].

Even though the scheme is unconditionally stable, errors can arise. The problem is that the scheme is inconsistent since the discretized equation leads to a hyperbolic equation which is not the original PDE:

τ • ∂ 2 u ∂t 2 + ∂u ∂t = α • ∂ 2 u ∂x 2 , (2.8)
where τ is defined in Equation (2.7):

τ def := α • ∆t 2 ∆x 2 .
However, for the DuFort-Frankel scheme to be consistent with the original equation, we must have: ∆t ∆x → 0 , so that the therm in the second derivative in time τ • ∂ 2 u ∂t 2 should be as small as possible to not disturb the solution.

The stability of the scheme (2.8) was studied by [START_REF] Taylor | The stability of the Du Fort-Frankel method for the diffusion equation with boundary conditions involving space derivatives[END_REF] and by [START_REF] Chetverushkin | Explicit schemes and numerical simulation using ultrahigh-performance computer systems[END_REF], in which the following relation is obtained so the scheme can be consistent: ∆t ∆x τ α .

(2.9) Therefore, τ is an important parameter in the scheme, particularly since it depends on the discretization parameters.

Boundary conditions. According to [START_REF] Taylor | The stability of the Du Fort-Frankel method for the diffusion equation with boundary conditions involving space derivatives[END_REF], when using the DuFort-Frankel scheme, the spatial derivatives of the boundary conditions should be discretized considering a truncation error of the second order O(∆x 2 ), to preserve the accuracy and stability properties. In this way, Equation (2.3) gives:

k n 0 • -3 • u n+1 0 + 4 • u n+1 1 -u n+1 2 2 • ∆x + O(∆x 2 ) = Bi L • u n+1 0 -u L ∞ t n+1 .
(2.10) Thus, u n+1 0 can be deduced by rearranging Equation (2.10):

u n+1 0 = 4 • u n+1 1 -u n+1 2 + 2 • b • u L ∞ t n+1 3 + 2 • b , (2.11)
where b = ∆x • Bi L /k n 0 . A similar approach is adopted for the last node j = N at the other boundary condition:

k n N • 3 • u n+1 N -4 • u n+1 N-1 + u n+1 N-2 2 • ∆x = Bi R • u n+1 N -u R ∞ t n+1 .
(2.12) which gives:

u n+1 N = 4 • u n+1 N-1 -u n+1 N-2 + 2 • a • u R ∞ t n+1 3 + 2 • a , (2.13) where a = ∆x • Bi R k n N .
Initial condition. To initialize simulations two initial conditions are required since this method has a three level time discretization. Thus, the initial condition u 0 is used twice to start the simulation. If the simulation starts at iteration n = 1 then u 0 is used for n = 1 and for n = 0 : u x, t 0 = u 0 and u x, t 1 = u 0 .

Important Features

The DuFort-Frankel scheme has the benefits of being unconditionally stable, second-order accurate in time O(∆t 2 ) and first-or second-order accurate in space, depending on the choice of ∆t. In addition, it has the advantage to explicitly computing the solution at each time step, avoiding costly sub-iterations. All these characteristics make the method to be an interesting option in co-simulation and in parallel simulations [START_REF] Gasparin | An improved explicit scheme for whole-building hygrothermal simulation[END_REF].

Attention should be paid for every unconditionally stable scheme because the choice of the time discretisation ∆t, which is an important issue to represent accurately the physical phenomena. As mentioned by [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF], the unconditional stability property does not imply that a physically realistic solution will be produced, no matter how large the time step is. The stability does not guarantee physically plausible solutions. Despite of being unconditionally stable, it does not mean that any value of ∆t can be used; it must respect the characteristic time of the problem (Gasparin et al., 2018b). For the coupled heat and moisture transfer we have used the same value of ∆t, according to the characteristic time of the slower diffusion process.

Extended Application

The paper Stable explicit schemes for simulation of nonlinear moisture transfer in porous materials (Gasparin et al., 2018b) is devoted to explore the use of improved explicit schemes to overcome the stability limitation of the standard explicit scheme. The proposed schemes are evaluated to solve nonlinear moisture transfer through porous materials. The first case study considers a linear diffusion transfer through a porous material. The DuFort-Frankel scheme was compared to (i) the Crank-Nicolson, (ii) the hyperbolisation, (iii) the classical Euler explicit scheme and (iv) the reference solution. For ∆t 10 -3 , the error of the DuFort-Frankel scheme is first-order accurate in ∆t, and for ∆t 10 -3 second-order accurate. The second case study focuses on non-linear transfer, with material properties strongly dependent on the vapour pressure field. Both DuFort-Frankel and Crank-Nicolson schemes were used to compute the solution of the problem. Results have shown that the error is proportional to O(∆t). The Crank-Nicolson is twice longer than the DuFort-Frankel to compute the solution, due to the supplementary operations required to compute the implicit and explicit parts of the scheme.

Since the DuFort-Frankel method showed promising results, we extended its application to whole-building hygrothermal simulations, to the scale of the wall, and, the scale of one-and two-zone models. This approach may reduce the computational cost by a factor of twenty, as well as it enables perfect synchronism for whole-building simulation and co-simulation. These results are shown in An improved explicit scheme for whole-building hygrothermal simulation [START_REF] Gasparin | An improved explicit scheme for whole-building hygrothermal simulation[END_REF]. For the first case study, of heat and moisture transfer through a wall, the Euler implicit scheme required around 3 sub-iterations, making it three times more costly than the DuFort-Frankel scheme. When coupling the wall and zone models using implicit schemes, a nonlinear system of equations has to be solved. By using the DuFort-Frankel explicit scheme, the system of equations becomes linear and no sub-iterations are needed. Therefore, with the DuFort-Frankel approach, the algorithm requires only 9% of the CPU time of the implicit approach. When considering the nonlinearities of the wall material properties and long-wave radiative heat transfer among the room surfaces, the computational savings rise to 95%. These results encourage to apply the DuFort-Frankel approach in building simulation tools. The computational gains should further increase with the number of rooms, walls, partitions, and furniture.

The Method of Horizontal Lines

The Method of Horizontal Lines is one of the methods that is adaptive in space. In building physics this type of methods is not common. What is normally used is the non-uniform grids as in [START_REF] Janssen | The influence of soil moisture in the unsaturated zone on the heat loss from buildings via the ground[END_REF] but it does not change over time.

Nowadays, the most popular approach to solve Equation (2.2) numerically is the socalled Method Of Lines (MOL) [START_REF] Hamdi | Method of lines[END_REF][START_REF] Kreiss | Method of lines for hyperbolic differential equations[END_REF][START_REF] Reddy | Stability of the method of lines[END_REF][START_REF] Schiesser | Method of lines solution of the Korteweg-de Vries equation[END_REF][START_REF] Shampine | ODE solvers and the method of lines[END_REF]. By using this technique, the first step is to semi-discretize the diffusion equation in space. Thus, the classical secondorder central finite-difference scheme yields to the following system of coupled ordinary differential equations [START_REF] Richtmyer | Difference Methods for Initial-Value Problems[END_REF]:

du j dt = α ∆x 2 • u j-1 -2 • u j + u j+1 , j = 1, 2, . . . , N -1 ,
remembering that u j = u (x j , t) is the definition for the semi-discretization in space.

For simplicity, the discretization is given only for the transport equation, there are still three remaining equations to complement the problem, two for the boundary conditions and one for the initial condition. Once the PDE (2.2) is discretized with the MOL, it yields to a large system of equations to solve. In typical accurate numerical simulations it has a size of approximatively N ≈ 10 2 . . . 10 3 or even larger. Then, the semi-discretized system can be solved using various explicit or implicit time-marching schemes [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF]. This approach is refereed to as the Method of Vertical Lines (MOVL) which is shown schematically in Figure 2.6(a). In the field of heat and mass transfer the dominant majority of numerical simulations follow this philosophy. However, if one looks at Equation (2.2) from a different perspective, for example, as the Cauchy-Kovalevskaya form, it can be written as [START_REF] Evans | Partial Differential Equations[END_REF]:

α • ∂ 2 u ∂x 2 = ∂u ∂t , (2.14)
which can be solved with a totally different numerical technique. By semi-discretizing Equation (2.14) in time, it yields to the following system:

α • ∂ 2 u (x, t n+1 ) ∂x 2 = u (x, t n+1 ) -u (x, t n ) ∆t .
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 2.6(b), which will be called in this work as the Method of Horizontal Lines (MOHL). Figure 2.6 presents the discretization difference between the MOVL and the MOHL. Next section aims at detailing the formulation and the appropriate method to solve Equation (2.3). 

Problem reformulation as a BVP

Consider the dimensionless nonlinear diffusion Equation (2.1). For our purposes it will be more advantageous to rewrite it in a non-conservative form. Thus, after taking one derivative, Equation (2.1) becomes:

∂u ∂t -λ (u, u x ) • ∂u ∂x -ν (u) • ∂ 2 u ∂x 2 = 0 , (2.15)
where coefficients λ and ν are defined as:

λ (u, u x ) def := Fo • k (u) • u x c (u) , ν (u) def := Fo • k (u) c (u) .
Equation (2.15) is equivalent to Equation (2.1) for smooth solutions. Indeed, under assumptions that coefficients c (u) and k (u) are positive defined, the solutions are smooth. Equation (2.15) is an evolution equation written in time-space: first, we write the time derivative, then, the ones related to space. According to the philosophy of the Method of Horizontal Lines, let us rewrite it in space-time:

ν (u) • ∂ 2 u ∂x 2 + λ (u, u x ) • ∂u ∂x = ∂u ∂t .
In other words, we see that the time derivative u t as a source term for our convenience here. Finally, we rewrite the last equation as a system of first-order differential equations in the spatial variable x :

       u x = ϑ , ϑ x = u t ν (u) - λ (u, ϑ) ν (u) • ϑ .
To have a true BVP, we need to eliminate the dependence of the solution regarding the time. For this reason, the solution u is semi-discretize in time, by replacing the surface u (x, t) by a sequence of horizontal lines:

u n (x) def := u (x, t n ) .
Thus, at each time layer t n we have to solve a true BVP in space:

       u n x = ϑ n , ϑ n x = (ũ t ) n ν (u n ) - λ (u n ) • ϑ n ν (u n ) • ϑ n ,
where (ũ t ) n is an approximation of the time derivative using some difference formulae. Depending on the desired accuracy the following backward finite-difference formulas can be considered:

ũ t n = u n -u n-1 ∆t + O ∆t , ũ t n = 3 • u n -4 • u n-1 + u n-2 2 • ∆t + O ∆t 2 .
In our work, the second-order accuracy in time is implemented in the calculations to have a more precise solution, so that we can use ∆t 10 -1 . The final scheme is unconditionally stable according to the construction of the method.

Implementation

As the MOHL is implemented in the Matlab TM environment, the solver bvp4c [START_REF] Shampine | Solving boundary value problems for ordinary differential equations in MATLAB with bvp4c[END_REF] is employed to solve the BVP problem at every time step. There are also other solvers that can be employed, such as the bvp5c [START_REF] Shampine | Solving boundary value problems for ordinary differential equations in MATLAB with bvp4c[END_REF] and even bvp6c [START_REF] Hale | A Sixth-Order Extension to the MATLAB Package bvp4c of J. Kierzenka and L. Shampine[END_REF]. The codes for bvp4c & bvp5c are available within any standard Matlab TM distribution, while the code bvp6c was developed by Dr. Nick Hale (Stellenbosh University, South Africa) 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 (I) approximation of the uniform fourth, fifth or sixth orders of accuracy in ∆x, respectively. Special attention should be given to the fact that other implementation details are significantly different among the solvers, the order is not the only difference among them. In our numerical simulations we use fourth-order adaptive methods, which are enough in most practical applications.

Let us consider a system of ordinary differential equations of the form u (x) = f (x, u) , within the interval [ a, b ] subjected to two boundary conditions ψ (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 (x). Algorithm 1 presents the procedure of implementation of the MOHL. 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, readers may refer to [START_REF] Shampine | Solving boundary value problems for ordinary differential equations in MATLAB with bvp4c[END_REF] and [START_REF] Hale | A Sixth-Order Extension to the MATLAB Package bvp4c of J. Kierzenka and L. Shampine[END_REF].

The essential feature of these algorithms is the adaptive distribution of collocation nodes. The grid adaptation and error controls are based on the residuals of the 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 BVPsolution 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 provided as one of the inputs of the solver. No special treatment needs to be carried out before supplying it.

Important Features

Mesh refinement is a great issue in the field of numerical simulations of thermal science problems. Therefore, the present method shows that the unsteady problem of nonlinear diffusion transfer can be solved from a different perspective. Namely, we interchange the order of discretization by proposing a numerical model accurate to the orders Algorithm 1 MOHL's algorithm for the weakly coupled heat and moisture transfer problems. 

[ u n , u n x ] = bvp4c (ODE of u n , BC of u n , u 0 ); 8:
Compute refined solution;

9:

Update: u 0 ;

10:

Increment: t := t + ∆t and n := n + 1 ;

11: return u (x, t) and u x (x, t) .

O(∆x 4 ) . . . O(∆x 6 ) with a moderate increase of computational efforts. Moreover, we took advantage of the fact that today there are well-tested robust adaptive numerical methods to tackle BVPs in one space dimension [START_REF] Hale | A Sixth-Order Extension to the MATLAB Package bvp4c of J. Kierzenka and L. Shampine[END_REF][START_REF] Shampine | Solving boundary value problems for ordinary differential equations in MATLAB with bvp4c[END_REF]. 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 [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF]. 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).

Extended Application

The work An adaptive simulation of nonlinear heat and moisture transfer as a boundary value problem (Gasparin et al., 2018a) presents the diffusion processes applied to the heat and moisture transfer through multilayered porous building materials. To demonstrate the benefits of this approach, three case studies are presented.

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 diffusion 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 O(10 -4 ) or less. The advantage of the proposed method is the adaptive spatial 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-differences 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.

The Quasi-Uniform Nonlinear Transformation

The Quasi-Uniform Nonlinear Transformation (QUNT) is a numerical method with an adaptive spatial grid which provides a high resolution where it is required [START_REF] Blom | Algorithm 731: A moving-grid interface for systems of onedimensional time-dependent partial differential equations[END_REF][START_REF] Khakimzyanov | On supraconvergence phenomenon for second order centered finite differences on non-uniform grids[END_REF][START_REF] Khakimzyanov | Numerical simulation of conservation laws with moving grid nodes[END_REF]. However, this approach is completely different from the Method of Horizontal lines. The adaptivity is obtained by using a monitor function which indicates the location of the spatial grid points during the simulation. The adapted grid is obtained by means of the equidistribution principle. This moving grid approach has the advantage of being conservative in space and also second-order accurate (or more: 4th order).

A finite-difference scheme on a fixed uniform mesh

Before solving the problem on a moving mesh, one has to chose a robust scheme on a fixed grid. This scheme will be then generalized to incorporate the motion of the mesh. Thus, the dimensionless diffusion equation is fully discretized with a uniform grid by using the IMEX scheme:

c n j • u n+1 j -u n j ∆t = Fo ∆x •   k n j+ 1 2 • u n+1 j+1 -u n+1 j ∆x -k n j-1 2 • u n+1 j -u n+1 j-1 ∆x   .
(2.16)

This scheme considers an IMplicit discretization for the field u n+1 and an EXplicit discretization for the storage and diffusivity coefficients, c n and k n . This approach is also known as the semi-implicit scheme and it approximates the continuous operator to order O(∆x 2 + ∆t). The advantage of the semi-implicit scheme over the fully implicit is to avoid sub-iterations in the solution procedure and, at the same time, being stable and consistent. For simplicity, the boundary conditions are not treated here. In what follows, the method based on this scheme with the moving grid is described.

Finite-difference scheme on a moving grid

Consider the computational domain I = [ 0 , l ], the reference domain Q = [ 0 , 1 ] and a bijective time-dependent mapping from Q to I:

x : Q × IR 0 -→ I .
This represents the quasi-uniform nonlinear transformation. In addition, the boundary points are required to map into each other:

x (0 , t) = 0 and x (1 , t) = l .
The reference domain Q is uniformly discretized into N elements, with q j = j h, j ∈ Z 0 representing the nodes of the grid, which is equally spaced h = 1/N . However, only the image of nodes q j under the map x (q , t) are needed, since they constitute the nodes of the moving mesh:

x (q j , t n ) = x n j .
The diffusion Equation (2.1) is rewritten on the domain Q, with the help of the composed function v (q , t)

def := (u • x) (q , t) ≡ u x (q , t) , t : J (q , t) • c (v) • ∂v ∂t = Fo • ∂ ∂q   k (v) J (q , t) • ∂v ∂q   + c (v) • ∂v ∂q ∂x ∂t ,
where

J def := ∂x ∂q
is the Jacobian of the transformation x (q , t). We assume that J > 0 , ∀ (q, t) . Therefore, the fully discrete form of the diffusion heat equation on a moving mesh is:

J n j • c n j • v n+1 j -v n j ∆t = Fo h •   k n j+ 1 2 J n j+ 1 2 •   v n+1 j+1 -v n+1 j h   - k n j-1 2 J n j-1 2 •   v n+1 j -v n+1 j-1 h     + + c n j •   v n+1 j+1 -v n+1 j-1 2 • h   •   x n+1 j -x n j ∆t   .
In summary, this is the parametrization process of the mesh motion by a bijective mapping.

In the sequence, the construction of this map is explained.

Construction and motion of the grid

The equidistribution method is used to construct the adaptive grid. To control the distribution of the nodes, a positive valued function w (x, t) has to be chosen. This function is called the monitor function. The choice of the monitor function and its parameters are important for the accuracy of the scheme. In particular, for unsteady nonlinear problems, there are several options for w. One can find more information, for example, in [START_REF] Stockie | A moving mesh method for one-dimensional hyperbolic conservation laws[END_REF]. In this work, the function chosen is based on [START_REF] Khakimzyanov | Numerical simulation of conservation laws with moving grid nodes[END_REF] and can be written as:

w (x, t) = 1 + α 1 |u| β 1 + α 2 • ∂u ∂x β 2 , (2.17)
where α 1 , α 2 , β 1 and β 2 are positive real parameters. The quantity w (x , t) assumes large positive values in areas where the solution has important gradients, inducing the elements of the grid mesh to concentrate there. The discrete formulation of the monitor function presented in Equation (2.17) is written as:

w n j = 1 + α 1 u n j β 1 + α 2 • u n j+1 -u n j x n j+1 -x n j β 2 .
A non-uniform grid I is given if we construct the mapping x (q , t) : Q -→ I and evaluate it in the nodes of the uniform grid. To this end, the equidistribution method proposes that the desired mapping x (q , t) be obtained as a solution of a nonlinear parabolic problem, which is written as:

∂ ∂q w (x , t) • ∂x ∂q = β • ∂x ∂t , ( 2.18) 
with x (0, t) = 0 and x (1, t) = l as the Dirichlet-type boundary conditions. Equation (2.18) can be written in the discrete form as:

1 h •   w n j+ 1 2 • x n+1 j+1 -x n+1 j h -w n j-1 2 • x n+1 j+1 -x n+1 j h   = β • x n+1 j -x n j ∆t ,
with the same boundary conditions x n+1 0 = 0 and x n+1 N = l. The parameter β > 0 plays the role of the inverse diffusion coefficient and it controls the smoothness of nodes trajectories.

Initial grid generation.

As the problem depends also on time, it is of capital importance to obtain a high-quality initial mesh. The initial condition must be adapted to the new grid before starting the dynamical simulation. At t = 0 we compute the monitor function w (x, 0) of the initial condition u (x, 0). Then, the mapping x (q, 0) is determined as the solution of Equation (2.18), which is reduced to a simple second-order ordinary differential equation:

∂ ∂q w (x , 0) • ∂x ∂q = 0 ,
supplemented with the Dirichlet boundary condition. Thus, the finite-difference approximation of this latter differential equation is:

1 h •   w 0 j+ 1 2 • x 0 j+1 -x 0 j h -w 0 j-1 2 • x 0 j+1 -x 0 j h   = 0 ,
with the discrete boundary conditions x 0 0 = 0, x 0 N = l. The previous nonlinear system of equations is solved iteratively, with iterations initialized with a uniform grid as the first guess. Its solution satisfies the equidistribution principle: in areas where w 0 j+ 1 2 takes large values, the space between two neighboring nodes x 0 j+1 and x 0 j has to be inversely proportionally small. One can find more information in [START_REF] Khakimzyanov | On supraconvergence phenomenon for second order centered finite differences on non-uniform grids[END_REF].

Smoothing step

To ensure the absolute smoothness of the mesh motion, the monitor function w can be further filtered. This process, called smoothing step enable to enlarge the set of acceptable values of the parameters α and β of the monitor function. According to [START_REF] Khakimzyanov | Numerical simulation of conservation laws with moving grid nodes[END_REF], the following filter inspired by an implicit scheme can produce robust results:

wj+ 1 2 = w j+ 1 2 -σ • wj+ 1 2 + σ 2 • wj-1 2 + wj+ 3 2 , (2.19)
where σ is a positive smoothing parameter. To complete Equation (2.19), boundary conditions are taken as in the original problem:

w 1 2 = w 1 2 , wN-1 2 = w N -1 2 .
The smoothed discrete monitor function wj+

1 2 N -1 j=0
is then obtained by solving the linear System (2.19). Note that the smoothing operator is applied to the monitor function and not to the nodes.

On the choice of the parameters

The parameters α 1 , α 2 , β 1 , β 2 , β and σ are chosen according to the problem and the solution under consideration. To determine the optimal values one has to make some numerical experiments. Interested readers may consult [START_REF] Khakimzyanov | Numerical simulation of conservation laws with moving grid nodes[END_REF] for more details on this aspect.

Important Features

In this work, the Quasi-Uniform Nonlinear Transformation (QUNT) method has been presented as an innovative method to perform the unsteady diffusion transfer simulations. The method is based on a non-uniform adaptive grid technique that identifies were the spatial nodes must be placed. This method has the advantage of being conservative in space and also second-order accurate. A satisfactory accuracy is shown by using much less spatial nodes than traditional methods, by moving its nodes where the gradients of the field are higher. Another remarkable point of this method is that it is easy to be implemented, by adding extra differential equations without being too computational costly. Simulations have been efficiently carried out for different wall and roof configurations, showing that the innovative method efficiently provides a gain of ≈ 25% on the computer run time if compared with the Euler IMplicit-EXplicit (IMEX) scheme.

Extended Application

Another adaptive method is the Quasi-Uniform Nonlinear Transformation (QUNT), which has been presented in An innovative method to determine optimum insulation thickness based on non-uniform adaptive moving grids (Gasparin et al., 2019b). The optimum insulation thickness of buildings walls is determined by taking into account the wall orientation and the position of the insulation in Brazilian buildings, which is determined through a parametric study. The annual heating/cooling transmission loads have been calculated using weather data and transient heat flux, simulated over one year.

First, the method is compared with a reference solution to properly understand its features. The high accuracy of the solution induced the analysis of the physical behaviour. Then, we studied a case focused on comparing a wall with insulation and another one with no insulation, applied to the climatic conditions of one Brazilian city. This case proved the importance of the insulation in reducing the transmission loads. Then, a second case study is simulated to verify if the insulation is better inside or outside of the building wall. The last case considered analyzing the effect of wall orientation and the roof assembly. As Brazil is mostly located in the Southern Hemisphere, the North facade is the one receiving more solar radiation. In consequence, for the hottest climate cities, the Northfacing wall is the one that requires the thicker insulation layers, while, in the South region -Curitiba -, the North-facing wall requires thinner layers of insulation.

A basic economic analysis was presented, demonstrating that the implementation of an insulation layer appears to be a cost-effective measure to save energy consumption. The optimum insulation thickness for all the wall and roof configurations studied is to be between 2 cm and 7 cm. The city that requires less insulation is São Paulo, due to its mild climate, and the city that needs more insulation is Salvador, due to its temperature to be high above the thermal comfort values.

The Reduced Spectral method

In this section the reduced Spectral method is described. While the previous methods were classified as complete-original models the spectral method is so efficient that is treated as a reduced-order model. The method is presented for three different applications (i) one-dimensional transient diffusion, (ii) one-dimensional transient diffusion with an extra parameter and (iii) two-dimensional transient diffusion.

One-dimensional space

For convenience, the nonlinear diffusion Equation (2.1) is written in the non-conservative form and it has its spatial domain transformed to the canonical interval x ∈ [0, 1]

x ∈ [-1, 1]. Thus, the diffusion equation is rewritten as follows:

∂u ∂t -ν (u) • ∂ 2 u ∂ x 2 -λ (u) • ∂u ∂ x = 0 , (2.20)
where,

ν (u) def := 4 • Fo • k (u) c (u) , λ (u) def := 4 • Fo c (u) • ∂ k (u) ∂ x = 4 • Fo • k (u) c (u) • ∂u ∂x .
The Spectral method assumes that the unknown u (x, t) can be accurately represented as a finite sum [START_REF] Fornberg | A pratical guide to pseudospectral methods[END_REF]:

u (x, t) ≈ u n (x, t) = n i = 0 a i (t) ϕ i (x) , i = 0, 1, 2, . . . , n . (2.21) Here, {ϕ i ( x )} n i = 0 is a set of basis functions that remain invariant in time, {a i ( t )} n i = 0
are the corresponding time-dependent spectral coefficients and n represents the number of degrees of freedom of the solution. Equation (2.21) can be seen as a series truncation after N = n + 1 modes. The Chebyshev polynomials are chosen as the basis functions since they are optimal in L ∞ approximation norm [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF]. It should be noted that other bases can be used, such as the Fourier and Legendre polynomials but for the problems of this work the Chebyshev polynomials ϕ i ( x ) ≡ T i ( x ) are the best choice [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF]. A detailed information about the Chebyshev polynomials can be found in Appendix A.

Therefore, the expression of the derivatives in the Chebyshev basis are:

∂u n ∂ x = n i = 0 a i (t) ∂T i ∂ x (x) = n i = 0 ã i (t) T i (x) , (2.22a) ∂ 2 u n ∂ x 2 = n i = 0 a i (t) ∂ 2 T i ∂ x 2 (x) = n i = 0 ã i (t) T i (x) , (2.22b) ∂u n ∂t = n i = 0 ȧ i (t) T i (x) , (2.22c)
where ȧ def := ∂a ∂t according to Newton's notation. Note that the derivatives are reexpanded in the same basis function. As a result, coefficients {ã i (t)} and { ã i (t)} must be re-expressed in terms of coefficients {a i (t)}. The connection is given explicitly from the recurrence relation of the Chebyshev polynomial derivatives [START_REF] Peyret | Spectral methods for incompressible viscous flow[END_REF]:

ã i = 2 c i n-1 p = i + 1 p + i odd p a p , i = 0, 1, . . . , n -1, (2.23) ã n ≡ 0 , ã i = 1 c i n-2 p = i + 2 p + i even p p 2 -i 2 a p , i = 0, 1, . . . , n -2, (2.24) ã n-1 ≡ ã n ≡ 0 ,
with,

c i =    2 , if i = 0 , 1 , if i > 0 .
These relations can be written in a matrix form as:

à = D • A , and à = D • A ,
where

A = a 0 , a 1 , • • • , a n-1 , a n T , Ã = ã 0 , ã 1 , • • • , ã n-1 , ã n T , Ã = ã 0 , ã 1 , • • • , ã n-1 , ã n T .
Moreover, the elements of matrices D and D are deduced from equations (2.23) and (2.24). Thus, with the derivatives defined in Equation (2.22), the series can be substituted into Equation (2.20) to provide the residual:

R = n i = 0 ȧ i (t) T i (x) -ν n i = 0 a i (t) T i (x) • n i = 0 ã i (t) T i (x) -λ n i = 0 a i (t) T i (x) • n i = 0 ã i (t) T i (x) , (2.25)
which is considered as a misfit of the approximate solution. The purpose is to minimize the residual, so that the solution satisfies the governing equations. To this end, the residual is minimized via the Tau-Galerkin method, which requires that Equation (2.25) be orthogonal to the Chebyshev basis functions R , T j = 0 :

R , T j = 1 -1 R ( x , t ) • T j (x) √ 1 -x 2 dx = 0 , j = 0, 1, 2, . . . , n -2 .
(2.26)

The residual can be minimized by other spectral approaches as well. For more details, one can consult (Mendes et al., 2017, Chap. 6) As a result, the projected residuals are:

M • ȧ i (t) = G i, j • ã i (t) + Λ i, j • ã i (t) , (2.27)
where, M is a diagonal and the singular matrix (rank ( M ) = N -2) which contains the coefficients of the Chebyshev weighted orthogonal system. The matrix M has the following form:

M =               π π 2 0 . . . π 2 0 0 0              
, and, the elements with indices (i, j) are written as:

G i, j = 1 -1 ν n i=0 a i (t) T i (x) • T i (x) • T j (x) √ 1 -x 2 dx , Λ i, j = 1 -1 λ n i=0 a i (t) T i (x) • T i (x) • T j (x) √ 1 -x 2 dx ,
with indices i, j being the ones defined in Eqs. (2.21) and (2.26). By using the Chebyshev-Gauß quadrature, the integrals are also approximated by a finite sum [START_REF] Peyret | Spectral methods for incompressible viscous flow[END_REF] :

G i, j ≈ π m m k = 1 ν k T i (x k ) T j (x k ) , Λ i, j ≈ π m m k = 1 λ k T i (x k ) T j (x k ) ,
where,

ν k def := ν   n i = 0 a i (t) T i (x k )   , and λ k def := λ   n i = 0 a i (t) T i (x k )   .
In addition, x k are the Chebyshev nodes:

x k = cos   2 • k -1 2 • m • π   , k ∈ {1, 2, . . . , m} .
The value of m was determined empirically according to numerical investigations, which is considered as m = N + 5 in our computations (Gasparin et al., 2019a).

To complete the problem, the boundary conditions from Equations (2.3) and (2.4) are also written in the form of the residuals:

ω 1 = k n i = 0 a i (t) T i (-1) • n i = 0 ã i (t) T i (-1) - Bi L 2 • n i = 0 a i (t) T i (-1) -u L ∞ (t) , ω 2 = k n i = 0 a i (t) T i (1) • n i = 0 ã i (t) T i (1) + Bi R 2 • n i = 0 a i (t) T i (1) -u R ∞ (t) ,
with T i (-1) = (-1) i and T i ( 1 ) ≡ 1 [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF]. The division of the Biot number (Bi) by 2 comes from the domain transformation, which can be seen in details in Appendix B.1.

In this way, it is possible to compose the system of ODEs to be solved plus the two additional algebraic expressions regarding the boundary conditions. Finally, the system of Differential-Algebraic Equations (DAEs) has the following form:

M • Ȧ = G i, j • D • A + Λ i, j • D • A + b , (2.28)
where b is the vector containing the boundary conditions, previously defined by ω 1 and

ω 2 : b = 0, 0, 0, • • • , ω 1 , ω 2 T
Initial values of the coefficients {a i (t = 0)} are calculated by the orthogonal projection of the initial condition [START_REF] Canuto | Spectral Methods -Fundamentals in Single Domains[END_REF]:

a 0, i ≡ a i (0) := 2 π • c i • 1 -1 u 0 (x) • T i (x) √ 1 -x 2 dx , i ∈ {0, 1, . . . , n} , (2.29)
where, u 0 (x) is the dimensionless initial condition. Therefore, the reduced system of ODEs composed from Eqs. (2.28) and (2.29) can be solved for the time-dependent coefficients:

   M • ȧ i = f (a i , t) , a i ( 0 ) = a 0, i .
(2.30)

Different approaches can be used to solve the system of ODEs (2.30). The most straightforward solution is to apply a numerical integration scheme, with moderate accuracy. So, with an embedded error control and not so stringent tolerances, it can be done very efficiently. In this work, the Matlab TM environment was used to perform simulations, and the solvers ODE15s or ODE23t were used to solve the differential-algebraic system of equations.

The output are the vectors of spectral coefficients {a i (t)} n i=0 at different successive instants of time. Then, it enables to reconstruct the solution thanks to spectral representations.

Parametric approach

For the sake of simplicity and without loosing the generality, this method is first explained considering the linear diffusion equation:

∂u ∂t = α • ∂ 2 u ∂ x 2 , (2.31)
for t ∈ [ 0, τ ] and x ∈ -1 , 1 . The diffusivity coefficient is consider as another coordinate with α ∈ [α min , α max ]. Thus, instead of propagating a solution in time for a fixed α, we propagate a family of solutions that are parametrized by α , as illustrated in Figure 2.7. It is almost like solving two-dimensional problems. The difference is that we do not have the partial derivatives with respect to the third space parameter α in Equation (2.31). Therefore, we seek for a solution of u that depends on three parameters: space, time and diffusivity:

u : (x, t, α) ∈ [ -1, 1 ] × [ 0, τ ] × [ α min , α max ] -→ IR .
Before applying the spectral method, the interval [ α min , α max ] must be transformed to the canonical one [ -1, 1 ]. The description of the transformation can be found in the Appendix B.2. Thus, ᾱ represents the transformation:

[ α min , α max ] ᾱ -→ [ -1, 1 ] .
The Spectral method assumes that the unknown u which depends on (x, t, ᾱ) can be accurately represented as a finite sum:

u (x, t, ᾱ) ≈ u n m (x, t, ᾱ) = n i = 1 m j = 1 a ij (t) T i-1 (x) P j-1 (ᾱ) .
(2.32)

Here, {T i-1 (x)} n i = 1 and {P j-1 (x)} m i = 1 are sets of basis functions that remains constant in time, {a i j (t)} n, m i,j = 1 are the corresponding time-dependent spectral coefficients and n • m represents the number of degrees of freedom of the solution. Equation (2.32) can be seen as a series truncation after N = n and M = n modes. The basis function P is chosen to be the Chebyshev polynomials since they proved to be practical and efficient in our applications but others basis such as Legendre polynomials could be used.

Using the expression of the derivatives provided in Appendix A, the residual of the diffusion Equation (2.31) is:

R (x, t, ᾱ) = n i = 1 m j = 1 ȧ i j (t) -ᾱ • ã i j (t) T i-1 (x) P j-1 (ᾱ) , (2.33)
which is considered as a misfit of the approximate solution. The purpose is to minimize the residual:

R (x, t, ᾱ) 2 -→ min ,
which is realized via two methods, the Tau-Galerkin and the Collocation methods. For this, we require that the residual of Equation (2.33) be orthogonal to the Chebyshev basis functions R , T k = 0 at the collocation points ᾱ p :

R , T k = 1 -1 R (x, t, ᾱ p ) • T k (x) √ 1 -x 2 dx = 0 , with, k ∈ {0, 1, 2, . . . , n -2} , and p ∈ {1, 2, . . . , m} .
Namely,

1 -1   n i = 1 m j = 1 ȧ i j (t) -ᾱ p • ã i j (t) P j-1 (ᾱ p ) T i-1 (x) T k (x) √ 1 -x 2   dx = 0 . (2.34)
This gives us (n -2)•m equations. The extra coefficients are obtained by substituting the derivative (A.1a) into the boundary conditions (2.3) and (2.4):

n i = 1 m j = 1 ã ij (t) T i-1 (-1) P j-1 (ᾱ p ) - Bi L 2 • n i = 1 m j = 1 a ij (t) T i-1 (-1) P j-1 (ᾱ p ) + Bi L • u L ∞ (t) 2 = 0 , (2.35a) - n i = 1 m j = 1 ã ij (t) T i-1 (1) P j-1 (ᾱ p ) - Bi R 2 • n i = 1 m j = 1 a ij (t) T i-1 (1) P j-1 (ᾱ p ) + Bi R • u R ∞ (t) 2 = 0 , (2.35b)
Equations (2.35a) and (2.35b) supplies the 2 • m equations to complete the system of a total of n • m equations. Therefore, the original problem is reduced to a differential-algebraic system of equations, which has the following form:

M (t, ᾱ) • ȧ l = f (t, ᾱ, a l ) , l ∈ {1, 2, . . . , n • m} , (2.36)
where, M ∈ Mat nm×nm (IR) is the mass matrix and f (t, ν, a l ) ∈ IR nm is a vector. The matrix {a i j } is transformed into a vector a i j (t) -→ a l = (i-1)n+j (t) so that the System (2.36) can be easily solved. Initial values of the coefficients {a ij (t = 0)} n,m i,j=1 are calculated by the projection of the initial condition as in Equation (2.29). Thus, the time-dependent coefficients {a ij (t)} n,m i,j=1 are computed by solving the following system:

   M (t, ᾱ) • ȧ i j = f (t, ᾱ, a i j ) , a ij (0) = a 0, ij .
(2.37)

This system is similar to Equation (2.30) and can be solved in the same way with the standard Matlab TM solvers.

Two-dimensional space

Consider the extension of the linear isotropic diffusion Equation (2.2) for the twodimensional space already transformed to the canonical interval:

∂u ∂t = α x • ∂ 2 u ∂ x 2 + α y • ∂ 2 u ∂ ȳ 2 , (2.38)
where α x = 4 α/l 2 x and α y = 4 α/l 2 y . For the two-dimensional case, four extra equations are given for the boundary conditions:

∂u ∂ ȳ ȳ = -1 = 0 , (2.39) ∂u ∂ ȳ ȳ = 1 = 0 , (2.40) u| x = -1 = u L ∞ (t) ,
(2.41)

u| x = 1 = u R ∞ (t) , (2.42)
in this case, Neumann-and Dirichlet-type. To complete the problem, an initial condition must be provided:

u (x, ȳ, t = 0) = u 0 (x, ȳ) .
(2.43) By using the spectral approach, the solution is approximated by the sum [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]:

u (x, ȳ, t) ≈ u n m (x, ȳ, t) = n i = 1 m j = 1 a i j (t) T i-1 (x) T j-1 (ȳ) .
(2.44)

Here, {T i-1 (x)} n i = 1 and {T j-1 (ȳ)} m i = 1 are sets of basis functions (Chebyshev polynomials) and {a i j (t)} n, m i,j = 1 are the corresponding time-dependent spectral coefficients which are the unknowns of the problem. Then, the residual is composed by substituting Equation (2.44) into Equation (2.38):

R (x, ȳ, t) = n i = 1 m j = 1 ȧ i j (t) T i-1 (x) T j-1 (ȳ) -α x • n i = 1 m j = 1 ã x i j (t) T i-1 (x) T j-1 (ȳ) -α y • n i = 1 m j = 1 ã y i j (t) T i-1 (x) T j-1 (ȳ) , (2.45)
The residual R (x, ȳ, t) is minimized via the collocation method:

R (x k , ȳ l , t) = 0 ,
which requires that the residual be equal to zero at the Chebyshev-Gauss-Lobatto points:

x k = -cos π • k n -1 , k ∈ {0, 1, 2, . . . , n -1} , (2.46a) ȳ l = -cos π • l m -1 , l ∈ {0, 1, 2, . . . , m -1} . (2.46b)
These points are chosen in order to minimize the error and to prevent against Runge's phenomenon [START_REF] Fornberg | A pratical guide to pseudospectral methods[END_REF]. Figure 2.8 presents the distribution of equally spaced points. However, their distribution follows Equation (2.46). To solve the problem, the residual can be written in a matricial form. For this, consider the spatial derivatives approximated as in [START_REF] Liu | Efficient Chebyshev spectral method for solving linear elliptic PDEs using quasi-inverse technique[END_REF]:

∂ 2 u n m ∂ x 2 = n i = 1 m j = 1 T i-1 (x k ) T j-1 (ȳ l ) ã x i j (t) = C (x k ) ⊗ C (ȳ l ) • D x ⊗ Id y • A i j (t) ,
and

∂ 2 u n m ∂ ȳ 2 = n i = 1 m j = 1 T i-1 (x k ) T j-1 (ȳ l ) ã y i j (t) = C(x k ) ⊗ C(ȳ l ) • Id x ⊗ D y • A i j ,
where Id x is the identity matrix of size n × n with respect to x, Id y is the identity matrix of size m × m with respect to ȳ and the operation denoted by ⊗ is the Kronecker product [START_REF] Kolda | Tensor decompositions and applications[END_REF].

The matrix of spectral coefficients {a i j } n,m i,j=1 was transformed into a vector:

A i j = a 11 , a 12 , . . . , a 1n , a 21 , a 22 , . . . , a 2n , a m1 , a m2 , . . . , a mn (nm×1) .

In addition, D x and D y correspond to the second order differentiation matrices with respect to x and y respectively such that D x ∈ Mat n×n (IR) and D y ∈ Mat m×m (IR), which are deduced from the recurrence Relation (A.2b). Finally, the vector of the Chebyshev polynomial values are:

C (x k ) = T 0 (x k ), T 1 (x k ), T 2 (x k ), . . . , T n (x k ) (1×n) , C (ȳ l ) = T 0 (ȳ l ), T 1 (ȳ l ), T 2 (ȳ l ), . . . , T m (ȳ l ) (1×m)
.

Therefore, the original problem is reduced to a system of ordinary differential equations, which has the following form:

C(x k ) ⊗ C(ȳ l ) • Ȧ i j = C(x k ) ⊗ C(ȳ l ) • α y • Id x ⊗ D y + α x • D x ⊗ Id y • A i j , (2.47)
which is valid only for the internal collocation points:

k = 1, 2, . . . , n -2 and l = 1, 2, . . . , n -2 .
For the boundary conditions, four extra equations appear to complete the system:

n i = 1 m j = 1 ã y i j (t) T i-1 (x k ) T j-1 (1) = 0 , (2.48a) n i = 1 m j = 1 ã y i j (t) T i-1 (x k ) T j-1 (-1) = 0 , (2.48b) n i = 1 m j = 1 a i j (t) T i-1 (1) T j-1 (ȳ l ) -u L ∞ (t) = 0 , (2.48c) n i = 1 m j = 1 a i j (t) T i-1 (-1) T j-1 (ȳ l ) -u R ∞ (t) = 0 , (2.48d)
which has the following matrix compact form:

C(x k ) ⊗ C(-1) • Id x ⊗ D y • A i j = 0 , (2.49a) C(x k ) ⊗ C(1) • Id x ⊗ D y • A i j = 0 , (2.49b) C(-1) ⊗ C(ȳ l ) • A i j -u L ∞ (t) = 0 , (2.49c) C(1) ⊗ C(ȳ l ) • A i j -u R ∞ (t) = 0 , (2.49d)
where D y corresponds to the first order differentiation matrix with respect to y, such that D y ∈ Mat (m × m), which is deduced from the recurrence Relation (A.2a).

Initial values of the coefficients {a ij (t = 0)} def := b ij are calculated by approximating the initial condition u 0 (x, ȳ) as:

u 0 (x, ȳ) = n i = 1 m j = 1 b ij T i-1 (x) T j-1 (ȳ) , (2.50)
which is expanded and solved for b ij using the collocation points. Therefore, the time-dependent coefficients {a i j } n,m i,j=1 are computed by solving the following the ODE system:

   M • Ȧ i j = f (t, a i j ) , a ij (0) = b ij . (2.51)
where, M ∈ Mat nm×nm (IR) is the mass matrix and f (t, a i j ) ∈ IR nm is a vector. To solve the System (2.51), the technique described in (Boyd, 2000, Chap. 10) is employed. We compose the matrix M and the vector f (t, a i j ) by the lines of the grid. Algorithm 2 describes the process to compose f (t, a i j ) which is similar to M . The integration in time is performed with the use of the solver ode15s as in the previous spectral applications.

Algorithm 2 Spectral-Collocation's algorithm for composing f (t, a i j ) and M.

1: f := zeros (n • m, 1) 2: M := zeros (n • m, n • m) 3: for k = 0 do Collocation point x 0 = -1 4:
for l = 0 : m -1 do 5:

f [l] := C(-1) ⊗ C(ȳ l ) • A i j -u L ∞ (t) Equation (2.49c) 6: for k = 1, 2, . . . , n -2 do Internal x k collocation point 7:
for l = 0 do 8:

f [(m -1)k + 1] := C(x k ) ⊗ C(-1) • Id x ⊗ D y • A i j Equation (2.49a) 9: for l = 1, 2, . . . , m -2 do 10: f [(m -1)k + l] := C(x k ) ⊗ C(ȳ l ) • α y • Id x ⊗ D y + α x • D x ⊗ Id y • A i j 11:
Right side of Equation (2.47)

12: M [(m -1)k + l] := C(x k ) ⊗ C(ȳ l )
Left side of Equation (2.47)

13:

for l = m -1 do 14: f [(m -1)k + (m -1)] := C(x k ) ⊗ C(+1) • Id x ⊗ D y • A i j Equation (2.49b) 15: for k = n -1 do Collocation point x m-1 = 1 16: for l = 0 : m -1 do 17: f [l + (m -1) • (n -2)] := C(+1) ⊗ C(ȳ l ) • A i j -u R ∞ (t) Equation (2.49d
) 18: return f (t, a i j ) and M.

Important Features

The main advantage of the Spectral-ROM is that the lower number of degrees of freedom needed to solve problem (2.31) is much lower than the ones needed to solve by means of conventional methods such as finite differences, finite elements and finite volumes. The Spectral approach assumes a separated tensorial representation of the solution by a finite sum of function products. It fixes a set of spatial basis functions to be the Chebyshev polynomials and then, a system of ordinary differential equations is built to compute the temporal coefficients of the solution using the Tau-Galerkin technique. The sum of polynomials generated suits for the whole domain, providing a high approximation of the solution which does not depend on the number of spatial nodes. The smoother the function is, the faster is the convergence of its spectral series [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF]. Actually, for considerably smooth problems, the error decreases exponentially, making the solution with the same order of accuracy of other methods but with a much lower number of degrees of freedom. As a result, this method has a low memory usage, allowing to store and operate a lower number of variables [START_REF] Trefethen | Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations[END_REF].

Nonetheless, spectral methods still have some constraints. For example, complex geometries are one of their main drawbacks as they work better when the geometry of the problem is fairly smooth and regular [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF], while finite-element methods are particularly well suited to problems in very complex geometries. Moreover, spectral methods can offer higher accuracy in geometries like boxes and spheres, which can be combined into more complex shapes [START_REF] Canuto | Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics[END_REF][START_REF] Fornberg | A pratical guide to pseudospectral methods[END_REF]. In applications where geometry-related disadvantages are not present, the classic finite-element, finite-volume and finite-difference methods do not come close in terms of efficiency.

Extended Application

In the first part of the work with the Spectral method Solving nonlinear diffusion problems in buildings by means of a Spectral Reduced-Order Model (Gasparin et al., 2019a) we simulate the moisture diffusion in one-dimension through porous materials as a Reduced-Order Model (ROM). The efficiency of the Spectral approach is demonstrated for simple and multilayered domains with highly nonlinear properties with sharp boundary conditions and profiles of solutions. Results show that the Spectral reduced-order model approach enables to simulate accurately the field of interest. Furthermore, numerical gains become particularly interesting for nonlinear cases since the proposed method can drastically reduce the computer run time, e.g. by a factor of 100, when compared to the traditional Crank-Nicolson scheme for one-dimensional applications.

Moreover, in A comparative study of two reduced order models for moisture diffusion problems in building physics [START_REF] Gasparin | Advanced reduced-order models for moisture diffusion in porous media[END_REF], we explore in details the capabilities of two model-reduction techniques -the Spectral Reduced-Order Model (Spectral-ROM) and the Proper Generalised Decomposition (PGD) -to numerically solve moisture diffusion transfer through porous materials. Both approaches are applied to three different problems to provide clear examples of the construction and use of these reduced-order models. Linear and non-linear unsteady behaviours of one-dimensional moisture diffusion are investigated and also a parametric problem in which the solution depends on space, time and diffusivity.

Results have highlighted that both model-reduction techniques supply accurate solutions and enable to reduce significantly the order of the model but the convergence of the spectral method is more efficient.

To complement the studies in one dimension with the Spectral-ROM we extended it to a coupled heat and moisture transfer problem. In A Spectral method for solving heat and moisture transfer through consolidated porous media (Gasparin et al., 2018e), the method is used to compute one-dimensional heat and moisture diffusion transfer in porous materials. Some adaptations in the method must be carried out due to the coupling between equations and the nonlinear coefficients. To demonstrate the benefits of this approach, three case studies were presented. The first one considers nonlinear heat and moisture transfer through one material layer. The second case includes the rain effect over two layers domain, while the last one compares the numerical prediction against experimental observations to appreciate the reliability of the model. Results show how the nonlinearities and the interface between materials are easily treated with the Spectral reduced-order model. For the reliability part, predictions show a good agreement with experimental results, which confirm robustness, computationally efficiency and high accuracy of the proposed approach for predicting the simultaneous heat and moisture transfer in porous materials. In this part, the methods described in Chapter 2 are applied to solve nonlinear partial differential equations that represent some transport phenomena in building physics. This study compares those numerical methods among each other and against the classical implicit and explicit schemes with the central finite-difference approaches, as well as with the reference solution. The performance of those numerical methods is presented in terms of accuracy, stability, convergence and computational time. Three case studies are considered on the comparison campaign: (i) one-dimensional nonlinear transfer, (ii) a parametric problem and, (iii) a two-dimensional linear transfer.

Numerical experiment 1: 1D nonlinear transport

This case study considers an isothermal one-dimensional nonlinear moisture transfer through a single layer material, as illustrated in Figure 3.1. The mathematical model is first defined, followed by the parameters of the problem and the results. 

Mathematical Model

The physical problem considers one-dimensional moisture transfer through a porous material defined in the spatial domain x ∈ I x ≡ [ 0, ] and in time domain

t ∈ I t ≡ [ 0, τ ] ,
where [m] is the length of the material and τ [s] is the final time of simulation. The moisture transfer occurs due to vapour diffusion, which mathematically formulated as (Gasparin et al., 2018b;[START_REF] Rouchier | Identification of the hygrothermal properties of a building envelope material by the covariance matrix adaptation evolution strategy[END_REF]:

c m • ∂P v ∂t = ∂ ∂x   k m • ∂P v ∂x   , (3.1)
whose field of interest is the vapour pressure is the the total moisture transfer coefficient under vapour pressure gradient defined as:

P v : (x, t) -→ P v (x,
k m : (P v ) -→ ρ l • R v • T P v • k l + δ v ,
where

δ v : P v -→ δ v (P v ) [s] and k l : P v -→ k l (P v ) [s]
are the vapour and liquid permeabilities of the material. The water vapour gas constant R v [J/(kg • K)], the temperature T [K] and the water specific mass ρ l [kg/m 3 ] are considered as constants. Table 3.1 presents the values of the water properties considered in this work. The saturation pressure expression was interpolated for the range of temperature between [ 273.15, 298.15 ] K, following the Antoine law to ensure the maximum accuracy.

The coefficient c m [kg/(m 3 • Pa)] is the the moisture storage coefficient and it is defined as:

c m : P v -→ w P s ,
where w : P v -→ w (P v ) [kg/m 3 ] is the material moisture content. The relation between the moisture content w and the relative humidity φ is given by the sorption isotherm, which depends on the material properties. In addition, the following assumptions are adopted in this study: (i) no hysteresis effect; (ii) no temperature dependency on the mass balance equation and (iii) properties are dependent only on the vapour pressure field. Boundary conditions. The vapour transport at the boundaries occurs due to the vapour pressure difference between the ambient air at the wall surface. The mathematical model is written as a Robin-type boundary conditions:

k m • ∂P v ∂x = h L m • P v -P L v, ∞ -g l, L , for x = 0, t 0 , -k m • ∂P v ∂x = h R m • P v -P R v, ∞ , for x = , t 0 ,
where P v, ∞ [Pa] stands for the vapour pressure of the ambient air, h m [s/m] is the convective transfer coefficient, g l [kg/(m 2 • s)] is the liquid flow (driving rain). In addition, the superscript L represents the left boundary side and the superscript R represents the right boundary. Furthermore, the ambient fields are considered as time dependent:

P v, ∞ : t -→ P v, ∞ (t) , g l : t -→ g l (t) ,
and the convective transfer coefficients are positive constants h m = const > 0 .

Initial condition. The initial distribution of the vapour pressure is considered uniform within the material:

P v (x, t = 0) = P v, i .
Flow. The density of the moisture flow g [kg/(m 2 • s)] is an important output and it is computed as:

g : t -→ -k m • ∂P v ∂x x 0 ,
where x 0 ∈ I x . Equation (3.1) with the initial and boundary conditions pose an important difficulty in dealing with the nonlinearities of the storage c m and diffusion k m coefficients: both of them depend on the moisture content field. These coefficients are usually given by empirical functions inferred from experimental data.

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 us also to estimate the relative magnitude of various terms in governing equations, and thus, eventually to simplify the problem using asymptotic methods [START_REF] Nayfeh | Perturbation Methods[END_REF]. In addition, the dimensionless form enables to manipulate numerically the quantities at the order of O(1) where the floating point arithmetics is designed to have minimal rounding errors [START_REF] Kahan | On a proposed floating-point standard[END_REF]. Considering the temperature range of interest in building applications, all transport coefficients are considered as moisture content dependent. Thus, the time and the space domains are transformed into a dimen-sionless representation:

x : I x -→ [ 0, 1 ] , t : I t -→ [ 0, τ ] , x -→ x , t -→ t t 0 ,
where the subscript 0 represents a reference value which is also valid for the next transformations and the superscript represents a dimensionless quantity of the same variable.

The dimensionless quantity of the vapour pressure is then defined as:

u : (x , t ) -→ P v P v, 0 ,
with the dimensionless properties given by:

c m : u -→ c m c m, 0 , k m : u -→ k m k m, 0 .
Therefore, the governing Equation (3.1) can be written in a dimensionless form as:

c m • ∂u ∂t = Fo m • ∂ ∂x   k m • ∂u ∂x   , ( 3.2) 
where Fo m = const > 0 is known as the Fourier number, which is defined as:

Fo m def := t 0 • k m, 0 2 • c m, 0 .
For the the boundary conditions, the dimensionless formulation is:

k m • ∂u ∂x = Bi L m • u -u L ∞ -g l, L , for x = 0, t 0 , (3.3) -k m • ∂u ∂x = Bi R m • u -u R ∞ , for x = 1, t 0 , (3.4)
with the Biot numbers Bi m = const > 0 being defined as:

Bi m def := h m • k m, 0 .
In addition, the ambient variables u ∞ and g l are written as:

u ∞ : t -→ P v, ∞ P v, 0 , g l : t -→ g l • k m, 0 • P v, 0 .
Finally, the dimensionless initial condition is:

u (x , t = 0) = 1 ,
for the reference vapour pressure chosen as the value of the constant initial condition

P v, 0 = P v, i .

Description of the case study

Moisture transfer is strongly nonlinear due to the variation of the material properties. For this reason, this case study is defined to investigate these effects. The case considers a 15-cm monolithic material, which properties are given in Table 3.2 and are graphically represented in Figure 3.2. The expressions of the properties are inspired from the wood fiber material, present in the paper of [START_REF] Rouchier | Identification of the hygrothermal properties of a building envelope material by the covariance matrix adaptation evolution strategy[END_REF]. The properties are given for a constant temperature of 23 • C on a isothermal case, which leads to a saturation pressure of 2804.4 Pa .

Property Expression Unit

Sorption isotherm, w (φ) 931 φ 5 -1575 φ 4 + 956.3 φ 3 -268.1 φ 2 + 55.32 φ -0.01074 kg/m 3

Vapour permeability, δ v (φ) 6.8 • 10 -10 + 5.446 • 10 -10 • exp (9 φ 4 -1.5) s

Table 3.2: Properties of the wood fiber material [START_REF] Rouchier | Identification of the hygrothermal properties of a building envelope material by the covariance matrix adaptation evolution strategy[END_REF]. The moisture storage and the moisture diffusivity coefficients are shown in Figures 3.3(a) and 3.3(b), respectively. Since the measurements were performed for a relative humidity on the range of φ = [ 0.25, 0.75 ] , the liquid transport is neglected. For the other values of relative humidity we have approximated the solution. The nonlinearities are considerably high in this case, as the increase of the moisture content increases of one order of magnitude for c m and three orders of magnitude for k m , correspondingly. Initial conditions are considered uniform over the spatial domain, with an initial vapour pressure of P v, i = 1168.9 Pa , corresponding to a relative humidity of 50 % . The boundary conditions, represented by the relative humidity φ are shown in Figure 3.4. They oscillate during the 168 hours of simulation as follows:

φ L ∞ (t) = 0.5 + 0.35 • sin 2 2 • π • t 48 • 3600 + 0.05 • sin π t 120 • 3600 , φ R ∞ (t) = 0.5 + 0.45 • sin 11 π • t + 70 168 • 3600 .
The convective mass transfer coefficients are set to h L m = 2•10 -7 s/m at the left boundary, and to h R m = 5 • 10 -7 s/m at the right boundary. This case aims at exciting the nonlinear properties of the material to induce sharp profiles of the vapour pressure. Dimensionless parameters. The reference time is t 0 = 1 h, thus, the final simulation time is fixed to τ = 168 . The reference vapour pressure was taken the same as the initial condition. The problem is considered with Fo m = 0.042 and P v, 0 = P v, i . At the boundaries, the Biot numbers assume the following values: Bi L m = 33.58 and Bi R m = 83.96 . The unitless vapour pressure varies sinusoidally over the time according to the following expressions:

u L ∞ (t ) = 1 + 0.7 • sin 2 2 • π • t /48 + 0.1 • sin π t /120 , u R ∞ (t ) = 1 + 0.9 • sin 11 (π • t + 70)/168 .
The properties are dimensionalized with the following reference values: c m, 0 = 0.0034 and k m, 0 = 8.9327 • 10 -10 .

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

c m (u) = 36.91 • u 4 -99.9 • u 3 + 90.99 • u 2 -34.01 • u + 7.018 , k m (u) = 1.368 + 9.478 • 10 -3 • exp (1.457 • u 3 ) .

Assessment performance for the numerical solution

All the numerical results in this thesis are computed using the Matlab TM environment (Matlab, 2017b). To analyse the accuracy of the proposed method, the error between the solution u num , obtained by one of the numerical methods, and the reference solutions u ref , are computed as functions of x by using the following Euclidean norm:

ε 2 ( x ) def := 1 N t • Nt n = 1 u n num ( x , t n ) -u n ref ( x , t n ) 2 ,
where N t is the number of temporal steps and u is the vapour pressure dimensionless field. The reference solution u ref (x, t) is computed by using the Matlab TM open source toolbox Chebfun [START_REF] Driscoll | Chebfun Guide[END_REF].

It is an open-source software system used for numerical computing. Chebfun is capable of dealing with linear and nonlinear differential and integral operators. By using pde solvers, it enables to compute a numerical solution of a partial derivative equation with the Chebyshev polynomials adaptive spectral methods.

The Chebfun solution is used as reference because it is not an efficient method. It has a good solution accuracy but for solving the type of problem of this work it has a low computational speed. Moreover, the global L ∞ error is given by the maximal values of ε 2 ( x ) :

ε ∞ def := sup x ∈ [ 0 , ] ε 2 ( x ) .
We state that if the error is above the maximum limit of ε ∞ 10 -2 , the solution provided by the numerical scheme is not considered to be accurate enough.

Another criterion of evaluation is the computational (CPU) time t cpu [s] spent by the numerical model to compute the solution. The measurement is carried out using the

Matlab

TM platform and a computer with a processor Intel R Core TM i5 @ 2.80GHz×4 and 16GB of RAM. The ratio R cpu is defined by: R cpu def

:= t cpu t • ,
where t • is a reference time. One should notice that the codes are not optimized, and the slight differences between the CPU time should not be considered.

For the flow, the Euclidean error (distance) is computed as:

ξ ( t ) def := g ref (x i , t) -g num (x i , t) 2 ,
where g ref is the reference flow computed with Chebfun solution and g num is the flow computed with other solutions. The global L ∞ error of the flow is given by the maximum value over time of ξ ( t ) :

ξ ∞ def := sup t ∈ [ 0 , τ ] ξ ( t ) .
In this document, when we refer to Implicit scheme that means the Euler method of finite-difference approximation with backward-time and centred-space (BTCS) discretization. For nonlinear problems, the Implicit scheme uses the fixed-point algorithm to perform sub-iterations at each time step to converge the solution. Here, N denotes the number of iterations performed by the algorithm. In addition, when we refer to Explicit scheme that means the Euler method of finite-difference approximation with forward-time and centred-space (FTCS) discretization.

Results and discussion

Simulations with the numerical methods described in the previous chapter are presented here. They are compared with the Implicit scheme and with the reference solution so that the accuracy can be evaluated.

DuFort-Frankel scheme

The classical Explicit scheme is generally avoided when dealing with building physics simulations, due to the Courant-Friedrichs-Lewy (CFL) [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF] stability condition that restrains the temporal discretization to prohibitively small time steps. By employing the Explicit scheme, the time step must satisfy:

∆t ∆x 2 2 • min u    c m (u) Fo m • k m (u)    .
(3.5)

With the given parameters of the case study, the CFL condition is ∆t 8.5 • 10 -4 , corresponding to a physical time step of 3.1 s , which is considerably low.

The DuFort-Frankel scheme is also an explicit scheme. However, it is unconditionally stable according to the Von-Neumann stability condition, as shown in preceding chapter on Section 2.2. To asses numerically this affirmation, the first simulation is computed using the DuFort-Frankel scheme, considering a spatial discretization parameter of ∆x = 10 -2 and a time step of ∆t = 10 -2 , which is 12 times larger than the CFL stability condition. Figure 3.5(a) presents the vapour pressure profiles at different times during the simulation. A good agreement between the reference solution and the DuFort-Frankel scheme is observed. It shows that the method provides a stable and reliable solution even for higher time steps than one given by the CFL restriction.

The time evolution of the fields at the surfaces of the material is presented in Figure 3.5(b). The vapour pressure at x = 0 m varies according to the frequency of the relative humidity at the left side, and the vapour pressure at x = 0.15 m varies according to the frequency at the right side. It is noticed that moisture diffuses faster when its content is higher within the material. The errors between the reference solution and the ones computed with the DuFort-Frankel and the Implicit schemes are given in Figure 3.6. For the Implicit scheme, the same discretization parameters where considered. For the iterative process, a tolerance of η 10 -5 has been used for the convergence of the sub-iterations, using a fixedpoint algorithm. It confirms that both numerical schemes enable to compute an accurate solution, at the order of 10 -3 , with the Implicit being slightly more accurate than the DuFort-Frankel scheme by considering the same spatial and temporal discretizations. However, the time necessary for simulating the solution with the Implicit scheme is 20 times higher than the one with the DuFort-Frankel approach, as reported in Table 3.3. The implicit scheme requires around 9 sub-iterations per time step to treat the nonlinearities of the problem. The DuFort-Frankel approach computes directly the solution and therefore has a reduced computational cost, which is only 5 % of the Implicit scheme based code.

A parametric analysis of the behaviour of the three numerical schemes namely DuFort- Frankel, Implicit and Explicit schemes is carried out for different values of the temporal discretization ∆t , while the spatial discretization is maintained to ∆x = 10 -2 . Results of the error ε ∞ are shown in Figure 3.7(a) and the respective computational times for each simulation is given in Figure 3.7(b). As expected, the Explicit scheme was not able to compute a solution when the CFL stability condition is not respected. The value computed from Equation (3.5), around ∆t 8•10 -5 is in accordance with the results from the convergence study. It also confirms that the DuFort-Frankel scheme is unconditionally stable, as it computes a solution for any discretization parameter ∆t . However, the error grows faster than those for the Implicit scheme. By using the DuFort-Frankel method, the maximum value for the time step that provides an accurate solution is ∆t = 5 • 10 -2 , equivalent to 3 min. For time steps greater than this, the error of the solution is not acceptable. Even though, if the Implicit is used with a time step of ∆t = 2 • 10 -1 ≡ 12 min the computational time to perform this simulation is still higher than the one used by the DuFort-Frankel scheme with a smaller time step of ∆t = 2 • 10 -2 ≡ 1.2 min, which provides a solution with the same order of accuracy ε ∞ O 5 • 10 -3 . These values are indicated by double-arrow in the Figures 3.7(a) and 3.7(b). The problem with the Implicit scheme is that when we increase the time step, the number of sub-iterations also grows, making the cost at each time step higher. This phenomenon is illustrated in Figure 3.8, presenting the average number of sub-iterations required by the fixed-point algorithm to perform the simulation for the respective time step. A second parametric analysis is also performed to evaluate the behaviour of the DuFort-Frankel scheme, by modifying the discretization values of the temporal and spatial grids, as presented in Figure 3.9(a). When ∆t is higher than ∆x , the inconsistency in the DuFort-Frankel scheme becomes apparent and the error starts to grow. This happens, since the DuFort-Frankel scheme modifies the original parabolic equation to a hyperbolic equation, in which the term τ • ∂ 2 u ∂t 2 is added to the original equation. To reduce this term influences as much as possible the solution, the parameter τ must be small, according to the consistency analysis in Equation (2.7). For this reason, when τ ∆t > 1 , the error begins to grow, as shown in Figure 3.9(b). This figure presents the error as a function of τ ∆t , which was computed with a constant value for the diffusivity α := Fo m = 0.042 . If large time steps are required, then one must also use large spatial grid. However, a balance is required, since there is a limit on the accuracy related to large time and spatial steps. For example, for ∆t = 5 • 10 -2 , equivalent to 3 min , the ideal spatial step would be around ∆x = 4 • 10 -2 in order to have an error of ε ∞ O 10 -3 . As illustrated, the choice of the discretization parameters has strong influences on the computed solution. When analysing the error of the DuFort-Frankel scheme as a function of the time discretization ∆t , three regions have been highlighted. The first one corresponds to small discretization parameter ∆t 10 -3 , where all schemes provide accurate solutions with different computational costs. The second region, ∆t ∈ [1 • 10 -3 , 5 • 10 -2 ], where the DuFort-Frankel scheme is highly recommended due to the relation between accuracy and computational time. Finally, the last region, where ∆t 5 • 10 -2 , which includes larger values of ∆t , in which the Implicit scheme is suggested for its accuracy but does not necessarily imply a reduction on the computational 58 cost.

Another feature regarding too large time steps ∆t is that the schemes do not succeed in representing the physical phenomena. As mentioned in (Gasparin et al., 2018b), the time step has to be carefully chosen in accordance with the characteristic time of the physical phenomena. Moreover, in building simulation, the recommended time step is not higher than 10 min (Dos [START_REF] Santos | Analysis of numerical methods and simulation time step effects on the prediction of building thermal performance[END_REF], which can make the DuFort-Frankel scheme interesting for such applications.

QUNT

The second numerical method applied in this study is the Quasi-Uniform Nonlinear Transformation (QUNT), in which the spatial grid is no longer uniform but adaptive. Simulations using the QUNT approach have considered the following values for the parameters of the monitor function: α 1 = 0.5 , β 1 = 0.8 , α 2 = 0.7 and β 2 = 1.5 . These parameters have been determined by our numerical investigations and correspond to accurate predictions for this case study. These parameters varies between ( 0, 2 ) for simulations of this kind of problem. For the grid diffusion parameter, an acceptable value tested was β = 1 ∈ ( 0, 100 ) and for the smoothing parameter was σ = 0.01 ∈ ( 0, 10 ). Simulations are performed for a time discretization of ∆t = 1 • 10 -2 and a reduced number of spatial points N x = 26 the equivalent to ∆x = 4 • 10 -2 . For initial comparisons, these parameters are fixed and, then, a parametric study will be performed. The QUNT method is first compared to the reference solution and with the Implicit scheme, since it is the main classical approach used in the literature. For the Implicit scheme, the same discretization values are considered.

The sample trajectory of the adaptive spatial nodes can be observed in Figure 3.10. The spatial points concentrate where the values of the gradient are higher, particularly at the boundaries, where important variations are imposed. As the vapour pressure within the material rises significantly around the 3rd day it forces the grid nodes to significantly move due to the high gradients. In the rest of the simulation, the motion of the modes is determined by the oscillation frequency associated to the left boundary condition.

Figure 3.11(a) presents the computed solutions of the vapour pressure profile for the time instants of t = 24 h, t = 48 h and t = 84 h. As it can be observed, the spatial nodes of the QUNT method does not coincide with the uniform grid but both methods provide stable and accurate solutions. In the Figure 3.11(b), the evolution of the vapour pressure fields is presented for x = 0 m and x = 0.15 m, the boundaries of the material. The grid points located at these places are the only ones not to move. They are fixed during the whole simulation time. The results for the error ε 2 over the time are presented in Figure 3.12, for the QUNT, the IMEX and the Implicit approaches. The IMEX was added to the comparison since the QUNT method is based on the IMEX, but with an adaptive spatial grid. The moving grid modifications have improved the method, which makes it to handle better with the high gradients than the original IMEX method, providing more accurate results, especially at the boundaries. This figure shows that the QUNT and IMEX methods have a spatial accuracy of order of O(10 -4 ), while the Implicit scheme has an accuracy of order of O(10 -3 ), for the same spatial and time discretization. The maximal error of each solution computed is reported in Table 3.4. Thus, for the same number of spatial nodes (N x = 26), the QUNT can be 4 times more accurate than the Implicit scheme and 2 times faster. However, the modifications in the adaptivity of the grid has increased the computational time. For this simulation, the IMEX method takes half of the simulation time required by the QUNT method. A convergence study was carried out to compare the QUNT, the IMEX and the Implicit schemes to evaluate the accuracy and the computational time of the solutions. Thus, for different values of time steps ∆t = 10 -3 , 10 -2 , 10 -1 , the global error ε ∞ and the computational time t cpu are computed as a function of the number of spatial nodes N x ∈ [ 6 , 101 ]. Figures 3.13, 3.14 and 3.15 display the results.

Figures 3.13(a), 3.14(a) and 3.15(a) present the error as a function of the number of spatial nodes, showing that the QUNT method is more accurate among these methods, followed by the IMEX and then by the Implicit one. For the QUNT method, as the number of spatial nodes increases, the global error decreases faster when compared to the other methods. Besides, for N x = 6 all methods provide solutions of the field with approximately the same order of accuracy, which is insufficient to get good approximations of the solution. The QUNT approach proved to be very sensible regarding variations of ∆t.

As ∆t decreases, the moving grid adapts better because the solution becomes practically smoother.

The difference between the IMEX and the Implicit schemes relies on the computation of the nonlinear properties at each time iteration. For the IMEX scheme, the coefficients c m and k m are computed at the time t n while for the Implicit scheme, they are computed at the time t n+1 . This difference avoids the sub-iterations in the IMEX scheme. Thus, the simulations with the Implicit method are slower than the one obtained by the IMEX as noticed in Figures 3.13(b), 3.14(b) and 3.15(b). The only situation in which the Implicit method comes close to the QUNT in terms of efficiency is when there are smaller spatial and time steps, such as ∆t = 10 -3 and N x = 11. In which the Implicit solution is 1.5 times faster than the QUNT and just 1.4 times less accurate, but still acceptable. These small discretization values make the Implicit scheme to require less sub-iterations.

One interesting feature of the QUNT method is that it provides accurate solutions for a lower number of spatial nodes than the other two methods. For example, when ∆t = 10 -1 only 20 spatial nodes are enough to provide a fast solution. Furthermore, the solution is more accurate than the Implicit scheme even with 100 spatial nodes. However, with small quantities of spatial nodes, the computation of the flow may lack the accuracy, as represented in Figure 3.16. A minimum quantity of points is needed to have a good approximation of the fluxes, which is around N x 20 for the QUNT method and N x 30 for the other two methods.

Figure 3.16 shows that the computation of the flow is less accurate than the computation of the field. This happens since the error increases when the derivative has to be computed. Depending on how the derivative is approximated we can lose one order of accuracy. The QUNT method also provides solutions for the flux more accurate than the IMEX method. To obtain even higher accuracy on the flow, it is necessary to increase the derivative approximation order. In this case, the error of the flux is relatively high due to the nonlinearity of the case study. For a less nonlinear case, the error of the flux density can be considerably reduced. The horizontal dashed line indicates the acceptable limit for the accuracy of the flow. Above this line, the error is not acceptable. Therefore, it is clearly shown that the QUNT approach computes an accurate solution of the problem with more efficiency than the Implicit scheme. Moreover, the adaptivity of the spatial grid increases the accuracy of the solution, making the method to converge faster with fewer spatial nodes than methods with uniform grids. The only problem with this method is that it requires the determination of some parameter values, that are not evident at the first moment. However, this problem can be overcame with gain of experience in using the method.

MOHL

The Method of Horizontal Lines is also adaptive in space but unlike the Quasi-Uniform Nonlinear Transformation method the quantity of spatial nodes may vary during the simulation time, because a boundary value problem is solved at each time step. Simulations for this case study with the Method of Horizontal Lines (MOHL) are performed using the solver bvp4c with the relative and absolute tolerances both set to tol = 10 -4 . The time step discretization is of ∆t = 10 -1 while the space step is adaptive and depends on the solver. The initial guess of the vapour pressure field is composed of 20 spatial nodes. For the Implicit solution, simulations are carried with the same time step ∆t = 10 -1 but with a fixed spatial discretization of ∆x = 10 -2 .

The vapour pressure distribution profiles for t = { 24, 48, 84 } h are presented in Figures 3.17(a), for the selected mesh, which is determined by the solver, for the reference solution and for the refined solution. This last one is an evaluation of the solution that covers all the interval [ 0, l ] in the chosen points, which are not the same as the selected mesh of the solver. A good agreement can be observed between the MOHL method and the reference solution. Simulation results of the vapour pressure evolution on the boundaries of the building component (x = 0 m and x = 0.15 m ) are presented in Figure 3.17(b). The proposed method has been able to follow successfully the variations 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 evolution. The vapour pressure diffuses slowly over the material, since the low influence from the opposite boundary can be experienced. However, as the material adsorbs the moisture, faster it is diffused, following the indicated properties.

Note that at each time layer, the number of spatial grid points of the selected mesh varies, since the distribution of spatial nodes is adaptive, with a minimum of 20 and a maximum of 178 points, as can be observed in Figure 3.18(a). If the tolerance of the solver is modified, the number of nodes of the spatial grid would change accordingly. 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). Figure 3.18(b) shows the location of the points during the simulation. As with the QUNT method, the spatial nodes are concentrated close to the left boundary as the vapour pressure inside the material is increased. Then, after the third day, the number of spatial nodes and its location stabilises as the step at the right boundary decreases. The error ε 2 between the MOHL and the Implicit methods against the reference solution are given in Figure 3.19, which is to the order of O 10 -3 for both methods. The error of the solution with the MOHL is consistent with the chosen tolerances that was set to tol = 10 -4 and to the discretization on time, that is second-order accurate O(∆t 2 ) . It is admissible that the actual error is larger than the prescribed tolerance. The solver estimates the residual and not the actual error stricto sensu (since the exact solution is unavailable). It is expected, however, that the error decreases when the tolerance decreases as shown in Table 3.5.

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 3.5 presents the error ε ∞ , as well as the number of average mesh points used by the solver and the computational time to perform the simulations. As the tolerances become more stringent, the solver adapts the solution making the number of mesh points to increase, in order to provide a consistent solution with the specified accuracy. The initial number of points chosen was 20 in all cases. For the tolerance tol 10 -3 , the solutions are considered to be accurate. However, for tol 10 -5 the computational time required to perform simulations rises significantly and becomes prohibitive. Thus the tolerances recommended would be 10 -4 tol 10 -3 . The computational time required to compute the solution of this case study is provided in Table 3.6. For the Implicit solution, it takes on average 11 iterations at each time step to converge to the solution which increases the global computational time of the method. However, the MOHL method could not be faster than the Implicit approach. In fact, to solve the boundary value problem at each time step, makes the method slower but more accurate and robust. In addition, there is the interpolation of the solution for the desired grid. Thus, small time steps are not recommended since the computational time would be too high.

MOHL Implicit

∆t 2 • 10 -1 10 -1 N x [20, 58] 101 ε ∞ 3.78 • 10 -3 3.09 • 10 -3 t cpu 27 s 7 s
Table 3.6: Features and computational time of the MOHL and Implicit scheme.

One of the advantages of using the solver bvp4c is that it has as output the respective derivative of the field, making the computation of the flow straightforward. Thus, Figure 3.20 presents the total moisture flow (liquid plus vapour) at both boundaries. 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 oppose boundary. The step of relative humidity at the right boundary can also be observed on the moisture flow which suddenly changes. The advantages of the proposed method are the adaptive spatial mesh grid according to the solicitations of the physical phenomena and, as the resulting discretization is fully implicit, the method is unconditionally stable. This feature is highlighted in Figure 3.21. The MOHL is more stable for large time steps than the Implicit method. In other words, the semi-discretization in time is not subject to any kind of restriction regarding the time step. Thus, the time step can be chosen based on the accuracy considerations solely. Typically, the MOHL can compute a solution with the same order of accuracy as the Implicit approach with a time step three times larger.

Spectral

The previous numerical methods tested needed to solve large systems of equations to compute the solution. The spectral approach considers a reduced system of equations which can drastically reduce the computational time and it also provides very accurate solutions.

The discretized method provides a reduced system of ordinary differential equations (ODEs). If before the methods had to solve a system with approximatively 100 degrees of freedom per time step, with the Spectral method it is reduced to approximatively 10. The spectral coefficients {a n ( t )} are calculated for every intermediate time instant by the solver ode15s. This solver is set with the absolute and relative tolerances to tol = 10 -4 . The inputs are the initial time, the final time and the time step (optional) and then the solver supplies the integration at the given time.

Similarly to what has been presented for the previous methods, the profiles of the vapour pressure and its evolution at the boundaries are given in Figures 3.22(a) and 3.22(b), in order to verify the convergence of the solution. It can be seen the physical phenomenon is well represented and the Spectral solution is in good agreement with the reference one. The number of modes used in the discretization was N = 11.

The error ε 2 of the Spectral and of the Implicit methods is given in Figure 3.23. Both solutions of the problem have been computed with the same discretization parameters, with ∆x = 10 -2 and ∆t = 10 -2 , which provide an order of accuracy around O(10 -4 ) for the two solutions. However, the Implicit scheme needs more time to perform simulations 32.3 s while the Spectral method needed only 1.6 s. To compare the efficiency of the methods employed, Table 3.7 presents some numerical features: (i) the computational time, (ii) the average number of sub-iterations and (iii) its respective error. The Spectral method is substantially faster than the other methods. It represents only 5 % of the CPU time needed for the Implicit approach, to obtain a solution with the same order of accuracy.

The first three and last three coefficients are shown in Figures 3.24 final solution with the Chebyshev polynomials. It can be seen that the first coefficients have the most significant values with the highest magnitudes, and as the order of the coefficients rises, their value decreases. The last coefficients are always the smallest ones.

When dealing with Chebyshev polynomial, the magnitude of the last spectral coefficient determines the error upper limit, acting as an error estimator to the solution [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF].

The spectral solution has been projected over a grid with N x = 101 spatial nodes so the results can be analysed. However, as the Spectral method provides a global solution to the problem, the number of spatial nodes does not influence the accuracy of the final Spectral N = 11 1.6 5 -6.36 • 10 -4

Table 3.7: Computational time required to perform simulations (∆x = 10 -2 and ∆t = 10 -2 ). result. Figure 3.25(a) can corroborate this statement. The spectral solution remains with the same accuracy as the number of spatial nodes increases, while the Implicit scheme gets more accurate, since it is a method based on a local approximation. In fact, to get more accurate results with the spectral method, one needs to increase the number of modes or to set more strict tolerances on the ODE solver. For each simulation, the computational time has been measured and is given in Figure 3.25(b), which proves once more the Spectral method is very efficient. The computational time almost does not change with the increase of spatial modes which compose the spectral solution.

To analyse the influence of the number of modes on the spectral solution, Figure 3.26 presents the error ε 2 for the Spectral solution computed using different number of modes, for the discretisation parameters set to ∆x = 10 -2 and ∆t = 10 -2 . As the number of modes increases, the solution of the Spectral method becomes more accurate. With 9 modes we already have a satisfactory solution to the problem, with the absolute error of the order of O(10 -3 ) , while the solution with 8 modes is not enough. However, if the 71 problem is linear, the number of modes required can be even further reduced. In Gasparin et al. (2019a), in the linear case, 6 modes were enough to compose a sufficiently accurate solution. The other factor influencing the Spectral solution is the tolerance of the solver. For this reason, a convergence study is performed to verify the efficiency of the proposed method. The discretization parameters are the same as in the previous study, while the number of modes N of the Spectral solution and of the tolerance of the solver vary. Figure 3.27(a) presents the maximum error ε ∞ as a function of the number of spectral modes. As we increase the number of modes, the solution is more accurate, converging exponentially as expected by the theoretical results [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF]. The smoother a problem is, the faster is its convergence. One can observe that the solution converges to a minimum value after a certain number of modes, which is related to the tolerance of the ODE solver. In addition, the computational time to perform each spectral simulation is presented in Figure 3.27(b). The computational effort to perform the simulation increases with the number of modes. However, it remains extremely low. Figure 3.28 displays the Fourier power spectrum of the signal frequency per unit of time, at the left and right boundaries and at the middle of the material. The picks of energy are higher on the left side of the material, which generates the highest gradients. Some oscillations due to aliasing errors may occur, as the ones around 5•10 -1 Hz. However, the power of this frequency is very low compared to the highest peak. It corresponds to a difference of 7 orders of magnitude, which does not cause problems to the solution.

The analysis of error ε ∞ as a function of the time step ∆t is shown in Figure 3.29. As shown in the previous study, the implicit scheme is first-order accurate in time. However, the way the spectral method is implemented does not allow it to be evaluated. Since the ODE15s was employed, it does not matter how large the time step is, because the solver will use the adapted one that is convenient to meet the tolerance and then supply the solution at the time instant required by the user.

Therefore, the spectral reduced-order model has shown a high accuracy and perfect agreement with the respective reference solutions. The advantage of the proposed method is the low computational burden which makes it very efficient.

Reduction of computational efforts

When different numerical methods must be compared, one option is to have a solution with the same order of accuracy and, then, evaluate the computational time required. Therefore, Table 3.8 provides the computational time required to compute each solution for the numerical method indicated. Simulations have been performed with different discretization values for each numerical method, in order to give a solution with the same order of accuracy, in this case ε ∞ O(10 -3 ). The Reference solution is used to evaluate the accuracy of the numerical methods while the Implicit scheme is used as a reference to the computational run time since it is the most widely applied method in the building physics domain.

Reference Implicit DuFort-Frankel QUNT MOHL Spectral N = 10 ∆t 10 -1 10 -1 3 • 10 -2 2 • 10 -1 2 • 10 -1 2 • 10 -1 ∆t [
To tell which method is the best, it depends on the circumstances and on the demand of the user. For the case presented, to get this order of accuracy, three methods have shown a significant reduction of computational efforts: the DuFort-Frankel method, the QUNT and the Spectral method. Even if the time step of the DuFort-Frankel method is not as large as the other methods, it provides a solution with the same order of accuracy and still reduces computational costs. However, if one still needs to use large time steps, other methods such as the QUNT, the MOHL and the Spectral one can achieve the goal with a satisfactory accuracy. It is important to point out that the MOHL requires a great computational effort compared to the other methods. Due to the sharp gradients of this problem, at each time step, the bvp solver requires much more time to calculate the converged solution.

The efficiency of a numerical method is defined here as the combination of accurate solutions with the minimum of computational efforts. If one wants to reduce computational costs by simply increasing the time step value, one must pay attention to the minimal accuracy required. Reductions on computational costs due to large time steps are penalized with loss of accuracy. In this way, the methods presented here can be attractive alternatives to reduce computational costs without losing accuracy.

Numerical experiment 2: a parametric problem

The issue of solving a parametric problem is to compute the solution of P v depending not only on the time and space coordinates x and t, but also depending on another parameter of the problem, such as the moisture diffusivity coefficient α. Thus, the family of solutions sought are:

P v : I x × I t × I α -→ IR , (x, t, α) -→ P v (x, t, α) ,
where α is defined as a coordinate of the problem within a given interval α ∈ I α ≡ [ α min , α max ] , in which N α is noted as the number of elements (the cardinal) of the domain.

By using the spectral method, two approaches can be adopted to solve this kind of problem. The first one is to compute the solution P v (x, t) for each value of the parameter α , by means of a loop in the algorithm. Other classical numerical methods may also use this approach. The second option is to see the problem as if it were a two dimensional space problem, but without any boundary conditions or derivatives regarding this parameter in the governing equation. Thus, the solution is sought as a decomposition on each coordinate of the problem (as presented in Section 2.5.2).

Mathematical Model

Assume that the transport coefficients of Equation (3.1) are positive constants:

c m = const > 0 and k m = const > 0 .
Then, the moisture transfer equation can be written in a simplified way:

∂P v ∂t = α • ∂ 2 P v ∂x 2 , (3.6)
where

α def := k m /c m [m 2 /s] is the diffusivity coefficient.
For the parametric problem, α ∈ I α and the vapour pressure field is computed as

P v : (x, t, α) -→ P v (x, t, α) .
Boundary conditions. The boundary conditions that complete the problem are expressed as:

k m • ∂P v ∂x = h L m • P v -P L v, ∞ , for x = 0, t 0 and ∀ α ∈ I α , -k m • ∂P v ∂x = 0 , for x = , t 0 and ∀ α ∈ I α ,
recalling that h m [s/m] is the convective transfer coefficient, which is constant in time h m = const > 0, that L represents the left boundary side and that P L v, ∞ [Pa] is the vapour pressure of the ambient air which varies with time as a known function:

P L v, ∞ : t -→ P L v, ∞ (t) .

Initial condition.

A uniform vapour pressure distribution is considered as the initial condition within the material:

P v (x, t = 0, α) = P v, i .

Dimensionless representation

As discussed in the previous case study, it is important to get a dimensionless formulation of the problem to be simulated. Therefore, let us consider the following dimensionless quantities:

u def := P v P v, 0 , u ∞ def := P v, ∞ P v, 0 x def := x , t def := t t 0 , α def := α • t 0 2 , Bi m def := h m • k m ,
where the subscript 0 denotes a reference value, chosen according to the problem under consideration and the superscript represents a dimensionless quantity of the same variable. Therefore, Equation (3.6) can be written in a dimensionless form as:

∂u ∂t = α • ∂ 2 u ∂x 2 , for x ∈ ] 0, 1 [ , t > 0 and α ∈ I α * = [α min , α max ] .
The dimensionless formulation of the boundary conditions are:

∂u ∂x = Bi L m • u -u L ∞ , for x = 0, t 0 and ∀ α ∈ I α * , - ∂u ∂x = 0 , for x = 1, t 0 and ∀ α ∈ I α * ,
where we have the Biot number and the dimensionless ambient vapour pressure given by the mappings:

Bi L m : α -→ Bi L m (α ) and u L ∞ : t -→ u L ∞ (t ) .
Lastly, the dimensionless initial condition is:

u (x , t = 0, α ) = 1 .
Therefore, the dimensionless problem is solved by means of the spectral method in the case study described next.

Description of the case study

For this case, we seek for a parametric solution of Equation (3.6). The vapour pressure is computed as a function of time t , space x and the moisture diffusivity coefficient α of different materials. The left boundary is exposed to cyclic changes of relative humidity between 33 % and 75 % , with a 24 h period as presented in Figure 3.30. The total time of simulation is four days, the equivalent to 96 h . The convective vapour coefficient is set to h v = 2 • 10 -8 s/m . The right boundary is set impermeable, with a null flow. The initial condition is homogeneous among all materials with a vapour pressure value of P v, i = 925.46 Pa , the equivalent to 33 % of relative humidity. Simulations undergo at a constant 23 • C temperature that leads to a saturation pressure of 2804.4 Pa. All materials have the same 10-cm length with the properties given in Table 3.9. The moisture diffusivity varies within the segment Ω α = [1.97•10 -10 , 2.19•10 -8 ] m 2 /s. This case study is inspired in the one presented by [START_REF] Rode | The Concept of Moisture Buffer Value of Building Materials and its Application in Building Design[END_REF], including the material properties, to study the so-called moisture buffer effects. Table 3.9: Properties of the materials used in the parametric case [START_REF] Rode | Moisture Buffering of Building Materials[END_REF]. [ 200, 133.33, 133.33, 80, 28.57, 66.67 ], in this order from α 1 to α 6 . The vapour pressure varies sinusoidally over the time following the mapping:

u L ∞ (t ) = 1 + 1.3 • sin 2 (2 • π • t /48) .
The dimensionless properties of each material are:

α 1 = 7.084 • 10 -5 , α 2 = 1.514 • 10 -4 , α 3 = 1.848 • 10 -4 , α 4 , = 8.141 • 10 -4 , α 5 = 1.979 • 10 -3 , α 6 = 7.887 • 10 -3 .

Results and discussion

Simulations are performed for 6 different values of moisture diffusivity coefficients in the segment Ω α , representing the materials of Table 3.9. The Spectral-Parametric model is compared with the Implicit scheme, with the Spectral method using a loop over the parameter values and also with the reference solution, obtained with the Chebfun toolbox.

Figure 3.31(a) shows the profiles of vapour pressure within each material at the time instant t = 48 h. It is possible to observe significant variations of the vapour pressure within each material as the value of the moisture diffusivity coefficient increases. Furthermore, Figure 3.31(b) presents the evolution of vapour pressure near the left boundary, at x = 1 cm . In these figures, the Spectral-Parametric approach is in a good agreement with the reference solution. However, some oscillations can be observed in Figure 3.31(a) only for α 1 . This occurs because the diffusivity coefficient is very low and the number of modes was not enough to support the stiffness of the problem.

In addition, with the simulated solution, it is possible to verify the variation of moisture content of each material. For this, Figure 3.32 presents the quantity of water in grams per kilogram of material, over the time. There are some materials that retain more moisture than others, which can help to determine the Moisture Buffer Effect (MBV). The material of the curve on the graphic corresponding to α min retains very low quantity of mass while the material of α max accumulates more the moisture according to the increase of ambient relative humidity.

To perform the parametric simulation, the Spectral-Parametric needed around N = 19 modes for the spatial base functions and M = 6 modes for the parametric base functions, with the tolerance of the solver ode15s set to tol = 10 -5 . The solution is projected on a grid composed by ∆x = 1 • 10 -2 and ∆t = 1 • 10 -2 . To have the same order of accuracy, the Implicit scheme was discretized with a time step equal to ∆t = 1 • 10 -3 and a spatial step of ∆x = 5 • 10 -3 . One can observe that the number of modes of the Spectral-Parametric approach is higher. Moreover, the discretization of the implicit scheme is very refined. This happened because the Aggregate concrete material has a very slow diffusion process with a low value of the diffusivity coefficient, which concentrates high gradients near the left boundary and slows down the moisture diffusion through the material. The Spectral-parametric model applies two different techniques of spectral methods, the Tau-Galerkin approach for the spatial basis function and the Collocation method for the parameter basis function. Both basis functions are the Chebyshev polynomials. The classical collocation points for the Chebyshev basis are the extrema points [START_REF] Fornberg | A pratical guide to pseudospectral methods[END_REF]. However, as the application is not the classical one, different collocation points have been tested and are presented in Table 3.10. As a matter of fact, the best collocation points are the value of the respective parameter to be simulated because during the simulation these actual values of the parameter are required. The error ε 2 is shown as a function of space for each material in Figure 3.33(a) for the Spectral-Parametric solution and in Figure 3.33(b) for the Implicit solution. The methods are built to give a minimum order of accuracy, around O( 10 -3 ) . As the value of the moisture diffusivity coefficient increases, the error decreases because the diffusion process becomes less stiff to solve. It can be noticed that the Implicit scheme faces difficulties to treat the oscillations at the left boundary while the Spectral parametric has a solution more homogeneous over space. The Spectral-parametric is also compared with the Spectral solution using a loop over the values of the parameter. Both solutions are computed using the same quantity of modes for the spatial basis N = 19 and the same tolerance of the solver tol = 10 -5 . The speed of calculation is higher for the spectral-parametric approach as reported in Table 3.11. These results are also commented in [START_REF] Gasparin | Advanced reduced-order models for moisture diffusion in porous media[END_REF]. Indeed, the number of modes in the Spectral-parametric approach increases by considering a new parameter as another coordinate, making the integration in time slower than computing the solution for each parameter. However, both approaches are able to compute the solution of the parametric problem faster than the reference solution and the Implicit scheme, for the same order of accuracy. 3.34(a) presents the last spectral coefficients of the solutions. The magnitude of the spectral coefficients is higher than the error of the solution O(10 -6 ) for α max and O(10 -3 ) for α min . It happens because the collocation points do not correspond for the optimal convergence of the solution. They are chosen according to the specific value of the parameter α k . Thus, the Theorem 6 (Boyd, 2000, Page 47) is not valid for this specific case. Nevertheless, if one works with an interval of variation of the parameter (and not determined values), the collocation points can be chosen in their optimal way. In this case, the order of magnitude of the last coefficients would be closer to the order of the error.

Figure 3.34(b) shows the Fourier power spectrum function of the signal frequency per unit of time, generated by the fast Fourier transform. The peak on the energy is due to the sinusoidal signal at the boundary, since it is the only parameter that interferes on the transient solution. The convergence of the basis function of the parameters is presented in Figure 3.35(a). As it can be observed, it drops at different rates with the increase of the number of modes. Up to 8 parameters the ideal number of modes or collocation points is the quantity of parameters itself. For more than 8 parameters, it does not worth to increase the number of modes because the computational efforts are not saved. The computational run time of each simulation is presented in model. One interesting result is that the Spectral-parametric model becomes more interesting when the solution has to be computed for more than 30 values of the parameter. The computational time of the Spectral-parametric and the Spectral-loop are almost the same. This occurs because the number of modes M stops growing after 8 parameters, which are enough to give a solution with a minimum order of accuracy around O(10 -3 ). Therefore, if one needs to compute a parametric solution with a large number of these parameters the Spectral-parametric model would be strongly recommended. The advantage of the Spectral-parametric method is related to its ability to compute at once the solution depending on the three coordinates, whereas the Spectral and the Implicit approaches compute the solution for each value of the desired parameter independently one-by-one with a loop. It should be noted that the Implicit approach is a low-order approximation that provides a solution less accurate than the one obtained by the spectral ones.

Tensor representation of the solution

As the solution of the problem depends on more than two parameters, an interesting option is to have a tensorial representation of the solution. By using a 3 rd order tensor, it can represent the solution with discrete points, as the one presented in Figure 3.37. In this figure, we have four slices of α and four slices of x for all values of t . The problem with the discrete tensor representation is that they demand a high CPU memory to work and to storage. Thus, a solution to this problem would be to decompose the tensor with low-rank approximations to compress the high-dimensional data. Remark that this process is performed after obtaining the solution using one of the numerical methods.

The most popular tensor decomposition method is the Canonical Polyadic Decomposition (CPD), also known as CANDECOMP-PARAFAC decomposition (CANonical DE-COMPosition or PARAllel FACtors model) [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of "Eckart-Young" decomposition[END_REF][START_REF] Harshman | Foundations of the PARAFAC procedure: Models and conditions for an "explanatory" multi-modal factor analysis[END_REF], which is similar to the Singular Value Decomposition (SVD) method for matrices. The SVD method decomposes the original matrix into a product of three other matrices, two unitary matrices and a diagonal matrix, reshaping the original matrix.

The idea of the Canonical Polyadic (CP) method is to decompose approximatively the third-order tensor U of size N x × N t × N α :

U = u i j k ∈ IR Nx×Nt×Nα ,
into the sum of rank-one tensors, such as a, b and c:

U ≈ Ũ = R r=1 λ r a r • b r • c r ,
where • denotes the vector outer product, and a r ∈ IR Nx , b r ∈ IR Nt and c r ∈ IR Nα for r = 1, 2, . . . , R [START_REF] Kolda | Tensor decompositions and applications[END_REF].

In this work, the toolbox developed by [START_REF] Bader | Matlab tensor toolbox version 3.0[END_REF] has been used to decompose the tensor of the parametric solution, with the function cp_als. This function is based on an alternating least-squares algorithm [START_REF] Phan | On fast computation of gradients for CANDE-COMP/PARAFAC algorithms[END_REF]. The original tensor file containing U has a size of ≈ 44.39 MB. After performing the CP decomposition, the size of the new file containing the decomposed tensor Ũ is of 1.11 MB, which is 40 times smaller, as can be seen in Table 3.12. The number of ranks used is R = 15, which provides a compressed tensor with errors ε ∞ = 1.52 • 10 -2 and ε 2 = 2.37 • 10 -3 . We can say that Ũ is a rank-15 approximation to U. If one wants to compress even more the tensor, it is possible to reduce the ranks value but the accuracy of the compressed tensor would be degraded. The computational effort to decompose the tensor is very low, in this case, for rank 15, it was of 0.91 s. The space required to store the solution computed with the implicit method is important, representing a file size of 44.39 MB. However, the time to recover one value of the solution u (x 0 , y 0 , α 0 ) is equal to 3 • 10 -4 s, which is almost instantaneous. By using the spectral method, the reduced matrix only needs 8.53 MB of storage. The computational time to compose one value u (x 0 , y 0 , α 0 ) scales with 2•10 -3 s. It is only six times greater than the Implicit one but the important advantage is the reduction on the storage. With the tensor obtained by the CP method, the storage of the solution is reduced to 1.11 MB and the computational time to compose the solution u (x 0 , y 0 , α 0 ) is 1.2 • 10 -3 s. Thus, in terms of post-processing the solution, the decomposed tensor is more advantageous to the Spectral method since the storage and the time to access an element are lower. The number of ranks to be chosen depends on the error of the compressed tensor, on the computational time to perform the compression and also on the size of the compressed data to be stored. According to Figure 3.38(a), the minimal number of ranks to have an acceptable accuracy would be 10. If the computational time to compress the tensor is taken into account, 10 ranks would be also good, because for lower rank, the solver needs more iterations to reach the desired tolerance. However, as the number of ranks is increased, the space needed to store the decomposed tensor also grows but it is still lower than the original tensor. Thus, based on the three graphics of Figure 3 .38, for R ∈ [10, 20] all values are good options, and the choice of the rank will depend rather on the desired accuracy.

Numerical experiment 3: 2D transport

For the one-dimensional cases, the numerical solutions are validated against a reference solution computed with the Chebfun package. Unfortunately, the Chebfun does not simulate yet two-dimensional transient phenomena, so that we need to find another way to validate the results. For this reason, two case studies are presented in this section. One simple case, to validate the codes with an analytical solution and a second case with more complex boundary conditions, which represents better the physics for a real case.

Mathematical Model

Consider the partial differential equation of the moisture diffusion transfer: 

c m • ∂P v ∂t = ∇ k m • ∇P v , ( 3 
c m = const > 0 and k m = const > 0 .
Thus, the two-dimensional linear diffusion over a plate is given by expanding Equation (3.7) to become:

∂P v ∂t = α • ∂ 2 P v ∂x 2 + α • ∂ 2 P v ∂y 2 , (3.8)
where

α def := k m /c m [m 2 /
s] is a constant diffusivity coefficient and the spatial domains are given by x ∈ [ 0, ] and y ∈ [ 0, h ]. The diffusion process is considered isotropic and the vapour pressure field of the plate, at position (x, y) and time t is denoted by P v (x, y, t) , which is depicted in Figure 3.39.

For the validation case, the vapour pressure value is imposed at the boundaries of the spatial domain using Dirichlet-type conditions: On the other hand, for the other case study, mixed boundary conditions are imposed:

P v (0, y, t) = P L v, ∞ , P v ( , y, t) = P R v, ∞ , P v (x, 0, t) = P B v, ∞ , P v (y, h, t) = P T v, ∞ , with P v, ∞ = const .
k m • ∂P v ∂x x = 0 = h L m • P v -P L v, ∞ (y, t) , P v ( , y, t) = P R v, ∞ (t) , ∂P v ∂y y = 0 = 0 , ∂P v ∂y y = h = 0 ,
with vapour exchange at the left boundary, null flow at the top and bottom boundaries and the prescribed value at the right boundary. The initial state of the plate is given by: P v (x, y, 0) = P v, 0 (x, y) , which must be compatible with the boundaries to seek for an analytical solution, so the problem can be strongly well-posed [START_REF] Kreiss | Initial-Boundary Value Problems and the Navier-Stokes Equations[END_REF].

Dimensionless representation

For the transformation of the original mathematical model, consider the following dimensionless quantities:

u def := P v P v, 0 , u ∞ def := P v, ∞ P v, 0 x def := x l , y def := y h , t def := t t 0 , α x def := α • t 0 l 2 , α y def := α • t 0 h 2 , Bi L m def := h L m • l k m ,
where the subscript 0 represents a reference value and the superscript represents the dimensionless quantity of the same variable. Thus, Equation (3.8) can be written in a dimensionless form as:

∂u ∂t = α x • ∂ 2 u ∂x 2 + α y • ∂ 2 u ∂y 2 . (3.9)
The dimensionless formulation of the boundary conditions for the validation case are:

u (0, y , t ) = u L ∞ , u (1, y , t ) = u R ∞ , u (x , 0, t ) = u B ∞ , u (x , 1, t ) = u T ∞
and for the other case are:

∂u ∂x x = 0 = Bi L m • u -u L ∞ (y , t ) , u (1, y, t) = u R ∞ (t ) , ∂u ∂y y = 0 = 0 , ∂u ∂y y = 1 = 0 .
Last, the dimensionless form of the initial condition is:

u (x , y , t = 0) = u 0 (x , y ) .

Assessment performance for the numerical solution

The error between the solution u num , obtained by one of the numerical methods, and the reference solutions u ref , are computed as functions of x and y by using the following root mean square error:

2 ( x, y ) def := 1 N t • Nt n = 1 u n num ( x, y, t n ) -u n ref ( x, y, t n ) 2 ,
where N t is the number of temporal nodes and u is the vapour pressure dimensionless field. The reference solution u ref (x, t) in this case is computed using an analytical solution.

Moreover, the global L ∞ error is given by the maximal values of ε 2 ( x, y ) :

∞ def := sup (x, y) ∈ [ 0, l ]×[ 0, h ] ε 2 ( x, y ) .
If the error does not respect the limit of ∞ 10 -2 , the solution provided by the numerical scheme is not accurate enough for our applications.

Another criterion of evaluation is the computational time t cpu [ s ] spent by the numer-ical model to compute the solution. The ratio R cpu is defined by: R cpu def

:= t cpu t • ,
where t • is the reference time.

Validation

In this section, the Spectral-2D is validated with an analytical solution since the Chebfun package does not simulate two-dimensional unsteady transport. Analytical solutions are very restrictive and complex problems cannot be solved by them. Therefore, this part considers a simplified two-dimensional transfer to validate the algorithms. Three types of numerical approaches are also tested to compare against the Spectral approach, (i) the Implicit scheme, (ii) the Alternating Direction Implicit (ADI) method [START_REF] Peaceman | The numerical solution of parabolic and elliptic differential equations[END_REF] and (iii) the Explicit scheme. The analytical solution. The analytical solution is generated by the Maple software, considering a the spatial domain [ 0, 1 ] × [ 0, 1 ], with a initial condition given by the mapping:

u 0 (x , y ) = x • (1 -x ) • y • (1 -y ) ,
and constant Dirichlet boundary condition:

u L ∞ = u R ∞ = u T ∞ = u B ∞ = 0 .
Figure 3.40 shows the initial condition u 0 (x , y ) . Therefore, the analytical solution is given by the expression:

u (x , y , t ) = 1 π 6 ∞ p = 1 ∞ q = 1 16 • (-1) p + (-1) q -(-1) p+q -1 • sin (qπx ) • sin (pπy ) • e -αt π 2 (q 2 +p 2 ) q 3 • p 3 .
To have a solution with 8 digits of precision, the analytical solution needs a truncation on the sum of around 10 modes for p and q.

Simulation tests. Simulations are performed with a dimensionless diffusivity coefficient of α = 0.02 for a time horizon of τ = 16, with a time discretization of ∆t = 10 -2 for the Implicit, ADI and Spectral method. Due to the CFL condition, the Explicit scheme is simulated with ∆t = 10 -3 . The spatial discretization parameters are the same for both directions, ∆x = 10 -2 and ∆y = 10 -2 . The Spectral-Collocation method is implemented with the same number of modes for both spatial bases, N = 11 and M = 11, making a total of N • M = 121 modes. In addition, the tolerance of the solver ode15s is set to tol = 10 -6 . This case study allows to see the drying process of the moisture occurring within the material. For this end, Figure 3.41(a) presents the value of the field u in the middle of the material (x , y ) = (0.5, 0.5) . As expected, with the course of time, the value of the dimensionless field tends to merge to the prescribed value at the boundary. All the numerical methods seem to be in accordance with the analytical solution at this location. In addition, Figure 3.41(b) presents the field u as a function of x , at y = 0.5 and for different instants t = [ 0, 2, 4, 16 ] . In fact, the diffusion process occurs until t = 16, where it reaches the steady state.

To better evaluate the solutions, the distribution of the error 2 among the x and y axes is presented in Figure 3.42, for each numerical solution implemented. The Implicit method provided the less accurate solution, with an order of O(10 -5 ) while the other methods provided a solution with an order of O(10 -6 ), which are satisfactory results. The solution of the implicit scheme looks different from the explicit because the system of equations to be solved is constructed from a different way. Normally, in building physics applications, the solutions do not need to be that accurate. In this case, as we are val- idating the numerical model, we sought very precise solution since the variation of u is on the order of O(10 -2 ). Therefore, it can be concluded that all numerical methods have converged to the appropriate solution.

The tensor corresponding to the Analytical solution was decomposed so it can be stored and used for the comparisons. The original tensor has a dimension of 101 × 101 × 1601 with a size of 124.6 MB. By decomposing it with an 8-rank approximation its dimension reduce to {8} {101 × 8} {101 × 8} {16016 × 8} with amounts to only 0.11 MB . The decompositions took 0.3 s to be computed and if one needs to compose the original tensor it is instantaneously done.

Figure 3.43(a) presents the last spectral coefficients of the solutions a i,m with i = 0, 2, . . . , n . In this case, the order of the last spectral coefficients is compatible with the error, which is of the order of ∼ O(10 -5 ∼ 10 -6 ). Figure 3.43(b) indicates the Fourier power spectrum function of the signal frequency per unit of time. It describes the distribution of power in frequency space that composes the signal. The spectrum of u (t) can contain some information about the nature of u . It is possible to verify that there is no extra signal forcing from the boundaries and that the problem is linear since there are no other oscillations on the signal frequency.

A convergence study is also performed for the two-dimensional methods. The spatial steps are fixed to ∆x = 2 • 10 -2 and ∆y = 2 • 10 -2 , while the time step varies. The error ∞ for the Explicit, Implicit and the ADI approaches as function of ∆t is presented in Figure 3.44. As expected, the Explicit scheme cannot provide a solution for ∆t higher than the CFL stability condition. In addition, the error of the Implicit scheme starts 94 to grow before the ADI, but both methods have the same order of accuracy O(∆t ) . Another interesting feature of the ADI method is that it is faster than the Implicit one, making the ADI method a better choice between the two-dimensional finite-differences schemes. In the next case study, only the ADI method will be used in comparison with the Spectral-Collocation. A convergence study is also performed for the Spectral-Collocation method but in relation to the total number of modes N •M . Since the collocation method is implemented, the minimum number of modes for each basis function must be 3, which gives a total of 9 modes. Thus, Figure 3.45 gives the error ∞ as a function of the number of total modes, considering different tolerances for the ode solver. The two-dimensional Spectral-Collocation method converges exponentially, at an order of ∼ O exp (C • N • M ) , where C is a constant. Between the tolerances, N • M = 10 • 10 modes is more than sufficient to have very accurate solutions and N • M = 7 • 7 modes can provide solutions of order of O (10 -4 ), which are appropriate for the type of problem that we are dealing with.

The computational time of the solutions to provide a solution with the same order of accuracy O 10 -6 is given in Table 3.13. The Implicit scheme is the slowest method between the ones implemented in our study. Indeed, it is slowed down by the inversion of the matrix A (NxNy×NxNy) on the solution of A • U = b at each time step. The modification of the ADI method makes it faster than the Implicit one. Thus, instead of inverting a very large matrix, it inverses two smaller matrices C (Ny×Ny) and D (Nx×Nx) to obtain the solution at each time step. Moreover, it is important to mention that the ADI method is also more accurate for the same discretization parameters. However, the Spectral-Collocation method is still more efficient compared to the other approaches. The number of degrees of freedom of the Spectral-Collocation is N • M , while for the finite-difference methods is N x • N y . Thus, the degrees of freedom of the Implicit method is 2.5 • 10 3 and of the Spectral collocation is 81, representing a reduction of 2 orders of magnitude at each 96 time step. Therefore, the next case study will be solved only by the Spectral-Collocation and the ADI methods since their efficiency proved to be the most satisfactory.

Analytic Explicit Implicit ADI Spectral (N M = 9 • 9) ∆t 10 -3 10 -3 10 -3 10 -2 10 -2 ∆x = ∆y 10 -2 10 -2 2 • 10 -2 2 • 10 -2 2 • 10 -2 tol - - - - 10 -6 ε ∞ - 3.

Description of the case study

For this case, the left boundary is considered as Robin-type, which increases the difficulty of the problem. The vapour pressure is computed as a function of time t within a plate of size 0.1×0.5 m . The material of the plate is the cellular concrete whose properties are given in Table 3.9. The left boundary is exposed to cyclic changes of relative humidity between 33 % and 89 % , with a 24 h period, which also changes according to the position y. The top and bottom boundaries are set impermeable, with null flow and the left boundary varies smoothly in a step. Figure 3.46(a) and 3.46(b) presents the variations of the ambient relative humidity at the left and right boundaries of the material, respectively. Due to the complexity of the boundary conditions, analytical solutions for this kind of problem are not simple and sometimes not possible to generate.

The total time of simulation is four days, the equivalent to 96 h ≡ 345600 s . The convective vapour coefficient at the left boundary is set to h v = 7 • 10 -8 s/m . The initial condition is homogeneous within the material with the a vapour pressure value of 925.46 Pa , the equivalent to 33 % of relative humidity. Simulations undergo at a constant 23 • C temperature which leads to a saturation pressure of 2804.4 Pa. Dimensionless parameters. The reference time is t 0 = 1 h ≡ 3600 s, thus, the final simulation time is fixed to τ = 96 . The reference vapour pressure was taken the same as the initial condition u 0 = 1. At the left boundary, the Biot number assumes the value Bi L m = 100. The vapour pressure varies sinusoidally over the time following the mapping:

u L ∞ ( y , t ) = 1 + 1.3 • sin π • t /24 2 + 0.4 • y , u R ∞ ( t ) = 1 + 1.9 • sin π • t /96 2 .
Furthermore, the dimensionless diffusivity coefficients are α x = 1.97 • 10 -3 and α y = 7.91 • 10 -5 . They have distinct values because the length of each domain is not the same.

that the problem is linear. The difference from the validation case is that two peaks with low frequency and strong energy are present. x and y and not only in one direction. If the ambient vapour pressure P v, ∞ were not dependent on the coordinated y, the solution would be a false two-dimensional, with variations only in the x coordinate. The variation of the boundary condition regarding y can be also considered in the convective transfer coefficient h m (y) and on an incoming liquid flow g l (y) .

The final solution provided by each model is a discrete tensor of dimension 101 × 101 × 9601 . To view the values of this tensor in a global way, Figure 3.50 presents several slices. In this figure, we have three slices of y and four slices of x, for all values of t . This complete tensor has a size of 747.22 MB and it only contains information for a simulation performed for a small period of time. For this reason, the tensor decomposition becomes more important when the problem dimension grows. By decomposing it with a 12-rank CP-approximation, the tensor is reduced to 0.9 MB . The decompositions took 2.2 s to be run and if one needs to recompose the original tensor it takes approximatively the same amount of time. 101

Chapter 4

Conclusion

Most of the numerical methods applied to mathematical models used in building physics are commonly based on classical methods to compute the solution of parabolic problems, such as finite differences and finite volumes approaches. The time discretization is chosen according to appropriate stability conditions. Explicit schemes are not used because they have limitations regarding the CFL stability condition, which leads to prohibitively small time steps that require a high overall computational run time. Thus, implicit schemes as the Crank-Nicolson and Euler Implicit schemes are routinely used. However, even with large time steps, the nonlinearities of the problem demand sub-iterations at each time step that are also computationally costly. For this reason, the following question has been raised: how to simulate efficiently, i.e., accurately and rapidly, nonlinear problems of heat and moisture transfer in buildings?

To answer that question, the manuscript presented a synthesis of the application of innovative numerical methods for moisture diffusion in porous material. The purpose is to highlight the strength of each method on simple but realistic cases. More complex cases of coupled heat and mass transfers are considered in the articles in Appendix C.

Chapter 2 details the methodology to build innovative numerical models with the DuFort-Frankel method, the Quasi-Uniform Nonlinear Transformation, the Method of Horizontal Lines and the reduced Spectral method. The first three methods are considered as improved complete original models because they need to solve large systems of equations but with a reduced cost compared to the classical approaches used in the literature. The spectral method is considered as a reduced order model, since the size of the system to be solved is significantly reduced when compared to the original models. Chapter 3 presents the features of the methods on a real case application. First, they were compared in a simple case of one-dimensional nonlinear moisture transfer so the important features can be highlighted. To summarize, Table 4.1 presents the numerical features of each method so one can compare the methods presented in this work. Regarding the grid mesh, the solution and the precision in space, the precision in time and the possibility to expand to higher dimensions, the reduced Spectral method highlights from the others. As spectral method presented the best features, it is evaluated in a parametric problem and in two-dimensional transfer problems, which are the case studies of the sequence. Therefore, the main contribution of this work is the detailed and systematic study of innovative methods to treat the complexities of nonlinear diffusion transport phenomena in building physics. The efficiency of the numerical methods depends on different requirements, which difficult the choice of the best method, as each one has their strengths and disadvantages. To illustrate the different applicabilities, Figure 4.1 presents some of the main criteria that can be used to evaluate the numerical methods: (i) spatial complexities, (ii) time scales and (iii) reduction of the computational cost. For example, if one aims at dealing with spatial complexities and accurate modelling of the phenomena at interfaces between two materials, the MOHL and QUNT are good candidates (Gasparin et al., 2018a(Gasparin et al., , 2019b)). For coupling the model of heat and moisture transfer with others through co-simulation approaches [START_REF] Wetter | Co-simulation of building energy and control systems with the building controls virtual test bed[END_REF], the DuFort-Frankel is a very promising approach. It is explicit and avoids extra sub-iterations due to model coupling [START_REF] Gasparin | An improved explicit scheme for whole-building hygrothermal simulation[END_REF]. For strong reduction of the computational cost with a very satisfactory accuracy of the solution, the Spectral has been shown as the best option (Gasparin et al., 2018e). It is particularly efficient to treat parametric problems where one aims to compute the solution depending on an input parameter [START_REF] Gasparin | Advanced reduced-order models for moisture diffusion in porous media[END_REF]. Furthermore, for 2D transport phenomena the Spectral method is also an interesting ap-proach. Not to mention that the DuFort-Frankel and the QUNT approaches can also be extended to two-dimensional problems.

Perspectives

For further work, many aspects can be considered in order to apply the presented methods in the building simulation field:

• Complexity of the physical phenomena: In this work, the diffusion transfer was the main physical phenomena studied. From a physical point of view, a lot of phenomena can be added to improve the mathematical model. For example, the hysteresis is an important phenomenon that can influence in the final solution but was not explored [START_REF] Kwiatkowski | Modelling of hysteresis influence on mass transfer in building materials[END_REF][START_REF] Pedersen | Prediction of moisture transfer in building constructions[END_REF]. We have considered the coupling between heat and moisture transfer in a week way, in which we neglected the heat gradient in the moisture balance equation. At the interface between two materials, we have considered continuous field, flux, and flow. However, when other driving potentials are used, we may have resistance at the interface, which generates a discontinuity of the solution.

• Multi-dimensional simulation: In this study, the one-dimensional heat and moisture transfer is widely explored with different numerical methods. The results provide important perspectives for integration in building simulation programs such as Domus, EnergyPlus, etc. However, our study on two-dimensional transfer problems is introductory. To go further in this direction, we need to investigate deeply the nonlinear transfer in two-and three-dimensional geometries, with multilayered materials. In addition, simulate the moisture transfer in geometries like corners and foundations are of fundamental importance since the moisture accumulation in these parts can bring serious durability problems. In addition, it would be interesting to analyze the difference between the computation of the heat balance in zones considering onedimensional transfers against the heat balance in zones considering two-dimensional transfers.

• Validation: We have compared our solutions of the numerical methods against experimental data in only one case, the one-dimensional heat and moisture transfer through one material. Such efficient numerical approaches can be used for comparison with experimental data for more complex model. Indeed, the MOHL could be used to compare the numerical prediction with experimental data at the interface between two porous material. It would require global measurement using for instance neutron beam [START_REF] Sedighi-Gilani | Visualization and quantification of liquid water transport in softwood by means of neutron radiography[END_REF]. Moreover, a experimental benchmark two-dimensional heat and mass transfer could be carried out using the Spectral ROM as a numerical model.

• Large spatial domains: As the spectral method is a type of global method, it can compute the solution for all the spatial domain independently on the number of spatial nodes. This feature makes this method very suitable for problems with a large spatial domain, such as in ground simulation, since we do not need a refined grid.

• Coupling with simulation software: All the methods and case studies have been implemented in the Matlab environmental. We would like to implement the DuFort-Frankel and the reduced Spectral methods in a simulation tool so we can evaluate the gains of the direct implementation. Another alternative would be to perform a co-simulation alongside with the software, in which the program and the numerical method are simulated separately but exchanging information at each time step [START_REF] Mazuroski | An artificial intelligencebased method to efficiently bring cfd to building simulation[END_REF].

• Inverse problems: The reduced spectral method can be a good option to perform inverse problems. When one needs to estimate a parameter, the direct model of heat and moisture transfer need to be simulated thousands of times which are computationally expensive. With the reduced order model, the estimation can be performed much faster [START_REF] Berger | Proper generalized decomposition model reduction in the bayesian framework for solving inverse heat transfer problems[END_REF].

• Optimisation of the algorithms: The numerical methods are implemented in the algorithms in the best form we could. However, they are not optimized. Currently, in the Spectral method algorithm, the number of modes is determined before starting simulations. It would be interesting if the algorithm could do it automatically.

Besides, when we simulate the heat and moisture transfer, we compute first the equation for the moisture transfer and then we compute the equation for the heat transfer. We would like to perform both equations at the same time.

This list of perspectives opens many fields of investigation and it encourages to go further on the research. so that the orthogonality property of the Chebyshev polynomials is written as:

T i , T j = 1 -1 T i (x) • T j (x) √ 1 -x 2 dx =          0 , if i = j , π , if i = j = 0 , π 2 , if i = j = 0 .
The integrals can be approximated by a finite sum using the Chebyshev-Gauss quadrature:

1 -1 f (x) • g (x) √ 1 -x 2 dx ≈ π m • m k = 1 f (x k ) g (x k ) ,
with x k being the Gauss points.

The spatial derivatives are re-expanded in the same Chebyshev basis function:

∂u n ∂ x = n i = 0 a i (t) ∂T i ∂ x (x) = n i = 0 ã i (t) T i (x) , (A.1a) ∂ 2 u n ∂ x 2 = n i = 0 a i (t) ∂ 2 T i ∂ x 2 (x) = n i = 0 ã i (t) T i (x) , (A.1b)
As a result, coefficients {ã i (t)} n i=0 and { ã i (t)} n i=0 must be re-expressed in terms of coefficients {a i (t)} n i=0 . The connection is given explicitly from the recurrence relation of the Chebyshev polynomial derivatives: 121

Figure 2 . 1 :

 21 Figure 2.1: Scheme of the numerical methods presented in this thesis.

Figure 2 . 2 :

 22 Figure 2.2: Cartesian representation of the space-time discretization.

  Figure 2.3(b) gives the Cartesian representation of semi-discretization in time scheme.

Figure 2 . 3 :

 23 Figure 2.3: Cartesian representation of the semi-descritization in space (a) and in time (b).

Figure 2 . 4 :

 24 Figure 2.4: Representation of a cell with the location of the points.

Figure 2 . 5 :

 25 Figure 2.5: Stencil of the DuFort-Frankel numerical scheme.

Figure 2 . 6 :

 26 Figure 2.6: The classical numerical solution strategy based on the MOL (a) and the proposed approach view based on BVPs (b).

Figure 2 . 7 :

 27 Figure 2.7: Representation of the parametric solution.

Figure 2 . 8 :

 28 Figure 2.8: Representation of the two-dimensional solution collocation points.

  Model . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Description of the case study . . . . . . . . . . . . . . . . . . . . 3.1.3 Assessment performance for the numerical solution . . . . . . . . 3.1.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Numerical experiment 2: a parametric problem . . . . . . . . 3.2.1 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Description of the case study . . . . . . . . . . . . . . . . . . . . 3.2.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4 Tensor representation of the solution . . . . . . . . . . . . . . . . 3.3 Numerical experiment 3: 2D transport . . . . . . . . . . . . . 3.3.1 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Assessment performance for the numerical solution . . . . . . . . 3.3.3 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Description of the case study . . . . . . . . . . . . . . . . . . . . 3.3.5 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: One-dimensional physical wall model.

  t) [Pa]. It comes from the relation: φ : (P v , P s ) -→ P v P s where φ [ -] is the relative humidity expressed without physical dimension, P s : T -→ P s (T ) [Pa] is the saturation pressure and T [K] is the temperature. The coefficient k m [s]

Figure 3 . 2 :

 32 Figure 3.2: Sorption isotherm (a) and vapour permeability (b) for a temperature of 23 • C.

Figure 3 . 3 :

 33 Figure 3.3: Moisture storage coefficient (a) and total moisture transfer coefficient (b).

Figure 3 . 4 :

 34 Figure 3.4: Ambient relative humidity φ ∞ .

Figure 3 . 5 :

 35 Figure 3.5: Vapour pressure profiles at t = 24 h, t = 48 h and t = 84 h (a) and its evolution at the edges of the material (b).

Figure 3 . 7 :

 37 Figure 3.7: Error ε ∞ computed between the reference solution and the other numerical schemes as a function of ∆t (a) and the respective ratio of the computational time R cpu to perform each simulation with relation to t • = 449.3 s (the maximal computational time t cpu of the Implicit simulation) (b).

Figure 3 . 8 :

 38 Figure 3.8: Average number of sub-iterations of the Implicit scheme as a function of time discretization ∆t .

Figure 3 . 9 :

 39 Figure 3.9: Error ε ∞ computed for different ∆x as a function of ∆t (a) and as a function of τ /∆t (b).

Figure 3 Figure 3 . 11 :

 3311 Figure 3.10: Trajectories of the grid nodes in space-time.

Figure 3 . 12 :

 312 Figure 3.12: Error ε 2 as a function of the spatial domain x .

Figure 3 Figure 3 . 14 :

 3314 Figure 3.13: Global error ε ∞ computed in function of the number of spatial nodes, for ∆t = 10 -3 (a) and its respective CPU time necessary to perform the simulation (b).

Figure 3

 3 Figure 3.15: Global error ε ∞ computed in function of the number of spatial nodes, for ∆t = 10 -1 (a) and its respective CPU time necessary to perform the simulation (b).

Figure 3

 3 Figure 3.16: Global error of the flow ξ ∞ computed as a function of the number of spatial nodes.

Figure 3 . 17 :

 317 Figure3.17: Three vapour pressure profiles with the selected mesh determined by theMOHL (a) and the evolution of the vapour pressure on the edges of the material (b).

Figure 3 . 18 :

 318 Figure 3.18: Number of spatial grid points over time (a) and the respected selected grid (b).

Figure 3 . 20 :

 320 Figure 3.20: Moisture flow at the boundaries computed wit the derivative provided by the MOHL.

Figure 3

 3 Figure 3.21: Global error ε ∞ computed as a function of the time discretization.

Figure 3 Figure 3 . 23 :

 3323 Figure 3.22: Vapour pressure profiles at different times instants, at t = {24; 48; 84} h (a) and its evolution on the edges of the material (b).

Figure 3 . 24 :

 324 Figure 3.24: Evolution of the first three (a) and of the last three (b) spectral coefficients a n .

Figure 3

 3 Figure 3.25: Error ε ∞ as a function of the number of spatial nodes N x (a) and the computational time spent in each simulation (b).

Figure 3 . 26 :

 326 Figure 3.26: Error ε 2 computed for different numbers of modes N .

Figure 3 . 27 :

 327 Figure 3.27: Maximum absolute error as a function of the number of spectral modes (a) and their respective CPU simulation time (b).

Figure 3 Figure 3 . 29 :

 3329 Figure 3.28: Fourier power spectrum of the Spectral solution computed at the boundaries and at the middle of the material x = { 0; 0.075; 0.15 } .

  α 6 = 4.20 • 10 -10

Figure 3 . 30 :

 330 Figure 3.30: Evolution of the relative humidity at the left boundary on the parametric case.

Figure 3

 3 Figure 3.31: Vapour pressure profiles for each one of the materials at the time instant of t = 48 h (a) and the vapour pressure evolution near the left boundary at x = 1 cm (b).

Figure 3

 3 Figure 3.32: Evolution of the mass variation of each material.

Figure 3

 3 Figure 3.33: Error ε 2 computed for the Spectral-parametric solution (a) and also computed for the Implicit solution.

Figure

  Figure3.34(a) presents the last spectral coefficients of the solutions. The magnitude of the spectral coefficients is higher than the error of the solution O(10 -6 ) for α max and

Figure 3 . 34 :

 334 Figure 3.34: Last spectral coefficients (a) and the Fourier Power Spectrum at the left boundary (b) for each material.

Figure 3

 3 Figure 3.35: Convergence of the Spectral modes M (a) and its computational run time (b) for different quantities of parameters α k .

Figure 3

 3 Figure 3.36: Computational time t cpu as a function of the number of elements N α .

Figure 3 .

 3 Figure 3.37: Tensor representation of the solution P v (x, t, α).

Figure 3 . 38 :

 338 Figure 3.38: Error ε 2 (a), computational time t cpu (b) and memory size of the compressed tensor (c) as function of the number of rank R .

  .7) with the moisture capacity coefficient c m [kg/m 3 Pa] and the moisture diffusivity coefficient k m [s] are considered as positive constants:

Figure 3 . 39 :

 339 Figure 3.39: Two dimensional wall model.

Figure 3 . 40 :

 340 Figure 3.40: Initial condition u 0 (x , y ) of the validation case.

Figure 3

 3 Figure 3.41: Evolution of the dimensionless moisture field u at the middle of the material (a) and the error of each numerical solution (b).

Figure 3 . 42 :

 342 Figure 3.42: Error 2 among the x and y axes.

Figure 3 . 43 :

 343 Figure 3.43: Evolution of the spectral coefficients a i,m (a) and the Fourier Power Spectrum for (x , 0.5) (b).

Figure 3

 3 Figure 3.44: Error ∞ as function of ∆t .

Figure 3

 3 Figure 3.45: Error ∞ as a function of the total number of modes N • M .

Figure 3 . 46 :

 346 Figure 3.46: Evolution of the relative humidity at the left (a) and right boundaries (b).

Figure 3 . 48 :

 348 Figure 3.48: Evolution of the spectral coefficients {a i,m } n i=0 (a) and the Fourier Power Spectrum for (x, y = 0.12) (b).

Figure 3 .

 3 Figure 3.49 presents the vapour pressure as function of x and y , at four different time instants. Each one of the graphics represents a slice of the discrete tensor solution. It can be observed that the solutions of the vapour pressure vary according to both directions

Figure 3 . 49 :

 349 Figure 3.49: Slices of the solution for different time instants t = {32, 50, 80, 96} h.

Figure 3 .

 3 Figure 3.50: Tensor representation of solution the solution P v (x, y, t).

Figure 4 . 1 :

 41 Figure 4.1: Synthesis of the numerical methods advantages.

A

  = a 0 , a 1 , • • • , a n-1 , a n , Ã = ã 0 , ã 1 , • • • , ã n-1 , ã n , Ã = ã 0 , ã 1 , • • • , ã n-1 , ã n ,and the elements of matrices D and D are deduced from Equations (2.23) and (2.24).Moreover, the temporal derivative is also written with the same Chebyshev basis function:

  

  , Stable explicit schemes for simulation of nonlinear moisture transfer in porous materials, J. Buil. Perf. Sim., Vol 11(2), pp 129-144, Doi https:

	org/10.1007/s11242-018-1106-2.
	7. Gasparin et al. (2019), A spectral method for solving heat and moisture transfer
	through consolidated porous media. Int. J. Numer. Methods Eng., Vol 117(11), pp
	1143-1170, Doi https://doi.org/10.1002/nme.5994.

//doi.org/10.

1080/19401493.2017.1298669. 

2.

Gasparin et al. (2018)

, An improved explicit scheme for whole-building hygrothermal simulation, Buil. Sim., Vol 11(3), pp 465-481, Doi https://doi.org/10.1007/ s12273-017-0419-3. 3. Gasparin et al. (2018), An adaptive simulation of nonlinear heat and moisture transfer as a boundary value problem, Int. J. Ther. Sci., Vol 133, pp 120-139, Doi https://doi.org/10.1016/j.ijthermalsci.2018.07.013. 4. Gasparin et al. (2019), An innovative method to determine optimum insulation thickness based on non-uniform adaptive moving grid J. Braz. Soc. Mech. Sci. Eng., Vol 41(4), pp 173, Doi https://doi.org/10.1007/s40430-019-1670-6. 5. Gasparin et al. (2019), Solving nonlinear diffusion problems in buildings by means of a Spectral reduced-order model, J. Buil. Perf. Sim., Vol 12(1), pp 17-36, Doi https://doi.org/10.1080/19401493.2018.1458905. 6. Gasparin et al. (2018), Advanced Reduced-Order Models for Moisture Diffusion in Porous Media, Transp. Porous Med., Vol 124(3), pp 965-994, Doi https://doi.

  Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 The mathematical problem . . . . . . . . . . . . . . . . . . . . . 2.1.3 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Chapter 2
	Numerical Methods

Contents 2.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 2.2 The DuFort-Frankel scheme . . . . . . . . . . . . . . . . . . . . 2.2.1 Important Features . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Extended Application . . . . . . . . . . . . . . . . . . . . . . . . 2.3 The Method of Horizontal Lines . . . . . . . . . . . . . . . . . 2.3.1 Important Features . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Extended Application . . . . . . . . . . . . . . . . . . . . . . . . 2.4 The Quasi-Uniform Nonlinear Transformation . . . . . . . . . 2.4.1 Important Features . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Extended Application . . . . . . . . . . . . . . . . . . . . . . . . 2.5 The Reduced Spectral method . . . . . . . . . . . . . . . . . . 2.5.1 One-dimensional space . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Parametric approach . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.3 Two-dimensional space . . . . . . . . . . . . . . . . . . . . . . . . 2.5.4 Important Features . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.5 Extended Application . . . . . . . . . . . . . . . . . . . . . . . .

Table 3

 3 Figure 3.6: Error ε 2 as a function of the spatial domain x .

	10 -2			
	10 -3			
	10 -4			
	0		0.05	0.1	0.15
	Reference	Implicit	DuFort-Frankel
	ε ∞	-	9.15 • 10 -4	1.35 • 10 -3
	t cpu	43.94 s	28.05 s	1.37 s
	R cpu	160 %	100 %	5 %

.3: Error ε ∞ and computational time of the solutions simulated with the same discretization values ∆t = 10 -2 and ∆x = 10 -2 .

Table 3

 3 

	98 • 10 -4

.4: Computational time and accuracy of the solutions simulated with the same discretization values ∆t = 10 -2 and ∆x = 4 • 10 -2 .

Table 3

 3 Figure 3.19: Error ε 2 computed over the time with the MOHL and with the Implicit approach.

	10 -2			
	10 -3			
	10 -4			
	0	0.05	0.1	0.15
	tol	ε ∞	Average N x	t cpu
	10 -6 7.22 • 10 -5	722	586.6 s
	10 -5 1.44 • 10 -4	396	228.8 s
	10 -4 2.33 • 10 -3	111	89.43 s
	10 -3 4.97 • 10 -3	61	57.35 s
	10 -2 2.38 • 10 -2	25	36.6 s
	10 -1 3.47 • 10 -2	20	32.2 s

.5: Maximum error ε ∞ computed for different tolerances with the solver bvp4c, considering ∆t = 10 -1 .

  Numerical Scheme CPU time (s) CPU time (%) Average iterations Error ε ∞

	Reference Chebfun	44.1	139	-	-
	Implicit	31.7	100	8	9.27 • 10 -4

  Table 3.8: Computational time of the solutions with the same order of accuracy O 3 • 10 -3 .

	min]	6	6	1.8	12	12	12
	∆x	10 -2	10 -2	2 • 10 -2	-	-	10 -2
	∆x [cm]	0.15	0.15	0.3	-	-	0.15
	N x	101	101	51	41	[ 20, 58 ]	101
	tol	-	-	-	-	10 -3	2.7 • 10 -3
	ε ∞	-	3.13 • 10 -3	3.22 • 10 -3	3.39 • 10 -3 3.78 • 10 -3	3.11 • 10 -3
	t cpu [ s ]	9.4	5.3	0.3	0.3	27	0.3
	R cpu	1.8	1	0.06	0.06	5.1	0.06

  Extrema3.79 • 10 -3 6.12 • 10 -2 7.22 • 10 -2 3.47• 10 -2 1.22 • 10 -1 1.29 • -5 Equi-spaced 3.79 • 10 -3 7.27 • 10 -2 8.79 • 10 -2 5.55 • 10 -2 1.13 • 10 -1 1.35 • -5Table3.10: Error ε ∞ for the solution of each parameter for different collocation points.

	Collocation points			ε ∞			
	α k	α 1	α 2	α 3	α 4	α 5	α 6
	Itself	3.79 • 10 -3 4.53 • 10 -4 2.34 • 10 -4 1.95 • 10 -6 1.67 • 10 -6 1.32 • -6

Table 3 .

 3 11: Computational time required to compute the parametric problem and the error ε ∞ of solution of α min and α max .

	Method	ε ∞ of α min ε ∞ of α max	t cpu	R cpu
	Reference	-	-	132.35 s 5.27
	Implicit	2.71 • 10 -3 2.67 • 10 -5	25.10 s	1
	Spectral-Parametric 3.79 • 10 -3 1.33 • 10 -6	5.29 s	0.21
	Spectral	3.79 • 10 -3 1.11 • 10 -5	2.06 s	0.08

Table 3 .

 3 12: Features of the original and of the decomposed tensors.

Table 3 .

 3 13: Computational time of the solutions with the same order of accuracy O 10 -6 .

  Solution in space Precision in space Grid mesh Precision in time DimensionTable 4.1: Global comparison of the numerical methods.

	DF	local	medium	fine	low	n
	QUNT	local	medium	adaptive	medium	n
	MOHL	global	high	adaptive	medium	1-2
	Spectral	global	high	coarse	high	n

Acknowledgements

First, I would like to thank the external members of the board, Prof. Jean-Jacques Roux, Prof. Aristeu da Silveira Neto, Prof. Carsten Rode and Prof. Monika Woloszyn,

The financial support came from the Brazilian agency CAPES of the Ministry of Education and from the COFECUB project. It consists of an interdisciplinary work involving the fields of Mechanical Engineering and Applied Mathematics. Scholarships 2016 Scholarship CAPES-Cofecub/ PhD student exchange to France. 2010 Marcelino Champagnat Prize -best student of the class/ Free scholarship for a master's degree at Pontifical Catholic University of Paraná (PUCPR).

Abbreviation

ADI

Alternating Direction Implicit method

xvi Chapter 1

Results and discussion

The Spectral-Collocation method is simulated considering 9 modes for both basis functions, namely N • M = 9 • 9 = 81 total modes. The tolerance of the solver ode15s is set to tol = 10 -4 which integrates the spectral coefficients in a temporal grid of ∆t = 10 -2 . The final solution is composed for a spatial grid consisting of N x = N y = 101 nodes. For the ADI method, the discretization parameters are defined as ∆t = 10 -2 and ∆x = ∆y = 10 -2 , the same values of the solution computed for the Spectral approach. As expected, the Spectral method is able to compute the solution in 16.8 s, which is almost twice faster than the ADI method, which requires about 25.9 s.

Figure 3.47(a) presents different profiles (also called the fibers of the tensor) of vapour pressure as a function of the space x . The profiles are located at y = 0.24 m and they are presented for the instants t = {32, 50, 80, 96} . On the other hand, Figure 3.47(b) presents the evolution of vapour pressure at the locations (x, y) = (0, 0); (0, 0.25);

(0, 0.5) m . Both ADI and Spectral-Collocation methods have converged to the same solution, apparently. The error cannot be calculated, since there is no reference solution. For the Spectral method, one of the indicative of error is the last spectral coefficient. In this simulation, the value of the last coefficients {a i,m } n i=0 are given in Figure 3.48(a), which oscillates at an order of O(10 -9 ) . As the results match the ADI solution, it is assumed that the number of modes is enough to provide a satisfactory solution. In Figure 3.48(b), one can find the Fourier power spectrum of the solution. The information given by this graphic is that there are two strong signals, that come from the boundary conditions and

Appendix A

Chebyshev polynomials

Here one can find some important information about the Chebyshev polynomials used in this work, which are based mainly in [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF]; [START_REF] Canuto | Spectral Methods -Fundamentals in Single Domains[END_REF]; [START_REF] Fornberg | A pratical guide to pseudospectral methods[END_REF]; [START_REF] Peyret | Spectral methods for incompressible viscous flow[END_REF]; [START_REF] Trefethen | Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations[END_REF].

The Chebyshev polynomials are a family of orthogonal polynomials that are a practical basis for approximating smooth functions on a segment. The Chebyshev polynomials of first kind are denoted by T n (x) , which is a polynomial of degree n on the segment x ∈ [-1, 1] that is written as:

For instance, the first few polynomials are:

Higher order polynomials can be constructed using the following recursive relation :

The values of the Chebyshev polynomials at the boundaries, at x = ±1, are given by:

In addition, the roots of the polynomial T n (x) = 0 are given by the Gauß points:

The scalar product f, g with respect to the weight 1 -x 2 -1/2 is defined as:

Rescaling

To apply the orthogonal polynomials used in this work, some domains must be changed to fit the segment [-1, 1].

B.1 Change of variables on the spatial domain

For this reason, a change in the spatial domain segment x : [ 0, ] -→ [ -1, 1 ] is carried out on the following system:

First, we consider a linear relation between x and s such that:

where a and b are real constants. To determine their values, we relate the extrema points:

which leads to:

Then, the function x is written as:

Substituting Equation (B.2) in the system of Equation (B.1) leads to the system: Consider a linear relation between α and α such that:

where a and b are real constants. To determine their values, we relate the extrema points:

Then, the function α is written as:

and, the function α is written as:

Substituting Equation (B.5) in the system of Equation (B.3) leads to the problem:

∂ 2 u ∂s 2 , s, α ∈ ( -1, 1 ) and t > , u = f (s) , s ∈ ( -1, 1 ) and t = ,