Yufei Han

Gridpocket Xiaolan

Sha Eurecom

Etta Grover

Silva Gridpocket

Pietro Michiardi Eurecom

On the Impact of Socio-economic Factors on Power Load Forecasting

In this work, we study the importance of socioeconomic factors of residential customers for estimating daily peak and total load, at fine temporal granularity. We do so by generating a compact set of heuristic rules using random forests. Compared with black-box time series prediction models in previous works, the rule set we obtain is highly interpretable and makes it easy to fuse human experts domain knowledge. In addition, we quantify and rank the importance of socio-economic features in the rule set for the forecast task. Our experimental analysis, which uses a publicly available dataset of over 3,800 households, providing consumption data for 1.5 year, highlights the superiority of tree-based models over state-of-the-art techniques that use support vector machines.

Introduction

Recent development of smart meter technologies have enabled energy providers to collect fine-grained (e.g., hourly) customer usage records. Analyzing large amounts of consumption data helps energy utilities to build user profiles, which opens opportunities to develop new customer oriented services. The liberalization of energy markets and added transparency of energy data motivates providers to offer energy tariff schemes better suited to customers' personal preferences: several energy providers have already put in place dynamic billing strategies and provide energysaving services based on consumption data analytics [START_REF] Wood | Dynamic energyconsumption indicators for domestic appliances: environment, behaviour and design[END_REF], [START_REF] Mcloughlin | Characterising domestic electricity consumption patterns by dwelling and occupant socioeconomic variables: An Irish case study[END_REF]. From the perspective of grid management systems, collecting and processing large amount of consumption records allow utilities to leverage energy production flexibly, to cope with peak usage, and monitor the health of the distribution grid at fine temporal granularities.

In this context, load forecast is a critical operation that involves the prediction of electricity consumption characteristics, e.g., demanded power level of grid substations or house-holds, over a variety of time horizons (minutely, hourly, daily or even yearly) [START_REF]Energy Vision 2020: Integrated Resource Plan/Environmental Impact Statement[END_REF], [START_REF] Kolter | A large-scale study on predicting and contextualizing building energy usage[END_REF], [START_REF] Centra | Hourly Electricity Load Forecasting: An Empirical Application to the Italian Railways[END_REF], [START_REF] Hesham | Electric load forecasting: literature survey and classification of methods[END_REF], [START_REF] Humeau | Electricity Load Forecasting for Residential Customers: Exploiting Aggregation and Correlation between Households[END_REF]. Utilities greatly benefit from load forecasts for the managament of the energy supplydemand balance and the interaction of flexible load switching schemes [START_REF] Aman | Improving Energy Use Forecast for Campus Micro-grids using Indirect Indicators[END_REF]. In energy science, many studies focus on how to improve forecast accuracy using statistical time series processing algorithms [START_REF]Energy Vision 2020: Integrated Resource Plan/Environmental Impact Statement[END_REF], [START_REF] Hesham | Electric load forecasting: literature survey and classification of methods[END_REF], [START_REF] Aung | Towards Accurate Electricity Load Forecasting in Smart Grids[END_REF], [START_REF] Miranda | Fuzzy Inference Applied to Spatial Load Forecasting[END_REF], [START_REF] Fu | Models for long-term energy forecasting[END_REF], [START_REF] Kermanshahi | Up to year 2020 load forecasting using neural nets[END_REF]. In contrast, little attention has been devoted to studying the impact of socioeconomic factors in building customer usage profiles. Like any other private consumption records, e.g., online shopping histories, electricity usage is also a consumer oriented process, which can be inferred from customers' behavioral habits and householding characteristics [START_REF] Mcloughlin | Characterising domestic electricity consumption patterns by dwelling and occupant socioeconomic variables: An Irish case study[END_REF], [START_REF] Beckel | Towards automatic classification of private households using electricity consumption data[END_REF], [START_REF] Beckel | Automatic socio-economic classification of households using electricity consumption data[END_REF]. As opposed to historical consumption associated with proprietary meters linked to a specific house, socio-economic factors are bound to customers and are expressed independently of a specific provider. The association between such factors and energy consumption patterns allows new incumbents to estimate energy demand profiles, even if few and incomplete historical consumption records are available. Moreover, although smart meter technology is emerging in the landscape of energy-related industry, energy providers still stand at the very early stage of deploying new generation meters. For most end-users, energy providers are unable to collect fine-grained historical consumption records.Therefore socio-economic customer profiles turn out to be a quintessential asset that should be exploited.

The aim of our work is two-fold: i) we propose a heuristic forecast model of customers' daily energy usage patterns based on their socio-economic factors and, ii) we identify the most important factors for forecasting purposes. As such, the results of our work can help utilities to investigate characteristics of energy usage of each end-user, avoiding verbose questionnaires and reducing the intrusiveness of private information collection. Our empirical study uses a large-scale public dataset of electricity usage patterns of residential consumers. We proceed by defining two behavioral indicators, namely peak load and total electric consumption as the forecasting targets. Then, we design a forecasting model based on random forests [START_REF] Eiber | Pruning decision trees and lists[END_REF], [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] to model the association between customers' socio-economic factors and the two indicators. Essentially, our ultimate goal is to estimate energy usage profiles based on their socio-economic attributes. Applications of the methods we study in this work are numerous. Based on the two selected indicators, energy providers can set a proper power limit for each customer and protect the grid from risk of overloading during daily peak hours [START_REF]Energy Vision 2020: Integrated Resource Plan/Environmental Impact Statement[END_REF], [START_REF] Miranda | Fuzzy Inference Applied to Spatial Load Forecasting[END_REF], [START_REF] Fu | Models for long-term energy forecasting[END_REF]. Utilities can use the forecasting results to tune billing prices and energy production dynamically, according to variation of energy usage [START_REF] Kolter | A large-scale study on predicting and contextualizing building energy usage[END_REF], [START_REF] Centra | Hourly Electricity Load Forecasting: An Empirical Application to the Italian Railways[END_REF]. For individual customers, the forecasting model can help them adjust their energy usage and reduce their bills [START_REF] Abreu | Household Electricity Consumption Routines and Tailored Feedback[END_REF].

The remainder of the paper is organized as follows. A review of the related research progress in load forecasting and a description of the necessary background on random forests is given in Section 2. In Section 3, a brief overview the statistical properties of the dataset used is performed along with a quantification of the information conveyed by the socio-economic factors extracted from the dataset. In addition, the feasibility of using data-driven models to estimate energy usage patterns through socio-economic factors is discussed. Section 4 describes the experimental configuration in detail. Moreover,the forecasting performance of this approach is compared to the current state-of-the-art. In doing so,the stability of our approach is investigated and a quantified ranking of the importance of the socio-economic features is defined for the forecast task.Section 5 concludes the paper.

Background and Related Work

In this Section, current efforts toward the design and operation of load forecast mechanisms is examined andthe necessary background on the techniques used in this paper is provided, namely decision trees and random forests.

Current Approaches to Load Forecasting. Load forecasting has received a lot of attention lately: most research work in this domain can be categorized into three categories.

Firstly, time-series based methods are used to model temporal causality of energy demand between the past and future [START_REF] Centra | Hourly Electricity Load Forecasting: An Empirical Application to the Italian Railways[END_REF], [START_REF] Hesham | Electric load forecasting: literature survey and classification of methods[END_REF], [START_REF] Humeau | Electricity Load Forecasting for Residential Customers: Exploiting Aggregation and Correlation between Households[END_REF], [START_REF] Aung | Towards Accurate Electricity Load Forecasting in Smart Grids[END_REF], [START_REF] Kermanshahi | Up to year 2020 load forecasting using neural nets[END_REF], [START_REF] Taylor | Neural Network Load Forecasting With Weather Ensemble Prediction[END_REF], [START_REF] Mohandes | Support vector machines for shortterm electrical load forecasting[END_REF], [START_REF] El Desouky | Hybrid adaptive techniques for electric-load forecast using ANN and ARIMA[END_REF]. Popular machine learning algorithms, such as ARIMA [START_REF] El Desouky | Hybrid adaptive techniques for electric-load forecast using ANN and ARIMA[END_REF], SVM [START_REF] Centra | Hourly Electricity Load Forecasting: An Empirical Application to the Italian Railways[END_REF], [START_REF] Mohandes | Support vector machines for shortterm electrical load forecasting[END_REF] and neural networks [START_REF] Kermanshahi | Up to year 2020 load forecasting using neural nets[END_REF], [START_REF] Taylor | Neural Network Load Forecasting With Weather Ensemble Prediction[END_REF], build black-box models of temporal dynamics. Thus, prediction is achieved following the learned functional mapping between past and future energy usage patterns. Accurate as they are, they offer no insights on the underlying phenomenon behind the time evolution of energy demand, and the predicitive models cannot be easily generalized.End user models are commonly userd as an alternative to black-box methods [START_REF] Wood | Dynamic energyconsumption indicators for domestic appliances: environment, behaviour and design[END_REF], [START_REF]Energy Vision 2020: Integrated Resource Plan/Environmental Impact Statement[END_REF], [START_REF] Kolter | A large-scale study on predicting and contextualizing building energy usage[END_REF], [START_REF] Hesham | Electric load forecasting: literature survey and classification of methods[END_REF], [START_REF] Aman | Improving Energy Use Forecast for Campus Micro-grids using Indirect Indicators[END_REF], [START_REF] Fu | Models for long-term energy forecasting[END_REF], [START_REF] Abreu | Household Electricity Consumption Routines and Tailored Feedback[END_REF].This model requires information about housing conditions, electrical appliance usage and environmental factors. The key idea of the this model is to disaggregate daily energy consumption into elementary components -including heating/cooling, water usage, cooking and other behaviors -which are used to interpret the temporal variation of clients' energy demand. Thus, forecasting models use relations between the collected information and energy consumption profiles. The shortcoming of enduse models is that forecasting performance depends heavily on the quality of available information, which makes them sensitive to noise.

Econometric methods [START_REF]Energy Vision 2020: Integrated Resource Plan/Environmental Impact Statement[END_REF], [START_REF] Kolter | A large-scale study on predicting and contextualizing building energy usage[END_REF], [START_REF] Hesham | Electric load forecasting: literature survey and classification of methods[END_REF], [START_REF] Beckel | Towards automatic classification of private households using electricity consumption data[END_REF], [START_REF] Beckel | Automatic socio-economic classification of households using electricity consumption data[END_REF], [START_REF] Abreu | Household Electricity Consumption Routines and Tailored Feedback[END_REF] combine the above two techniques. Such models estimate the relationship between energy consumption profiles and the factors influencing consumption behavior, using the least fitting error criterion as in timeseries based methods. Econometric models are built by learning the mapping from pairs of input factors and output energy consumption profiles automatically, which is appealing for realistic application deployments. Hence, this category has gained popularity recently.

Although energy demand-supply does not closely obey typical market laws, the purchase-consumption relation between utilities and end-users indicates energy usage is intrinsically a type of consumption behavior similar to on-line purchasing in e-commerce. The basic motivations driving clients to purchase more electricity is to satisfy their living requirements. Therefore, clients' comfort preferences, income levels, family structure, residential status, electrical appliances and environmental features, or, in short, their socioeconomic status, can therefore change electricity usage behavior dramatically, which in turn determines consumption profiles.

As discussed above, end-use model explicitly integrate such factors in the forecasting model. However, they require human experts' interference to tail the input factor set. To the best of our knowledge, the integration and ranking of socio-economic features for load forecasting through automatic learning from energy usage data has remained elusive in the literature. Using the methodology described in this work, we study how to appropriately, and automatically select socio-economic factors to build an energy forecast model. In doing so, we unveil which are the most relevant features that contribute to accurate predictions, and which can be safely dismissed as they contribute little to the predictive power of our model.

Random Forest Models.

The forecasting model in this work is based on random forests. Its concept dates back from the indenpendent work of Tin Kam Ho and Amit and Geman [START_REF] Eiber | Pruning decision trees and lists[END_REF], [START_REF] Breiman | Random Forests[END_REF] as a an ensemble variant of decision tree [START_REF] Eiber | Pruning decision trees and lists[END_REF], [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. The central idea is to construct a set of decision trees independently by bootstrapping training data. The output of a random forest is obtained, e.g., through majority voting or with a simple average of the output of each individual tree model. Each decision tree in a random forest is a tree-like rule chain based on input variables for classification or regression. Each rule is a branch-split operation comparing one input variable with a predefined threshold. An input sample is then forwarded to different decision branches for finer analysis. The hierarchical split-branch operations form a coarse-to-fine white box model, which can explain explicitly how input factors are combined to achieve a complete decision making procedure.

The main motivation of random forests stems from the fact that a single tree model is sensitive to noise and produces predictions with large variance. Random forests mitigate the issue by injecting randomness into the tree structures and finally aggregating the output of the decision trees built independently. First, bootstrap sampling helps produce different randomly sampled training data subsets. Then, a decision tree is trained using each subset. In addition, input attributes are also bootstrapped, and the selected attribute subsets are used to learn the tree structure: as such, individual decision trees are largely independent one from each other. Averaging the outputs over several learned trees can substantially reduce the variance and improve stability in the final decision output. Despite a seemingly complex procedure, random forests can be used to assess the importance for each input attribute.

The properties of random forest makes itself a perfect fit for the goal of our work. One further aspect that has not been addressed so far, pertain to the scalable implementation of a random forest algorithm. Although out of the scope of this paperand thus excluded from our contributions -significant development was dedicated to parallel design of a random forest algorithm. We did so by relying on the functional paradigm offered by data-intensive scalable computing systems such as Apache Spark [START_REF]Apache Spark[END_REF]. Essentially, the proposed algorithm follows the MapReduce programming model, and the forecasting models are built by executing on a parallel platform, particular attention was given to the iterative nature of the random forest process: intermediate data materialization was avoided (which contributed to a large extent to poor runtime performance) using in-memory, reliable data structures. We omit further details of our algorithm to avoid distracting readers from the goal of this work.

Load Forecasting using Socioeconomic Factors

A detailed description of the dataset used in this work is provided, which identifies the features used to learn the proposed model. The statistical properties of such features are analyzed with emphasis on socioeconomic indicators. As a general remark, the large volume of the available training samples offers a large coverage of energy and socio-economic factors, which in turns improve model generability.

3.1.The dataset

The CER ISSDA dataset is a publicly available energy consumption trace, collected by the Irish Commission for Energy Regulation (CER) in a smart meter study [START_REF]Eletricity Customer Behaviour Trial[END_REF]: it contains electricity consumption data of 4,225 private households and 485 small / medium enterprises; the trace covers 1.5 years (from July 2009 to December 2010). For each customer, the daily load curve is sampled every 30 minutes: energy data can be thought of as a series of timestamps and energy readings. In addition to energy data, the dataset includes a series of survey sheets and answers for each consumer, describing their housing condition, occupancy, employment status, income level, social class, appliance usage information and other socioeconomic factors. 41 survey questions belonging to five categories were carefully selected, that allowed to build consumer profiles based on heating and lighting behavior, hot water and other electrical appliances usage: a complete list of features is shown in Tab.

3.2.Predictors used in the forecasting task

In addition to the socio-economic factors outlined in Sec. 3.1, we complement the predictors we use in the forecasting task with environmental features, which are significant for energy consumption profiling [START_REF] Taylor | Neural Network Load Forecasting With Weather Ensemble Prediction[END_REF]. The first two are heating degree days and cooling degree days, evaluating quantitatively the needs to start heating or air conditioning appliances to keep an adequate environmental temperature. We also consider the daily average humidity, which represents the average air humidity level during a one-day interval: indeed, days with the same temperature but higher humidity level generally needs more energy to keep warm in winter or cool in summer. Finally, we include month of year, day of week and holiday index: the first feature represents the seasonal weather change, directly and strongly affecting energy usage profiles of all involved customers; the others features differentiate occupancy patterns of residential customers.

The input of our data-driven forecasting model includes a total of 41 factors, that are summarized in Tab. 1. Note that our predictors are of mixed formats, including both numeric and categorical types. As a consequence, a random forest model -which treats predictors uniformly in the splitting operation -is a more appropriate choice as compared to previously used statistical regressors, such as Support Vector Machine (SVM) [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF], for which categorical predictors are problematic.

3.3.Statistical profile of the load forecast task

The statistical profile of the dataset can be constructed using well-known feature ranking statistics to provide formal grounding to the key idea of this work: the goal is to show that socio-economic factors play an important role as predictors for a forecasting model.

As we observe in Fig. 1, the distributions of daily peak loads and daily total consumption are light-tailed: the majority of the measurements locate around the mean, and there are no outliers. 1 mutual information criterion [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF], [START_REF] Witten | Data Mining: Practical Machine Learning Tools and Techniques[END_REF] was then used to illustrate the correlation between each input factor and the forecasting target: the larger the mutual information, the higher the correlation between the two random variables. Specifically, Kullbeck-Leiber divergence of the product of marginal distributions of two random variables x and y from the joint distribution of them was used, as illustrated in Eq. 1:

I(X; Y ) = y∈Y x∈X p(x, y) log p(x, y) p(x) p(y) (1)
Fig. 2a and 2b report the mutual information for all the 41 socio-economic factors extracted from our dataset. More specifically, Tab. 2 reports the top-15 factors (sorted in descending order of mutual information), that are mostly related to daily peak loads and daily total energy usage respectively. We emphasize that here, the mutual information is not used to rank predictors for the purpose of feature selection. From Tab. 2, occupancy attributes, the number of bedrooms, age of the clients, employment status of the clients, usage of tumble dryer and usage of dishwashers, are the common key elements influencing both peak loads and total amount of energy consumption, which are related to clients consuming capacities, potential electricity needs to keep proper temperatures and daily housework.

1. In general, outliers are often caused by either mistakes in measurement or faulty meters, and they usually increase the bias of data-driven models, leading to performance degradation. For the purpose of the forecasting task, the analysis indicates that cooking appliances -because they consume a lot of energy -largely contribute to the magnitude of the daily peak energy consumption. Instead, such appliances only marginally affect the the daily total energy consumption, as their usage is restricted to a small interval of time. In contrast, our analysis indicates that house category and usage of a standalone freezer device are key factors for the daily total energy consumption, but are weakly related to the daily peak load. Such results follow common sense: the type of the house is related to the thermal performance and inertia of the house, which crucially affects the total amount of energy -required to maintain a desirable temperature in the house. It has been shown that a datadriven model based on socio-economic factors may indeed find associations between user profiles and their consumption behavior. This method will be compared with the state of the art for validation.

Experimental Evaluation

The performance of the proposed random forest model is evaluated for forecasting daily peak loads and total electricity consumption of each residential user. In addition, the proposed method is compared to a state-of-the-art approach of time series forecasting, Support Vector Machine (SVM), for validation purposes. We conclude with an analysis of the importance ranking of input socio-economic and environmental factors.

4.1.Comparison of Forecasting Models

The experimental methodology is as follows. For each of the 3,822 users and for the 1.5 year duration of our dataset, daily peak load measurement (in KiloWatt) and daily total electricity consumption (in KiloWattHour) are extracted to build the prediction targets. Input predictors are the 41 features described in Sec. 3. Therefore, we obtain 1,615,541 training samples, containing the input-output pairs for each of the forecasting tasks.

Metrics. For all models involved in our comparative analysis, we measure the forecasting performanceprecisely, the accuracy of model fiiting -using coefficient of determination R 2 [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF], [START_REF] Witten | Data Mining: Practical Machine Learning Tools and Techniques[END_REF]: it is defined by the ratio between the sum of square regression residual error and total sum of squares of the forecasting target. Large values of R 2 indicate an agreement between the model and the underlying real output, which translates in higher accuracy. Note that for really biased predictions, the coefficient can be zero or negative, while it is upper bounded by 1.

In this work, the relative ratio coefficient is used instead of the sum of square error, such that the variance of underlying forecasting target is taken into consideration to define an accuracy indicator. Higher variance implies that it is more difficult to construct a model to explain the overall data distribution. Thus, the R 2 coefficient is a more fair evaluation criterion of the model forecasting power than the sum of square errors. Eq.2 formalizes the performance metric, where y i and f i are the ground truth and the corresponding predicted value and ȳ is the empirical expectation of the ground truth.

R 2 = 1 -i (f i -y i ) 2 i (y i -ȳ) (2)
Forecasting models. The main parameter of random forests is the number of trees constructed in the ensemble model to perform a final vote. In the experiment, this parameter is varied in the range {50, 100, 200, 300, 400}, and the stability of forecasting performance is investigated. Once the random forest model is constructed, the out-of-bag error [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF], [START_REF] Witten | Data Mining: Practical Machine Learning Tools and Techniques[END_REF] Model of the random forest is used directly as the estimation of the model generalization error. Note that the splitbranch operation in the random forest constructs a piecewise linear regression model, which can handle non-linearly distributed data without introducing additional algorithmic components.

Since the training set consists in more than 1.8 million records, a standard approach such as SVM can not afford the construction of a huge kernel matrix to enforce a non-linear regression in memory: hence, we use a linear kernel SVM. The training configuration of SVM is carefully selected through a 5-fold crossvalidation. After fixing the configuration parameter, the SVM training ends by providing a linear regressor: 5-fold cross-validation is employed again to estimate its generalization error. Overall, the computational complexity of building the random forest model is O(mn log(n)), where m is the number of trees, n is the number of bootstrapped training samples to construct each tree in the model. Considering that m is usually much smaller than n, the complexity of random forest is comparable to or smaller than the quadratic programming cost of SVM [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF], [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF], which is between O(n 2 ) and O(n 3 ). It is noted that the implementation of the random forest uses the MapReduce programming model, which is compiled to be executed on the Apache Spark platform. In the interest of space, the algorithmic details of the approach is not provided. However its capability goes beyond the embarassingly parallel nature of the random forest model, by also using parallel algorithms to build individual trees. As such, the method can scale to large datasets.

Results. Tab. 3 summarizes the generalization error for the propsed method and SVM. In general, this method outperforms SVM: precisely, the highest accuracy is achieved with an ensemble of 100 trees, for both forecasting tasks. Indeed, although more trees can reduce variance, they might increase the model bias. There are two main reasons for the superiority of the approach. First, SVM treats categorical inputs as numerical variables, which hinders the task of defining an appropriate similarity measure: therefore, categorical information is essentially lost during the construction of the SVM model. Furthermore, due to the computational complexity and memory requirements, the SVM model is linear: thus, it cannot model the non-linear relation between socio-economic features and energy consumption profiles. In contrast, the random forest method does not suffer from such limitations. For illustrative purposes, Fig. 3 and4 show two instances of the forecasting task for daily peak load and daily total consumption respectively: the solid line corresponds to the ground truth of daily peak load or daily total consumption measurements, while the dashed line represents the estimated values using the random forest model. These figures reveal that our method achieves higher predictive power for the total daily consumption than for the daily peak load, which confirms our findings in Tab. 3. The random forest model obtains distinctively higher R 2 in daily total consumption prediction. The results are explained as follows. Peak loads, as an instant power level measurement, are easily affected by householders' occupancy duration, appliances' working status, and in general, time-of-day related behavior. Such effects inject randomness into daily peak load patterns, which in turn are hardly represented comprehensively by the questionnaire used to derive input features of our model. In contrast, the daily total consumption, as a cumulated sum of energy consumption measurements, is less prone to random fluctuations due to time-ofday behavior. As a consequence, the model built with the approach contains a strong association to the static socio-economic features available in the dataset.

4.2.Importance Ranking of the Socio-economic and Environmental Factors

An important by-product of the proposed random forest model is the inherent ranking of each input factor toward the forecasting task. The quantitative scores produced by our method are used to evaluate the importance of association between the socio-economic factors and daily energy usage profiles.

The ordering rank of input factors were recorded, which were derived through bootstrap sampling the training data during the construction of the random for-est model. Therefore, the mean and standard deviation of the feature importance score of each input factor can be estimated. Note the standard deviations are at least three times smaller than the mean values for all top 15 factors. It indicates the ranking result is statistically stable with respect to the random forest parameter, i.e., the number of trees used to build the model and training data sampling bias. Hence the factor importance ranking is beyond only data-driven output, but represents underlying association rules contained in the data. Thus the model built with 400 trees was used as a reference in the following discussion.

The experimental results are summarized in Tab. 4, which reports the top-15 input factors for predicting peak loads and total energy consumption. The ranking obtained is consistent to that obtained using the mutual information criterion in Sec. 3. The number of occupants is ranked as the most important feature, affecting both daily peak load level and the total amount of energy consumption, which follows the intuition of a direct relation between household size and energy consumption. General appliances and age are common features for the two forecasting tasks, but they are ranked higher for the daily peak load than for total consumption prediction task. Indeed, appliances contribute to a large extent to instantaneous peaks in energy consumption, and are generally responsible for aggregate energy consumption as well. Appliances that have a daily cycle (e.g. freezers) and average daily heating degree days are also shared features in two forecasting tasks. However, they are ranked much higher for predicting the total daily energy consumption, as they contribute to a steady energy consumption, rather than representative of peak demand. Some important distinctive traits between the ranking obtained by the proposed method and that based on mutual information are examined. For peak load forecasting, the random forest method boosts the importance of the water heating appliances, while it downgrades indirect factors, such as the number of bedrooms, and completely neglect some economic factors, like social class of the chief income and employment status. For total daily energy consumption, the random forests method assigns higher importance scores to factors related to heating usage -including environmental heating degree index, month of year, the year when the house is built, the number of bedrooms -while it decreases the importance of social class, age and employment status of the owner. Overall, the ranking produced by the random forest method is in line to that obtained through a manual inspection of the data by domain experts. It should be noted that environmental temperature factor does not play an important role in both forecasting tasks. The possible reason is high heating loads at all times due to Irish weather status, which limits the capability to trend the consumption data with varying temperatures. Furthermore, the collected temperature data is lack of precision for specific region, coupled with lack of exact knowledge of the geographic location of each household.

Conclusion

In this paper, the issue of forecasting daily electricity usage patterns and evaluating the influence of clients' socio-economic factors on energy consumption has been addressed.

A forecasting method based on random forest that seamlessly incorporates both clients' socio-economic factors and environmental factors has been proposed, resulting in an ensemble of split-branch decision rule chains, whereby a voting mechanism is used to achieve a stable prediction. The rule-chain based structure enabled the explaination of how each input features contributes to energy consumption forecasting, thus unveiling the underlying physical association between predictors and the daily energy usage pattern. These learned associations can either be used to discover unknown factors of energy consumption behaviors, or they can be applied as a complementary decision support to human experts.

Experimental results based on large-scale energy consumption records showed that the proposed random forest method dramatically outperforms a state-of-theart SVM-based method, both in terms of prediction accuracy and in terms of scalability. Additionally, our work evaluated automatically importance of socioeconomic factors in the forecasting tasks. The result can guide the construction of socio-economic surveys without requiring human intervention, leading to a decreased intrusiveness of measurement campaigns.

As shown in these experiments, daily peak load forecasting is a difficult task, compared to aggregate energy consumption. Future work intends to build detailed usage patterns of large power appliances that play a key role in determining peak loads and integrate historical peak load information to extend the input set of the forecasting model. In addition, we will provide details on our parallel algorithm design of the random forest model, and focus on its scalability properties by conducting a large-scale experimental campaign.

Figure 1 :

 1 Figure 1: Probability density function of the two forecasting targets.

Figure 2 :

 2 Figure 2: Mutual information of each feature in the prediction of peak load consumption and daily consumption.

Figure 3 :

 3 Figure 3: Prediction of peak load level of one user during one month.

Figure

  Figure Prediction of total consumption usage of one user during one month.

Table 1 :

 1 1.For this work, users with less than 10% of survey Breakdown of features extracted from surveys in the dataset and additional environmental factors.

	Category	Features
		Total number of occupants
	Household Profiles	Number of occupants > 15 years old
		Number of occupants < 15 years old
		Gender
		Age
	Chief Income Owner	Employment status
		Social class of the chief income
		Education level
		Interests in reducing bills by installing smart
		meters
		Interests in helping protect environments
	Behavior Related	Interests in helping others who live with you
		to reduce bills
		Having succeeded in reducing the energy cost
		Having made changes in the energy consump-
		tion
		Willing to reduce energy consumption to pro-
		tect the environment
		Having the knowledge to reduce the energy cost
		Type of the house
	House Characteristics	The year when the house is built Number of bedrooms
		The proportion of double gazed windows in the
		house
		Type of the heating appliances in the house
		Type of the appliances to heat water in the
	Main Electrical Appliances	house Type of the cooking appliances
		Type of the heating energy control system
		Type of the heating water control system
		Whether the house is always kept with adequate
		temperature
		Whether there is a washing machine
		Whether there is a tumble dryer
		Whether there is a dishwasher
		Whether there is an instant electric shower
	Other Electrical Appliances	facility Whether there is an electric shower which
		pumps from hot water tank
		Whether there is an electric cooker
		Whether there is an electric heater
		Whether there is a stand alone freezer
		Whether there is a water pump or electric well
		pump or pressurized water system
		Whether there is an immersion
		Heating Degree Days
	Weather	Cooling Degree Days
		Daily Average Humidity
		Is Holiday
	Temporal Features	Month of the Year
		Day of the Week
	coverage are eliminated, which results in a dataset of
	3,822 households.

Table 2 :

 2 Top 15 factors in prediction tasks using mutual information criterion

	Prediction of daily peak loads	Prediction of daily total en-
		ergy usage
	Total number of occupants	Total number of occupants
	Number of occupants > 15	Number of occupants > 15
	years old	years old
	Whether there is a dishwasher	Whether there is a dishwasher
	Whether there is a tumble dryer	Number of bedrooms
	Number of occupants < 15	Whether there is a tumble dryer
	years old
	Whether there is an electric	Number of occupants < 15
	cooker	years old
	Number of bedrooms	Whether there is a stand-by
		freezer
	Type of the cooking appliances	Age
	Age	Type of the house
	Employment status	Employment status
	Whether there is an instant elec-	Whether there is an electric
	tric shower facility	cooker
	Willing to reduce energy con-	Month of year
	sumption to protect the environ-
	ment
	Social class of the chief income	Type of the heating appliances
		in the house
	Type of the heating appliances	Type of the appliances to heat
	in the house	water in the house
	Type of the appliances to heat	Social class of the chief income
	water in the house

Table 3 :

 3 Generalization error of peak load and total electricity consumption prediction, with the Random Forest Model (RF) and SVM. Results for the RF model include values for the tree parameter (in parenthesis) we used.

		Peak Load: R 2	Total Load: R 2
	RF(50)	0.5214	0.7151
	RF(100)	0.5240	0.7261
	RF(200)	0.5127
	RF(300)	0.5208	0.7100
	RF(400)	0.5157	0.7070
	SVM	0.4816	0.6572

Table 4 :

 4 Top 15 factors in prediction tasks.

	Prediction of daily peak loads	Prediction of daily total en-
		ergy usage
	Total number of occupants	Total number of occupants
	Number of occupants > 15	Number of occupants > 15
	years old	years old
	Type of the cooking appliances	Average daily heating degree
		days
	Whether there is an electric	Number of bedrooms
	cooker
	Average daily heating degree	Whether there is a dishwasher
	days
	Whether there is a tumble dryer	Month of year
	Whether there is a dishwasher	Whether there is a tumble dryer
	Type of the appliances to heat	Whether there is an electric
	water in the house	heater
	Whether there is an instant elec-	Whether there is a stand-by
	tric shower facility	freezer
	Whether there is an electric	Type of the appliances to heat
	heater	water in the house
	Age	Whether there is an electric
		cooker
	Whether there is a stand-by	The year when the house is built
	freezer
	Education level	Whether there is an electric
		shower which pumps from hot
		water tank
	Having made changes in the en-	Employment status
	ergy consumption
	Number of bedrooms	Age


Stéphane Etienne

M J Pettigrew

C Charreton

C B Éguin

A Ross

S Étienne

D Épartement De

Two-phase damping for internal flow : physical mechanism and effect of excitation parameters

Keywords: Two-phase, damping, flow pattern, bubble motion diphasique, amortissement, régime d'écoulement, mouvement de bulle

Two-phase damping for internal flow: Physical mechanism and effect of excitation parameters

C Charreton, Cédric Béguin, A Ross,

Introduction

Two-phase flow induced vibration in steam generators is well documented in the literature ( [START_REF] Pettigrew | Fluidelastic instability of heat exchanger tube bundles; review and design recommendations[END_REF], [START_REF] Weaver | A review of cross-flow induced vibrations in heat exchanger tube arrays[END_REF], [START_REF] Au-Yang | Flowinduced vibrations in power and process plant components -progress and prospects[END_REF]). Extensive experimentations have been carried out over the last forty years to get a better hand on the several excitation mechanisms involved, such as quasi-periodic forces or fluidelastic instability. Review and design guidelines for heat exchangers constructors have been proposed, improving nuclear power plant safety and reliability. The duality of two-phase flows from a vibration point-of-view lies in the fact that they bring about destructive phenomena while causing significant damping on the structure. Damping is a crucial input parameter to predict vibration effects in steam generators. However, the nature of the damping is not well understood. A better knowledge of the physical mechanism involved would lead to improved modeling of vibration effects in the near future.

The first experimental studies on two-phase damping were performed by [START_REF] Carlucci | Damping and hydrodynamic mass of a cylinder in simulated two-phase flow[END_REF] and [START_REF] Carlucci | Experimental studies of damping and hydrodynamic mass of a cylinder in confined two-phase flow[END_REF]. For a cylinder confined in axial two-phase flow, they found that total damping is strongly dependent on void fraction. Moreover, the two-phase damping component is much higher than the damping due to fluid viscosity for single-phase flow. It can reach up to 3%. [START_REF] Hara | Analytical model for evaluating added mass and damping of a vibrating circular rod in two-phase fluid[END_REF] derived an analytical model for a cylinder confined in axial-two-phase flow. They modeled the gas phase as columns having no mass nor stiffness. Cylinder and gas motions were described by beam equations, coupled by the fluid forces. Coupling coefficients were extracted from potential flow theory. The eigenvalue problem was solved to find the damping coefficients which compared well with the experiments. However, this approach does not provide a physical explanation for the mechanism. More recently, [START_REF] Uchiyama | Numerical prediction of added mass and damping for a cylinder oscillating in confined incompressible gas-liquid two-phase mixture[END_REF] proposed a numerical simulation, also for a cylinder confined in two-phase flow, assuming a bubbly flow. Damping values of the same order of magnitude as in [START_REF] Carlucci | Experimental studies of damping and hydrodynamic mass of a cylinder in confined two-phase flow[END_REF] were observed. However, the damping ratio goes down to zero for void fractions higher than 60%, which is not verified experimentally. This fact raises questions about the applicability of numerical codes for high void fraction. Indeed, these codes assume bubbly flow but usually do not take flow pattern transitions into account. Nevertheless, [START_REF] Uchiyama | Numerical prediction of added mass and damping for a cylinder oscillating in confined incompressible gas-liquid two-phase mixture[END_REF] explains damping by "the phase lag of the drag force acting on the cylinder behind the cylinder displacement". This introduces a notion of relative displacement inherent to the two-phase mixture.

Two-phase damping has also been measured for cross-flow. It is undoubtedly the most important flow configuration since most vibration mechanisms are critical in the U-bend region of the steam generator. Semi-empirical relations for design purposes are given by [START_REF] Pettigrew | Damping of heat exchanger tubes in two-phase flow: review and design guidelines[END_REF]. They compiled a considerable amount of data to identify the most influent parameters. It was shown that flow velocity and tube frequency have minor influence, contrary to confinement and surface tension. Therefore, several phenomena are potential dissipative mechanisms that could be responsible for two-phase damping. Mainly, flow structure, relative motion of gas phase and liquid phase and coalescence/breakup of bubbles are suspected.

Other studies were steered towards the influence of fluid properties on two-phase damping at Ecole Polytechnique of Montreal. These were performed for internal axial flow on clamped-clamped tubes. This configuration is less interesting from a practical point-of-view since in CANDU nuclear plants, pressurized heavy water is supposed to be almost liquid inside steam generator tubes. Still, it is interesting to notice that two-phase damping vs void fraction curves are oddly the same for the three flow configurations reported : annular, internal axial and cross-flow (see [START_REF] Carlucci | Experimental studies of damping and hydrodynamic mass of a cylinder in confined two-phase flow[END_REF], [START_REF] Gravelle | Damping of tubes due to internal twophase flow[END_REF] and [START_REF] Pettigrew | Damping of heat exchanger tubes in two-phase flow: review and design guidelines[END_REF] respectively). This suggests that the mechanism involved is the same for each case. From a design point-of-view, internal axial flow configuration is the simplest. [START_REF] Gravelle | Damping of tubes due to internal twophase flow[END_REF] reported damping measurements in 20 mm tubes. The decrease of two-phase damping at the transition between bubbly and slug flow is explained by the decrease of interface surface area when slugs appear. This is somehow contradictory with Pettigrew's observations: ζ was found to increase with surface tension σ. A possible explanation given was that bigger bubbles are more prompt to dissipate energy. This work was pursued by [START_REF] Béguin | Influence of viscosity, density and surface tension on two-phase damping[END_REF] who tested several air-liquid mixtures, in order to assess the effect of viscosity and density on two-phase damping. [START_REF] Béguin | Two-phase damping and interface surface area in tubes with vertical internal flow[END_REF] also performed two-phase damping experiments with rigid spheres in sedimentation in stagnant liquids. It appeared that density difference between phases has a major effect, contrary to the viscosity. Damping values with rigid spheres were somehow smaller than in air-liquid mixtures by a factor 2, but proportionality with respect to interface surface area was confirmed. For large number of spheres, interaction occurring between spheres (e.g. onset of coalescence in case of a gas phase) seems to modify this trend. [START_REF] Béguin | Two-phase damping and interface surface area in tubes with vertical internal flow[END_REF] also presented a 2D model of a bubble in an oscillating tube filled with liquid, and solved the Navier-Stokes equations analytically. They showed that viscous dissipation due to the presence of a bubble can be related to the relative motion of the bubble with respect to the structure.

These conclusions motivated the design of a new test section which would not only allow damping measurements but also let us observe the gas phase behavior. Also, damping values have only been extracted at the natural frequency of the considered system so far. Thus, the objective of this project is to measure two-phase damping accurately, so as to relate it to the relative motion of the gas phase that we also measured.

In the next section, a new test rig is presented. It allows to command excitation parameters, such as frequency, on a structure subjected to internal two-phase flow. This leads to interesting information on the nature of the two-phase flow energy dissipation. Then, the experimental parameters involved in the system are described. In 4, the technique to extract the two-phase damping component of the oscillating structure is presented. Results on damping and the relation with flow patterns are described afterwards. Then, in 6, the motion of the gas phase is characterized with the processing of high-speed videos. It is related to the two-phase values through an analytical model of the forces exerting on the bubbles in 7.

Experimental setup

The experimental setup is comprised of several features that will be described separately for the sake of clarity.

Sliding tube: The section itself is a stiff aluminum square tube. The hydraulic diameter D h is 76.2 mm and the length L is 1 m. It is vertically mounted and supported by four linear bearings installed on two parallel shafts. This feature allows the tube assembly to slide in the y-direction only, transverse to the flow (cf. 1). The linear bearings are self-aligning to minimize friction in the system. Three panels of the tube are transparent, and the square section prevents refraction effects, allowing a good visualization of the flow.

Compression spring

Square tube

Excitation wheel

Extension spring

Two-phase flow inlet

Variable-speed motor z x y

Cable pulley Compression system: Four compression springs of stiffness k c and concentric with the rods retain the test section in the y-direction. One of their ends is in contact with the tube assembly and the other with what we refer to as compression plates. These are bolted with the sturdy bench structure via threaded rods (see 2). They can be displaced along those threaded rods to set an initial compression to the springs. This is an important feature since it prevents the tube assembly from impacting the two sets of springs when oscillating. The four springs are always working. However, the maximum amplitude that the system can reach is equal to the length of the initial compression. Beyond that point, the tube assembly would disconnect from a pair of springs. The compression springs used in the present study have a free length of 152 mm and the initial compression is 45 mm. The tube assembly behaves like a classic mass-spring system, from a structural point-of-view.

Excitation system: The "mass-spring system" is forced with a simulated sinusoidal excitation, provided through extension springs of stiffness k 0 . A motor drives an excitation wheel (cf. 1) via a pulley. A cable is attached to a shoulder-bearing system screwed eccentrically on the excitation wheel. The rotary motion is converted via the cable-pulley system into a sinusoidal motion of the extremity of extension springs (cf. 1). Since the extension springs are always stretched, the tube undergoes a sinusoidal force. We can control its magnitude F 0 by changing the eccentricity e on the excitation wheel. The latter is drilled with multiple tapped holes at different positions. The excitation frequency ω can be changed using the variable-frequency drive that pilots the motor. Notice the symmetry of the system with respect to the yz-plane: it is excited through two springs. Two additional extension springs are supporting the tube on the left side of the bench on 1, to ensure symmetry with respect xz-plane. All the springs, including compression springs, were tested on a traction-compression machine in order to verify their linearity over their range of operation. Dynamic stiffnesses have not been tested but are neglected since the excitation frequencies are very low (around 6Hz).

Hydraulic loop: In order to allow the transverse motion of the oscillating tube under internal flow, we use two flexible hoses that are fitted on both ends of the tube, and to the rest of the rigid pipe system on the other end. The flexible hoses are corrugated to ensure local pressure resistance and global flexibility, so as to not affect the tube motion. Also, their lengths were chosen so as to avoid an unbalancing mass effect at mid-section induced by water presence and make sure that they "follow" the tube smoothly.

The hydraulic loop works as follows: a centrifugal pump takes water from a reservoir. Air is injected just upstream the test section, as illustrated in 3. A mixer placed in the circular to square expansion homogenizes the air-water mixture. The mixture flows upward through the test section and then back to the tank at atmospheric pressure, where the phases separate.

Measurement system: The position of the tube is measured in time with laser sensors. They yield a resolution lower than 12 µm. The acquisitions are performed at a sampling rate of 2048 Hz.

Water

Air injection

Flexible tube

Experimental parameters

Two-phase flow parameters

Two-phase flow in a tube can be characterized either by the void fraction ε in a portion ∆L or by the volumetric fraction β [ [START_REF] Collier | Convective boiling and condensation[END_REF]]. The void fraction represents the proportion of gas volume over the total volume:

ε = V g V g + V l = A g ∆L A g ∆L + A l ∆L = A g A g + A l (1)
where A k is the area occupied by phase k in a given section of the tube. The void fraction would require specific instrumentation such as capacitance or fiber optic probes to be accurately determined.

On the other hand, the volumetric fraction only requires the volume flow rates of each phase:

β = Q g Q g + Q l = A g u g A g u g + A l u l = A g A g + A l /s (2)
The volume and area void fractions are equal for fully developed two-phase flows. where s is the slip ratio between the average velocities of the phases: s = u g / u l . Note that ε = β if s = 1, which is the definition of a homogeneous flow. We characterize the mixture velocity using the definition of superficial velocity:

j = Q g + Q l A (3)
The experiments were performed at constant j, using the fact that j g = βj and j l = (1 -β)j.

Volumetric flow rates are both measured with appropriate flow meters.

Sources of damping

The presence of fluid around a structure will affect its damping. [START_REF] Carlucci | Experimental studies of damping and hydrodynamic mass of a cylinder in confined two-phase flow[END_REF] identified several sources of damping:

ζ t = ζ s + ζ f + ζ ν + ζ 2ϕ (4)
(i) ζ s : Structural damping is caused by the energy losses inherent to the motion of a structure.

In our case, it is due to the friction in the linear bearings and inner losses in the springs. Its value is obtained by determining the oscillating characteristics of the system while no fluid is flowing through (j = 0). (ii) ζ f : When water is flowing through a slender pipe, the latter can warp and oscillate under fluid force. This force causes a damping effect at the elbows: it is the fluid damping mechanism. Since our structure is practically stiff, ζ f can be neglected.

(iii) ζ ν : Fluid friction on the sides of the test section induces viscous damping, which depends on fluid properties, fluid velocity and channel geometry. As can be seen on 4, it is negligible at 100% void fraction, and is monotonic until 0% (water only). Even there, ζ ν is small compared to ζ 2ϕ . We determine ζ ν at ε = 0%. Then, it is assumed that ζ ν decreases linearly with ε and is null at ε = 100% for a given superficial velocity of the mixture. (iv) ζ 2ϕ : The two-phase component is obtained by subtracting the aforementioned components from the measured total damping ζ t . This formulation suggests that all the components are purely velocity dependent. This is not the case for this system, as dry friction occurs in the linear bearings. 4 aims at describing how the two-phase component is therefore extracted.

Experimental technique 4.1 Structural damping

Information on the nature of structural damping is determined by analyzing the free vibrations of the system. The empty tube was initially displaced by 20 mm and then released. The result is shown on 5. The linear envelope is typical of Coulomb friction, occurring in the linear bearings. It is explained by the fact that a constant amount of energy is extracted from the system at each cycle, meaning that it is not velocity dependent. Since in this case, we know the weight of the tube assembly, we can calculate a friction coefficient µ of 0.06 per bearing. This value is inside the expected range [0.04, 0.07] for linear bearings [ [START_REF] Van Beek | Advanced engineering design: lifetime performance and reliability[END_REF]], where friction is indubitably higher than in roller bearings. In order to test the effect of mass and frequency, the tube assembly was ballasted with several masses or stagnant water, and different springs were used. As a result, total mass m over structural mass m s was varied within [1, 1.6] and natural frequency within [START_REF] Béguin | Influence of viscosity, density and surface tension on two-phase damping[END_REF][START_REF] Zuber | On the dispersed two-phase flow in the laminar flow regime[END_REF] rad/s. Coulomb friction coefficient was found to vary within [0.04, 0.08].

In theory, pure Coulomb friction allows an infinite response at resonance for a mass-spring system, for low values of µ such as in our case [ [START_REF] Den | Mechanical vibrations[END_REF]]. Therefore, a combination of Coulomb friction and viscous damping is commonly used. A more elaborate friction model [ [START_REF] Olsson | Friction models and friction compensation[END_REF]] did not prove to be necessary, as the tube is constantly sliding (no stick-slip phenomena, for instance, are involved). We then introduce ζ s , the viscous component due to structural losses. The determination of µ and ζ s , along with the other sources of damping, requires an appropriate modeling of the system for operating conditions.

System modeling

A top view of the test section is provided on Fig. 6(a). We can identify the four compression springs and the two pairs of extension springs, clamped on the right hand side and excited on the left hand 6(b) shows the full system model from a structural point-of-view. The mass of the tube assembly is noted m s and its position y s . The springs are considered to be linear, based on the tractioncompression tests mentionned earlier in this paper. The excitation springs are subject to a sinusoidal displacement. We isolate the forces (including the different sources of damping and friction) on the free body diagram of the tube assembly. Remembering that because of the initial compression, all the springs are always working, the equivalent variables of the system under operating conditions (6(c)) are:

Force magnitude : F 0 = 2k 0 e (5)
Stiffness :

k t = 4k 0 + 4k c (6)
Damping :

c t = c s + c ν + c 2ϕ (7)
Mass :

m t = m s + m 2ϕ + m a (8)
The total damping coefficient c t is the sum of the three aforementioned components. m 2ϕ is the mass of fluid vibrating with the structure and can be deduced from ε and tube geometry. m a is the added mass, due to a relative motion between the liquid and the structure. The equation of motion of the structure in direction y, for a given excitation frequency ω, is:

m t ÿs + c t ẏs + k t y s = F µ ( ẏs ) + F e (t) (9)
= -µm t g sgn( ẏs ) + 2k 0 e cos(ωt) [START_REF] Béguin | Influence of viscosity, density and surface tension on two-phase damping[END_REF] It can be rewritten as:

ÿs + 2ζ t ω n ẏs + ω 2 n y s = -µg sgn( ẏs ) + 2k 0 eω n k t cos(ωt) (11)
with the damping ratio ζ t and the natural frequency ω n defined as usual as:

ζ t = c t 2 √ k t m t and ω n = k t m t (12)
The total damping ratio is non-dimensionalized by the total mass, including the added mass. In this case, the added mass is negligible compared to m s and m 2ϕ . The relative error is estimated to be lower than 3% on ζ t . 11 is non-linear because of the friction force. Thus, finding an analytical expression for ζ t is far from straightforward. The next paragraph describes the method to retrieve variables of 12 from the frequency response function of the system.

Protocol

The position of the tube is acquired with the laser sensors. The time domain samples have a duration of 20 seconds. A typical tube response is shown on 7. The RMS amplitudes of tube Y rms s and excitation frequencies ω are extracted from these samples. This operation is repeated by increasing ω with the variable frequency drive in order to cover the resonance peak of the structure in the frequency spectrum. The frequency resolution that we get on the tube (including power transmission ratios) is 0.03 Hz. This frequency sweep allows to construct the frequency response function of the tube, as illustrated on Fig 8.

The black and gray curves calculated for a linear system (ζ = 2.2% and ζ = 6.0% respectively) have been overlaid to illustrate the slight non-linearity of the system. Around resonance, the experimental points are closer to the black curve whereas when excited farther away from the natural frequency, they are closer to the grey curve.

The natural frequency of the system is extracted from a polynomial fit around the resonance point. Then, we use the least-squares method to find values for µ and ζ t of 11. Indeed, we do not know the weight supported by the bearings when there is flow. The bottom flexible tube may pull the tube assembly down whereas the flow going up might reduce the force on the bearings. Since the peak of resonance is well defined (ζ t is expected to be smaller than 6%), only a narrow range of excitation frequency is required to accurately determine the damping. Therefore, we consider µ to be constant over a frequency response function for given conditions. So, 11 is solved for all the frequencies tested, the RMS amplitudes are calculated and compared with the experimental values. Note on 10 that the model suggests that the total mass m t is supported by the bearings. Since a smaller value for the mass is expected in reality (as explained in the previous paragraph), the compensation is done directly on the friction coefficient. The values of µ for the conditions tested lie within [0.025, 0.045], just below the range expected from the free vibration tests, where the total mass of the system was supported by the bearings in this case. For the experiments with two-phase flow, where the reaction supported by the bearings is unknown, the parameter µ should be interpreted as the ratio of friction forces on the overestimated reaction = f µ /(m t g). In fact, the values of µ vary depending on the flow conditions. For the case presented on Fig. 7, if we force µ rather than fitting it, we obtain ζ t = 1.6% underestimating the friction force with µ = 0.025 and ζ t = 1% overestimating the friction force with µ = 0.045 presented in doted and dashed line. The best fit is ζ t = 1.4% for µ = 0.033 closest to the experimental condition. We can therefore estimated the absolute experimental error as ±0.3%

This confirms that the non-linearity of the system can be explained by Coulomb friction induced by the bearings only. Furthermore, the good agreement between the model and experiments shows that two-phase damping can be modeled as a velocity dependent damping.

Results

Influence of fluid velocity

Total damping results are shown on Fig. 9 as function of volumetric fraction β, for the three superficial velocities of the mixture tested. Each point is the average over three measurements. First-of-all, we can see for β = 0% that ζ s + ζ ν is fairly small, of the order of 0.5 to 1%. Then, damping seems proportional to volumetric fraction until a change in slope occurs (vertical dashed line). This transition corresponds to the transition between bubbly and churn regimes, as will be discussed in 5.2. The breakdown in damping values was related to the change in interface surface area by [START_REF] Gravelle | Damping of tubes due to internal twophase flow[END_REF]. As anticipated, the two-phase damping values are very high and reach 3%. However, damping is expected to go down for high volumetric fraction, since the flow becomes air single-phase. This behavior can be explained with 10.

10 shows the evolution of the natural frequency of the system. When β is increased, the mass of the vibrating fluid decreases. Thus, the natural frequency of the total system increases. The dotted line represents the expected natural frequency in the context of a perfectly homogeneous flow, neglecting pressure and added mass effects. It can be observed that when the volumetric fraction is increased in the experiments, the response gets farther away from the line. On the other hand, when superficial velocity is increased, results get closer to the line. This proves that the higher the velocity, the more homogeneous the flow (void fraction is closer to volumetric fraction). The relation between volumetric fraction and void fraction can be summarized by this equation: When s tends to 1, β becomes closer to ε. Physically, in vertical co-current upward flow, gas is faster than liquid due to drag-buoyancy equilibrium. This equilibrium creates a velocity difference that is not much affected by liquid velocity. Therefore, as the liquid velocity gets higher, s tends to 1. For higher void fractions, flow is intermittent because there is no definite continuous phase. The large diameter of our tube amplifies this effect. Stagnant volumes of water are observed and are subject to sloshing which causes an increase in slip ratio. A superficial velocity of 0.7 m/s is the maximum velocity that can be reached with the current pump.

s 1 β -1 = 1 ε -1 (13)
In the meantime, ω n is a good qualitative indicator of the actual void fraction in the test section, as described in 5.3. The flow conditions at β = 90% are estimated to correspond approximately to ε = 60%.

Flow patterns

Videos of the flow were taken with a high-speed camera. Selected photographs of each volumetric fraction are shown on 11 for a superficial velocity of 0.7 m/s. Bubbly flow can clearly be observed up to β = 40%. The breakdown in damping values occurs during the transition to churn regime. Perturbation due to larger bubbles can be made out for β = 50% -60%.

We note that no slug flow regime is observed. That is because the length to diameter ratio is L/D h = 13, whereas it should be higher than 20 for slugs to "develop". This had been observed by [START_REF] Ohnuki | Experimental study on transition of flow pattern and phase distribution in upward air-water two-phase flow along a large vertical pipe[END_REF] who performed void fraction measurements on a very large circular tube of 200 mm inner diameter, and flow conditions comparable to ours. It was also concluded that churn flow was dominant in large pipes for conditions where slug flow exists in smaller pipes. They also reported a phenomenon that can be clearly observed at 80% volumetric fraction. Large Taylor bubbles intermittently perturb a liquid film still filled with many small bubbles. [START_REF] Prasser | Influence of the pipe diameter on the structure of the gas-liquid interface in a vertical two-phase pipe flow[END_REF] note that those large bubbles can however freely move and deform in three dimensions in the large tube and are thus far from ideal Taylor bubbles.

As reviewed by [START_REF] Coleman | Characterization of two-phase flow patterns in small diameter round and rectangular tubes[END_REF], flow patterns in large circular tubes have been experimented over the years, mostly for horizontal ducts. Flow patterns maps [ [START_REF] Mandhane | A flow pattern map for gas-liquid flow in horizontal pipes[END_REF]] were derived and compared reasonably well with theoretical ones [ [START_REF] Taitel | A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow[END_REF]]. Rectangular microchannels and minichannels have also been studied and reviewed thoroughly by [START_REF] Cheng | Two-phase flow patterns and flow-pattern maps: fundamentals and applications[END_REF]. Unfortunately, very little information exists on large rectangular channels. Aspect ratio, surface tension, hydraulic diameter and pressure have an effect on transitions in small (around 5 mm) to micro channels. However, these conclusions may not be verified for large vertical square ducts.

An extensive study of flow pattern transition determination would require additional instrumentation. It is beyond the scope of this study. Our observations of the flow are meant to support damping results and give away the lack of studies in this range of diameter for vertical rectangular ducts. 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%

Influence of excitation force

The magnitude of excitation force was doubled for j = 0.7 m/s. In order to ensure repeatability with respect to volumetric fraction, natural frequency is used to compare the results. Indeed, if we neglect added mass, ω n is related to hydrodynamic mass m 2ϕ which depends directly on void fraction ε:

     ω n = k t m s + m 2ϕ m 2ϕ = (AL + V corr )(ερ g + (1 -ε)ρ l ) (14)
where V corr is a correction volume that corresponds to the additional volume of vibrating fluid within the flexible tubing. It was determined for ε = 0%. ρ g depends on average pressure in the tube. The pressure upstream was found to be below 1.6 bar for all flow conditions. Given the high pressure loss due to the mixer on 3 and the short horizontal length of tubing downstream, the average pressure is close to atmospheric conditions in the oscillating tube. Therefore, ρ g was determined at atmospheric conditions, for 23 • C. The results are shown on 12.

Clearly, the excitation force has an effect on transition since the two curves separate after the transition. The higher the force, the higher the confinement which is likely to affect flow pattern and the damping values for moderately high void fraction. However, it is interesting to notice that the curves match perfectly before transition, for bubbly flow conditions. This independence of ζ 2ϕ with respect to F 0 suggests that the damping mechanism can be modeled as a purely velocity dependent mechanism for bubbly flow. In other words, ζ 2ϕ is independent of excitation amplitude and frequency for bubbly flow. Furthermore, we have depicted the direct dependence of two-phase damping with interface surface area, so we did not expect ζ 2ϕ to depend on force magnitude. This motivated us to steer the study towards the modeling of two-phase damping mechanism for bubbly flow in 7.

Gas phase behavior 6.1 Video processing

The direct relation between interface surface area and two-phase damping is not sufficient to explain the two-phase damping mechanism. We suspect that relative motion between the phases induces dissipation, through the work of the forces exerting on the bubbles.

We used a high-speed camera at 1000 fps to characterize the gas phase motion for bubbly flow within the oscillating structure. The objective is to correlate the gas phase motion to the damping values measured. To do so, β should be sufficiently high (for lower volumetric fraction, two-phase damping would be too small to be measured accurately). Unfortunately, the higher the void fraction, the harder the characterization. For instance, [START_REF] Lindken | A novel PIV technique for measurements in multiphase flows and its application to two-phase bubbly flows[END_REF] presented an accurate measurement of the flow field with bubbles using a combined PIV/shadowgraphy technique. However, the local void fraction around the cluster of bubbles does not exceed 2.5%. The challenge in our case is that many interbubble motions were observed, including overlapping. For this reason, individual tracking of the bubbles with image segmentation techniques is laborious.

Therefore, we resorted to a global tracking of the gas phase motion. The method consists in selecting a thumbnail of an image at instant t, and looking for the same thumbnail in the next image at time t + ∆t. We note f (t) the full image matrices taken with the high speed camera, of constant size (M, N ) at instant t. The indices of the thumbnail position are noted (i(t), j(t)). The new coordinates indices (i(t + ∆t), j(t + ∆t)) are calculated in order to minimize the Mean Square Error (MSE), expressed as:

M SE = M/2<I≤M/2 -N/2<J≤N/2 [f (i(t) + I, j(t) + J) -f (i(t + ∆t) + I, j(t + ∆t) + J)] 2 (15)
The high frame rate was set so that the minimization could be performed in a region of pixels close to the original thumbnail. The thumbnail should of course be initialized in the middle of the gas phase. At each time step, the new thumbnail is used to look for the next one. This evolution of f (i(t), j(t)) tends to average the inter-bubble motion and gives a good representation of the bulk gas phase motion.

The resolution is around 1 pixel = 0.2 mm. While very simple and robust, the drawback of this technique is that its accuracy can only be appreciated qualitatively.

On 13, one can directly observe the gas phase relative motion with respect to the structure. The images are shown in the tube reference frame. On the first image, the test section is accelerating towards the left, and the gas column is already compressed on the left wall of the tube. This suggests that bubbles are in phase lead with respect to the structure. On the second frame, the tube reaches a zero velocity point, so the gas phase is uniformly dispersed, before being compressed on the right. Notice the white rectangle, standing for the tracked thumbnail. It has moved not only in the transverse direction but also upward. Indeed, upward co-current air-water mixtures are studied.

Characterization of gas phase motion

Experiments for three volumetric fractions β = 5, 10 and 15%, at constant superficial velocity j = 0.6 m/s. Since the shape of the thumbnail does not change, only one point of its coordinates (i(t), j(t)) is required to track the bulk gas phase motion with respect to time. i(t) and j(t) correspond respectively to the vertical and transverse motions. Vertical relative velocity with respect to the liquid u gz -u lz , ranges from 0.22 to 0.28 m/s. This is around the expected value of 0.25 m/s (given by the difference between buoyancy and drag force on the bubbles). Those values give a slip ratio s of approximately 1.5. The three corresponding void fractions ε are thus roughly 3, 7 and 11%.

A typical graph of transverse relative motion of the bubbles is shown on 14. The signal is the combination of a sine wave and a linear deviation, meaning the tracking rectangle is oscillating while drifting towards one side of the tube. Note that only bubbles close to the front window are observed. We attribute the drift behavior to a swirling motion of the bubbles. When the structure is oscillating, a radial instability is onset (the swirling motion was not observed for the tube at rest). When bubbles come into contact with the side walls, they tend to slide in a preferential direction instead of simply bouncing back and most likely coalesce with other bubbles. This is attested by the photographs in 13 where structural oscillations do not seem to prematurely onset the formation of larger bubbles. A similar behavior had been reported by [START_REF] Béguin | Two-phase damping and interface surface area in tubes with vertical internal flow[END_REF]. When their clamped-clamped tube was released to measure the free vibrations, they observed elliptical motion instead of oscillations in a plane. This phenomenon was attributed to the cylindrical shape of the tube causing swirling. Our observations suggest that swirling is generated by tube oscillations only, and not by tube geometry. The swirling velocities on the side wall were around 0.05 m/s (slope of the dashed line on 14), which is negligible compared to up to 1 m/s for structure peak velocity. Thus, the swirling motion is not dissipative. However, it prevents coalescence of the bubbles, keeping a constant interface surface area. It could explain why tube oscillations do not seem to affect flow pattern transitions.

Therefore, only the harmonic motion of the gas phase is extracted. The tube amplitude is noted Y s and the relative amplitude of motion of the gas phase with respect to the tube is Y b|s . The corresponding time-dependent motions are noted y s and y b|s . Fig. 15 shows the different movements of the system involved. Note that the gas phase is in phase lead with respect to the structure.

Amplitude of motion of the gas phase

The relative amplitude Y b|s was extracted for several excitation amplitudes. Y s was varied by changing the frequency close to the natural frequency of the system for a given void fraction. Since the peak of resonance is very well defined, the excitation frequency can be considered constant at 33.75 ± 1.5 rad/s (6 ± 0.23 Hz) over all our experiments. Results are shown on 16. Y b|s /Y s represents the gain of the gas in terms of amplitude. For a constant excitation amplitude, the gain of the gas tends to decrease when void fraction increases. This is in accordance with the drag relations that take into account the effect of void fraction [e.g. [START_REF] Ishii | Drag coefficient and relative velocity in bubbly, droplet or particulate flows[END_REF], [START_REF] Zuber | Steady state and transient void fraction of bubbling systems and their operating limits (part i, steady state operation)[END_REF]]. Indeed, it is well-established that C D increases dramatically with ε through confinement, therefore limiting the transverse amplitude. For high Y s , the gain seems to reach an asymptote. As explicitly shown in 7.1, important shape deformation for high excitation amplitudes cause an increase of the drag force on the bubbles. Also, it is well established that wall proximity causes a similar effect on the drag force [START_REF] Takemura | Drag and lift forces on a bubble rising near a vertical wall in a viscous liquid[END_REF]. In this case, the tube walls are moving and thus represent a physical barrier that cannot be crossed by the bubbles, leading to an important confinement effect. All these effects contribute to the existence of the limit cycle.

The relative motion can not be related to two-phase damping values without the forces on the bubbles. Therefore, a simple model has to be analytically derived.

Analytical model

We propose a simple model of a bubble in an oscillating structure subjected to internal two-phase flow. This sections aims at testing existing correlations and implement them in the model under certain hypotheses, to assess the extent to which the two-phase damping can be reproduced. Based on our observations, the transverse relative motion of the gas phase will be the main output of the calculations to explain damping, through the work of the forces exerting on the bubbles.

Equation of motion

We consider a deformable bubble of equivalent diameter a. It is immersed in an oscillating structure filled with water, as illustrated by 17. The bubble has the properties of a confined bubble in two-phase flow, of given volumetric fraction and superficial velocity. Thus, its movement is representative of the gas phase motion.

We note #" u b the relative velocity of the bubble with respect to the fluid:

#" u b = #" u g -#" u l (16)
We assume that bubbles do not have influence on the liquid. Thus, the liquid undergoes a solid body motion with the structure. Therefore:

#" u l = (-Y s ω sin(ωt), u lz ) (17)
where u lz is a function of j and ε, the void fraction still being calculated based on volumetric fraction and slip ratio of the bubble. The 2D equation of motion of the bubble takes the form: where #" F B is the buoyancy force:

m d #" u b dt = #" F B + #" F D + #" F M + #" F I (18)
#" F B = 4 3 πa 3 ∆ρg #" z (19)
#" F D is the drag force, based on bubble relative velocity:

#" F D = - 1 2 πa 2 ρ l C D u b #" u b (20)
We calculate the drag coefficient using a relation by [START_REF] Béguin | Modélisation des écoulements diphasiques: amortissement, forces interfaciales et turbulence diphasique[END_REF]. It is valid for up to Re = 300, and takes the effect of void fraction into account:

C D = 16 Re      1 + 2 2 + 3µ * 2 + 2µ * 2 1 + Re c Re      1 -ε (1 -ε 1/3 ) 3 P 1 + µ * P 2 P 3 + µ * P 4 G(χ)
with :

P 1 = 4 + 6ε 5/3 P 2 = 6 -6ε 5/3 P 3 = 4 + 6ε 1/3 + 6ε 2/3 + ε P 4 = 4 + 3ε 1/3 -3ε 2/3 -4ε Re c = 33 + 8600ε 2/3 G(χ) = 1 3 χ 4/3 (χ 2 -1) 3/2 χ 2 -1 -(2 -χ 2 )sec -1 χ (χ 2 sec -1 χ -χ 2 -1) 2 (21)
The oblateness χ of the bubble is calculated with a terminal velocity correlation by [START_REF] Kushch | Finite-weber-number motion of bubbles through a nearly inviscid liquid[END_REF] (first-order term only):

χ = 1 + 9 64 We ( 22
)
where the Weber number We is defined as:

We = 2aρ l u 2 b σ (23)
The relation of 22 is valid for a single bubble rising in stagnant liquid, so its application is questionable in our case. For instance, oblateness tends to increase with void fraction ( [START_REF] Duineveld | The rise velocity and shape of bubbles in pure water at high reynolds number[END_REF], [START_REF] Roig | Measurement of interstitial velocity of homogeneous bubbly flows at low to moderate void fraction[END_REF]). However, our observations show that bubbles change shape depending on structure position, and so that velocity is the major parameter controlling the bubble deformations.

#" F M is the added mass force:

#" F M = - 4 3 πa 3 ρ l C M d #" u b dt (24)
where C M is the added mass coefficient. A few studies propose the effect of void fraction on the added mass of a bubble. See for instance the works by [START_REF] Zuber | On the dispersed two-phase flow in the laminar flow regime[END_REF] or [START_REF] Cai | A more general cell model for added mass in two-phase flow[END_REF]. However, added mass could increase or decrease with ε depending on certain hypotheses. A more recent study by [START_REF] Béguin | Void fraction effect on added mass in bubbly flow[END_REF] proposes a new semi-empirical correlation based on potential flow theory for random clouds of bubbles. It appears that C M weakly increases with ε. Therefore, we will simply use the added mass coefficient C M = 0.5 of an isolated sphere as a first approximation. Finally, the bubble undergoes an inertia force #" F I because of the pressure gradient caused by the accelerating tube on the water, in the transverse direction. This force can be expressed as:

#" F I = 4 3 πa 3 ρ l #" y s (25)
To sum up, the only input parameters for the model are a, j and β. The motion of the bubble is governed by its added mass and fluid viscosity. Thus, it has no mass nor rigidity, as modeled by [START_REF] Hara | Analytical model for evaluating added mass and damping of a vibrating circular rod in two-phase fluid[END_REF].

Transverse amplitude

A few images of the oscillating structure at different instants were segmented in order to have an estimate of the bubble sizes and shapes, for Y s /D h = 0.22. The bubble radii a were found to range within [1.2, 2] mm and oblateness χ within [1.2, 1.9], assuming a revolution ellipsoid. For the same conditions and a = 1.4, the model predicts χ up to 2.3, which seems reasonable considering the fact that the correlation we use do not take the effect of acceleration into account.

The predicted bubble relative amplitude Y b|s is presented on 18(a). The model tends to overpredict the amplitude by a factor 2. Thus, some effects have obviously been overlooked. Still, the [START_REF] Coleman | Characterization of two-phase flow patterns in small diameter round and rectangular tubes[END_REF][START_REF] Mandhane | A flow pattern map for gas-liquid flow in horizontal pipes[END_REF][START_REF] Taitel | A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow[END_REF][START_REF] Cheng | Two-phase flow patterns and flow-pattern maps: fundamentals and applications[END_REF][START_REF] Lindken | A novel PIV technique for measurements in multiphase flows and its application to two-phase bubbly flows[END_REF][START_REF] Ishii | Drag coefficient and relative velocity in bubbly, droplet or particulate flows[END_REF][START_REF] Zuber | Steady state and transient void fraction of bubbling systems and their operating limits (part i, steady state operation)[END_REF][START_REF] Takemura | Drag and lift forces on a bubble rising near a vertical wall in a viscous liquid[END_REF] trend seems very well respected. The limit cycle reached is mostly due to the oblateness of the bubbles, as attested by 18(b).

Dissipated energy

It is also useful to compare the model directly to the damping values measured for bubbly flow. We decided to compare the energy dissipated by the two-phase damping force F 2ϕ over one cycle of oscillation:

E 2ϕ = 2π/ω 0 F 2ϕ dy s dt dt = 2π/ω 0 c 2ϕ ẏs 2 dt = πc 2ϕ ωY 2 s ( 26
)
This energy has to be compared with that dissipated by the work of the forces applied on the bubbles. We are only interested in the projection in the y direction of these forces since we try to explain two-phase damping with transverse motion of bubbles. The buoyancy force is vertical and as such is not accounted for. Added mass and inertia forces are purely inertial, hence they do not work over one cycle. Only the work of the drag force projected on the y-axis is contributing, and can be written as:

W D = N b 2π/ω 0 # " F D • #" y dy b|s dt dt (27)
N b stands for the number of bubbles in the test section. It is calculated based on the void fraction:

N b = εAL 4 3 πa 3 (28)
This equation is valid only for bubbly flow (up to β ≈ 50% experimentally), so no effect of coalescence is considered. The void fraction depends on volumetric fraction and slip ratio (previously defined in the z direction). The model compares fairly well with two-phase damping in terms of dissipated energy. It is odd to have a good agreement for the energy, and a poor one for relative amplitude of the bubbles as reported in 7.2. Assuming that the drag relation is correct, we believe this is due to the fact that we measured relative amplitude of the gas phase with respect to the structure, and not the liquid, as underlined in the next section.

Model calibration with glycerin experiments

Changing the Reynolds number is a good way to test the validity of the model for other flow conditions. Other videos were taken with a stagnant glycerin solution (j = 0). Glycerin density is 1.21. Its viscosity was tested with a rheometer and was found to be of 0.163 Pa.s at 22 • C. Single bubbles were injected with a needle. The Reynolds numbers based on bubble relative velocities were below 3. Using segmentation imaging techniques, we were able to measure bubble radii, transverse amplitude and vertical terminal velocity (cf. 20(a) and 20(b)).

The scatter in the data is explained by the fact that we do not control the bubble radius. This is why the points were gathered in several radius categories, thus collapsing the data. The structural transverse amplitude does not have a strong effect on Y b|s . Given the low Weber numbers of the experiments, the bubbles hardly deform (χ ≈ 1). Hence, no additional drag is induced. Also note that the graphs are basically identical on their respective scale. Experimentally, the bubbles would rise in glycerin in a straight line due to the very low Reynolds number (two path instabilities occur at higher Re, as reported extensively in the literature, e.g. [START_REF] Tomiyama | Terminal velocity of single bubbles in surface tension force dominant regime[END_REF] and [START_REF] Mougin | Path instability of a rising bubble[END_REF]). Hence, at this range of Re, transverse migration is only due to tube motion and depends directly on bubble radius, as attested by 21. The figure summarizes the twenty-five experiments in stagnant glycerin. Since all the points collapse on a straight line, the relative amplitude of the bubble is only a function of its radius. However, two bubbles with different radii have different terminal velocities in the z direction, as attested by Fig. 19(b). This confirms that no coupling exists between the y and z directions in this case.

Contrary to 7.2, the agreement between the model and experiments on 20 is very good, for both U t and Y b|s . This suggests that considering the relative amplitude with respect to the structure or the liquid is the same. Hence, assuming solid body motion for the liquid is correct in this case. The poor agreement reported before for air-water mixture experiments may be caused by a recirculation of water around bubbles or a wrong drag coefficient given the bubbles' high Reynolds number. Indeed, the higher the void fraction, the more important the recirculation to fill the space between the moving bubbles. This interstitial flow would affect both drag on bubbles and relative velocity. Furthermore, we stated in 6.2 that structural oscillation did not seem to onset bubble coalescence. The nature of the liquid film between bubbles would control their coalescence/break-up conditions [ [START_REF] Razzaque | Bubble size in coalescence dominant regime of turbulent air-water flow through horizontal pipes[END_REF]]. The local gas-liquid interaction due to the relative bubble motion is given by pseudoturbulence theory but is beyond the scope of this study. Bubble-induced liquid agitation was modeled [e.g. [START_REF] Elghobashi | A two-equation turbulence model for two-phase flows[END_REF] and [START_REF] Béguin | Modélisation des écoulements diphasiques: amortissement, forces interfaciales et turbulence diphasique[END_REF]] and characterized experimentally at moderate void fractions [ [START_REF] Lance | Turbulence in the liquid phase of a uniform bubbly air-water flow[END_REF] and [START_REF] Riboux | Experimental characterization of the agitation generated by bubbles rising at high reynolds number[END_REF]]. Information on local interstitial velocity in an oscillating structure with gas bubbles would allow to model turbulence forces on bubbles as well as lift forces (through information on local vorticity) which were neglected in the present model. This effect is important even for two-phase flow in tubes at rest for void fractions larger than 2% [ [START_REF] Roig | Measurement of interstitial velocity of homogeneous bubbly flows at low to moderate void fraction[END_REF]]. Thus, it is likely that structural oscillations cause additional pseudo-turbulence when inducing the relative motion of the gas phase. Characterizing this phenomenon might help to understand the two-phase damping mechanism more thoroughly.

Conclusion

In this paper, a new test section offering control over the excitation parameters in order to determine two-phase damping experimentally was presented. Observations as well as processing of high-speed videos gave novel information on the gas phase motion, especially on its relative motion with respect to the structure. A simple analytical model fed with correlations was derived. Although not entirely complete, the model gives useful information on the physical dissipative mechanism. It also underlines the missing information to catch the full nature of the phenomenon. We suspect that the latter lies in the complexity of the liquid phase motion. So far, the following conclusions and perspectives can be brought out:

(i) Two-phase damping can reach 3% in this square 76.2 mm tube. It is fairly proportional to void fraction, until a change in slope at transition between bubbly and churn flow regime occurs. (ii) The frequency response function of the tube subjected to internal two-phase flow confirms that two-phase damping is a viscous damping (velocity dependent) mechanism. This is supported by the fact that ζ 2ϕ seems independent of the excitation force magnitude F 0 for bubbly flow. For higher void fraction, the mechanism is different. We suspect that energy is extracted by sloshed liquid phase within the continuous gas phase. Damping is affected by the stochastic nature of the flow at this regime, causing an increase of standard deviation in the results. (iii) There is definite relative motion of the gas phase with respect to the liquid for bubbly flow.

The bubbles have a bulk body motion due to the tube oscillations, and are in phase lead with respect to the structure. The oscillations seem combined with a swirling motion. It is not very dissipative but helps to prevent a premature coalescence of the bubbles, keeping a maximum interface surface area prompt to transfer energy between the phases and the structure. (iv) For bubbly flow, the power dissipated by two-phase damping is equivalent to that dissipated by the drag forces on the bubbles. However, the discrepancy between bubble amplitudes predicted by the model and those measured suggests a complex interstitial liquid flow between the bubbles under tube oscillations, affecting both drag and relative velocity. It seems like even though two-phase damping is observed on a large scale, small scale effects cannot be overlooked for ζ 2ϕ to be accurately modeled.

Figure 1 .

 1 Figure 1. Test section (pipe system not shown).

Figure 2 .

 2 Figure 2. Sliding system.

Figure 3 .

 3 Figure 3. Air injection upstream to the test section.

Figure 4 .

 4 Figure 4. Components of total damping in two-phase axial flow and experimental apparatus. Total and structural damping are measured, fluid and viscous damping modeled [4].

Figure 5 .Figure 6 .

 56 Figure 5. Free response of the empty tube after an initial displacement.

Figure 7 .

 7 Figure 7. Measured position signal of the tube.

  µ and ζ t are found to best fit the experimental data in the frequency domain. The result is summarized by the dotted line on Fig 8.

Figure 8 .

 8 Figure 8. Frequency response function of the system (F0 = 23 N, j = 0.7 m/s, β = 30%).

Figure 9 .

 9 Figure 9. Influence of fluid velocity on total damping.

Figure 10 .

 10 Figure 10. Natural frequency of the system for the flow conditions tested. The dotted line represents the expected natural frequency with the homogeneous model, obtained with 14 (fn = ωn/2π).

Figure 11 .

 11 Figure 11. Flow patterns.

Figure 12 .Figure 13 .

 1213 Figure 12.Influence of force magnitude. Natural frequency is taken as reference for comparison, since it represents the actual void fraction more accurately (cf. 10). Note that ζ2ϕ is independent of F0 for bubbly flow.

Figure 14 .

 14 Figure 14. Typical unfiltered transverse position of the tracking rectangle in time. It represents the gas phase motion on the front tube wall. The linear deviation is attributed to a global swirling motion of the bubbles.

Figure 15 .

 15 Figure 15. Different movements of the system involved. Light dashed lines are the the raw signals extracted from the videos, and solid lines are the corresponding sinusoidal fits.

Figure 16 .

 16 Figure 16. Relative amplitude of the gas phase with respect to the structure.

Figure 17 .

 17 Figure 17. Schematic of the bubble model in an oscillating structure subjected to two-phase flow.

  Comparison between model and experimental values for gas phase transverse amplitude.

Figure 18 .

 18 Figure18. Analytical model results solving eq.[START_REF] Coleman | Characterization of two-phase flow patterns in small diameter round and rectangular tubes[END_REF][START_REF] Mandhane | A flow pattern map for gas-liquid flow in horizontal pipes[END_REF][START_REF] Taitel | A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow[END_REF][START_REF] Cheng | Two-phase flow patterns and flow-pattern maps: fundamentals and applications[END_REF][START_REF] Lindken | A novel PIV technique for measurements in multiphase flows and its application to two-phase bubbly flows[END_REF][START_REF] Ishii | Drag coefficient and relative velocity in bubbly, droplet or particulate flows[END_REF][START_REF] Zuber | Steady state and transient void fraction of bubbling systems and their operating limits (part i, steady state operation)[END_REF][START_REF] Takemura | Drag and lift forces on a bubble rising near a vertical wall in a viscous liquid[END_REF]

Figure 19 .

 19 Figure[START_REF] Mandhane | A flow pattern map for gas-liquid flow in horizontal pipes[END_REF]. Energy dissipated over one oscillation cycle. -: by two-phase damping using 26; --: calculated with analytical model using 27. Only the power of bubble drag forces is significant in the model. Calculations are presented for structure amplitude Ys = 18.8 mm, superficial velocity j = 0.6 m/s and bubble radius a = 1.4 mm.

  Model: a = 0.9mm Model: a = 1.1mm Model: a = 1.3mm a < 1mm 1 ≤ a < 1.2mm a ≤ 1.2mm (b) Terminal velocity.

Figure 20 .

 20 Figure 20. Comparison between model and experiments with single bubbles rising in stagnant glycerin.

Figure 21 .

 21 Figure 21. Gain of the bubbles as a function of bubble radius, in stagnant glycerin (experimental).

Table 1 .

 1 Characteristics of the system.

	Property	Notation Value/Range
	Hydraulic diameter	D h	76.2 mm
	Tube area	A	D h × D h (square tube)
	Tube length	L	1 m
	Mass of tube assembly 2	m s	17 kg
	Compression spring constant 1 k c	6545 N/m
	Extension spring constant 1	k 0	474 N/m
	Natural frequency	f n	6.45 Hz
	Excitation force	F 0	23 -69 N
	1 Average of measured values over the four springs.
	2 Calculated with 12 for the empty tube.

  Comparison between the model and experiments is shown on 19.

				Two-phase damping Work rate of drag force on bubbles
		1.2
	Dissipated energy E (J)	0.4 0.6 0.8 1
		0.2
		0	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5
								β

Acknowledgments

This work was sponsored by the Natural Sciences and Engineering Research Council, Babcock & Wilcox Canada and Atomic Energy of Canada Ltd., through the BWC/AECL/NSERC research chair in fluid-structure interactions.