Kashif Saeed
email: kashif.saeed-5@etudiants.ensam.fr

Nazih Mechbal

Gérard Coffignal

Michel Vergé

Recursive modal parameter estimation using output-only subspace identification for structural health monitoring

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

stabilization diagrams [START_REF] Van Der Auweraer | Modal parameter estimation from inconsistent data sets[END_REF], [START_REF] Allemang | Correlation coefficient for modal vector analysis[END_REF], which exhibit the variation of modal parameters relative to increments in modal order, making them an effective tool for spurious mode identification. However, they require extensive user interaction as some numerical spurious modes behave as physical ones, when the model order is changed making them very difficult to identify on stabilization diagrams.

In order to reduce computation complexity, recursive SubID algorithms have been proposed. Mercère et al [START_REF] Mercère | Recursive subspace identification based on instrumental variable unconstrained quadratic optimization[END_REF] and Delgado et al [START_REF] Delgado | Recursive canonical variate subspace algorithm[END_REF] have developed methods that identify system matrices without extracting state sequence information. Whereas, Pongpairoj and Pourboghrat [START_REF] Pongpairoj | Real-time optimal control of flexible structures using subspace techniques[END_REF] estimated state sequences from approximated subspace intersection of sequential updated vector. Methods [START_REF] Fan | Data processing in subspace identification and modal parameter identification of an arch bridge[END_REF], [START_REF] Cauberghe | On the influence of the parameter constraint on the stability of the poles and the discrimination capabilities of the stabilisation diagrams[END_REF] have been proposed to reinforce stabilization diagrams but they do not eliminate spurious modes automatically. Therefore, the problems of improving the accuracy of modal parameter estimates and the automatic elimination of spurious modes using covariance driven output-only SubID methods are addressed in this paper.

We propose a two step iterative procedure for improving the accuracy of modal parameter estimation on active structures, while using covariance based output-only SubID. In the first step, excitation force is iteratively modified to have a better estimate of the modal parameters. In the second step, an attempt is made to eliminate the spurious modes by using an original stabilization histogram and by resampling of measured output data. In order to reduce the complexity of modal testing and to increase the flexibility of implementation, output-only SubID method was applied.

The paper is organized as follows. A brief introduction to output-only covariance driven SubID algorithm is given in section 2. The proposed algorithm is introduced in section 3. Results from two case studies are demonstrated in section 4 and 5. Section 6 concludes the paper with future prospects.

II. MODAL ANALYSIS BY SUBSPACE IDENTIFICATION

A. State Space Model

It is generally convenient to convert the second order differential equations of a flexible structure in state-space form. If 1/τ is the sampling rate then the discrete state-space model is given by

x k+1 = A d x k + B d u k + w k (1)
y k = C x k + v k (2)

I. INTRODUCTION

Detecting and localizing damage in smart structures is a topic of growing interest in Structural Health Monitoring (SHM) [START_REF] Lestari | Curvature mode shape-based damage assessment of carbon/epoxy composite beams[END_REF], [START_REF] Hurlebaus | Smart structure dynamics[END_REF], [START_REF] Mechbal | Application of a combined active control and fault detection scheme to an active composite flexible structure[END_REF]. These structures have built in piezoelectric sensors and actuators, which can be used to identify the changes caused by damage in structural parameters, such as modal properties. The sensitivity and use of modal properties for damage detection have been extensively studied [START_REF]An overview of intelligent fault detection in system and structures[END_REF], [START_REF] Sohn | A review of structural health monitoring literature: 1996-2001[END_REF] and many modal parameter identification methods have been proposed. One of the modal parameter identification methods, which has gained popularity in recent years, is Subspace Based Identification (SubID) [START_REF] Katayama | Subspace methods for system identification[END_REF].

SubID methods estimate the state-space sequence using input/output signals. The original SubID method is noniterative and is restricted only to offline implementation because of computational complexity of robust matrix operations such as singular value decomposition (SVD). Moreover, like many other identification methods, SubID depends heavily on excitation force. Weak ambient excitations along with ambient noise have strong influence on the accuracy of estimates. In order to benefit from ambient excitation, long samples are used with models of redundant order. However, using model of redundant order generates spurious modes, which mix intensively with physical modes, making them difficult to identify. It has been observed that as the model order is increased, the physical modes tend to stabilize as compared to the spurious ones. This is the core idea behind where, x k ∈ R ns , u k ∈ R nu and y k ∈ R ny are respectively the state-space, input and output vectors at time instant k. n d is the system order, n s is state-space order, n y is the number of outputs and n u is the number of inputs. A d , B d and C are the discrete system, input and the output matrices respectively. w k and v k are respectively the process and the measurement noise vectors.

Considering that the system behaves like a stationary linear dynamics system and that input forces can be modeled as a non-stationary white noise, the model can be modified as

x k+1 = A d x k + w k (3)
y k = C x k + v k (4)
Here, the process and the measurement noise are assumed to be zero-mean Gaussian white noise processes such as

E w k v k w k v k = Q S S T R δ(t) (5)
where, E is the expectation operator and δ(t) is the Kronecker delta. Q and R are process and measurement noise covariance matrices. S is the cross-covariance between process and measurement noise.

B. Output-only Covariance Driven Subspace Identification

Covariance driven SubID is based on the following steps [START_REF] Mevel | Input/output versus output-only data processing for structural identification-application to in-flight data analysis[END_REF]. Let Λ i be output covariance matrix given by

Λ i = E y k+i y T k (6)
A Hankel matrix is constructed using p + 1 block rows and q block columns of the output covariance matrix.

H p+1,q =       Λ 1 Λ 2 • • • Λ q Λ 2 Λ 3 • • • Λ q+1 . . . . . . . . . . . . Λ p+1 Λ p+2 • • • Λ p+q+1       (7)
Considering that the system is stationary and that the process and measurement noises are white, Gaussian with zero mean, we can rewrite [START_REF] Katayama | Subspace methods for system identification[END_REF] as

Λ i = C Σ C T + R i = 0 (8)
Λ i = C A i-1 d G i > 0 (9)
where, G is the cross-covariance between state and observed output. This leads to the factorization of Hankel matrix as

H p+1,q = O p+1 C q (10)
where

O p+1 = C C A d • • • C A p d T (11)
and

C q = G A d G • • • A q-1 d G (12)
are p + 1 order Observability and q order Controllability matrices of rank n d respectively.

In order to calculate the system state matrix A d the first block row of the observability matrix O p+1 is deleted to get a truncated observability matrix O ↑ p+1 given by

O ↑ p+1 = C A d C A 2 d • • • C A p d T (13)
Thus, the system state matrix A d can be calculated by

A d = O ⊕ p O ↑ p+1 (14)
where, O ⊕ p represents the pseudo-inverse of p order observability matrix. The output state matrix C can be obtained from the first block row of O p . Thus, the state-space model given in ( 3) and ( 4) can be identified.

Modal parameters are estimated by using the output matrix C and system matrix A d . Eigenvalues of A d yield the natural frequencies f i and damping ratios ζ i of the system. Modeshapes φ i are obtained in the dimension of the measured degrees of freedom from corresponding eigen-vectors.

C. Balanced Realization

For the time being, we use Balanced Realization [START_REF] Akaike | Stochastic theory of minimal realization[END_REF] for the implementation of SubID, which uses SVD to obtain observability and controllability matrices. In practice, only finite number of samples are available. As a result, the output covariance matrices are calculated by

Λi = 1 n τ -i nτ -i k=1 y k+i y T k (15)
where, n τ is the number of samples. Experimental Hankel matrix Ĥp+1,q is populated using the estimated output covariance matrices Λi . The SVD of Ĥp+1,q and its truncation at a desired model order, yields the estimated observability Ôp+1 and controllability Ĉp+1 matrices

Ĥp+1,q = U 1 U 2 S 1 0 0 S 2 V 1 V 2 T (16) Ôp+1 = U 1 S 1/2 1 (17) Ĉp+1 = S 1/2 1 V 1 (18)

III. RECURSIVE METHOD

For a damage detection algorithm, using changes in modal parameters as an indication of damage, underling modal parameter identification method must be precise and accurate. To achieve this goal, we propose an alternative recursive approach using a stabilization histogram.

A. Stabilization Histogram

Order stabilization diagram exhibits variation of modal parameters estimates with increments in model order. Though increasing model order reveals stability of physical modes but it also increases the interference of spurious modes. There are two types of spurious modes, the ones generated by noise/weak ambient excitation and others resulting from numerical inaccuracies. By far, stabilization diagram is the most commonly used method for physical mode identification [START_REF] Van Der Auweraer | Modal parameter estimation from inconsistent data sets[END_REF] but it requires extensive user interaction as some numerical spurious modes behave like physical ones. In order to automate modal parameter extraction from SubID method, an alternative approach using a histogram is proposed.

To construct a stabilization histogram, we propose first, to categorize the modes in n c categories. The modes are categorized according to the stability of natural frequencies, damping ratios and modes shapes, with respect to increments in model order. For example, here, we define three categories: stable modes, semi-stable modes and new modes. A mode calculated for a given order is said to be stable if frequency varies less than 1% from the previous order and damping ratio varies less than 5%. Semi-stable modes have stable frequencies but no stable damping coefficients. New modes have neither stable frequency nor stable damping coefficients. Adding a mode shape stability criteria will increase the number of catagories and might improve results.

Weights are associated with each category. The weight of a more stable category is selected to be comparatively higher than all other less stable categories. Considering the previous stated example, the weight of a stable mode is greater than that of a semi-stable mode and the weight of semi-stable mode is greater than the weight of a new mode.

The working frequency band is divided into discrete intervals. With each interval a counter is associated. For an estimated mode that falls in a given frequency interval, the counter is incremented by the weight of the associated category. As a result stable modes will have higher counts as compared to less stable modes. Plotting counts with respect to frequency, results in a histogram that shows higher peaks for more stable modes.

B. Excitation Signal

Results generated by output-only covariance driven SubID method depend heavily on the excitation force, which is considered here to be unknown. In order to improve the results and yet still retain the simplicity of this algorithm, it is in our interest to excite a structure with a well defined input signal but disregard it during the identification process. We propose an iterative procedure, actively changes this input based on the results of previous identification attempts. The procedure is given as follows:

1) Data acquisition from ambient excitation 2) SubID and construction of stabilization histogram 3) Generation of a pseudorandom binary signal having considerable energy around the identified frequencies 4) Data acquisition from forced excitation 5) SubID and construction of stabilization histogram 6) Going to step 3 till the identified frequencies stabilize We assume that on an active structure the excitation points are well optimized to excite the modes in the working frequency band. We also assume that during the period of algorithm convergence, the structural properties will not change. After the stabilization of input signal, spurious modes should be eliminated so that results can be used for damage detection algorithm.

C. Spurious Mode Rejection

We have observed that as spurious modes generated by noise do not stabilize while increasing the model order, however, the numerical spurious modes do tend to stabilize and can be considered as physical modes of the system. We have also observed that these numerical spurious modes depend on the sampling period. In the second step of spurious mode rejection, we propose to downsample the system output, which was originally measured at a high sampling rate and re-estimate the modal parameters. Injecting these results, from different sampling rates, in a single stabilization histogram reveals spurious modes because physical modes remain stable under both conditions of changing sampling rate and model order. Flow diagram of the proposed algorithms is shown in Fig. 1 IV

. CASE STUDIES -ACTIVE COMPOSITE BEAM

The proposed online modal parameter estimation method is applied to an active composite beam structure (Fig. 2). The beam consists of a composite filling covered with two external thin plates. There are three pairs of PZ29 piezoelectric ceramics, bounded symmetrically on both sides of this beam. They are positioned parallel to the mid plate surface and are polarized in a way that permits sensing or generating pure bending motion.

The aim is to accurately and automatically identify the first five natural modes of the structure, using measurements from a single piezoelectric pair and one excitation source. Using output measurements from a single piezoelectric pair restricts us to only consider the stability of natural frequencies and damping coefficients. The categories defined in the example stated in section III will be used for all the results shown in this paper. Weights selected for stable, semi-stable and new modes are 8, 2 and 1 respectively and the frequency interval set to 0.1Hz.

A. Simulation

A finite element model (FEM) of the beam is constructed using Bernoulli hypothesis, for simulations purposes. The modeled structure is damped satisfying Caughey's criteria [START_REF] Adhikari | Damping modelling using generalized proportional damping[END_REF]. Natural frequencies and the corresponding damping ratios are shown in Table I. The dynamic response of piezoelectric pair 2, to a random signal excitation by piezoelectric pair 1, is simulated using Newmark Integration method. White noise is added on the simulated response with signal to noise ratio of 10dB.

Stabilization diagram generated after applying SubID to the simulated response data, measure at 2500Hz, is shown in Fig. 3. The corresponding stabilization histogram is shown in Fig. 4. Both the diagrams show difficulties in identifying the first two modes of the beam. The histogram shows high peaks for the two well identified modes at 273Hz and 490Hz. However, counts for the first two modes are lower as compared to spurious modes. Using the first 15 frequencies, having the highest count, an input signal is generated by the proposed algorithm, in order to better excite the structure in the next iterations.

Five peaks on the stabilization histogram (Fig. 5), obtained after 3 rd iteration, correspond well with the first five modes of the beam. As there are no spurious modes, the final estimated values are shown in Table II.

To demonstrate the effectiveness of spurious mode rejec- tion, the proposed algorithm is implemented on an experimental setup.

B. Experimental Setup

The experimental setup is completed by a control loop: charge amplifiers for the conditioning of measurement signals, voltage amplifiers for the actuators and a specialized dSpace c card performing the real time measurements and control. Initially the beam is excited by random signal using the 1 st piezoelectric pair as actuator. Vibration data is measured from 2 nd piezoelectric sensor at a sampling rate of 10kHz for a duration of 2 minutes. This data is downsampled to 1250Hz and then used to generate the stabilization diagram and the corresponding histogram, given in Fig. 6 and Fig. 7 respectively. These figures reveal that it is difficult to identify the first and the second modes of the structure. Moreover, there are spurious modes, which stabilize with increasing model order.

The advantage of stabilization histogram is that it not only allows better identification of first mode but also assists in constructing a richer signal for excitation. Exciting the beam with a refined input, generated iteratively, yields the final stabilization diagram shown in Fig. 8 The 4 th mode at 257Hz is weakly identified and there rests a numerical spurious mode, which stabilizes with increasing model order at 430Hz. To eliminate the spurious mode, the response is measured at a higher sampling rate of 10kHz. It shifts the numerical spurious modes to higher frequencies but the task of identifying the low frequency modes becomes more difficult (Fig. 9).

Histograms allow the user to analyze extremely large data Counts of all the first five vibration modes are well above the selected threshold, allowing them to be automatically selected. Whereas, the effect of the numerical mode is cancelled by changing the sampling period. Final estimated values of the first five modes are shown in Table III.

V. CASE STUDIES -CACTUS

The proposed algorithm is also tested on a thin aluminum plate like structure called CACTUS (Fig. 11). CACTUS is 2mm thick and is clamped in vertical position. There are nine pairs of PZ29 piezoelectric ceramics, bounded symmetrically on both sides of the plate.

The aim is to identify the first seven natural modes of the structure up to 55Hz. Initially, the structure is excited by random signal using the 1 st piezoelectric pair as actuator.

Abstract-

  Precise identification of structural properties is a vital step towards detection and localization of damage in structures. In this paper the problems of improving accuracy of modal parameter estimates and automatic elimination of spurious modes in covariance driven output-only subspace identification method are addressed. An iterative procedure is proposed, which in the first step, actively modifies the excitation signal, resulting in improvements in identification results. In the second step, for spurious mode elimination, an alternative stabilization histogram is introduced to automatically combine and extract identified modal parameters. It is shown that use of measured output signals, of different sampling rates, along with combination of identified results on a single stabilization histogram, can enhance the effectiveness of spurious modes rejection. One numerical and two experimental examples, on the modal parameter estimation of a composite beam and an aluminum plate (CACTUS) are presented to demonstrate the efficacy of the iterative procedure. The proposed algorithm allows automatic and accurate generation of modal parameter residuals for structural health monitoring applications.

Fig. 1 .

 1 Fig. 1. Online Modal Parameter Estimation using Output-only Subspace Identification.

Fig. 2 .

 2 Fig. 2. The Active Composite Beam Structure.

Fig. 3 .

 3 Fig. 3. Stabilization diagram obtained after random excitation of modeled beam after 1 st Iteration. Stable, semi-stable and new modes are represented by square, circle and dot respectively. Right y-axis corresponds to Bode diagram obtained from the FEM model.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Stabilization histogram obtained after 1 st Iteration from random excitation of modeled beam.

Fig. 7 .

 7 Fig. 7. Stabilization histogram obtained from experimental data measured at a sampling rate of 1250Hz with random excitation.

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig. 8. Stabilization diagram obtained from experimental data measured at a sampling rate of 1250Hz with improved excitation

TABLE I FIRST

 I FIVE MODES OF THE MODELED BEAM

	Modes	1 st	2 nd	3 rd	4 th	5 th
	Natural Frequency (Hz)	6.7	44.8 126.2 273.5 490.5
	Damping Ratio ×10 -3	4	6	7	15	11

TABLE III FIRST

 III FIVE IDENTIFIED MODES OF THE COMPOSITE BEAM

	Modes	1 st	2 nd	3 rd	4 th	5 th
	Natural Frequency (Hz)	6.2	44.1 126.5 255.5 499.0
	Damping Ratio ×10 -3	4.2	5.5	6.9	98.4	26.7

Vibration data is measured from 2 nd piezoelectric sensor at a sampling rate of 1kHz for duration of 1 minute. Stabilization diagram (Fig. 12) thus obtained, reveals that 2 nd , 3 rd and 6 th modes at 7.7Hz, 13.7Hz and 26.1Hz are not identified. The stabilization histogram generated by the proposed algorithm obtained after cancellation of spurious modes on 4 th iteration is shown in Fig. 13. The distinct peaks allow the first seven modes to be easily distinguished from the spurious modes.

VI. CONCLUSIONS

Identifiability of modal parameters is enhanced by exciting the structure with a well defined input signal. Spurious modes are eliminated by changing sampling rate of measured output signal and by combining the results on a single stabilization histogram. The estimated results have show that the proposed algorithms not only allowed automatic selection of excitation signal in small number of iterations but also removed the spurious modes from the generated results. The present approach was applied using a SubID method, which is based on SVD of Hankel matrix, which is computationally complex. For a real-time SHM applications, the underlining identification algorithm should have reduced computational complexity. Next steps in this study are to observe the performance of the proposed algorithms using recursive SubID methods and to analyze their robustness in generating modal parameter residuals for SHM applications.


Fei Zheng

Stéphane Bonnet

Emma Villeneuve

Maeva Doron

Aurore Lepecq

Florence Forbes

Unannounced Meal Detection for Artificial Pancreas Systems Using Extended Isolation Forest

This study aims at developing an unannounced meal detection method for artificial pancreas, based on a recent extension of Isolation Forest. The proposed method makes use of features accounting for individual Continuous Glucose Monitoring (CGM) profiles and benefits from a two-threshold decision rule detection. The advantage of using Extended Isolation Forest (EIF) instead of the standard one is supported by experiments on data from virtual diabetic patients, showing good detection accuracy with acceptable detection delays.

I. INTRODUCTION

Currently, the most used Artificial Pancreas (AP) system implements an hybrid closed-loop requiring patients with Type 1 Diabetes (T1D) to declare both the time and the amount of their carbohydrates (CHO) intakes. Consequently, such systems are prone to CHO underestimations due to missing declarations and may lead to inadequate insulin dosing, and expose the subjects to hyperglycemia risk [START_REF] Gingras | The challenges of achieving postprandial glucose control using closed-loop systems in patients with type 1 diabetes[END_REF]. Unannounced meal detection and quantification is thus an essential task to reduce such risk as well as to achieve an accurate fully automated closed-loop AP system. Studies have been devoted to unannounced meal detection with a majority of them implementing physical model-based detection methods. The minimal model [START_REF] Bergman | Quantitative estimation of insulin sensitivity[END_REF] or Medtronic Virtual Patient (MVP) model [START_REF] Stocker | Virtual patient software system for educating and treating individuals with diabetes[END_REF], [START_REF] Kanderian | Identification of intraday metabolic profiles during closed-loop glucose control in individuals with type 1 diabetes[END_REF] are modified for the application of different prediction algorithms, for example Unscented Kalman Filter [START_REF] Turksoy | Meal detection in patients with type 1 diabetes: a new module for the multivariable adaptive artificial pancreas control system[END_REF], [START_REF] Ramkissoon | Unannounced meals in the artificial pancreas: detection using continuous glucose monitoring[END_REF], Moving Horizon Estimation [START_REF] Kölle | Pattern recognition reveals characteristic postprandial glucose changes: Non-individualized meal detection in diabetes mellitus type 1[END_REF] or Sliding Mode Observer [START_REF] Sala-Mira | Sliding-mode disturbance observers for an artificial pancreas without meal announcement[END_REF] to estimate the states (specifically glycemia and meal appearance rate). Some linear state-space models are adjusted also for describing the pharmacodynamic response of insulin and glucose appearance, so that Kalman Filter [START_REF] Fathi | An unannounced meal detection module for artificial pancreas control systems[END_REF] or Variable State Dimension [START_REF] Xie | A variable state dimension approach to meal detection and meal size estimation: in silico evaluation through basalbolus insulin therapy for type 1 diabetes[END_REF] can be applied for the states estimation. Then, detections are made by comparing the residuals between the estimated and observed states to a threshold [START_REF] Sala-Mira | Sliding-mode disturbance observers for an artificial pancreas without meal announcement[END_REF], by calculating their cross-covariance [START_REF] Ramkissoon | Unannounced meals in the artificial pancreas: detection using continuous glucose monitoring[END_REF], or by applying further decisions rules, for example hypothesis testing [START_REF] Fathi | An unannounced meal detection module for artificial pancreas control systems[END_REF] or statistical analysis [START_REF] Kölle | Pattern recognition reveals characteristic postprandial glucose changes: Non-individualized meal detection in diabetes mellitus type 1[END_REF].

Besides the model-based detection approaches, data-driven methods have the advantage that their performance is much less dependent on the identification accuracy of the applied physical model. Three data-driven methods are studied in [START_REF] Meneghetti | Data-driven anomaly recognition for unsupervised model-free fault detection in artificial pancreas[END_REF] for insulin pump fault detection, showing that Isolation Forest (IF) [START_REF] Liu | Isolation forest[END_REF] outperforms the other two density-based outlier detection methods. Inspired by [START_REF] Meneghetti | Data-driven anomaly recognition for unsupervised model-free fault detection in artificial pancreas[END_REF], our work adapts the application of IF to unannounced meal detection. We propose suitable features for individual Continuous Glucose Monitoring (CGM) profiles, and a two-threshold-based rule for risk evaluation and detection. More importantly, the replacement of the standard IF by its recent generalization named Extended Isolation Forest (EIF) [START_REF] Hariri | Extended isolation forest[END_REF] is explored and shows promising performance on a virtual data set for this specific unannounced meal detection task.

The remaining of the paper is organized as follows: Section II presents the proposed EIF based detection method, Section III shows the superiority of EIF compared to standard IF and evaluates the performance of the method on a virtual data set; finally, Section IV concludes the work and draws some perspectives.

II. METHODS

A. Isolation Forest & Extended Isolation Forest

Isolation forest (IF) is an approach to detect anomalies that explicitly isolates anomalies rather than profiles normal instances [START_REF] Liu | Isolation forest[END_REF].

Let X = {x 1 , • • • , x N } denote the training set, where x n ∈ R d . X ψ [m]
denotes a randomly sampled subset of size ψ from X with m ∈ {1, • • • , M } represents the subset size. An isolation tree (iTree) recursively divides X ψ [m] by randomly selecting at each node an attribute among d features and a split value. The process is stopped for a node when it contains only one instance or when the instances are with the same value.

We denote by h (x) the so-called "path length" of x which is the number of edges that x traverses from the root node to the external node of an iTree. In practice, the path length is limited to the depth l = ceil(log 2 ψ). Then, the Anomaly Score (AS) is defined by:

s (x) = 2 -E[h(x)] c(N ) . (1)
c (N ) is the theoretical average path length which is only related to the training set size N [START_REF] Liu | Isolation forest[END_REF]. E [h (x)] represents the expectation of path length of instance x, which can be approached by averaging over M iTrees. The AS takes its value in (0, 1]. A score close to one indicates that the corresponding instance x is easy to isolate and could then be considered as an outlier.

As a variant of IF, Extended Isolation Forest (EIF) creates nodes by randomly selecting an hyper-plane in the feature space. It can thus provide relatively more consistent anomaly scores, as well as more accurate detection of the feature structure [START_REF] Hariri | Extended isolation forest[END_REF]. When applying EIF, different extension levels can be chosen related to the number of feature axes that hyper-planes intersect, ranging from 1 (intersect 1 feature,

B. Features for Unannounced Meal Detection

Considering the unannounced meal detection problem, inspired by [START_REF] Meneghetti | Data-driven anomaly recognition for unsupervised model-free fault detection in artificial pancreas[END_REF], d = 3 features are chosen such that x = [G p (t) , DCOB (t) , DIOB (t)]. G p (t) represents plasma glucose which is approximated by CGM measurement, while DCOB, DIOB are defined as follows:

DCOB(t) = α exp γ Ġp (t) ωCOB (t) + 1 , DIOB(t) = β exp γ Ġp (t) υIOB (t) + 1 .
(2) Ġp (t) represents the time derivative of G p (t) [mmol/L], COB (t) [mmol] and IOB (t) [mU] are CHO on board and insulin on board respectively. COB and IOB are calculated using declared CHO (if any) and injected insulin through a modified MVP model (see Section V for details) identified on a training set. Doing so, the MVP model parameters as well as the features can adapt to specific profile. The description of the virtual patient simulator is in Section III. The hyper parameters are manually tuned and set as: α = β = 15, γ = 5, ω = 8 × 10 -4 and υ = 10 -4 .

An example of features from a virtual patient is illustrated in Fig. 1. It is worth noticing that the computation of Ġp (t) requires tke knowledge of G p (t + 1), if 2nd order accurate central difference is applied as in this work, which leads to a minimum detection delay of one sample interval, = 5 min.

C. Risk levels & Detection Decision

Outlier detection is usually performed using a single threshold. We propose a more reliable decision rule using two thresholds related to two unannounced meal risk levels, low and high, respectively. They are defined as: We then define the low risk and high risk thresholds using parameters η = {η l , η h }. An alarm of unannounced meal is then raised when either: 1) one high risk appears (AS > thr(η h )); 2) or three successive low risks are observed (AS > thr(η l )). See Fig. 3 for an illustration. To summarize, the meal detection process for individual subject consists of a training and testing step as follows:

thr(η) = F -1 (η) 0 ≤ η ≤ 1 1 -(2 -η) 1 -F -1 (1) 1 < η ≤ 2 , (3)
• Training: 1) a MVP model is first identified using a training dataset. 2) EIF is trained on the training feature set and AS thresholds are defined. • Testing: 1) The MVP initial state of the tested scenario is estimated; 2) for each coming sample, its AS is computed and the announced meal risk is evaluated with the following appropriate detection decision. For evaluation, the regular detection window is set to 3h (maximum detection delay) except when two meals are less than 3h apart (see irregular window in Fig. 3).

III. EXPERIMENTAL RESULTS

A. Virtual Patient Data Set

The proposed method is tested on four virtual patients simulated from SIM-HOV1 based on Hovorka model [START_REF] Hovorka | Partitioning glucose distribution/transport, disposal, and endogenous production during IVGTT[END_REF]. The data includes CGM measures, insulin doses and announced meals with a sample rate of 5 min. In addition, SIM-HOV considers that the data is contaminated by a Gaussian noise with a standard deviation (sd) of 0.012 × G p (0) mmol/L, where G p (0) is the initial glycemia value. The patients have various physiologies, such as different insulin sensitivity and digestion speed. Simulated scenarios correspond to an average of 3.4 meals per day, with time and size of the meals randomly sampled from the Diabeloop WP7 data set2 (median meal size: 40g CHO and interquartile 35g). Thus the virtual patients have similar eating behaviour to the real patients. Their glycemia was continuously regulated in closed-loop by an AP Controller for Simulation (APCS) designed by Algorithms Teams from Diabeloop S.A and CEA. APCS adjusts in real-time the insulin bolus and basal infusions according to a certain strategy in order to maintain glycemia in a healthy range. As we can see in Fig. 1, some compensation insulin doses are delivered after meals. However, APCS does not realise any meal detection.

In the experiment, for each virtual patient, a 5-day scenario (with all meals announced) of data is simulated for training, and two scenarios of 20 days are available for testing. One scenario includes 50% of declared meals, to mimick the case where meals are irregularly announced by patients. Another scenario has 0% meal declared to mimick a fully automatic meal detection case.

B. Performance

The procedure described in Section II-C is applied for unannounced meal detection to each virtual subject. The MVP model is identified through a non-linear least-square method. Then, EIF with a certain extension level is trained on the 5-day scenario training set. The parameters of EIF are set empirically to ψ = 256 and M = 100. Various thresholds are tried with η l in [0.6, 1.2] and η h in [η l , 1.2] with a 0.05 step. When η l = η h , the decision is degraded to the commonly used single threshold case. Testing is conducted on the two scenarios with 50% and 0% declared meal. The initial state of the MVP model is estimated using the initial 5-hour data. The detection is then processed "real-time" like, sample by sample, calculating the AS value through trained EIF, measuring the risk and taking a subsequent decision, as illustrated in Fig. 3.

ROC surfaces of the detection performance for different choices of η are illustrated in Fig. 4. From it we can conclude that under the proposed detection framework, EIF with the highest extension level outperforms the other extension levels as expected. The thresholds chosen by f 1 score [START_REF] Hand | A note on using the F-measure for evaluating record linkage algorithms[END_REF] in the test of 50% and 0% meal declaration are close to each other and both with η l < η h . This shows the interest of the two threshold decision.

To further assess the detection performance, we detail the fully extended IF case (extension level 2) with η = [0.950, 1.000]. Table I reports the performance matrix, which shows for both 50% and 0% meal declaration cases, that EIF reaches a high true positive detection rate around 90%. The false alarm rate of the all unannounced meal case is 11.47% while it is 6.2% for the 50% declaration case. This is probably due to the continuous abnormal behavior affected by dense unannounced meals. Around 39% announced meals are detected as unannounced, which may mean that, measured by the trained EIF, the delivered insulin regulated by APCS is not adequate to compensate the effect of the announced meal. Indeed, the proposed method is tested under a harder condition than the pure open-loop case, since less evidence is observed thanks to the compensation from APCS. The standard deviations shown in brackets confirm an acceptable performance variance in between subjects.

Table II reports the delays of detection which are within 40 min on average, and risk notification are within 31 min for all cases. It happens for the virtual patients that some unannounced meals raise risks that maintain longer than their irregular detection windows, and are detected as an alarm of the following announced meal. This phenomena leads to a shorter detection delay concerning announced meals as observed in Table II.

The detection result from a virtual patient with 50% declared meals is illustrated on its glycemia record in the bottom figure of Fig. 5. The corresponding meal and insulin intake, estimated AS values and risk thresholds are depicted in the figure above. The good performance of the proposed detection method can be visualized. Running on a 3.6 GHz CPU, the MVP identification and EIF training take 197 s, then the risk evaluation of each sample takes only 0.005 s.

IV. CONCLUSIONS & PERSPECTIVES

An Extended Isolation Forest (EIF) based method is proposed to detect unannounced meals for artificial pancreas. The method captures differences in individual profiles, with real-time applicability. The detection performance is illustrated on virtual patients data simulated from SIM-HOV. It shows the superiority of using EIF over the standard Isolation Forest (IF) and the interest brought by the introduced two threshold decision rule.

Future work includes: 1) MVP update for adapting the physiological dynamics of a subject; 2) safe thresholds decision strategy for real patients individually; 3) meal size estimation and compensation bolus calculation for fully closed-loop AP, possibly using the estimated AS from EIF.

V. APPENDIX A. Modified Medtronic Virtual Patient Model

The modified MVP model [START_REF] Stocker | Virtual patient software system for educating and treating individuals with diabetes[END_REF], [START_REF] Kanderian | Identification of intraday metabolic profiles during closed-loop glucose control in individuals with type 1 diabetes[END_REF] with a complexity trades off between minimal models and maximal ones is defined by: CHO absorption Based on the MVP model, COB(t) = D 1 (t) + D 2 (t) and IOB(t) = S 1 (t) + S 2 (t) for features in [START_REF] Bergman | Quantitative estimation of insulin sensitivity[END_REF].

Ḋ1 (t) = A G U m (t) -D 1 (t)/τ D Ḋ2 (t) = D 1 (t)/τ
The estimated parameters of the subject in Fig. 5

Fig. 1 :

 1 Fig. 1: Feature time waveform. Green area in top figure shows the ideal glycemia range. Gray stems in middle figure show unannounced meal and orange stems are announced ones. Insulin doses (basal + bolus) are displayed by blue stems in bottom figure.

Fig. 2 :Fig. 3 : 1 ,

 231 Fig. 2: Anomaly Score (AS) threshold.

Fig. 4 :

 4 Fig. 4: ROC surfaces for different IF extension levels with respect to different thresholds (annotation shows the value of η which gives the best f 1 score).

Fig. 5 :

 5 Fig. 5: An example of detection result. Risks are shown on glycemia (bottom), with corresponding meal [g] and insulin (5× original value) [U/min] information, AS and AS thresholds (top).

  D -D 2 (t)/τ D insulin absorption    Ṡ1 (t) = U i (t) -S 1 (t)/τ S Ṡ2 (t) = S 1 (t)/τ S -S 2 (t)/τ S İP (t) = S 2 (t)/(V I τ S ) -k e I P (t) İeff (t) = -p 2 I ef f (t) + p 2 S I I p (t) Ġp (t) = -[p 1 + I ef f (t)] G p (t) + p 1 G b + D 2 (t)/(V G τ D ) The model considers 2 inputs, CHO intake U m (t) [mmol/min],insulin delivery U i (t) [mU/min], and 7 states, amount of glucose in remote compartments D 1 (t), D 2 (t) [mmol], amount of insulin in remote compartments S 1 (t), S 2 (t) [mU], plasma insulin concentration I p (t) [mU/L], insulin effect I ef f (t) [1/min] and plasma glucose G p (t) [mmol/L]. The 10 parameters are: CHO utilization A G , CHO and insulin absorption delay τ D , τ S [min], insulin elimination rate k e [1/min], insulin and glucose distribution volume V I , V G [L], insulin sensitivity S I [L/mU/min], glucose effectiveness at zero insulin p 1 [1/min], rate of disappearance of insulin effect p 2 [1/min], and basal glucose G b [mmol/L].

  are: A G = 0.86, τ D = 51.7 min, τ S = 97.6 min, k e = 0.09 /min, V I = 12.3 L, V G = 6.5 L, S I = 3.96 × 10 -4 L/mU/min, p 1 = 0.026 /min, p 2 = 0.01 /min, G b = 8.62 mmol/L.

TABLE I :

 I Performance (detection rate average (sd)) for the fully extended IF with η = [0.95, 1.00].

	Event	Unannounced meal	Announced meal	No event
	50% declared	90.8% (8.0%)	38.9% (17.5%)	6.2% (4.1%)
	0% declared	90.0% (2.1%)	-	11.47% (6.8%)

TABLE II :

 II Detection delay (average (sd) min)

	Delay	Unannounced meal 50% declared 0% declared	Announced meal
	Detection	39.0 (22.8)	35.1 (16)	31.1 (16.6)
	risk notification	30.2 (22.3)	26.3 (15)	22.4 (15.4)

INRIA, CNRS, LJK, Univ. Grenoble Alpes, Grenoble INP, F-38000 Grenoble, France.

Univ. Grenoble Alpes, CEA, LETI, DTBS, LS2P, F-38000 Grenoble, France.

SIM-HOV is a Python simulator developed by Diabeloop S.A.CONFIDENTIAL. Limited circulation. For review only.Manuscript 952 submitted to

42nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) in

Diabeloop WP7 study[START_REF] Benhamou | Closed-loop insulin delivery in adults with type 1 diabetes in real-life conditions: a 12-week multicentre, open-label randomised controlled crossover trial[END_REF] is a 12-weeks nationwide (France) randomized crossover trial study concerning 63 T1D patients.

Manuscript 952 submitted to 2020 42nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) in

ACKNOWLEDGMENT

The authors would like to thank Diabeloop S.A for providing the virtual patient data set.


Anne Daubresse

Mohend Evelyne Prat

D Kosiachevskyi
email: kosiachevskyi.dima@gmail.com

K Abahri

M Chaouche

Assessment Of The Hygrothermal, Microstructural And Chemical Evolution Of A Hemp-Based Cementitious Mortar Under Etics Total Weathering Aging Protocol

Keywords: Hemp mortar, ETICS, rain cycles, heat-cold cycles, microstructure, thermal conductivity, moisture buffer value, vapor permeability

Introduction

The building industry is a very energy consuming sector and is responsible for a large share of the world's CO2 emission. Besides large-scale activity, like transportation, it is due to the use of the energy intensive materials like cement, bricks, steel and others. Moreover, the construction of the buildings represents a small part of the energy consumption in their total life cycle. As for the use of residential and tertiary buildings, it consumed 45% of total final energy and emitted 22% of CO2 in France in 2016 [1]. The only viable solution for developing the sustainable building production system is to minimize the net CO2 emission and use alternative low-energy building materials and more energy efficient technologies. That is why the use of bio-based materials based on different types of agricultural products gained a particular interest in France and around the world in order to solve the problem of CO2 emission. In addition to their low carbon footprint, the use of these materials is promoted due to their hygrothermal properties and the availability of the raw materials [START_REF] Walker | Mechanical properties and durability of hemp-lime concretes[END_REF].

Moreover, the necessity of more thermal insulation to get better energy performance of buildings leads to the problem of air quality, and hemp mortar is the material that can regulate indoor conditions by reducing daily humidity variations. Different natural fibers in building insulation materials are well studied, especially in regard to mechanical, thermal and hygroscopic properties. There are a lot of works on the properties of wood [START_REF] Hachem | Fiber's hygromorphic effect on thermal conductivity of wooden fibrous insulation characterized by X-ray tomography[END_REF][START_REF] El Hachem | Original experimental and numerical approach for prediction of the microscopic hygro-mechanical behavior of spruce wood[END_REF], straw [START_REF] Ren | Study on the hygrothermal properties of a Chinese solar greenhouse with a straw block north wall[END_REF], cork [START_REF] Panesar | The mechanical, transport and thermal properties of mortar and concrete containing waste cork[END_REF][START_REF] Pacheco Menor | Granulated cork with bark characterised as environment-friendly lightweight aggregate for cement based materials[END_REF][START_REF] Jerónimo | Hydraulic lime mortars incorporating micro cork granules with antifungal properties[END_REF], jute [START_REF] Garikapati | Mechanical behavior of flax-lime concrete blocks made of waste flax shives and lime binder reinforced with jute fabric[END_REF], and hemp [START_REF] Walker | Mechanical properties and durability of hemp-lime concretes[END_REF]. However, the use of bio-based building materials is hampered by the lack of a database relating to its durability.

In this paper, the investigated material is hemp mortar which is used for its thermal, acoustic and hydric performances [START_REF] Bennai | Experimental characterization of thermal and hygric properties of hemp concrete with consideration of the material age evolution[END_REF]. Nevertheless, hemp mortar is a material of complex, heterogeneous and anisotropic composition. Like all composites based on natural fibers, hemp mortar must be investigated taking in account its multiscale nature related to its strong heterogeneity: from the macroand microscale where hemp mortar is represented as a constitutive part of the envelope or as a multiphase material with interphases between hemp particles and a cementitious matrix to the nanoscale of organic cell walls tissues or cement hydrates. Another feature that improves the performances (for example thermal properties) of the insulating hemp-based mortar is that it has a double porosity: intra-(hemp particles porosity) and inter-particular (hemp/matrix interface).

Consequently, such compositional and structural complexity complicates the investigation of hemp mortar when talking about the durability.

Another important point that must be taken into account is the hygroscopicity of hemp mortar and its sensitivity to thermal and hygric solicitations and to the environmental conditions of its life cycle. In fact, directly after the hemp mortar cast at early-age, cracks appear at the interfaces between the hemp particles and the cementitious matrix. Later, these cracks increase due to dimensional variations of hemp particles caused by the relative humidity solicitations. This modifies the hemp mortar characteristics and leads to the degradation of the material [START_REF] Bennai | Influence of hydric solicitations on the morphological behavior of hemp concrete[END_REF][START_REF] Bennai | Microscopic hydric characterization of hemp concrete by X-ray microtomography and digital volume correlation[END_REF]. The evolution of the microstructural properties must be investigated in order to understand the behavior and durability of the bio-based composite as regards to the hygrothermal stresses.

When talking about the durability of hemp mortars, few researchers studied the hygrothermal stresses impact on the hemp mortar properties. For example, Sentenac et al. [START_REF] Sentenac | Investigation on the performance and durability of treated hemp concrete with linseed oil[END_REF] studied the performance and durability of treated hemp mortar. To estimate the evolution of the different properties and physical characteristics, the full immersion test was used. The conditioning represents 12 to 15 cycles of water immersion at 20 °C during 48 hours and drying at 40 °C during 48 °C depending on the samples. This immersion protocol has often been criticized because it is highly aggressive and also not very representative of the real ageing conditions of the materials.

Other researchers used accelerated wetting and drying cycles to investigate the aging of hemp mortars. For example, Marceau et al. [START_REF] Marceau | Influence of accelerated aging on the properties of hemp concretes[END_REF] adopted the wet/dry cycles protocol at constant temperature of 30°C (each cycle lasts 8 hours). Each conditioning consists of 8 cycles: 4 cycles of 40% / 90% RH and 4 cycles of 40% / 98% RH. They studied the changes in the microstructure (density evolution, hemp particles / cementitious matrix interface) and physical properties (thermal and acoustical performance). Delannoy et al. [START_REF] Delannoy | Durability of hemp concretes exposed to accelerated environmental aging[END_REF] proposed another accelerated aging protocol consisting of 5 days at 98% RH and 2 days of 40% RH constantly at 30 °C. Their study examines the evolution of the different properties of hemp mortar and its chemical composition.

In all of the proposed literature works, non-standardized accelerated aging protocols were adopted that are not representative of the real-life cycle conditions of hemp mortar insulation. This is explained by the total water immersion in the case of the protocol of Sentenac et al. [START_REF] Sentenac | Investigation on the performance and durability of treated hemp concrete with linseed oil[END_REF] or high relative humidity in the works of Marceau et al. [START_REF] Marceau | Influence of accelerated aging on the properties of hemp concretes[END_REF] and Delannoy et al. [START_REF] Delannoy | Durability of hemp concretes exposed to accelerated environmental aging[END_REF], that do not appear under standard conditions of application and use recommended by Piot et al. [START_REF] Piot | Study of a hempcrete wall exposed to outdoor climate: Effects of the coating[END_REF]. That is why the main objectives of this paper are to adapt and propose the accelerated aging protocol basing on the standardized one used for External Thermal Insulation Composite Systems (ETICS) with rendering [START_REF] Eota | Guideline for European technical approval of external thermal insulation composite systems (ETICS) with rendering[END_REF] and to investigate the microstructure and properties evolution of hemp mortar due to proposed aging cycles. The main advantages of the proposed protocol are the use of heat-rain and heat-cold cycles and protection against the direct water access that allow to represent real conditions of use of biobased insulation materials with rendering. Experimentally, the evaluation of the hemp mortar characteristics was conducted before and after weathering cycles.

Experimental procedure

In this section, the studied material, the accelerated weathering aging cycles conditions and the experimental characterization methods are presented. The studied hemp mortar can be used for both external and internal insulation. Characterization of hemp mortar was conducted on two benchmarks: before and after the exposure to aging protocol. This allows us to compare the results and study the influence of the conditioning on hemp mortar microstructure, chemical composition and hygrothermal properties.

Materials

The hemp shives of France origin were used for manufacturing the hemp mortar samples. Their visual aspect and granulometric curve are presented in Figures 1 and2 respectively. The binder represents a commercial mix designed by ParexGroup SA. It is composed of hydrated lime, Portland cement and various additions to achieve the rheology suitable for spray application. Chemical composition of mineral binder is presented in Table 1. The hemp mortar proportion ratios are represented in the Table 2. This mix is suitable for material application by spraying. This requires special attention to the rheology that allows the industrial application of insulating hemp mortar. The studied samples were manufactured in laboratory using the laboratory mixer. The mixing was done in two stages: firstly, hemp shives were mixed with water during 1 minute at 140 rpm, secondly, the binder was added and it was mixed during 5 minutes at 140 rpm. After their laboratory manufacture, the specimens were dried for 28 days at laboratory conditions (23 ± 2 °C, 50 ± 5% RH) as shown in Figure 3. The adopted sample sizes depend on experimental procedure and are presented for every type of test in the paragraph 2.3.

Accelerated weathering aging cycles

After the 28 days of drying, the hemp mortar samples were conditioned in two different ambiences. Total weathering protocol is represented in Figure 4-A as is explained in the following:

• The first batch of specimens was subjected to the laboratory conditions (23 ± 2 °C, 50 ± 5% RH) for the whole period of the aging cycles. These samples were used later as a reference.

• The second part of the specimens was exposed to the accelerated aging cycles that consist of two parts: heat -rain and heat -cold cycles that are described below:

Heat -rain cycles: Samples are subjected to a series of 80 cycles, comprising the following phases: 1-heating to 70°C (rise for 1 hour) and maintaining at (70 ± 5)°C and 10 to 30% RH for 2 hours, 2-spraying for 1 hour (water temperature at + 15 ± 5 °C, and amount of water 1 l/m² min), 3 -leave for 2 hours (drainage). These cycles are represented in Figure 4-B.

Heat-cold cycles: After a subsequent conditioning of at least 48 hours at temperatures between 10 and 25 °C and a minimum relative humidity of 50 %, the same test rig is exposed to 5 heat/cold cycles of 24 hours comprising the following phases: 1-exposure to (50 ± 5)°C (rise for 1 hour) and maximum 30% RH for 7 hours (total of 8 hours), 2exposure to (-20 ± 5)°C (fall for 2 hours) for 14 hours (total of 16 hours) (Figure 4-C).

Samples of two batches (those used as reference and those subjected to the accelerated aging cycles) were characterized and compared later to assess the aging effect. As hemp mortar is used as insulation material with rendering, it is protected against direct access of water during the accelerated aging cycles. That is why in this study, we used a roofing to avoid direct spraying of water on the samples, as shown in Figure 5. The overall protocol schema is presented in Figure 6.

Microstructural characteristics and chemical analysis

The total porosity tests In order to measure the total porosity, the water porosity test was used accordingly the NFP 18-459 standard [START_REF]NF P18-459 -Concrete -testing hardened concrete -testing porosity and density[END_REF]. It is based on the hemp mortar saturation in water under vacuum. Then different weightings were done to calculate the total porosity of the samples. The scales sensitivity was 0.01 g. One of the weightings was the weighting of wet sample. As hemp mortar becomes fragile after wetting, we used a cotton tissue in order to protect the wet sample from damage (Figure 7). Three specimens (15x15x7 cm 3 ) were tested for each time of conditioning to validate the repeatability of the test.

Scanning Electron Microscopy (W-SEM) microstructural analysis

Scanning electron microscope (W-SEM) analyses were performed by using a Hitachi S-3400N

device. The surface of the hemp mortar and the interface zone between cementitious matrix and hemp particles were observed and analyzed before and after aging. Particular attention has been devoted to the identification of the same observation area for these two states.

Micrographs were recorded by using secondary electrons detectors at an acceleration voltage of 15 kV, in high vacuum mode (pressure <1Pa). Previously to observations, the specimens were covered by carbon in order to improve the image quality.

PH level evolution

There are a lot of different methods of the identification of a pH level. Some authors identify a pH level of interstitial solution by compression that allow to get the most reliable results as we don't modify the pH level adding water [START_REF] Alonso | Low-pH Cementitious Materials Design and Characterisation[END_REF]. Nevertheless, this method is a heavy procedure that requires special equipment and a large quantity of cement. Others propose another method -of cementitious suspensions [START_REF] Codina | Les bétons bas pH -Formulation, caractérisation et étude à long terme[END_REF]. It consists of suspending of some part of the mortar in water and measuring a pH level of obtained suspension. At the same time, the use of suspensions implies problems of dissolution of the solid phases but especially of dilution, moreover if talking about hygroscopic composite material. That is why in this work a surface electrode WTW Sentix Sur was used [START_REF] Tran | Influence of the intrinsic characteristics of mortars on their biofouling by pigmented organisms: Comparison between laboratory and field-scale experiments[END_REF]. It allowed us to get a reliable result and especially on different parts of samples to evaluate a chemical degradation of the mineral matrix.

In order to better understand the degradation of the samples, we have investigated two different areas for the measurement of pH: the surface and the middle of the sample. The idea is to understand whether or not the degradation caused by hygrothermal stresses has actually propagated inside the sample or not. In fact, most of the works in the literature studding hemp mortars do not pay attention to the position of pH measurement in the sample. Nevertheless, it is important as the chemical degradation starts from the surface and is a long process.

For this test the specimens of the same composition, made with the same raw materials, of the size 16x4x4 cm 3 were used. For each value three measures were performed at different positions on the specimen. As the pH level measurements are possible just in solution, the constant volume of demineralized water (0.2 ml) was deposed on the studied surfaces. Then, the electrode was kept on the sample surface, so that the electrode's membrane remains in contact with the hemp mortar (Figure 8). At equilibrium, the WTW Sentix Sur displays a constant value of PH which will be subsequently recorded. For the inner pH level measures the protocol is the same, but the measurements were conducted in the middle of the specimen.

Figure 8. Measurement of the surface pH level using a WTW Sentix Sur electrode

X-ray diffraction (XRD) crystallography analysis

Prior to XRD crystallography analysis, hemp mortar samples of the size 16x4x4 cm 3 were grinded and sieved to remove organic fibers. This step allowed us to decrease the influence of the dilution on the final result as in this case we analyze only mineral matrix. To link the results of chemical composition to the pH level, the crystallography at the surface and in the middle was analyzed.

Free water was removed using the acetone. Then, the specimens were dried at 50 °C in order not to modify the composition of the sample. XRD crystallography patterns were obtained using Cu Kα anode tube (λ= 1.54182 Å) radiation at 40 kV and 20 mA with a Bruker D8 diffractometer (Bruker, Karlsruhe, Germany). The diffractometer was scanned from 5° to 60° (2θ°) in step size of 0.015° and the counting time per step was 1.2 s.

Hygrothermal properties evolution

The water vapor permeability (WVP) (wet cup method)

The water vapor permeability is the ability of a material to transfer or retain water vapor due to a difference in the partial pressure of water vapor at the same atmospheric pressure on both sides of the material. It was identified according to the NF EN ISO-12572 standard using the wet cup method [START_REF]Hygrothermal performance of building materials and products -Determination of water vapour transmission[END_REF]. Firstly, specimens of diameter 11 cm and thickness of 2.4 cm (3 samples for each test) were dried at 23 °C, 50% RH during 28 days. Then, they were sealed on cups containing a saturated salt solution that ensures constant relative humidity in the cup (Figure 9). In the case of wet cup, ammonium dihydrogen phosphate was used that maintains the relative humidity at 93±1%. The cups are conditioned in a constant ambiance of 23 °C and 50% RH. The water vapor permeability value was determined in function of the vapor transmission through the specimen. To do so, the mass monitoring was conducted every day during 18 days for each cup to calculate the mass variation velocity:

∆ = • • (1)
where, ∆ is a mass variation velocity (kg/s); and are masses of the cups with hemp mortar samples at the time and respectively (kg); and are the successive times of weighting (s).

Then, the water vapor flux density was calculated:

g = (2)
where, G is an average value of five successive measurements of ∆ and A is an exposed area (m2).

The water vapor permeance was calculated as

W = •∆ (3)
where, ∆p is difference of water vapor during the test.

The water vapor resistance is an inverse value to permeance:

Z = (4)
In fact, there are several resistances in the experimental montage:

= + + + (5)
where, is an experimental calculated value of resistance, is a water vapor resistance of air, and are the water vapor resistances that take into account the interface diffusion at the sample/air and air/saline solution interfaces respectively.

Finally, the water vapor permeability value was calculated:

δ = (6)
where, d is a thickness of our sample (mm).

The thermal conductivity tests

For measuring thermal conductivity, there are two types of experimental methods: steady-state and non-steady-state ones. When talking about measurement setups, we can distinguish four main types: the guarded hot plate (GHP), the heat-flow meter (HFM), the hot wire, and laser flash diffusivity [START_REF] Yüksel | The Review of Some Commonly Used Methods and Techniques to Measure the Thermal Conductivity of Insulation Materials. Insulation Materials in Context of Sustainability[END_REF]. In this work we used the method of hot wire as it is well adapted for bio-based materials. Thermal conductivity measurements were taken on three different pairs of the 15x15x7 cm 3 samples under transitory conditions using the device CT-Metre (Figure 10). Before measurements, these samples were conditioned at 23 °C and 50 % RH until a constant mass, according to [START_REF][END_REF].

Measurements were made according to the device precautions on the measurements. The used protocol consists of placing a probe between two samples keeping the maximal parallelism of the surfaces to avoid the distortion of results. The thermal conductivity values were calculated by the CT-Metre automatically, coefficients of correlation of presented results were 0.999.

Moisture buffer value (MBV) tests

The moisture buffer value represents the ability of the material to mitigate the relative humidity changes in ambient air. It refers to the amount of water that could be absorbed/released when exposed to relative humidity change. In this paper, the Nordtest protocol was used to investigate the MBV of the studied hemp mortar. During this test the samples of hemp shives placed in the aluminum recipient of 10x10x3 cm 3 size and hemp mortar samples of 15x15x7 cm 3 size were used. To do this, variable cycles of relative humidity of 75% (for 8 hours) and 33% (for 16 hours) were applied. Samples were weathered at laboratory conditions (23 ± 2 °C, 50 ± 5% RH) beforehand. The MBV value was calculated:

! = ∆ "•($% &'( )$% &*+ ) (7)
where ∆m: sample mass variation [kg],

S: exposed area of sample [m 2 ],

-. /0 : maximal relative humidity [%],

-. 12 : minimal relative humidity [%].

According to [START_REF]BYG DTU R-126 Moisture buffering of building materials[END_REF], building materials are classified using the MBV range represented in Figure 11.

Results and discussions

In this section, the assessment of the microstructure, the chemical composition and its hygrothermal performance evolution under the accelerated aging cycles are presented.

Microstructural characteristics and chemical analysis

In order to study the microstructure evolution, firstly the total porosity was analyzed before and after accelerated aging cycles. Then, the interface between cement matrix and hemp particles was observed and the evolution of the chemical composition was studied. The analysis of these results was complemented by SEM.

We note that the studied hemp mortar has a porosity of 85.2% which is more important than the other mortars, presented in literature [START_REF] Bennai | Experimental characterization of thermal and hygric properties of hemp concrete with consideration of the material age evolution[END_REF]. This is explained by the organic fibers porosity and the properties of the cementitious matrix that ensure good insulation properties.

According to the results given in Table 3, weathering aging cycles affect very slightly the total porosity, increasing the value from 85.2% to 85.8%. This result remains in the range of a standard deviation. Such slight evolution of porosity could be explained by the dimensional variation of the hemp mortar due to hygrothermal solicitations. In fact, hygrothermal stresses cause swelling or shrinkage of the organic fiber causing cracking of the interface and, consequently, an increase in porosity. At the same time there is another effect -a carbonation of mineral matrix and such porosity change could be explained by competition between these two phenomena. Both effects are investigated in detail. In order to study the hemp mortar cracking before and after proposed weathering aging cycles, the W-SEM observation of the hemp particles / cementitious matrix interface was conducted. To do this, the surface areas of 2 samples for each type of conditioning were observed. For each sample, 10 interfaces were examined and the crack sizes were measured in 10 spots all along the crack.

example of the image before and after aging are present in Figure 12. We notice the cracking of the The carbonatation of the cementitious matrix was investigated by the X-ray diffraction crystallography (XRD). It was conducted for the samples before and after the weathering aging cycles which allowed estimate the chemical composition evolution. For each specimen the analysis of the composition on the surface and in the middle of the studied samples was performed.

The XRD patterns presented in Figure 13 showed the main mineralogical compounds of the hemp mortar as Mg4Si6O15(OH)2.6H2O (sepiolite), Ca6Al2(SO4)3(OH)12.26H2O (ettringite), Ca(OH)2

(portlandite), CaCO3 (calcite), SiO2 (quartz), Ca3SiO5 (tricalcium silicate), Ca2SiO4 (dicalcium silicate) and Ca3Al2O6 (tricalcium aluminate). Figures 13 (a) and 13 (c) present the crystallographic analysis of the hemp mortar specimen before the weathering aging cycles in the middle and on the surface, respectively. We note that after one month of drying at 23 °C and 50% RH there is nonhydrated Ca3SiO5 (tricalcium silicate), Ca2SiO4 (dicalcium silicate) and Ca3Al2O6 (tricalcium aluminate), which means that hydration is not completed. As for carbonatation of the cementitious matrix, we note that the quantity of the ettringite and portlandite is lower on the surface, is more visible in Figures 14 (I (a andc)) and 14 (II (a and c)) respectively. On contrary, Figure 14 (III-a and III-c) shows the increase of calcite. The carbonation of the calcium hydroxide is accompanied by a decrease in pH. That is why the measurements of the pH level were carried out to study the influence of the pH on the crystallographic composition. The results presented in Figure 15 show that the pH levels on the surface and in the middle of the sample were equal to 8.46 and 12.18 respectively before the weathering aging cycles and 8.71 and 8.69 -after. It demonstrates the effect of the carbonatation of the cementitious matrix.

In fact, the diffusion of the air dioxide carbon in the hemp mortar provokes its reaction with the portlandite and decrease the pH level. As the surface of the sample is exposed to the air, the quantity of portlandite at the surface is lover because of the carbonation. Later, that leads to the degradation of the ettringite, as it is very sensible to the ambient pH level.

Hygrothermal properties evolution

Hemp mortar is used for external or internal insulation applications. It is also considered as a hygroscopic material with a good ability to moderate variations in relative humidity of the external environment. The results of their evolution investigation are present below.

First, the water vapor permeability evolution due weathering aging cycles was investigated. This parameter expresses the quantity of water vapor that passes through sample's porosity. Thus, the evolution of this parameter also gives information on the evolution of the porosity and cracking caused by the hygrothermal dimensional variations after aging. The results of water vapor permeability before and after the weathering aging cycles are present in Figure 16. Firstly, one can notice that the results obtained from this study are much higher than the values measured for other formulation of the hemp mortar in the literature [START_REF] Bennai | Experimental characterization of thermal and hygric properties of hemp concrete with consideration of the material age evolution[END_REF]. The measured value at the age of 30 days for wet cup method was 2.02 E -10 kg.s -1 .m -1 .Pa -1 . That is explained by the composition, the particle size of hemp particles, and different binder types that allow to produce more aerated hemp mortar in case of studied formulation. Also, we notice the increase of the average water vapor permeability after aging from 2.02E -10 to 3.41E -10 kg.s -1 .m -1 .Pa -1 , that justify the results of the increase of the porosity. This is explained also by the cracking of the hemp mortar, which increases its vapor diffusion capacity and thus its permeability. Concerning the moisture buffer value, one can note that the MBV of hemp mortar is bigger than the MBV of hemp shives before aging cycles. It is explained by the fact that the studied hemp mortar is much more aerated and lightweight than conventional hemp mortars that allows it to preserve the hygric properties of hemp shives. Nevertheless, we notice that the accelerated aging cycles did not affect the hemp shives moisture buffer value. It stays excellent with 2.25 moisture buffer value according to Nordtest protocol (Figure 11) [START_REF]BYG DTU R-126 Moisture buffering of building materials[END_REF]. Unlike hemp shives, hemp mortar represents the degradation of its MBV value. It decreased from 2.22 g/(%RH.m 2 ) to 2 g/(%RH.m 2 ) (Figure 19). In fact, there is the competition between the interface damage of the hemp mortar samples and the carbonatation of the mineral matrix during the hygrothermal solicitations caused by the accelerated aging cycles. So, the carbonatation of the matrix caused the closing of the microporosity and complicated the air access to the hemp particles that finally led to the moisture buffer value decrease.

Despite the decrease of the moisture buffer value, we note that the MBV of the hemp mortar stays excellent.

Conclusion

In this work, an original accelerated weathering aging protocol inspired by the standardized one used for External Thermal Insulation Composite Systems (ETICS) with rendering was proposed. The particularity of this protocol consists in the representativeness of the real cycles of ageing applied to the building walls without damaging the material by total immersion in water. The proposed cycles of relative humidity and temperature have been studied and normalized to represent the climatic conditions of these insulation materials with rendering.

Experimentally, the hemp mortar microstructure, its composition and properties evolution were investigated. To do this, all characteristics were studied before and after the hemp mortar exposure to accelerated aging cycles. Firstly, the hemp mortar total water porosity, its water vapor permeability, its mineralogical composition and pH level evolutions were analyzed. Also, the SEM observation of the hemp particles / cementitious matrix interface was realized. Then, the characterization of evolution of the hygrothermal properties as thermal conductivity and moisture buffer value was conducted.

Results of the total porosity evaluation show a slight increase of the porosity but it remains in a standard deviation range. This is explained by a compensation of two processes that occurs and are in contest during the hygrothermal stresses: the cracking of the hemp particles / cementitious matrix interface that increase the hemp mortar porosity and the carbonatation that provoke the filing of the porosity by calcite. To study the cracking effect the evolution of water vapor permeability was analyzed. It increased from 2.02E -10 to 3.41E -10 kg.s -1 .m -1 .Pa -1 that shows the increase of the ability to transfer a water vapor. This change was also noticed using a SEM observation because of the increase of the interface cracks on the surface of the mortar. The XRD analysis demonstrated the carbonatation of the mineral matrix -the increase of the quantity of the calcite and decrease of the quantity of the portlandite and ettringite due to the aging cycles. The carbonatation was also demonstrated by the study of pH level evolution since the pH level of the samples decreased to 8.69 in the middle and to 8.71 on the surface after the weathering cycles. As for the hygrothermal properties, the value of thermal conductivity dependent on the porosity follows the same as total porosity pattern -it slightly increases at the range of the standard deviation. The moisture buffer value decreases from 2.22 g/(%RH.m 2 ) to 2 g/(%RH.m 2 ) that is explained by the closing of the hemp mortar microporosity and complication of the air access to the hemp particles. Obviously, both the interface cracking and the carbonatation processes present the considerable influence on the durability of the hemp insulating mortar. At the same time, we note that the studied hemp mortar saves their previously declared hygrothermal properties (as the thermal conductivity and the moisture buffer value). It is primordial when talking about the durability of the insulating material. As recommendation for using the proposed hemp mortar, it must be used according to the construction rules. The control of the continuity and the regularity of coating surface is a critical point because the access to water could lead to the molds proliferation.

This study provides literature with a less aggressive and more representative accelerated aging protocol with an improvement in knowledge when talking about the degradation and the hygrothermal behavior of the biosourced materials considered as very heterogeneous. These results represent a database for the durability of hemp mortars. Efforts are still needed to validate the accelerated aging tests with a comparison to an in situ experimental campaign.

Figure 1 .

 1 Figure 1. Hemp shives Figure 2. Granulometric curve of hemp mortar by image analysis

Figure 4 .Figure 5 .

 45 Figure 4. Diagram of the total accelerated aging cycles (A) with detailed heat-rain (B) and heat-cold (C) cycles

Figure 6 .

 6 Figure 6. Experimental protocol schema

Figure 7 .

 7 Figure 7. Hemp mortar samples without (A) and with (B) protective cotton tissue

Figure 9 .

 9 Figure 9. Hemp mortar sample sealed on the experimental cup

Figure 10 .

 10 Figure 10. CT-Metre measurement device

Figure 11 .

 11 Figure 11. Moisture buffer value [g/(%RH.m 2 )] standardization according to Nordtest protocol[START_REF]BYG DTU R-126 Moisture buffering of building materials[END_REF]

  hemp mortar and mineral matrix interface that show the nature of the hemp mortar damage. That explains the results of the water vapor permeability evolution. Moreover, the mean results of the crack size show the important increase of the crack size from 5.5±1.5 µm to 24.8±3.1 µm in the longitudinal direction of fibers. These results emphasize the influence of hygrothermal solicitations on the morphology of the hemp mortar.

Figure 12 .

 12 Figure 12. W-SEM observation of the surface of the specimens before (a) and after (b) the weathering aging cycles

Figure 13 .

 13 Figure 13. X-ray diffraction patterns of: (a) interior part of the hemp mortar samples before the weathering aging cycles; (b) interior part of the hemp mortar samples after the weathering aging cycles; (c) surface part of the hemp mortar samples before the weathering aging cycles; (d) surface part of the hemp mortar samples after the weathering aging cycles.

Figure 14 .

 14 Figure 14. Zoomed X-ray diffraction patterns for ettringite (I), portlandite (II) and calcite (III) evolution: (a) interior part of the hemp mortar samples before the weathering aging cycles; (b) interior part of the hemp mortar samples after weathering aging cycles; (c) surface part of the hemp mortar samples before the weathering aging cycles; (d) surface portion of the hemp mortar samples after the weathering aging cycles. The same pattern was noted for the sample after the weathering aging cycles. The results show that the ettringite and portlandite almost vanished in the middle and on the surface of the sample (figure 14 (I) and (II)), and calcite quantity increased (Figure 14 (III)). This shows that after the aging protocol the matrix is almost carbonated.

Figure 15 .

 15 Figure 15. Diagram of pH level evolution on the surface and in the middle of the samples before and after the weathering aging cycles

Figure 16 .Figure 17 .

 1617 Figure 16. Water vapor permeability (kg.s -1 .m -1 .Pa -1 ) evolution of hemp mortar samples before and after the weathering aging cycles

Figure 17 .

 17 Figure 17. Thermal conductivity [W/(m•K)] evolution of hemp mortar samples before and after the accelerated aging cycles

Figure 18

 18 Figure18shows the mass and relative humidity evolution that occurred during the first 5 MBV experimental cycles for the second sample of hemp mortar before the weathering aging cycles. The temperature was held constant at 23 °C.

Figure 18 .

 18 Figure 18. Masse (g) and relative humidity (%RH) evolution in function of the time

Figure 19 .

 19 Figure 19. Moisture buffer value g/(%RH.m 2 ) evolution of hemp shives and hemp mortar samples before and after the accelerated aging cycles

Table 1 .

 1 Oxide analysis of used lime-based binder with hydraulic additions

	Oxide	CaO	Al2O3	Na2O	K2O	Fe2O3	SiO2	MgO	SO3
	Percent
	by mass	57.40	2.20	0.15	0.27	0.92	25.40	0.50	3.09
	(wt%)

Table 2 .

 2 Mass ratios of hemp mortar composition

	Component	Hemp shives	Mineral binder	Water
	Percent by mass (wt%)	12.5 %	25 %	62.5 %
		Figure 3. Hemp mortar samples laboratory drying

Table 3 .

 3 Total water porosity values (%) evolution of hemp mortar samples before and after

	the weathering aging cycles

Acknowledgments

The authors gratefully acknowledge ParexGroup SA for financing this project.


D Soulat

A Cheruet

P Boisse

S Guinard

S Maison-Le Poec

Experimental and numerical analysis of interply porosities in composites thermoforming

Keywords: Thermoplastic composites, re-consolidation, shell finite elements

Effective forming of Continuous Fiber Reinforced Thermoplastic Composites requires a detailed understanding and modeling of the forming mechanisms and the development of computational techniques for process simulation. A re-compaction stage is performed at the end of the process in order to avoid interply porosities. The through the thickness behaviour of the composite during forming and consolidation is analysed by a shell finite element with pinching degrees of freedom. This element avoids a locking due to pinching by a modification of the constitutive relation.

INTRODUCTION

The CFRTP (composites with continuous fibres and thermoplastic matrix) forming starts from flat laminated plates made of several unidirectional or textile reinforcements which have been consolidated by a thermoplastic matrix in an autoclave. After heating at a temperature higher than melting (figure 1.1), the forming is made with a punch and die process usually using a rubber on the die (figure 1.2). Finally a reconsolidation is obtained by applying a pressure on the punch (figure 1.3). The objective is to avoid all residual porosity at the interface of the plies. This last stage is important and critical because the health requirements for loaded aeronautical parts are severe. If the use of CFRTP is increasing in the manufacturing of new aircrafts [START_REF] Maison | Technical developments in thermoplastic composites fuselages[END_REF][2] because the processing cycle is shorter than those of thermoset composites, there is an important need of simulation codes of the forming operations in order to define the forming conditions for a given shape and a given stack of plies. The present paper is focused on the experimental and numerical analysis of the compaction stage. Micrographs made on specimens formed at different moment of the process show the evolution of the porosity during forming. They show the porosities are large and consequently the importance of the compaction stage in order to guarantee compaction. 6). After one minute heating, many voids and porosities can be observed in the part (Figure 7).

T fusion

shell element is defined in which a through the thickness strain degree of freedom is introduced.

The stress through the thickness is taken into account and the element can model the compression during the re-compaction stage within a shell finite element analysis. used to stop the process ). There are still voids in flat parts but there is no porosity in curved parts of the reinforcement. For H=0 (complete process), (figure 10) all the voids have disappeared. This last stage corresponds to the end of the re-consolidation phase where a pressure is applied on the tools. That shows that this stage is very important for the final quality of the part. In the curved zones of the part, the reconsolidation is done during the forming accounting for the geometrical effect that leads to transverse stresses during forming.

ANALYSIS OF THE RE-CONSOLISATION STAGE

Experiments

In contrast of classical shell, a sixth degree of freedom appears to describe the displacement in the thickness.. If β =0 equation ( 2) leads to the classical shell kinematics without pinching. We consider β constant thought the thickness. β is an additional pinching degree of freedom. The strain tensor The re-consolidation has been studied in [START_REF] Lee | A model of the manufacturing process of thermoplastic matrix composites[END_REF] and some models have been proposed for the local consolidation. Some studies on consolidation have shown that it depends on the stress state in the laminate and mainly on the normal stress in the consolidation stage. This stress is not present in classical shell theory. Some finite elements with stress/strain through the thickness have been proposed [START_REF] Simo | On a stress resultant geometrically exact shell model. Part 4 : Variable thickness shells with trough-the-thickness stretching[END_REF][5][6] [START_REF] Bischoff | Shear deformable shell elements for large strains and rotations[END_REF]. In the present work a new shell element is used where a through the thickness strain degree of freedom is introduced [START_REF] Coquery | Modélisation d'un joint de culasse multifeuille[END_REF].

T 1 (u) (u) (u) 2  ε=∇+ ∇      θ  β can
The transverse shears are :

m ,1 13 z,1 y m ,2 23 z,2 x β 2ε u θ z β 2ε u θ  +    =+    -     (4)
and the through the thickness strain (pinching) is :

{} 33 ε =
(5) The displacement is obtained in subtracting the initial position from the current position. The transverse shear strains are modified by the pinching. In contrast of classical shell, the normal stress through the thickness is not zero. It is deduced from ε 33 =β using the compaction behavior law. with x z ~ the projection of this point on the midsurface and the position of this point on the normal . It has been shown [START_REF] Coquery | Modélisation d'un joint de culasse multifeuille[END_REF][9] that this element exhibits a "pinching" locking. To avoid this locking it is necessary to modify the constitutive relation in order to remove the coupling between pinching and bending [START_REF] Coquery | Modélisation d'un joint de culasse multifeuille[END_REF] [START_REF] Cheruet | Analyse et simulation de la mise en forme des composites thermoplastiques[END_REF].The example presented figure 12 and 13 (Clamped square plate under pressure) shows the pinching locking with is obtained with a complete behaviour law and the accurate result if pinching and bending are not coupled.

Let

the displacement of the point in the thickness and define β as the strain in the thickness,

0 t z z h - = ∆ 0 z ~h ∆ = β
. In the case of small rotations between two states the displacement expression becomes:

0t t ˆûuz X z X =- × θ + β ( 2
)
This element is used in order to compute the through the thickness stress during the forming stage as well as in the proper compaction phase. As an example the through the thickness stress computed after compaction of a cylindrical shell is shown in figure 14. The radius of the shells is 10 mm, the width is 3 mm and the thickness of each ply is 0.2 mm. In each of the two plies, the thought the thickness is proportional to the projection of the normal on the compaction direction. The simulation of a Z forming and the computation of the through the thickness stress and its comparison with experimental results can be found in [START_REF] Cheruet | Analysis of the interply porosities in thermoplastic composites forming processes[END_REF]. analysis of the forming that porosities induced by the process are important and must be suppressed. This finite element has a degree of freedom in the thickness and give the stress and strain through the thickness. The further work will be focused on the definition and the identification of the compaction constitutive relation.

Figure 6 :

 6 Figure 6 : Composite before heating

Figure 7 :Figure 8 :

 78 Figure 7 : Composite after heating

Figure 9 :

 9 Figure 9 : Composite in curved part for H=5.3mm

Figure 10 :Figure 5 .

 105 Figure 10 : Composite in flat part for H=0

Figure 8

 8 Figure 8 and 9 correspond to H= 5.3 mm (H is the height before the end of forming. Chocks have been

XFigure 12 .

 12 Figure 12. Clamped square plate under pressure

Figure 13 .

 13 Figure 13. Displacement at the center of a plate under pressure.

Figure 14 .

 14 Figure 14. Through the thickness stress in a cylindrical shell compaction