Alexis Fouilloy
email: fouilloy-a@univ-corse.fr

Cyril Voyant
email: voyant@univ-corse.fr

Gilles Notton
email: notton@univ-corse.fr

Marie Laure Nivet
email: nivet@univ-corse.fr

Jean Laurent Duchaud
email: duchaud_jl@univ-corse.fr

Machine Learning Methods for Solar Irradiation Forecasting: A Comparison in a Mediterranean Site

Keywords: Time Series forecasting, machine learning, variability, ARMA, ANN, Regression trees, Gaussian process, SVR

In this survey, several statistical and machine learning tools are analyzed and compared in view to forecast the solar irradiation in Ajaccio (Corsica, France, 41°55 N, 8°44 E, 4m asl). The forecasting horizon is from 1 to 6 hours with an hourly time granularity. Eleven forecasting models are compared: persistence, scaled persistence, ARMA, MLP, regression trees, boosted regression trees, bagged regression trees, pruned regression trees, random forest, Gaussian processes and support vector regression. The models are compared in term of error metrics: nRMSE (normalized root mean squared error), MAE (mean absolute error) and skill score related to the smart persistence.

INTRODUCTION

Solar radiation is one of the principal energy sources for physical, biological and chemical processes, occupying the most important role in many engineering applications [START_REF] Hoff | Modeling PV fleet output variability[END_REF]. The process of converting sunlight to electricity without combustion allows to create power without pollution. The major problem of such energy source is its intermittence and its stochastic character which make difficult their management into an electrical network [START_REF] Voyant | Multi-horizon solar radiation forecasting for Mediterranean locations using time series models[END_REF]. Thereby, the development of forecasting models is necessary to use ideally this technology [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF]. By considering their effectiveness, it will be possible for example to identify the most optimal locations for developing a solar power project or to maintain the grid stability and security of a power management system [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF]. Thus the solar energy forecasting is a process used to predict the amount of solar energy available for various time horizons [START_REF] Lorenz | Benchmarking of different approaches to forecast solar irradiance, 24th[END_REF]. Several methods have been developed by experts around the world and the mathematical formalism of Times Series (TS [START_REF] Join | Solar energy production: Short-term forecasting and risk management[END_REF]) has been often used for the short term forecasting (among 6 hours ahead) [START_REF] Lorenz | Benchmarking of different approaches to forecast solar irradiance, 24th[END_REF]. In this study, the solar forecasting is realized with an hourly time granularity for 1 hour to 6 hours' time horizon. The comparison between the models is evaluated by statistical indexes and graphical methods, the models are: persistence (P) [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF][START_REF] Voyant | Forecasting method for global radiation time series without training phase: Comparison with other well-known prediction methodologies[END_REF], scaled or smart persistence (SP) [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF][START_REF] Breiman | Classification and regression trees[END_REF], auto regressive moving average (ARMA) [START_REF] Breiman | Bagging predictors[END_REF], artificial neural network (ANN) [START_REF] Breiman | Classification and regression trees[END_REF][START_REF] Breiman | Random forests[END_REF], regression trees (RT) [START_REF] Rasmussen | Gaussian Processes in Machine Learning[END_REF], pruned regression trees (RT-pruned) [START_REF] Rasmussen | Gaussian Processes in Machine Learning[END_REF], boosted regression trees (RT-boosted), bagged regression trees (RT-bagged) [START_REF] Smola | A tutorial on support vector regression[END_REF], random forest (RF) [START_REF] Breiman | Random forests[END_REF], Gaussian processes (GP) [START_REF] Rasmussen | Gaussian Processes in Machine Learning[END_REF] and support vector regression (SVR) [START_REF] Smola | A tutorial on support vector regression[END_REF].

METHOD

Several preliminary stage must be realized before the forecasting phase:

-A quality control of the solar data: an automatic quality check used in the frame of GEOSS project (Group on Earth Observation System of System) (http://www.earthobservations.org/geoss.php) was applied to the data; -A data preprocessing: removing of the night hours;

sunset and the sunrise data induce some problems (mask effects and bad response of pyranometers at high zenith angle): a filter is applied on the datasets which remove all the data with a solar elevation angle up to 10°. The solar irradiation time series contains some properties like seasonality and periodicity which must be deleted to make the time series stationary; to do it, a clear sky model [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] is used and the clear sky index calculated (paragraph 2.1). -At last, the choice of the input data must be realized and is described in paragraph 2.2.

Clear Sky Irradiation (CSI)

The clear sky model is a theoretical model which determines the maximum theoretical global horizontal irradiation (GHI) for any geographical point at the ground level. The model used is SOLIS model [START_REF] Mueller | Rethinking satellitebased solar irradiance modelling: The SOLIS clear-sky module[END_REF][START_REF] Ineichen | Comparison of eight clear sky broadband models against 16 independent data banks[END_REF] based on radiative transfer models, Beer Lambert function, and their integration on the solar spectra uses the information's of news satellites ERS-2/ENVISAT. The calculation of the clear sly irradiation CS(t) is given by:

𝐶𝑆(𝑡) = 𝐻 0 . 𝑒 -𝜏 𝑠𝑖𝑛 𝑏 (ℎ(𝑡)) . 𝑠𝑖𝑛(ℎ(𝑡)) (1)
h(t) is the solar height in degrees and H0 is the extraterrestrial irradiation; τ is the global total atmospheric depth and the parameter b is a fitting parameter, both of them depending on the meteorological characteristics of the site [START_REF] Jiang | Mutual information algorithms[END_REF]. τ and b have different values for the three studied sites (https://aeronet.gsfc.nasa.gov/). The clear sky index (CSI) is calculated by:

𝐶𝑆𝐼(𝑡) = 𝐺𝐻𝐼(𝑡) 𝐶𝑆(𝑡) (2)
For all the machine learning tools used here, the stationary hypothesis is necessary and then the different time series, CSI are directly used as an input in the forecasting model. GHI is then obtained using Eq. 2 and all the errors metrics are given in term of GHI in absolute and relative values.

Dimension of the input matrix

We must choose the dimension of the input matrix of the forecasting tool. The approach used consists in predicting the future solar irradiation (at different time scales) based on the past observed data [START_REF] Diagne | Review of solar irradiance forecasting methods and a proposition for smallscale insular grids[END_REF]. Mathematically, the researched formulation is (ϵ(t+h) is a random white noise):

𝐺𝐻𝐼(𝑡 + ℎ) = 𝑓(𝐺𝐻𝐼(𝑡), 𝐺𝐻𝐼(𝑡 -1), … , 𝐺𝐻𝐼(𝑡 -𝑛)) + 𝜖(𝑡 + ℎ) (3)
The future time step (t+h) is forecasted based on the observed data at the times (t, t-1…, t-n). In other words, the objective is to calculate the value of n and to obtain ϵ as low as possible in absolute value. The choice of n, i.e. the dimension of the input matrix, is made by an auto mutual information method [START_REF] Jiang | Mutual information algorithms[END_REF][START_REF] Parviz | Mutual Information Based Input Variable Selection Algorithm and Wavelet Neural Network for Time Series Prediction[END_REF]. This auto mutual information is a property of the time series and depends on each dataset. It determinates the degree of statistical dependence of the variables specific to each site. Another step of preprocess is the k-fold sampling [START_REF] Salcedo-Sanz | Prediction of Daily Global Solar Irradiation Using Temporal Gaussian Processes[END_REF]: the dataset is divided in ten samples and each sample is used at least one time for the training and one time for the test, the simulation is repeated as many time as necessary. This method allows to have results independent of the set of data used for the training (only one data set being able to have some particularities that disturbs the robustness of the conclusions).

METEOROLOGICAL STATION AND DATA

The datasets used in this study are time series of horizontal global solar irradiation measurements (GHI) in Ajaccio (Corsica, France, 41°55 N, 8°44 E, 4m asl) at about 100 m from the Mediterranean sea (Figure 1). After the calculation of the auto-mutual information between the solar radiation data [START_REF] Jiang | Mutual information algorithms[END_REF][START_REF] Parviz | Mutual Information Based Input Variable Selection Algorithm and Wavelet Neural Network for Time Series Prediction[END_REF], it has been decided that the number of input data n in the forecasting model will be n=8.Thus, for estimating the solar irradiation at time t, the solar irradiations at time (t-1), (t-2), (t-3)… (t-n) will be used.

FORECASTING MODELS

Eleven models are used to forecast the global horizontal irradiations and are briefly described here. These models are classified in 3 categories: naïve models, classical machine learning models and regression trees based models. In this paper the symbol ̂ indicates that the value is predicted, without this symbol the value is measured.

Naïve models

The two next models are generally used only as a reference in view to compare it with more sophisticated models in terms of performances. The first model is the simplest model, the persistence, it is the repetition of the measure at the instant t to the instant t+h (h is the forecasting horizon) [START_REF] Diagne | Review of solar irradiance forecasting methods and a proposition for smallscale insular grids[END_REF]:

𝐺𝐻𝐼 ̂(𝑡 + ℎ) = 𝐺𝐻𝐼(𝑡) (4)
𝐺𝐻𝐼 ̂(𝑡) and 𝐺𝐻𝐼(𝑡) are respectively the predicted and measured hourly global horizontal solar irradiation at time t. This type of model is very simple to use but it gives results with very low accuracy.

In order to improve this model, the daily profile of the solar radiation (by clear sky) is taken into account and added in using the Solis clear sky model previously presented, then Eq. ( 4) becomes:

𝐺𝐻𝐼 ̂(𝑡 + ℎ) = 𝐺𝐻𝐼(𝑡). 𝐶𝑆(𝑡+ℎ) 𝐶𝑆(𝑡) (5)
This model is named scaled persistence (SP), a simple improvement of the persistence model [START_REF] Voyant | Statistical parameters as a means to a priori assess the accuracy of solar forecasting models[END_REF]. The advantages of naïve models is the fact that no historical data are required what makes these models easily usable in all cases and for all sites even if no solar radiation measurements exits for a long time. The accuracy of such models decreases rapidly with the time horizon and is generally not adequate for a more than one hour ahead.

Machine Learning models

A. Auto Regressive Mobile Average (ARMA)

The ARMA model [START_REF] Gooijer | 25 years of time series forecasting[END_REF] is a predictive model with two separate parts, the auto regressive one and the mobile average one. This model is able to predict the future values of time serie from a combination of past values of the time series and a white noise:

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = 𝜀 𝑡 + ∑ 𝜑 𝑖 𝑝 𝑖=0 . 𝐶𝑆𝐼(𝑡 -𝑖) + ∑ 𝜃 𝑖 . 𝜀(𝑡 -𝑖) 𝑞 𝑖=0 (6)
CSI (t+h) the clear sky index time series at time t+h, φ and θ the parameters of the autoregressive and moving average part deduced by a least square method, p and q are the orders of the model and ε(t) is an error term distributed like a Gaussian white noise for the time t.

B. Multi-Layer Perceptron (MLP)

The MLP is a type of artificial neural networks, here, a feed forward MLP with two layers (one hidden layer and one output layer) is used with, in input, the solar irradiation time series [START_REF] Kalogirou | Applications of artificial neural-networks for energy systems[END_REF]. The input matrix size was defined according to the results of the auto mutual method applied to the solar data.

The mathematical formula for a MLP with one hidden layer of m neurons, one output neuron and p input variables is a function described by:

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = ∑ 𝜔 𝑗,𝑠 𝑚 𝑗=1 . (𝑔. (∑ 𝜔 𝑖,𝑗 . 𝐶𝑆𝐼(𝑡 -𝑗) 𝑛-1 𝑖=0 + 𝑏 𝑗 )
+𝑏 𝑠 [START_REF] Voyant | Forecasting method for global radiation time series without training phase: Comparison with other well-known prediction methodologies[END_REF] with CSI the input vector of clear sky index constituted by the values of n variables, 𝐶𝑆𝐼 ̂(𝑡 + ℎ) the output vector with the predicted values of the model, bj the bias of the hidden neuron j and ω(i,j) the weight affected at the measured input CSI(t). g is the transfer function of the hidden neurons, bs the bias of the output neuron and ω(j,s) his weight affected to the output of the hidden neuron j. The optimization of the MLP is made by a classical technique: several configurations with a different number of hidden nodes are tested (with a number of hidden nodes varying between 3 and n+2) and the most efficient is selected.

C. Gaussian Process (GP)

The Gaussian process (GP) is a nonlinear modeling. It's a generalized Gaussian distribution with an infinity of variables [START_REF] Rasmussen | Gaussian Processes in Machine Learning[END_REF][START_REF] Salcedo-Sanz | Prediction of Daily Global Solar Irradiation Using Temporal Gaussian Processes[END_REF]

D. Support Vector Regression (SVR)

The Support vector regression (SVR) is a Kernel based model (as Gaussian process), it was firstly developed for regression problems and the application to time series forecasting is another successful evolution [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. Equation 10gives the formulation of the SVR applied on forecasting problems, with a training dataset 𝐷 = {𝐶𝑆𝐼(), 𝐶𝑆𝐼(𝑡 + ℎ)} with the test vector CSI(t+h) and h the prediction horizon:

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = ∑ 𝛼 𝑖 𝑡-1 𝜏=1
. 𝑘 𝑟𝑏𝑓 (𝐶𝑆𝐼(𝑡 + ℎ), 𝐶𝑆𝐼(𝑡 -𝜏) ) [START_REF] Breiman | Classification and regression trees[END_REF] and the Kernel radial basis function defined by:

𝑘 𝑟𝑏𝑓 (𝐶𝑆𝐼(𝑡 𝑝 ), 𝐶𝑆𝐼(𝑡 𝑞 ) ) = 𝑒𝑥𝑝 [ -(𝐶𝑆𝐼(𝑡 𝑝 )-𝐶𝑆𝐼(𝑡 𝑞 )) 2 2𝜎 𝑓 2 ] (11)
αi is the difference between Lagrange multipliers which are the solutions of a quadratic problem, b is the bias determined by specific conditions and Eq. ( 10).

E. Regression Trees

Decision trees based on "If-Then" rules are one of the most popular methods used in machine learning for classification, since they offer results that can be easily interpreted. Thus, this approach obtains ordered sets of If-Then rules for prediction that produces understandable models. In the last years, learning algorithms generating decision trees for classification problems have been extended for predicting values of attributes when they are numeric. These extensions have led to regression trees. A Regression Tree (RT) is a decision tree in which the leaf nodes have been set as regression models, and therefore, continuous numeric values can be predicted [35].

Standard and pruned regression trees (RT and RTpruned)

A model of regression trees is directly derived from the classification trees, Hastie and Tibshirani [START_REF] Hastie | Generalized additive models[END_REF] propose a formalization of the models:

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = ∑ 𝑘 𝑖 × 𝐼(𝐶𝑆𝐼(𝑡 -𝑖)) 𝑡-1 𝑖=1 (12)
with ki constant factors, I is a function which return 1 if input is used and 0 if it's not used; The trees are built by splitting the data based on the values of predictive attributes. Once the tree has been constructed, a regression model is computed for each node. The error for each node is the mean of the absolute difference between the predicted and the actual value of each instance of the training set that reaches the node. Pruning aspect of regression tree is operated with an elevation of the quadratic error tolerance per node. Splitting nodes stops when quadratic error per node drops below tolerance. For normal RT, the tolerance is close to zero, while for the pruned RT, a higher value is chosen using a heuristic method based on the minimizing of the global error of prediction. In pruned RT, the function I has more 0 than in the normal mode.

Boosted and bagged regression trees (RT-boosted and RT-bagged)

For RT, two types are boosted [START_REF] Breiman | Bagging Predictors[END_REF] and bagged of classification trees [START_REF] Friedman | Greedy function approximation: A gradient boosting machine[END_REF]. In boosting, successive trees give extra weight to points incorrectly predicted by earlier predictors. In the end, a weighted vote is taken for prediction. In bagging, successive trees do not depend on earlier trees each is independently constructed using a bootstrap sample of the data set. In the end, a simple majority vote is taken for prediction. The boosting method consists in assembly weak RT classifiers and take the average of predictions in order to improve the efficiency. In this case a weak predictor is a simple single split regression tree. The main technique is that the next trees give more weights to the data which have a bad prediction at the precedent point, at the end of the simulation, a weighted vote is taken for the prediction. The Bagging method is another version of the prediction with regression trees, it was described by Breiman [START_REF] Breiman | Bagging predictors[END_REF]. Bagging mean bootstrap aggregating, the model is an aggregation of regression trees which grown from samples of dataset. The subtrees are employed for the prediction and a vote is taken for the prediction.

Random Forest

Breiman [START_REF] Breiman | Random forests[END_REF] proposed the random forests, which add an additional layer of randomness to bagging. In addition to a construction for each tree using a different bootstrap sample of data, random forests change how the regression trees are constructed [START_REF] Breiman | Random forests[END_REF]. In standard trees, each node is split using the best split among all variables. In a random forest, the dataset is equally divided in samples but each regression tree grows differently, each node is split using the best among a subset of predictors randomly chosen at that node [START_REF] Breiman | Random forests[END_REF]. This improvement by randomness give robustness to the model and decreases the over-training risks.

ACCURACY PARAMETERS

In order to evaluate the accuracy of the different models, three different statistics are used. The mean absolute error (MAE), defines the absolute value of the gap between the observed and the predicted value:

𝑀𝐴𝐸 = ∑ |𝐺𝐻𝐼 ̂(𝑖)-𝐺𝐻𝐼(𝑖)| 𝑁 𝑖=1 𝑁 (13)
N is the number of variables.

The normalized root mean squared error (nRMSE) is more sensitive to big forecast errors, and hence is suitable for applications where small errors are more tolerable and larger errors cause disproportionately high costs, as for example in the case of utility applications. It is probably the reliability factor that is most appreciated and used: it's a good statistical index to evaluate the accuracy of a model, the aim of an operator is to minimize it in order to improve model performances:

𝑛𝑅𝑀𝑆𝐸 = √ 1 𝑁 ∑ (𝐺𝐻𝐼 ̂(𝑖)-𝐺𝐻𝐼(𝑖)) 2 𝑁 𝑖=1
𝐺𝐻𝐼 ̅̅̅̅̅̅ [START_REF] Smola | A tutorial on support vector regression[END_REF] where 𝐺𝐻𝐼 ̂ is the algebraic average of the observed values. The skill score is an index calculated in order to compare the performance of a given model with a reference model, here the reference model is the smart persistence model (SP) described in 4.1.).

𝑆𝑘𝑖𝑙𝑙 𝑆𝑐𝑜𝑟𝑒 = 𝑀𝑒𝑡𝑟𝑖𝑐 𝑓𝑜𝑟𝑒𝑐𝑎𝑠𝑡 -𝑀𝑒𝑡𝑟𝑖𝑐 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 𝑀𝑒𝑡𝑟𝑖𝑐 𝑝𝑒𝑟𝑓𝑒𝑐𝑡 𝑓𝑜𝑟𝑒𝑐𝑎𝑠𝑡 -𝑀𝑒𝑡𝑟𝑖𝑐 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 = 1 - 𝑛𝑅𝑀𝑆𝐸 𝑓𝑜𝑟𝑒𝑐𝑎𝑠𝑡 𝑛𝑅𝑀𝑆𝐸 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 (15)
The skill score is always inferior at 1, negative if the forecaster is less performant than the reference, 0 if the performances are similar, and positive if it is better.

FORECASTING MODELS PERFORMANCE

Tables 1 and2 show the results of the nRMSE and MAE calculations for the eleven models and a forecast horizon from 1 to 6 hours with an hourly resolution. The naïve models (persistence and smart persistence) give always bad performances in all the sites and is the worst model.

For Ajaccio with a low variability dataset the better models are ARMA and MLP, followed by scaled persistence and bagged regression trees; however the calculated error is similar for all the models excepted for the persistence and the simple regression tree. With such a variability, the machine learning models allow a good forecasting. Concerning the MAE, excepted for the naïve model, all the models are similar in term of absolute error, the results are in the same magnitude.

The skill score was computed in the same conditions for all the models, Figure 2 is the graphical representation of the skill scores (related to the smart persistence reference model).

CONCLUSION

Several statistical and machine learning tools (benchmark related to 11 forecasting models) were analyzed and compared in view to forecast the horizontal solar irradiation in Ajaccio. For a weak variability as Ajaccio, the two classical methods: ARMA and MLP seem to be the most efficient.

In a future work, we will integer solar data from other sites in order to improve the conclusions and maybe to modify the number of variability classes using clustering methods.

Figure 1 :

 1 Figure 1: Position of the meteorological station In this station, the period of available solar data are from 01/01/1998 to 12/31/2000 with 26280 validated data after the quality check.After the calculation of the auto-mutual information between the solar radiation data[START_REF] Jiang | Mutual information algorithms[END_REF][START_REF] Parviz | Mutual Information Based Input Variable Selection Algorithm and Wavelet Neural Network for Time Series Prediction[END_REF], it has been decided that the number of input data n in the forecasting model will be n=8.Thus, for estimating the solar irradiation at time t, the solar irradiations at time (t-1), (t-2), (t-3)… (t-n) will be used.

Figure 2 :

 2 Figure 2: Skill score vs forecast horizon for Ajaccio

Table 1 : nRMSE vs forecast horizon, the two best models are highlighted (results in %) Horizon h+1 h+2 h+3 h+4 h+5 h+6

 1

	Persistence 26.60 42.62 54.10 61.39 64.51 63.86
	SP 19.26 26.46 31.18 34.15 36.92 38.93
	ARMA 18.35 29.27 31.38 32.25 33.18 33.69
	MLP 18.26 29.26 31.31 32.47 32.98 33.84
	RT 24.64 36.88 38.47 39.74 39.95 41.24
	Boosted RT 18.75 29.55 31.89 32.51 33.55 33.98
	Bagged RT 18.76 29.80 31.10 32.17 33.35 34.02
	Pruned RT 18.72 30.88 32.27 33.76 34.01 35.00
	RF 18.97 29.63 31.62 32.38 33.37 33.91
	GP 18.97 30.08 31.96 33.29 33.55 34.44
	SVR 18.55 38.78 41.03 41.56 41.66 41.60

Table 2 : MAE values (in Wh/m²) vs forecast horizon, the two best models are highlighted Horizon h+1 h+2 h+3 h+4 h+5 h+6

 2

	Persistence 104.6 176.0 220.2 252.2 259.2 256.0
	SP 55.5 79.5 96.0 107.1 116.4 124.0
	ARMA 60.7 88.8 95.0 97.4 100.6 102.3
	MLP 60.6 89.4 95.2 98.9 100.7 103.3
	RT 73.3 117.7 125.0 129.8 132.0 135.6
	Boosted RT 61.4 91.8 98.1 100.1 103.3 104.6
	Bagged RT 61.4 93.3 97.3 101.3 104.6 106.0
	Pruned RT 60.8 95.4 100.0 105.0 106.1 108.9
	RF 61.1 92.9 99.1 102.0 105.0 106.5
	GP 61.9 96.1 101.1 104.8 106.3 108.3
	SVR 54.6 126.6 133.3 135.5 135.1 134.6


Bahloul Rachid

Bahaj Mohamed

Sidki Omar

Periodic solutions of differential equations with three variable in vector-valued space

Keywords: Differential equations, L p -multipliers, vector-valued function space

The aim of this work is to study the existence and uniqueness of a periodic solutions of differential partial functional differential equations with three variable:

Our approach is based on the L p -multipliers of linear operators.

Introduction

Motivated by the fact that neutral functional integro-differential equations (abbreviated, NFDE) with finite delay arise in many areas of applied mathematics, this type of equations has received much attention in recent years. In particular, the problem of existence of periodic solutions, has been considered by several authors. In this work, we study the existence and uniqueness of periodic solutions for the following integro-differential equation with two variable of the following form ∂ ∂x ∂ ∂y ∂ ∂w u(x, y, w) + Au(x, y, w) = f (x, y, w) for x, y, w ∈ R,

where (A, D(A)) is a closed linear operator on a Banach spase X, and f : R 3 → X is a locally p-integrable and 2π-periodic function for 1 ≤ p < ∞.

In [START_REF] Arend | The operator-valued Marcinkiewicz multiplier theorem and maximal regularity[END_REF] and [START_REF] Arendt | Operator-valued Fourier multipliers on periodic Besov spaces and applications[END_REF], Arendt gave necessary and sufficient conditions for the existence of periodic solutions of the following evolution equation in vecteur valued and Besov space.

In [START_REF] Bahloul | Periodic solutions of degenerate equations with finite delay in UMD space[END_REF], Bahloul et al established results on the existence for some degenerate differential equation.

1 : E-mail address : bahloul33r@hotmail.com 2: E-mail address : mohamedbahaj@gmail.com 1: E-mail address : osidki@hotmail.com

In [START_REF] Bahloul | Periodic solutions of differential equations with two variable in vector-valued function space[END_REF], Bahloul established results on the existence of periodic solution of Eq. (1.1) when a = 0, namely, for the following partial functional differential equation

∂ ∂x ∂ ∂y u(x, y) + Au(x, y) = f (x, y)
where (A, D(A)) is a linear operator. Similar characterizations were then obtained for the scale of degenerate integro-differential and second differential equations [START_REF] Bu | Mild well-posedness of second order differential equations on the real line[END_REF], [START_REF] Bu | Solutions of second order degenerate integro-differential equations in vector-valued function spaces[END_REF] and [START_REF] Bu | L p -Maximal Regularity of Degenerate delay Equations with Periodic Conditions[END_REF].

In [START_REF] Hernan | Prokopczyk Periodic Solutions of abstract neutral functional differential equations[END_REF], Hernan et al, studied the existence of periodic solution for the class of linear abstract neutral functional differential equation described in the following form:

d dt [x(t) -Bx(t -r)] = Ax(t) + G(x t ) + f (t) for t ∈ R
where A : D(A) → X and B : D(B) → X are closed linear operator such that D(A) ⊂ D(B) and G ∈ B(L p ([-2π, 0], X); X). We refer the readers to papers [START_REF] Favini | Degenerate Differential Equations in Banach Spaces[END_REF], [START_REF] Keyantuo | Fourier multipliers and integro-differential equations in Banach spaces[END_REF], [START_REF] Keyantuo | Maximal regularity for a class of integro-differential equations with infinite delay in Banach spaces[END_REF] and [START_REF] Lindenstrauss | Classical Banach Spaces II[END_REF] and the references listed therein for information on this subject.

In [START_REF] Lizama | Fourier multipliers and periodic solutions of delay equations in Banach spaces[END_REF], Lizama gave necessary and sufficient conditions for the existence of periodic solutions of the following equation.

d dt x(t) = Ax(t) + L(x t ) + f (t) for t ∈ R,
where A is a closed linear operator on an UMD-space Y . We note that, starting with the work [START_REF] Lizama | Fourier multipliers and periodic solutions of delay equations in Banach spaces[END_REF], the problem of characterization of maximal regularity for avolution equations with periodic initial conditions have been studied intensively in the last years. See e.g [START_REF] Lizama | Periodic solutions of fractional differential equations with delay[END_REF], [START_REF] Lizama | Maximal regularity of delay equations in Banach spaces[END_REF] and [START_REF] Hino | Almost periodic solution of Differential Equations in Banach Spaces[END_REF], see alsoo [START_REF] Hale | Theory of Functional Differential Equations[END_REF] and [START_REF] Poblete | Maximal regularity of second-order equations with delay[END_REF] and references therien. the main novelty in this paper relies in the presence of two variable. This work can be considered as a progress on the treatment of such kind of problems. This work is organized as follows : After a creterion in the second section, we are able to characterize in section 3 the existence and uniqueness of strong L p -solution for the problem (1.1) solely in terms of a property of L p -multiplier for the sequence of operators -ikzl(-iklzI + A) -1 . We optain that the following assertion are equivalent in general Banach spaces:

(1)

(-iklzI + A) is invertible and {-iklz(-iklzI + A) -1 , (k, l, k) ∈ Z 3 } is an L p -muliplier.
(2) For every f ∈ L p (T 3 ; X) there exist a unique function u ∈ H 1,p (T 3 ; X) such that u ∈ D(A) and equation (1.1) holds for a.e (x, y, w) ∈ [0, 2π] 3 .

In section 4, we give the main abstract result (theorem (4.2)) of this work.

In section 5, we propose an application to some integro-differential equation. In section 6, we give the conclusion.

A creterion for periodic solutions

Let X be a Banach Space. Firstly, we denote By T the group defined as the quotient R/2πZ. There is an identification between functions on T and 2π-periodic functions on R.

We consider the interval [0; 2π) as a model for T. Given 1 ≤ p < ∞, we denote by L p (T 3 ; X) the space of 2π-periodic locally p-integrable functions from R 3 into X, with the norm:

f p := ( 2π 0 2π 0 2π 0 f (x, y, w) p dxdydw) 1/p
For f ∈ L p (T 3 ; X), we denote by f (k, l, z), (k, l, z) ∈ Z the (k, l, z)-th Fourier coefficient of f that is defined by:

f 3 (k, l, z) = ( 1 2π ) 3 2π 0 2π 0 2π 0 e -i(kx+ly+zw) f (x, y, w)dxdydw for (k, l, z) ∈ Z × Z and (x, y, w) ∈ R × R. Notation 2.1. : Let (k, l, z) ∈ Z × Z and 1 ≤ p < ∞. Denote by H 1,p (T 3 ; X) = u ∈ L p (T 3 ; X) : ∃v ∈ L p (T 3 ; X), v 3 (k, z) = -iklz u 3 (k, z) Remark 2.1. [9] Let f ∈ L 1 (T; X). If g(t) = t 0 f (s)ds and k ∈ Z, k = 0, then ĝ(k) = i k f (0) -i k f (k) such that: f (k) = 1 2π 2π 0 e -ikx dx. Lemma 2.2. [4] Let f ∈ L 1 (T × T; X). If w(x, t) = x 0 t 0 f (s, ξ)dξds and (k, z) ∈ Z × Z -{(0, 0)} then w(k, z) = -1 kz f (0, 0) + 1 kz f (0, z) + 1 kz f (k, 0) -1 kz f (k, z) such that: f (k, z) = ( 1 2π ) 2 2π 0 2π 0 e -i(kx+zy) dxdy. Lemma 2.3. [4] Let 1 ≤ p < ∞ and let f, g ∈ L 1 (T × T; X).
The following are equivalent:

(1) v 3 (k, l, z) = 0 if k = 0, l = 0 or k = 0 and there exists a ∈ X such that u(x, y, w) = a + x 0 y 0 w 0 v(s, τ, ξ)dξdτ ds, (x, y, w) ∈ T 3 (2) v 3 (k, l, z) = -iklz u 3 (k, l, z), (k, l, z) ∈ Z 3 Proof. Its proof is inspired by that of [4, Lemma 2.3].
Remark 2.2. By Lemma (4.1) we have u ∈ H 1,p (T 3 ; X) ⇔ there exists v ∈ L p (T 3 ; X) such that v 3 (k, l, z) = 0 if k = 0, l = 0 or k = 0 and there exists a ∈ X such that u(x, y, w) = a +

x 0 y 0 w 0 v(s, τ, ξ)dξdτ ds, (x, y, w) ∈ T 3 . Now we can formulate the following multiplier definitions. Definition 2.1. For 1 ≤ p < ∞, a sequence {M k,l,z } (k,l,z)∈Z 3 is said to be an L p -multiplier (or (L p , L p )-multiplier) if for each f ∈ L p (T 3 , X), there exists u ∈ L p (T 3 , Y ) such that

g 3 (k, l, z) = M k,l,z f 3 (k, l, z) for all (k, l, z) ∈ Z 3 . Definition 2.2. : For 1 ≤ p < ∞ , a sequence {M k,l,z } (k,l,z)∈Z×Z is said to be an (L p , H 1,p )- multiplier if for each f ∈ L p (T 3 , X), there exists u ∈ H 1,p (T 3 , Y ) such that g 3 (k, l, z) = M k,l,z f 3 (k, l, z) for all (k, l, z) ∈ Z 3 .
We recall the following results Lemma 2.4. [START_REF] Bahloul | Periodic solutions of differential equations with two variable in vector-valued function space[END_REF] Let f, g ∈ L p ([0, 2π] 3 ; X), where 1 ≤ p < ∞. Then the following are equivalent.

(i) f (x, y, w) ∈ D(A) and Af (x, y, w) = g(x, y, w),

(ii) f 3 (k, l, z) ∈ D(A) and A f 3 (k, l, z) = g 3 (k, l, z) for all (k, l, z) ∈ Z 3 3 A caracterization of stong L p -solution Definition 3.1. : Let f ∈ L p (T 3 ; X). A function u ∈ H 1,p (T 3 ; X)
is said to be a 2π-periodic strong L p -solution of Eq. (1.1) if u ∈ D(A) and Eq.(1.1) holds almost every where. Proposition 3.1. Let X be a Banach space. Suppose that for every f ∈ L p (T 3 ; X) there exists a unique strong solution of Eq. (1.1) for 1 ≤ p < ∞. Then (1) for every (k, l, z) ∈ Z 3 the operator (-iklzI + A) has invertible.

(2) -iklz(-iklzI + A) -1 (k,z)∈Z×Z is an L p -multiplier.

Proof. 1. Let k, l, z ∈ Z and v ∈ X. Then for f (x, y, w) = e i(kx+ly+zw) v , there exists u ∈ H 1,p (T 3 ; X) such that:

∂ ∂x ∂ ∂y ∂ ∂w u(x, y, w) + Au(x, y, w) = +f (x, y, w)
Taking Fourier transform. Then we obtain kx+ly+zw) v is a 2π-periodic strong L p -solution of Eq. (1.1) corresponding to the function f = 0 Hence u(x, y, w) = 0 and v = 0 then (-iklzI + A) is injective.

(-iklzI + A) u 3 (k, l, z) = f 3 (k, l, z) = v ⇒ (-iklzI + A) is surjective. if (-iklzI + A)(v) = 0, then, u(x, y, w) = e i(
2) Let f ∈ L p (T 3 ; X). By hypothesis, there exists a unique u ∈ H 1,p (T 3 , X) such that the Eq.(1.1) is valid. Taking Fourier transforms, we deduce that

u 3 (k, l, z) = (-iklzI + A) -1 f 3 (k, l, z) for all k, l, z ∈ Z, then (-iklzI + A) -1 (k,l,z)∈Z 3 is an L p -multiplier. Hence -iklz u 3 (k, l, z) = -iklz(-iklzI + A) -1 f 3 (k, l, z) for all k, l, z ∈ Z. On the other hand, since u ∈ H 1,p (T 3 ; X) , there exists w ∈ L p (T 3 ; X) such that w 3 (k, l, z) = -iklz u 3 (k, l, z) = -iklz(-iklzI + A) -1 f 3 (k, l, z) i.e -iklz(-iklzI + A) -1 (k,l,z)∈Z 3 is an L p -multiplier.

Main result

The next Theorem is the main result of this section.

Lemma 4.1. Let f ∈ L 1 (T 3 ; X). If w(x, y, w) =

x 0 y 0 w 0 f (s, τ, ξ)dξdτ ds and (k, z) ∈ Z 3 -{(0, 0, 0)} then

w 3 (k, l, z) = i 3 klz [ f 3 (0, 0, 0)-f 3 (k, 0, 0)-f 3 (0, l, 0)-f 3 (0, 0, z)+ f 3 (k, l, 0)+ f 3 (0, l, z)+ f 3 (k, 0, z)-f 3 (k, l, z)]
Proof.

w 3 (k, l, z) = i z [ x 0 y 0 f (s, τ, 0)dτ ds] - i z [ x 0 y 0 f (s, τ, z)dτ ds] = i z [ i l x 0 f (s, 0, 0)ds - i l x 0 f (s, l, 0)ds] - i z [ i l x 0 f (s, 0, z)ds - i l x 0 f (s, l, z)ds] = i 3 zl [ x 0 f (s, 0, 0)ds - x 0 f (s, l, 0)ds - x 0 f (s, 0, z)ds + x 0 f (s, l, z)ds] = i 3 zl [( i k f 3 (0, 0, 0) - i k f 3 (k, 0, 0)) -( i k f 3 (0, l, 0) - i k f 3 (k, l, 0)) -( i k f 3 (0, 0, z) - i k f 3 (k, 0, z)) + ( i k f 3 (0, l, 0) - i k f 3 (k, l, 0))] = i 3 klz [ f 3 (0, 0, 0) -f 3 (k, 0, 0) -f 3 (0, l, 0) -f 3 (0, 0, z) + f 3 (k, l, 0) + f 3 (0, l, z) + f 3 (k, 0, z) -f 3 (k, l, z)].
Theorem 4.2. Let A be a closed linear operator and 1 ≤ p < ∞. If

(1) for every (k, l, z) ∈ Z 3 the operator M k,l,z = (-iklzI + A) has invertible.

(2) -iklz(-iklzI + A) -1 (k,l,z)∈Z 3 is an L p -multiplier. Then for every f ∈ L p (T 3 ; X) there exists a unique 2π-periodic strong L p -solution of Eq. (1.1). Lemma 4.3. [4] Let 1 ≤ p < ∞ and X be a Banach space. The following assertions are equivalent

(1) {M k,l,z } (k,l,z)∈Z 3 is an (L p , H p )-multiplier. (2) {-iklzM k,l,z } (k,l,z)∈Z 3 is an L p -multiplier. Proof. 1 ⇒ 2. Let f ∈ L p (T 3 , X), then there exists g ∈ H 1,p (T 3 ; X) ( in particular g ∈ L p (T 3 ; X)) such that g 3 (k, l, z) = -iklzM k,l,z f 3 (k, l, z), then {-iklzM k,l,z } (k,l,z)∈Z 3 is an L p -multiplier. 2 ⇒ 1. Let f ∈ L p (T 3 ; X), by hypothesis ∃v ∈ L p (T 3 ; X) such that v 3 (k, l, z) = -iklzM k,l,z f 3 (k, l, z).
In particular v 3 (k, l, z) = 0 if k = 0, l = 0 or z = 0. Let w(x, y, w) =

w 3 (k, l, z) = i 3 klz [ v 3 (0, 0, 0)-v 3 (k, 0, 0)-v 3 (0, l, 0)-v 3 (0, 0, z)+ v 3 (k, l, 0)+ v 3 (0, l, z)+ v 3 (k, 0, z)-v 3 (k, l, z)] w 3 (k, l, z) = -1 iklz v 3 (k, l, z) = M k,l,z f 3 (k, l, z
), for all k = 0, l = 0 and z = 0 Let u(x, y, z) = w(x, y, z) + M 0,0,0 f 3 (0, 0, 0) -w 3 (0, 0, 0) i.e u(x, y, z) = a+ x 0 y 0 w 0 v(s, τ, ξ)dξdτ ds, a = M 0,0,0 f 3 (0, 0, 0)-w 3 (0, 0, 0). Then by Remark (2.1) we have u ∈ H 1,p (T 3 ; X).

u 3 (k, l, z) = w 3 (k, l, z) + 0 = w 3 (k, l, z) = M k,l,z w 3 (k, l, z), for all (k, l, z) ∈ Z 3 , Proof of theorem (4.2). Existence: Let f ∈ L p (T 3 ; X). Define M k,z = (-iklzI + A) -1 by theorem (4.2), the family -iklz(-iklzI + A) -1
(k,l,z)∈Z 3 is an L p -multiplier and by Lemma (4.3), it is equivalent to the fact that the family {M k,l,z } (k,l,z)∈Z 3 is an (L p , H p )multiplier, i.e, there exists u ∈ H 1,p (T 3 ; X) such that

u 3 (k, l, z) = M k,l,z f 3 (k, l, z) = (-iklzI + A) -1 f 3 (k, l, z) (4.1)
In particuler, u ∈ L p (T 3 , X) and there exists v ∈ L p (T 3 ; X) such that For all (k, l, z) ∈ Z 3 . Then using that A is closed, we conclude that u ∈ D(A) and ∂ ∂x ∂ ∂y ∂ ∂w u(x, y, w) + Au(x, y, w) = f (x, y, w). Uniqueness : Suppose that the Eq. (1.1) has two strong L p -solution u 1 and u 2 , then u = u 1 -u 2 is a solution of Eq. (1.1) corresponding to the fonction f = 0. Taking Fourier transform on both sides of this equality we get (-iklzI + A) u 3 (k, l, z), which implies that u 3 (k, l, z) = 0 for all (k, l, z) ∈ Z 3 , consequently, u = 0. Then u 1 = u 2 .

v 3 (k, l, z) = -iklz u 3 (k, l, z) (4.2) Moreover, ∂ ∂x ∂ ∂y

Conclusion

For every function f ∈ L p (T 3 ; X), this paper deals with the existence and uniqueness of periodic solutions of the differential equation with three variable ∂ ∂x ∂ ∂y ∂ ∂w u(x, y, w) + Au(x, y, w) = f (x, y, w). Using the method of Fourier transform of equation (1.1) and L p -multipliers of linear operators. the results obtained in this paper generalizes the result of Bahloul [START_REF] Bahloul | Periodic solutions of differential equations with two variable in vector-valued function space[END_REF].

  , τ, ξ)dξdτ ds, then by Lemma (4.1) we have,

∂ ∂w u 3 (

 3 k, l, z) = -iklz v 3 (k, l, z). By Fejer's theorem one has in L p (T 3 ; X) u 3 (k, l, z) = lim n→∞ lim n →∞ lim m→∞ σ n,n ,m (u)(x, y, w) such that σ n,n ,m (u)(x, y, w) e i(kx+ly+zw) u 3 (k, l, z). By (4.1) and (4.2) we have -iklz u 3 (k, l, z) + A u 3 (k, l, z) = f 3 (k, l, z)


M P Pai'doussis

G X Li

CROSS-FLOW-INDUCED CHAOTIC VIBRATIONS OF HEAT-EXCHANGER TUBES IMPACTING ON LOOSE SUPPORTS

The chaotic dynamics of heat exchanger tubes impacting on the generally loose baffle plates is studied using an analytic model that involves delay differential equations. The critical flow velocity for the local instability of the flexible cylinder near the static equilibrium position is obtained by assuming a harmonic solution in the discretized linearized model and solving the resulting algebraic equations. Numerical solutions show that, with increasing flow beyond the critical, the amplitude of motion grows until impacting with the loose support occurs; more complex motions then arise, leading to chaos for a sufficiently high flow velocity. A Lyapunov exponent technique is developed for delay differential equations and is applied to the present model, showing definitely and unequivocally-for the first time for this system-that chaotic motions do occur. To understand the system behaviour better, a one-degree-of-freedom, forced, negatively damped impact oscillator is studied. The response is quite complex and results in chaos in certain parameter regions. The analysis is performed by finding periodic solutions and determining their stability and bifurcations with the Poincare map technique.

INTRODUCTION

Arrays of cylinders subject to cross-flow are known to be subject to fluid elastic instabilities at sufficiently high flow velocities. In tube-in-shell type heat exchangers, steam generators, boilers, condensers and similar equipment, these instabilities can cause severe damage to the tubes: the resultant high-amplitude vibration may sometimes cause inter-tube clashing and always wear with the tube-supporting baffle plates [START_REF] Pa | iooussis 1980 Practical Experiences with Flow-induced Vibrations[END_REF], eventually leading to breakdowns.

Two mechanisms are generally accepted as being responsible for fluid elastic instability [START_REF] Chen | [END_REF][3][START_REF] Price | [END_REF][START_REF] Paidoussis | [END_REF]: a velocity controlled, negative damping mechanism, and a stiffness controlled one. The first mechanism gives rise to single-degree-of-freedom flutter; hence array stability may be analyzed, approximately, by considering all cylinders of the array except one to be rigid, and considering motions of just that one cylinder [6; 4, 5]. The second mechanism involves at least two degrees of freedom and gives rise to coupled-mode flutter [START_REF] Paidoussis | [END_REF]. To clarify matters further by reference to well known aeroelastic phenomena, the first mechanism is similar to Den Hartog's for the galloping of iced transmission lines; the second is akin to wake flutter of aeroplane wings.

For cylinder arrays, both mechanisms are generally at work, although the contribution of one of the two may in some cases be so small as to be negligible; thus, for instance, allowing in some circumstances the analysis of a fully flexible array of cylinders by considering only one to be flexible, as mentioned above. For some array geometries,t it has been found that for low values of the mass-damping parameter,:j: m8 1 pD 2 , say for m8 1 pD 2 < 300 approximately, the negative damping instability is dominant. This is the case, for instance, for rotated triangular arrays and square arrays. However, this is not always true. For rotated square arrays, even for low values of m8 I pD 2 , the instability is of the two-degreeof-freedom variety; moreover, this type of array is prone to be subject to a recently discovered static instability, i.e. fluid-elastic divergence [7].

There are three main analytic or quasi-analytic (semi-empirical) models for predicting oscillatory fluid-elastic instability [START_REF] Chen | [END_REF][START_REF] Price | [END_REF]6,8,9], as well as others. The interested reader is referred to Chen's [10] and Paidoussis' [11] discussions of alternative classifications of these models, as well as to Paidoussis and Price's [START_REF] Paidoussis | [END_REF] discussion of the mechanisms giving rise to the instability in simple terms. All these models are linearized, both in terms of the fluid-dynamic forces acting on the cylinders and in terms of structural motions. Of course, in reality, non-linear fluid and structural forces come into play when the oscillation amplitude becomes substantial, limiting motion into a stable limit cycle. The effect of fluid-force non-linearities on the dynamics of arrays of cylinders in cross-flow has been studied by Price and Valerio [START_REF] Price | Flow-induced Vibration-1989[END_REF]. However, in the particular case of heat-exchanger type cylinder arrays, the amplitude limitation may be controlled by another, physical, motion-limiting constraint, as will be explained below.

In heat exchangers, the tubular cylinders are threaded through the baffle plates (the purpose of which is to direct the flow and to support the cylinders) through holes, which are oversize to permit easy threading and to allow for thermal expansion. After some use, crud may be deposited on the baffle plates, rendering some, but not all, loose-fitting holes into positive supports; also, the fact that the tubes are not perfectly straight originally helps explain how in a real heat exchanger some baffle plates provide positive support to any given cylinder, while others do not. Now, after the onset of fluid elastic instability, impacting will occur with the baffle plate at the loose supports, which is known to give rise to structural non-linear forces [START_REF] Axisa | [END_REF][START_REF] Fricker | Proceedings ASME International Symposium on Flow-induced Vibration[END_REF]. This non-linearity is considered to perhaps be more important for heat exchanger arrays than the fluid-force non-linearities, since impacting occurs even at relatively small amplitude (the gap to the baffle plate being generally small), when fluid non-linearities are still relatively unimportant.

The present paper is concerned with the dynamics of heat exchanger type cylinders in cross-flow, ofthe kind (in terms of geometry and values of 8 and m/ pD 2 ) where the singledegree-of-freedom, negative damping instability is dominant. Hence, the analysis is carried out by considering a single flexible cylinder, with all the other cylinders in the array rigid and immobile; the analytic model used is that of Price and Paidoussis [START_REF] Price | [END_REF]. For the reasons given in the previous paragraph, the only non-linearities considered here are those associated with impacting on the baffle plates.

Although the proposed model is in many ways close to an actual heat-exchanger type array, some simplifications and idealizations have been introduced in order to enable the analysis to bring out the essential dynamics, without the encumbrances that more realistic modelling would entail. Thus, although tubular cylinders in most heat exchangers are multi-span, here a two-span model will be considered, with the extremes of the cylinder being positively supported (clamped) and the middle support being loose; see Figure l(a).

-r Here "geometry" is characterized by the layout of the cylinders (e.g., into square or triangular patterns) and by their orientation vis-a-vis the upstream flow velocity vector (thus, "normal triangular", or "rotated square" arrays).

t Here m is the cylinder mass per unit length, 8 is the logarithmic decrement of damping, the p is the fluid density and D is the cylinder diameter. Also, although in a real heat exchanger the baffie plates are circular and hence post-impact sliding could involve quasi-orbital motion of the flexible cylinders, t it will here be assumed for simplicity that motions will be planar and in the direction of the single degree of freedom considered, in the cross-stream direction.

The main motivation for this work is to show definitively, in quantitative terms, whether chaotic oscillations are possible in such systems. Axisa et al. [START_REF] Axisa | [END_REF] obtained chaotic-looking oscillations with their model of an impacting heat exchanger tube; however, since they did not conduct calculations which would provide a quantitative measure of whether the oscillations were truly chaotic, the question of the existence of chaos in such systems was not unequivocally answered in that work. Nevertheless, a most significant finding of that work [START_REF] Axisa | [END_REF] was that increased wear rates are obtained with the higher frequencies associated with such oscillations. Thus, this problem is of both fundamental and practical interest.

The dynamics of the system are studied through simulation, obtaining bifurcation and phase flow diagrams and also Lyapunov exponents. In order to understand the dynamics better and explain some of the phenomena observed in the numerical results, an even more simplified model is examined in the last part of the paper, involving a single-degree-offreedom, negatively damped impact oscillator under the excitation of harmonic forcing. Such a model allows the performance of some non-linear analysis, which helps toward understanding the underlying mechanism for the observed chaotic behaviour. and D is the diameter of the cylinder, p and U are the fluid density and velocity, respectively, Cn and CL are the drag and lift coefficients, respectively, based on the flow velocity in the gap between the cylinders, Cma is the virtual or "added" mass coefficient of the fluid around the cylinder, lit is the time delay, and w•..,(9(1); see also Pai'doussis and Price [START_REF] Paidoussis | [END_REF].

The origin of the delay lit is in the assumption made in this quasi-steady model that there is a delay (phase 1ag) between cylinder motions and the fluid-dynamic forces generated thereby, associated with re-adjustment of the viscous wake to the altered position of the cylinder vis-a-vis its neighbours.

It should be mentioned here that if another model leading to a negative-damping type of fluid-elastic instability had been utilized instead of that of references [START_REF] Price | [END_REF]8], qualitatively similar results would have been obtained; thus, the conclusions of this study are not specific to the fluid-elastic instability model utilized. None of the available mathematical models are "perfect"; for a discussion of strengths and weaknesses of the available models, the reader is referred to references [10,11].

By introducing the following non-dimensional quantities, in which the cylinder length is denoted by L, w X -2g

1 7J=jj• t;=z• t =A-1 mL 4 t=Qit, s = -. mnl p = 1/ ( 1 + n ~ma), (3)
m -2nU J-f m=-- U=- -~, pD2' DQI' m1
and then dropping the tilde over timet, equation (1), combined with (2), may be written in the following dimensionless form: where A-1 is the dimensionless eigenvalue of the first mode for a clamped beam, and r = 2nj0 for Jl =I.

1 o 2 11 ( OcD) o17 1 o 4 7J
t Supports in real heat exchangers are somewhere between a simple support and a clamped one, and exact integral relationships between frequencies are rather rare.

The system is now discretized by utilizing the standard Galerkin expansion N 11 <~. t)= I t/J;(~)q;(t), i=l [START_REF] Paidoussis | [END_REF] where the t/J;( ~)are the dimensionless orthonormal set of eigenfunctions for the clampedclamped beam and have the properties:

J 1 {1 for i=j, t/J;(~)t/Jj{~)dx= 0 CO '..J.'• 0 10f l-r],
and Upon substituting equation ( 5) into (4), multiplying by <Pi~) and integrating over the domain [0, 1], one obtains q;(t)+ -+ -4 _ q;(t)+v;q;(t)-- 8 2 _--;-q;(t-r)+./(1'/b)t/J;(~b)=O, (6) 1tc nm TC m uT/ where V;= (A.;/ A. 1 ) 2 are the ratios of the non-dimensional natural frequencies to the first one, 1Jh is the displacement at~= ~b. and the viscous damping term has been replaced by the modal damping 8;V;jtc, in which it has been assumed that damping coefficients are small.

Notice that the only coupling term in the above equation is the force f( 1'/b) due to the constraint at the intermediate loose support. If there were no such constraint, then the resulting equations would be decoupled.

In what follows, equation ( 6) will be studied by modelling the constraint force by a cubic spring, /(1'/) = ICT/ 3 , and by a trilinear spring, /(

1'/) = 1C[1J-!<111 +di-11J -di)],
respectively, where 2d represents the gap at the constraint, non-dimensionalized with respect to D. For the cubic spring representation of the constraint force, symmetric and asymmetric cases are studied separately ( ~b = 0• 50 and 0•45, respectively, in sections 4 and 6). The cubic spring representation is less physically realistic than the trilinear one; the reason for utilizing it at all is that it is analytic and so permits calculations (e.g., of the Lyapunov exponents, section 5) which are not possible otherwise. Some of the calculations are time-consuming and hence they were all conducted for a single set of physical parameters. The system chosen was a square array with pitch-todiameter ratio pjD=I•5, for whiCh CD=0•26 and oCdoT/=-8•1 have been measured [START_REF] Price | [END_REF].t To avoid complications arising from multiple stability boundaries [START_REF] Price | [END_REF], a not-too-low value of m and 8 were chosen, namely m= 3 and 8;=0•06 for all i. Note that m= 3 signifies that, realistically, the fluid is a liquid (e.g., water flowing around brass tubes); with Cma= I• 2 this leads to f3 = 0• 24.

Concerning the impact stiffness, An tunes [START_REF] Antunes | de r Universite Paris VI. Contribution a I' etude des vibrations de faisceaux de tubes en ecoulement transversal[END_REF] calculated a value of 6• 6 x 10 6 N /m, while Yetisir and Weaver [START_REF] Yetisir | Fluid-Structure Interaction and Aerodynamic Damping[END_REF] give values of 10 6 -10 7 Njm for impact with anti-vibration bars (A VB). Adopting the range of 10 6 -10 7 N /m, this gives a non-dimensional stiffness 1( = 10 3 -10 4

• Since numerical simulation showed that the dynamics of the system are qualitatively insensitive to 1C in this range, 1C = 10 3 is used throughout.

In summary, unless otherwise stated, the fixed parameters in this study are CD= 0• 26, oCdo1J = -8•1, m=3, 8;=0•06, f3 =0•24 and 1(= 10 3

• The main parameter to be varied is the dimensionless flow velocity, 0.

In what follows, all quantities are understood to be dimensionless, so the tildes on 0, l and m are henceforth dropped.

THE HOPF BIFURCATION

To find the critical flow velocity for the onset of instability, all system parameters are fixed as defined in the previous section, and the dimensionless flow velocity U is gradually increased. Eventually, for high enough U, there is a point, UH, at which the stable equilibrium (at the undeformed position) of the cylinder is lost. Thereafter, for any U> UH, the motion of the cylinder is amplified, until the non-linear elements of the system counterbalance the growth, bringing the system to a steady oscillatory state. A Hopf bifurcation is defined as the loss of the stable equilibrium and onset of amplified oscillation.

The critical value,UH, may be found by solving the linearized system of equations ( 6) by letting/( 1Jb) = 0. Upon substituting the parameter values given in section 2, one obtains (7) where a1 =0•0145v;, a2=0•00524, a 3 =0•76v1, a4=0•026; and V; and thus a; are known for each mode. The above equations are decoupled and thus can be treated independently.

For the ith mode, assuming a solution of the form in which j is the imaginary unit and m; is the ith dimensionless eigenfrequency, and substituting it into equation ( 7), gives

-m7 + (aJ + a2U)jm;+ (a 3 + a 4 U 2 e-im,~) =0.
At the threshold of the Hopfbifurcation, U= UH, m; is purely real, and this equation leads to (8) where r = 2n-I U H. For fixed parameter values, U H and m; are the only unknowns for each mode, and are solved for numerically; their values are listed in Table 1.

Of all the critical flow velocities, the first one, U HI = 1• 785, is the most important, since thereafter an oscillatory state prevails for the double-span system of Figure 1. However, higher velocities are also of interest ( U> U HI), since in the design of such systems the value of UH is normally computed assuming a single span (i.e., of length L/2 in Figure 1). Referring to the definition of the dimension1ess U in equation ( 3), it is seen that, based on L/2, the "design" dimensional critical flow velocity is higher than even that of the second mode Hopf bifurcation found here, because the middle support is loose. If three or more spans had been considered, and they happened to be loose also, one can see that flow velocities many times larger than U HI are of interest. Of course, linear analysis, yielding the Hopf bifurcation, cannot give any information about limit-cycle motion; moreover, in the linear system the modes are uncoupled. Nevertheless, it is instructive to discuss here how limit-cycle motion can develop with the model in this paper. Once U exceeds U 8 , amplified motions lead to impact with the loose support, whereupon the modes do become coupled. Thereafter, the energy flows into the system (from the flowing fluid) in the mode or modes for which U> UH, while the higher modes that are still stable provide damping. Thus, a balance between the negatively damped unstable modes, dumping energy into the system, and the positively damped stable ones, dissipating it, leads to a state of dynamic equilibrium, i.e., to a limit cycle.

BIFURCATION DIAGRAM FOR THE SYMMETRIC CASE

Numerical simulations of the dynamics of certain types of cross-flow systems have been reported by many researchers (see Axisa et al. [START_REF] Axisa | [END_REF] and Fricker [START_REF] Fricker | Proceedings ASME International Symposium on Flow-induced Vibration[END_REF], for example). But all of them were presented in the form of either phase flows or time histories for certain individual values of the flow velocity, with other parameters fixed. With the fast computers now available, it is easy and helpful to display the system response for some continuous ranges of the flow velocity in order to obtain a full understanding of system dynamics; this is also helpful in deciding more definitely if the system becomes chaotic or not. Such a diagram also allows one to observe and pinpoint any sudden changes (bifurcations) in the system behaviour. The specific choice of the output signal to be processed in this way depends on the individual system, and usually the one which can provide most information is used.

For the present system, the simulation was carried out by integrating equation ( 6) using the standard fourth order Runge-Kutta algorithm. The parameter varied is U, and the output utilized to display bifurcations is the approximate dimensionless displacement at the middle point of the cylinder, q(O• 5, t), obtained through equation [START_REF] Paidoussis | [END_REF]. The triggering signal is oq(O• 5, t) I ot: when it is equal to zero and q(O• 5, t) > 0, then the value of q(O• 5, t) is recorded. In the bifurcation diagrams to be presented, solutions for the first 200 time units were discarded, so as to achieve steady state solutions, and then q(O• 5, t) was recorded for the next 50 time units.

Before proceeding with the computation, it is important first to decide on the number of modes necessary in the simulation to generate the essential dynamics of the physical system. For the cubic constraint model (symmetric when ~b = 0• 5 and asymmetric when ~b = 0•45) and the trilinear one, phase plots were generated for values of U before and after the critical value for instability. Study of these phase flows indicates that simulation results from models with five modes and higher are qualitatively the same both in magnitude and in shape when the integration step size is sufficiently small, and thus five-mode models (N = 5) are used in all the simulations to follow.

The bifurcation diagram for the symmetric cubic model is displayed in Figure 2. It is clearly seen that the Hopf bifurcation occurs at U 8 = 1• 78, as predicted earlier, and thereupon stability is lost, leading to the birth of a limit cycle; for U < U 8 , all oscillations, if any, die out as time becomes large. It can also be seen that the first post-Hopf bifurcation occurs at Ub = 3• 24, followed by another one at Uch = 4• 39. In order to have a clear picture of the system dynamics in the neighbourhood of these bifurcations, phase flows were computed for some specially chosen U, and are shown in Figure 3.

A regular limit cycle for U=3•0 (U 8 < U< Ub) is shown in Figure 3(a). After the first post-Hopf bifurcation, the cylinder vibrates slightly once, before it "leaves" the springt t Strictly speaking, the cylinder is always in contact with the spring for the cubic model of the constraint.

However, the spring force versus deflection characteristic has the followng feature: the spring force is very small for small deflections and then abruptly very high for large deflections; hence, one may use such words "leaving" and "hitting" the spring. and moves to the other side, as seen in Figure 3(b) for U = 3• 6 where the motion is symmetric. A further increase in U drives the cylinder to hit the springs harder, and the corresponding phase flow, as shown in Figure 3(c) for U=4•0, looks still periodic but may have higher period. In Figure 3( d) is shown the phase flow for U = 4• 3 before the occurrence of chaos; it is seen that sudden changes of the direction of the phase flow appear. As U increases further, chaotic response develops, as indicated by the orderless phase flow shown in Figure 3(e). Regardless of how long the simulation was carried out, no pattern could be observed. Qualitatively speaking, this phase flow is similar to the flow for a two-well potential oscillator, where homoclinic orbits exist and the homoclinic bifurcations lead the system to chaos, as seen in Figure 4 of reference [START_REF] Moon | [END_REF]. t

A common feature of chaotic systems is that, within the chaotic regimes, there exist windows in which the system response is periodic. In Figure 3(f) is shown such an example for the present system (U=4•82). This periodic window can also be seen fairly clearly in Figure 2, buried in the chaotic region. It should be commented that chaotic response in physical systems may become periodic (transient chaos), and periodic motion surrounded by chaos may also become chaotic with the increase in time. The phase flow in Figure 3(f) was recorded in a time period of 500-520 non-dimensional units.

L Y APUNOV EXPONENTS

There are many ways in which to characterize chaotic motions. Phase flows, Poincare maps and power spectra, to name but a few, are common techniques utilized for this purpose, but are only useful for each individual signal. To study the system dynamics over a range of a given parameter, for example the flow velocity, the appropriate tool is the Lyapunov exponent technique, especially useful in numerical simulation. Unlike some of the other aforementioned techniques, this one gives a quantitative measure for the onset of chaos. For each given parameter value, one can compute a Lyapunov exponent for the corresponding motion. The motion is attractive (a stable fixed point), periodic or chaotic, if the exponent is negative, zero or positive, respectively. This simple procedure can easily t The underlying mechanism which makes the phase flows here similar to those in Moon's non-autonomous problem [START_REF] Moon | [END_REF] could not be pinpointed. Nevertheless, simialr structures can be observed in other systems which, on the surface at least, are quite different: see, for example, Ueda's (18) and Plaut and Hsieh's [19] work. be repeated for different values of the parameter over the range of interest, and the results saved and looked at later.

Lyapunov exponents for a set of delay differential equations, such as those for the system under study, can be derived in the same spirit as that utilized for ordinary differential equations. To proceed with this derivation, the delay differential equation in the form

x(t) = f(x(t)) + g(x(t-r)) (9)
is considered, where x represents ann-dimensional vector function, andfand g are analytic functions of x at t and t-r, respectively. Denoting by l/J(t) a solution of equation ( 9), and introducing a variational vector function u(t) such that x(t) = l/J(t) + u(t), then substituting x(t) into equation ( 9) and Taylor expanding yields u(t) = Df( l/J (t))u(t) + Dg( 4J (t-r))u(t-r), (10) where Dfand Dg are the Jacobian matrices of/and g along the trajectory c/J(t), and the non-linear terms have been discarded. With cjJ (t) known, by numerically integrating equation (9) for given initial conditions, the solution of ( 10) gives information concerning the stability of c/J(t). For a bounded solution c/J(t), its Lyapunov exponents are defined by

a-+.!_ In iu(T)i T iu(O)i as T-+ +oo. ( 11
)
There are n values of a for each c/J, since n linearly independent u(O) may be chosen in !Rn. However, only the maximum exponent, denoted simply by a, is important for detecting chaos, and this a may be computed by assigning an arbitrary u(O) [20].

The computation of Lyapunov exponents for the cubic model was carried out by first writing equation ( 6) in its first order form and then deriving the variationals: tj;==p;, p;= -(al + azU)p;-a3qi-a4U 2 q;(tr)-K(7 cPihr c/J;, V;= -(a!+ azU)v;-a3u;-a4U 2 u;(tr)-K(7 c/J,4 r(7 c/Jjuj)c/J;, [START_REF] Price | Flow-induced Vibration-1989[END_REF] where the only delay terms are those associated with the coefficients 2 shows that the system indeed undergoes a Hopf bifurcation at U H and the transition to chaos occurs at Uch. This result also confirms the validity of the application of the Lyapunov exponent technique to delay differential systems. However, one may have noticed that the periodic window in the neighbourhood of U=4•8, shown in Figure 4(a), is considerably larger than the one shown by the bifurcation diagram in Figure 2; an alternative possibility is that this region contains numerous thin periodic regions interspersed between equally thin chaotic ones. Individual checks of some points in this neighbourhood showed that some seemingly chaotic motions (observed from the phase flows) do have near-zero Lyapunov exponents. Such motions may well be quasi-periodic, as will be suggested by an analysis for a simplified model in section 8. Notice that the original definition of the Lyapunov exponent, equation ( 11), is over an infinite time interval. This, of course, is not practical in numerical simulation and was never achieved. In practice, the result obtained from calculations over a finite time period is perfectly acceptable if convergence is reached. This averaging process is shown in Figure 4(b) over a time interval of 800 non-dimensional units (that over the first 200 was discarded) for U== I• 5, 3•0 and 4•61, respectively; these three cases correspond to attractive, periodic and chaotic motions.

There are many technical subtleties in the implementation of the computation of Lyapunov exponents for delay differential equations. The most important one is the process of renormalization during the computation. Notice that when a> 0 (in the chaotic region), the solution to the variational equation (10), which is linear, will blow up as time approaches infinity. Numerically, one corrects this problem by renormalizing u(t) within some finite time intervals. When doing so, one must also renormalize the history in (t-r, t), since this information is needed for the computation later. With the results discussed in this section it is considered that the main question raised in this study has been answered: chaotic oscillations are indeed possible for sufficiently high flow velocities.

BIFURCATION DIAGRAM FOR THE ASYMMETRIC CASE

So far, the necessary tools have been developed for characterizing chaotic responses and they have successfully been applied to the symmetric model ( ;b = 0• 5). Now it is of interest to investigate what changes would result if the stop constraint was placed off the midcylinder point. It is noticed that, when ;b = 0• 5, due to symmetry the even modes do not contribute, and thus do not affect the system behaviour; they can be entirely ignored. When ;b =!0.5, however, full participation of all modes is expected in the simulated dynamics.

In Figure 5 are• displayed the bifurcation diagram and the corresponding Lyapunov exponents in the velocity range of 3•0-4•5 for ;b=0•45. The Hopf bifurcation remains unaffected, at UH= 1•78. Then, at Ub=3•39, another bifurcation occurs, and the cylinder now begins to oscillate one more time "locally" before it moves to the other side of the spring, this being analogous with the symmetric case (cf. Figure 3(b)). At Uch=4 •12, chaotic motion begins to appear, and the solution blows up for U near 4• 5. A major difference between the two bifurcation diagrams, Figures 2 andS Lyapunov exponent diagrams, Figures 4(a) and 5(b), is that the periodic window just beyond the onset of chaos, observed in the symmetric case, did not materialize in the asymmetric case. A typical chaotic phase flow diagram for the asymmetric case is shown in Figure 6.

NUMERICAL RESULTS FOR THE TRILINEAR MODEL

Up to this point, the constraint force has been modelled by a cubic spring,/(17) = tC1J 3 , which represents a compromise: although this form off is mathematically convenient, it does not model the physical situation perfectly, since it does not recognize that in the free gap, the constraint force is zero, rather than just small. Therefore, in this section, the constraint is modelled by a trilinear spring,j=TC[1J-~(I1J+di-11J-dl)], which is obviously closer to physical reality when TC and d are chosen properly. It is noticed that the representation of the trilinear stiffness is continuous but not analytic at lxl = d (i.e., at the contact points), and thus no Lyapunov exponents can be computed from the governing differential equations.

The phase trajectory for U= 2• 5 (> UH) is shown in Figure 7. It was simulated by setting a small initial displacement ( 1J <d) and zero initial velocity. With increasing time, the phase trajectory begins to grow until it impacts with the stop constraints. Thereafter, the cylinder vibrates more irregularly back and forth within the gap, repeatedly striking the two sides of the loose support. In an attempt to see any possible convergence to a limit cycle, the simulation was performed for the rather long time period of 1000 non-dimensional time units; however, no periodic pattern was found.

In Figure 8(a) and (b) are shown the phase flows for U= 1•9 (just after the Hopf bifurcation) and U = 3, with the transient trajectories of the first 500 time units discarded. For the case of U= 1•9, the response appears to have a small chaotic component superimposed on a larger periodic component indicated by the main peak of its power spectrum, shown in Figure 8(c). For the case of larger flow velocity, U=3•0, the phase plot looks more chaotic (Figure 8(b)), and the corresponding power spectrum is broad-banded (Figure 8(d)). The same trends continue for higher flow velocities. The increase in response frequencies with flow is of particular importance, since that would result in increased wear of the tubes in heat exchangers.

Although the constraint force is modelled by a trilinear stiffness with TC = 1000, the system essentially behaves as if it were impacting on two rigid stops. This may be seen from the phase flows in Figures 7 and8: when the cylinder hits one side of the constraint, it immediately bounces back. The resulting motion appears to be chaotic, without any pattern. This behaviour is similar to that of a zero-stiffness impact oscillator where chaotic motions are possible even for very small forcing [21]. A similar phenomenon was also observed for a pipe conveying fluid when the motion-limiting constraint was modelled by a trilinear spring [22]. In general, for a multiple-degree-of-freedom non-linear system with delay terms, the theoretical analysis is difficult, and may be practically impossible, and one is usually limited to numerical simulation. For the present system with a trilinear spring model, it may be useful to simplify the system further, to a point which will allow some 3•0 ,;;-

1•5 'g ~ 1!-0•0 u 0 ~ -1•5 -3•0 -2 -1 I 2
Displace men! ( x 10-3 ) meaningful analysis and could provide some insight regarding the underlying mechanism leading to chaos.

Considering all the results obtained in sections 4-7, the situation may be summarized as follows. With both models of the impact stiffness and for the symmetric and asymmetric loose-baffle location, it has been shown that, for sufficiently high cross-flow, chaotic oscillations are possible. Hence, perhaps the erratic, chaotic-looking oscillations reported by Axisa et al. [START_REF] Axisa | [END_REF] were chaotic also.

Nevertheless, there remain several unanswered questions. First, the route to chaos for the present models is still unclear, although qualitative changes in the bifurcation diagrams can be clearly observed. The second question is the significance of the apparently chaotic behaviour obtained for some parameter ranges with the symmetric, cubic-spring model (Figure 2 for u~4•7-5•0), which are shown to be periodic by the Lyapunov exponent calculations (Figure 4(a)). The third question is related to the result obtained with the trilinear spring, where it was found that chaotic motions appear to occur at any U large enough to cause impacting. It will be attempted to provide possible answers to these remaining questions in the next section, through an even further simplified, lower-dimensional model.

ANALYSIS WITH A SIMPLIFIED MODEL

In its simplest form, the analytic model utilizeld in the foregoing has only two modes: one which for high enough U develops negative damping, and the other with positive damping; the two are coupled via the impact-related non-linear stiffness terms. The modes concerned are the lowest two interactive modes: for the symmetric case, the first and the third; for the asymmetric case, the first and the second. Here, the simplest possible model is proposed, namely a one-degree-of-freedom model involving a negatively damped, forced impact oscillator. The negative damping feature is carried over directly from the previous model, and is associated with the lowest mode behaviour therein. The positively damped higher modes are represented in the model by the forcing function, where it is understood that the forcing frequency ought to be substantially higher than that of the oscillator; finally, the damping associated with the higher mode(s) is represented by a restitution coefficient less than unity. This forced impact oscillator cannot pretend to model the physics of the real system; nevertheless, it does possess some similar key features, and it is in this sense that it is pursued.

The mathematical model for the impact oscillator described above is a second order piecewise linear differential equation for which there exist explicit solutions for the time intervals between any two consecutive impacts. The non-dimensional equation of the oscillator takes the form x+2S"x+x=ycosmt, for lxl<l,

.X-+ -rx, for lxl =I, (13)
where S" ( < 0) is the negative damping coefficient, and r ( < I) is the restitution coefficient, representing the energy dissipation in the system, m (> I) and r represent the forcing frequency and amplitude, respectively; the system is non-dimensionalized such that impacting occurs at lxl = I rather than at the real physical stops.

It should be remembered that there have been many studies on piecewise linear impact oscillators, notably by Shaw [START_REF] Shaw | The dynamics of a harmonically excited system having rigid amplitude constraints[END_REF] and Li et al. [21]. The system (13) looks rather simple, but its dynamics are extremely rich and complicated, and could themselves be the subject of a separate study (the case of S" = 0 is reported by Shaw [START_REF] Shaw | The dynamics of a harmonically excited system having rigid amplitude constraints[END_REF]). Here only some key results are presented which may be relevant to this problem, and which may provide possible explanations to some of the unanswered questions referred to at the end of the previous section. Before proceeding, the Poincare map technique for system (13) is introduced, as it is needed in later analysis.

THE POINCARE MAP

Equation ( 13) is three-dimensional, and may be written into a set of first order equations by introducing a phase angle 4>:

y=-2S'y-x+ycosm4J, for lxl<l, ~=1, y-+-ry, forlxl=l, x=y, (14)
where 4> = t(mod 2tr I m) is simply the time t, except that it only takes values from 0 to 2tr I m. The central idea of the Poincare map technique is to choose a cross-section and convert a continuous system into a discrete one on this surface, reducing the system dimension by one. Equation ( 14) is discontinuous at lxl =I, and thus it is convenient to choose such a cross-section at x = 1 (or -1):

L = { (x, y, </J ) E fiX~+ X §: x= + 1, y > 0, </JE(0, 2tr jco ]},
where U = [ -1, 1], and § is the circle of 2trjm. The variables in the Poincare section are 4>

and y. The Poincare map is a mapping from I: onto I:, defined symbolically by

( </>, y) = P( </Jo, Yo). ( 15
)
In words, P takes the point ( f/1 0 , y 0 ) in I: back to another point ( f/1, y) still in I: via equations [START_REF] Fricker | Proceedings ASME International Symposium on Flow-induced Vibration[END_REF]. If the orbit is periodic, then (f/J,y)=(f/Jo,Yo).

The advantage of employing this technique for system ( 13) is that thereby one is effectively dealing with a two-dimensional discrete map rather than a three-dimensional discontinuous flow, and the stability information of a continuous periodic solution of ( 13) may be obtained by computing eigenvalues of P at the corresponding fixed point. The solution is stable if all IA.;I < 1. The computation of Pis made through obtaining solutions [START_REF] Fricker | Proceedings ASME International Symposium on Flow-induced Vibration[END_REF] in intervals between two consecutive impacts and using the rebound condition at the two boundaries.

DETERMINATION OF PERIODIC SOLUTIONS

Equation ( 14) contains many different types of periodic solutions of order n, depending on the values of the system parameters and the initial conditions. The simplest type of such solutions is the one in which the system impacts only once with each of the two boundaries, and spends a period of 2ntc I ro to complete a whole cycle. This is the type of solution considered in the following.

To find such a periodic solution, it is assumed that the oscillator starts at a pont Po on the Poincare section Po = (xo = I, Yo > 0, c/Jo) e :r.. A complete cycle is as follows: from Po the solution jumps to point p, = (x, = 1, y, < 0, c/J,H:r. via the impact rule, and then it leaves p 1 and flows toward the left until it hits the left boundary at point Pz = (xz = -1, Yz < O,cpz); at that moment, the oscillator once again undergoes a jump from pz to another point P3 = (x 3 = -1, y 3 > 0, cp 3 ), and after that it continues its journey toward the right until it returns to the Poincare section at P4 = ( X4 = 1, y 4 > 0, cp4) e :r.. The motion is periodic when cp4 = c/Jo and y4=yo.

In the following analysis, all parameters will be fixed except for the forcing amplitude y, which will be varied to investigate possible bifurcations. To avoid confusion, t is used still, instead of cp, while keeping in mind that cp = t(mod 2tc I ro ). A general solution of equation ( 14 From the above equations, one can in principle solve for t4 and Y4 in terms of to and Yo, step by step from p 0 to p 4 . It is immediately recognized that, because of the involvement of trigonometric and transcendental functions, such an implicit solution is impossible. It is fortunate that to fulfill the purpose of determining stability and bifurcations, an implicit form of the Poincare map P is sufficient. After all, the values of r and initial conditions (t 0 , x 0 , y 0 ) for periodic solutions can be computed numerically through solving the above non-linear algebraic equations.

In order to determine the stability of a periodic solution emanating from (to, Yo), the Jacobian matrix of the Poincare map P at (to, Yo) must be computed:

a(t, y) DP(to, Yo) = . a( to, Yo) (18)
Because of the discontinuity in y, the computation of the derivative (18) must be divided into four parts: (19) All the individual derivative terms can easily be determined using the impact rule and the functions f and g defined in equations ( 16): and the derivative iJ(t 4 ,y 4 )ja(t 3 ,y 3 ) can be obtained from equation (20b) by shifting the numbers from 1 to 3 and 2 to 4. Multiplying the four matrices, in order, yields DP(t 0 , y 0 ) and because of the assumption that (to, Yo) is periodic, Y4 should be replaced by Yo in the final result. The eigenvalues of DP can be determined from

A. 2 -tr A.+det=O, (21)
where (22) is the determinant of DP, with Ll21 and Ll43 being the time periods spent from p 1 to Pz and p 3 to P4, and tr is the trace of the matrix, which is too long to be given here. For a given y, the initial conditions for a periodic solution can be computed numerically from equations [START_REF] Moon | [END_REF]. In turn, the characteristic equation ( 21), involving the initial conditions (t 0 , x 0 , y 0 ) and y, may be solved to give the eigenvalues A..

RESULTS

All calculations in this section have been conducted with m= 3• 5. !; = -0•05 and r=0•75.

The main objective here was to find whether a period-doubling sequence exists, which is one of the common paths leading to chaos for many physical systems. For a discrete map, such a bifurcation sequence is often predicted by a ( -1) eigenvalue of the map. To test whether this phenomenon was occurring, A was first set at }.,=-0•75 in equation ( 21) and then the resulting equation was solved, yielding y=0•9010. The initial conditions corresponding to this periodic solution were t 0 = 0• 8697, x 0 = + 1 and y 0 = 1•0336. The phase flow for such a periodic solution is shown in Figure 9(a). Since IAI < 1, the solution was stable and, regardless of how long one observed the simulation, it never changed. By then setting }., = -1• 25 and repeating the same procedure, y = 0•9329, to= 0• 8728, x 0 = + 1 and y 0 = 1• 0378 were obtained, corresponding to an unstable solution. In this case, the phase flows over a short time period looked the same as that of the previous case, but as time was increased the flow gradually diverged to a new steady state. The result is shown in Figure 9(b) after discarding the transient solution. It is clear that a period-doubling bifurcation has indeed occurred at y c = 0• 9133, obtained directly by setting }.,= -1 in equation ( 21) and solving it.

In order to determine whether this single period-doubling bifurcation event is the first in a doubling sequence, a bifurcation diagram was systematically simulated, displaying the steady state positive velocities (y > 0) at x = 0 versus the parameter y, over a fairly large range, 0•05< y< 1•2. The result, plotted in Figure lO(a), is in remarkably accurate agreement with the prediction of Ye• However, as was commented at the beginning of this section, the periodic solutions being investigated in the foregoing analysis are but a small portion of the whole group of periodic solutions. The whole picture of the system dynamics is far more complex. This is clearly seen in Figure lO(a): even before the predicted bifurcation at Ye, the system response has already become chaotic. In Figure lO(b) is shown an enlarged portion in the neighbourhood of y = 0• 67, and one sees that there exist some higher order period-doubling sequences involving impacting more than once with each The importance of these results is in their application, by analogy, to the original problem. The last set of results (Figure 11) shows clearly that very-long-period motions can arise which can easily be mistaken for chaos (cf. results of section 4). The other point is that despite the apparent disorder in Figure 10, the underlying mechanism leading to chaos is through a cascade of period-doubling bifurcations. The qualitative similarity of Figures 10, 2 and 5(a) is evident, and hence the same conclusion with regard to the route to chaos may apply to the earlier model as well. Finally, the results obtained by the trilinear model, giving chaos immediately after impacting begins, may now be seen as being associated with the immediacy of multiple impacting and quasi-periodic motions. In principle, there may be a period-doubling sequence leading to this, but the interval of the varied parameter ( U) over which this occurs may be too small to be pinpointed.

CONCLUSIONS

The purpose of this paper was to decide whether chaotic oscillations are possible in multi-span heat exchanger tubes subjected to cross-flow, when there exists at least one loose support, a common occurrence in practice. This question was addressed with the aid of a simple model involving one two-span flexible tube in the midst of an array of rigid cylinders, with the middle support on the flexible tube being loose, i.e., with a gap between the tube and the supporting baffle plate. To make matters simple, motion was supposed to occur exclusively in the cross-stream direction, and impact dynamics was modelled via cubic and trilinear spring representations.

The study involved numerical computation of bifurcation and phase flow diagrams and Lyapunov exponents. It was shown that, as the flow velocity was increased to beyond the first Hopf bifurcation (negative-damping fluid elastic instability), limit-cycle motion resulted; beyond that, impact of the cylinder on the supports gives birth to chaos when the flow velocity becomes sufficiently high. With the trilinear spring model of the impact constraint, chaotic-looking results are obtained, immediately after (in terms of increasing flow velocity) the occurrence of impacting on the constraint.

The first and foremost contribution of this paper is to show, for the first time and in quantitative terms, that chaotic oscillations are indeed possible for this system; i.e., for loosely supported heat-exchanger tubes subjected to cross-flow. This result is of practical importance also, since the higher frequencies that are associated with such motions have been known to cause accelerated wear at the loose support, and hence to signal premature breakdown in heat exchangers.

The second contribution made by this paper is that, to the authors' best knowledge, this is the first application of the Lyapunov exponent technique to delay differential equations; it was through the results obtained thereby that the existence of chaos could be decided upon unequivocally. Of course, the importance of this technique has long been recognized; its implementation to delay differential equations broadens its applicability even further. In the present problem, the computation of the Lyapunov exponents for the case of symmetric location of the constraint clarified an additional point: whereas the system response in phase flow and bifurcation diagrams over a small range of U appeared to be aperiodic and chaotic-looking, Lyapunov exponent calculations showed this to be in fact quasi-periodic motion.

A number of small points remained unanswered, and the 10-dimensional (five-degreeof-freedom) model utilized in the calculations did not lend itself to analytic or semi-analytic methods for clarifying these points. Among these points were (i) the determination of the route to chaos and (ii) the finding that for the trilinear spring case chaos apparently occurred at any flow velocity at which impact could take place.t For that purpose, a highly simplified one-degree-of-freedom mathematical model was implemented in Section 8, involving a forced impact oscillator with negative damping. This being a three-dimensional model only, it was possible to determine semi-analytically the mechanism underlying the onset of chaotic oscillations: a period-doubling cascade was predicted and confirmed numerically. The results suggest that even though the actual route to chaos for the physical system was not clearly demonstrated, period-doubling cascades may well be the underlying structure for chaos. Also, the results obtained with this model suggest that the chaoticlooking response obtained with the trilinear model of the impact constraint beyond, but close to, the Hopf bifurcation point may in fact be quasi-periodic.

Figure

  Figure I. (a) The system under consideration, showing the flexible cylinder positively supported at its extremities and with a loose support at xb, subjected to cross-flow. (b) A cross-sectional view of the idealized system, in which the cylinder vibrates only in the cross-stream direction and the loose support (gap df(..D) is modelled as shown.

2 . 1 B=

 21 THE EQUATION OF MOTION The equation of motion of the cylindrical beam subject to cross-flow may be written in the form a\v aw o 2 w E l -4 +c-+m-2 +8(x-xb)f(w)=F(w, w, w), ax ar at (I) t The dynamics of generalized impacting of a circular cylinder in a circular hole can be very complex, involving both Coulomb and sliding frictional effects. where w is the cross-stream lateral displacement of the cylinder, El its flexural rigidity, m its mass per unit length, and c is the damping coefficient;fis the non-linear force due to the constraint, with o(x-xb) being the Dirac delta function and xb the constraint location (see Figure l(a)), and Fis the flow-induced force. The cylinder is supposed to be clamped at its two extreme supports (Figure l(a)) rather than simply supported, to disallow natural frequencies which could be integral multiples of the first. t The flow-induced force F may generally be expressed as . .. (fw ow F(w, w, w)=M-2 (x, t)+ B-(x, t)+ Cw(x, t-M), ot otwhere, in accordance with the Price and Pai'doussis[START_REF] Price | [END_REF]8] model,

Figure 2 .

 2 Figure 2. The bifurcation diagram for the system representing an in-line square array with p/ D= 1•5 (CD= 0•26, iJCdil11 = -8•1), m= 3, 81=0•06, J3 =0•24 and K= 10 3 , showing the dirnensionless mid-point displacement amplitude in terms of the dimensionless fluid velocity, U. The loose support, positioned at the cylinder midpoint ("symmetric" system, xb/ L = 0• 5), is modelled here as a cubic spring. U" denotes the critical U for the Hopf bifurcation, and Ub is the first post-Hopf bifurcation; Uch denotes the onset of chaos.

Figure 3 .

 3 Figure 3. Phase flow diagrams associated with the bifurcation diagram of Figure 2, for the followng values of U: (a) U=3•0; (b) U=3-6; (c) U=4•0; (d) U=4•3; (e) U=4•6!; (f) U=4•82.

  a 4 ; the eigenfunctions c/Ji c;) are evaluated at c; = c;b. When an N-mode model is employed in the computation, then i,j== 1, 2, ... , N. The exponents for the symmetric case of section 4 are shown in Figure 4(a), where c;b = 0• 5. It is clearly shown that for U < U H= I• 78, the exponents are indeed negative, and are nearly zero for periodic motions, until the point at which the motion becomes chaotic, at Uch• Direct comparison of the results in Figure 4(a) with the bifurcation diagram in Figure

~TimeFigure 4 .

 4 Figure 4. (a) The largest Lyapunov exponent for the system of Figure 2 as a function of the dimensionless flow velocity, U; (b) the rate of convergence of the calculation for three specific values of U as a function of dimensionless time.

Figure 5 .Figure 6 .

 56 Figure 5. (a) The bifurcation diagram for the asymmetric version (xb/ £=0•45) of the system of Figure 2; (b) the corresponding largest Lyapunov exponents.

Figure 7 .FrequencyFigure 8 .

 78 Figure 7. Phase flow diagram for the system with trilinear modelling of the constraint stiffness at the loose support (at xh;L=0•5). for U=2•5, showing the growth of the oscillation in response to impact with the loose supports, and the subsequent motion.

  ) between two consecutive impacts has the form x = e -w-<ol[a cos roJ.tt 0 ) + b sin rod(tto)]+ {3 cos (rot-'!')= f(t; to, Xo, Yo),[START_REF] Yetisir | Fluid-Structure Interaction and Aerodynamic Damping[END_REF] y=.X=g(t; to, Xo, Yo))-tc + arctan --2 1-(I) for ro > 1,and a, b are integration constants depending on the initial conditions to, Xo and Yo.In order to determine the Poincare map, it is necessary to find a mathematical solution from Po to Y4 e :r.. This can be done as follows: from Po top,: y, = -ryo, t, =to; fromp, topz: -l=f(t2; +I,y,), Yz=g(t2; t,, +l,y,); frompz top3: y3=-ryz, t3=tz;[START_REF] Moon | [END_REF] fromp 3 to p4: +1 =j(t4; t3, -1,y3), y4=g(t4; t3, -l,yJ).

  au1, YI) = act3, y3) [1 o ]. a( to, Yo) a(tz, Yz) 0 -r a(tz,yz) -::; :~ -:~;:~ act~, yJ) ag atz ag --+ -ag at 2 ag --+atz at! at! atz ay! ay!

DisplacementFigure 9 .

 9 Figure 9. Phase flow diagrams for the simplified model of the negatively damped, forced impact oscillator (S'=-0•05, w=3•5, r=0•75): (a) a periodic solution for a forcing amplitude of y=0•9010; (b) a period-two solution for y = 0•9323.

Figure 10 .

 10 Figure 10. (a) Bifurcation diagram for the simplified, three-dimensional model of the negatively damped, forced impact oscillator(~= -0•05, w = 3• 5, r = 0• 75), with the forcing amplitude, y, as the bifurcation parameter; (b) a portion of the bifurcation diagram, 0•64<y<0•70, in enlarged form.

0- 2 '

 2 Figure 11. (a-c) Phase flow diagrams and (d-f) the corrsponding Poincare maps for the simplified threedimensional model of the negatively damped. forced impact oscillator (S'=-0•05, ro=3•5, r=0•75): (a. d) for r=0•5; (b, e) for r=0•7; (c, f) for r= 1•0.

TABLE 1 Critical

 1

	Mode	UH	OJ;
	1	1•785	0•824
	2	4•920	2•272
	3	9•645	4•455
	4	19•944	7•364
	5	23•818	11•000

flow velocity values, UH, and the corresponding non-dimensional natural frequencies

t Note that the values of drag, lift and their derivatives given by Price and Paidoussis[START_REF] Price | [END_REF] have been converted here from those corresponding to far-u~stream flow velocity to values associated with the so-called reference gap velocity, by dividing by [pj(p-D)], equal to 9 in this case.

t Since the trilinear model is not analytic, the computation of Lyapunov exponents, which would have decided the issue, is not possible.

ACKNOWLEDGMENTS The authors gratefully acknowledge the support provided by the Natural Sciences and Engineering Research Council of Canada, Le Fonds FCAR of Quebec and the Faculty of Graduate Studies and Research of McGill University.