M H Benzaama

L H Rajaoarisoa
email: lala.rajaoarisoa@imt-lille-douai.fr

B Ajib

S Lecoeuche
email: stephane.lecoeuche@imt-lille-douai.fr

A data-driven methodology to predict thermal behaviour of residential buildings using piecewise linear models

Keywords: Thermal behavior, Smart-metering, Data-driven models, PWARX, experiments on real buildings

Nowadays, data-driven approaches are a good way to estimate very efficient black box models for different engineering systems. This class of model is well recognized by its outstanding performances to describe the overall behavior of one system based on its input-output relationships without any physical knowledge. In the context of building modeling, this approach is particularly well suited to predict future temperatures or energy consumption in a building. This paper presents an innovative method that uses input-output data to establish reliable and suitable thermal behavior models for residential buildings, especially for existing buildings where only measurements are available and no numerical models are at the disposal of the facility managers. The main paper contributions consist in the design of a new methodology based on the adaptation of a switched model estimation technique and in its validation to model accurately building thermal behaviors. The paper describes different stages needed to reproduce faithfully complex behaviors : data collection, PieceWise affine Auto-Regressive eXogenous (PWARX) identification technique, sensitivity analysis... It also explains how the procedure and the data-driven estimation algorithm are efficient in extracting submodel parameters and sequence that give an outstanding ability to reproduce thermal dynamics of buildings, requiring the only collection of available data. The effectiveness of our methodology is discussed through experiments on different buildings located in the North of France. Indeed, through a comparative study between the piecewise ARX model and other existing models such as nonlinear ARX, indexed ARX and ARX models, the PWARX model gives good results in terms of indoor temperature estimation with 78.48% accuracy.

Introduction

Control of energy consumption continues to be the primary concern in all areas of research [9; 11; 16]. Indeed, the optimization of energy consumption allows for improving the energy performance of the building. Furthermore, to design a set of optimal control, a thermal dynamics modeling step is necessary [START_REF] Hubert | Rules for modeling energy systems in low consumption buildings[END_REF]. However, this task is complicated due to various factors influencing the thermal behavior, particularly by [START_REF] Hiroshi | Total energy use in buildingsanalysis and evaluation methods[END_REF]: (i) climate, (ii) building envelope, (iii) building services and energy systems, (iv) building operation and maintenance, (v) occupant activities and behavior and (vi) indoor environmental quality provided.

Having an accurate thermal model plays a vital role in improving prediction and evaluating energy performance. In the literature, three main categories of modeling approaches have been considered [START_REF] Royer | A procedure for modeling buildings and their thermal zones usingco-simulation and system identification[END_REF]:

• white box models which are based on physical knowledge of the system and thermal balance equations. These are often obtained through energy simulation software like EnergyPlus [START_REF] Al-Janabi | Comparison of energyplus and ies to model a complex university building using three scenarios: Free-floating, ideal air load system, and detailed[END_REF], TRNSYS [START_REF] Benzaama | Investigation of the thermal behavior of a combined geothermal system for cooling with regards to algeria's climate[END_REF], etc;

• black box models which use only measured input/output data and statistical estimation methods (e.g. [24; 26]);

• grey box models, a mix of the first two categories above. They use input/output data as well as some a priori knowledge on the system. A popular grey-box model is the equivalent RC networks [24; 34; 36].

A comparative study between these different models was done in [START_REF] Abdul | Development and performance comparison of low-order black-box models for a residential hvac system[END_REF]. The main conclusions of this study are that: (i) the use of white box model often requires important set-up and computation time; (ii) it also involves a large number of inputs to define the model, such as the composition of the building envelope for example [START_REF] Crawley | Contrasting the capabilities of building energy performance simulation programs[END_REF]. In some studies, it is difficult, if not impossible, to recover this input [START_REF] Rick | Inverse modeling of simplified hygrothermal building models to predict and characterize indoor climates[END_REF]. To overcome this problem, data-driven methods have emerged in building framework. We find from the literature that the most commonly used data-driven techniques for building thermal modeling and energy performances prediction are based on the ARX (AutoRegressive eXogeneous) model [21; 31]. However, despite its performances, such as quick implementation, a good accuracy..., this model has important limits, mainly due to the estimation of a unique thermal model. A unique model may not consider dynamic changes due to usages, equipment configurations or external factors, such as wind and solar radiation, that influence the thermal dynamics of the building. Today, we can find several data-driven models derived from the ARX model. One can read the following references for further information [2; 15; 17; 22; 27; 29; 37]. Each of them is mainly differentiated by the parameter estimation techniques and the structure of the model. For instance, in [START_REF] Lin | Fractional order models for system identification of thermal dynamics of buildings[END_REF], the authors tested the Fractional order Auto-Regression with eXogenous variable (FARX) model on building integrated energy systems. The model has been validated using the input-output data retrieve of a residential building simulated with software IES<VE>. As a result, it has been shown that the FARX model gives better accuracy than ARX one. Also, the results in [START_REF] Mechaqrane | A comparison of linear and neural network arx models applied to a prediction of the indoor temperature of a building[END_REF] and [START_REF] Soleimani-Mohseni | Estimation of operative temperature in buildings using artificial neural networks[END_REF] show that the Non-linear ARX (NARX) model performance was significantly greater than the one of ARX model. However, each technique has its own advantages and inconveniences, and one of our motivation is to discuss on how to recover the best data-driven model structure for simulating the thermal behavior of the building. Issues that one may remark into the following papers [33; 37] and [1; 22].The first response to this discussion has been presented in [START_REF] Guo | Hourly cooling load forecasting using time-indexed arx models with two-stage weighted least squares regression[END_REF] where the best prediction method includes a combination of two separate time-indexed ARX models to improve the prediction accuracy of the cooling load over different forecasting periods. In the following, we will also contribute to this discussion by making a comparison of existing ARX-derived techniques.

Paper contribution

In this paper, a black-box identification approach by means of a PieceWise affine auto-Regressive eXogenous (PWARX) model will be developed in the framework of building thermal modeling. This new methodology takes profit from recent advances in the hybrid system identification community. Hybrid and PWARX systems are heterogeneous dynamic systems that combine simultaneously continuous and discrete dynamics. These systems are helpful to introduce expert knowledge in the data-driven models, especially when various behaviors or uses have to be explained. They can be represented by switching models, i.e. by a set of continuous-submodels indexed by a discrete mode or a specific building management system setting. However, the idea of using partially observed regime switching models for building thermal modeling has not received satisfactory solutions apart from a few attempts which do not combine real world learning applications such as: missing data, dependent time series and noisy observations, see [14; 28]. Also, these works mainly highlight the necessity of using different models to represent different sorts of dynamics in a building. It was argued too in [START_REF] Ajib | Data-driven building thermal modeling using system identification for hybrid systems[END_REF] and has been proved in the context of the prediction of buildings energy consumption into [START_REF] Ajib | Prediction of standardized energy consumption of existing buildings based on hybrid systems modeling and control[END_REF]. Then, the purpose of this paper will be to highlight the interest of using PWARX systems to model building thermal dynamics, in its 3 aspects: presentation of its scientific foundations, with regard to conventional ARX techniques; discussion around its interest to give explanations on different thermal dynamics; experimental comparison of the results obtained with different state-of-the-art methods.

Besides, this paper also discusses the selection of suitable inputs to define the best model structure. This will be done by performing sensitivity analysis and testing several configurations. This step aims to state on the influence of each input on the accuracy of the model and on the quality of identified parameters. Thus, the PWARX model should be able to explain the true thermal and energy behaviors of the building by identifying the use scenarios. Finally, the effectiveness of our methodology will be shown by presenting thermal modeling results for different building architecture located in the north of France. In particular, we assess the performance of our model through a comparative study between the piecewise ARX model and other existing models such as nonlinear ARX, indexed ARX and ARX models.

Paper outline

The rest of the paper is organized as follows. The formal definitions of the ARX model and its derivatives, as well as PWARX model are introduced in Section 2. The system identification technique used to identify and validate the PWARX model is detailed in Section 3. To illustrate the effectiveness of the proposed method, experimental results obtained from a student residential are provided in Section 4. Finally, conclusions and further works are given in Section 5.

Background

ARX model

Let us first consider an ARX model, using an input/output (y, u) representation, writing as follows:

y(t) = -a 1 y(t-1)-• • •-a n a y(t-n a )+b 1 u(t-1)+• • •+b n b u(t-n b )+e(t)
(1) where n a and n b are the model order, a i and b i are the model coefficients, and e(t) ∈ R n e is a white noise process. In other terms, the ARX model can be defined by the following relations:

A(z)y(t) = B(z)u(t) + e(t) (2)
with

A(z) = 1 + a 1 z -1 + • • • + a n a z -n a (3)
and

B(z) = b 1 z -1 + • • • + b n b z -n b ( 4
)
where z is a backward shift operator. So, for available input measurements, we can estimate the output (ŷ) at each time t by the ARX model as follows:

ŷ(t) = B(z)u(t) + e(t) + (1 -A(z))y(t). (5)
Finally, in compact form the estimate of the output can be writing as:

ŷ(t) = ϕ T (t)θ (6)
with

ϕ(t) = [y(t -1) • • • y(t -n a ) u(t -1) • • • u(t -n b )] T (7)
and

θ = [a 1 • • • a n a b 1 • • • b n b ] T . (8)
In these equations, ϕ represents the regression vector and θ the parameter vector. In the context of the design of a building thermal model, y generally represents the indoor temperature, and u represents the measured factors that influence the temperature evolution.

Taking the thermal ARX model example in [START_REF] Amasyali | A review of data driven building energy consumption prediction studies[END_REF], the indoor thermal behavior (T i ) depends on the following inputs: the outdoor air temperature (T o ), solar gains (Ra), internal gains (Q u ) and thermal gains (Q h ). Thus, the relationship between the input-output is expressed by:

T i (t) = - n a j=1 a j T i (t -j) + n b r=1 (b 1,r T o (t -r) + b 2,r Ra(t -r) + b 3,r Q h (t -r) + b 4,r Q u (t -r) + v(t)). (9)

NARX model

NARX is a nonlinear auto-regressive network with exogenous inputs. It has a recurrent dynamic nature and is commonly used in time-series modeling [START_REF] Asgar | Narx models for simulation of the start-up operation of a single shaft gas turbine[END_REF]. The formal definition of the NARX model is given by the following expression [START_REF] Beale | Neural network tool boxtm user's guide[END_REF]: [START_REF] Benzaama | Investigation of the thermal behavior of a combined geothermal system for cooling with regards to algeria's climate[END_REF] where y is the output variable, u is the deterministic externally variable, and f represents the nonlinear function. n a and n b represent the model order. Thus, the next value of the output signal depends on its previous values and those of the input signal (exogenous input) [START_REF] Asgar | Narx models for simulation of the start-up operation of a single shaft gas turbine[END_REF].

y(t) = f [u(t-1) • • • u(t-n b ) y(t-1) • • • y(t-n a )] + e(t)
Basically, a NARX model can be implemented by using a feed-forward neural network to approximate the function f . This network can be used for nonlinear filtering of noisy signals or for predicting future behavior [START_REF] Lin | Learning long-term dependencies in narx recurrent neural networks[END_REF]. Moreover, it converges faster, and is less susceptible to variations and dependencies. It has typically better generalization abilities than other networks [START_REF] Ferracuti | Data-driven models for short-term thermal behaviour prediction in real buildings[END_REF].

As an example, let us consider the NARX structure described in [START_REF] Ferracuti | Data-driven models for short-term thermal behaviour prediction in real buildings[END_REF]. It uses as output the indoor temperature (T i ) and as input variables the outdoor air temperature (T o ), solar gains (Ra), internal gains (Q u ) and thermal gains (Q h ). They are necessary to predict the internal temperature as follows:

T i (t) = f [T i (t -1), • • • , T i (t -n a ); T o (t -n k ), • • • , T o (t -n k -n b ); Ra(t -n k ), • • • , Ra(t -n k -n b ); Q h (t -n k ), • • • , Q h (t -n k -n b ); Q u (t -n k ), • • • , Q u (t -n k -n b )] + v(t), (11)

Indexed ARX model

Several researchers introduce a more efficient ARX model by including physical understanding of the building in the model structure. This kind of model is called indexed ARX model. This model is more accurate than typical regression models and more efficient, in terms of computation time for instance. So, the ARX model is indexed with respect to different time intervals, different modes or configurations and different inputs to improve accuracy and efficiency. The (nq-m)-indexed model is then defined as a ARX model while considering for example n past times for the output (y : {y(t), y(t -1), ..., y(tn)}), a q modes or configurations and m input numbers (u : {u 1 , u 2 , ..., u m }). The input-output relationship can be written by the following expression:

y(t) = -a 1,q y(t -1) -• • • -a n,q y(t -n) + b 1,q u 1 (t) + • • • + b m,q u m (t) + e(t). (12)
Practically, a 4-3-5 ARX indexed model has been defined in [START_REF] Mechaqrane | A comparison of linear and neural network arx models applied to a prediction of the indoor temperature of a building[END_REF] by the following configurations: the load prediction Load as output and temperature T o , humidity H o , solar radiation Ra, wind speed W d and people occupancy O c as inputs, can be written as follows:

Load(t) =b 1,q T o (t) + b 2,q H o (t) + b 3,q W d (t) + b 4,q Ra(t) + b 5,q O c (t)
+ a 1,q Load(t -1) + a 2,q Load(t -2) + a 3,q Load(t -3)

+ a 4,q Load(t -4) + a 5,q .

Remark 1. Here three modes describe respectively : the daytime (8 AM to 9 PM), the transition time (6 AM to 8 AM), and the nighttime (9 PM to 6 AM).

PWARX model

PieceWise affine ARX model is a special class of non-linear models, obtained by splitting the state-input domain into a finite number of polyhedral regions. It can model a large number of physical processes such as residential building systems with several inherent operating modes or configurations. The general structure of a PWARX model is:

y(t) = f (ϕ(t)) + e(t) (14)
with f as a piecewise affine map of the following form:

f (ϕ(t)) =              θ T 1 φ(t) if σ(t) = 1 . . . θ T q φ(t) if σ(t) = q, (15)
where φ = ϕ T 1 T is the extended regression vector. σ(t) is the switching rule defined by:

σ(t) = i iff ϕ(t) ∈ R i , for i = 1, • • • , q; (16)
and {θ i } q i=1 are the parameter vectors that define the sub models. {R i } q i=1 represent a complete partition of the region R ⊂ R n , with n = n e n a + (n b + 1), and each region is a convex polyhedron with:

R i = {ϕ ∈ R n : H i φ 0} (17)
where H i and 0 are respectively a matrix of appropriate dimensions defining the limit of the region partitioning the set of regression vector and the null vector. To complete and give more details about the definition given above, figure 1 shows the partitioning into {R i } 3 i=1 region of the regression space according to scenario described by each parameter vector {θ i } 3 i=1 (see [START_REF] Boukharouba | Identification of piecewise affine systems based on dempster-shafer theory[END_REF] for more details). So PWARX model is defined by a set of ARX sub-models that help to analyze and predict the various thermal behaviors of the building. Each sub-ARX model is assumed to describe one mode or configuration [START_REF] Boukharouba | Identification of piecewise affine systems based on dempster-shafer theory[END_REF].

Basically, the physical presentation of thermal behavior of indoor air by considering the q th mode and with m input variables is defined by:

T i (t) = n a j=1 a j,q T i (t -j) + m v=1 n b j=1 b j,q u v (t -j), (18)
and in terms of outdoor air temperature (T o ), outdoor humidity (H o ), solar radiation (Ra) and heating power (Pw), it can be presented as:

T i (t) = n a j=1 a j,q T j (t -j) + n b r=1 b 1 r,q T 0 (t -r) + b 2 r,q H 0 (t -r) + b 3 r,q Pw(t -r) + b 4 r,q Ra(t -r) . (19)
So, based on the input/output pair generated by the system described by relation [START_REF] Herie | Modeling of a building system and its parameter identification[END_REF], the identification of the PWARX model consists in (i) determining the number of sub-models in order to have the best fitting between the measured and estimated output, (ii) associating the data of each affine sub-model respectively, (iii) and finally, estimating each parameter vector θ i associated with these sub-models. Details about this identification procedure is given in the following section.

PWARX thermal model identification methodology

Iin this paper, we present an original data-based method to identify and estimate a set of models that are able to reproduce the variety of building thermal behaviors (Figure 2). So, to estimate all parameters of the PWARX model, we adopt an identification procedure defined by the following steps. Firstly, we design the experiment and the system necessary to collect the data. Then, we acquire some measurements from different room parts, without requesting specific usages or functioning conditions. Our case study consists of student residential buildings located in Douai (France). Different scenarios will be considered during the measurement campaigns, such as the orientation (rooms situated at East/West), a different heating sequence and various periods of acquisition in the objective to provide sufficient data for discussing on the validation of our new methodology.

Secondly, we proceed to the determination of the model structure by a sensitivity study according to several parameters. Indeed, the sensitivity analysis makes it possible to identify the parameters of influence of the model for the prediction and to choose well (or validate) the model structure for a maximum precision and a minimum computation time in the context of building thermal modeling. The main objective of this study is then to give a good approach to estimate switched thermal dynamics and then to reproduce complex building thermal behaviors based on a data-driven modeling. Details are given in the next paragraphs.

Smart metering

Our objective is to compute a thermal model for a building without any particular prior knowledge on the system. To this end, more extensive and uniform data is required, to provide more general information about the thermal behavior of each building studied. Thus, a smart supervisory, control and data acquisition (2SCADA) system is designed. Figure 3 illustrates its overall architecture. The system is based on wireless sensors (developed by CLEODE TM company) that allow measuring the indoor and outdoor temperature (T i , T o ), relative humidity (H i , H o ), solar radiation and diffuse solar radiation in the building (R a , R ai ).

Figure 4 shows examples of different trends for the indoor and outdoor temperature (on the top), and the indoor and outdoor humidity (on the bottom). Each sensor is defined with a sampling time equals to T s = 2.4min.

Complementary to the smart-metering system, a control system is used to control (ON/OFF) the heating operation for a few scenarios. Table 1 below defines each sensor characteristics used for collecting data. These components are based on wireless system network (WSN) technology which is today the most popular technology used in smart-metering systems. Moreover, many WSN platforms have been developed [START_REF] Badreddine | Optimized energy management for a smart building multi-sources multi-charges: different principles of validations[END_REF], ZigBee/IEEE 802.15.4 protocols are often part of them. These protocols are a global hardware and software standard designed for WSN requiring high reliability, low cost, low power, scalability, and low data rate [START_REF] Markus | Buildings, climate and energy[END_REF].

The data are collected, using the smart-metering and the control systems, in order to analyse and predict the thermal behaviors of the building for several non-defined operating modes.

Estimation of the PWARX model

Here, we justify the choice to model the building thermal behavior as a switching system by its ability to detect changes in thermal dynamics resulting from a given operating mode, in particular an operation associated with the heating system, human activities, weather conditions, etc. Thus, we can implement a set of sub-models that reproduces more faithfully the reality, by comparison with the standard ARX model.

The procedure used to estimate the different parameters of each sub-model is described into four steps. The first one is the initialization step. Knowing that the number of sub-models is a priori not known, this step allows to proceed to the creation of Λ singleton classes with the Λ regression available data. The second step is the re-affectation of data followed by estimating parameters of each sub-model. Let notice that this step consists mainly in grouping the regression data (initially divided into Λ classes) into a minimum number of disjoint classes. The third one is a verification step that is necessary for the stabilization of the data classification. And the last one is the validation phase of the model by applying a structural analysis. Each step is recovered and detailed in the following paragraphs.

A-Initialization. In this part, the number of sub-models is set to max(q) = Λ, i.e. Λ classes are obtained,

C = {C 1 , • • • , C Λ }.
Each class corresponds to one sub-model defined by the parameter vector θ i as shown in figure 1. The initial value of each parameter vector per class is noted by θ0 1 , • • • , θ0 Λ . So, to compute θ0 i , we use the data x(i) of the class C i and its c nearest neighbor (c-nn) and finally apply the least-square estimation technique.

Thus, for an input/output data {(ϕ(t), y(t))} N i t=1 generated by the affine model [START_REF] Herie | Modeling of a building system and its parameter identification[END_REF], the parameter vector θi is computed by least square technique while minimizing:

V(θ i ) = 1 N i N i t=1 (y(t) -θ T i φ(t)) 2 (20)
where N i is the sample number.

B-Data re-affectation and model estimation. The aim is to reduce the number of clusters by reclassifying each data point and to estimate the parameters of each sub-model. To do that, we try to minimize the Euclidean distance (Figure 6) between each pair of data (x(i), x( j)) by:

d j i = x(i) -x( j) (21)
as well as the error between the measured output and the output of each sub-model is minimal. Moreover, data x(i) is ranked according to the regression vector and the measured output as the following expression:

x(i) = [ϕ(i) T , y(i)] T . (22)
Thus, the re-affectation of data x(i) is done relative to their nearest neighbors (c-cnn). Formally, data will migrate towards the most representative clusters according to a specific decision rule:

ψ i j = exp(-α σ( j) x(i) -x( j) 2 -β σ( j) (y(i) -θT σ( j) φ(i))) 2 (23)
with θT σ( j) is the parameter vector corresponding to the class C σ( j) . α σ( j) and β σ(i) are positive parameters. Thereafter, we need to update the class number C = {C 1 , • • • , C q} where q ≤ Λ. And, to decide on the re-affectation of the data x(i) to the remaining class q, we use the membership criterion of x(i) to the class C p , defined by: [START_REF] Liu | The re-optimization strategy of multi-layer hybrid building's cooling and heating load soft sensing technology research based on temperature interval and hierarchical modeling techniques[END_REF] where q determine the persistent number of class in order to have the best fitting between the measured and estimated output. Let notice that P(x(i) [START_REF] Boukharouba | Identification of piecewise affine systems based on dempster-shafer theory[END_REF].

P(x(i) ∈ C p ) = j|x( j)∈{Γ c (x(i)∩C q )} ψ i j c j=1 ψ i j , p ∈ {1, • • • , q}
∈ C p ) = 1 if the elements of Γ(x(i)) ∈ C p and P(x(i) ∈ C p ) = 0 if the element of Γ(x(i)) C p
C-Convergence criterion. The convergence is done by comparing the earlier parameter vectors θ

(r) = θ(r) 1 , • • • , θ(r)
q and the posterior parameter vectors θ (r+1) following:

θ (r+1) -θ (r) ≤ υ ( 25
)
where r is the index of the iteration and υ is an arbitrary threshold defined by the user. Generally, one use υ = 10 -5 .

D-Sensitivity analysis.

A study of the influence of each input onto the model structure is presented in this article by comparing the FIT criterion value for different configurations. The FIT represents the similarity between the measured output y and the output ŷ predicted by the model. It is given by the following equation:

FIT = (1 - ŷ -y y -ȳ ) × 100% (26)
where ȳ and ŷ are respectively the mean and the estimate of the measures y(t). Once a model with a satisfactory FIT has been estimated, it is possible to evaluate the influence of each parameter on the model and on the prediction accuracy. To this end, we analyze the sensitivity of each parameter by identifying and characterizing the: (i) input uncertainties, (ii) model orders, and (iii) class number.

Brief, this analysis allows us to appreciate the contribution of each input variable onto the FIT. This contribution is computed by using the following influence index defined by:

W ui = ρi m i=1 ρi × 100%, ( 27
)
where ρi is the contribution estimated of each input obtained from the computation of F IT k . This last is the estimation of the FIT subject to the influence of inputs considered for the k th model structure.

It can be expressed by:

F IT k = N i t=1 m i=1 ξ i (k) ρi u i (t), for each model structure k (28)
with

ξ i (k) =         
1 if the i th input is considered for the model structure k, 0 otherwise.

(29) Concerning the other terms, u i (t) represents each input variable used for the identification of PWARX model. Model structure k is then the configuration considered, in terms of input number, to estimate thermal model parameters. For instance scenario where only two inputs between m are considered.

Given the influence index W ui for each input signal, we can realize an immediate analysis of how much the uncertainty on the measurements can affect the indoor air temperature prediction. Thus, this procedure should allow to properly structure the model in terms of input data in order to increase the accuracy of the model and reduce the calculation time.

Experimental results and discussion

Case studies

To show the effectiveness of our methodology, experiments on the Lavoisier student residential building, located in Douai, in the north of France (Figure 6) have been made. The total area of the building is approximately 3500m 2 and it is subdivided into 2 sub-building parts. The first part is composed by 4 floors and the second one by 5 floors. Each floor has respectively for each sub-building 10 and 32 rooms of 11m 2 with a double glazed window size (135cm × 110cm). Notice that for this building, the envelop is composed by brick walls with a thickness 22cm and a glass wool insulation with a thickness 14cm.

So, we launched experiments on the second sub-building during this study. Herein, half of the rooms are located in To recover a persistent database for identifying the underlying causal thermal model of the building, several experiment scenarios have been taken into account: two different orientations "East-West", heating sequence and different floor levels (1 st and 3 rd floors). In addition, rooms operate on different thermal excitations driven by a heating system supply and are influenced by natural solar radiation. Precisely, the heating supply will be controlled on two different sequences such as a controlled sequence and a random heating sequence.

Also, the same experimental protocol was renewed for Condorcet student residence located in Douai (Figure 7). The difference of two buildings is mainly the envelop. Indeed, for Condorcet we have a block wall with a thickness 20cm combined with polystyrene insulation with thickness 10cm. These two buildings are considered with the aim of validating the model for its ability to reproduce the building thermal behaviors faithfully only using a data collection.

So, for each scenario we evaluate the thermal environment thanks to different sensors and actuators placed in each room as it is illustrated onto Figure 8 above (example for west oriented room). Once again, readers can refer to Table 1 for sensor and actuator characteristics. Moreover, the location of the temperature sensors makes it possible to measure the indoor air temperature of the irradiated zone (T i1 ) and the shaded zone (T i2 ) for the two orientations (i.e. East-West), that we add the solar radiation sensor (noted Lux ext and Lux int on the figure), the heating power consumption measure (noted Ch 1 and Ch 2 on the figure), the outdoor humidity and temperature (H e and T e ).

Besides, during the measurement campaigns each collected data associated with each scenario is saved in a database in CSV format. This is possible thanks to an in-house developed application based on LABVIEW TM and NetBeans TM software's. The historical data is used after that to identify the PWARX-thermal model according to the experimental planning detailed after.

Experimental design

During this experience we would firstly like to compare two rooms of the Lavoisier residence according to its orientation (i.e. east and west). This stage aims to find a model that validates the indoor thermal behavior, for these two orientations under the same conditions.

East side

The experimental study was made from February 07 th -12 th February, 2018. The evolution of the outdoor temperature varies between 0 and 10 • C for an outdoor humidity ranging between 70 and 100% as indicated respectively in Figure 9a and 9b. Moreover, Figure 10a and 10b show the evolution of indoor air temperature (black solid line) for two different sensor positions (T i1 and T i2 ) and the heating power demand (black dot line), during the period defined before. The temperature varies between 19 and 25 • C and the electrical power consumed between 0 and 1400W.

In the other hand, the solar radiation received on the east wall and the radiation incident on the room is shown in figure 11a and11b. We can see that solar radiation received by the wall reaches a maximum value of 11.5W/m 2 , while the incident solar radiation on the floor is approximately 0.45W/m 2 . And finally, two scenarios are considered in this study: the first one is the free evolution scenario (i.e. without heating supply), and the second one consists to apply a random activation of the heating source. This is renewed also for the west side.

West side

For the same period, i.e. from 07-02-2018 to 12-02-2018, the solar radiation received on the west wall reaches a maximum value of 95W/m 2 while the solar radiation received on the floor takes a maximum value of 42.5W/m 2 (see Figure 12). Moreover, Figure 13 below shows that the indoor air temperature of the shaded area varies between 19 and 26 • C, however for the irradiated area, the indoor temperature varies between 19 and 31 • C, which generates a difference of 7 • C between the two positions.

Analysis for different heating sequences

In this part, experiments are made from 13/02 to 17/02, with the following conditions: (i) two west oriented rooms respectively at the 1 st and 3 rd floor are considered, (ii) weather conditions are the same in terms of outside temperature, outdoor As example, we can observe in Figure 14, the west side receives a solar radiation power which rises to 90W/m 2 . On the other hand, to meet the third condition cited before, the 3 rd floor heating system is controlled in a random manner, which causes a random variation of the indoor thermal behavior (Figure 16b). While, the room on the 1 st floor evolves with two different scenarios: free evolution (i.e. heating system OFF) and heating system activated (Figure 15b). shown in Figure 16a.

PWARX-thermal model validation

The objective of this part is to propose a model validation process in order to present the precision of estimated parameters and to evaluate how the choice of hyper-parameters (that define the model structure) could influence the results. For example, in [START_REF] Mechaqrane | A comparison of linear and neural network arx models applied to a prediction of the indoor temperature of a building[END_REF], authors have taken only the outside temperature as inputs using an order number equals to n a = n b = 4. In contrast, authors of [START_REF] Sarwar | Field validation study of a time and temperature indexed autoregressive with exogenous (arx) model for building thermal load prediction[END_REF] used the outdoor temperature, outdoor humidity, solar radiation and wind speed as inputs for the same order number n a = n b = 4. However, the choice of initial model structure hyper-parameters was not justified in most of related works. Approach for analyzing and validating the structure of the model will be presented in the next paragraphs.

Sensitivity analysis

Sensitivity analysis is used to determine which is the best structure for the model to faithfully reproduce the thermal behavior on a variety of measurement campaigns. Thus, the structural analysis serves to regulate the algorithm by making an analysis on the hyper-parameters such as n a , n b , class number (C) and the input choices.

The determination of the inputs and their influence on the structure of the model are first presented in this section. To this end, we estimate a model derived by following configurations:

• n a = 2, n b = 2, n k = 0, C = 200, β 0 = 0.5, γ 0 = 0.5.
• Inputs measured: Outdoor temperature (T o ), outdoor humidity (H o ), solar radiation (Ra), heating power (Pw).

• Output measured: Indoor Temperature T i .

By varying the system order n a and n b in Table 2, we can see several model structures (1 to 4) and the FIT values obtained with each of them. We observe that the best estimation is reached when n a and n b correspond to an order of 5 with a FIT equals to 78.48%. this model (resp. accuracy of the parameters) has been also investigated. Results in the bottom part of Table 2 show the impact of inputs on the FIT's value. We can note that the FIT decreases or increases according to the inputs selected. The worst model corresponds to a structure composed only of T o and H o (FIT=15.27%). By adding the solar radiation to this structure, the FIT increases and becomes 24.61%. The model structure with Ra and Pw as inputs gives a better FIT compared to the model structure with T o and H o with a difference of 40.41%. By adding the outdoor humidity to this structure, the FIT increases and becomes 70.38%. To summarize, Figure 17 represents the effect of the system orders n a and n b (varying from 2 to 5) and c-nearest-neighbors c (varying from 100 to 300). The results show that, with all inputs, the best model corresponds to n a = n b = 5 and c = 200 neighbors.

On the other hand, Figure 18 and 19 show a comparison between the model structures " 4 D , 4 A and 4 C " as defined in table 2 (indoor estimated temperature by red solid line, indoor measured temperature by black solid line). Results indicate that solar radiation and heating power are essential in order to estimate models with accuracy and precision. The absence of the heating power as input (Model 4 C ) reduces the FIT by 53.87% compared to model 4 (Figure 18b), while the absence of the solar radiation and heating power (Model 4 A ) reduces the FIT by 63.21% compared to model 4 (Figure 18a). We have also noted that the absence of the outdoor temperature (Model 4 D ) reduces the FIT by 8.1% compared to model 4 as shown in Figure 18b.

Figure 20 shows the influence index of each input variable on the model structure. It is mainly based on the sensitivity analysis conducted before and computed in order to classify the set of inputs (see equation ( 27)). The results show that the heating power (Pw) is the most influential input with a relative influence index of 70%. This is due particularly to the impact of the heating power on the indoor thermal behavior. For the next, we remark that solar radiation (Ra) affects with an influence index of 20%. Let note that the thermal response time of solar radiation is greater than of heating power. Then, the indoor thermal behavior is influenced by the diffuse solar radiation which depends on the building orientation, glazing type and weather conditions. At the end, outside humidity (H o ) and temperature (T o ) are less influential, with computed indexes of 8% and 4% respectively.

To conclude this analysis, modeling in efficient way the thermal behavior of the building requires at least the heating power and the solar radiation as inputs of the model. The other terms can be used simply to add further information to the model.

PWARX and ARX derived models

This section presents the validation of the PWARX-thermal model by comparing it with ARX and its derived models. We choose the PWARX structure model equivalent to the model type number 4 (see Table 2). Moreover, for the building investigated, the algorithm estimates 3 operating modes describing the thermal behavior of the building. Figure 21 shows the estimated switching evolution of these operating modes (black dashed line) based on data collected from the eastfacing room.

Precisely, three operating modes (i.e. discrete state) were detected with respect solar radiation (Ra) and the use of the heating system (Pw); which is simply the logical consequence of results obtained during the sensitivity analysis above. For instance, the first sub-model < S M1 > (resp. first operating mode) corresponds to the activation of the heating system (ON). The second sub-model < S M2 > (resp. second operating mode) defines the situation where heating system is (OFF) and solar radiation is greater than 2.19W/m 2 , while the third sub-model < S M3 > (resp. third operating mode) corresponds to the situation where heating system is (OFF) and solar radiation less than 2.19W/m 2 . This is summarized in table 3 below.

Modes

Pw Ra Mode 1 ON -W/m 2 Mode 2 OFF ¿2.9 W/m 2 Mode 3 OFF ¡2.9 W/m 2 Table 3: Switching conditions Furthermore, we can see on Figure 21, the validation of these switching conditions by comparing the estimated operating modes (black dashed line) with the predicted (red dot line) by the rules defined in table 3. And finally, each vector of parameters associated with each operating mode is given in Table 4, as numerical application example for the PWARXthermal model identified. Basically, for the q th mode and with m input variables, we have: resentation for mode q = 1 (resp. S M1): We remark here, on the one hand, that the building thermal behavior is not well recovered by the ARX model, and on the other hand, the operating points defined for the indexed-ARX model described in [START_REF] Mechaqrane | A comparison of linear and neural network arx models applied to a prediction of the indoor temperature of a building[END_REF] are not really suitable to explain building configurations.

T i (t) = n a j=1 a j,q T i (t -j) + m v=1 n b j=1 b j,q u v (t -j). ( 30
)
T i (t) = -1.1562T i (t -1) + 0.1606T i (t -2) -0.0031T o (t -1) + 0.0031T o (t -2) -0.0017H o (t -1) + 0.0019H o (t -2) + 2.508e-06Ra(t -1) -8.99e-07Ra(t -2) + 3.32e-05Pw(t -1) + 4.15e-05Pw(t -2) + 0.03 (31)
So, by confronting with PWARX model, the NARX model can give good results (Figure 22 -top and left part) in terms of indoor temperature estimation (pink solid line), compared with the measurements (black solid line). However it remains a challenging task to interpret without proper knowledge of the internal structure of the model. In other terms, it is difficult for building experts for instance to understand physical phenomena, especially when various behaviors or uses have to be explained. As we can see, PWARX-thermal model makes possible the understanding of the system behavior in a much more intuitive and easy way than other models.

Finally, Table 5 gives a numerical quantification of the performance of each model based on the FIT criteria. The results confirm that the PWARX model structure is the best for prediction purposes, being able to estimate accurately the real temperature according to different scenarios or configurations, influencing the thermal behavior of the building. February, 2018, for the east orientation and using training data as presented in section 4.1. In this section, we want to validate our model with other scenarios, firstly by considering data collected from west oriented room for the period February 14 th -17 th February, 2018. Figure 23 shows the indoor temperature estimation results (red solid line) with a FIT = 60.02%. Thus, we observe with satisfaction that the structure of the estimated model is adapted to east and west orientations.

Model

In the next, the same PWARX-thermal model was tested with data collected from the room on the 3 rd floor and located in the west-facing, from February 14 th to 17 th , 2018. The thermal behavior of this room is mostly different compared to the reference room one due to the orientation, floor level, heating scenarios and weather conditions. Once again the result in Figure 24 shows there is good agreement between the experimental measure (black solid line) and the prediction (red solid line) with a FIT = 77.66%.

PWARX-thermal building methodology tested on an-

other building Authors of [START_REF] Herie | Modeling of a building system and its parameter identification[END_REF] is to show that the reference model structure may also be suitable for any building's archetype. Besides, it is important to note that the thermal behavior conditions are not also similar due to different heating scenarios and weather conditions. So, Figure 25 below gives us any ideas on indoor temperature dynamics. Indeed, herein the indoor thermal behavior, respectively for rooms in east and west side, can be estimated by the PWARX model with a FIT of 57.40% and 77.66% for both rooms. These results are acceptable and in line with configurations given before.

Comparison between models and evaluation of thermal behavior similarities

Figure 26 below shows the comparison procedure that could be implemented in order to compare the thermal behaviors of different rooms or buildings, using our PWARX methodology, in general manner. In fact, this procedure consists in testing a reference PWARX-thermal model on other measurements derived for other rooms to be compared. In this example, five rooms belonging to both the Condorcet and Lavoisier residences are used, one as a referenced model (Room 1) and 4 others to be compared. The objective is to analyse the ability of the reference model to recover another thermal behavior identified, based on a particular case study (examples: heating mode, orientation, weather conditions), but also to detect possible true different thermal dynamics. Briefly, it will allow us to verify the similarity of the thermal behavior of each room with respect to the reference room from Lavoisier building (noted here Room 1). The FIT's difference between the Room 1 and the other rooms will be considered as a similarity rate (S R) and makes it possible to quantify the thermal behavior discrepancy. This last can be computed by the following expression: S R = 100 -(FIT room 1 -FIT room i ), with i = 2, ..., 5. [START_REF] Royer | A procedure for modeling buildings and their thermal zones usingco-simulation and system identification[END_REF] where FIT room i is computed using the measurements collected on room i and the PWARX model estimated on room 1. The FIT results in Table 6 show that Rooms 1 and 5 probably have the same behavior (similarity rate equal to 99.18% as shown in Figure 27). Of course, both rooms are located in two different residences (Condorcet, Lavoisier) but this can be explained by the fact they are in the same orientation and level floor, and have also the same heating equipment. Furthermore, the comparison between Room 1 and Room 2 shows the floor level difference impacts on the thermal behavior of the rooms. For this case, the similarity of the thermal behavior of both is around 81.54%, as illustrated in Figure 27. Now, if we consider the thermal behavior of the Rooms 1 and 4, it is not really similar related to the orientation. Indeed, the data collected from Room 4 have lower solar energy values than those from Room 1. As a result, the similarity rate decreases significantly and reaches 78.92% (see Figure 27). Finally, Room 3 can be considered as the worst case for this validation. In fact, as the heating operating equipment and floor level are different from the initial scenario, the similarity rate obtained is around 70.67% compared to the reference Room 1. Overall proportionality of the similarity rate is illustrated in Figure 27. In summary, experimental data allows to interpret easily the thermal behavior of different rooms (resp. buildings). In other terms, it is possible to compare each room (resp. building) to another reference one by deducing in particular the orientation and/or the floor level are/is different or similar of the reference for instance. On the other hand, we can use the similarity rate value to state if a reference model is suitable to also explain the thermal behavior of any building's archetype, even if the configurations are not identically exact due to different heating scenarios and/or weather conditions.

Conclusions

In this paper, we present a novel methodology for modeling the thermal behavior of residential buildings by using the PWARX model. The identification procedure is detailed to obtain the model and its parameters. Also, several analyses have been conducted to have a suitable structure for the model, as well as to validate it. The first of them is the sensitivity analysis which is necessary to determine what inputs will be required in order to guarantee the accuracy of the model and the quality of identified parameters. We can conclude that solar radiation and heating power are the main influential inputs for modeling thermal behavior with an influence index of 20% and 70% respectively. On the other hand, we also remark that model parameters are impacted by the location and the orientation of the building. Once again, it has been shown that the PWARX model is more suitable for any configuration than the ARX and its derived models. Indeed, through a comparative study between the piecewise ARX model and other existing models such as nonlinear ARX, indexed ARX and ARX models, the PWARX model gives good results in terms of indoor temperature estimation with 78.48% accuracy. Furthermore, a comparison procedure details how the PWARX model could be used to estimate and to compare the thermal behavior of different archetype buildings.

Hence, the future works will be directed to extend this approach for predicting energy performance certificates for existing buildings, for helping to derive strategies for better energy efficiency of the Heating, Ventilation and Air-Conditioning system (HVAC) while taking into consideration the occupant's activities and comforts, and usage modes.

Figure 1 :

 1 Figure 1: Representation of a PWARX system with 3 sub-models [12].

Figure 2 :Figure 3 :

 23 Figure 2: Thermal behavior model identification

Figure 4 :

 4 Figure 4: Example of collected data during one week: Temperature and hygrometry curves

Figure 5 :

 5 Figure 5: Data re-affectation into three affine sub-classes in the regression space [12].

Figure 6 :

 6 Figure 6: (a) Lavoisier student Residence in Douai (b) Lavoisier student Residence 3D model

8

 8

Figure 7 :Figure 8 :Figure 9 :Figure 10 : 2 08Figure 11 :Figure 12 :

 7891021112 Figure 7: (a) Condorcet student Residence in Douai. (b) Condorcet student Residence 3D model

Figure 13 :Figure 14 :

 1314 Figure 13: (a) Indoor air temperature (T i1 and T i2 ). (b) Solar radiation effect on indoor air temperature.

Figure 15 :

 15 Figure 15: (a) Indoor air temperature for the 1 st floor. (b) Heating scenario for the 1 st floor.

FigureFigure 16 :

 16 Figure15ashows that the indoor temperature for the first stage varies between 20 and 31 • C for an outside temperature

FIT

Figure 17 :

 17 Figure 17: Effects of c-nearest-neighbors (c) and the system order (n a , n b ) on the FIT.

Figure 18 :Figure 19 :

 1819 Figure 18: Training data: (a) Model Structure "4". (b) Model structure "4 D " without outside temperature

Figure 20

 20 Figure 20: Influence index

Figure 21 :

 21 Figure 21: Operating mode validation.

Moreover, regarding Figure 22 ,

 22 one can see the comparison of performance of the PWARX-thermal model with other ARX-derived models. The two bottom sub-figures in Figure 22 illustrate the estimation of the indoor temperature with ARX model (green solid line) and indexed-ARX (black solid line) one, compared with the measurements (black solid line).

Figure 22 :

 22 Figure 22: Indoor air temperature validation: Comparison between PWARX, NARX, ARX and Indexed-ARX models.

Figure 23 :

 23 Figure 23: Validation of Temperature (West Side), FIT = 60.02% (14-02 / 17-02).

Figure 24 :

 24 Figure 24: Validation of Temperature (West Side, 3 rd Floor), FIT = 77.66% (14-02 / 17-02).

Figure 25 :

 25 Figure 25: (a) Validation of Temperature (East Side), FIT = 57.40% (14-02/17-02). (b) Validation of Temperature (West Side), FIT = 77.66% (09-02/12-02)

Figure 26 :

 26 Figure 26: Thermal behavior comparison procedure.

Figure 27 :

 27 Figure 27: Similarity rate.

Table 1 :

 1 Wireless sensor characteristics

	Measured parameter Measuring range	Accuracy	Type of sensor	Units
	Temperature	-40 to 125 o C	±0, 4 o C (max), at -10 to 85 o C	Thermistor	[ o C]
	Heat flux	0 -2000W/m 2	40µV/(W/m 2 )	HFP	[W/m 2 ]
	Relative humidity	0 to 100%	±3%(max), at 0 -80%	Capacitive polymer	[%]
	Solar radiation	1 -65536Lux	±100Lux	Light sensor	[Lux] to [W/m 2 ]
	Electrical power	up to 3500W	up to 16A on 220V	Power meter	[W]

  This is achieved with T o , H o , Ra and Pw as inputs. The influence of each input on the accuracy of

	Model	Model structure Inputs selected n a n b C	FIT
	1	2	2	200	T o , H o ,Ra, Pw 71.41%
	2	3	3	200	T o , H o ,Ra, Pw 64.70%
	3	4	4	200	T o , H o ,Ra, Pw 51.60%
	4	5	5	200	T o , H o , Ra, Pw 78.48%
	4 A	5	5	200	T o , H o ,	15.27%
	4 B	5	5	200	Ra, Pw,	55.68%
	4 C	5	5	200	T o , H o , Ra	24.61%
	4 D	5	5	200	Ra, H o , Pw	70.38%

Table 2 :

 2 Best fitting per model structure: (Top) according to system order, (Bottom) according to selected inputs for the 4 th model

Table 4 :

 4 Each behavioral model parameters.

	As an example, let us considerate the following numerical rep-

Table 5 :

 5 Performance comparison between PWARX model and ARX-derived models4.3.3. Impact of the orientations and floor levelsThe parameter's estimation of the PWARX model was made with the measurement campaign from February 07 th -12 th

  mentioned in their study, the same sets of ARX model inputs may not apply directly to a different building type and location due to measurement properties. From this context, hereafter we verified if the reference PWARX model structure estimated for Lavoisier building can reproduced the indoor thermal behavior of two east-west oriented rooms, parts of Condorcet building. The purpose of this step

		28	Temperature Curves		28	Temperature Curves
	[°C]	26		Measurement PWARX Model	[°C]	26	Measurement NARX Model
	Temperature	22 24			Temperature	22 24
	Internal	20			Internal	20
		18				18
	08/02 8:32	09/02 1:34	09/02 18:35	8/02 8:32	09/02 1:34	09/02 18:35
		28	Temperature Curves
	[°C]	26		Measurement ARX Model
	Temperature	22 24
	Internal	20
		18
	08/02 8:32	09/02 1:34	09/02 18:35	08/02 8:32	09/02 1:34	09/02 18:35
			Date-Hour			Date-Hour

Table 6 :

 6 Cross-validation results

Acknowledgements This work was supported by the European project "SHINE: Sustainable Houses in Inclusive Neighborhoods". A project granted by Interreg 2 Seas and the European Regional Development Fund.


MD Gautier Chene
email: chenegautier@yahoo.fr

A Atallah

E Cerruto

D Hallez

Y Castori

E Nohuz

Easy mesh fixation to the promontory during sacrocolpopexy procedure

Keywords: sacral promontory, sacrocolpopexy, ureter, vascular anatomy, mesh, laparoscopy

We propose a simplified technique to approach the promontory and to fix the mesh on the anterior longitudinal ligament during the sacrocolpopexy procedure by using an instrument initially designed for the vaginal approach of the sacrospinous ligament. Using this easy technique in case of anatomic variations on the sacral promontory may avoid complications.

Introduction:

The approach to the promontory during laparoscopic sacrocolpopexy may sometimes be particularly delicate with vascular and / or ureteral wound rates of up to 4.4% and 1% respectively [START_REF] Nygaard | Abdominal sacocolpopexy: a comprehensive review[END_REF][START_REF] Good | Vascular and ureteral anatomy relative to the midsacral promontory[END_REF]. Its dissection is however essential to fix the mesh. The anatomical studies have made it possible to demonstrate the proximity of the vascular and ureteral structures to the site of fixation, particularly in three dangerous structures (because of their proximity): the left common iliac vein, the right ureter and the right internal iliac artery [START_REF] Good | Vascular and ureteral anatomy relative to the midsacral promontory[END_REF].

The average distance between the left common iliac vein and the promontory is 27 mm (extremes vary between 9.5 and 47.5 mm). The average distance between the right ureter and the promontory is 27 mm (the extremes vary between 16 and 48 mm). The average distance between the right internal iliac artery and the promontory is 25 mm (the extremes vary between 14 and 39 mm) [START_REF] Good | Vascular and ureteral anatomy relative to the midsacral promontory[END_REF]. There is also the proximity of the hypogastric nerve and the sacral plexus (between 15 to 30 mm from the promontory) [START_REF] Good | Vascular and ureteral anatomy relative to the midsacral promontory[END_REF][START_REF] Giraudet | The anatomy of the sacral promontory: How to avoid complications of the sacrocolpopexy procedure[END_REF]. Moreover, the abundance of pre-sacral fatty tissue could complicate the precise identification of these different vascular and ureteral elements.

We had the idea of using the Capio ™ Slim device (Boston Scientific ™, Marlborough, USA) which allows sacrospinous fixations to be carried out vaginally (Fig 1A and1B). It is easy to introduce this device through the suprapubic incision or via a 15 mm trocar. The radius of curvature of the Capio ™ Slim head makes it very easy to pass a suture through the anterior longitudinal ligament on the sacral promontory in one step.

Surgical technique: Figures 1 to 3

The patient is set in supine position, in Trendelenburg position and the digestive loops displaced, under general anesthesia, with oro-tracheal intubation, the legs apart in flexion of 30 °, the two arms along the body, the buttocks slightly overhanging the operating table, the catheterized bladder. device is completely removed and the 10 mm trocar replaced. We can then check that the grip is satisfactory by pulling the two strands (Fig 2D). The strands are then inserted into the mesh (Fig 3A and3B). The thread is tied internally or extra-corporeal with the knot pusher (Fig 3C and3D). The rest of the intervention is similar to a conventional sacrocolpopexy [START_REF] Chene | How I do… easily to reduce the operative time of laparoscopic sacrocolpopexy[END_REF].

Conclusion:

Using the Capio ™ Slim device is simple, quick and will effectively fix the mesh on the anterior longitudinal ligament of the sacral promontory. Its use in this indication is being evaluated in our department and may be of service when the dissection of the headland is limited in size because of anatomic variations. It could therefore avoid vascular or ureteral complications. We recall the need for a perfect knowledge of the anatomy of the lumbosacral region and its variability.

Caption:

  Four trocars are required: a 12 mm umbilical trocar for the laparoscope, two 5 mm trocars in the right and left iliac fossa two fingers across the antero-superior iliac spine, and a 10 mm suprapubic trocar or a 15 mm suprapubic trocar (in case of the surgeon wants to put the Capio ™ Slim device through the trocar). Dissection of the promontory requires exposing perfectly the anterior vertebral ligament by first identifying the vascular and ureteral structures to be preserved. The dissection space is limited in width, corresponding to the dimensions of the head of the Capio ™ Slim device. The needle is inserted into the Capio ™ Slim device (Fig 1B). The 10 mm suprapubic trocar is removed. The Capio ™ Slim device is introduced through the same skin incision (Fig 1C and 1D). Another possibility is the use of a 15 mm suprapubic trocar in which the Capio ™ Slim device could be introduced. The assistant in surgery must place his finger with a gauze on the skin incision to avoid losing the pneumoperitoneum. The Capio ™ Slim device is applied to the anterior longitudinal ligament (Fig 2A). The operator presses the plunger to insert the needle into the ligament (Fig 2B). Then surgical forceps removes the crimped wire from the throat of the device head (Fig 2C). The Capio ™ Slim

Figure 1 :Figure 2 :

 12 Figure 1:

Figure 3 :

 3 Figure 3:

Conflict of interest: none


R Shakeri

A Mesgouez
email: arnaud.mesgouez@univ-avignon.fr

G Lefeuve-Mesgouez

Transient response of a concrete tunnel in an elastic rock with imperfect contact

Keywords: Dynamic response, Imperfect interface, Circular tunnel, Semi-analytical approach, Transient wave propagation

The two-dimensional transient response of an imperfect bonded circular lined pipeline lying in an elastic infinite medium is investigated. Imperfect boundary conditions between the surrounding elastic rock and the tunnel are modelled with a two-linear-spring design. The novelty of the manuscript consists in studying at the same time transient regimes and imperfect bonded interfaces for simulating the dynamic response of a tunnel embedded in an elastic infinite rock. Wave propagation fields in tunnel and rock are expressed in terms of infinite Bessel and Hankel series. To solve the transient problem, the Laplace transform and the associated Durbin's algorithm are performed. To exhibit the dynamic responses, influences of various parameters such as the quality of the interface conditions and the thickness of the lining are presented. The dynamic hoop stresses and the solid displacements of both the tunnel and the rock are also proposed.

Introduction

The effect of underground structures on surrounding environments and the effect of surrounding environments on underground structures are of great importance in various fields such as civil engineering constructions, mining or geophysics. More specifically, studying wave propagation in complex heterogeneous media e.g. ground structures is useful not only for geotechnical and mechanical applications or seismic events, but also for medium characterisation since wave propagation signals contain information on the medium itself. Various techniques use for instance the determination of electrical resistivity properties or/and mechanical wave speeds from first arrival times to characterise heterogeneous grounds. Such methods have been applied in the specific underground environment of the LSBB (Low Noise Inter-Disciplinary Underground Science and Technology Laboratory, situated in Rustrel, Vaucluse, southern France), which is constituted of a system of galleries with or without concrete layers [START_REF] Beres | Characterisation of anisotropy in a karstified carbonate plateform using seismic and electrical resistivity methods: a joint approach[END_REF]. It belongs to the Albion plateau, a karstic limestone unit of the southern French Alps. With the aim of better characterising the configuration of the LSBB, a new strategy is proposed and the methodology, based on a mechanical wave propagation modelling and especially the study of the full waveform signal, will follow three main steps in a long term process: (1) implementation of a robust and efficient forward problem, in relation with the LSBB configuration geometry that focuses on the transient responses in terms of stress and displacement fields (short term); (2) a sensitivity analysis [START_REF] Mesgouez | Use of global sensitivity analysis to assess the soil poroelastic parameter influence[END_REF] to extract the main influential mechanical parameters in the LSBB configuration (medium term); and (3) an adequate inversion process related to the latest mechanical parameters and with a full waveform inversion approach to take into account the entire content of the transient response and thus as much information as possible (long term) [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF][START_REF] Lefeuve-Mesgouez | Retrieval of the physical properties of an anelastic solid half space from seismic data[END_REF][START_REF] Mesgouez | Characterization of a viscoelastic heterogeneous object with an effective model by nonlinear full waveform inversion[END_REF]. The paper focuses here on the first part of the project. Resulting from different recent observations, it has been noticed that the concrete layer does not present everywhere a perfect contact with the surrounding rock. Moreover, the digging of the tunnel has weakened some parts of the massif, and a damaged zone has appeared in several areas. It appears then to be important to take into account in the model different kinds of contacts, and to study their effects. A bibliographical review shows that wave propagation modelling can still be improved by taking into account more realistic boundary conditions between the massif and the tunnel for transient excitations. Beyond the specific objective of our study, the work proposed here can also be considered for other applications such as modelling seismic wave propagation in fluid-filled boreholes [START_REF] Sidler | A pseudo-spectral method for the simulation of poro-elastic seismic wave propagation in 2D polar coordinates using domain decomposition[END_REF], deeply buried pipelines subjected to hydraulic fracture initiation or internal blast loads [START_REF] Zhou | Dynamic response of a circular pipeline in a poroelastic medium[END_REF][START_REF] Gao | An analytical solution for the transient response of a cylindrical lined cavity in a poroelastic medium[END_REF], safety of tunnels during earthquakes [START_REF] Yi | Effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves[END_REF].

In the transient regime, researchers have employed numerical methods such as a combination of finite element method (FEM) and boundary element method (BEM) [START_REF] Manolis | Seismic analysis of buried pipeline in a 3D soil continuum[END_REF][START_REF] Clouteau | Freefield vibrations due to dynamic loading on a tunnel embedded in a stratified medium[END_REF] or BEM [START_REF] Tadeu | Scattering of waves by subterranean structures via the boundary element method[END_REF][START_REF] Karakostas | Dynamic response of unlined tunnels in soil with random properties[END_REF][START_REF] Tadeu | 3D scattering of waves by a cylindrical irregular cavity of infinite length in a homogeneous elastic medium[END_REF][START_REF] Panji | Transient analysis of wave propagation problems by half-plane BEM[END_REF][START_REF] Alielahi | Seismic ground amplification by unlined tunnels subjected to vertically propagating SV and P waves using BEM[END_REF] to study the transient wave scattering of cavities or tunnels in half and full space environments and its effect on ground surface and tunnel's responses. Concerning the semi-analytical developments for dynamic studies, Kouretzis et al. developed an analytical methodology to calculate the strain field in buried pipelines due to surface blast load [START_REF] Kouretzis | Analytical calculation of blast-induced strains to buried pipelines[END_REF]. The 3D cylindrical thin shell theory was considered and the authors focused on the analytical results of strain analysis of P-and Rayleigh waves. A good comparison with experimental results and numerical FEM results was achieved. The 3D and transient response of an explosion occurring inside a buried pipeline was also studied by Feldgun et al. who used the Timoshenko elastic plastic shell theory for the pipeline and a Grigoryan model for the soil [START_REF] Feldgun | Internal blast loading in a buried lined tunnel[END_REF]. Smerzini et al. focused on the effects of underground structures, such as circular cavity or elastic shell, on the surface motion of a 2D homogeneous, isotropic and linear viscoelastic half-space subjected to seismic SH excitation in frequency and time domains [START_REF] Smerzini | Effect of underground cavities on surface earthquake ground motion under SH wave propagation[END_REF]. The theoretical approach is based on the expansion of wave functions in terms of Bessel and Hankel functions. Nonlinear dynamic response of functionally graded shells resting on elastic foundations for various geometries [START_REF] Anh | Nonlinear stability analysis of thin FGM annular spherical shells on elastic foundations under external pressure and thermal loads[END_REF][START_REF] Duc | Nonlinear dynamic analysis and vibration of eccentrically stiffened S-FGM elliptical cylindrical shells surrounded on elastic foundations in thermal environments[END_REF][START_REF] Kim | Nonlinear vibration and dynamic buckling of eccentrically oblique stiffened FGM plates resting on elastic foundations in thermal environment[END_REF] has been performed with several classical shell theories, studying particularly the elastic foundation influence on the vibration of functionally graded materials. Kouretzis et al. established a closed-form solution to calculate tunnel forces due to P-wave propagation in an elastic medium under full-slip conditions and showed that Pwaves can lead to higher axial hoop forces compared to S-waves [START_REF] Kouretzis | Analysis of circular tunnels due to seismic P-wave propagation, with emphasis on unreinforced concrete liners[END_REF]. Zhai et al. studied the effect of a cylindrical cavity buried in an elastic soil subjected to an anti-plane Heaviside step load [START_REF] Zhai | Dynamic response of cylindrical cavity to antiplane impact load by using analytical approach[END_REF]. They solved the equations in the Laplace transform domain and the inverse transform was calculated using contour integration techniques along an appropriate integration path. Comparisons both with static results and with Durbin's method were proposed. Fu et al. used indirect boundary element method with Green's function of distributed loads to study soil-tunnel interaction of a rigid tunnel buried in an elastic multilayered soil (an elastic layer over an elastic bedrock) [START_REF] Fu | Dynamic soil-tunnel interaction in layered half-space for incident P-and SV-waves[END_REF].

Beside the transient aspect of the response, one of the most influential parameters for the dynamic response of the mechanical tunnel is the quality of the interface between the rock and the tunnel. Modelling the interface conditions should be close to real considerations. Actually the interface between tunnel and rock is imperfect due to interstitial defaults like foliation, microcracks, mesoscale cavities, faults or joints at the interface and wall roughness. Moreover these discontinuities can be partially saturated with air or water and can thus modify the dissipation of the incoming energy. A few papers modelled imperfect interface with a linear spring [START_REF] Yi | Effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves[END_REF][START_REF] Yi | Dynamic response of a circular lined tunnel with an imperfect interface subjected to cylindrical P-waves[END_REF] or a parallel linear dashpot and spring designs [START_REF] Fang | Visco-elastic imperfect bonding effect on dynamic response of a non-circular lined tunnel subjected to P and SV waves[END_REF]. Indeed, imperfect bonding means here that the stress field is assumed to be continuous whereas the linear displacement field is not at the interface. The displacement jump is assumed to depend linearly on the traction vector, with potential normal and tangential component differences. The dependence can be real and frequency-independent (elastic spring case) or it may be complex and/or frequency-dependent (viscoelastic approach) [START_REF] Schoenberg | Elastic wave behavior across linear slip interfaces[END_REF]. Yi et al. investigated the frequency response of a circular tunnel with imperfect bonded interface in a full space elastic medium modelled with spring constants [START_REF] Yi | Effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves[END_REF][START_REF] Yi | Dynamic response of a circular lined tunnel with an imperfect interface subjected to cylindrical P-waves[END_REF]. Fang and Jin proposed a viscoelastic interface model to study harmonic stress field around a noncircular tunnel in an elastic medium [START_REF] Fang | Visco-elastic imperfect bonding effect on dynamic response of a non-circular lined tunnel subjected to P and SV waves[END_REF]. The viscoelastic imperfect boundary was modelled with a linear spring and a dashpot design.

The above literature review shows that the coupling of transient regimes and imperfect boundary conditions of lined tunnels seems to be absent. The novelty of this paper consists in studying at the same time transient regimes and imperfect bonded interfaces for simulating the dynamic response of a tunnel embedded in an elastic infinite rock. The subject of this paper thus focuses on the 2D transient response of a circular pipeline lying with an imperfect bonded contact in an elastic medium.

This paper is organized as follows. In section 2, a brief statement of the semi-analytical development is recalled. Details on the imperfect boundary conditions are precised. In section 3, the transient hoop stresses and displacements of both the tunnel and the rock are presented. Particularly, 2D snapshots of wave propagation in the rock and the lining are obtained for some parameters for which the effects of the imperfect interface can be observed and quantified. Indeed, the authors have noted that the displacement field in the tunnel is strongly modified and globally lower in comparison with the perfect case response.

Formulation

The two-dimensional problem configuration is illustrated in Fig. 1. A long cylindrical tunnel is embedded in an infinite rock medium. The outer radius of the tunnel is b and its inner radius is a. A transient incident Pwave propagates in the rock in the decreasing x axis direction and acts on the outer surface of the tunnel, as shown in Fig. 1.

Imperfect boundary conditions

Here an imperfect boundary condition hypothesis at the interface between the rock and the tunnel is assumed. Actually the interface zone can present various partially saturated discontinuities like foliation, microcracks, mesoscale cavities, faults or joints. These discontinuities can significantly impact the dynamic response of the tunnel. Even if this inhomogeneous area is quite difficult to apprehend, because it is mainly phenomenological and due to the intrinsic nature of the rock and the way the tunnel itself has been built, one of the most well-known models taking into account the imperfect boundary conditions is a two-linear-spring design: displacement field is assumed to be discontinuous along the entire length of the interface whereas traction field is assumed to be continuous. Applications of this concept do not exclusively concern seismic wave propagation but also non-destructive tests and natural or industrial composite materials. Relationships between experimental observations and the two-linear spring model can be found in [START_REF] Honarvar | Nondestructive Evaluation of Cylindrical Components by Resonance Acoustic Spectroscopy[END_REF][START_REF] Lourme | Etude des assemblages collés par ondes guides ultrasonores: étude expérimentale et modélisation par éléments finis[END_REF][START_REF] Vlasie | Mechanical and acoustical study of a structural bond: comparison theory/numerical simulations/experiment[END_REF][START_REF] Valier-Brasier | Scaled behavior of interface waves at an imperfect solid-solid interface[END_REF], among others.

Based on this concept, the imperfect boundary conditions at the interface surface between the lining and the rock at r=b can be written as follows [START_REF] Yi | Effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves[END_REF][START_REF] Schoenberg | Elastic wave behavior across linear slip interfaces[END_REF][START_REF] Achenbach | Wave propagation in elastic solids[END_REF].

1 2 1 2 1 1 2 1 1 2 , , ( ), ( )
rr rr r r rr r r r r k u u k u u θ θ θ θ θ θ σ σ σ σ σ σ = = = - = - (1
) And the boundary conditions at the inner surface of the lining at r=a are written as a free surface. (2) where kr and kθ are the normal and tangential spring constants, respectively; u the displacement vector; and σ the stress tensor. Actually kr and kθ represent the roughness in radial and tangential directions respectively between rock and tunnel. Dividing the third equation in Eq. ( 1) by kr or the fourth equation in Eq. ( 1) by kθ, when these parameters tend to infinity, the left-hand sides tend to be zero, implying no jump in displacement between the tunnel and rock. This limit case represents thus the perfect continuity conditions; consequently, the perfect boundary conditions mean hereafter the largest values of kr and kθ. On the contrary, very low values of kr and kθ mean that traction efforts at the interface do not exist anymore.

Waves in rock and tunnel

The usual gradient and curve decomposition of the displacement vector,

φ Ψ = ∇ + ∇ × u
, is used here leading to the following Helmholtz equations. (1 )( ), (1 )( ) (where µ and λ are the Lamé constants and ρ is the material density). The displacement field u does not depend on the z-coordinate due to the 2D studied configuration, so the vector potential Ψ only has one single component ψz. For simplicity it will be written as ψ. The Laplace transform defined by 0 ( , , ) ( , ,

p s c t c t φ φ ∇ = ∂ ∂ ∇ = ∂ ∂ Ψ Ψ (3) where ( 2 )
) st f r s f r t e dt θ θ ∞ - = ∫ (4)
is then applied to Eq. ( 3). The transient plane P-wave propagates in the rock in the negative x-direction, as shown in Fig. 1. The time variable is adjusted so that the incident wave front strikes the tunnel at r=b and θ=0 at t=0. The total scalar potential t 1 φ , in the rock medium (subscript 1 denotes the wave features in the rock medium) and in the Laplace domain, is then composed of incident and refracted components as written in the following.

1 0 i t 1 1 1 0 0 refracted P-wave , incident wave ( ) ( ) cos ( ) ( ) cos x n n n n n n e F s I r n A s K r n α φ φ ε α θ α θ ∞ ∞ - = = = + ∑ ∑ 14444 4244444 3 1444 424444 3 (5)
where 1 α =s/cp1 is the wavenumber of P-wave in the Laplace domain in the rock; x0 the distance between the tunnel and incident wave before t=0; F(s) an arbitrary excitation term expressed in the Laplace domain; εn=1 for n=0 and εn=2 for n≥1; Kn the modified Bessel function of the second kind of order n; and In the modified Bessel function of the first kind of order n.

Due to the interface between the rock and the tunnel, two refracted waves (P-wave as indicated in Eq. ( 5) and SV-wave in Eq. ( 6)) in the rock medium propagate when the incident wave acts on the outer surface of the tunnel. Since the scattering field vanishes when r→∞, they are represented by modified Bessel function of the second kind Kn [START_REF] Zhai | Dynamic response of cylindrical cavity to antiplane impact load by using analytical approach[END_REF]. Consequently, the total vector potential t 1 ψ , in the rock medium and in the Laplace domain, can be written as

t 1 1 0 refracted SV-wave ( ) ( )sin n n n B s K r n ψ β θ ∞ = = ∑ 1444 424444 3 (6)
where An(s) and Bn(s) are unknown coefficients; and 1 β =s/cs1. Also, for the tunnel medium (subscript 2) the scalar and vector potentials can be written as a combination of refracted and reflected waves as follows.

t 2 2 2 0 [ ( ) ( ) ( ) ( )]cos n n n n n C s I r D s K r n φ α α θ ∞ = = + ∑ (7) t 2 2 2 0 [ ( ) ( ) ( ) ( )]sin n n n n n M s I r N s K r n ψ β β θ ∞ = = + ∑ (8)
where Cn(s), Dn(s), Mn(s) and Nn(s) are unknown coefficients. From Eqs. ( 5)-( 8), the potential decomposition of displacements, and the stress-strain relations, the displacement and stress components can be obtained everywhere in terms of the unknown coefficients. Then combining with the interface and boundary conditions, and with the use of the orthogonality properties, a matrix system can be written as

1 0 ( ) ( ) ( ) ( ) x n n n n s s e F s s α - = ∈ E Q G ( 9
) where Qn(s) includes the unknown coefficients as follows.

( ) [ ( ), ( ), ( ), ( ), ( ), ( )]

t n n n n n n n s A s B s C s D s M s N s = Q (10
) And equations for Gn(s) and En(s) are shown as follows (see more details in [START_REF] Shakeri | Dynamic response of a circular tunnel with imperfect surface interaction embedded in an elastic medium[END_REF]).

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (2 ( )) ( ) ( ) (2 ( )) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (2
µ α µ α α α α α µ α µ α α α α α µ α α µ α α α α α α µ - - - - - - = E 51 1
( )) ( ) ( ) (2 ( )) ( ) ( ) ( ) ( ) ( ) ( ) 0 0 ( ) ( ) ( ) ( ) 0 0 ( ) ( ) ( ) ( )
α µ α α α α α α α α α α α α α α           - -           (12)
where µ =μ1/μ2; and eij are given as follows.

= - + = - -- 2 2 2 1 2 2 2 2 2 41 1 42 1 2 2 2 2 2 43 1 44 1 5 2 ) ( ) ( ) ( ) ( 2 ) ( ) ( ) ( ) [( 1) ( ) ( )] ( ) ( 2 ) ( ) ( ) ( ) [( 1) ( ) ( )] ( ) [( 1)
= -+ + = - -- + + = -+ - = - - 1 51 1 2 2 2 52 1 53 1 2 2 2 54 1 ) ( )] ( ) [( 1) ( ) ( )] ( ) ( 2) ( ) ( ) ( ) [( 1) ( ) ( )] ( ) ( 2)
Solving Eq. ( 9) for different values of

[0, ] n N ∈
, with N being the truncation constant of the infinite series, gives coefficients An(s) … Nn(s) and thus the displacement and stress fields everywhere in the Laplace domain.

Durbin's method

The inversion of Laplace transforms to obtain the transient responses is performed numerically by using the following Durbin's formula in the interval [0,2T0] [START_REF] Durbin | Numerical Inversion of Laplace Transforms: An Efficient Improvement to Dubner and Abate's Method[END_REF][START_REF] Hasheminejad | Vibro-acoustic response of an elliptical plate-cavity coupled system to external shock loads[END_REF].

0 0 0 0 0 1 ( ) (2 )[1 [2 Re[ ( )]]] [Re[ ( 2 )]cos(2 ) Im[ ( 2 )]sin(2 )]] M ct k t e T c c ik T kt T c ik T kt T π π π π = Λ = Λ + Λ + - Λ + ∑ (14)
where M is the truncation parameter; and c an arbitrary real number. Λ(s) can be any of ( )

r u s , ( ) u s θ , ( ) rr s σ , ( )
s θθ σ
. The suggested value of cT0 and M are between 5-10, and 50-5,000 for sufficient accuracy, respectively [START_REF] Durbin | Numerical Inversion of Laplace Transforms: An Efficient Improvement to Dubner and Abate's Method[END_REF].

Numerical results

Validation

To validate our development, and due to the lack of transient response when taking into account an imperfect boundary scenario, a comparison is proposed with the frequency approach of Yi et al. [START_REF] Yi | Effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves[END_REF]. The comparison concerns the harmonic response of a tunnel embedded in an elastic rock subjected to incident plane compressional wave. Second, in order to show the reliability of transient response, a cylindrical cavity in the elastic medium under a step load acting on the inner surface of the cavity is considered and compared to [START_REF] Zakout | Transient response of a cylindrical cavity with and without a bonded shell in an infinite elastic medium[END_REF]. The following nondimensional parameters have been chosen: v1=v2=0.3, μ1/μ2=1, cp1/cs1=cp2/cs2=√3.5, b/a=1.2, cs1/cs2=1 and kr=kθ=100μ1/b. A big value for the spring constant at the interface between the rock and the lining changes the imperfect boundary condition to a perfect interface, as discussed previously. Good agreement is obtained, as shown in Fig. 3, for the mode of the displacement (n=0) and at points (r=a, θ=0).

Main results

Mechanical values (tunnel and rock parts) are given in Table 1 from Fang and Jin [START_REF] Fang | Visco-elastic imperfect bonding effect on dynamic response of a non-circular lined tunnel subjected to P and SV waves[END_REF]. The outer and inner radius of the tunnel are b=1 m and a=0.8 m. Results are obtained for the following dimensionless quantities.

* * 2 * * * * * 1 1 1 1 2 2 1 1 1 1 2 2 1 , , , 1, , , p p p
= = = = = = = (16)
The incident plane wave in the Laplace domain in Eq. ( 5) is chosen as F(s)=1/(s+3), and the corresponding time expression is f(t * )=e -3t * . Durbin's formula is evaluated with M=1,200, c=5/T0 and T0=60. Fig. 4 displays the dimensionless hoop stress of the rock * 1 θθ σ at (r=b, θ=0 and π/2) and the tunnel * 2 θθ σ at (r=a, θ=0 and π/2) for the same three different values of imperfect boundary condition coefficients. Recall that perfect conditions are modelled with the highest values of kr and kθ as explained previously, corresponding here to the black curves in Fig. 4. The dimensionless hoop stress in the rock at (r=b, θ=0) is increased when tending towards the perfect bond configuration. On the contrary, for (r=b, θ=π/2) the dimensionless hoop stress is slightly larger for the weak bond (blue curves). Actually this fact has been observed in the harmonic response studied by Yi et al. who showed the complicated effects of kr and kθ with respect to θ [START_REF] Yi | Effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves[END_REF]. Consequently it cannot be asserted that the displacement or stress fields are larger in the perfect case for every θ, even if the overall response for perfect boundary conditions clearly presents higher amplitudes. By increasing the value of kr and kθ the dimensionless hoop stresses of the tunnel in both points are strongly increased (lower figures) due to the higher level of wave transmission. The most important point here concerns the oscillating feature of the dynamic response of the tunnel, and to a lower extent of the rock, when the perfectly bonded interface conditions are considered. On the contrary, for very low values of kr and kθ, or highly imperfect boundary conditions, the materials are detached from each other and tunnel does not affect rock response significantly. Actually, the response of the rock is similar to the case without tunnel (dashed lines).

Values of kr=kθ=1µ1/b lead to intermediate responses for both the tunnel and the rock (see red curves in Fig. 4). Since the dimensionless P-wave velocities in rock and tunnel are * 1 p c =1 and * 2 p c =0.91 respectively, the elastic wave arrives sooner in the rock in comparison with the tunnel. Also, a short delay in the response of the tunnel at (r=a, θ=0) can be seen due to the thickness of the tunnel. Observations mentioned for Fig. 4 are confirmed for the displacement analysis concerning the influence of the boundary conditions on the amplitude of the response and on the oscillating behaviour. For instance, for θ=0 the nondimensional displacement in tunnel (left lower sub-figure) shows differences up to several hundred percent between cases kr=kθ=10µ1/b and kr=kθ=1µ1/b (black and red curves respectively). This illustrates the fact that the incoming energy for weak contact is consumed between tunnel and rock due to rough interface, and that the tunnel receives less energy. A pseudo-periodicity around T * ≈2×(0.2/0.91)≈0.44 due to the multiple reflections in the tunnel and specifically for the perfect bonded case is also observed.

Figs. 6 and 7 are analogous to Figs. 4 and 5 respectively but for a 0.4 m tunnel lining thickness (b=1 m, a=0.6 m). The main difference concerns the oscillations: the pseudo-period is doubled due to the thickness of the lining. The amplitudes are not so much perturbed and the boundary condition type impacts the mechanical fields in a similar trend.

Figs. 8-11 exhibit the nondimensional displacement snapshots of the rock and the tunnel for various selected times, for perfect and imperfect boundary conditions, for a 0.2 m lining thickness (Figs. 8 and9) and a 0.4 m lining thickness (Figs. 10 and11). Fig. 12 presents the snapshots relative to the case of the rock without lining. The plane wave front is clearly seen and is deformed when it impacts the tunnel surface. In the perfect boundary condition cases, the nondimensional displacement of the tunnel is larger in comparison with the imperfect cases. As mentioned in section 2.2, concerning the imperfect boundary conditions, discontinuities and faults between the tunnel and the rock behave like a barrier against the input energy of incident wave. The multiple reflections in the tunnel, discussed in Figs. 4 and5, are thus observed in the perfect contact configurations, leading to disturbed fields. In the imperfect cases, the displacement field in the tunnel is smoother. No notable difference in rock displacement is noticed.

Conclusion

In this paper, the 2D transient response of an imperfect bonded circular lined pipeline lying in an elastic infinite medium has been obtained. Imperfect interaction surface between elastic medium and tunnel has been considered as a two-linear-spring model. The elastodynamic wave equations of P-and S-waves have been solved in the Laplace domain, using Helmholtz decompositions and Bessel and Hankel series expressions. Durbin's numerical Laplace transform inversion approach has been used to obtain transient responses.

The global dynamic semi-analytical development has been validated by comparison with previous results. Indeed, the closed form solution in this study paves the way for future to optimize the mid-layer properties in order to get the acceptable dynamic stress around the underground structures for any type of wave loading. The influences of imperfect conditions on both the dynamic hoop stresses and displacements of the tunnel and rock have been analyzed. The tunnel behavior is greatly impacted by the nature of the interface. Perfect bond interface configurations have shown oscillating features in the responses of the tunnel. Weak bond interface configurations have led to low overall level of dynamic responses. No significant effect on the displacement of the rock has been observed.

Further works will consider from the theoretical point of view both the implementation and the sensitivity study of a damaged zone area and its specific contact with the tunnel and the rock, and the development of anisotropic and/or porous characteristics of both the lining and the external environment. From the experimental point of view, specific results extracted directly from the Low Noise Inter-Disciplinary Underground Science and Technology Laboratory will be compared, after data treatment, to simulation results.

Fig. 1 .

 1 Fig. 1. Problem geometry.

=

  are the P-and S-wave velocities in the elastic medium, respectively

Fig. 2

 2 presents the dynamic stress concentration factor * θθ σ of the rock at r=b, recalled by the amplitude of the incident wave. The hoop stress is considered here since, in cylindrical tunnel configurations, it is bigger than radial or axial stresses. The rock and lining parameters in Fig. 2 are taken directly from Yi et al.: v1=0.25, v2=0.3, μ1/μ2=0.31, cp1/cp2=0.7, b/a=1.2 and α1a=2 [9]. The same three different values of kr and kθ are also considered. Excellent agreement between present results and those from [9] is obtained.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Comparison of dimensionless hoop stress of the present results and results of [9] in the rock for different values of kr and kθ.

Fig. 5

 5 exhibits the dimensionless radial displacements of the rock * 1 r u at (r=b, θ=0 and π/2) and the tunnel * 2 r u at (r=a, θ=0 and π/2) for the same three different values of imperfect boundary condition coefficients.

Fig. 4 .

 4 Fig. 4. Dimensionless hoop stresses of rock and tunnel for θ=0 and π/2, different ratios of kr and kθ and a 0.2 m lining thickness.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Dimensionless radial displacements of rock and tunnel for θ=0 and π/2, different ratios of kr and kθ and a 0.2 m lining thickness.

Fig. 7 .

 7 Fig. 7. Dimensionless radial displacements of rock and tunnel for θ=0 and π/2, different ratios of kr and kθ and a 0.4 m lining thickness.

Fig. 8 .

 8 Fig. 8. Dimensionless displacement snapshots of rock and tunnel for large value of kr and kθ (perfect interface case) for selected times and for a 0.2 m lining thickness.

Fig. 9 .

 9 Fig. 9. Dimensionless displacement snapshots of rock and tunnel for small value of kr and kθ (imperfect interface case) for selected times and for a 0.2 m lining thickness.

Fig. 10 .

 10 Fig. 10. Dimensionless displacement snapshots of rock and tunnel for large value of kr and kθ (perfect interface case) for selected times and for a 0.4 m lining thickness.

Fig. 11 .

 11 Fig. 11. Dimensionless displacement snapshots of rock and tunnel for small value of kr and kθ (imperfect interface case) for selected times and for a 0.4 m lining thickness.

Fig. 12 .

 12 Fig. 12. Dimensionless displacement snapshots of rock without lining for selected times.

Table 1

 1 Material properties for the rock and tunnel.

	Parameters	Rock	Tunnel
	ρ (kg/m 3 )	2,730	2,500
	E (GPa)	35	30
	ν		0.3	0.25
	µ (GPa)	13.462	12
	λ (GPa)	20.192	12
	cp (m/s)	4,154.3	3,794.7
	cs (m/s)	2,220.6	2,190.9
	* p c	1	0.9134
	c	* s	0.5345	0.5274

Acknowledgements

This work has been carried out under a PhD grant from Conseil Général Provence-Alpes-Côte d'Azur, France. This support is gratefully acknowledged. The authors also wish to express their gratitude to Stéphane Gaffet, Head of the Low Noise Inter-Disciplinary Underground Science and Technology Laboratory for his constructive comments during this project.