Julien Berger
email: julien.berger@univ-lr.fr

Benjamin Kadoch

Estimation of the thermal properties of an historic building wall by combining Modal Identification Method and Optimal Experiment Design

Keywords: Parameter estimation problem, Model Calibration, Inverse Heat Conduction Problem, Modal Identification Method, Reduced Order Model

The estimation of wall thermal properties by in situ measurement enables to increase the reliability of the model predictions for building energy efficiency. Nevertheless, retrieving the unknown parameters has an important computational cost. Indeed, several computations of the heat transfer problem are required to identify these thermal properties. To handle this drawback, an innovative approach is investigated. The first step is to search the optimal experiment design among the sequence of observation of several months. A reduced sequence of observations of three days is identified which guarantees to estimate the parameter with the maximum accuracy. Moreover, the inverse problem is only solved for this short sequence. To decrease further the computational efforts, a reduced order model based on the modal identification method is employed. This a posteriori model reduction method approximates the solution with a lower degree of freedom. The whole methodology is illustrated to estimate the thermal diffusivity of an historical building that has been monitored with temperature sensors for several months. The computational efforts is cut by five. The estimated parameter improves the reliability of the predictions of the wall thermal efficiency.

Introduction

Within the environmental context, requirements on building energy efficiency becomes more and more important. In France, the building stock increases slowly [START_REF] Stephan | Experimental assessment of thermal inertia in insulated and non-insulated old limestone buildings[END_REF] with an average around 1% . Thus, there is a crucial issue in retrofitting existing building. To this end, models such as building simulation programs are employed to define the best actions. However, there are some uncertainties on the input parameters of the model [2,3].Particularly, for such old buildings, the material properties may differ strongly from one site to another [4]. Indeed, due to the vernacular traditional architecture, the material were generally taken from locals around the building. Furthermore, the aging of materials may have changed the properties [5]. Parameter such as the thermal properties of the walls require to be known precisely since they play a crucial role on the assessment of the building energy efficiency [START_REF] Jumabekova | Sensitivity analysis in the framework of parameter estimation problem for building energy performance: a continuous derivative based approach[END_REF]7].

The unknown properties can be estimated by solving inverse problems or more precisely parameter estimation problem [8] to calibrate the model. This procedure aims at minimizing a so-called cost function between the model numerical predictions and the experimental observations obtained by on-site measurements.Such approaches have been successfully applied for real walls and a detailed state-of-the-art is proposed in [9]. However, the determination of thermophysical properties using in situ dynamic measurements faces important challenge in terms of computational effort. This challenge arises from two reasons. The first one is the duration of the experimental observations. If long measurement duration are better for the accuracy of the retrieved parameter, it implies important computational costs to solve the inverse problem. Note that in situ measurement depends also on occupants' acceptance and long period can be constraining. The standard techniques ISO 9869 requires more than a month of measurement periods [START_REF] Rasooli | In-situ characterization of walls' thermal resistance: An extension to the iso 9869 standard method[END_REF]. Thus, various alternative techniques have been developed to reduced the monitoring time. In [11] an experiment of 7 days is carried. In [START_REF] Cornaro | Dynamic simulation and on-site measurements for energy retrofit of complex historic buildings: Villa mondragone case study[END_REF], the accuracy of the parameter estimation is investigated according to the measurement period. For their case study, the authors suggest to carry at least three days of monitoring to have an accurate estimation. In [9], a fast estimation method is proposed based excitation pulse method and less than 1 h of measurement.

The second reason is the mathematical model to represent the physical phenomena. Indeed, the computational efforts required to solve the inverse problem are directly proportional to the so-called direct model (as well as the optimization algorithm). Thus, numbers of investigations aimed at developing direct model with reduced computational cost. Lumped model such as RC approaches have been proposed due to their very small computational cost. However, their reliability to estimate accurately the thermal diffusivity is questioned in [13]. The model has to be based on detailed heat diffusion process. In [14] an equivalent homogeneous model is proposed for the description of multi-layer walls with the issue of reducing the cost of parameter estimation as shown in [5]. When possible, analytical solution are also developed as for instance in [15].

Thus, the objective of this article is to answer both issues in order to propose a fast and accurate parameter estimation method. First, the length of the measurement observations is determined with the optimal experiment design methodology [16,17]. The optimal design is searched according to the conditions of the experiment. Carrying the measurement for the optimal experiment design ensure to estimate the unknown parameter with the highest accuracy. This approach has been used in [18,[START_REF] Rasooli | In-situ rapid determination of walls' thermal conductivity, volumetric heat capacity, and thermal resistance, using response factors[END_REF] to determine the optimal boundary conditions of laboratory experiments to retrieve material properties of heat and mass transfer. In [START_REF] De Monte | Optimal experiment design for thermal property estimation using a boundary condition of the fourth kind with a time-limited heating period[END_REF], an optimal design is found to retrieve the thermal conductivity and volumetric heat capacity of a material using laboratory experiments. This approach can be extended for in-situ measurement to determine the optimum sequence of measurement observations. Then, to decrease the computational effort of the direct numerical model, model reduction techniques can be employed. Several methods are reported in the literature to model the physical phenomena in building walls [START_REF] Berger | Review of Reduced Order Models for Heat and Moisture Transfer in Building Physics with Emphasis in PGD Approaches[END_REF]. Among all, the Modal Identification Method (MIM) demonstrated successful applications for inverse problem. The primer works of Girault, Petit, and Videcoq introduced the MIM approach for direct simulations in [START_REF] Petit | A Modal Identification Method To Reduce A High-Order Model: Application To Heat Conduction Modelling[END_REF]23]. Then, it has been applied for the identification of boundary fluxes through sequential procedure in [24][25][START_REF] Girault | Identification methods in nonlinear heat conduction. part ii: inverse problem using a reduced model[END_REF] or surface temperature as in [START_REF] Gérardin | Réduction de modèle thermique par Méthode d'Identification Modale (MIM) pour déterminer la température de surface des composants de machine de fusion[END_REF]. To our best knowledge, it has never been employed for the estimation of thermal diffusivity in building walls under climatic varying conditions.The MIM ROM is employed since in the state-space representation of the MIM model, both the field and its sensitivity to the unknown parameter are computed with a reduced computational cost. It can be combined with a gradient-based algorithm to retrieve the thermal properties of the wall.

This article presents and evaluate the proposed methodology to estimate the thermal diffusivity of a wall monitored in an old building. First the mathematical model is described in Section 2. The MIM model reduction method is detailed in Section 3. Then, the search optimal experiment design combined with the solution of the parameter estimation problem is presented in Section 4. A brief synthesis of the whole methodology and the metrics to evaluate its efficiency are given in Section 5. Subsequently, a first case study is introduced to validate the MIM method for estimating the thermal diffusivity of a wall. Last, a real case study is investigated in Section 7.

Mathematical model

First, the mathematical model is presented. It is based on two equations. The first corresponds to the diffusion one occurring in the wall. In addition, the sensitivity equation is also described. It enables to compute the sensitivity coefficient which qualifies the sensitivity of the field of temperature relative to change in the unknown diffusivity. This is required in the algorithm to solve the inverse problem.

Physical formulation

The field of interest is the temperature T evolving in a building wall materialaccording to a diffusion process. The space domain is defined by x ∈ 0 L , where L is the length of the wall. The time domain is defined by t ∈ Ω t with Ω t def := 0 , t f . The temperature is computed using the heat diffusion equation [START_REF] Ozisik | Inverse Heat Transfer: Fundamentals and Applications[END_REF]:

∂T ∂t = α • ∂ 2 T ∂x 2 , ( 1
)
where

α def := k c
is the heat diffusion, k is the heat conductivity and c is the volumetric heat capacity. At the interface between ambient air and the wall, Dirichlet type boundary conditions are assumed.The temperature is prescribed by the one measured:

T = T ∞ , L t , x = 0 , T = T ∞ , R t , x = L ,
where T ∞ are the temperatures measured by sensors and depending on time:

T ∞ : t -→ T ∞ t .
At t = 0 , the temperature in the material is known:

T = T ini x , t = 0 .
with T ini a function of x corresponding to the initial field in the wall:

T ini : x -→ T ini x .
The so-called observable field T s is defined as follows:

T s : ( x s , t ) -→ T x s , t , s ∈ 1 , . . . , N s ,
where N s is the number of points of interests. In our case, it corresponds to the sensor locations.Thus, the observable fields correspond to the temperature computed with the model at the sensor locations. It enables to compare the model predictions with the experimental observations. The direct problem is now transformed into a dimensionless formulation to be solved by the numerical model.

Dimensionless formulation

As discussed and thoroughly motivated in [START_REF] Nayfeh | Perturbation Methods[END_REF]30], it is of capital importance to obtain a dimensionless problem before elaborating a numerical model.It allows to define scaling parameters such as Fourier number that may decrease the number of unknown parameters in the inverse problem. For this, dimensionless fields are defined:

u def := T -T min T ref , u ∞ , L def := T ∞ , L -T min T ref , u ∞ , R def := T ∞ , R -T min T ref , u ini def := T ini -T min T ref ,
where T ref and T min are user-defined reference temperature. The space and time coordinates are also transformed into dimensionless variables:

t def := t t ref , x def := x L .
The Fourier dimensionless number is defined:

Fo def := t ref • α L 2 .
It corresponds to the thermal diffusivity in the dimensionless representation and quantifies the magnitude of diffusion inside the material. A high number indicates a fast heat transfer process inside the material. With these transformations, the dimensionless problem is written as:

∂u ∂t = Fo • ∂ 2 u ∂ x 2 (2)
with the boundary condition:

u = u ∞ , L t , x = 0 (3a) u = u ∞ , R t , x = 1 (3b)
and the initial condition:

u = u ini x , t = 0 . ( 4
)
To have a well-posed problem, initial and boundary conditions must be compatible. The observable field is also transformed according to u s def

:= T s -T min T ref , s ∈ 1 , . . . , N s .

Sensitivity equation

The issue is to estimate the unknown thermal diffusivity α of the wall or its equivalent in the dimensionless representation, the Fourier number Fo . The solution of the inverse problem requires the so-called sensitivity function θ :

θ : ( Fo , x , t ) -→ ∂u ∂Fo .
It quantifies the sensitivity of the field of temperature according to the unknown parameter. A small magnitude reveals that large change in the parameter induces small change in the field. Thus, the parameter cannot be estimated with accuracy. Inversely, high magnitude of sensitivity coefficient are favorable conditions to estimate the unknown parameter. One can obtain can compute the sensitivity equations by differentiating Eq. (2) relatively to the parameter:

∂θ ∂t = Fo • ∂ 2 θ ∂x 2 + ∂ 2 u ∂x 2 , ( 5
)
with the following boundary conditions obtained by differentiating Eq. (3):

θ = 0 , x = 0 , θ = 0 , x = 1 ,
and the initial condition arising from Eq. (4):

θ = 0 , t = 0 .
Particularly, the sensitivity function at the location x s of the sensors are needed. It is computed by:

θ s : ( x s , t ) -→ θ x s , t , s ∈ 1 , . . . , N s .

MIM reduced order model of the field and its sensitivity

In this Section, the numerical models used to perform simulations are described starting first with the large original model. Then the reduced MIM approach is detailed. A uniform discretisation based on finite differences, is adopted for the space domain, with the discretisation parameter denoted by ∆x . The discrete value of u ( x , t ) and θ ( x , t ) becomes u j def := u ( x j , t ) and θ j def := θ ( x j , t ) , respectively, with j ∈ 1 , . . . , N . Note that other discretisation methds (finite volumes, finite elements, ...) could be used and will give the same MIM formulation.

The state space representation of the large original model

The straightforward semi-discretization of Eqs. (2) and (5) using central finite differences yields for j ∈ 2 , . . . , N -1 to:

∂u ∂t = Fo ∆x 2 • u j+1 -2 • u j + u j-1 , ( 6a
)
∂θ ∂t = Fo ∆x 2 • θ j+1 -2 • θ j + θ j-1 + 1 ∆x 2 • u j+1 -2 • u j + u j-1 . ( 6b
)
The problem is governed by Dirichlet boundary conditions for both fields u and θ :

u 1 = u ∞ , L ( t ) , u N = u ∞ , R ( t ) , θ 1 = 0 , θ N = 0 .
Thus, the semi-discrete equations ( 6) can be written in the matrix form:

U = Fo • A • U + Fo • B • Q , Θ = Fo • A • Θ + A • U + B • Q ,
where A ∈ M R N ×N and B ∈ M R N ×2 are matrices defined to ensure a second order accuracy in space:

A def := 1 ∆x 2 •            -2 1 0 . . . 0 1 -2 1 . . . 0 0 . . . . . . . . . 0 0 . . . 1 -2 1 0 . . . 0 1 -2            , B def :=            1 0 0 0 . . . . . . 0 0 0 1            , Q def := u ∞ , L u ∞ , T
The computation of the observable fields u s and θ s can also be formulated in a matrix form:

Y u = C • U , Y θ = C • Θ .
where

C ∈ M R N s×N s , so that Y u ∈ M R N s×1 and Y θ ∈ M R N s×1 are vectors.
In the end, the state space representation of the problem is formulated as:

                         U = Fo • A • U + Fo • B • Q , Θ = Fo • A • Θ + A • U + B • Q , Y u = C • U , Y θ = C • Θ . ( 7
)
The model is solved using a Runge-Kutta order 2 solver. It computes the field and its sensitivity at the points of interests x s , corresponding to the sensor locations. It is denoted as Large Original Model (LOM).

The MIM reduced order model

The MIM reduced order model enables to compute the field ũ s and θ s , s ∈ 1 , . . . , N s . They are gathered in the vectors Ỹ u and Ỹ θ . Both are computed from the unknown reduced vectors

X u ∈ M R N r ×1 and X θ ∈ M R N r ×1 .
The reduced model equations is assumed as:

                         Ẋ u = Fo • F • X u + Fo • G • Q , Ẋ θ = Fo • F • X θ + F • X u + G • Q , Ỹ u = H u • X u , Ỹ θ = H θ • X θ , ( 8
)
where the matrix F ∈ M R N r ×N r is assumed as diagonal. And the matrices

G ∈ M R N r ×2 , H u ∈ M R N s×N r and H θ ∈ M R N s×N r
are fulled. Indeed the general formulation of MIM (Eq. 8) is obtained by considering the transformation U = M X u applied to the LOM equations (Eq. 7). It can be noted that the columns of the matrix M ∈ M R N ×N r contains the eigenvectors of A and the diagonal matrix

F ∈ M R N r ×N r the N s eigenvalues of A.
The reduced model will thus give at time t the reduced vectors X u and X θ . To come back to the fields at the interest points, the operations Ỹ u = H u • X u and Ỹ θ = H θ • X θ are performed. These observable fields are obtained by the reduced model and thus denoted with the super script ∼ . It is interesting to remark that knowing matrices F , G, we have straightforwardly the reduced model for the sensitivity X θ . Here, the MIM model is solved using a Runge-Kutta order 2 method.

The construction of the MIM reduced model is performed by defining the differences between the reduced model and the complete model as a square residues functional J rom :

J rom N r , F , G , H u , H θ def := N s i=1 t f 0 ũ s ( t ) -u s ( t ) u s ( t ) 2 + θ s ( t ) -θ s ( t ) θ s ( t ) 2 dt .
The functional J rom is minimized to determined the order of the reduced model N r and matrices F , G , H u and H θ . The minimization procedure is carried out preliminary to the solution of the parameter estimation problem, using any optimizer such as quasi-Newton,PSO or genetic algorithm. This procedure is denoted as the learning step. Here the construction is performed for several values of the unknown diffusivity α . The set of thermal diffusivity used for the construction of the ROM is denoted Ω α .

Parameter estimation problem with reduced order model

The reduced order model is used in the framework of parameter estimation problem. For the sake of clarity, the methodology is explained for the unknown parameter denoted by:

p = α .
It is also distinguished the solution of the parameter estimation problem p • . It is different from the exact solution of the problem, noted p r , when known. The so-called a priori parameter, used as initial guess in the algorithm to solve the inverse problem, is denoted by p apr . The set of searched parameter is denoted Ω p . The experimental observations of the field u are written with the super-script m . Thus, the vector of observation is denoted

Y m ∈ M R N s×1 .

The optimal sequence of observations

To use the MIM ROM in the framework of parameter estimation problem, it requires a learning step to build a posteriori the matrices. Since this step is time consuming, the optimal experiment design (OED) methodology [START_REF] Rasooli | In-situ rapid determination of walls' thermal conductivity, volumetric heat capacity, and thermal resistance, using response factors[END_REF]31] is used to select a reduced measurement sequence of 3 days, denoted:

Ω oed t def := t ini , t ini + 3 ,
where t ini d is the beginning of the reduced sequence. Note that the reduced sequence respects the condition Ω oed t ⊂ Ω t . The measurement plan Π is introduced:

Π = t ini ,
The search of the optimal measurement plan Π • is carried out through the maximization of the D-criteria [16, 17]

Π • = arg max Π Ψ ,
where Ψ is the determinant of the Fisher matrix F:

Ψ def := det F( Π ) ,
and

F def := 1 σ 2 N s i = 1 Ω oed t θ s t 2 dt ,
with σ is the measurement uncertainty. This methodology is performed using the LOM and the a priori value of the unknown parameter. It enables to determine the optimal reduced sequence of measurement that maximize the parameter estimation accuracy. With this approach, the learning step of the MIM and the solution of the inverse problem are performed over a reduced sequence Ω oed t ⊂ Ω t .

Estimating the unknown parameter

The parameter estimation problem is solved over the reduced sequence Ω oed t to cut the computational efforts. It aims at determining the unknown parameter p verifying:

p = min Ω p J ( p ) ,
where J is the cost function defined using the least square estimator and the reduced model:

J : p -→ N s i=1 t f 0 ũ s ( p , t ) -u m s ( t ) 2 dt . ( 9
)
The minimisation of the cost function is realized using the Gauß algorithm [17,[START_REF] Ozisik | Inverse Heat Transfer: Fundamentals and Applications[END_REF]. It assumes the following condition on J :

∇ p J = 0 . ( 10
)
Thus, Eq. ( 10) is equivalent to:

N s i=1 t f 0 θ s ( p , t ) • ũ s ( p , t ) -u m s ( t ) = 0 . ( 11
)
Using the Taylor expansion of ũ around a parameter p k , equation (11) gives a system of N s equations:

ũ s ( p k , t ) + θ s ( p k , t ) • p -p k -u m s = 0 .
It enables to compute the next candidate p = p k+1 better than p k in an iterative procedure, written with the matrix formulation:

p k+1 = p k + Y θ ( p k ) • Y θ ( p k ) -1 • Y θ ( p k ) • Y m -Ỹ u ( p k ) . ( 12
)
In terms of computational time, the cost to compute the new candidate p k+1 is strongly reduced due to the use of MIM reduced order model. Indeed, it only manipulates the reduced vectors X u and X θ of size N r N s . Moreover, the MIM model provides straightforwardly the sensitivity functions of the observable fields.

The Algorithm 1 describes the main steps to solve the inverse problem. At the iteration k , the field and its sensitivity are computed using the MIM ROM and a candidate p k for the unknown parameter. With these results, the cost function between measurement and model predictions is evaluated. Moreover, using the sensitivity Y θ , the candidate p k+1 for the next iteration is computed. The algorithm runs until a maximum number of iterations N k is reached or one of the two following criteria is satisfied:

γ 1 η 1 & γ 2 η 2 ,
where η 1 and η 2 are user defined values. The criteria γ 1 and γ 2 evaluates the relative magnitude changes of the cost function and candidate parameter, respectively:

γ 1 def := J( p k+1 ) -J( p k ) 2 J( p k ) 2 , γ 2 def := p k+1 -p k 2 p k 2 .
Algorithm 1 Gauß algorithm using the MIM ROM.

N k & γ 1 η 1 & γ 2 η 2 do 5:
Compute the direct problem Y u and Y θ using reduced order model with Eq. ( 8) and a candidate p k 6:

Evaluate the cost function J using Eq. (9)

7:

Compute the next candidate for the unknown parameter p k+1 using Eq. [START_REF] Cornaro | Dynamic simulation and on-site measurements for energy retrofit of complex historic buildings: Villa mondragone case study[END_REF] 8:

Increment: k = k + 1 , p k ← p k+1 9:
Compute stopping criteria γ 1 and γ 2 with Eqs.

10: end

Synthesis of the methodology and metrics of its efficiency

The global methodology to solve the parameter estimation problem of heat transfer using MIM ROM is synthesized in Fig. 1. First, the OED methodology is used to define an optimal reduced sequence Ω oed t of the observations of the field. This results have two advantages noticeable at steps two and three. At step two, the learning step of the MIM is carried only for 3 days and not for the whole signal. It saves important computational efforts to determine the matrices of the ROM MIM. Then, at step three, the inverse problem is solved for a reduced sequence. Moreover, the sequence is optimal since it ensures the maximum accuracy of the retrieved parameters. Last, when the unknown parameter is estimated, the MIM ROM is computed considering the whole signal. The error between the numerical predictions and the whole sequence of experimental observations is computed to evaluate the reliability of the model. Several metrics are defined to evaluate the efficiency of the proposed method. The L 2 and L ∞ errors are denoted by ε 2 and ε ∞ , respectively. They are computed to evaluate the accuracy of the predictions of the numerical models:

ε 2 def := 1 N s N t N t n = 1 N s s = 1 u ( x s , t n ) -u ref ( x s , t n 2 , ε ∞ def := max x s , t n u ( x s , t n ) -u ref ( x s , t n ,
where u ref is a so-called reference solution, obtained generally by the LOM model, and N t the number of time step. In addition, the relative error ε r is calculated to evaluate the accuracy of the retrieved parameter:

ε r def := p • -p r p r .
Another important criteria is the computational efforts evaluated in the Matlab TM environment with a computer equipped with Intel i7 CPU and 32 GB of RAM. The measured computational time is denoted by t cpu s . A dimensionless version is also computed:

t cpu def := t cpu t cpu , ref ,
where t cpu , ref corresponds to the computational time required using the LOM.

Validation of the MIM model for parameter estimation problem

The construction of the MIM ROM and its use to solve parameter estimation problem is first explained on a simple case with simulated experimental data. Since the computational costs are reduced, step one of the methodology illustrated in Figure 1 is not carried out. Thus, next section presents the learning step to build the MIM ROM for the field and its sensitivity to the unknown parameter. Then, the model is used for parameter estimation problem.

Construction of the MIM ROM

The MIM ROM is built during the so-called learning step. For this, Eq. ( 1) is solved for a wall of length L = 0.5 m , an horizon of simulation t f = 24 h and the following initial and right boundary condition:

T ini = T ∞ , R = 20 • C .
The solution is generated for four cases at five points of interests defined as:

x s = 0.1 , 0.15 , 0.25 , 0.3 , 0.45 . ( 13
)
Three cases corresponds to the following diffusivity Ω α = 2.03 , 5.21 , 5.81 • 10 -7 m 2 . s -1 , corresponding to the ones of concrete, gypsum and mineral wool, combined with the left boundary condition:

T 1 ∞ , L = T ini + 5 sin 2 π 24 • 3600 t . ( 14
)
The last case is carried out for the diffusivity of concrete α = 5.21•10 -7 m 2 . s -1 and a slightly modified left boundary condition:

T 2 ∞ , L = T ini + 20 sin 6 π 24 • 3600 t . ( 15
)
To obtain each order of the MIM, PSO is used to perform the learning step with 100 iterations and a size of 50 populations. It can be noted that taking 100 iterations is more than sufficient to ensure the convergence.

Figures 2(a) and 2(b) show the error of the reduced order model with the order N . After the order 10 , the accuracy of the model becomes stable. Thus, a MIM ROM with an order N scaling with O( 5 ) enables to cut the computational complexity as well as ensure a sufficient accuracy of the predictions.

Accuracy for direct simulations

To evaluate the accuracy of the built model, the solutions y u and y θ are computed with the MIM ROM and the LOM for 50 values of α in the interval 2.03 , 6.37 • 10 -7 m 2 . s -1 . The following initial and boundary conditions are defined:

T ini = T ∞ , R = 20 • C , T ∞ , L = T ini + 5 sin 2 π 24 • 3600 t .
The final time is t f = 24 h and the wall of length L = 0.5 m . The points of interest are the same as defined in Eq. (13). The discretisation parameters is ∆t = 30 s for both models. For the LOM, 75 points of space discretisation are used.

The time evolution of the fields and the sensitivity at the five points of interests is shown in Figures 3(a) and 3(b) for two values of diffusivity, obtained with the MIM order 5 . It highlights that the model enables to simulate the whole dynamic of the observable outputs for thermal diffusivity different from the one used in the learning step.

Figure 4 shows the variation of the error according to the values of the diffusivity and for two different orders. The ROM of order 5 has a satisfactory error, lower than 10 -1 for both the field and its sensitivity. Moreover, the error is stable for the whole range of thermal diffusivity remaining inside the domain used during the learning step. As noted for α > 6 • 10 -7 m 2 . s -1 , the accuracy of the model is not satisfying anymore. Indeed, this diffusivity is out of the range used during the learning step.

Table 1 provides a synthesis of the effiency of the MIM ROM for a direct simulation. The accuracy of the model is almost similar between order 5 and 17 . In terms of computational time, the MIM ROM of order 5 enables to cut the computational cost by 87% compared to the LOM. Since the model reduction method is a posteriori, there is an inherent computational cost to build the MIM model. As a synthesis, a learning step composed of 4 signals is sufficient to build a model of low order and reliable for a wide range of parameter α.

Parameter estimation problem

The reliability of the MIM ROM is now evaluated for solving parameter estimation problem. For this, experimental observations are numerically generated using the LOM for a diffusivity

T ini = T ∞ , R = 20 • C , T ∞ , L = T ini + 3 • 1 -cos 2 π 12 • 3600 t .
It can be remarked that the boundary condition at T ∞ , L is very different from the one used during the learning step of the MIM ROM ( 14) and (15). The LOM is computed using ∆t = 30 s and ∆x = 6.8 mm . Then, the observations are projected at the points of interest considering a time step of 1 hour and adding a random normally distributed error. The standard deviation of the measurement is reported in Table 2. Both the LOM and the MIM ROM of order 5 are used to solve the parameter estimation problem. The initialguess is α apr = 2.1 • 10 -7 m 2 . s -1 , corresponding to a different value from those used during the learning step. The stopping criteria are set as η 1 = η 2 = 10 -14 and N k = 100 . The convergence of the algorithm to solve the parameter estimation problem is illustrated in Figures 5(a) and 5(b). The algorithm searches the unknown parameter until the criteria γ 2 on the cost function is lower than 10 -14 . It requires 8 and 10 iterations when using the LOM and the ROM, respectively. A comparison of the model prediction with the estimated parameters against the measurement is shown in Figures 6(a) and 6(b). A satisfying accuracy is observed for both models. The residuals have no specific signature as illustrated in Figures 6(c) and 6(d). Table 2 provides a synthesis of the results. The MIM ROM enables to cut the computational time to solve the inverse problem with a satisfying relative error of 4% on the estimated parameter.

Further numerical tests are carried for 50 values of thermal diffusivity. Figure 7(a) shows the evolution of the estimated parameter using the MIM ROM of order 5 . The algorithm requires around 5 iterations to estimate the unknown parameters with an relative accuracy of 4% . Figure 7(b) gives the variation of the computational time ratio between the algorithm using the LOM and the MIM ROM. The use of the reduced order model enables to reduce the computational effort by 60% .

Solving the parameter estimation problem for a real case study

Now the methodology illustrated in Figure 1 is applied to a realistic case study presented in the next section. Then, the reduced sequence of three days is determined using the OED methodology.

Description of the case study

The studied building is an old house, built at the end of the XIX th century and located in Saint Julien de Crempse, France. Further information on the architectural aspect of the building can be found in [START_REF] Cantin | Field assessment of thermal behaviour of historical dwellings in france[END_REF]. The investigated wall is oriented East and has a length of L = 43 cm .The wall is composed of a mixture of lime stones and clay mortar as shown in the picture Figure 8(c). Due to the vernacular traditional architecture, the stones arises from local around the house.

There is no information on the exact position of the stones.The issue is to estimate the global thermal diffusivity α of the wall by an equivalent model [14]. Note that the structural identifiability of this parameter has demonstrated in [13]. For this aim, the wall are monitored using five intrusive calibrated sensors HOBO TMC-6-HA as shown in Figures 8(a) to 8(c). Three sensors are installed inside the wall by obliquely drilling a whole of 13 mm. In the hole, the contact between the material and the sensor is ensured by a thermal conductive paste. The hole is then filled using an insulation material. The sensor located at the inside and outside surfaces are protected by insulation materials from incident radiation. The location of the sensor is set to assume that the heat transfer occurs mainly from the outside to the inside. An illustration of the experimental design is given in Figure 8(d).

The monitoring is carried out from January until April 2009. Experimental measures of temperature are available at x = 5 , 23 , 40 cm , as illustrated in Figure 9(b). Thus, N s = 3 are hold. It is not known if the sensors are located in the stone or in the mortar. The temperature T ∞ , L and T ∞ , R are also measured and given in Figure 9(a). Measures are available for Ω t = 0 , 2670 h , corresponding to four months. The measurement time step is 1 h . The experimental uncertainties σis computed using [START_REF] Taylor | The stability of the Dufort-Frankel method for the diffusion equation with boundary conditions involving space derivatives[END_REF]: where σ m = 0.1 • C is the sensor measurement uncertainty and σ x is the sensor position uncertainty. Due to experimental design based on drilling a large wall, the latter can be important. It is computed according to:

σ = 2 m + σ 2 x ,
σ x = ∂T ∂x x = x s δ x ,
where δ x = 2 cm is the uncertainty on the sensor position. The partial derivative ∂T ∂x is computed using the numerical model. Figure 10(a) shows the probability of the measurement uncertainty at the three points of observation. The uncertainty increases for the sensor located near the outside air conditions. For the initial conditions, a first order polynomial is fitted using the experiment:

T ini : x -→ T 0 • 1 + x 1 ,
where T 0 = 11.75 • C and 1 = -0.406 m . As shown in Figure 10(b), the error fitting is satisfactory with ε 2 = 0.09 • C 2 . The a priori knowledge from standards gives α apr = 2.8 • 10 -7 m 2 . s -1 .In addition, Table 3 gives the diffusivity of the different elements composing the wall.

Determining the reduced sequence of observations using the OED methodology

The OED methodology is now carried out to determine the optimal reduced sequence of three days Ω oed t of observations. The sensitivity function at the point of observations x s , s ∈ 1 , 2 , 3 are computed and shown in Figure 11(a). Then, the criteria Ψ is for each sequence of three days in Ω t . Figure 11(b) shows the variation of the dimensionless criteria according to the choice of t ini . By comparing Figure 11(a) and 11(b), it can be remarked that the optimal sequences occurs when the sensitivity of the parameter have high magnitudes of variations. The criteria verifies Ψ 0.98 for three initial days. According to these results, the following reduced sequence is adopted for the construction of the MIM ROM and then the estimation of the unknown parameters:

Ω oed t = 72 , 75 d .
The chosen monitoring time of 3 days arises from a compromise among several opposite criteria. Indeed, it could be reduced compared to traditional methods. From an pure inverse problem point of view, the longest is the measurement period, the highest is the accuracy of the parameter estimation. This can be proven by looking at the variation of the OED criteria Ψ according to the length of the sequence in Figure 12. The criteria is monotonously increasing according to Ω oed t . However, one has to deal with experimental constraints such as the acceptation of the building occupants. Furthermore, the longer is the measurement period, the higher is the computational time of the direct model and as a consequence the inverse problem algorithm. Figure 12 shows the variation the criteria Ψ added to the inverse of the computational cost of the algorithm. For this case, a good compromise is in three days of monitoring. Note that the influence of the length of the measurement period is investigated empirically in [START_REF] Cornaro | Dynamic simulation and on-site measurements for energy retrofit of complex historic buildings: Villa mondragone case study[END_REF]. For their case study, the authors suggest to carry at least three days of monitoring to have an accurate estimation.

Building the MIM ROM

For the reduced sequence Ω oed t , the boundary conditions are given in Figure 13(a). The initial condition is the following:

T ini : x -→ T 0 • 1 + x 1 - x 2 2
, where T 0 = 16.87 • C , 1 = 1.764 m and 2 = 0.5 m . The error fitting of the initial condition is satisfactory (ε 2 = 0.24 • C 2 ) as shown in Figure 13(b). The learning step is carrying considering four values of diffusivity Ω α = 1 , 5 , 10 , 70 • 10 -7 m 2 . s -1 . The training is carried until the order N = 20 taking around t cpu = 3 d . Figures 14(a

Parameter estimation

After the learning step the MIM ROM matrices are known. Thus, the model is now used to estimate the thermal diffusivity of the material. Figure 13(d) shows the variation of the sensitivity functions using the a priori value of the diffusivity. The functions have high magnitude of algorithm. As presented in Table 4, the value of the estimated parameter α • remains stable from N = 7 . The computational time to solve the inverse problem is cut by 5 compared to the LOM. Slightly differences in the computational time are remarked for the MIM ROM due to the number of iterations to converge. The model of order 20 is more accurate and requires less iterations to retrieve the unknown parameter.It is important to note that the experimental design relies on embedded sensors inside the wall to obtain the experimental observations. This may be a restriction of the methodology for other case study. Other experimental devices could be used such as heat flux sensor or infrared camera. The methodology to search the optimal experimental design could be extended for such approaches. Furthermore, a MIM ROM needs to be constructed according to the new type of observable field.

Reliability of the model

To evaluate the reliability of the model with the estimated thermal diffusivity α • = 1.06 • 10 -6 m 2 . s -1 , the numerical predictions computed with the MIM ROM N = 10 are com- give the probability of the error with observation over the whole sequence of measurement. The error is compared with the probability of the measurement. The prediction of the model are more precised using the estimated parameter. Indeed the probability has a smaller standard deviation. However, at the point x 3 , errors of larger magnitude occur for probabilities of both models. It reveals that some phenomena are not considered in the model and the mathematical formulation should be improved by considering for instance heterogeneous diffusivity in the wall.

To verify the reliability of the proposed model, the inside heat flux j and the thermal loads E are computed according to:

j def := -k ∂T ∂x x = L , E def := Ω t j q dτ .

Conclusion

The uncertain thermal properties of existing building walls can be determined by solving parameter estimation problems. However, such problems have important computational cost for at least two reasons. The inverse problem algorithm requires to solve the direct problem several times. In addition, the experimental measurement are carried on-site implying simulations over several month under climatic variations. These complexities increase the computational requirements.

To answer this issue, innovative methodology is proposed based on two concepts. First, a reduced order model based on MIM method is used to reduce significantly the computational time of the direct problem without loosing accuracy. The methodology is presented in Section 3. The reduced model is based on a state space representation where the matrices are built during a learning step. The latter is based on a minimization procedure between the predictions of the reduced and complete models. An interesting point is that the model is built to compute the field of interest and its sensitivity to the unknown parameter. The sensitivity is known straightforwardly using the matrices of the reduced model.

The second concept is the optimal experiment design methodology, described in Section 4. It is employed to determine a reduced sequence of observations. Three advantages are enhanced with this approach. First, it reduces the inherent computational cost of a posteriori model reduction methods since the learning step is carried for a reduced sequence. Then, the chosen sequence is optimal to estimate the parameter with accuracy. In addition, the inverse problem is solved only for a few days of observations. It reduced again the computational effort to retrieve the parameter. The solution of the inverse problem is obtained with the Gauß algorithm since the reduced order model computes the field and its sensitivity.

A first case study is proposed in Section 6 to validate the MIM model reduction method. This step is important to evaluate the efficiency in terms of accuracy and computational cuts. The Regardless the inherent cost of the learning step, the model cut by 5 the computational cost of the direct problem. The model is also evaluated in the framework of inverse problem with simulated experimental observations. The unknown parameter is estimated with an error lower than the measurement uncertainty and a reduced computational time.

After this simple validation case, the whole methodology is applied to a more realistic one. The issue is to estimate the thermal diffusivity of an old building wall. The latter is monitored during four months with three sensors drilled inside. With the OED methodology, a sequence of three days is identified as optimal. Then, the MIM ROM is built for the reduced sequence and a signal composed of 4 values of thermal diffusivity. Then, a model of order 10 is chosen to solve the inverse problem. The estimated diffusivity is three times higher than the one provided by standards. In terms of computational efficiency, the parameter is estimated with an algorithm 5 times faster. Eventually, the reliability of the model is evaluated by comparing the predictions of the MIM model with the experimental observations of the complete period. It highlights that the global methodology is efficient to calibrate the model with a reduced computational effort.

Future works should focus on improving the definition of the mathematical model by dealing with anisotropic thermal diffusivity in the building wall. In addition, to evaluate the real efficiency of the methodology, the Gauß algorithm combined with the MIM ROM should be installed in embedded system for fast estimation of building material properties.

Nomenclature and symbols

Figure 1 .

 1 Figure 1. illustration of the methodology.

Figure 2 .

 2 Figure 2. Evolution of the error of the MIM ROM model according to the order N .

1 Figure 3 .

 13 Figure 3. Evolution of the field and the sensitivity at the points of interest for two values of diffusivity.

Figure 4 .

 4 Figure 4. Variation of the error of the MIM ROM, Orders 5 (a,b) and 17 (c,d), for the whole interval of thermal diffusivity, for the field (a,c) and its sensitivity (b,d).

Figure 5 .

 5 Figure 5. Evolution of the estimated parameter (a) and of the convergence criteria γ 2 according to the iteration.

Figure 6 .

 6 Figure 6. Comparison of the (simulated) measurement and the model prediction using the estimated parameter for the LOM (a,c) and the MIM ROM Order 5 (b,d).

Figure 7 .

 7 Figure 7. Evolution of the estimated parameter according to the iteration (a) and of the computational time ratio of the algorithm using the MIM ROM Order 5 (b) for 50 values of initial diffusivity.

Figure 8 .

 8 Figure 8. Illustration of the experimental design with picture of the inside sensors (a,b), the outside sensors (c) and the experimental design (d).

  ) and14(b) shows the evolution of the error according to the order of the model. It can be remarked that a model of order N = 10 provides a satisfactory accuracy.

Figure 9 .Figure 10 .

 910 Figure 9. Time variations of the boundary conditions (a) and the measured observation data (b).For the sake of clarity, the measurement uncertainty is not presented.

Figure 11 .

 11 Figure 11. Time variations of the sensitivity functions of the problem (a) and variation of the OED criteria according to the choice of the initial time of the reduced sequence Ω oed t (b).

Figure 12 .Figure 13 .

 1213 Figure 12. Variation of the OED criteria according to the length of the reduced sequence Ω oed t

Figure 14 .Figure 15 .

 1415 Figure 14. Evolution of the error of the MIM ROM model according to the order N trained over the reduced sequence Ω oed t .

Figure 17 (

 17 Figure17(a) shows the variation of the thermal flux. It has important magnitude inducing high thermal loads as noted in Figure17(b). The thermal loads are computed for perturbations of 20% of the estimated diffusivity. It can be noted that the relative error on the loads remains stable around 5% . It confirms that the approximations introduced from the whole methodology (uncertainty measurement, reduced order model or estimation on a reduced sequence) are minor compared to the ones in the mathematical model itself.

Figure 16 .

 16 Figure 16. Comparison of the numerical predictions for the last week (a,c,e) and probability of the error between the predictions and the experimental observations for the whole period Ω t (b,d,f).

Figure 17 .

 17 Figure 17. Time evolution of the inside heat flux (a) and variation of the thermal loads according to small perturbation of α • (b).

Table 1 .

 1 Efficiency of the MIM ROM.

	Model Order	ε ∞		t cpu s		t cpu -
			u s	θ s	Learning step Direct simulation Learning step Direct simulation
	LOM	75	-	-	-	0.60	-	1
	ROM	17	0.08 0.06	2.4	0.10	4	0.16
	ROM	5	0.07 0.06	2.4	0.08	4	0.13
	α

r = 5.21 • 10 -7 m 2 . s -1 . The initial and boundary conditions are the following:

Table 2 .

 2 Results of the parameter estimation problem.

		Point of observation x s m				Error ε r -on
		0.1	0.15	0.25	0.3	0.45		Computational time	parameter estimation
	Model			Error ε 2			Order t cpu s	t cpu -	α apr	α •
	LOM	0.059 0.092 0.073 0.10 0.073	75	4	1	-0.59	-0.01
	ROM 0.068	0.11	0.083 0.10 0.074	5	0.9	0.23	-0.59	-0.04
		Measurement uncertainty • C
		0.94	0.70	0.46	0.30	0.13

Table 3 .

 3 Wall stratigraphy.

	Layer	Thermal diffusivity , m -2 . s -1	Thickness cm
	Lime stones	2.7 • 10 -7	O( 20 )
	Mortar	5.5 • 10 -7	O( 4 )
	Equivalent homogeneous wall	2.8 • 10 -7	43

Table 4 .

 4 Results of the parameter estimation problem, noting that α apr = 2.8 • 10 -6 m 2 . s -1 .

			Estimated parameter	Computational time	Iterations to
	Model Order	α • • 10 7 m 2 . s -1	t cpu s	t cpu -	to converged
	LOM	100	1.064	572	1	19
	ROM	5	1.097	114	0.2	24
	ROM	7	1.064	112	0.19	23
	ROM	10	1.062	101	0.17	21
	ROM	15	1.065	102	0.18	21
	ROM	20	1.067	94	0.16	19

Acknowledgments

BK thanks Johnathan Gerardin for the fruitful discussions on MIM.


Mathias Cuny

Arnaud Lapertot

Jian Lin

Benjamin Kadoch
email: benjamin.kadoch@univ-amu.fr

Olivier Le Metayer

Multi-criteria optimization of an earth-air heat exchanger for dierent French climates

Keywords: earth-air heat exchanger, multiple-criteria decision-making, optimization multi-criteria, renewable energy, sensitivity analysis

The Earth-Air Heat Exchanger (EAHE) is a system for cooling or preheating the blown air into a building. A modeling of this system is proposed from a model for the ground and a model for the exchanger. The results are in agreement between the outside air temperature of the numerical model and the experimental temperature measurements. A sensitivity analysis with factorial plans is carried out to determine the most impactful parameters, taking into account two energy criteria and one economic criterion. The results show that tube radius, tube length, air velocity, burial depth and soil nature are the most impactful parameters. A multi-criteria optimization study with genetic algorithms is then performed to determine the Pareto front. The criteria do not evolve in the same trend because when the cost of the energy recovered decreases, the coecient of performance and the EAHE eciency will deteriorate, and reciprocally. Finally, a multiple-criteria decision-making is carried out using the TOPSIS method to provide the optimal pipe conguration.

The optimal solution provides a large tube length, an intermediate burial depth, and a small air velocity and tube radius. The EAHE can achieve strong energy performance in any French climate.

Introduction

The residential sector accounts for 40 % of total energy consumption in developed countries [START_REF] Pérez-Lombard | A review on buildings energy consumption information[END_REF].

To ght global warming, one of the European Union's strategies is to reduce energy consumption in buildings [START_REF]The European Parliament and the Council of the European Union[END_REF]. New regulations impose thermally insulated buildings using renewable energy systems and high-eciency energy systems. In this context, ventilation heating systems are favored.

For example, in Europe, heating, ventilation and air conditioning (HVAC) accounts for 68 % of the energy consumed in the residential sector [START_REF] Pérez-Lombard | A review on buildings energy consumption information[END_REF].

The earth-air heat exchangers (EAHE) are renewable energy systems using ground temperature to preheat or cool the building. This system also allows the ventilation of the exhaust air of a building and the injection of renewed air from outside. The EAHE can improve the performance of a controlled mechanical ventilation for the heating system. The global energy performance of the system is thus improved by 29 % in winter and between 36 % and 46 % in summer [START_REF] Ascione | Earth-to-air heat exchanger for nzeb in mediterranean climate[END_REF]. In addition, the EAHE are very well suited to provide air ventilated at comfort temperature during the summer [START_REF] Yusof | Thermal analysis of earth-to-air heat exchanger using laboratory simulator[END_REF]. Such energy system is particularly suitable for contrasting climates, either for preheating [START_REF] Li | Feasibility study on applications of an earth-air heat exchanger (eahe) for preheating fresh air in severe cold regions[END_REF], or for cooling [START_REF] Fazlikhani | Numerical analysis of the eciency of earth to air heat exchange systems in cold and hot-arid climates[END_REF]. This energy system is also suitable for temperate climates such as [START_REF] Cuny | Inuence of coating soil types on the energy of earth-air heat exchanger[END_REF], [START_REF] Lin | Impact of soil moisture on the longterm energy performance of an earth-air heat exchanger system[END_REF] and [START_REF] Thiers | Thermal and environmental assessment of a passive building equipped with an earth-to-air heat exchanger in France[END_REF]. The use of an EAHE and its energy performance are highly dependent on climatic conditions [START_REF] Li | Experimental study of a ground sink direct cooling system in cold areas[END_REF]. The rst objective is to investigate if this system can be applied to temperate climates such as France. The soil type also allows improving the performance of the system. Indeed, the thermal performance can reach 15.9 % with a soil composed of sand [START_REF] Cuny | Inuence of coating soil types on the energy of earth-air heat exchanger[END_REF]. Combined with a dual-ow controlled mechanical ventilation, the system's coecient of performance (COP ) is dened by the ratio between the heating capacity and the system consumption power and the COP can reach 16.3 [START_REF] Li | Feasibility study on applications of an earth-air heat exchanger (eahe) for preheating fresh air in severe cold regions[END_REF].

The second objective is to apply a methodology of sensitivity analysis, multi-criteria optimization and multi-criteria decision-making. The aim is to determine the optimal design parameters.

Sensitivity analysis quanties the impact of parameter variation on dened criteria and makes it possible to classify the inuence of the parameters. This method can be used to selected main pa-rameters. This EAHE energy system can involve a large number of parameters and the sensitivity analysis could be used to reduce the number of parameters retained for the optimisation process [START_REF] Saltelli | Global Sensitivity Analysis[END_REF]. Moreover, a multi-criteria optimization study is carried out with the main parameters. The aim is to determine the optimal design parameters that gives maximum performance and minimum cost. Indeed, the energy system must be optimized to be ecient and cheap. The main design 40 parameters of an earth-air heat exchanger are: the material, diameter and length of the tube, the air ow velocity and the burial depth of the tube [START_REF] Bordoloi | An intense review on the latest advancements of earth air heat exchangers[END_REF]. Concerning the COP of an EAHE, it can be optimal for a given length. Indeed, the dierence between the air temperature owing through the EAHE and the ground temperature decreases but pressure losses increase [START_REF] Benhammou | Parametric study on thermal performance of earth-to-air heat exchanger used for cooling of buildings[END_REF]. It is therefore not necessary to use exchangers of large length. 45 The third objective is to apply this energy system and this methodology for dierent French climates. Five French climatic zones are studied (H1a, H1b, H2b, H2c, H3), as illustrated in gure 1. Each climate is composed of its own specicity. Indeed, France is surrounded by a continental Figure 1: French climatic zones [START_REF]Les diérentes zones climatiques françaises[END_REF]. climate in the east, an oceanic climate in the west and north, and a Mediterranean climate in the south. Some French regions is constituted by a mountain climate due to the Pyrenees, the Alps, the Massif Central and the Vosges. Indeed, the continental climate is the most contrasting climate because it is characterized by a cold winter and a hot summer. Thus, the oceanic climate is inuenced by the marine currents of the Atlantic Ocean, which provides a mild winter and a cool summer. Afterwards, the Mediterranean climate is composed of a hot and a dry summer and a mild winter. This climate has a lot of sunshine and wind. Then, the mountain climate is constituted by a mild and short summer, and a long and very cold winter. The dierence between the climatic zones is based on meteorological data. Indeed, the outside air temperature and coecients used to model the ground temperature are specic to each zone. Figure 2 shows the evolution of the outdoor air temperature for the dierent climatic zones in France [START_REF]Les données météorologiques RT[END_REF]. Zone H3 has a stronger average temperature and a smaller temperature amplitude than the other zones.

In addition, the H1b and H2b zones can have a very low air temperature in winter. Finally, zone H2c gives the most contrasted climate since the lowest and highest temperature are obtained for winter and summer, respectively. Indeed, the H1a and H2b zones have an oceanic and continental climate. The H1b and H2c zones not only have a continental climate but also a mountain climate.

Zone H3 has a predominantly Mediterranean climate and a bit of mountain climate.

In this article, a multi-criteria optimization methodology is performed for dierent French climates.

In section 2, the energy system and the dierent methods used will be described. Indeed, the EAHE is implemented from two distinct models. The rst one simulates the thermal behaviour of the ground by computing the ground temperature, which is a function of depth and time. The second one models the air temperature between the inlet and the outlet of the EAHE. This annual dynamic modeling includes also the calculation of the power required by the fan to impose an air velocity.

The experimental measurements and the numerical estimates of the outlet air temperature are compared in order to validate the model. In this study, several geometric parameters and the criteria are taken into account: two energy criteria (the EAHE eciency and the COP ) and an

Time [months]

H 1a H 1b H 2b H 2c H 3
Figure 2: Temporal evolution of the outside air temperature for zone H1a (blue curve), zone H1b (red curve), zone H2b (orange curve), zone H2c (purple curve) and zone H3 (green curve) [START_REF]Les données météorologiques RT[END_REF]. economic criterion (the cost of the energy recovered). Section 3 describes the results obtained by the sensitivity analysis, multi-criteria optimization and multi-criteria decision-making methods. Indeed, a sensitivity analysis is carried out to nd the most impactful variables. Then, the design parameters are optimized according to the criteria to nd the set of solutions that minimize the cost of the energy recovered, and maximize the EAHE eciency and the COP . Finally, the optimal system pipe conguration is provided for ve climatic zones studied (H1a, H1b, H2b, H2c, H3).

Modeling description and Methods

Modeling description

The models used for the EAHE system is based on studies [START_REF] Thiers | Thermal and environmental assessment of a passive building equipped with an earth-to-air heat exchanger in France[END_REF] and [START_REF] Benhammou | Parametric study on thermal performance of earth-to-air heat exchanger used for cooling of buildings[END_REF]. The energy system is described in this subsection.

Description of the global system

The EAHE is composed of a pipe buried in the ground at a depth z. The tube is characterized by: a length L, an inner radius r and a thickness e. The air velocity v air owing through the EAHE is imposed by a fan located at its outlet (gure 3). The EAHE is the area of thermal contact between two dierent sources: the air owing through the tube, the soil surrounding the tube. In winter, the ground temperature is warmer than the outside air temperature. The air owing through the EAHE heats up thanks to the warmth of the soil. The outlet air is warmer than the inlet air of the EAHE. In summer, the opposite phenomenon occurs and the outlet air is colder than the inlet air. The EAHE is illustrated by a tube and the ground by a volume. The inlet air temperature of the tube T air, in corresponds to the outside air temperature T air, ext . This temperature is easily measurable from hourly meteorological databases. However, the ground temperature at the soiltube interface is dicult to obtain because the installation of a measuring device at the burial depth of the tube is problematic. In addition, these measurements are only valid for a single EAHE conguration. In this context, a system modeling is necessary. The hypothesis used for the EAHE model is taken from [START_REF] Thiers | Thermal and environmental assessment of a passive building equipped with an earth-to-air heat exchanger in France[END_REF] and [START_REF] Benhammou | Parametric study on thermal performance of earth-to-air heat exchanger used for cooling of buildings[END_REF]:

1. a one-dimensional problem is considered; 2. the soil is considered as a homogeneous and isotropic medium; 3. the air and ground thermophysical properties are constant and independent of temperature.

The system is divided into one model for the ground and one for the EAHE.

The soil thermal model

The soil thermal model provides the ground temperature for any depth and at any time without taking into account the EAHE operation. The eects of mass and heat transfer are considered negligible due to small dierences in ground temperatures [START_REF] Saix | Physical Model for the Study of Mass and Energy Transfers in the Non-Saturated Layer of Soil Located Above a Solar Energy Storage Zone[END_REF]. The soil is considered as a semi-innite homogeneous medium with constant thermophysical properties: thermal conductivity λ ground , density ρ ground and mass thermal capacity c p ground [START_REF] Thiers | Thermal and environmental assessment of a passive building equipped with an earth-to-air heat exchanger in France[END_REF]. In this study, the soil thermal model takes only into account thermal conduction.

Considering these hypotheses, the propagation of a sinusoidal temperature signal in a semi-innite medium has an analytical solution [START_REF] Sacadura | Fourth Printing. Techniques et Documentation[END_REF]. In the investigated case, the sinusoidal temperature is the outside air temperature, it is written as follows:

T air, ext = T mean + i=3 i=1 [T i sin (ω i t -φ i )] (1)
with T mean the average temperature over the tube, T i the amplitude temperature, ω i frequency oscillations, φ i the phase shift. The oscillations consider an annual frequency ω 1 , a monthly frequency ω 2 and a daily frequency ω 3 to take into account the annual, monthly and daily meteorological uctuations. The amplitudes and the frequency are determined by the least square's method applied to meteorological data. The values of the amplitude and the phase shift are given in table 1.

These values depend only on climatic zones. By applying this signal to the surface of a semi-innite medium, the ground temperature is therefore given by:

T ground (t, z) = T mean + i=3 i=1 T i exp -z ω i 2 α ground * sin ω i t -φ i -z ω i 2 α ground (2)
Zone T mean T 1 with α ground the thermal diusivity of the ground, which is dened by:

T 2 T 3 ω 1 ω 2 ω 3 φ 1 φ 2 φ 3 Unit • C • C • C • C s -1 s -1 s -1 - - - * 10 -7 *
α ground = λ ground ρ ground * c p ground (3)

The EAHE model

The EAHE modeling gives the outlet air temperature. The air owing through the EAHE is dened by a mass ow rate ṁair and a constant mass heat capacity c p air . Depending on the direction of air 125 ow, the tube is decomposed into several cylinders with the same length ∆x (gure 4). For each cylinder, the air and ground temperatures are considered uniform. The power balance is written by:

ṁair c p air * dT e dx = U A air-tube * (T ground -T e ) (4)
The mass of the tube is negligible compared to the mass of the ground. Consequently, the tube is characterized by its thermal conductance U A air-tube . This variable is dened by:

1 U A air-tube = 1 h c * (2πr ∆x) + 1 2πr ∆x λ tube * ln r + e r (5)
where λ tube is the thermal conductivity of the tube and h c is the convective exchange coecient between the air and the tube. h c is computed using formulations in forced convection [START_REF] Ozisik | Heat Transfer: A basic approach[END_REF] and it is written by: At each time step, the outlet air temperature T s is considered as the new inlet air temperature T e for the following cylinder. The electrical power consumed by the fan at the output of the EAHE is dened by:

h c = 0.036 λ air 2r * Re 0.8 P r 1 3 * 2r ∆x 0.055 (6)
P f an = 1 η f an ṁair ρ air ∆p (7)
with η f an the electromechanical eciency of the fan. Pressure losses of the system ∆p are calculated from [START_REF] Benhammou | Parametric study on thermal performance of earth-to-air heat exchanger used for cooling of buildings[END_REF] and are described by:

∆p = f L 2r + 2 * 1.3 ρ air v 2 air 2 (8)
where the regular pressure losses coecient f is dened by:

f = (1.82 log(Re) -1.6) -2 (9)

Validation

The results from the EAHE model proposed in this study are compared to the experimental data of an EAHE installed on the IUT Robert Schuman geothermal platform (University of Strasbourg), including the meteorological data [START_REF] Lin | Impact of soil moisture on the longterm energy performance of an earth-air heat exchanger system[END_REF]. The procedure of model identication and boundary conditions are explains and details in [START_REF] Cuny | Inuence of coating soil types on the energy of earth-air heat exchanger[END_REF].

The experimental EAHE is buried at 1.03 m in a multilayer soil: 10 cm of plant material, 60 cm of natural backll and 50 cm of sand surrounding the tube. In addition, the EAHE is composed of a tube with a length and diameter equal to 17.5 m and 20 cm, respectively. Concerning the operating mode of the experimental measurements, the fan extracts the outside air at a speed of 0.51 ± 0.02 m.s -1 at the tube outlet. Outside and outlet air temperatures are measured every 20 minutes with sensors PT100 with a accuracy of ± 0.1 • C. The geothermal platform is schematized by the gure 5.

Equations 2 and 4 are solved numerically at each time step during a full year using Runge-Kutta algorithm 4 th order with OpenModelica software [START_REF] Consortium | Openmodelica[END_REF]. In this study, the time step is 10 minutes and the hourly meteorological data are linearly interpolated. Figure 6 shows the evolution of the output air temperatures obtained from experimental data T exp air out and the numerical model T num air out .

Thermal behaviour of the EAHE model is similar to the measurements. The relative average error over the full year is 5.7 %. The modelled data are in good agreement with the experimental data.

Methods

The methods used for sensitivity analysis, multi-criteria optimization and multi-criteria decisionmaking is dened in this subsection.

Denition of variables and criteria

The To quantify the EAHE performance, three criteria are dened: two energy criteria such as the average eciency η EAHE and the average coecient of performance COP of the exchanger, and an economic criterion like the cost of the energy recovered (CER).

The average eciency of an EAHE η EAHE is the ratio between the energy recovered of the EAHE and the recoverable energy over a one-year period T [START_REF] Benhammou | Parametric study on thermal performance of earth-to-air heat exchanger used for cooling of buildings[END_REF]:

η EAHE = T 0 (T air, in -T air, out ) dt T 0 (T air, in -T ground ) dt (10)
The average COP is expressed by the ratio between the energy recovered and the energy consumed by the EAHE [START_REF] Benhammou | Parametric study on thermal performance of earth-to-air heat exchanger used for cooling of buildings[END_REF]:

COP = ṁair c p air * T 0 (T air, in -T air, out ) dt T 0 (P f an ) dt (11)
The cost of the energy recovered CER is dened by the ratio between the total cost of an EAHE (installation cost, maintenance cost and operating cost) over its lifetime N year = 20 years and the energy recovered by the EAHE [START_REF] Launay | Analysis strategy for multi-criteria optimization: Application to inter-seasonal solar heat storage for residential building needs[END_REF].

CER = cost installation + cost operation + cost maintenance N year ṁair c p air * T 0 (T air, in -T air, out ) dt (12)
The installation cost estimates the cost of workmanship such as earthworks and the cost of materials such as ventilators, air intake terminals and collectors. The installation cost is a function of tube length and burial depth. The operating cost is determined by the product of the electrical energy consumed by fans and the electricity prices. The maintenance cost quanties the change of the lters three times per year.

Sensitivity analysis

The model validity and parameter sensitivity are studied in this section. This study consists in identifying and ranking the most impactful inputs, determining the non-impactful inputs and mapping the behaviour of the outputs in relation to the inputs [START_REF] Saltelli | Global Sensitivity Analysis[END_REF]. To carry out this study, the factorial plan method is used as it can quantify the contribution of input parameters and their interactions. The two-level factorial plan requires 2 k calculations with k the number of inputs [START_REF] Montgomery | Design and analysis of experiments[END_REF].

Therefore, factorial plans require more calculation than methods that vary one parameter at a time, but they do provide information on the interactions between variables.

The factorial plan procedure is to approximate the objective function y by a linear system [START_REF] Launay | Analysis strategy for multi-criteria optimization: Application to inter-seasonal solar heat storage for residential building needs[END_REF].

y = γ 0 + k i 1 =1 γ i 1 X i 1 + k i 1 =1 k i 1 >i 2 γ i 1 i 2 X i 1 i 2 + ... + γ i 1 i 2 . . . i k X i 1 i 2 . . . i k ( 13
)
with γ 0 the mean of the objective function, γ i the vector of the main eects of the parameters, γ i 1 i 2 the vector of the interaction eects between parameters. X i is composed of the normalized parameters vector.

∀ i ∈ [1, k], X i = x i -x ref i (x max i -x min i )/2 (14)
where x i is characterized by the parameter value, x ref i the reference value, x max i the maximum value,

x min i the minimum value. For the two-level factorial plan method, x i is equal to its maximum or minimum value. Then the coecients γ are obtained by solving the linear system 13. To obtain the contribution of the dierent parameters, the sum of squares must be dened.

∀ i ∈ [1, 2 k -1], SS i = 2 k γ 2 i ( 15
)
Finally, the contribution percentage of each parameter and their interactions is thus determined by the following equation. P C which quantify the importance of parameters and couplings for each objective. The coupling corresponds to the interaction between parameters.

∀ i ∈ [1, 2 k -1], P C i = SS i k i=1 SS i (16)
In the investigated study, inputs correspond to variables or parameters (k = 12) and outputs referred to criteria or objectives. The most impactful parameters is determined using the factorial plan method [START_REF] Montgomery | Design and analysis of experiments[END_REF].

Multi-criteria optimization

A multi-criteria optimization problem seeks to optimize several objective functions simultaneously.

However, it does not provide an optimal solution because some objectives are not necessarily comparable. The works of Vilfredo Pareto [START_REF] Ehrgott | Multicriteria Optimization[END_REF] show that it is impossible to improve all criteria simultaneously, which means that the improvement of some objectives will necessarily degrade the other objectives. The result of the multi-objective optimization gives a set of solutions that have a good compromise for the dierent criteria, including the best solution of each criterion. These solutions come together on a curve called the Pareto front, which is composed of non-dominated solutions. A solution X is dominated by another solution Y if only, for any objective f i , f i (X) is less than or equal to f i (Y ) with at least one strict inequality. Many optimization methods exist as Nelder-Mead method, simplex method, particle swarm algorithms and genetic algorithms [START_REF] Arora | Optimization: Algorithms and applications[END_REF].

Genetic algorithms are suitable for constraints and multi-objectives problem. Genetic algorithms are probabilistic optimization methods based on the evolution of the species. They were developed originally by the works of John Holland [START_REF] Holland | Adaptation In Natural And Articial Systems[END_REF]. Thanks to selection, crossing and mutation processes, individuals change over generations and the best performing species tend to survive in their natural environment.

In this case, genetic algorithms NSGA II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] are performed with a percentage of mutation and crossing of 50 %, and with 100 individuals and 1000 generations, which corresponds to 100,000 evaluations per zone.

Multiple-criteria decision-making

Decision aid is a technique to choose on the best solution from a set of possible solutions. There are a number of multiple-criteria decision-making methods. The best known and used are the weighted sums and TOPSIS [START_REF] Tzeng | Multiple attribute decision making: Methods and applications[END_REF]. The Technique for Order Preference by Similarity to Ideal Solution TOPSIS is a decision aid technique developed originally in [START_REF] Tzeng | Multiple attribute decision making: Methods and applications[END_REF]. It consists in selecting the solution that has the shortest distance from the ideal solution.

In the investigated case, the best pipe conguration of the EAHE is selected using the TOPSIS method. The ideal point is the maximum of the η EAHE and the COP , and the minimum of the CER.

In this article, the optimization procedure is carried out in three steps. First, sensitivity analysis is used to determine the most impactful parameters in order to reduce the number of parameters.

Then, multi-criteria optimization gives the set of solution that has a good compromise for the dierent criteria. Finally, the multi-criteria decision-making method selects one solution in the Pareto front. This selected solution is called the optimal solution.

3. Results and discussions

Sensitivity analysis

This section shows the results of the sensitivity study in order to select the most impactful parameters for multi-criteria optimization. The minimum x min and maximum x max limits for each variable are dened according to their reference values x ref and a variable β varying from 5 % to 35 %: Furthermore, ground thermal conductivity λ ground has a percentage, which increases with β and reaches 3 %. Other parameters and couplings have no impact for the η EAHE . The second graph in gure 7 shows the results for the COP . Tube radius is predominant for the COP and decreases between 75 % and 60 % when β increases. Moreover, air velocity contributes 20 % of the general evolution of the COP . Similarly, parameter L and the coupling r -v air are small but increase with β, they exhibits 12 % and 2 %, respectively. Then, the burial depth has a small contribution, which is dependent on β. Finally, other variables have no impact for the COP . The graph in gure 7 exhibits the contribution percentages for the CER. Tube radius r is the dominant variable for the CER. This parameter decreases with β, because it goes from 70 % to 60 % when β increases.

     x min = x ref * (1 -β) x max = x ref * (1 + β)

Parameter sensitivity

In addition, the contribution of air velocity and tube length are important independently of β.

The inuences of v air and L are equal to 20 % and 6 %, respectively. Similarly, the parameter ρ ground * c p ground and the coupling between r and v air have a low contribution percentage but they increase with β. The contribution of the other parameters is negligible for the CER.

Physical interpretation

In summary, the most impactful parameters are tube radius, tube length, air velocity, burial depth and soil nature for three criteria considered. The impact of each parameter on the objectives are explained by carrying out a physical interpretation.

Tube radius has an inuence on the EAHE eciency, on the coecient of performance and on the cost of the energy recovered. Indeed, tube radius is correlated to the outlet air temperature of the EAHE (eq 10-12) and has an impact on the pressure drop and the fan power (eq 7-8). As a consequence, when r increases, ∆p and ∆T =|| T air out -T air in || decrease, so the heat exchanges are lower quality. ∆p and ∆T decrease independently. Therefore, when r increases, the CER increases, the COP and η EAHE decrease.

In addition, air velocity has a contribution to the three objectives because v air is directly related to the air mass ow rate. Indeed, the fan power (eq. 7) and the power balance (equation 4)

depend on the mass ow rate of air. Finally, the air velocity is correlated to the operating cost of the system which depends on the fan power. Therefore, when v air increases, the heat exchanges between soil and air, ∆p and cost operation increase. Thus when v air increases, the CER increases, the η EAHE and COP decrease.

Similarly, tube length has an inuence on the criteria. Indeed, pressure drops (eq. 8) and the outlet air temperature of the EAHE depend on tube length. It aects also the installation cost which is also related to the installation cost of the system. When L increases, the ground-air heat exchanges, ∆p and cost installation increase. Consequently, when L increases, the η EAHE and the COP increases, the CER decreases.

Furthermore, burial depth z has an inuence on the three objectives. According to the equation 2, the ground temperature depends on z and has a contribution on the outlet air temperature of the EAHE (eq. 4) and on the EAHE eciency (eq. 10). Burial depth is also correlated to the installation cost of the system. So, when z increases, cost installation increases, and T ground , T air out increase in winter and decrease in summer, so the heat exchanges are better. Therefore, when z increases, the COP increases, the CER and η EAHE decrease.

Moreover, the soil thermophysical properties λ ground and ρ * c p ground have an impact on the η EAHE , the COP and the CER. Indeed, the ground temperature is dened by a function of the ground thermal diusivity (eq. 2). As previously for burial depth, the ground temperature is related to the outlet air temperature of the EAHE (eq. 4) and thus to the three objectives. When α ground increases, the heat exchanges between the soil and the air are lower quality so the COP decreases and the CER and the η EAHE increase.

Choice of parameters for optimization

Based on the sensitivity study, the variables selected for the optimization of the EAHE are: tube radius r, tube length L, burial depth z, air velocity v air and soil nature characterized by its thermal diusivity α ground . These parameters have a signicant inuence on the performance of an EAHE. The results allowed obtaining the same design parameters given by [START_REF] Bordoloi | An intense review on the latest advancements of earth air heat exchangers[END_REF] and the soil nature provided by [START_REF] Sipio | Factors inuencing the thermal eciency of horizontal ground heat exchangers[END_REF]. Other parameters are considered negligible. The reference values are used for these parameters for the rest of the study. In this case, the boundary values are summarized in table 3 and 3: Range of values for the design parameters and for the nature of the soil characterized by its thermal diusivity.

Three types of ground [START_REF] Thiers | Thermal and environmental assessment of a passive building equipped with an earth-to-air heat exchanger in France[END_REF] are commonly present when installing an EAHE. These values are selected for the study, see table 4.

Soil nature

I ground λ ground c p ground ρ ground α ground Units - W.m -1 .K -1 J.kg -1 .K -1 kg.m -3 * 10 -6 m 2 .s -

Multi-criteria optimization

The selected parameters for optimization are: tube radius r, tube length L, burial depth z, air velocity v air and three soil types characterized by the ground index I ground . These variables vary independently over the range shown in tables 3 and 4. This optimization procedure is carried out for each climatic zone with dierent boundary conditions [START_REF]Les données météorologiques RT[END_REF]. The optimization results are described using two formats: Pareto front and mapping. Due to the large number of graphs generated, this section only presents the results for zone H1a. Other areas have the same trend.

Figure 8 illustrates the Pareto front, which is a function of each objective. The best solutions vary from 0.12 to 0.26 e/kW h for the CER, 2.9 to 14.8 for the COP and 0.90 to 0.98 for the η EAHE . The criteria do not evolve in the same trend because if the cost of the energy recovered decreases, the coecient of performance and the EAHE eciency will deteriorate, and inversely.

Figure 8a refers to the ground index. Two forms can be distinguished: blue points correspond to a clayey-sand silt (I ground = 2) and red points to sand (I ground = 3). The CER is small with silty soil and the η EAHE is strong with sand. The solutions obtained with clay soil are dominated, which is why no points are in gure. Figure 8b shows the best solutions are obtained for tube length varying from 90 to 100 meters. When the tube length increases, the CER and the COP decrease because pressure losses increase. Indeed, when L increases, the η EAHE increases due to a better heat exchange between the ground and the air. Figure 8c exhibits the Pareto front for tube radius between 0.12 and 0.20 meters. When tube radius increases, the COP , the CER and the η EAHE decrease. Indeed, when r increases, the exchanges between soil and air are lower quality and pressure losses are smaller. Figure 8d shows the dominant solutions for air velocity varying from 0.6 to 1.6 m.s -1 . When air velocity increases, the COP and the CER decrease, and the η EAHE increases. The same observation is obtained for tube length. Figure 8e exhibits the good compromises for a burial depth ranging from 2.4 to 4 meters. The burial depth varies according to the type of soil: between 2.4 and 3.4 meters for the sand and 3 to 4 meters for the clay-sand silt. For any type of soil, the CER decreases for decreasing z, this result is consistent since the system cost depends on its depth. Mapping allow exhibiting objectives according to parameters for all solutions obtained during the process. Each point corresponds to one simulation.

Multiple-criteria decision-making

The optimization method made possible to obtain the Pareto front. In order to install this system, it is essential to select a single solution for the optimal pipe conguration of the EAHE. For these reasons, decision aid methods are used.

Table 5 summarizes the best solutions for dierent French climates. The optimal values are between 0.17 and 0.23 e/kW h for the CER, 11.4 to 14.4 for the COP and between 0.96 and 0.97 for the η EAHE . Furthermore, the ground index corresponds to sand (I ground = 3). The value of parameters varies between 14 and 16 cm for r, from 97 to 98 meters for L, from 2.9 to 3.3 meters for z, from 0.7 to 0.8 m.s -1 for v air . As a result, the disparity on values is small.

The cost of the energy recovered is optimal when the total cost is minimal and when the energy recovered is maximum. Indeed, a small air velocity and a small burial depth will generate little cost. In addition, a large tube length and a small tube radius will increase the dierence between the output and input air temperature of the EAHE. It is necessary to obtain the strongest possible temperature dierence to maximize the energy recovered. Then, the EAHE eciency is stronger when the outlet air temperature is very close to the ground temperature. The latter is obtained where the tube length is very large and the tube radius is small. Finally, the coecient of performance is maximum when the energy recovered is maximum and when the energy consumed is minimal. For the same reasons, the energy recovered is maximum for a large tube length and a small tube radius. Moreover, a small air velocity will cause a minimal energy consumption.

The mean and the relative average error are dened by:

x = 1 5 Hi x Hi σ = |x Hi -x| x (17)
with x Hi the value of parameters or objectives for the climate zone Hi.

Therefore, to achieve optimal pipe conguration, the EAHE must be buried in sandy soil at a burial depth of 3. for the CER, the COP , the η EAHE , respectively. The total cost of the system is about 11400 e over 20 years including 58 % for the installation cost, 36 % for the maintenance cost and 5 % for the operation cost. The energy recovered by the system is approximately 2920 kW h and the energy consumed by fans is around 240 kW h. The relative average error on the geometric parameters are less than 10 % so the EAHE can achieve strong energy performance in any French climate.

Conclusions

In this study, a methodology for optimizing the design parameters of an EAHE was carried out for dierent French climates. The proposed modeling of the system was divided into two models, one for soil and one for the EAHE. The model was validated using experimental measurements from the Robert Schuman Institute of Technology (University of Strasbourg). A sensitivity analysis was carried out on the design parameters of the EAHE according to two energy criteria (the average eciency and the coecient of performance) and an economic criterion (the cost of the energy recovered). The results show that tube radius, tube length, burial depth, ventilated air velocity and soil nature are the most impactful parameters. These ones were selected for multi-criteria optimization. Multi-criteria optimization was performed on the EAHE system. The criteria do not evolve in the same trend. Indeed, if the cost of the energy recovered decreases, then the coecient of performance and the EAHE eciency will deteriorate, and inversely. Pareto front solutions give a small tube radius, a large length, an intermediate burial depth and a small air velocity. Finally, a multi-criteria decision-making method was applied to determine the optimal pipe conguration. The EAHE should be buried in a sandy ground, at a burial depth of 3.2 meters, a tube radius of 15 cm, a tube length of 97 meters, a ventilated air velocity of 0.7 m.s -1 . In this case, the value of objectives is 0.20 e/kW h, 12.6, 96.4 % for the CER, the COP , the η EAHE , respectively. In addition, the values obtained are of the same order of magnitude for each climate zone. In complement to this study, it would be interesting to study the sensitivity and optimization of an EAHE associated with controlled mechanical ventilation with double ow and a building. This kind of system is able to renew air inside the home and could satisfy the demand for heating and cooling. Finally, this methodology can be applied for any climate. Similar results can be expected for comparable climates but the results will not necessarily be the same for contrasting climates.

Figure 3 :

 3 Figure 3: Structure of the ventilation system integrated into the building.

Figure 4 :

 4 Figure 4: Schematization of the mesh size of an EAHE.

  air v air 2r µ air P r = µ air c p air λ air the Reynolds number and the Prandtl number. µ air is the air dynamic viscosity and λ air is the air thermal conductivity.

Figure 5 :Figure 6 :

 56 Figure 5: Schematization of the EAHE installed on the IUT Robert Schuman geothermal platform.

Figure 7

 7 Figure 7 illustrates the contribution percentages P C. The rst graph in gure 7 shows contribution percentages of the variables for the η EAHE . Tube radius, tube length and air velocity are the three most impactful parameters because they correspond respectively to 53 %, 42 % and 2 % of the general variation of the η EAHE . However, r, L and v air have a small contribution for dierent β.

Figure 7 :

 7 Figure 7: Percentage contribution of variables for the dierent objectives.

  correspond to the design boundary values of an EAHE for a single-family house of 450 m 3 [12]. Air renewal in the house requires a minimum volume ow rate of 150 m 3 .h -1 . A constraint on the volume ow rate at the outlet of the EAHE ensures this minimum air change.

Figure 8 :

 8 Figure 8: Pareto front for the dierent parameters.

Figure 9

 9 Figure 9 shows the EAHE eciency as a function of air velocity v air and tube length L (left), and as a function of tube radius r and soil index I ground (right). The EAHE eciency values are included between 0 and 0.98. The results show that the η EAHE increases when L and I ground increases, and when v air and r decrease. The points on the Pareto front, illustrate by the pink dots, are located for a large length, small air velocity and tube radius, and for soil index I ground = 2, 3.Figure 9

Figure 9 on

 9 the right shows that there is indeed no solution for a clay soil (I ground = 1).

Figure 9 :

 9 Figure 9: Evolution of the EAHE eciency η EAHE according to the air velocity and burial depth (left) and according to the tube radius and the ground index (right). The Pareto front is constituted by the pink dots.

Figure 10

 10 Figure10shows the COP for burial depth z, air velocity v air and for tube length L, tube radius r. The coecient of performance values are between 0 and 14.8. The coecient of performance increases when v air and r decrease, and when z and L increase. In addition, the COP achieves an optimal value for a burial depth, an air velocity, a length and a radius equal to 2.3 m, 0.8 m.s -1 , 100 m, 0.1 m, respectively. The Pareto front is located for small air velocity and tube radius, intermediate burial depth and large length.

Figure 10 :

 10 Figure 10: Evolution of the COP according to air velocity and burial depth (left) and according to the length and radius of the tube (right). The Pareto front is constituted by the pink dots.

Figure 11

 11 Figure 11 illustrates the CER as a function of air velocity and tube radius. According to this gure, the cost of the energy recovered can vary between 0.12 and 2000 e/kW h. Moreover, air 355

Figure 11 :

 11 Figure 11: Evolution of the CER according to the air velocity and the tube radius (left). Zoom of the Pareto front (right). The Pareto front is constituted by the pink dots.

Table 1 :

 1 10 -6 * 10 -5 Amplitude, frequency and phase shift values for the dierent French climates.

	H1a 12.1 7.19 0.38 3.04 1.99	2.02	7.27 10.7 3.53 10.2
	H1b	11.8 8.44 -0.80 3.34 1.99	2.02	7.27 10.7 3.53 10.4
	H2b	12.3 6.61 -0.54 3.35 1.99	2.02	7.27 10.6 3.94 10.3
	H2c	13.5 7.68 -0.61 4.13 1.99	2.02	7.27 10.6 4.77 10.2
	H3	16.0 8.01 -0.34 2.89 1.99	2.02	7.27 10.5 -0.62 10.5

Table 2 :

 2 tube parameters are characterized by a: radius r, length L, thickness e, rugosity tube , thermal conductivity λ tube , mass thermal capacity at constant pressure c tube and density ρ tube . Burial depth z, air velocity v air , ground thermal conductivity λ ground , mass thermal capacity at constant pressure of the ground c p ground and ground density ρ ground are also studied. Finally, the fan eciency η f an and the electricity price price electricity are taken into account. Parameters considered are summarized in table 2. The reference value corresponds to the experimental data. EAHE design parameters.

	Parameters	References	Units
	r	0.1	m
	L	17.5	m
	e	0.005	m
	tube	0.015	mm
	z	1.03	m
	v air	1.0	m.s -1
	λ tube	0.6	W.m.K -1
	c tube ρ tube	2097600	J.m -3 .K -1
	λ ground	2.7	W.m.K -1
	c p ground ρ ground	1210000	J.m -3 .K -1
	η f an	0.63	-
	price electricity	0.12	e/kW h

Table 4 :

 4 Thermophysical properties for the three types of soil.

	1

Table 5 :

 5 2 meters, a tube radius of about 15 cm and a tube length of 97 meters. The Results of the multi-criteria optimization methodology for dierent French climates. air velocity should be xed at 0.7 m.s -1 , which is equivalent to a volume ow of 175 m 3 .h -1 ; the minimum volume ow constraint is thus respected. The criteria values are 0.20 e/kW h, 12.6, 0.96

	/	Objectives				Variables
	Name	CER	COP η EAHE I ground	r	L	z	v air
	Units e/kW h	-	-	-	m	m	m m.s -1
	H1a	0.21	11.4	0.96	sand 0.15 96.9 2.91 0.72
	H1b	0.17	13.4	0.96	sand 0.16 97.6 3.22 0.67
	H2b	0.23	12.0	0.97	sand 0.14 97.5 3.26 0.71
	H2c	0.20	14.4	0.97	sand 0.14 96.6 3.05 0.69
	H3	0.19	11.8	0.96	sand 0.15 98.2 3.41 0.76
	x	0.20	12.6	0.96	/	0.15 97.4 3.17 0.71
	σ [%]	8.0	8.3	0.3	/	4.3	0.5	4.8	3.4

Acknowledgement

The authors express their gratitude to M. Siroux and C. Fond for the use of experimental data from the geothermal platform of the University of Strasbourg.