domain: math phys

text:
Alain Jégat

ONE-DIMENSIONAL ELASTIC COLLISIONS IN A PLATONIC QUADRIDIMENSIONAL SPACE

whether they are published or not. The documents may come

The geometrical framework

This modeling is based on the Platonic space outlined in the following articles:   , , , , O i j k h is a frame for the four-dimensional Euclidean space whose axes are denoted   OX ,   OY ,   OZ ,  

Ow ; the direction of the projection is that of the vector h .

Following the hal-01207447 v1 and hal-01213447 v1 articles, the notion of relativistic mass of a particle is described here as a result of its interaction with a stratification of the four-dimensional Platonic space by a sequence of hyperplanes

Introduction of some absolute concepts

This section is intended to introduce some concepts, described as absolute, which are independent of any moving reference frame in the Platonic space.

The application of these absolute concepts will allow us to obtain the relations between the initial and final trajectories of two relativistic particles after a one-dimensional elastic collision.

Absolute time concept

Given the regular motion of the objects observed in the Platonic universe (« between two events, the distances traveled by all the observed objects are equal, regardless of their trajectory »), we can introduce a notion of absolute time in the following terms :

« Absolute time between two events is the distance traveled by (all) the mobiles between these two events .» Thus, in the above diagram, if the event E1 occurs while the mobiles M, N and P are respectively at M1, N1 and P1 and if the event E2 occurs while the mobiles M, N and P are respectively at M2, N2 and P2, we can define the absolute time T  between these two events by the relation Tr    .

Note that in this Platonic space, one can thus introduce a notion of absolute time which is measured in meters.

Absolute velocity

Into the frame   ,, O i h , let us consider a ponctual mobile M whose straight trajectory is directed by a vector i  , with   Note that here the absolute velocity of M is a dimensionless quantity.

Absolute frequency of De Broglie waves

Let us consider a ponctual particle P moving along a straight line led by a vector i  , with  

, ii    , linked to a stratification of the Platonic space whose distance is 0 w  .

To this particle may be associated De Broglie waves (see the articles cited in reference).

Consider a line whose equation is Xa  .

We can define the absolute frequency of these De Broglie waves considering two of its successive occurrences Note that this absolute frequency is expressed in

1 m  .

Absolute energy and absolute mass of a particle

 By adopting for absolute quantum of action the quantity abs h hc  (in .

Jm), associated to the absolute frequency of the De Broglie waves

0 1 sin abs w     (in 1
m  ), we are led to define as absolute energy of the particle P (in J ) the quantity :

0 sin abs abs abs hc Eh w     .
The link between the energy abs E and the energy E of the particle measured in a reference frame R  (where  is arbitrary, not a multiple of ) is given by :  Into the reference frame R  , is so obtained the following coherence :

2 0 1 cos cos sin sin hc E h mc w          , with 2 0 2 sin 1 cos cos sin sin 1 abs c w v c           and 2 0 2 sin 1 cos cos sin sin 1 abs m h m cw v c          .

Absolute linear momentum of a particle

By analogy with the classical notion of momentum of a particle, we define its absolute momentum as being equal to the product of its absolute velocity by its absolute mass, (Note that this absolute momentum is expressed in kg.)

Application to the one-dimensional elastic collisions

Two hypotheses

h h h h c w c c w c w            .
With the relations (see articles cited by reference), these two assumptions will lead to the results described in the following paragraph.

One-dimensional elastic collisions -Diagram 1 --Diagram 2 --Diagram 3 -

Let us consider two particles P1 and P2 whose rest masses are respectivly m01 and m02 coming into elastic collision in the Platonic space at the point « Impact » (see Diagram 2).

This impact is perceived by the observer I1 into the reference frame R  ; the trajectory of this observer is represented by a blue dotted line.

(Note that in the Platonic space, this collision can be represented with distinct values of w(P1) and w(P2). Here , they are chosen equal in Diagram 2 by pure graphical convenience.)

Before the collision, in Diagram 1, the trajectory of P1 is directed at an angle  , the trajectory of P2 at an angle  . The particles P1 and P2 are observed into R  at P1 and P2.

After the collision, the trajectory of P1 is directed at an angle '  , that of P2 at an angle '  .

The particles P1 and P2 , designated by P'1 et P'2 , are observed into R  at P'1 and P'2.

From the two assumptions made in the previous paragraph , these two angles are given, with 0   and 0  , by the following equations (cf. expressions with the rest masses at the end of this article):

     2 2 2
                            2 2 2
                       .
It is noteworthy that on the one hand, these relations are absolute to the extent that they do not depend on the angle of  of the reference frame.

On the other hand, they lead, regardless of the value of the angle of the reference frame, to the conservation in this reference frame of the global energy and global momentum, namely:

' ' 01 1 02 2 01 1 02 2 2 2 '2 '2 1 2 1 2 2 2 2 2 2 2 2 2 01 02 01 02 2 2 '2 '2 1 2 1 2 2 2 2 2 1 1 1 1 1 1 1 1 m v m v m v m v v v v v c c c c m c m c m c m c v v v v c c c c                        i.e.
          .   If 0   (i.e. if   1 v P c  ), then         22 02 01 22 02 01 1 cos 1 cos cos ' 1 cos 1 cos ww ww               and '0   (i.e.   2 ' v P c  ).  If   (i.e. if   1 v P c  ), then         22

Expressions of cos(') and cos(') based on the rest masses

«

  UN MODÈLE PLATONICIEN (EUCLIDIEN-PROJECTIF) POUR LA THÉORIE DE LA RELATIVITÉ RESTREINTE » (pré-publication hal-01081576, version 1). « A PLATONIC (EUCLIDEAN-PROJECTIVE) MODEL FOR THE SPECIAL THEORY OF RELATIVITY » (pre-publication hal-01165196, version 1).

H

  which are orthogonal to the direction of the projection h , regularly spaced by a distance

  distinct lines Xa  and Xb  . With X b a    , we can define the absolute velocity of M by :



  From this result, we can define the notion of absolute mass (in kg): wavelength of a particle, m0 being the rest mass of this particle P.)

  Let us consider two ponctual particles P1 and P2 whose straight trajectories are directed respectively by the vectors i  et i  , associated to the stratifications whose distances are

	respectively	w 	and	w  .
		01		02
				w 	sin hc			w 	sin hc			w 	sin ' hc 		w 	sin ' hc 	;
				01			02			01		02
	 the sum of absolute momentum is conserved,
		i.e.
				tan				tan				tan '		tan '
				01				02				01		02

After the elastic collision of these two particles P1 and P2, if we call ' i  and ' i  the direction vectors of their new trajectories , we will choose as hypotheses :

 the sum of absolute energy is conserved, i.e.

. Study of some particular cases

	 If	01 w  or if 0	02 w 	  , then '   and	cos ' 		 cos 1 cos 2 cos cos 2     2  2 cos 1 cos 	 		.
	This last result, of abstruse appearance, leads just to the following fact: in any reference
	frame R  where the particle	1 P is at rest, the speed of	2 P and	2 ' P are opposite .
	Indeed, in this reference frame, we have (cf. hal-01165196, v1) :	  2 vP c 		cos 1 cos cos cos    
	01 01 sin sin '  02 01 02 02 01 02 cos cos cos cos cos ' cos cos ' cos sin sin sin sin sin sin '                 and there is obtained, after simplification:   2 ' cos ' cos cos 1 cos cos ' 1 cos cos cos vP c              	.
	3.3 First, if	ww   	(i.e.
			01	02
	 Symmetrically, if	w    or if	0  , we obtain w	'   and
						01				02
			2 cos		 cos 1 cos	2	
	cos '
				2
			1 cos			2 cos cos

01 02 mm  ) we obtain, independently of the reference frame R  , a geometrical symmetry whose simplicity is noteworthy . Indeed, in this case, the above formulas lead to '   and '   .

Conclusion

Through this study, this article aims to show that the introduction of absolute concepts (energy, momentum , ...) can easily model relativistic phenomena in the Platonic framework .

From this, it seems promising to continue their application (eg to study the elastic twodimensional relativistic collisions or inelastic collisions), or even to increase their number (eg by introducing the concept of absolute angular momentum) to solve the modeling of other problems, as the two-body problem .

And beyond is needed of course the search for deeper understanding of the links between these concepts of absolute energy, momentum, angular momentum, … and the geometrical structures proposed into the Platonic model.


L Morelli
email: luca.morelli@bo.infn.it

M Bruno

M. D 'agostino

G Baiocco

F Gulminelli

M Cinausero

M Degerlier

D Fabris

F Gramegna

T Marchi

S Barlini

M Bini

G Casini

N Gelli

G Pasquali

S Piantelli

S Valdrè

Clustering effects in fusion evaporation reactions with light even-even N=Z nuclei. The 24 Mg and 28 Si cases

Clustering effects in fusion evaporation reactions with light even-even N=Z nuclei. The 24 Mg and 28 Si cases

INTRODUCTION

The NUCL-EX collaboration has recently started an experimental campaign of exclusive measurements of fusion-evaporation reactions with light nuclei as interacting partners. The aim is to progress in the understanding of statistical properties of the decay of light nuclei at excitation energies above particle emission thresholds and to measure observables linked to the presence of cluster structures in nuclear excited levels.

The statistical decay of hot nuclei and nuclear clustering are very active research topics in nuclear physics [START_REF] Valdrè | Clusters in Nuclei[END_REF]. Using the statistical theory of Compound Nucleus (CN) decay, the detailed output of a fusion-evaporation reaction is uniquely predicted under the knowledge of nuclear ground state properties and level densities. The knowledge of level densities is not only important for the understanding of nuclear structure [5], but it is also required for different applications of nuclear physics, from nucleosynthesis calculations to reactor science.

Despite the interest of this issue, few studies, mainly based on inclusive experiments, exist on the evaporation of very light nuclei (A ∼ 20 region) at relatively high excitation energy (ε * ≈ 3 A.MeV). The interest on this mass -excitation energy region is easily justified: from the experimental point of view, measuring light nuclei means low multiplicity events, which, together with a high detection coverage and high energy and angular resolution, leads to the possibility of achieving a complete event reconstruction, thus having a full control on the reaction mechanism.

In addition, the knowledge of the decay mechanism could help in the backtracing procedure in multifragmentation events. Indeed in the statistical framework fragments formed in the freeze-out volume are mainly in this energy and mass region.

Nuclear structure signatures are especially evident in light nuclei, even at high excitation energy. By using an exclusive channel selection and a highly constrained statistical code, it is possible to put into evidence deviations from a statistical behavior in the decay of the hot sources formed in the collision. Indeed some excited states of different nuclei in this mass region are known to present pronounced cluster structures. These correlations may persist in the ground state along some selected isotopic chains [6] and, according to the Ikeda diagrams [2], α-clustered excited states are sizeably expected at high excitation energies close to the multi-alpha decay threshold in all even-even N = Z nuclei.

The subject of α-clustering has been a central issue in nuclear physics and has witnessed a gain of interest in recent years. On the theoretical side, highly sophisticated ab-initio calculations have shown pronounced cluster features in the ground state of a large number of light nuclei [6]. Concerning experimental research, rotational bands consistent with α-cluster structures have been identified in different even-even light nuclei and shown to persist even along their isotopic chains. Such effects can be experimentally seen as an excess of cluster production with respect to the prediction of the statistical model, provided that the ingredients of the latter are sufficiently constrained via experimental data. A campaign of exclusive measurements of fusion-evaporation reactions with light nuclei as interacting partners is currently carried on by our collaboration. In the framework of this campaign, the 12 C+ 12 C and 14 N + 10 B reactions have been measured in order to study the decay of the same 24 Mg compound nucleus, populated at the same excitation energy but through different entrance channels. The results obtained in our first measurements are briefly reported in the following sections.

A very interesting clustering effect is for instance that of the so-called Hoyle state, i.e. the excited state 0 + at 7.65 MeV of 12 C, which plays a decisive role in stellar nucleosynthesis of 12 C.

Recently different experiments have studied the low-lying 12 C excited states and in particular the Hoyle state at 7.65 MeV excitation energy [7] and several interesting features of these states have been pointed out [8][9][10][11][12][13][START_REF] Quattrocchi | Contribution to Nucleus nucleus 2015 NN2015[END_REF].

In several cases experimental results have interpreted the Hoyle state as mainly sequential through an intermediate 8 Be gs decay [8][9][10][11][12]; other experimental data show sizeable amount of instantaneous three α-particle decay [13,[START_REF] Quattrocchi | Contribution to Nucleus nucleus 2015 NN2015[END_REF] where the three α particles are characterized by low kinetic energies and low kinetic energy dispersion. Most of the results indicate a nearly complete agreement with a sequential decay with a very small amount (0.2 % [12]) of the instantaneous one.

In this work we show some preliminary results from the analysis of semiperipheral 12 C + 12 C reactions where the projectile decays in three α-particles passing through the Hoyle state.

The results compatible with indications for non-sequential decay have been obtained in reactions involving heavy ions [13,[START_REF] Quattrocchi | Contribution to Nucleus nucleus 2015 NN2015[END_REF]. We have therefore investigated the Hoyle state when obtained in the decay chain of 24 Mg* formed in central collisions and in particular we have analyzed six α-particles in the final stage of the 24 Mg* decay. We will compare the experimental results both for central and peripheral collisions with the prediction of a model based on Hauser-Feshbach formalism which includes light nuclei ecxcited states. The main features of this model have been extensively described in [4]. The predictions have been filtered in order to be compared to experimental data, taking into account detector characteristics such as thresholds, energy and angular resolutions. Data have been also compared to filtered simple model calculations for simultaneous decay (DDE), i.e. three α-particles simultaneously emitted from the excited Carbon and for a linear chain (DDL), i.e. two α-particles emitted back to back with the third α-particle at rest in the 12 C frame.

Similary to the 24 Mg*, the compound nucleus 28 Si* can be formed through fusion of α-cluster stable nuclei in the entrance channel and in the future of light radioactive beams provided by SPES. It is interesting for such systems that they can be studied at different beam energies, thus extracting information on how (and to what extent) structure effects are still at play in the decay of hot nuclei at different excitation regimes.

THE EXPERIMENTS

The experiments were performed at the LNL (Laboratori Nazionali di Legnaro), with 12 C, 14 N and 16 O beams provided by the XTU TANDEM accelerator. We have studied the two systems 12 C + 12 C and 14 N + 10 B, leading to the same compound nucleus 24 Mg* at the same excitation energy. The decay of the 28 Si* compound nucleus has been investigated using 16 O + 12 C reaction at three different energies. For these last reactions the analysis has been focused up to now on the necessary preliminary checks on the quality of collected data. Charge (and mass, where the information is available) identification and energy calibration are in progress.

The apparatus has been described in [START_REF] Bruno | [END_REF]. Here we recall the main features. The experimental setup is composed of the GARFIELD apparatus and the Ring-Counter (RCo) annular detector, fully equipped with digital electronics [16]. The forward polar angles (7 o ÷17 o ) are covered by a three stage apparatus (RCo), ionization chamber (IC), silicon strip s (Si) and CsI(Tl) scintillators. The angular resolution is ∼ ± 0.7 o for the polar angle and ∼ ± 11 o for the azimuthal one for particles detected in the scintillators. The performances of the RCo allow to reach energy determination with the accuracy of percent, charge identification of particles and fragment with threshold as low as 0.8 ÷ 1 AMeV and mass identification of light isotopes with ∼ 6 AMeV threshold through ΔE -E technique in Si-CsI(Tl) and/or pulse shape analysis (PSA)in the CsI [START_REF] Bruno | [END_REF]. The remaning polar angles (θ = 30 o ÷ 170 o ) are covered by a two stage, drift ionization chamber -CsI(Tl) scintillator apparatus (GARFIELD). Energy is determined with accuracy of the order of few percent and reaction products are identified in charge and mass (for light products) with a threshold around 1 AMeV [START_REF] Bruno | [END_REF]. The combination of the two devices allows a nearly-4π coverage, which, combined with a high granularity, permits to measure the charge, the energy and the emission angles of nearly all the charged reaction products, with an excellent discrimination of the different reaction mechanisms.

RESULTS

Central 12 C + 12 C and 14 N + 10 B reaction

The selection of the fusion-evaporation mechanism is based on the coincidence between LCP's and a fragment detected at forward angles (RCo). At the same time we ask the complete charge detection and longitudinal momentum conservation. We compare experimental data to the predictions of a dedicated Monte Carlo Hauser-Feshbach (HF ) code for the decay of the compound nucleus [3,4], explicitly including all the experimentally measured particle unstable levels from the archive NUDAT2(http://www.nndc.bnl.gov/nudat2/). For both reactions the best reproduction of the systematics of the fusion cross sections is obtained assuming two different maximum values (J 0max ) for the angular momentum distribution of the hot fused source [17]. In particular we use J 0max =18 for 12 C+ 12 C reaction and J 0max =15 for 14 N + 10 B reaction and we adopt the same diffusness parameter ΔJ = 2 for both reactions. Energy spectra for protons and α particles detected at GARFIELD angles are plotted in figure 1, for residues of different charge, for the 14 N reaction, and compared to HF [4] calculations and to data for the 12 C+ 12 C reaction. For both reactions, a good reproduction of proton and α energy spectra is achieved in all channels. Some discrepancies are present in events with an Oxygen residue, where the energy tails for α particles are not reproduced by the model, especially for the 12 C+ 12 C case. It has been already shown in [18] that this discrepancy is mostly due to the an extra experimental cross section for channels of the type (2α, 16 O gs/ * ) 1 . Such outgoing channels populated in the 12 C+ 12 C collisions can be attributed to an entrance channel effect, given the α-like structure of reaction partners and produced fragments. Indeed, for the 14 N + 10 B case inclusive distributions are better reproduced by calculations. This seems not to exhaust the total discrepancy with HF calculations, as we infer from the energy spectrum of α particles in coincidence with Oxygen and with 20 Ne for the 14 N reaction. Thanks to the completeness of event reconstruction, an estimate of the dissipated energy [18] can be extracted to further investigate all the channels involving alpha particles in both reactions:

Q kin = 2 i=1 E α i + E Z Res -E beam
, where E α i and E Z Res are, respectively, the laboratory energy of α particles and residues, and E beam is the energy of the incident projectile. Figure 2 displays the obtained Q kin distributions for 12 C+ 12 C (black dots) and 14 N+ 10 B (blue dots) reactions, taking into account the different initial Q-value for the two reactions. In particular we plot only the channels with the maximum α multiplicity associated to the residue of charge Z. In the Q kin distribution for the channel (2α, 16 O gs/ * figure 2d), the two peaks correspond to α-decay chains, starting from the 24 Mg * compound nucleus and leaving a 16 O residue either in its ground state 12 C and Q kin = -0.8 MeV for 14 N reactions correspond to the opening of the 4-body channel 15 O+n + α + α. Neutrons are not detected in the experiments, and the broader distribution observed for lower Q kin values is due to events in which neutron(s) emission has taken place. A difference in the relative population of less dissipative events is evident between the 12 C+ 12 C (black dots) and 14 N+ 10 B (blue dots) reactions in figure 2 2 . In particular, a much higher percentage of (2α, 16 O) events populates the less dissipative Q-value region in the 12 C sample. This larger deviations are evident in the relative population of the different regions for all the even-Z residues, while we observe a very good agreement in Q kin distributions for the odd-Z residue. This difference between the two data-sets confirms a possible larger contribution of direct reactions for the 12 C+ 12 C experiment. Since a residual deviation is observed in figure 2 for α particles emitted in coincidence with an Oxygen in the 14 N reaction, we turn now to an estimation of the possible α clustering effects for both reactions, both in the entrance channel and in the excited 24 Mg. In order to put into evidence the differences in terms of branching ratio, not only between data samples, but also considering the HF predictions for both reactions we present, in Table 1, the most populated channel in the experimental sample for each residue. We can see that the BR (branching ratio) of the dominant decay channels is reasonably well reproduced by the statistical model for odd Z residues, while discrepancies can be seen for even-Z ones.

The evaporation chains leading to a final Carbon or Oxygen or Neon residue show a preferential α decay in both reactions. A possible interpretation of this α excess could be the presence of residual α correlations in the excited 24 Mg or in its daughter nucleus 20 Ne, populated irrespective of the entrance channel of the reaction.

12 C* Hoyle state in 24 Mg* decay: the 6-α channel

We now turn to examine events where the whole available mass and charge is found as α particles. This selection has been performed in [4] and the inclusive results show compatibility with a sequential α-particle emission. To improve the statistics we have analyzed here also events where only five αparticles have been detected and extracted the properties of the sixth α-particle from momentum and energy conservation. We have reconstructed the energy dissipated by the quantity Q kin previously defined (Q kin = 6 i=1 E i -E beam ). In figure 3 (left panel) the Q kin distribution is shown, together with the cut used for the analysis.

In order to reconstruct the intermediate 12 C* disassembling into 6 alphas of the 24 Mg* formed in central collisions we have selected the three (out of six) α-particles corresponding to the minimum The lowest energy peak in the lower panel of figure 3 corresponds to the very well known Hoyle state, which has been studied in past years by several authors [8-11, 13, 14]. The main debate consists in the interpretation of the decay of this state as sequential, via the 8 Be gs formation, or if there is a contribution of instantaneous breakup. Different observables have been proposed to clarify this point and we will show our results for some of these observables.

The first observable is the minimum relative energy (see figure 4) of two out of the three α-particles assigned to an excited 12 C* produced along the 24 Mg* decay chain. Despite of the low statistics, the results are compatible with the ones obtained for HF , indicating that the decay mainly proceeds through the intermediate formation of a 8 Be gs . The same is true for the Dalitz plot where data and HF predictions show the same configuration (see figure 5). As usual one defines the two coordinates of the Dalitz plot as x d = √ 3(e 1e 2 ) and y d = 2e 3e 1e 2 , where e i, j,k = E i, j,k /(E i + E j + E k ) are the α-particle energies in the 12 All the results for central collisions are compared to HF predictions, showing a very good agreement. All the results obtained in inelastic channels in 12 C+ 12 C reaction [START_REF] Morelli | annual report 2014 LNL[END_REF] and in various reactions at different energies [8][9][10][11][12] indicate a very low limit of this contribution of the order of permil. On the other side reactions with heavier ions at higher energies, involving more complex systems, seem to be compatible with a simultaneous decay of some percent [13,[START_REF] Quattrocchi | Contribution to Nucleus nucleus 2015 NN2015[END_REF]. Our results for central reactions with the formation of 12

The Hoyle state in peripheral collisions

With a different selection we repeated the same analysis performed for central collisions, with the aim of investigating the decay mechanism of the Hoyle state in different conditions. We have selected events where three α-particles are detected in the forward cone (RCo) and nothing is detected in the rest of the solid angle. If we want to unambiguously identify "true" α-particles we have to exploit the ΔE-E Si-CsI or the pulse shape analysis in CsI(Tl) scintillators [16]. This implies that the threshold for the identification is about 24 MeV and consequently that the kinetic energy of 12 C nucleus emitting three α-particles is relatively high. The recoil nucleus 12 C is therefore very slow and its energy not sufficient to reach the detector, even with a very thin target. In figure 6(left panel) the energy spectrum of the reconstructed 12 C quasi-projectile is shown, as calculated from the sum of the energy of the three α-particles:

E tot = 3 i=1 E i + E rec
where E rec is the recoil energy of the 12 C quasi-target. As it can be easily seen most of the kinetic energy is peaked at ≈88 MeV and a small peak is present at about 84 MeV, corresponding to the 12 C quasi-target at about 4.4 MeV excitation energy, the first 12 C excited level. The excitation energy of the 12 C quasi-projectile can be calculated as:

E * = 3 i=1 E i - p 2 C 2m C + E th E * is presented in Figure 6 (right panel) (E th = -Q) where Q =-7.
272 MeV is the Q-value for the decay of 12 C in three α-particles. Three excited 12 C quasi-projectile levels are present for "true" αparticles (≈ 7.7, 9.6 and 10.8 MeV). We have analyzed data corresponding to the first 7.65 MeV 12 C* excited state, the first observable is the minimum relative energy of two out of the three α-particles (figure 7). The results are compatible with the ones obtained for central collisions, indicating that the decay mainly proceeds through the intermediate formation of a 8 Be. The other observable proposed in the previous section is the Dalitz-plot calculated for the αparticle energies in the 12 C frame, normalized to the total energy. The results from data and HF calculations are shown in figure 8. We would expect that the simultaneous decay should result in three equal kinetic energies, i.e. in an enhancement of the central part of the plot which is not the case of our data. On the contrary, it shows enhancements in regions where the energy of two α-particles are very close one to the other and the third α-energy is far from the two others. This can be considered as a signal of the predominance of the sequential decay, as confirmed by HF predictions (see figure 8 right panel). Other observables, proposed in [12] have been analyzed such as the mean energy < E α >, the max- imum energy normalized to the total energy max and the deviation from the average center-of-mass energy, defined as

E rms = < E 2 α > -< E α > 2 .
The results show a good agreement with HF predictions, i.e. with a sequential decay [START_REF] Morelli | [END_REF].

The 16 O+ 12 C reaction

For the new performed experiment 16 O+ 12 C at three different beam energies, namely E beam = 90, 110 and 130 MeV, the analysis has been focused up to now on the preliminary checks of collected data. In the case of complete fusion, such reactions lead to a fused 28 Si* system respectively at 55, 62 and 74 MeV excitation energy. From our previous results on the decay of excited 24 Mg we expect cluster correlations to possibly persist in the 28 Si* system, even at such high excitation energies. Extrapolating the information from Ikeda diagrams, we also expect the cluster degree of freedom be differently expressed varying the E* of the system, thus leading to different decay patterns deviating from statistical expectations. For the three different beam energies, the charge distributions corresponding to the preliminary selection of the fusion-evaporation channel out of the entire dataset is shown in figure 9, together with the α multiplicity spectra detected at GARFIELD angles in coincidence with an evaporation residue in the RCo. As expected, we note a quick increase of the alpha multiplicity when the beam energy goes up.

Conclusions and perspectives

In this work we have discussed the clustering effects in nuclear reactions involving N=Z nuclei. In particular results from the measurements with a 12 C beam on a 12 C target at 95 MeV and the 14 N+ 10 B reaction at 80 MeV beam energy are presented.

We have compared experimental data for the two reactions among them and, moreover, each experiment with the results of HF statistical calculations for the decay of the 24 Mg * source. The measured data are compatible with the expected behavior of a complete fusion-evaporation reaction, with some exception of specific channels corresponding to the emission α particles in coincidence with a Oxygen, Carbon or Neon residue. The experimental branching ratio excess for α particle emission has been quantified for both reactions, putting into evidence an effect due to the cluster nature of projectile and target in 12 C+ 12 C reaction but, at the same time, an indication of the persistence of cluster correlations also in the hot fused 24 Mg, as suggested by data from 14 N+ 10 B collisions.

To investigate the importance of discrete levels and the clustering effects in the CN decay, we have studied the decay of the well known Hoyle state in three α-particles. We have isolated this level in the 6-alpha decay channel in central collisions and, for the first time, we have compared its features with those of the same Hoyle level populated for the 12 C* quasi-projectile. It has been debated if the decay mechanism is sequential, through a 8 Be g.s. formation, or if a simultaneous decay in three α-particles could be contribute. Our results obtained comparing peripheral and central reactions with the formation of 12 C* Hoyle state in the decay chain of 24 Mg*, give no indications of deviations from the sequential decay mechanisms. Further measurements are needed, increasing the energy and the number of nucleons involved, in order to investigate if in medium effects may distort the decay of the Hoyle state. A natural extension of the same experimental scheme will be to compare the decay of a compound system formed in collisions between stable α-clustered reaction partners, or induced by an unstable light projectile on a chosen target in the future SPES facility at Legnaro National Laboratory. Together with the previously studied 24 Mg, the compound nucleus 28 Si can be formed with light radioactive beams provided by SPES.

DOI: 10 Figure 1 .

 101 Figure 1. Proton (left part) and α (right part) energy spectra in complete Z det =12 events detected in coincidence with a residue of charge Z res , indicated in each figure column. Data (black and blue for 12 C+ 12 C and 14 N+ 10 B respectively) are compared to model calculations (lines) for 12 C+ 12 C (see text for details). All distributions are normalized to unitary area.

Figure 2 .

 2 Figure 2. Q kin -values distribution when the residue of charge from 5 to 10, indicated in each panel, is detected in coincidence with the maximum possible number of α-particles. Data are from 12 C+ 12 C (black dots) and from 14 N+ 10 B (blue dots) reactions.The Q kin is normalized to the Q-value of 14 N+ 10 B data. For more details see text.

Figure 3 .

 3 Figure 3. Upper panel: Q kin distribution for the six α-particles events. The vertical lines indicate the events retained for the analysis. Lower panel: Reconstructed 12 C* excitation energy from the three α-particles with the minimum total energy in central collisions.

Figure 4 .Figure 5 .

 45 Figure 4. Central collisions: minimum relative α-α energy (MeV).The continuous lines are the predictions of HF predictions.

Figure 6 .

 6 Figure 6.Total energy spectrum for "true" α-particles (left panel). The 12 C excitation energy spectrum (right panel) is presented for the "true" α-particles selecting the12 C target in the ground state.

Figure 7 .

 7 Figure 7. Peripheral collisions: minimum relative α-α energy (MeV).The continuous lines are the predictions of HF predictions.

Figure 8 .

 8 Figure 8. Peripheral collisions: Left panel: experimental energy Dalitz plot. Right panel: HF Dalitz plot.

Figure 9 .

 9 Figure 9. Preliminary results of inclusive charge distribution (left part) and alpha particle multiplicity (right part) of events selected for the reaction16 O + 12 C at three different beam energies.

Table 1 .

 1 gs/ * indicate the16 O either in the ground or in an excited state Most probable experimental channel and its branching ratio together with the value predicted by the HF calculations. Errors on the experimental value (about 5%) take into account both statistical error and the possible particle contamination

	DOI: 10.1051/

1

O or in one of its excited bound states. The Q-values Q kin = -15.78 MeV for

  2 with less dissipative we indicate events with Q kin greater than the neutron emission threshold

	EPJ Web of Conferences 12	2	,	11002 (2016)				DOI: 10.1051/	epjconf/2016	122	11002
	CNR * 15
				-30	-20	-10	kin Q 0	MeV

  C frame, normalized to the total energy.

	EPJ Web of Conferences 12	2	,	11002 (2016)	DOI: 10.1051/	epjconf/2016	122	11002
	CNR * 15

  C* Hoyle state in the decay chain of 24 Mg*, give no indications of deviations from the sequential decay mechanisms.

	EPJ Web of Conferences 12	2	,	11002 (2016)	DOI: 10.1051/	epjconf/2016	122	11002
	CNR * 15

C* excitation energy. This is shown in figure3(upper panel). The two lowest12 C* excited states are present at 7.65 and 9.64 MeV.


Simon Jones

Mathias Legrand

Forced vibrations of a turbine blade undergoing regularized unilateral contact conditions through the wavelet balance method

Keywords: wavelet analysis, Petrov-Galerkin method, weak formulation, unilateral contact conditions, nonsmooth dynamics

The method of weighted residuals can efficiently enforce time-periodic solutions of flexible structures experiencing unilateral contact. The Harmonic Balance Method (HBM) based on Fourier expansion of the sought solution is a common formulation, though wavelet bases that can sparsely define nonsmooth solutions may be superior. This hypothesis is investigated using a full three-dimensional blade with unilateral contact conditions on a set of N c discrete contact points located at its tip. The unilateral contact conditions are first regularized and a distributional formulation in time is introduced allowing L 2 .S 1 / N trial functions to properly approximate in the time-domain the solution to the governing equations. The mixed wavelet Petrov-Galerkin solutions are found to yield consistent or better results than HBM, with higher convergence rates and seemingly more accurate contact force prediction.

INTRODUCTION

Efficiently predicting the vibratory responses of flexible structures which experience unilateral contact is becoming of high engineering importance primarily because of possible subsequent mechanical failure originated by fatigue. This type of response is increasingly common in industrial applications due to implementation of light materials and thin designs involving larger displacements together with tighter operating clearances between components. For example, consider aircraft engines where slender, twisted blades rotate at high rotational velocity within stationary casings where minimal clearance is desirable for turbine energy efficiency. Simulating contact between the blades and casing is not a trivial exercise since unilateral contact is usually described by inequalities and complementary conditions [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. In the time domain, structural displacements and velocities which satisfy these non-penetration Signorini conditions are known to respectively feature absolute continuity and bounded variation only [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: Definition and outlook[END_REF]. This implies displacements are not necessarily differentiable everywhere in the defined domain and velocities may exhibit jumps; these types of problems are generally referred to as nonsmooth. Broadly speaking, existence of solutions is still an open problem and intense research is devoted to the derivation of efficient timestepping solution methods [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF][START_REF] Leine | Dynamics and bifurcations of non-smooth mechanical systems[END_REF].

This class of unilateral problems can also be approached using periodic vibration theory. This allows the original initial-value formulations to be transformed into formulations that are periodic in time [START_REF] Laxalde | Nonlinear modal analysis of mechanical systems with frictionless contact interfaces[END_REF]. In structural mechanics, two families of numerical techniques can efficiently describe approximated periodic solutions while maintaining computational efficiency. The first is commonly known as the shooting method which consists of finding the initial conditions that realize a periodic motion; for free responses, the period of the motion can also be an unknown. The overall approach relies on time integration of the governing equations over one period using a nonlinear solver which iterates on the initial conditions, and possibly the period. The second technique is based on weighted residual formulations which are of interest in the current investigation. This method involves approximating the solution using a set of time-dependent basis functions, called trial functions, and enforcing the respective residual error to be orthogonal to an independent set of weighting functions [START_REF] Lapidus | Numerical solution of ordinary differential equations[END_REF][START_REF] Reddien | Projection methods for two-point boundary value problems[END_REF]. The orthogonality is enforced using an inner product in the sense of a scalar product adequately defined on a functional Hilbert space. Unlike the shooting method which can become numerically sensitive to possible jumps in the velocity field, weighted residual techniques directly enforce the periodicity conditions while the remaining possibly regularized unilateral contact constraints and governing local equations of motion are satisfied in a weak integral sense [START_REF] Wang | Dynamics of unsymmetric piecewise-linear/non-linear systems using finite elements in time[END_REF]. It is worth noting that the well-known Harmonic Balance Method (HBM) is a weighted residual approach where both the trial and weighting functions are Fourier series. The main goal of the current work is to explore relevant basis functions whose order of smoothness can be adapted to a particular system to attain accurate approximations and rapid convergence.

The paper is broken into five sections: Section 2 introduces the system of interest; Section 3 reviews the weighted residual formulation in time; Section 4 provides an overview of the different basis functions under investigation; Section 5 specifies the unilateral contact model; and Section 6 presents and compares results for a number of basis function combinations.

SYSTEM OF INTEREST

The current investigation is motivated by the contact interaction between blade tips and the surrounding casing in modern turbine engines. Structural contact is a consequence of improved energy efficient technology resulting in tight operating clearances between the rotor and casing; slight out-of-roundness of the casing can potentially yield severe damage [START_REF] Padova | Experimental results from controlled blade tip/shroud rubs at engine speed[END_REF][START_REF] Park | Analysis of a J69-T-25 engine turbine blade fracture[END_REF].

The model is simplified to a single blade, shown in Fig. 1; details of the model can be found in Section 5. The blade is subjected to unilateral contact by a rigid wall harmonically moving; the harmonic displacement is described by a time-dependent T -periodic function y.t/ which simulates the blade rotating within an out-of-round casing. Accordingly, the domain As a first approach, the contact interface c is simplified to N c D 7 contact points denoted x c;i 2 R 3 , i D 1; : : : N c and located on the tip of the blade, see Fig. 2: they are used to enforce unilateral contact conditions in a pointwise fashion where friction is ignored (i.e. contact is purely in the radial direction). This pointwise contact approximates the pressure distribution over the blade tip; this is deemed an acceptable simplification given the other approximations made in the current investigation.

The discrete gap which exists between the tip and the oscillating casing boundary is uniquely defined by N c gap functions as g i .x c;i ; t / D g 0;i n .u.x c;i ; t/ y.t//, i D 1; : : : ; N c where n is the normal to the surface at the considered contact point, and g 0 D .g 0;i / i D1;:::;N c , the initial gap distances. These gap distances are stored in a vector g.x c ; t / D .g i .x c;i ; t// i D1;:::;N c . When the gap closes to 0, the blade tip comes into contact with the boundary and unilateral contact conditions must be satisfied. The existence of periodic solutions of period T is assumed. The unknown displacement u.x; t/ has to satisfy the following complementary boundary value problem:

1. local equation of motion R u.x; t / div D f.x; t /; 8x 2 ; 8t

2. linear strain and linear elasticity constitutive law

" D 1 2 ru C r > u I D D W "; 8x 2 N ; 8t (2)
3. conditions of periodicity in time

boundary condition in force

n D 0; 8x 2 ; 8t (5)
6. unilateral contact conditions (to be satisfied for each contact point)

g i .
In the above equations, signifies density, f.x; t/ the potential external body forces acting on the blade, the superscript > denotes a transpose, and the dot superscript represents a temporal derivative. Structural damping is incorporated later by assuming ˇ-damping of the elastic modulus, as detailed in Section 3. It is assumed in Eq. ( 5) that the potential other external loadings on (e.g. aerodynamic or pressure) are negligible compared to the unilateral contact loads. The quantity .t / D . i .t // i D1;:::;N c , that mechanically is a contact force for node-wise contact conditions, stems from the enforced non-penetration conditions g.x c ; t/ 0 (to be read componentwise) and is necessarily non-negative (by convention). The complementarity condition i .t/g i .x c;i ; t/ D 0, i D 1; : : : ; N c states that the contact forces .t/ and the corresponding distances g.x c ; t/ separating the blade tip contact points from the oscillating boundary may not be zero at the same time. These three conditions are such that the mathematical object pairing the contact force to the displacement is not a function in the usual sense. It is also known that the displacement field may be non-differentiable at given instants within the period of motion and the velocity field may be discontinuous. This motivates the derivation of numerical techniques capable of efficiently handling this nonsmoothness. As a first approach, the unilateral contact inequalities (6) are simplified and replaced by a penalty function; an exponential spring K p is used to approximate the contact forces. For a single contact location, the gap g i .t/ in Fig. 1 separates the tip of the blade and the exponential spring. The penalty function is of the form f c;i .u.x c;i ; t // D max a c .e ˛gi .x c;i ;t/ 1/; 0 :

Specific values for the penalty function are provided with the other model parameters in Section 5.

WEIGHTED RESIDUAL FORMULATIONS IN TIME

The method of weighted residuals is a classic method of obtaining numerical solutions to boundary value problems by expanding the sought solution as a finite sequence of time-dependent basis functions in a proper functional space [START_REF] Lapidus | Numerical solution of ordinary differential equations[END_REF], commonly referred to as trial functions. The periodicity conditions result in a problem which can be formulated on a circle in time denoted S 1 where the period of the steady-state solution T has been normalized to 1 [START_REF] Berberian | Introduction to Hilbert space[END_REF]. Accordingly, these basis functions are taken from the L 2 .S 1 / N Hilbert space [START_REF] Holschneider | Wavelet analysis on the circle[END_REF]. The subsequent residual is rendered orthogonal to a set a linearly independent functions of the same space, referred to as weighting functions, through an inner product. Trial and weighting functions may contain the same basis function but not necessarily: the Galerkin method is a special case where the weighting functions are taken from the same functional space as that of the trial functions; the Petrov-Galerkin method involves the selection of weighting function comprising a basis which is independent of the trial function [START_REF] Durran | Numerical Methods for Fluid Dynamics[END_REF]. The purpose of this study is to seek various Hilbert spaces of square integrable functions defined on the circle S 1 which efficiently predict the displacement of the blade.

To solve Eqs. ( 1) through [START_REF] Beylkin | On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases[END_REF], the unknown displacement is expanded into a truncated series of N functions separated in space and time:

u.x; t / D N nod X i D1 # i .x/ Q u i .t / (8)
where Q u i .t / are the nodal displacements at the N nod nodes of the mesh shown in Fig. 2, and # i .x/, the shape functions in space within the standard Finite Element Method which is here implemented for the spatial variable using quadratic hexahedral elements [START_REF] Bathe | Finite element procedures in engineering analysis[END_REF][START_REF] Liu | Finite element method: a practical course[END_REF]. The corresponding dynamics is then reduced Figure 2. Finite element mesh of rotor blade; red dot shows unilateral contact nodes through a slightly adapted version of the Craig-Bampton method [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF][START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF]. The blade itself is a considered as the unique substructure of the system, its root defining d is clamped and the N c nodal displacements in the radial direction at the contact points (red nodes in Fig. 2) form the boundary degrees-of-freedom. The remaining internal degrees-of-freedom are synthesized. The resulting truncated set of component modes and full set of static modes, which capture the dynamics of the contact interface in our study, are then calculated. This yields the following reduced-order vector ordinary differential equation of size N in the Craig-Bampton basis:

M R Q u.t / C C P Q u.t / C K Q u.t / D f ext . Q u.t/; t/ (9)
where N N nod , together with the remaining periodicity conditions in time and the unilateral contact conditions. Here M and K are the base-constrained mass and stiffness matrices for the blade, and ˇdamping is enforced such that C D ˇK to account for structural damping. The displacement vector Q u.t / in Eq. ( 9) stores the temporal unknowns Q u i .t/, i D 1; : : : ; N . Similarly, f ext . Q u.t/; t/ stores the reduced external forcing functions, if any, as well as the reduced contact forces stemming from the regularized contact force introduced in Eq. [START_REF] Blatter | Wavelets a primer[END_REF].

The remainder of the derivation involves a weighted residual formulation in the time dimension. Three forms of the weighted residual method are discussed: the strong integral form, the weak form, and the distributional formulation. For this discussion, it is understood that Q u (the trial functions) and Q v (the weighting functions) are time-dependent vectors, i.e. Q u.t/ and Q v.t/, but for the sake of clarity the .t/ is omitted below.

Strong integral form

The standard weighted residual formulation of a differential equation is commonly termed the strong form. Taking the inner product of Eq. ( 9) with a weighting function stored column-wise in v results in the strong integral form of the equation: find

Q u 2 H 2 .S 1 / N such that 8Q v 2 L 2 .S 1 / N ; Z S 1 Q v > M R Q u C Q v > C P Q u C Q v > K Q u Q v > f ext . Q u; t/ dt D 0 (10)
This strong form of the equation is not necessarily the best framework for obtaining a solution [START_REF] Liu | Finite element method: a practical course[END_REF]; for this example involving a vector ordinary differential equation of order 2, the solution must be at least H 2 , limiting the permissible basis of trial functions .

Weak form

The respective weak form of the weighted residual statement can be obtained by performing one integration by parts over the domain S 1 for all terms containing a double time derivative in Eq. ( 10). This results in:

find Q u 2 H 1 .S 1 / N such that 8Q v 2 H 1 .S 1 / N ; Z S 1 P Q v > M P Q u C Q v > C P Q u C Q v > K Q u Q v > f ext . Q u; t/ dt D 0 (11)
The integral form of the weak formulation offers the advantage of shifting a portion of the functional smoothness requirement from the trial functions onto the weighting functions. More precisely, both the trial and weighting functions must now be H 1 . This allows the trial functions to be chosen from a wider permissible space [START_REF] Liu | Finite element method: a practical course[END_REF].

Formulation in a distributional sense

The above procedure of obtaining the solution from the weak formulation by performing integration by parts on the weighted residual statement can be extended one step further. Theoretically, this extension is not necessary since it is known that displacements should be absolutely continuous in time, though it may assist in the numerical derivations and allow very simple basis functions to be considered [START_REF] Studer | Numerics of unilateral contacts and friction. Modeling and numerical time integration in non-smooth dynamics[END_REF]. A weaker formulation is proposed by integrating again the terms involving time derivatives of the trial functions. This formulation can be understood in the sense of distributions, also known as generalized functions, ie:

find Q u 2 L 2 .S 1 / N such that 8Q v 2 H 2 .S 1 / N ; Z S 1 R Q v > M Q u P Q v > C Q u C Q v > K Q u Q v > f ext . Q u; t/ dt D 0 (12)
Here the double time differential on the field variable is transferred to the weighting function and the continuity requirement on the trial function is reduced. As discussed later, the desired displacement functions Q u can now be described using a series of constant piecewise functions for instance.

Time discretization

Let each nodal displacement Q u i .t / be approximated by a combination of M linearly independent trial functions q k .t /, k D 1; : : : ; M defined on t 2 S 1 . Similarly, let each Q v i .t/ be approximated by a combination of linearly independent weighting functions p k .t/, k D 1; : : : ; M also defined on t 2 S 1 .

A very brief description of the Sobolev spaces used herein can be found in the Appendix.

Collectively, this leads to

Q u i .t / D M X kD1 a ik q k .t / and Q v i .t/ D M X kD1 b ik p k .t/; i D 1; : : : ; N (13)
or, equivalently in a vector form:

Q u.t / D '.t /a and Q v.t / D .t/b (14)
In Eq. ( 14), '.t / and .t / are rectangular matrices of dimension N NM ; a and b are vectors of size NM 1.

Each of the strong, weak, and distributional formulations can now be discretized with proper timedependent basis functions. The overall goal of this paper is to discuss the numerical properties of certain wavelet bases to efficiently solve the problem of interest. Depending on the selected formulation, the order of derivation acting on Q u.t / and Q v.t/ (i.e. '.t/ and .t/ respectively) will affect the admissible basis functions. The arbitrariness of the weighing function in Eq. ( 10), [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF], and ( 12) is reflected by the arbitrariness of vector b and the corresponding discretized versions are: strong integral form:

Z S 1 > M R ' C > C P ' C > K' a > f ext .'a; t/ Á dt D 0 (15)
weak form:

Z S 1 P > M P ' C > C P ' C > K' a > f ext .'a; t/ Á dt D 0 (16)
distributional form:

Z S 1 R > M' P > C' C > K' a > f ext .'a; t/ Á dt D 0 (17)
The resulting nonlinear equations can generically be recast in the following form Ga f.a/ D 0

where Ga and f.a/ respectively stand as the linear internal and the nonlinear external contributions in Eq. ( 15), [START_REF] Durran | Numerical Methods for Fluid Dynamics[END_REF], and (17).

Matrix condensation

It is useful to note that after discretization in space and time, the penalty-like contact force is applied to only seven degree-of-freedoms (i.e. the radial components of the contact nodes). This allows the equations of motion to be condensed to seven active degrees-of-freedom. Consider the partitioning of Eq. ( 18) into passive p and active a degrees-of-freedom as follows [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF]:

Ä G pp G pa G ap G aa  a p a a à  0 f a .a a / à D  0 0 à : (19)
In this form f a .a a / is the displacement-dependent nonlinear contact force. The first block row of Eq. ( 19) is linear and gives

a p D G 1 pp G pa a a ( 20
)
resulting in the condensation equation

G aa G ap G 1 pp G pa a a f a .a a / D 0 (21)
or

N Ga a N f.a a / D 0 (22)
for simplicity. The square matrix and vectors of Eq. ( 22) are of size 7M . Eq. ( 22) is nonlinear and must be solved using an appropriate solver. Condensing the equations of motion reduces the number of unknowns in the nonlinear solution, hence reducing computational times and generally increasing the rate of convergence.

TRIAL AND WEIGHTING FUNCTION BASES

The selection of functional bases to be used in the above approaches is an important factor in approximation accuracy and computational efficiency [START_REF] Bathe | Finite element procedures in engineering analysis[END_REF][START_REF] Liu | Finite element method: a practical course[END_REF]. A priori selection of the optimal bases for unilateral contact problems is not a trivial matter due to potential nonsmoothness of the response. In the current investigation a number of bases are investigated to compare the quality of approximation, including Fourier functions, B-spline wavelets, Daubechies wavelets, and Haar wavelets.

Harmonic balance and Fourier series

The harmonic balance method (HBM) is a special case of the weighted residual method where the Fourier basis is used for both the trial and weighting functions [START_REF] Mickens | Comments on the method of harmonic balance[END_REF]. This technique is particularly effective when dealing with smooth nonlinear systems; convergence is often reached with very few terms. It has been used to study steady state response of turbine engine blades with friction dampers using a multiterm approximation [START_REF] Wang | Investigation of the vibration of a blade with friction damper by HBM[END_REF]. This approach was extended to unilateral contact and friction conditions [START_REF] Legrand | n-dimensional Harmonic Balance Method extended to non-explicit nonlinearities[END_REF] through an Alternating Frequency/Time domain strategy proposed by Cameron et al. [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF] and Pilipchuk [START_REF] Pilipchuk | Application of special nonsmooth temporal transformations to linear and nonlinear systems under discontinuous and impulsive excitation[END_REF]. It is worthy to note that the HBM formulations of Eqs. ( 15), [START_REF] Durran | Numerical Methods for Fluid Dynamics[END_REF], and ( 17) are identical. The Fourier basis for L 2 .S 1 / is defined as f1g [ fcos.2m t/; sin.2m t/jm 2 N g where m signifies the harmonic number of the function. The first six functions are shown in Fig. 3. Fourier basis functions feature an infinite degree of smoothness. While this property can be beneficial in some cases, it is unclear whether Fourier functions are optimal when simulating potentially nonsmooth problems. It is possible that a large number of harmonics would be required to accurately capture the nonsmooth response, or the approximation may exhibit Gibbs phenomenon at localized discontinuities in the sought displacement and velocity fields. Accordingly, other basis functions featuring a lesser degree of functional smoothness are explored.

0 1 -1 0 1 (a) constant term 0 1 -1 0 1 (b) cos.2 t/ 0 1 -1 0 1 (c) sin.2 t/ 0 1 -1 0 1 (d) cos.4 t / 0 1 -1 0 1 (e) sin.4 t/ 0 1 -1 0 1 (f) cos.6 t/

Brief review of discrete orthogonal wavelets

Discrete orthogonal wavelet families are composed of highly localized, oscillatory functions which provide a basis of L 2 .R/ and can be adapted to the periodic domain L 2 .S 1 / [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. These localized characteristics, or compact support, allow sparse representation of piecewise signals including transients and singularities [START_REF] Burrus | Introduction to wavelets and wavelet transforms: a primer[END_REF]. This makes them useful functions for use in the Galerkin approach when nonsmooth solutions are predicted [START_REF] Strang | Wavelets and filter banks[END_REF]. There are a large number of wavelet families and definitions thereof; both Mallat [START_REF] Mallat | A wavelet tour of signal processing[END_REF] and Strang [START_REF] Strang | Wavelets and filter banks[END_REF] provide excellent introductions to wavelet theory and history. Galerkin methods using appropriate discrete wavelet families as the trial functions have been shown to accurately approximate the solutions to both ordinary and partial differential equations [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF][START_REF] Beylkin | On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases[END_REF][START_REF] Pernot | A wavelet-Galerkin procedure to investigate time-periodic systems: transient vibration and stability analysis[END_REF][START_REF] Williams | Introduction to wavelets in engineering[END_REF][START_REF] Christon | The numerical performance of wavelets for PDEs: the multi-scale finite element[END_REF][START_REF] Burgos | Direct solution of differential equations using the wavelet-Galerkin method[END_REF].

The discrete wavelet family is built from scaling functions .t/ and wavelet functions .t/. These functions are analogous to the low-pass and high-pass filters of a filter bank [START_REF] Strang | Wavelets and filter banks[END_REF]; decomposition using scaling functions will give a "smoothed" approximation of the original signal, while decomposition using wavelet functions provides the details of the signal, or high-frequency content.

The exact decomposition of a continuous time signal Q u i .t/ can be written

Q u i .t / D X k g k J;k .t / C m X j DJ X `hj;` j;`. t/ (23)
and

J;k .t / D 2 J =2 .2 J t k/ (24) J;k .t / D 2 J =2 .2 J t k/ (25)
where J; k 2 Z; J is the dilation parameter (i.e. level), k is the translation parameter, and m is the maximum resolution given by the sampling rate of the function Q u i .t/. The span of the scaling functions at level J is commonly denoted V J , while the wavelet span is denoted W J . For orthogonal wavelet families, W J is the orthogonal complement to

V J in V J C1 V J C1 D V J ˚WJ : (26)
Provided the wavelet family is orthogonal [START_REF] Strang | Wavelets and filter banks[END_REF], the space of square-integrable functions on the real line L 2 .R/ can be decomposed using multiresolution analysis as a nested sequence of closed subspaces [START_REF] Mallat | A wavelet tour of signal processing[END_REF]

V 2 V 1 V 0 V 1 V 2 L 2 .R/ (27)
such that lim

J !1 V J D L 2 .R/: (28)
This implies that V J , the subspace composed of the set of scaling functions at level J , can approximate any function in L 2 .R/. A reduced orthonormal basis of L 2 .R/ is constructed by truncating the wavelet terms, resulting in

Q u i .t / D X k u k J;k .t / (29)
The accuracy of this approximation increases as the level J is increased. This property of the scaling functions makes them excellent trial and weighting functions in weighted residual methods because they can be adapted to the accuracy level required. The reduced orthonormal basis of scaling functions is used in the current investigation. This approach of increasing the level of scaling functions, rather than including wavelet functions to increase approximation accuracy, is a common approach in the literature [START_REF] Nielsen | Wavelets in Scientific Computing[END_REF] and used in the bulk of the applicable references cited in this paper.

;k .t `/ 0 Ä k Ä 2 J 1: (30)
This is equivalent to "wrapping around" the support R on S 1 through summation. The finite size of the interval results in the condition J 0

V 0 V 1 V 2 V 3 L 2 .S 1 / (31)
such that lim

J !1 V J D L 2 .S 1 /: (32)
A number of periodic discrete wavelet families exist [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The investigation considers three families to determine how they perform in unilateral, nonsmooth contact problems: B-spline, Daubechies, and Haar.

Orthogonal cubic B-spline scaling functions

The scaling function for the orthogonal cubic B-spline wavelet family built on R is given as [START_REF] Blatter | Wavelets a primer[END_REF] .

t / D X k c k B 3 .t k/ (33)
where cubic B-spline B 3 .t / can be written using the following formulas

B 3 .t / D 8 < : t 3 =6 .0 Ä t Ä 1/ t 3 =2 C 2t 2 2t C 2=3 .1 Ä t Ä 2/ t 3 =2 4t 2 C 10t 22=3 .2 Ä t Ä 3/ t 3 =6 C 2t 2 8t C 32=3 .3 Ä t Ä 4/ 0 .otherwise/: (34)
The coefficients c k can be determined using [START_REF] Blatter | Wavelets a primer[END_REF]

c k D c k D 1 Z 0 cos.k / p p 3 .cos / d .k 0/ ( 35
)
where k is an integer. The cubic polynomial p 3 is given as

p 3 .cos. // D 1 630 .cos 3 . / C 60 cos 2 . / C 297 cos. / C 272/: (36)
Numerical simulations showed that truncating the summation in Eq. ( 33) at 50 Ä k Ä 50 is sufficient; larger k terms add negligibly to the summation. After periodization to S 1 and normalization by R 1 0 J;k .t / dt D 1 [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF], sample orthogonal cubic B-spline scaling functions for J D 0; 1; : : : ; 5 are shown in Fig. 4 4.4. Daubechies Wavelets The Daubechies wavelet family is defined by a set of L filter coefficients fp `W `D 0; 1; : : : ; L 1g, where L is an even integer. The scaling function is defined by the fundamental two-scale equation [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF] . which has fundamental support over the finite intervals OE0; L 1. This equation can be used to determine the value of the scaling function at dyadic points t D n=2 J , n D 0; 1; : : : using the algorithm provided by Chen et al. [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]. The exact solution for derivatives of the Daubechies scaling functions are also available up to derivative order L=2 1.

t / D L 1 X `D0 p ` .2t `/ (37)
The wavelet filter coefficients p `were derived by Daubechies to produce scaling and wavelet functions with specific properties [START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF][START_REF] Daubechies | Ten lectures on wavelets[END_REF], some of which include: the area under the scaling function is unity

Z 1 1 .t / dt D 1 (38)
the coefficients sum to two

L 1 X `D0 p `D 2 (39)
the scaling function and its translates are orthogonal

Z 1 1 .t / .t k/ dt D ı 0;k k 2 Z (40)
The corresponding scaling functions are highly nonsmooth and fractal in nature: as one increases the resolution of the functions, the shape does not converge but rather continues to increase in complexity. This makes accurately estimating the inner products of such scaling functions with each other prone to error when numerical integration is used [START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF].

When Daubechies scaling functions are used in a Galerkin approach, it is necessary to derive the inner products of the scaling function with itself or derivatives of itself. The exact solution to these inner products can be found by using the recursive nature of the fundamental equation on L 2 .R/ [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF][START_REF] Chen | The computation of wavelet-Galerkin approximation on a bounded interval[END_REF]; the solution to these inner products are commonly referred to as connection coefficients. When Daubechies scaling functions are periodized on S 1 , the wrapping procedure results in functions which are no longer scale-invariant at low J values (scale-invariance requires wavelets at any scale to be a pure dilation of the mother-wavelet). Fig. 5 provides examples of the periodized scaling functions for L D 6. Interestingly, the lack of scale-invariance for small J values does not invalidate the connection coefficient algorithms derived for unbounded domains; the connection coefficients can simply be wrapped around the periodic domain as necessary [START_REF] Restrepo | Inner product computations using periodized Daubechies wavelets[END_REF].

Since the compact support of the father scaling function is S 1 , the periodized function is equivalent. Example Haar scaling functions for J D 0; 1; : : : ; 5, normalized by Eq [START_REF] Pilipchuk | Application of special nonsmooth temporal transformations to linear and nonlinear systems under discontinuous and impulsive excitation[END_REF], are shown in Fig. 6. The geometry [START_REF] Strazisar | Laser anemometer measurements in a transoic axial-flow fan rotor[END_REF] is commonly used for validation of aerodynamic simulation for rotating machinery codes. Limited documentation is available regarding structural analysis of the blade; Doi [START_REF] Doi | Fluid/structure coupled aeroelastic computations for transonic flows in turbomachinery[END_REF] reports the blade is manufactured from a titanium alloy with the following properties: elastic modulus of 117:2 GPa, Poisson's ratio of 0.3, density of 4539:5 kg m 3 . The mesh, shown in Fig. 2, is constructed in ANSYS using 20-node quadratic-hexahedral elements; the blade-foot is rigidly constrained. The mesh was adjusted until the first 30 eigenmodes converged to within 1 %. Seven nodes along the blade tip (roughly equally spaced) are selected as contact nodes (21 degreesof-freedom). This mesh is then further synthesized through an adapted version of the Craig-Bampton method [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF][START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF] using 20 component modes, resulting in a final model size of 41 degrees-of-freedom. Validation of this reduced model showed mode convergence of 1 % or better for the first 20 modes, thus the mesh of the finite element model is deemed acceptable for use in the unilateral contact simulation.

The contact penalty function parameters a c and ˛in Eq. ( 7) are selected to simulate contact with a rigid wall: a c D 10 4 N, ˛D 5 10 4 cm 1 . To select these values, penalty function parameters were increased until the results converged such that small increases in the penalty function parameters negligibly affected the predicted response; further increase in the parameters did not significantly decrease wall penetration and only served to diminish the conditioning of the system. The beta-damping ˇin Eq. ( 9) is set to ˇD 10 5 s, which results in light damping .< 0:5 %/. Note, the penalty function of Eq. ( 7) assumes elastic contact. It is possible incorporating a restitution coefficient to simulate the dissipation found in inelastic contact could help improve convergence in the time-marching solution [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF][START_REF] Flores | On the continuous contact force models for soft materials in multibody dynamics[END_REF].

To act as a comparison solution, the unilateral contact finite-element equations detailed above are solved using a variable-order numerical differentiation formula (NDF) time-stepping algorithm [START_REF] Shampine | The Matlab ODE Suite[END_REF]. As the time-stepping solution does not solve for the periodic response directly, the solution is deemed to have converged to its periodic state when the rms of the relative error between the tip displacements for the i and i 1 periods is below 10 5 %. For boundary oscillation frequencies between 50 1000 Hz this generally required simulation of 100 200 periods of vibration; this approach is possible because the structural damping causes the transient response to decay away, albeit slowly.

Weighted residual formulation

Table I lists the combination of trial and weighting functions employed in the current investigation; where wavelet families are listed (B-spline, Haar, Daubechies), it is implied the scaling functions are used according to Eq. ( 29). The integrals involved in Eq. ( 18) are computed as discrete inner products using 1024 points for all but the Galerkin formulation using Daubechies functions. For the Daubechies Galerkin cases, the inner products are derived using the connection coefficients as mentioned in Section 4.4.

RESULTS

The six combinations of trial and weighting functions listed in Table I are employed to solve the condensed, nonlinear weighted residual formulation given in Eq. ( 22) for prescribed boundary frequencies ranging between 50 Hz to 1000 Hz. The computed evolution of the blade tip displacement vs. boundary frequency for the center contact node is shown in Fig. 7, normalized over one period of vibration. I; the tip displacements at this frequency are still relatively large thus large contact forces were expected. This typically results in numerical convergence difficulties thus comparison at this frequency should highlight the strengths and weakness of the proposed formulations.

Tip displacement

Samples of the predicted tip displacement responses at 450 Hz using 32 (J D 5) and 128 (J D 7) basis functions are provided in Fig. 8. All results shown are for the center contact node. All contact nodes were found to interact with casing wall and produce unique responses; the center node response was deemed a representative response thus only these results are shown for brevity.

The HBM method (Fourier:Fourier) approximates the tip displacements well using both 32 and 128 basis functions compared to the time-stepping solution. This is expected as the displacement curves are relatively smooth and well suited to the HBM method. The B-spline Galerkin formulation appears to have difficulty accurately computing the response using fewer basis functions but improves as the number of functions increases. The Daubechies Galerkin formulation (DB6:DB6) appears to perform well for both cases shown, though a quantitative analysis of the results in Sec. 6.5 show it is slightly less accurate than HBM.

The Haar:Fourier and Haar:B-spline formulations result in a stepped approximation of the tip displacement due to the blocky nature of the Haar scaling functions, particularly noticeable for lower numbers of basis functions. As the number of basis functions is increased, the step-size of the Haar function within the period is reduced, and at 128 basis functions the approximation follows the time-stepping solution closely. The benefits of using Haar scaling functions as trial functions become more evident when viewing the predicted contact forces, as discussed in Section 6.3.

Tip velocity

Consideration is also given to the accuracy with which the trial functions can predict velocities. For all displacement approximations where the trial basis functions can be differentiated pointwise to approximate velocities using Eq. ( 13), the following equality is used: The only cases where Eq. ( 42) is not valid are those involving Haar trial functions; Haar scaling functions are piecewise constant functions thus they cannot be directly differentiated in the usual sense (i.e. pointwise) to determine velocities. Velocities for these cases can be approximated by assuming the derivative of a Haar scaling function is the combination of a positive delta function at t 1 and a negative delta function at t 2 , where t 1 corresponds to the positive jump discontinuity in the Haar scaling function and t 2 corresponds to the negative jump. Using a set of these functions as P q k .t/ in Eq. ( 42) provides a reasonable approximation to the corresponding velocity field for cases involving Haar functions.

P Q u i .t / D M X kD1 a ik P q k .t /: (42)
The approximate center-node tip velocity response at 450 Hz is plotted in Fig. 9 for all cases using both 32 and 128 basis functions. The delta function representation discussed above can be seen in Figs. 9(e) and 9(f) for the 32 basis function case; only discrete points at the tips of the delta functions are shown for the 128 basis function case for clarity of the figure.

It is visible from these approximations that there is a sharp change in velocity due to the contact condition about halfway through the period. This type of response is known to typically produce in ringing due to Gibbs phenomenon at the sharp gradient change using Fourier trial functions. While this phenomena was found to exist, it was slight and not readily visible on the plot. Interestingly, the ringing is more pronounced in the Galerkin B-spline case, suggesting this formulation may not be an optimal basis for this particular simulation.

The 128 function DB6:DB6 combination closely follows the time-stepping solution, more accurately predicting the sharp gradient changes which the HBM solution smoothed out. This is a result of a few DB6 scaling function characteristics: vanishing moments up to order 2 allowing accurate approximation of smooth curves; and the highly localized and fractal nature of the function which capture the nonsmooth changes in the response.

For cases involving Haar trial functions the functional velocity envelope is reasonably approximated by the delta function representation, especially as the number of basis functions is increased. They appear to do particularly well in capturing the rapid gradient changes.

Tip contact force

The contact forces at the blade tip, calculated using the penalty function provided in Eq. [START_REF] Blatter | Wavelets a primer[END_REF], are computed at the seven contact nodes simultaneously. The contact force magnitudes for the seven nodes predicted at 450 Hz using the time-stepping method are presented in Fig. 10. The contact force for the middle node predicted at 450 Hz are presented for the six combinations in Fig. 11. Notice the time duration of the contact forces coincide with the interval where the tip displacement equals the boundary displacement of Fig. 8 as this is the only period during which the blade tip is in contact with the oscillating boundary. The effect of Gibbs phenomenon can be seen in the HBM and Galerkin B-spline cases at the nonsmooth transition between a free and forced blade tip. When simulating contact using a rigid contact law, this ringing can be detrimental to numerical convergence as it artificially causes the blade tip to bounce repeatedly off the boundary. This can cause the nonlinear solver to diverge from the desired solution.

The DB6 Galerkin formulation does well in capturing both the smooth sections and the nonsmooth transitions, particularly using 128 basis functions. The use of Haar trial functions does especially well in locating the nonsmooth transition. It is hypothesized that if a rigid contact law is enforced the Daubechies or Haar scaling functions will perform better relative to the other functions; this will be investigated in future work.

Energy norm convergence

Convergence of the normalized energy norms for each basis combination at 450 Hz are given in Fig. 12(a). For these results the distributional form of the governing equations discussed in Section 3.3 are used in all cases. Results using the strong and weak forms of the governing equation are discussed further in Section 6.7.

The system energy is calculated as

E D 1 2 P Q u T M P Q u C 1 2 Q u T K Q u C E penalty : (43)
where E penalty is the energy stored in the penalty springs used to enforce contact. The scalar energy norm is calculated as the rms value of the system energy E over a single period; this is then normalized by the energy norm from the time-stepping solution. The velocity terms P Q u are calculated as detailed in Section 6.2. As shown in Fig. 12(a), all combinations of trial and weighting functions converge to the time-stepping solution (i.e. energy norms are approximately unity using 256 basis functions). The HBM, Galerkin B-spline, and B-spline:Fourier formulations show energy norms close to unity even using relatively few basis functions. The combinations utilizing Haar trial functions tend to under-predict the system energy when few trial functions are used. This is predictable considering the blocky nature the response (e.g Fig. 8(e)) which gives a coarse representation of the periodic displacement and velocity responses.

Also notice the large energy norms predicted for the DB6:DB6 formulation using 32 or few basis functions. This is due to poor representation of the response compared to the time-stepping solution. Whilst the DB6:DB6 formulation has many good attributes for solving nonsmooth contact problems, these results suggest a greater number of basis functions may be necessary to predict globally accurate solutions as compared to the HBM method (see Section 6.7 for further discussion).

Tip displacement and wall force relative error

It is also desirable to quantify the global approximation error for the blade tip displacement and resulting wall force as predicted using the weighted residual combinations. The time-stepping solution is again used as a comparison solution and the results for the center contact node are shown. The rms displacement error over one period is plotted for the six basis function combinations in Fig. 12(b), while the rms wall force error is plotted in Fig. 12(c).

As shown in Fig. 12(b), the HBM and B-spline:Fourier approaches have the lowest error when 128 or few basis functions are used. This quantifies the visual evidence reported in Section 6.1. The relative error using Haar trial functions is largest out of the combinations investigated. This is expected due to the blocky nature of the Haar scaling functions; however the mid-points of the Haar functions match very well to the trend of the time-stepping solution. This suggests that if only a low resolution approximation of the response is required, the Haar scaling functions still perform well.

The Galerkin DB6 case performs reasonably well throughout, though is always outperformed by the Fourier cases. This is attributed to the fact the displacement curves are relatively smooth thus well suited to HBM. The benefits of the nonsmooth DB6 scaling functions to accurately predict sharp gradient changes is not beneficial here. However, as shown in the force error plot Fig. 12(c), the Daubechies wavelets show promise in predicting the non-smooth contact force.

The relative errors in the predicted contact forces are presented in Fig. 12(c). The HBM and Bspline:Fourier formulations perform the best when fewer basis functions are employed, but the DB6:DB6 case produce the most accurate prediction of the contact force once 256 functions are employed. This is again attributed to the ability of the Daubechies basis to represent sharp gradient changes, such as the transition from free vibration to contact with the casing. The Haar:Fourier formulation does not perform as well, but not as poorly as expected. This is attributed to the Haar functions ability to accurately capture the flat portion of the force response, and adequately represent the zone of contact.

Relative computational cost

While each formulation requires computation of the same number of unknowns for a given basis size, the computational cost of each method was found to differ. The difference between the methods shown in Fig. 12(d) is mainly due to the number of iterations required for the non-linear solver to converge. The Galerkin B-spline and methods involving Haar functions were found to require significantly more iterations in general compared to the other methods.

The DB6 Galerkin formulation required the least amount of computing time for the non-linear solution to converge at all basis sizes. In general, this formulation required significantly fewer iterations for the non-linear solver to reach the convergence criteria. It is hypothesized that this is due to the more accurate representation of the contact forces which prevented numerical "bouncing" or "ringing".

The HBM solution also performed very well. Note, fast-Fourier transforms were not employed in this investigation to allow a more balanced comparison between the methods. If the FFT method were incorporated into the code this would significantly reduce run times for HBM cases. Fast-wavelet transforms are also available for Daubechies wavelets and will be investigated further to determine the benefits to computational time.

Principal Angles Between Subspaces

The solution predicted using a weighted residual technique must exist in the approximate solution space defined by the trial and weighting functions chosen. The Galerkin and Petrov-Galerkin approaches require the projection of the trial function basis (and/or its derivative) onto the weighting function basis (and/or its derivative); it is within this projected space that the approximate solution is determined. Increasing the dimensionality of either the trial or weighting functional bases only improves the prediction accuracy if the inner product of the bases expands the approximate solution space to encapsulate more of the L 2 .S 1 / space.

Principal angles provide a measure of the commonality or "overlap" of the trial and weighting function spaces, where Shonkwiler [START_REF] Jiang | Angles between Euclidean subspaces[END_REF] shows that there are an equal number of principal angles as trial (or weighting) functions. If the trial and weighting functions span the same space, all the principal angles will equal zero and the approximate solution space will be large. If the trial and weighting functions are orthogonal to each other, all the principal angles will equal 90 ı and the approximate solution space will be null. The more principal angles with values close to 0 implies the projection of the trial functions onto the weighting functions span a larger approximate solution space and thus can potentially provide more accurate predictions.

The principal angles for the six Galerkin and Petrov-Galerkin combinations with 64 degrees-of-freedom are provided in Fig. 13 for the three projections required in the weaker form (distributional approach) Principal Angle TF (0) → WF (0) TF (0) → WF (1) TF (0) → WF (2) TF (1) → WF (1) (a) Fourier:Fourier outlined in Section 3.3: TF .0/ ! WF .0/ ; TF .0/ ! WF .1/ ; TF .0/ ! WF .2/ . For example, TF .0/ ! WF .1/ signifies the trial function basis being projected onto the first derivative of the weighting function basis. Also included are the principal angles for TF .1/ ! WF .1/ which is required if the weak form is used. As expected the Fourier:Fourier combination has principal angles equal to zero for all but a couple of dimensions. Consider TF .0/ ! WF .2/ for example: the trial function basis contains the 0 th Fourier term (i.e. the constant offset vector), whereas in the second derivative of the weighting functions this vector is null. Thus when TF .0/ is projected onto WF .2/ the information contained in the 0 th Fourier term is lost, hence the single 90 ı principal angle value. All other basis vectors are shared by both spaces so the remaining 63 principal angles are nil as illustrated in Fig. 13(a). This property is not shared for all Galerkin bases. In Fig. 13(b) the principal angles for the B-Spline Galerkin method are shown and have a greater range than the Fourier:Fourier case. This implies the projected solution space is somewhat smaller for the B-Spline Galerkin approach. This aides in the explaining why the relative errors shown in Fig. 12(b) and 12(c) are larger for B-Spline:B-Spline than for Fourier:Fourier.

Note that for all the Galerkin combinations TF .0/ ! WF .0/ and TF .1/ ! WF .1/ the principal angles are all zero. This is because for these cases the trial and weighting function bases are the same. This implies the well-known fact that the weak formulation may provide better prediction accuracy than the distributional formulation, where TF .0/ ! WF .2/ is required. Numerical simulation verifies this hypothesis, but the improvement in accuracy is orders of magnitude less than the relative error thus results for the weak form are not included in the discussion.

For the Petrov-Galerkin cases shown in Fig. 13(d) to 13(f), all the required projections result in a reduction of the solution space. The combinations involving Haar functions as the trial basis tend to have large mean principal angles, which attests to the Haar:Fourier and Haar:B-Spline combinations performing poorly according to Fig. 12(b) and 12(c).

Fig. 13(c) shows the DB6:DB6 formulation also has a large mean for its principal angles. This aids in explaining why the numerical convergence of this formulation was poor using 64 basis functions or fewer; the approximate solution space is relatively small. As the total number of trial and weighting functions is increased, this increases the number of degrees-of-freedom with low principal angle values making possible an accurately predicted solution.

FURTHER BLADE RESULTS

While this investigation focuses on determining effective functional basis for predicting blade response under unilateral contact conditions, a few post-processed FEA results are included below to highlight the practicality of the proposed approach.

Full-field result profiles

Once the periodic displacements from the reduced FEA model are computed using the wavelet-Galerkin method, the full-field displacements of the blade can be reproduced using the adapted Craig-Bampton method [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF][START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF]; Fig. 14 shows the displacement field during the periodic response at 320 Hz. From these results, stress and strain contours can be created using standard finite element methods.

Forced response

Another result of interest is the ability to construct the frequency response function for the unilateral contact condition. The penalty function simulating contact is effectively a stiffening condition, thus the resonant peaks for the system are expected to increase in frequency relative to the no-contact response.

To capture this nonlinear response an arc-length continuation method [START_REF] Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF] is utilized. The frequency response function over one of the system's resonant frequencies is shown in Fig. 15 using both the HBM and Galerkin DB6 combinations with 256 basis functions. To compute the forced response, the position of the rigid boundary is kept constant and a harmonic point load is introduced to the contact node.

As shown, the stiffening behaviour of the unilateral contact condition is captured by both trial:weighting function combinations; other combinations showed similar approximations. These results are promising as it implies the arc-length continuation method can be accurately applied to both Galerkin and Petrov-Galerkin approaches using the distributional formulation for a number of different bases. The method of weighted residuals for capturing the periodic responses of unilateral contact problems has been investigated on an industrial blade for a number of trial:weighting function combinations. The contact condition is simulated using an exponential penalty function approach. To extend the allowable trial function bases, a weak and a distributional formulations are presented which transfer trial function continuity requirements to the weighting functions. This allowed piecewise constant Haar scaling functions to be used as a trial basis in the current investigation. Results show that a number of trial:weighting function combinations produce accurate solutions which rapidly converge as the size of the discrete spaces is increased. As expected, Fourier functions perform well as a trial basis, though nonsmooth functions such as Haar and Daubechies scaling functions are also attractive since they provide comparable prediction accuracy and even out perform the Fourier functions in some measures.

Research is continuing on this subject. The goal is to improve on the penalty approach by implementing an exact unilateral contact Lagrange multiplier method. The hypothesis is that nonsmooth trial functions will be able to approximate these flat sections better than Fourier functions, and thus could provide better contact force estimation and improved convergence for nonlinear solution algorithms.

Furthermore, the results from this investigation show that the prediction accuracy using wavelet scaling functions may increase at a greater rate than Fourier function when the global scale is increased. However, this global scale increase may not be optimal. In most cases the coarse wavelet functions can reasonably predict general trends in the response, it is only locations of sharp gradient changes where fine scale functions are required (i.e. very short compact support); this local increase in resolution is not possible with Fourier trial functions [START_REF] Li | A wavelet balance approach for steady-state analysis of nonlinear circuits[END_REF]. Future work will focus on resolving only these high gradient areas using wavelets and leaving coarse scaling functions to describe the smooth sections. Increased accuracy shall be met while minimizing the number of basis functions required. L 2 .S 1 / is the space of square integrable functions on the circle. H 1 .S 1 / is the space of square integrable functions on the circle in which the first derivative (in a distributional sense) is also square integrable. H 2 .S 1 / is the space of square integrable functions on the circle in which the first and second derivatives (in a distributional sense) are also square integrable. L 2 .S 1 / N is the space containing N square integrable functions on the circle, where N is the number of spatial degrees-of-freedom in the current investigation-this is analogous for H 1 .S 1 / N and H 2 .S 1 / N .

Figure 1 .

 1 Figure 1. Schematic of the rotor blade undergoing unilateral contact conditions. The continuous contact interface c is approximated by N c D 7 contact points.

Figure 3 .

 3 Figure 3. First six functions of the Fourier basis

Figure 4 .

 4 Figure 4. Examples of the first six levels of periodized orthogonal cubic B-spline scaling functions

Figure 5 .

 5 Figure 5. Examples of the first six levels of periodized 6 coefficient Daubechies scaling functions (DB6)

17 Figure 6 .

 176 Figure 6. Examples of the first six levels of periodized Haar scaling functions

Figure 7 .

 7 Figure 7. Evolution of the blade tip center displacement for 50 Hz to 1000 Hz boundary frequencies

Figure 8 .

 8 Figure 8. Blade tip center displacement over one period T at 450 Hz forcing frequency

Figure 9 .

 9 Figure 9. Blade tip center velocity over one period T at 450 Hz forcing frequency

Figure 10 .

 10 Figure 10. Contact force magnitude across all contact nodes (node 1: leading edge; node 7: trailing edge)

Figure 11 .

 11 Figure 11. Blade tip center contact force over one period T at 450 Hz forcing frequency

Figure 12 .

 12 Figure 12. Convergence for center contact node at 450 Hz as the number of basis functions are increased: (a) energy norm; (b) blade tip displacement error; (c) contact force error; (d) run-time

Figure 13 .

 13 Figure 13. Principal angles for Galerkin and Petrov-Galerkin subspaces with 64 basis functions. TF signifies trial function, WF signifies the weighting function, and the superscript .n/ signifies the derivative order.

Figure 14 .

 14 Figure 14. Displacement field in deformed blade during a single time-step of periodic response at 320 Hz

Figure 15 .

 15 Figure 15. Frequency response function over the 1 st fundamental frequency of the unilateral contact system

  4.2.1. Periodization of wavelet familiesStandard wavelet definitions (i.e. scaling and wavelet functions) are commonly built on the real line. The functions can be adapted to periodic functions of L 2 .S 1 / by utilizing a standard periodization technique[START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Erlebacher | Wavelets: Theory and applications[END_REF][START_REF] Meyer | Ondelettes et opérateurs[END_REF]. Let .p/ .t/ be the periodized form of the scaling function .t / defined on R

	.p/ J;k .t / D	X	J
		`2Z

Table I .

 I List of trial and weighting function combinations

	Trial function Weighting function Fourier Fourier B-spline B-spline Daubechies 6 Daubechies 6 B-spline Fourier Haar Fourier Haar B-spline

ACKNOWLEDGEMENTS

The authors would like to thank Professor Gantumur Tsogtgerel from McGill University and Professor Stéphane Junca from University of Nice Sophia-Antipolis for their insightful discussions regarding wavelets and functional analysis. Thanks also go to Alain Batailly and Markus Meingast for support concerning the Craig-Bampton reduction of the blade. Jones would also like to acknowledge the Tomlinson Postdoctoral Fellowship for partially funding this research.

A. SOBOLEV SPACES

In the paper, S 1 denotes the circle compact support that can be identified to the interval .0; 1/ on the real line. Due to the periodicity conditions in time that should be satisfied by the solution, integral forms are readily derived on this set.

Below is a list of the Sobolev spaces defined on S 1 with brief descriptions on how they relate to the current investigation: