Myriam Foucras

Jérôme Leclère

Cyril Botteron

Olivier Julien

Christophe Macabiau

Pierre-André Farine

Bertrand Ekambi

Study on the cross-correlation of GNSS signals and typical approximations

Keywords: Acquisition, Correlation, Cross ambiguity function, Doppler, Galileo, GPS

In global navigation satellite system (GNSS) receivers, the first signal processing stage is the acquisition, which consists of detecting the received GNSS signals and determining the associated code delay and Doppler frequency by means of correlations with a code and a carrier replicas. These codes, as part of the GNSS signal, were chosen to have very good correlation properties without considering the effect of a potential received Doppler frequency. In the literature, it is often admitted that the maximum GPS L1 C/A code cross-correlation is about -24 dB. We show that this maximum can be as high as -19.2 dB when considering a Doppler frequency in a typical range of [-5, 5] kHz. We also show the positive impact of the coherent integration time on the cross-correlation, and that even a satellite with Doppler outside the frequency search space of a receiver impacts the crosscorrelation. In addition, the expression of the correlation is often provided in the continuous time domain while its implementation is typically made in the discrete domain. It is then legitimate to ask the validity of this approximation. Therefore, the purpose of this research is twofold. First, we discuss typical approximations and evaluate their regions of validity. Second, we provide characteristic values such as maximums and quantiles of the auto and cross-correlation of the GPS L1 C/A and Galileo E1 OS codes in presence of Doppler, for frequency ranges up to 50 kHz, and for different integration times.

Introduction

The first stage of a global navigation satellite system (GNSS) receiver is acquisition, which consists in determining the Doppler frequency and the code delay of the received GNSS signals (Van Diggelen 2009, pp 127-224). As shown in Fig. 1, this is done by multiplying the received signal with local replicas of the carrier and of the code and integrating the result during a certain time, called coherent integration time and denoted . The output value then provides the degree of similarity between the replicas and the received signal. Since the receiver does not know a priori the Doppler frequency and the code delay, the different possibilities have to be tested. The output of the acquisition, which is a two-dimensional function of the code delay and Doppler frequency, is called the cross ambiguity function, or CAF (Ipatov 2005, pp. 7-76;Motella et al. 2010). An illustration of a CAF is given in Fig. 1, with a GPS L1 C/A signal where the incoming signal has a Doppler frequency of 2600 Hz and a code delay of 766 chips. The evaluation of the CAF can be performed sequentially by computing one point after the other, or in parallel using methods based on fast Fourier transforms (FFT) [START_REF] Akopian | Fast FFT based GPS satellite acquisition methods[END_REF]Foucras et al. 2012;Leclère et al. 2013). Quite often in the literature, approximate models are used for the expression of the CAF. However, the validity of such approximations has not been completely justified. For example, the model provided in Spilker (1996, pp. 57-120), andHolmes (2007, pp. 349-472) considers that the result of the codes correlation is independent from the carrier Doppler frequency, which is actually not the case. Indeed, the real auto and crosscorrelations of the GNSS codes can be affected by a Doppler frequency. For example, Spilker (1996, pp. 57-120) and Kaplan et al. (2005, pp. 113-152) provide the probability that the cross-correlation of C/A codes reaches a certain level, but only for Doppler frequencies multiples of 1 kHz. However, as shown later, these frequencies do not necessarily correspond to the worst cases, because they depend on the coherent integration time used. Some papers already discussed the impact of the cross-correlation with a Doppler on the carrier-tonoise ratio or on the tracking discriminator error (Raghavan et al. 1999,Van Dierendonck et al. 2002[START_REF] Balaei | Cross Correlation Impacts and Observations in GNSS Receivers[END_REF], Lestarquit and Nouvel 2012, Margaria et al. 2012). Other papers discussed the performance of several families of GNSS codes, looking at the auto and crosscorrelations affected by the carrier Doppler and by the data bit transitions (Soualle et al. 2005, Wallner et al. 2007, Soualle 2009, Qaisar and Dempster 2010). However, the impact of the Doppler frequency is essentially studied only for several frequency candidates or for a small frequency range. Qaisar and Dempster (2007) also presents some results regarding the impact of the sampling frequency and of the Doppler on the code correlation, but only for some specific values.

Therefore, what is missing in the literature is a more exhaustive study regarding the impact of the Doppler frequency for an extensive set of values and the coherent integration time used. Here, we give the auto and cross-correlation results for Doppler frequencies in different frequency ranges corresponding to specific GNSS applications, and the impact of the coherent integration time is also discussed.

The first section presents the signals definitions and the assumptions used. The second section deals with the study of the usual approximations on the CAF. The main objective is to provide the associated region of validity of the approximations. The results regarding the impact of the Doppler frequency in the CAF are presented in the third section, for the GPS L1 C/A and Galileo E1 OS codes. This allows us to compare the auto and cross-correlation properties for each signal. Finally, conclusions are provided.

Signals definitions

The GPS L1 C/A signal transmitted by a satellite can be defined as

(1) where is the time, is the amplitude of the signal, is the pseudo random noise (PRN) code, is the data, and is the carrier frequency used by the satellite (Navstar 2014). The GPS L1 C/A codes contain 1023 chips and are generated with a chipping rate of 1.023 MHz, the code length is thus 1 ms.

At the receiver antenna, the signal is the sum of signals coming from satellites, where each signal is the transmitted signal that has been respectively attenuated, delayed and affected by the Doppler effect. After the antenna, there is the front-end where the signal is amplified, filtered, down-converted to baseband or to an intermediate frequency , and sampled at a frequency . At the output of the front-end, considering a complex sampling, the signal received can be defined as

∑ ( ( )
) )

(2) where is the sampling period ( ), is the delay of the code, is the Doppler frequency, is the carrier phase, and is a white Gaussian noise of zero mean coming from the thermal noise. Note that this model does not take into account some effects, such as the Doppler effect on the code, since it can be neglected regarding the small coherent integration times considered later, or the impact of the oscillator on the sampling frequency (Leclère 2014). Note also that we consider a complex sampling because with a real sampling there would be a term with the sum of frequencies after the mixing with the local carrier that would impact the correlation result (Motella et al. 2010). Since we want to investigate specifically the impact of the Doppler frequency and of the coherent integration time, this double frequency term is not wanted. However, the results considering a real or complex sampling would be relatively similar.

In the same way, considering the Galileo E1 OS signal which has a data component and a data free component called pilot, the signal after the frontend can be written as

∑ ( ( ) ( ( )
) ) (3)
where is the secondary code on the pilot component, and and denote the data and pilot subcarriers due to CBOC modulation as defined in European Union (2015). The Galileo E1 OS codes contain 4092 chips and are generated with a chipping rate of 1.023 MHz, the codes length is thus 4 ms.

In this study, it is assumed that the navigation data and secondary code are constant during the integration. See Wallner et al. (2007) and Soualle (2009) for some results considering a bit sign transition during the integration. Also, the filtering and quantization performed by the frontend are not taken into account. See [START_REF] Borre | Reducing Front-End Bandwidth May Improve Digital GNSS Receiver Performance[END_REF] for more details about that topic. Finally, the noise term is not considered since the focus of this research is on the correlation properties.

Study of usual approximations

In this section, we perform a study on two typical approximations, the first uses a continuous time domain model although actual processing is performed in discrete time domain, and the second models the CAF assuming that the code and Doppler are independent. For simplicity, we develop the results for the GPS L1 C/A signal, i.e.

, since they can be easily extended to the Galileo E1 OS signal and other GNSS signals.

Approximation continuous/discrete

In the literature, e.g. Kaplan et al. (2005, pp. 113-152) and Holmes (2007, pp. 349-472), the output of the CAF is often modeled by the following continuous time version:

( ̂ ̂ ) (4) |∫ ̂ ̂ | | ∫ ̂ ( ( ̂ ) ) |
where the superscript stands for continuous, denotes the L1 C/A code number of the code received, denotes the L1 C/A code number of the code replica, ̂ is the local estimate of the code delay, and ̂ is the local estimate of the carrier frequency. This expression is an approximation of what is really performed in a GNSS receiver since the signals are actually discrete. Thus, the actual output of the CAF is

( ̂ ̂ ) (5) |∑ ̂ ̂ | | ∑ ̂ ( ( ̂ ) ) |
where the superscript stands for discrete, and is the number of samples during the coherent integration time, i.e. . Remember that (5) is still an approximation, since it does not take into account the filtering of the front-end and the other factors mentioned in the previous section.

Therefore, the question is "are the results obtained with (4) close to those obtained with (5) ?", which is intuitively believed to be the case.

To check this, both equations have been evaluated for all the GPS L1 C/A codes, with a step of 10 Hz for ̂ over a range of 50 kHz and a step of 1 chip for , and considering one sample per chip. The details of how to compute exactly the integral of (4) are given in the Appendix.

Fig. 2 shows the difference between ( ̂ ̂ ) and ( ̂ ̂ ) as function of the code delay and Doppler frequency, and shows the distribution of this difference. It can be seen that the higher the Doppler frequency, the higher is the difference between the two results, but this difference is clearly negligible since the maximum difference is 0.0342 dB for a Doppler of 50 kHz. Therefore, we can conclude that using (4) as approximation of ( 5) is a valid hypothesis. ̂ is the difference of carrier frequencies, and is the cross-correlation between the codes and , defined as

∫ (7) If ,
is the autocorrelation of the code . The approximation given by ( 6) is the result of the following operation:

̅̅̅̅ ( ̂ ̂ ) (8) |∫ ̂ ∫ ( ( ̂ ) ) |
which is clearly not equal to (4). In Holmes (2007, pp. 349-472), it is mentioned that in (6) has been factored out of the integral as an approximation assuming that is small compared to the chip rate. Therefore, the question is "when this approximation is valid, i.e. when the results obtained with (6) are equal or very close to those obtained with (4) ?". First, it can be easily checked that if , i.e. ̂ , or if , i.e. ̂ , equation ( 6) is identical to (4), i.e.

̅̅̅̅ ( ̂ ) ( ̂ ) (9) ̅̅̅̅ ( ̂ ) ( ̂ ) (10
) This is illustrated by Fig. 3with ms. As a note, the CAF is symmetric with respect to , but not to . Note that the model used does not consider any filtering; in practice the bandwidth is limited and the correlation peak is rounded. chip, the approximation is very close to the exact value when kHz, i.e. for the main lobe of the sinc. Outside the main lobe, the approximation differs from the exact value. The higher | |, the higher is the difference. However, when | | chip, there are two cases as shown in Fig. 5 and Fig. 6. If the code autocorrelation without Doppler is very low, i.e. near -60 dB, a small error in the Doppler such as dozens of hertz will imply few dB of difference rapidly. However, if the code autocorrelation without Doppler is not too low, i.e. near -24 dB, the difference will be seen for a higher error in the Doppler, from hundreds of hertz, i.e. for frequencies within the main lobe. In order to have a global representation of the previous discussions, Fig. 8 shows in blue the areas where the difference between the exact computation and the approximation is lower than 1 dB. In conclusion, it has been verified that the continuous-time expression of the CAF is a very good approximation of the discrete-time one. Furthermore, it has been observed that the CAF approximation that considers the code and Doppler independently given by ( 6) is:

 Exact for a right estimation of the code delay or of the Doppler frequency (Fig. 3).

 Quite exact for a code delay less than 1 chip and a Doppler frequency less than 1 kHz, the inverse of the code period (Fig. 4), corresponding to the main peak.

 Quite exact for a code delay leading to a not too low correlation value such as -24 dB and a Doppler frequency less than about 200 Hz (Fig. 5 top and Fig. 6).

 Not valid for a code delay leading to not too low correlation value such as -24 dB and a Doppler frequency higher than about 200 Hz (Fig. 5 top and Fig. 6)  Not valid for a code delay leading to a very low correlation value such as -60 dB as soon as there is a very small Doppler frequency (Fig. 5 middle and bottom and Fig. 6).

Fig. 8 Difference between exact computation and approximation of the CAF for the L1 C/A code 1. Blue dots correspond to a not null difference smaller than 1 dB.

Evaluation of the correlation properties

In the previous section, it has been shown that the usual approximation that considers separately the code correlation and the sinc function is not representative of reality. Therefore, it is interesting to evaluate the level of the CAF for various Doppler frequencies. In this section, we provide statistics related to the values of the CAF for the auto and cross-correlation for the GPS L1 C/A and Galileo E1 OS codes; in particular, we give maximum values, distribution, i.e. probability of occurrence, and cumulative distribution through some quantiles.

Simulation scheme

The presented results come from exhaustive simulation runs for all the codes or couples of codes, for all code delays multiple of 1 chip, and for Doppler frequencies within different ranges using a step of 10 Hz. In order to speed up the computation time, the CAFs are computed following gives exactly the same correlation result as with frequencies of 2500 and 1000 Hz respectively or with frequencies of 1500 and 3000 Hz respectively. This is why Fig. 10 shows only the positive frequencies on the x-axis. 18 dB, i.e. there is only one point of the CAF over 10 that has a higher value. In the same way, still for the band [-5; 5] kHz, looking at the quantile 99.9 % only one point over 1000 has a value higher than -21.5 dB; and looking at the quantile 100 % the maximum CAF value is -19.18 dB. As it can be observed, for all the frequency bands, there is at least 1 % of the points (99%) of the CAF that are above the 0 Hz Doppler worst case of -23.94 dB.

( ̂ ) | ( ( ̂ ) )| (

Fig. 11

Level of CAF values excluding the correct alignment according to the selected quantiles and for different frequency bands, considering the 32 L1 C/A codes and an integration time of 1 ms

The same study has been done considering longer coherent integration times. In Fig. 12, we report the maximum of the CAF in the band [-5; 5] kHz for different coherent integration times (1, 2, 4, 10 and 20 ms). It can be seen that with longer integration times, the values on multiples of 1 kHz are the same while the values on not multiples of 1 kHz decrease. Therefore, when the coherent integration time is longer than the code period, the CAF maximums are on multiples of 1 kHz. This can be intuitively explained looking at the spectrum of L1 C/A codes. Since the C/A codes are periodic with a period of 1 ms, their spectrum is discrete and the spectral lines are spaced by 1 kHz (Spilker 1996, pp. 57-120), as illustrated in Fig. 13 (left) which zooms around 0 Hz, and where the height of each line depends on the code. The product of two L1 C/A codes is another L1 C/A code of the same period (Spilker 1996, pp. 57-120) therefore its spectrum is similar. A C/A code with a carrier Doppler has a spectrum almost similar to the one of a C/A code except that it is shifted, as illustrated Fig. 13 (right). Therefore, the product of a received C/A code with a carrier Doppler by a local C/A code will have a spectrum as shown in Fig. 13 (bottom). Such signal is then integrated; the transfer function of the integration is a sinc function as illustrated in Fig. 13 (bottom). Therefore, considering first an integration time of 1 ms, if the incoming signal has no Doppler, we have the spectrum shown in Fig. 14 (top left), where it can be seen that only the line at 0 Hz is kept after the integration. If the incoming signal has a Doppler of 500 Hz, we have the spectrum shown in Fig. 14 (middle left), and all the spectrum lines will be attenuated but conserved. If the incoming signal has a Doppler of 1000 Hz, we have the spectrum shown inFig. 14 (bottom left), and again only the line at 0 Hz is kept after the integration. Now, if we consider an integration time of 4 ms, the width of the sinc in the transfer function will be four times smaller. If the incoming signal has a Doppler that is a multiple of 1 kHz, see Fig. 14 (top and bottom right), only the line at 0 Hz is kept after the integration, which gives exactly the same result as with an integration of 1 ms. This explains why the correlation value is the same for Doppler that are multiple of 1 kHz whatever the integration time is. But now, with a Doppler that is a multiple of , all the spectrum lines will be cancelled as shown in Fig. 14 (middle right). This explains the nulls that appear in Fig. 12 when the integration time is increased. Note that these nulls can also be found analytically using ( 17), given in the Appendix. Finally, for other Doppler shifts all the spectrum lines will be attenuated and conserved, but the attenuation is much more compared to an integration time of 1 ms, which explains the lower CAF values.

f (kHz) 1 2 0 1 2 0 T T f (kHz) f (kHz) 1 2 0
Fig. 13 Spectrum of a C/A code without carrier Doppler (left), spectrum of a C/A code with a carrier Doppler (right), transfer function of an integration of duration ms (bottom) and 1 kHz (bottom), and transfer function using 1 ms of integration time (left) and using 4 ms of integration time (right)

1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 Doppler of 0 Hz f (kHz) f (kHz) f (kHz) f (kHz) f (kHz) f (kHz)
These observations are completed with the distributions of the CAF presented in Fig. 15. For example, with an integration time of 1 ms, 90 % of the values are below -26.35 dB, as already shown in Fig. 11, whereas with an integration time of 20 ms, 90 % of the values are below -44.56 dB (more than 18 dB). In conclusion, the use of a coherent integration time longer than the code period is much better in terms of cross-correlation protection.

Galileo E1 OS

The same study is now presented for the Galileo E1 OS signal. In this case, only the first 3 Doppler frequency ranges are considered. Indeed, because of the higher number of codes and number of code delays, to be exhaustive the quantiles should be computed on several tens of billions of points for a Doppler range of [-50; 50] kHz, which is not possible due to Matlab memory limits (using Matlab R2014a version).

Unlike GPS L1 C/A codes, for Galileo E1 OS codes, when not considering Doppler frequency, the number of possible autocorrelation or crosscorrelation values is not limited to 4 values. However, thanks to their longer lengths, 4 ms instead of 1 ms, the maximum autocorrelation value excluding the correct alignment of Galileo E1 codes is -25.39 dB without Doppler, which is lower than the maximum of -23.94 dB with the GPS L1 C/A codes. Nonetheless, with a Doppler up to 10 kHz, the maximum can be up to -23.45 dB.

Fig. 16 shows the values for some probabilities of occurrence for the different frequency bands and for the data and pilot codes. It can be observed that the autocorrelation properties of the Galileo E1 OS codes assigned to the data component are equivalent to those of the codes assigned to the pilot component in terms of distribution, and that the autocorrelation without any Doppler has indeed a multitude of different levels. It can also be observed that there is a global improvement of 2 to 8 dB compared to the GPS L1 codes, thanks to a code length that is four times longer. Foucras 2015). In this study, it is assumed that there is no bit sign transition.

The maximum of the CAF is shown in Fig. 17 for the band [-5; 5] kHz and different integration times. It can be seen that the values on the multiples of 250 Hz are the same whatever is the integration time, which was expected given the 4 ms period of the primary codes. Concluding this section on the autocorrelation, we present Table 1, which allows us to compare the autocorrelation values for the presented quantiles. When comparing the GPS L1 C/A and Galileo E1 OS code isolations, it can be observed that Galileo E1 OS codes have better correlation performance for an integration time equal to the spreading code period, with an improvement between 4 and 6 dB. Even if we compare both code families with the same integration time of 4 ms, the maximum is at least 2 dB lower for Galileo E1 OS codes. However, for a coherent integration time of 20 ms, which corresponds to 20 L1 C/A code periods or 5 Galileo E1 OS code periods, the lowest quantiles (90 % and 95 %) are lower for GPS L1 C/A. But things are different for the other quantiles (99 % and above) for which Galileo E1 codes perform better by at least 3 dB. This can be understood looking at Fig. 12 and Fig. 17, where we can see that there are more high lobes with the Galileo E1 OS signal, but the local maxima are still below the local maxima of the GPS L1 C/A signal. Note that only the first 6 Galileo E1 OS data codes have been used due to memory constraints. It is important to note that these results do not take into account the square sub-carrier modulation of BOC modulation. Cross-correlation In this section, we present statistics related to the values of the CAF for the cross-correlation for the GPS L1 C/A and Galileo E1 OS codes. As it will be shown, the results are relatively similar to those of the autocorrelation.

GPS L1 C/A

In presence of Doppler, the cross-correlation is similar to the autocorrelation, i.e. for an integration time of 1 ms, the maximum crosscorrelation value can be up to -19.1 dB, as shown in Fig. 18. The statistics for a Doppler range of [-10; 10] kHz are provided in Table 2. For longer integration times, a behavior similar to the one shown in Fig. 12 will be observed, and the maximum is obtained for a Doppler that is a multiple of 1 kHz. Van Diggelen (2009, p.219) states regarding detection problem due to cross-correlations that "both the code delay and the frequency offset of the incorrect (strong) satellite would have to be in the same search zone as the intended (weak) satellite". However, this statement is not correct, because as demonstrated in the Appendix, what matters is the difference between the received and local frequencies For example, let us consider a receiver with a frequency search space of ±500 Hz, which receives two GPS L1 C/A signals, one weak with a Doppler of 300 Hz (PRN 1) and one strong with a Doppler of 4100 Hz (PRN 2).

When the receiver will search PRN 1 and compute the correlation with a local Doppler of 100 Hz, there will be a high cross-correlation with the incoming PRN 2 signal, because the difference between its frequency and the local one is 4 kHz, which is a multiple of 1 kHz. 2 allows us to compare the maximum of each cross-correlation. The best isolation is for a code couple assigned to the same satellite. Indeed, they are chosen as orthogonal as possible to not interfere for the correlation of the local component with the received Galileo E1 OS signal.

Conclusion

We have presented various results regarding the correlations of the GPS L1 C/A and Galileo E1 OS codes. Some results are relatively well-known by people working in the GNSS field but have been clarified, while some of them are less known. In a first part, we have studied the validity of some approximations generally used. First, it was shown that the discrete summations done in GNSS receivers due to the sampling can be approximated by a continuous integration, as intuitively expected. The second point was about the classical correlator output expression, which considers the code and the Doppler frequency as independent. The region of validity of this expression has been clarified, and it is not only around the main peak as it is often thought, but also for relatively low Doppler frequencies, such as one or two hundreds of hertz with the GPS L1 C/A signal, when the theoretical code correlation is not extremely low, e.g. near -24 dB.

In a second part, we have provided a statistical description of the autocorrelation and crosscorrelation for GPS L1 C/A and Galileo E1 OS codes, considering different carrier frequency bands. For example, it has been shown that Galileo E1 OS codes have a lower crosscorrelation than the GPS L1 C/A codes for any Doppler frequency between -10 and 10 kHz considering an integration time of one code period, namely -23.44 dB against -19.18 dB. When considering all potential Doppler frequencies, comparing to no Doppler frequency the isolation is degraded by a maximum of about 2 dB and 4.75 dB for Galileo E1 OS and GPS L1 C/A respectively. The same trends were shown for the cross-correlation. For GPS L1 C/A, the isolation is ensured to be at the minimum -19.08 dB for any Doppler frequency between -10 and 10 kHz, which is 2 dB lower than the classical value given in the literature which is only valid for frequencies that are multiple of 1 kHz. For Galileo E1 OS, the isolation between data and pilot codes affected to the same satellite is -22.82 dB, which is 2.34 dB higher than any couple of Galileo E1 OS codes. We have also highlighted the significant impact of the coherent integration time, and shown that having a coherent integration time longer than the code period provides much better correlation performance. We have also shown that even the satellites which are outside the frequency search space impact the correlation result, since what matters is the difference between the frequencies, and not the frequencies themselves.

This study could be extended in many ways. First, other GNSS signals could be considered. Then, it could be extended by considering different sampling frequencies, since the sampling frequency also changes the code correlation properties. Finally, in the same way, the impact of the code Doppler could be studied, especially with signals having a chipping rate of 10.23 MHz or having pilot channels which enable long coherent integration. The remaining integral is the cross ambiguity function for a coherent integration over one code period, and is given as

∫ (16) ∫ ∑ ∑ [ ] ( ) ∫ ∑ [ ] ∑ ∫ [ ] ∑ [ ] ∫ ∑ [ ] [ ] ∑ [ ] ∑ [ ] ∑ [ ] ( ) ∑ [ ] ( ) ∑ [ ]
Therefore, the magnitude of the cross ambiguity function is which can be computed numerically exactly. In the same way, it is possible to numerically compute the cross ambiguity function when the delay between the codes is not a multiple of one chip.

From ( 14) and ( 15), it can be seen that the magnitude of the CAF does not depend on the starting point of the integration. Therefore, performing a non-coherent integration, i.e. averaging the magnitude or the magnitude squared of consecutive CAF results, will give the same result as the magnitude of one CAF result. Consequently, the distribution of the CAF magnitude does not change with non-coherent integrations.

Fig. 1

 1 Fig. 1 Principle of the acquisition (top), Illustration of the cross ambiguity function for a GPS L1 C/A signal (bottom)

Fig. 2

 2 Fig. 2 Difference between the CAF in the discrete and continuous cases, i.e. ( ̂ ̂ ) ( ̂ ̂ ), assuming receiving and generating L1 C/A code 1

Fig. 3

 3 Fig. 3 Comparison of the CAF of L1 C/A code 1 computed using (4) and (6) when chip (top), and Hz (bottom)

Fig. 4 Fig. 5

 45 Fig. 4 Comparison of the CAF of L1 C/A code 1 computed using (4) and (6) when | | chip. chip (top), chip (middle), chip (bottom)

Fig. 6

 6 Fig. 6 Comparison of the CAF of L1 C/A code 1 computed using (4) and (6) when | | Hz. Hz (top), Hz (middle), Hz (bottom)

Fig. 9

 9 Fig. 9 Simulation scheme

Fig. 10

 10 Fig. 10 Maximum of the CAF excluding the correct alignment as a function of the Doppler frequency, for an integration time of 1 ms, i.e. one code period. Each color represents a different L1 C/A code. The horizontal blue dotted line represents the CAF for a 0 Hz Doppler, i.e. -23.94 dB. The horizontal magenta line represents the maximum CAF for a not null Doppler, i.e. -19.18 dB.

Fig. 12

 12 Fig. 12 Maximum of the CAF excluding the correct alignment of the GPS L1 C/A code 1 as a function of the Doppler frequency, for integration times of 1, 2, 4, 10 and 20 ms, respectively

Fig. 14

 14 Fig.14Spectrum of a C/A code shifted by a Doppler frequency of 0 Hz (top), 500 Hz (middle) and 1 kHz (bottom), and transfer function using 1 ms of integration time (left) and using 4 ms of integration time (right)

Fig. 15

 15 Fig. 15 Cumulative distribution of the CAF excluding the correct alignment of L1 C/A code 1, for different coherent integration times, and a search space of [-10; 10] kHz

Fig. 16

 16 Fig. 16 Level of CAF values excluding the correct alignment according to the selected quantiles and for different frequency bands, considering the 100 Galileo E1 OS codes and an integration time of 4 ms. Galileo E1 OS data codes (top), Galileo E1 OS pilot codes (bottom)

Fig. 17

 17 Fig. 17 Maximum of the CAF excluding the correct alignment of the Galileo E1 OS data code 1 as function of the Doppler frequency, for integration times of 4, 8, 12, 16 and 20 ms, respectively

Table 1

 1 Levels of CAF values in dB according to the distribution considering [-10; 10] kHz as range for the Doppler frequency. For the integration time of 20 ms, only the first 6 Galileo E1 OS data codes have been used due to memory constraints.

Fig. 18

 18 Fig. 18 Distribution of the maximum of the crosscorrelation function per GPS L1 C/A codes couple considering [-10; 10] kHz as range for the Doppler frequency

Fig. 19

 19 Fig. 19 Level of CAF values according to the selected quantiles and for different frequency bands, considering the 50 Galileo E1 OS code data/pilot couples of each satellite and an integration time of 4 ms

Fig. 20

 20 Fig. 20 Illustration of a PRN code

Table 2

 2 Levels of CAF values in dB according to the distribution considering [-10; 10] kHz as range for the Doppler frequency. Blank elements for Galileo E1 OS in the table were not computed due to the high number of combinations.

				Galileo E1 OS (4 ms)
		GPS
		L1 C/A (1 ms)	Data/Pilot	Data/Data	Pilot/Pilot	Same	satellite
	90 %	-26.2
	95 %	-24.9
	99 %	-23.5
	99.9 %	-21.6
	99.99 %	-20.5
	99.999 %	-19.8
	99.9999 %	-19.3
	Max	-19.1	-22.8 -22.6 -23.0 -23.6

Acknowledgments

The authors would like to thank Chris J. Hegarty, Jack. K. Holmes, Francis Soualle and Phillip W. Ward, for their valuable comments and advice.

European Union (2015) European GNSS (Galileo) Open Service Signal. In Space Interface Control Document (OS SIS ICD) Issue 1.2. Foucras M (2015) Performance Analysis of Modernized GNSS Signal Acquisition. Ph.D. thesis, Institut National Polytechnique de Toulouse Foucras M, Julien O, Macabiau C, Ekambi B (2012) A Novel Computationally Efficient Galileo E1 OS Acquisition Method for GNSS Software Receiver. Proc. ION GNSS 2012, Institute of Navigation. Nashville, TN (USA), September, pp 365 -383 Holmes JK (2007) Spread Spectrum Systems for GNSS and Wireless Communications. Artech House Ipatov VP (2005) Spread Spectrum and CDMA: Principles and Applications. Wiley Kaplan ED, Hegarty CJ, Ward P (2005) Understanding GPS: Principles and Applications, 2nd edn. Artech House Leclère J (2014) Resource-efficient Parallel Acquisition Architectures for Modernized GNSS Signals. Ph.D. thesis, École Polytechnique Fédérale de Lausanne (EPFL) Leclère J, Botteron C, Farine P-A (2013) Comparison Framework of FPGA-based GNSS Signals Acquisition Architectures. IEEE Trans Aerosp Electron Syst 49:1497 -1518. Lestarquit L, Nouvel O (2012) Determining and measuring the true impact of C/A code cross-correlation on tracking-Application to SBAS georanging. In: Position Location and Navigation Symposium (PLANS), 2012 IEEE/ION. Myrtle Beach, SC, USA, pp 1134 -1140 Margaria D, Motella B, Dovis F On the Impact of Channel Cross-Correlations in High-Sensitivity Receivers for Galileo E1 OS and GPS L1C Signals, International Journal of Navigation and Observation, vol. 2012, Article ID 132078 Motella B, Lo Presti L (2010) The Math of Ambiguity: What is the acquisition ambiguity function and how is it expressed mathematically? Inside GNSS 5:20 -28. Navstar (2014) GPS Space Segment/Navigation User Interfaces (IS-GPS-200H). Qaisar SU, Dempster AG (2010) Cross-correlation performance assessment of global positioning system (GPS) L1 and L2 civil codes for signal acquisition. Radar Sonar Navig IET 5:195 -203. Qaisar SU, Dempster AG (2007) An Analysis of L1-C/A Cross Correlation & Acquisition Effort in Weak Signal Environments. In: Proceedings of International Global Navigation Satellite Systems Society, IGNSS Symposium. Sydney, Australia Raghavan SH, Kumar R, Lazar S, et al (1999) The CDMA Limits of C/A Codes in GPS Applications -Analysis and Laboratory Test Results. Proc. ION GPS 1999, Institute of Navigation, Nashville, TN, USA, September, pp 569 -580. Soualle F (2009) Correlation and Randomness Properties of the Spreading Coding Families for the current and Fututre GNSSs. In: Fourth European Workshop on GNSS Signals and Signal Processing. Oberpfaffenhofen, Germany Soualle F, Soellner M, Wallner S, et al (2005) Spreading Code Selection Criteria for the future GNSS Galileo. In: Proc. of the European Navigation Conf. GNSS (2005). Munich, Germany Spilker JJ (1996) Gold Code Cross-Correlation Properties. In: GPS Signal Structure and Theoretical Performance. American Institute of Aeronautics and Astronautics Tsui JB-Y (2005) Fundamentals of Global Positioning System Receivers: A Software Approach, 2nd edn. Wiley Series in Microwave and Optical Engineering Van Dierendonck AJ, Erlandson RJ, McGraw GA, Coker R (2002) Determination of C/A Code Self-Interference Using Cross-Correlation Simulations and Receiver Bench Tests. Proc. ION GPS 2002, Institute of Navigation, Portland, OR, USA, September, pp 630 -642 Van Diggelen FST (2009) A-GPS: Assisted GPS, GNSS, and SBAS. Artech House Wallner S, Avila-Rodriguez J-A, Hein GW, Rushanan JJ (2007) Galileo E1 OS and GPS L1C Pseudo Random Noise Codes -

Appendix: Computation of the continuous-time CAF

This appendix shows how to compute numerically the continuous-time expression given by (4). A PRN code is a succession of chips of duration , whose value can be 1 or -1, that repeats each chips; thus the duration of one code period is . A PRN code is illustrated in Fig. 20

Since the code is periodic with a period , ( ) ( ). Therefore,


Youssef K Hamidi

Gorkem E Guloglu

Jacob P Anderson

M Cengiz Altan

Modeling and Experimental Validation of Anomalous Moisture Absorption of Micro and Nanocomposite Laminates

Keywords: Composites, Moisture Absorption, Anomalous Diffusion, Hindered Diffusion Model, Anisotropy

published or not. The documents may come

Modeling and Experimental Validation of Anomalous Moisture Absorption of Micro and Nanocomposite Laminates

Youssef

  K Hamidi, Gorkem E Guloglu, Jacob P Anderson, M Cengiz Altan To cite this version: Youssef K Hamidi, Gorkem E Guloglu, Jacob P Anderson, M Cengiz Altan. Modeling and Experimental Validation of Anomalous Moisture Absorption of Micro and Nanocomposite Laminates. 2016. hal-01354014


Roques Sébastien
email: sebastien.roques@ec-nantes.fr

Gem

Carlo M Stoisser
email: carlo.stoisser@edf.fr

Patrice Cartraud
email: patrice.cartraud@ec-nantes.fr

Mathias Legrand
email: mathias.legrand@ec-nantes.fr

Bernard Peseux
email: bernard.peseux@ec-nantes.fr

Christophe Pierre
email: christophe.pierre@mcgill.ca

Modelling of rotor speed transient with rotor-to stator contact

Keywords: Speed transient, heavy rotor-to-stator contact, rubbing, Lagrange constraints

A rotor-stator model of a turbogenerator has been developed to solve rotational velocity transients: to this end, the angular position of a beam cross section is considered as a variable, allowing for the angular deceleration due to rotorto-stator friction to be calculated in a more realistic representation of the interaction. The equations of motion derived from an energetic formulation are highly nonlinear and require a specific time integration procedure. Calculations were run considering different types of diaphragms.

NOMENCLATURE

D R A F T - → x - → y - → z u v w - → X - → X - → Y - → Y - → Z - → Z -→ x1 - → y1 - → z1 -→ x2 - → y2

INTRODUCTION

In nuclear power plant turbosets, the design basis accident consists of the loss of a blade on the low pressure turbine last stage. During the accidental shut down, a severe rotor-casing interaction may occur at critical speeds due to large shaft line displacements originated by a high unbalance excitation. The contact between the shaft and the diaphragm induces an important deceleration rate and modifies the turbogenerator dynamics including the amplitude of the loads in the bearings. Risks of failure of the contacting components may also occur due to the heavy friction torque. Therefore the main objective is to verify that the designed turbine is capable of going through critical speeds without catastrophic consequences for the shaft line.

In a general manner, contact dynamic is highly nonlinear and many studies have been carried out to understand the rubbing phenomenon in order to prevent serious damages. Most of papers usually make the assumption that the angular velocity is known, even during contact [START_REF] Muszynska | Rotordynamics. Taylor & Francis edition[END_REF]. As far as we know, a few models have considered the angular position as a variable and the rotor was assumed to be rigid ( [START_REF] Dai | Dynamic behavior of the full rotor/stop rubbing: numerical simulation and experimental verification[END_REF] and [START_REF] Bartha | Dry friction backward whirl of rotors[END_REF]).

In this paper, only the shaft-diaphragm interaction is considered: the blades-to-diaphragm contact [START_REF] Legrand | Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion[END_REF] is assumed to be negligible in comparison with the shaft-stator interaction to compute the angular deceleration due to friction. Moreover the nonlinear behavior of the bearings is linearized. The Lagrange multipliers method in conjunction with an explicit time-stepping procedure is used to solve equations of motion [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF].

ROTOR SPEED TRANSIENT MODELLING 2.1 Shaft line components

A 1-D model has been developed considering the angular position of a cross section as a variable. Flexural vibrations in two orthogonal planes, torsional vibrations of a straight shaft and axial vibrations are considered. The following assumptions are used throughout the paper:

1. The shaft has a uniform circular cross section along its length. 2. Timoshenko theory is used. 3. The shaft is balanced (no mass eccentricity). 4. Gyroscopic terms are included. 5. External torques with constant direction along the undeformed centroidal axis are applied at each extremity of the shaft.

Rigid disk

Considering a rigid circular disk with an unknown angular speed, and due to its geometric symmetry, the kinetic energy E D c takes the form:

2E D c = M D u2 + v2 + ẇ2 + I DX θ2 u cos 2 θ v + θ2 v + I DZ θ2 u sin 2 θ v + φ + β 2 + 2 φ + β θu sin θ v (1)
D R A F T where M D ,I DX ,I DZ respectively stand for the mass, the second moment of inertia and the polar moment of inertia of the disk. Assuming small rotations and small displacements, the kinetic energy is rewritten as follows:

2E D c = M D u2 + v2 + ẇ2 + I DX θ2 u + θ2 v + I DZ φ + β 2 + 2 φ + β θu θ v (2)
This formula is exact up to the third order and equations of motion will be written up to the second order.

Shaft

The kinetic energy E c of a rotating beam is written by considering that the energy of the shaft corresponds to the integration of the disk energy along the longitudinal direction:

2E c = ρ l 0 S u2 + v2 + ẇ2 + I x θu 2 + θv 2 + I p φ + β 2 + 2 φ + β θu θ v dz (3)
The potential energy E d of a spinning Timoshenko beam is equal to:

2E d = l 0 ESw 2 ,z + EI x θ 2 u,z + θ 2 v,z + k x GS (u ,z -θ v ) 2 + (v ,z + θ u ) 2 dz (4)
Bending and torsion are coupled together because of gyroscopic terms and to the unknown angular position. The equations of motion derived from Hamilton's principle are nonlinear even without contact.

Linearized bearing

The dynamical study of a shaft line requires to describe the nonlinear behavior of the oil film. Assuming small displacements of the rotor, stiffness and damping coefficients of the oil film can be calculated by linearizing Reynolds' equations with respect to the equilibrium position. Calculating the virtual work δW of external forces acting on the shaft yields:

δW = F u F v • δu δv
where F u and F v are the components of the generalized forces acting in the bearings. After linearizing, they can be written in a matrix form:

F u F v = - k xx k xy k yx k yy • u v - d xx d xy d yx d yy • u v

Imbalance

The imbalance is a concentrated mass. Then writing the kinetic energy, it comes:

2E B c = m b u2 + v2 + ẇ2 + r 2 b φ2 + 2r b φ ( v cos ϕ -u sin ϕ) (5)

Diaphragm

Different types of casing have been considered.

FE modelling

Shaft line Traction and torsion are discretized using the usual shape functions. Writing ξ = z l and q = w ,β yields:

q(z) = 1 -ξ ξ q 1 q 2
The modified Hermite shape functions for bending used in [START_REF] Chen | Dynamic stability of a Cantilever shaft-disk system[END_REF] are chosen. The degrees of freedom u(z), θ v (z), v(z) and θ u (z) are written as follow:

                         u(z) = N 11 (z) N 12 (z) N 13 (z) N 14 (z)        u 1 θ v1 u 2 θ v2        θ v (z) = N 21 (z) N 22 (z) N 23 (z) N 24 (z)        u 1 θ v1 u 2 θ v2        and                          v(z) = N 11 (z) -N 12 (z) N 13 (z) -N 14 (z)        v 1 θ u1 v 2 θ u2        θ u (z) = N 21 (z) -N 22 (z) N 23 (z) -N 24 (z)        v 1 θ u1 v 2 θ u2        Diaphragm C3
The inner ring of casing C3 is discretized with cubic polynomials in u and v. Defining ζ = s l c the shape functions are:

N 1 (s) = 1 -3ζ 2 + 2ζ 3 N 2 (s) = l c ζ 1 -2ζ + ζ 2 N 3 (s) = ζ λ + 3ζ -2ζ 2 N 4 (s) = l c ζ -ζ + ζ 2
where l c denotes the length of a finite element. The discretized displacement field becomes for an element whose nodes are denoted 1 and 2:

u c (s) = N 1 (s)u 1 + N 2 (s)u 1,s + N 3 (s)u 2 + N 4 (s)u 2,s v c (s) = N 1 (s)v 1 + N 2 (s)v 1,s + N 3 (s)v 2 + N 4 (s)v 2,s
Equations of motion are then derived using Hamilton's principle and the stiffness, gyroscopic and imbalance matrices depend on the angular speed and acceleration.

D R A F T

Contact forces

The forces of particular interest in this study are the contact forces acting between the shaft and the diaphragm where friction is considered. As a first approach, contact is assumed to be punctual at the hot spot, i.e. highest value of rotor eccentricity (point C, see figure 3). With the assumption of small rotations and displacements, the radial gap is equal to the sum of the eccentricity of the cross section's center (distance OR) and the radius of the shaft. The radial contact force component F N is obtained by using Lagrange multipliers method, which guarantees the impenetrability condition. The tangential contact force F T is deduced using the Coulomb friction law F T = µF N . Finally the friction torque takes the form C f ric = µRF N .

Boundary conditions and external loads

As the angular speed is a variable, boundary conditions are necessary for the angular speed law:

-Driven torque: assumed to be constant and equal to C max ; -Resisting torque (alternator): equal to -C max φ Ω nom (model use to obtain an exponential shape).

During the accidental shut down, forces acting on the turbogenerator correspond to aerodynamical and fluids forces:

-Newtonian fluid friction torque:

C newt = -A newt φ; -Aerodynamical friction torque: C aero = -A aero φ2 .
A newt , A aero are positive constants. These friction coefficients have been identified by a EDF's study [START_REF] Fortin | Modélisation de ralentissement de GTA[END_REF]. The accidental shut down is only of interest for our studies: as the shaft line is disconnected from the rest of the network there is no driven torque (C max = 0)

General algorithm

A previous study [START_REF] Legrand | Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion[END_REF] shows that the penalty method is not adapted for the contact detection. A sensitivity study [START_REF] Giraudon-Guilloteau | Modélisation du contact en dynamique: Construction d'un élément simplifié de contact et application à l'interaction rotor/stator[END_REF] on implicit time scheme with contact constraints shows a dependency on Newmark parameters. It turns out that the central finite differences based on Carpenter's work [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF] (prediction-correction algorithm) exactly satisfies the contact detection and ensures the compatibility of the speed and the acceleration. Consequently, equations of motion are discretized in time according to central finite differences. Because of the nonlinear system of equations, the following time procedure is developed as follows: D R A F T Initialization of variables X and ϕ for t 0 and t 1 : X 0 , X 1 and ϕ 0 , ϕ 1 are estimated For n from n 2 to n end do Prediction of (X n+1 , ϕ n+1 ) assuming that the acceleration is locally linear. end For * : when the angular speed become negative, which means the shaft is counter rotating, the calculation is stopped because this case is not of interest.

Search of the solution

* * : the Lagrange multipliers is calculated assuming that the angular position, which is estimated in the neighborhood of the converged solution, is correct. At each time step, the rotor displacements are predicted. The matrices depending on the spinning speed and acceleration are updated. If a contact is detected, the rotating speed and acceleration change according to the friction torque.

SIMPLIFIED MODEL OF SHAFT LINE

In order to check the developed algorithm, a simplified turbine model has been considered (see figure 4(a)). FE models for the turbogenerator and the diaphragm are chosen so that the first eigenfrequencies match the corresponding ones of the real structures. Campbell diagram of the rotor is shown in figure 4(b).

TIME INTEGRATION RESULTS

The developed algorithm was first checked and validated simulating a speed transient. The speed law was given and results have been compared to the ones obtained from a commercial code.

Simulations with casing C1

Case without torsion As a preliminary study case, the mathematical diaphragm C1 is considered. The imbalance mass is set to 5 kg in order to excite the rotor. Figure 5(a) depicts the computed spinning velocity. A zoom of this figure during the contact phase is depicted in Figure 5(b)) and the angular acceleration is plotted in Figure 5(c). According to simulations, convergence of the displacements is achieved. The eccentricity of the node in front of the diaphragm (see figure 6(a)) shows that the contact is well considered: even for an high value of time step that ensures convergence, the residual penetration is negligible (lower than 0.2µm, value to be compared with the gap 8mm). The rotor response spectrum (see figure 6(b)) shows components at the rotating frequency and at the natural frequencies of the model.

The influence of the friction coefficient is shown in figures 7(a) and 7(b). For a small value, just the forward whirl motion is excited. With µ greater than 0.5, the reverse whirl phenomenon is observed. Backward whirl phenomena are pointed out with a lower friction coefficient with a higher imbalance mass too. These two figures show the trivial dependency on the friction coefficient which is directly linked with the contact friction torque. Case with torsion When torsion β is added to the model, results become more complex. Figure 8(a) shows that there is convergence except during contact: during interaction, results locally depend on the time step, which may be due to high frequency phenomena. The spatial discretization of the shaft might be too poor and could explain it. Even with an high value of time step that ensures convergence, the residual penetration is negligible (lower than 0.2µm, value to be compared with the gap 8mm). Further analysis must be performed to understand the results. Figure 8(c) depicts the angle of torsion β. The frequency spectrum of the latter is complex to analyze due to the nonlinear coupling between bending and torsion and the contact nonlinearity. As a intermediate conclusion, the developed procedure has been validated simulating rotor-stator interactions.It has been shown also that additional analysis concerning the d.o.f. of torsion have to be done. However flexibility and damping of the diaphragm are not included in such a casing. The unbalance mass is now set to 45kg and only simulations without considering the torsion d.o.f. are considered.

D R A F T

D R A F T

Simulations with casing C2

Simulations are performed with the diaphragm C2, which is more realistic but simplified. Analogously figure 9(a) depicts the calculated rotating speed and figure 9(b) the corresponding angular deceleration. The results can be compared to computations with casing C1 (see -on figures 9(a) and 9(c)): without flexibility (casing C1), the contact is harder. The contact detection is still well considered as figure 9(c) shows. Figure 10 diaphragm displacements during interaction. Some frequencies of the stator are common with the frequency content of the rotor: it seems that these frequencies are harmonics of the rotating speed and related to the number of impacts per revolution. Loads on bearings directly depend on the casing stiffness. If coupled with a nonlinear solver Reynolds equations solver, this algorithm could compute the real effects of nonlinear bearings and verify that the turbine passes critical speeds without catastrophic consequences.

Authors start to study a 3-D or pseudo 2-D casing FE model in order to add axial deformations of the diaphragm.

CONCLUSION AND PROSPECTS

A numerical tool has been developed to analyze speed transients of a turbogenerator with rotor-stator interactions. The developed contact algorithm have been checked according to the first time integration results. With the use of Lagrange multipliers, the computations don't depend on time step and the contact detection is well considered so that there is just numerical residual penetration.

Figure 1 :

 1 Figure 1: Nomenclature

Figure 2 :

 2 Figure 2: Schematic of diaphragm models

Figure 3 :

 3 Figure 3: Schematic of contact detection

  (X n+1 , ϕ n+1 ) satisfying the governing equations (nonlinear solver) While ( Residual vector of the discretized equations of motion 2 ≤ ε) do If (No penetration) then λ N = 0 else Calculation of the Lagrange multipliers λ N = λ N (predicted solution) * * Correction by radial and tangential (according to Coulomb law) contact forces and then the friction torque C f ric = µRλ N end If New solution guess (X n+1 , ϕ n+1 ) done If (ϕ n+1 ¡ ϕ n ) then Computation stopped * end If Results are recorded Next time step increment a n-1 ← a n and a n ← a n+1

Figure 4 :

 4 Figure 4: Model of shaft line

  time (s) angular speed (rad/s) dt = 10 -5 s dt = 5 • 10 -6 s dt = 10 -6 s dt = 5 • 10 -7 s angular speed (rad/s) dt = 10 -5 s dt = 5 • 10 -6 s dt = 10 -6 s dt = 5 • 10 -7 s

Figure 5 :

 5 Figure 5: Rotating speed and angular deceleration with casing C1

  Frequency content of the rotor without contact

Figure 6 :Figure 7 :

 67 Figure 6: Eccentricity of the node 4 and frequency content of the shaft line

  time (s) angular speed (rad/s) dt = 10 -5 s dt = 5 • 10 -6 s dt = 10 -6 s d.o.f. of torsion β (node 4)

Figure 8 :

 8 Figure 8: Rotating speed, deceleration and the torsion d.o.f. with casing C1

of node 4 Figure 9 :

 49 Figure 9: Rotating speed, deceleration and eccentricity with casing C2

Figure 10 : 4 . 3

 1043 Figure 10: Frequency spectrum when contact occurs

Dof node 4 Figure 11 :

 411 Figure 11: Rotating speed and eccentricity of node 4 with casing C3

  , E d) respectively kinetic and strain energies (k ij , d ij )

	w	traction (Z axis)	β		torsion (Z direction)
	(u, θ v )	bending in (X,Z) plane	(v, θ u )	bending in (Y,Z) plane
	ϕ	angular position	s		curvilinear abscissa
	(E c stiffness and damping coefficients
	φ	angular spinning velocity	Ω nom		nominal rotating speed
	Disk nomenclature
	M D	disk mass	I DX		moment of inertia
	I DZ	polar moment of inertia
	Beam nomenclature
	E	Young's modulus	ν		Poisson's ratio
	G	shear modulus	S		cross section area
	k x I x	transverse shear form factor moment of inertia	λ = I p	12EI x k x GSl 2 transverse shear effect polar moment of inertia
	ρ	mass density	l		element length
	z	axial distance
	Subscripts
	D	disk	b		imbalance
	c	casing

Table 1 :

 1 Radius of the inner ring = outer radius of the shaft + gap (gap=8mm) Numerical characteristics of the shaft line

	Disk
	Inner radius	0.5m	Thickness	0.1m
	Outer radius	1m	Mass density 7860kg/m 3
	Shaft
	E	200GPa	ν	0.3
	Total length Radius Imbalance Imbalance mass 5kg for C1 and 45kg for C2 and C3 Phase shift 10 m Mass density 7860kg/m 3 0.5m k x 6/7 0rad Distance between imbalance and center of the disk 1m Bearing β 0.0002 (0.02% of the stiffness coefficients ) k xx 2 • 10 5 N/m k xy 0N/m k yx 0N/m k yy 5 • 10 5 N/m c xx β × 2 • 10 5 N.s/m c xy 0N.s/m c yx 0N.s/m c yy β × 5 • 10 5 N.s/m Casing C3 E 200GPa ν 0.3 D R A F T Mass density 7860kg/m 3

Paper ABC-123


S W Shaw

An Invariant Manifold Approach to Nonlinear Normal Modes of Oscillation

Keywords: nonlinear vibrations, nonlinear oscillations, invariant manifolds, normal modes, Galerkin

A method for determining the amplitude-dependent mode shapes and the corresponding modal dynamics of weakly nonlinear vibratory systems is de• scribed. The method is a combination of a Galerkin projection and invariant manifold techniques and is applied to a class of distributed parameter vibratory systems. In this paper the general theory for a class of conservative systems is outlined and applied to determine the nonlinear mode shapes and modal dynamics of a linear Euler-Bernoulli team attached to a nonlinear elastic foundation.

Introduction

Traditional methods for determining the effects of nonlinearities on the mode shapes of vibrating systems typically involve an assumed form of harmonic time behavior and series expansions for the frequency of oscillation and the mode shapes in terms of the amplitude of vibration (see, for example, Bennouna and White [1], Benamar et al. [2], or Szemplinska [START_REF] Szemplinska | The Behaviour of Non-Linear Vibrating Systems[END_REF]). A more formal approach is found in the classical work of Rosenberg [START_REF] Rosenberg | [END_REF], which was followed by many subsequent efforts, for example, those of Rand [5][6][7] and Vakakis [START_REF] Vakakis | Analysis and identification of linear and non-linear normal modes in vibrating systems[END_REF][START_REF] Vakakis | [END_REF][10] and their co-workers, all of which used definitions of normal modes in terms of periodic motions for finite-dimensional, conservative, nonlinear systems, typically with two degrees-of-freedom. Recent work by Jezequel and Lamarque [11] uses normal form theory to generate nonlinear normal modes, again for systems with two degrees of freedom. Their approach does allow for the inclusion of special types of damping. Vakakis has used nonlinear normal mode methods to analyze nonlinear effects on mode localization in multi degree of freedom, cyclic systems [12].

The method described herein takes a point of view, which is fundamentally different from that found in previous engineering literature on the subject. Rather than view nonlinear normal modes as synchronous motions of conservative systems, they are viewed as motions on invariant manifolds which are tangent to, and of the same dimension as, the linear eigenspaces in the system phase space. These manifolds are the standard foliation of the center-stable manifold of the equilibrium point of interest, and their existence depends crucially on certain nonresonance conditions being satisfied.

These manifolds are familiar to someone trained in dynamical systems theory, but they have only recently been exploited for the construction of nonlinear normal mode shapes and modal dynamics. Shaw and Pierre developed methods of constructing nonlinear normal modes for both lumped-parameter [13] and distributed parameter [14] vibratory systems in this manner. Such an approach allows for the definition and construction of nonlinear normal modes even in the presence of gyroscopic and/ or dissipative terms. These efforts were followed by the work of Nayfeh and co-workers [15,16], who use both perturbation and invariant manifold techniques to generate nonlinear normal modes, including the investigation of a case with internal resonance [17]. King and Vakakis [18] have used a combination of Rosenberg's ideas with the method presented in [14] in order to generate normal modes for continuous systems. This work also includes an investigation of the stability properties of nonlinear modes.

Along similar lines, much work has been done on the existence and smoothness of these manifolds for finite-dimensional, conservative, Hamiltonian systems. The main result in the nonresonant case is the Lyapunov center theorem which states that if the Hamiltonian is C 1 and the linearized system has purely imaginary eigenvalues +fij (j = l, ... ,m) such that for j =I= k, nk;nj =I= n (n = 1,2,3, ... ), then there exist m two-dimensional, local C 1 invariant manifolds. (This result is due to Lyapunov, and a simple proof is given by Kelly [19).) These are the nonlinear normal mode manifolds of interest here, and each manifold contains a oneparameter family of periodic solutions. Other work in this area includes that of Weinstein and Moser on bounds for the number of periodic orbits which can exist near an equilibrium point [START_REF] Weinstein | Lagrangian submanifolds and hamiltonian systems[END_REF][START_REF] Moser | Periodic orbits near an equilibrium and a theorem[END_REF], and extensive work on the resonant cases; see [22][23][START_REF] Schmidt | A unifying theory in determining periodic families for[END_REF][25], for example. It is important to point out that in the strong resonance cases (n = 1, 2, 3), the number of families of periodic solutions is not generally the same as the number of pairs of complex conjugate eigenvalues. The nonlinear normal modes defined in the present work do not generally exist in such cases.

The method developed by Shaw and Pierre [14] is a systematic means for generating nonlinear normal mode manifolds for nonresonant partial differential equations which arise in structural vibrations. It involves some rather cumbersome calculation steps and the results depend on the selection of a suitable base point. The method described herein is aimed at providing an alternative approach which avoids some of these difficulties. The first step is to apply a usual Galerkin projection of the partial differential equations of motion onto a set of basis functions which consist of the mode shapes of the linearized system. This leads to a set of ordinary differential equations which describe the dynamics of the nonlinear system in terms of the linear modal amplitudes. This standard first step leads to equations of motion which are uncoupled at linear order. The invariant manifold approach developed in [13] is then applied in order to determine how the nonlinear coupling distorts the mode shapes and dynamics in the weakly nonlinear system. At this step it is important to treat the linear modal amplitudes and velocities as independent variables, as they provide the coordinates which parameterize the manifolds. Overall, this procedure provides a straightforward means for determining the amplitudedependent shapes of natural vibration for weakly nonlinear, distributed parameter systems. The resonance cases are flagged directly by singularities in the equations which are to be solved for the coefficients describing the manifolds.

One important conclusion of the results obtained is that one must be careful when generating reduced order models for the vibration of continua using modal projections. The typical step of projecting the partial differential equation onto a subset of the linear modes must be done with care, especially if one requires more than first-order nonlinear information about the system at hand. This shortcoming was observed by Nayfeh et al. [26] in the analysis of the dynamic response of a valve system, and is systematically explored here.

The paper is arranged as follows. In Section 2 the general method is described for a class of weakly nonlinear partial differential equations which arise in mechanical vibrations of distributed parameter systems. The method is applied to an example problem in Section 3 and the paper closes with some remarks in Section 4. It should be noted that the methods presented herein were motivated by the analysis of nonlinear mechanical structures, but they may have• a wider range of applicability.

The Method

Equations of Motion

Consider one-dimensional vibrations of a one-dimensional continuum which occu~ pies the region described by the spatial variable s E (0, 1). The equations of motion are assumed to be of the form a 2 w(s, t)

--at--::-2 - + L(w(s~ t)) + N(w(s, t)) = 0 (1)
with linear boundary conditions at s = 0, 1 of

B 1 ( W ( S, t)) is= 0 = 0, B 2 (w(s, t))ls=l = 0. ( 2
)
The equilibrium configuration about which the system vibrates is taken to be w(s, t) = 0, that is, L(O) = N(O) = 0. The operator L is taken to be linear and self-adjoint, while the nonlinear operator N is assumed to possess smoothness properties that allow it to be expanded to any order desired. The linear system, N = 0, is assumed to have nonzero natural frequencies flj with nj;nk::;:: 1, 2, 3, ... and associated normal mode shapes cfJ/s) for j = 1, 2, 3, ... , which are normalized according to the definition given in equation ( 4) below. The system eigenvalues are thus ±i!lr The inner product between real-valued functions on the unit interval is defined in the usual manner as

(f(s),g(s)) = l 1 f(~)g(~)d~. 0 (3)
The following general properties then hold for the linear system:

(¢i(s),¢/s)) = sij' (4) ( ¢i(s), L(¢/s))) = 11J5ij'
where oij is the Kronecker delta.

The solution of the equation of motion is now assumed to be of the form of an infinite series involving the linear normal modes and their amplitudes as functions of time as w(s,t) = L 4>/s)q/t). j ... l (5) Here qi(t) represents the contribution of the jth linear mode to the response. This solution form is substituted into the equation of motion, which is then projected onto the linear modes by using the inner product. Using the properties of the linear system, this results in the following set of differential equations of motion in terms of the linear modal amplitudes: i = 1' 2, 3, ... ' (6) where an overdot denotes a time derivative, q is the infinite-dimensional vector comprised of the q;'s, and the nonlinear coupling terms in Gi are determined from i = 1,2,3, ... . (7) These can, to cubic order in the q/s, be expressed as 00 00 00 00 co

1-1 m = 1 /=1 m = 1 n-1 i = 1,2,3, ... , (8)
where the vilm and J.Lumn coefficients depend on the nonlinearity and the mode number i.

In order to determine the nonlinear normal modes in terms of invariant manifolds by the method developed in [13], these equations of motion are written in first-Jrder form as qi(t) = Pi(t), [START_REF] Vakakis | [END_REF] pi(t) = -OfqJt) -Gt(q(t)), i = 1, 2, 3, . . .. We restrict our attention to systems for which OJOi is not equal to an integer for i :f:.j.

In the following, the present approach will be compared with the more commonly used method of "direct projection" onto a single mode. That approach employs a single mode, say c/Jj, in the expansion (5) and projects the equations of motion onto that mode. The i = j equation of motion is the only one of interest in that case, and it has the same form as equation ( 9) with i = j, except that the multiple sums in Gj are reduced to the terms vjjjq](t) + J.Ljjjq}(t) + •• • . The effects of the missing terms are considered in detail below.

The Normal Mode Invariant Manifold Equations

The key observation made in [13} is that for either a linear or a nonlinear system a normal mode motion is one of a family of motions for which the system behaves like a single degree of freedom oscillator; that is, like a second-order nonlinear oscillator. Such a motion will take place on a two-dimensional invariant manifold. For linear oscillatory systems these manifolds are the planar eigenspaces associated with the complex conjugate pairs of eigenvalues. The existence of these eigenplanes for the linear system is guaranteed for this class of systems with the given assumptions on the frequencies. For weakly nonlinear systems of the type being considered here, the curved manifolds are necessarily tangent to the linear eigenspaces at the equilibrium (w(s, t) = 0). Once these manifolds have been determined, the manner in which the linear modes combine to form the nonlinear normal modes will be evident. In addition, the dynamics on these manifolds, which represent the associated dynamics of the nonlinear normal mode, are easily recovered.

The construction of the modal manifolds is carried out in the usual manner. First, one must choose coordinates which describe the manifold; in this case (qk, pk) are a natural choice for the k th normal mode. Next, it is assumed that if the motion of (qk, pk) is known, then the motion of all other (qi, p) pairs can be described in terms of (qk, pk); this reduces the dynamics of the entire system to that of a single degree of freedom system or, equivalently, the dynamics are restricted to a twodimensional manifold. Also, the equations of motion must be satisfied for such a specific type of motion. These ingredients lead to a set of equations for the normal mode manifolds which, although not solvable in general, easily admit solutions in series form in the variables (qk, pk), just as is done in the construction of other types of invariant manifolds (see, for example, [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] or [START_REF] Gilsin | [END_REF]). The solution of the manifolds contains the desired information regarding the composition of the nonlinear normal modes in terms of the linear normal modes.

By defining (qk, pk) = (uk, vk), a normal mode motion is assumed to exist and be (at least locally) expressible in the form

qi(t) = Qik(uk(t), uk(t)), Pi(t) = Pik(uk(t) , uk(t)),
where, by definition,

i = 1,2,3, ... , Qkk(uk(t), vk(t)) = uk(t), Pkk(uk(t),vk(t)) = vk(t).
(10) (11) In this way, the dynamics of the entire structure is specified by the dynamics of the kth linear modal amplitude and velocity.

The equations for the Qik's and Pik's are obtained by first taking a time derivative of each item in equation (10) to obtain . aQik . aQ;k .

qi = --uk + --uk , auk avk (12) .
aPtk . aPik .

Pi= --uk + -v k , auk avk i = 1,2,3, ... .
The following substitutions from the equations of motion and equation (10) are then made:

where Qk is the vector of Qik's with i = 1, 2, ... . This process eliminates all time derivatives from equation (12), resulting in the following equations for the Qik's and P;k's in terms of uk and vk:

i = 1,2,3, ... , (13)
i = 1, 2, 3, .. . . (14) Note that the i = k equations are trivially satisfied.

A closed form solution of these equations is not generally attainable (except in the case of so-called similar normal modes [START_REF] Vakakis | Analysis and identification of linear and non-linear normal modes in vibrating systems[END_REF]), but a local solution near the origin can be obtained in terms of a series expansion. Matters are simplified by noting that for conservative, nongyroscopic systems, only those terms which are consistent with constant amplitude, standing wave normal modes are required. Thus, to cubic nonlinear order, the series is given by

Qik(ub vk) = alikuk + a3iku~ + asikvl + a6iku~ + astkukvl + •••, i = 1,2,3, ... , (15)
i=1,2, 3, ... , where the following identities for i = k are obvious from equation (11):

alkk = bzkk = 1, afkk = blkk = 0, j * 1, l * 2.
(16) (17) (The missmg afik and b 1 ;k coefficients are identically zero in the conservative, nongyroscopic case; see [13] for examples where they are nonzero. This is easily confirmed by including them in the above and subsequent steps.)

Substituting the expansions (15) into equations ( 13) and ( 14), expanding in terms of uk and vk, and gathering the coefficients of the u'f: 1 uf:' 1 terms, where m 1 , m 2 = 0, 1, 2, 3 and m 1 + m 2 < 4 (we are working only to cubic order here), yields a set of equations in terms of the unknowns afik and buk• These equations are as follows.

At linear order: from equation ( 13), the uk coefficient (18) from equation ( 14), the uk coefficient (19) At quadratic order: from equation ( 13), the ukuk coefficient [START_REF] Weinstein | Lagrangian submanifolds and hamiltonian systems[END_REF] from equation ( 14), the u~ coefficient [START_REF] Moser | Periodic orbits near an equilibrium and a theorem[END_REF] from equation ( 14), the v'f coefficient (22) At cubic order: from equation ( 13), the u~uk coefficient (23) from equation ( 13), the vt coefficient from equation ( 14), the u~ coefficient

-a6ik0f -L L Vum(a3lm + a3mt) -ILtkkk l= 1m= 1 00 00 (24)
= -b7ikO%-bz;k L L Vktm(a3lm + a3ml)-1Lkkkkb2ik-b4ikvkkk; (25) l=lm=l from equation ( 14), the ukv~ coefficient

-aSikD.f -L L 11 um<astm + asmt) l=lm=l 00 00 = 2b7ik -3b9ik0~ -b2ik L L vklm(asim + asmt). ( 26
) l= 1m= 1
Several features of these equations are worth noting. First, they are linear in the unknowns, and nonhomogeneous terms arise from the nonlinearity. Also, the known solution for the i = k case [see equation (11)] is consistent with these equations. The solution for the coefficients of the linear, quadratic, and cubic terms are now considered in turn.

Regarding the linear terms, it is immediately clear from equations ( 18) and ( 21) that [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] since nk :fo ni and using identity (16) above. This is not unexpected, since it simply says that the only contribution at the linear order to the kth nonlinear normal mode comes only from the kth linear mode. As described in detail below, this implies that the linear modes are sufficient for correctly capturing the leading order nonlinear dynamics, as nonlinear corrections to the mode shapes will riot affect the leading order nonlinear terms in the equations of motion. However, the modal distortions must be accounted for if one considers higher order dynamics.

Regarding the quadratic terms, these represent three linear equations in three unknowns, and are easily solved when the equations are nonsingular. The solution is given by vikkca

; -2nD tlzik -2vikk tlzik ' 2vikknf b4ik = --- 112ik (28) Vi =P k.
The coefficients for i = k are given in equation ( 16) above. Note that these equations are singular in two cases: (1) If any ni = 2flk, an internal resonance condition arises in which the i and k modes cannot be uncoupled. ( 2) If any ni = 0, a condition not allowed by the original restrictions on the natural frequencies.

The solution for the cubic terms involves four linear equations in four unknowns, with nonhomogeneous terms which arise from the cubic order nonlinearity in the equation of motion and from the solution at quadratic order. The appearance of terms from lower order at the current order is typical of such procedures and is well known in, for example, perturbation and normal form calculations. The solution for the cubic coefficients, for i -:f:. k, are given by

(7!1~-flf)( 1-Likkk-vkkkb4kk + o-i3) + 2vkkka5kkD~(3!1~-Df)-2Dko-is tl3ik 00 00 o-ij = L L. 11 um(ajim + ajml)' l= 1 m = 1 ( 29
)
j = 3,5.
The coefficients for i = k are given in equation ( 17) above. Note that these equations are singular whenever fl; = 3flk or nj = nk, resonance conditions not allowed by the original restrictions on the natural frequencies.

The uij terms can be simplified by noting that, although the a 11 m terms are not symmetric in the l and m subscripts, they are summed over both l and m, so that

00 co o-ij = L L Vilm(ajlm + ajml) l= 1 m= 1 00 00 = 2 L L Vumajlm' {= 1 m= 1
This is useful for calculations in specific examples .

.

Nonlinear Normal Mode Shapes and Dynamics

j = 3, 5. ( 30
)
The nonlinear mode shapes are now reconstituted by using equation ( 5) with the q/s replaced by their modal versions, the Q/s.

where the reader is reminded that u k and u k represent the displacement and velocity of the kth linear modal amplitude during the kth nonlinear normal mode motion. These normal modes represent standing wave motions of the system as can be seen by the following consideration. At any instant, say t = t *, at which the linear modal velocity is zero, that is, vk(t*) = 0, the entire velocity field is zero: (awk(s, t*))/ at = 0 'Vs E [0, 1]. Thus every point of the system reaches its maximum displacement at the same instant. This fact allows mode shapes, denoted here by Wk(s ), to be plotted in a manner just as is done for linear systems. They are simply the zero-velocity configurations of the system, which are given by 00

Wk(s) = wk(s, t*) = c/>k(s)u't + E c/>;(s)( a 31 ku'k 2 + a 6 ;ku'k 3 + ••• ), (33) i=l, i=Pk
where u\ = uk(t*) is the maximum value that uk(t) achieves during the motion.

Note that at any instant, say t = t 0 , at which the linear modal amplitude is zero, uk(t 0 ) = 0, the entire displacement field is generally not zero, wk(s, t 0 ) =I= 0. This is due to the presence of the nonsymmetric, quadratic nonlinearities, and in particular the term a 5 ;kvf(t 0 ). For systems which do not possess quadratic nonlinearities, wk(s, t 0 ) = 0 to leading order, since all a 5 ik terms are zero. In that case, every point of the system passes through zero at the same instant.

Note that, in contrast to the linear mode shapes, the nonlinear mode shapes depend on the peak amplitude of motion. The distortion from nonlinearities is contained in the quadratic and cubic terms, and in particular is expressed here in terms of the linear modes shapes with the contribution of the ith linear mode shape on the kth nonlinear mode shape being captured in the coefficients a 3 ;k and a 6 ik• In addition, during a given modal motion, the shape of the structure varies in time, as can be seen by considering equations ( 30) and ( 31) and assuming that uk(t) and vk (t) are periodic in time. This aspect of the modal response is not allowed by the methods used in [1][2][START_REF] Szemplinska | The Behaviour of Non-Linear Vibrating Systems[END_REF], which essentially assume a fixed shape whose amplitude varies according to a second-order nonlinear equation of motion.

The dynamics of a nonlinear normal mode are obtained by simply restricting the equation of motion to the modal manifold. This is accomplished by substituting the modal expressions for q and p given in equation ( 15) into the equations of motion (6). The result is

ilk + n~uk + akuz + f3kuf + 'Ykukit~ + ... = 0, ak = v kkk, f3k = P..kkkk + a-k3' /' k = a-ks' (34)
which, due to the aki terms, differs from what one would obtain by simply projecting the equation of motion onto the kth linear mode, <f>k(s). However, a projection onto c/>k(s) gives exactly the same dynamics as this result if one is •working only to first nonlinear order; this is discussed in more detail below. The amplitude-dependent frequency of nonlinear oscillation can be determined form equation (34) by using a perturbation method (here Lindstedt's method is used). It is found to be

9{3k0~ + 3yk0k -lOaf 2 n + u* + ... k 24.0.~ k ' ( 35
)
where uZ is the peak amplitude of uk during the motion.

Remarks

At this point it is interesting to examine some special cases. First, consider the general system with quadratic and cubic nonlinearities, but for which only the first-order nonlinear modal dynamics are of interest. These are trivially obtained by simply using a projection onto the single linear mode c/>k(s), which directly yields the quadratic coefficient ak = vkkk in equation (34). In this case there is no need to compute the nonlinear mode shapes, even to first nonlinear order. However, note that the frequency information will not be correct to leading nonlinear order; see equation ( 35). If one is interested in second-order nonlinear modal dynamics, which are required to obtain the correct frequency estimate in this case, it is clear from equation (34) that the first-order nonlinear corrections for the mode shapes must be known. These are expressed in terms of a 3 ;k and a 5 ;k , which combine with the vilm to give ak 3 and aks• For a system which possesses symmetric nonlinearities, in which case Pum = a 3 ;k = a 5 ik = b 4 ik = ak 3 = aks = 0, a projection onto the single mode </>k(s) will give the correct nonlinear modal dynamics to cubic order, via the required coefficient {3k = P..kkkk• Again, the correct higher (e.g., fifth-) order dynamics will require the appropriate nonlinear mode shape terms. Continuing this line of reasoning, a generalization of the development above provides the cubic order nonlinear terms for the mode shapes, and this information can be used to generate the correct nonlinear modal dynamics which includes terms up through and including quintic order in uk and ilk. Boivin et al. have carried out such calculations and have demonstrated the benefits of including such higher order terms in simulation studies [START_REF] Boivin | Normal normal modes, invariance, and modal dynamics approximations of non-linear systems[END_REF]. The fact that shape corrections are required to generate correct higher order dynamics has been used in estimating wave speed-wave amplitude relationships for water waves (see, e.g., Lamb [START_REF] Lamb | Hydrodynamics[END_REF], article 250), but has been overlooked in the structural dynamics community, where reduced order models for systems with quadratic and cubic nonlinearities are typically obtained by direct projections onto the linear modes of interest.

Example

The example considered here is a simply supported, linear, Euler-Bernoulli beam attached to a nonlinear elastic foundation. Figure 1 depicts the physical system. The linear problem is well known and can be found in textbooks (e.g., [START_REF] Weaver | Vibration Problems in Engineering[END_REF]). This problem was selected since the linear mode shapes are simple functions, making closed form solutions possible. The general approach will be made clear from this relatively simple example, and it will be evident where one will be required to carry out numerical computations in less straightforward cases.

The nondimensional equation of motion, expressed out to cubic order, is a 2 w(s,t) a 4 w(s , t)

2 3 at2 + as4 + KW(S , t) + p(w(s, t)) + 7](W(S, t)) + ••• = 0 (36)
with boundary conditions

w(O, t) = w(l, t) = 0, a 2 w a 2 w -2 (0, t) = . -2 (1, t) = 0.
as as

(37)
Here K, p, and TJ represent the linear, quadratic, and cubic coefficients of the foundation, respectively. The linear problem ( p = TJ = 0) has natural frequencies and normal mode shapes, respectively, of

n k = v K + (k1T ) 4 , cpk(s) = fi sin(k1Ts ).
w(s,t) (

) 38
Two cases will be considered. In Case I the system with both quadratic and cubic nonlinearities will be studied. First-order nonlinear terms for the mode shapes will be generated, and from these the second-order dynamics will be obtained. For Case II, the system with symmetric nonlinearities, p = 0, will be studied, and for this system the first -order nonlinear terms will be computed for both the mode shapes and the modal dynamics.

Case I

The coefficients vilm and ILilmn need to be generated; this is carried out as follows for this system. First, the nonlinear operator acting on the assumed form of the solution, equation ( 5), is written as

N(w(s,t)) = p(w(s,t)) 2 + 1J(w(s,t)) 3 + ••• N( .£: 4>/s)q/t)) = p( .£: ¢/s)q/t)) 2 + 11( _[ 4>/s)q/t)) 3 + ••• j=l J=l J=l (39) = p L L 4> 1 (s)<f>m(s)q 1 (t)qm(t) l= 1 m = 1 00 CQ 00 + 1] L L L 4>£(s )<J>m(s )<J>/s )q 1 (t )qm(t )qn(t) + ••• . l=-1 m=l n=l
Next, this is projected onto the ith linear mode according to equation (7) in order to obtain Gi: 00 00 G;(q(t)) = L 1 lf>Jt)p L L <J>,(t)lf>m(t)q,(t)qm(t)dt+ L 1 lf>;(~)1J

0 l=lm=l 0 • L L L <f>t(t)l/Jm(g)lf>n(g)ql(t)qm(t)qn(t)dg+ •••. ( 40
)
1=1 m=l n=l

The order of summations and integration is switched in order to obtain the coefficients vilm and J.Lamn from equation [START_REF] Vakakis | Analysis and identification of linear and non-linear normal modes in vibrating systems[END_REF]. The result is (41) given in equation ( 38). The nontrivial a 3 ik and a 5 ik simplify to

4/2 pk 2 (1 -( -l)i)(-K + i 4 7T 4 -2k 4 7T 4 ) a3ik = l7T(4k2 -j2 )( K + i47T4 )(3K -l47T4 + 4k47T4) ' (44) i * k.
Here the internal resonance condition is obvious: K = l1r 4 (i 4 -4k 4 )/31 is exactly the condition for which ai = 2Gk. At this point it is also seen that the contributions of the ith linear mode shape to the kth nonlinear normal mode decreases at the rate of i-7 for large i. It is also worthwhile to note that a 3 ik and a 5 ik are zero for all even values of i, including i = 2k, and for i = k (even or odd). Thus, only the odd, that is, the symmetric, linear modes contribute to the nonlinear modes at this order.

The summation terms uk 3 and uk 5 are obtained by direct substitution of the above expressions for vilm' a 3 ik, and a 5 ik into the expression for uij given in equation ( 28), and making the appropriate subscript assignments. These equations are rather unenlightening and are not presented explicitly here. It is sufficient to note that they converge since the coefficients of uk 3 decay in l and m like 1-3 m -a for large l, m while for uk 5 the decay is at the very rapid rate of z-12 m -7 .

Convergence is also aided by the fact that many of the terms in the sums are zero. For plotting purposes, it is convenient to define the following function

j = 2,3, (45)
which, for each mode number k, are functions of only the linear foundation coefficient K. Figure 2 depicts these functions for the first three modes. In each case the resulting plots indicate a close approximation to the exact value of the doubly infinite sum. The number of terms was varied and convergence was inspected visually; the typical required number of nonzero terms was in the range 3-6. Note that the sums are less significant for beams with a large linear stiffness, and it is easily shown that for large K, Sk 3 -K-1 and Sks -K-2 • These decay rates are interrupted by values of K at which internal resonances arise, at which points the S kJ are singular. In the K range shown in Figure 2, only the 2:1 resonance between the first and the third linear modes appears [at K = 7T 4 (4 X 3 4 -1) /3 :::: 10, 488].

All terms required for the mode shapes and modal dynamics are now in hand.

In this case the peak-amplitude nonlinear normal mode shapes are easily expressed to first nonlinear order by using equation ( 33), where the a 3 ik are given in equation ( 43) above, the cubic order terms are ignored, and six nonzero terms were used in each summation. Plots of the peak-amplitude mode shapes for the cases k = 1, 2 are depicted in Figure 3 for various amp1itudes, specified here by u't. In each case the zero-velocity beam configurations for both the "up" and the "down" positions are shown, indicating the lack of symmetry which arises since the foundation is softer (relative to the linear foundation) in the "down" direction and stiffer in the "up" direction. (See the Appendix for the amplitudes used for the "down" configuration.) It is worth noting that the effect of the elastic foundation does not have a decreasing effect on higher modes. This is in contrast with the linear case [cf. In order to determine the modal dynamics, the coefficients vkkk' J.Lkkkk' uk 3 , and ak 5 are needed. They are given above for the general case and are used to produce the oscillator (46) which governs the dynamics of the kth nonlinear normal mode. The amplitudedependent frequency of nonlinear oscillation can be written as (47) where ut is the peak amplitude of uk during the motion and the coefficient Ak is obtained directly from equation ( 35) and is given by A.CK, p, 1Jl = (9(K + e"•l( ~ 1J + s.,CK)p 2 )

(24( K + k47T4)3/ 2),
where the terms S k/ K) are given in Figure 2.

(48)

For a fixed amplitude of motion, the frequency of oscillation is seen to be linear in the cubic foundation coefficient TJ, quadratic in the quadratic foundation coefficient p, and has a complicated dependence on the linear foundation coefficient K. In order to plot Ak and examine the difference between the usual direct projection result and the present result, it is convenient to express A k as

Ak(K, p, i]) = 17/f(K) + p 2 (/t(K) + /ff 2 (K)), 9 jf(K) = 16flk' (49)
where /{ and ft are the terms one would obtain by a direct projection, and fff 2 involves the S kj terms and arises from the nonlinear distortions to the mode shapes. The net effect of the nonlinear mode shape depends on the relative sizes of iJ and p and the size of fff 2 relative to ff and ft. Plots of these three functions for the first two modes are given in Figure 5.

The importance of including the modal distortions in the modal dynamics can be demonstrated by a single case: that for the second mode, k = 2, with 17 = 0. Since f!P = 0, the only nonzero term in A 2 is proportional to f!J 2 , which is precisely the term that is neglected if one ignores the coupling between the linear modes. The magnitude of A 2 in this case is given by p 2 fi_ 2 , where f!J 2 is shown in Figure 5. It is seen that, while the correct result of zero frequency shift is approached as K ~ co, in the lower range of K the direct projection method can lead to significant discrepancies.

It is thus confirmed that, although the calculations are nontrivial, the first-order nonlinear corrections to the mode shapes provide the correct modal dynamics to second order, and that these are nontrivial in some cases. We now turn our attention to a simpler case, in which the foundation contains only cubic nonlinearities. Note that fi. 1 = 0.

Case II

The calculations here follow the above with p set equal to zero. Since only first-order information regarding mode shapes and dynamics is to be obtained, the calculations are much less involved. The modal dynamics are immediately obtained from the results of Case I with p set equal to zero, and are precisely what one would obtain from a simple projection of the equations of motion onto <f>is ).

For the nonlinear mode shapes, the following quantities are required and are easily obtained using equation ( 41): It is clear that the above procedure can be generalized in a straightforward manner to other situations. In fact, in [13] the complete theory has been completely worked out for nonresonant systems with a finite number of degrees of freedom possessing all possible quadratic and cubic nonlinearities. This includes cases which are nonconservative and/or gyroscopic in nature.

1-Likkk = 0, -11 1-Likkk = 2' Vi =f:. k,3k, (50) i = 3k.
In [14] an alternative formulation is used for defining nonlinear normal modes of continuous systems in which the dynamics of the entire structure is specified by the dynamics of a single point on the structure, s = s 0 • While the method in [14] is more naturally suited to handle some problems, such as those with nonlinearities in boundary conditions, the solution procedure is rather cumbersome and the results are sensitive to the choice of the base point s 0 • The present formulation is equivalent, yet significantly simpler to implement. The present approach is capable of handling discrete nonlinear elements. For example, a nonlinear spring attached at some point along a beam can be handled by simply using delta functions in the spatial variable (see [START_REF] Boivin | Normal normal modes, invariance, and modal dynamics approximations of non-linear systems[END_REF]). If nonlinear effects exist at the boundaries, the method described in [14] may be easier to apply, although a suitable change of variables may render the equation of motion in the form of equation (1).

A possible extension of these ideas is in the area of nonlinear waves. Here, one would seek special solutions of the governing partial differential equation which represent the dynamics of some low dimensional dynamical system. The method described in [14] is well suited for this. For example, for determining the "single mode" dynamics of a nonlinear wave equation, one simply says that if the displacement and velocity of a single point are known, then the entire displacement and velocity field are determined by these.

One issue which has remained unexplored in the current work is a consideration of the stability of nonlinear normal modes. While much work has been done in this area using traditional methods for two degree-of-freedom systems, it has not been considered in the context of the present method. This has important implications with regard to whether or not these nonlinear normal modes will be observable. The usual parametric instabilities between coupled modes are flagged by the present method in the form of singularities. However, other instabilities may occur, especially at large amplitudes (see [START_REF] Vakakis | [END_REF], for example).

The method described herein can be thought of as a means of generating the "best" single-mode model of the given system. A generalization of the procedure has been developed which can be used to generate invariant, reduced-order models of arbitrary order. In order to achieve this, one expresses the response in a form in which the unmodeled modes are expressed as functions of the modeled modes. This idea is of potential use in developing low-order models which may be helpful in computational algorithms or in the development of control schemes. These ideas have some relation to those of the nonlinear Galerkin methods which are used to approximate inertial manifolds [START_REF] Marion | [END_REF], although there is an important difference since that method is applicable to systems with dissipation whereas the systems of interest here are conservative.

A final note: While the above method is in some respects similar to standard perturbation techniques, a comparison between the results obtained by this method and those obtained using a combination of harmonic balance and eigenfunction expansions (see, e.g., [1][2][START_REF] Szemplinska | The Behaviour of Non-Linear Vibrating Systems[END_REF]) shows that different results are obtained for the nonlinear mode shapes, even for the simple example presented here. The approxi-mations for the frequency of nonlinear modal oscillations agree to first order, but the discrepancies in the mode shape approximations at first order will lead to differences in higher order frequency estimates. Since the present method is based on the fundamental principle of dynamic invariance, the results obtained here are the correct ones, and naive perturbation strategies should be applied with caution.

Fig. 1 .

 1 Fig. 1. Beam attached to a nonlinear elastic foundation.

Fig. 2 .

 2 Fig. 2. The functions Ski versus 1<, fork = 1, 2 and j = 3, 5.

Fig. 2 .Fig. 3 .

 23 Fig. 2. (Continued).

Fig. 3 .Fig. 4 .

 34 Fig. 3. (Continued). (d) k = 2, second mode in "down" state; (e) k = 3, third mode in "up" state; (f) k = 3, third mode in "down" state.

Fig. 5 .

 5 Fig. 5. The functions ff, ft, and f% 2 versus K for k = 1, 2.

Fig. 5 -

 5 Fig. 5-Continued.

Fig. 6 .

 6 Fig. 6. Nonlinear normal mode shapes at peak deflection, shown as solid lines, for 1< = 10 4 , p = 0, and 'TJ = 10 6 . The mode shapes are shown for two amplitudes: u'% = 0.02 and 0.04. The linear mode shapes for the same uZ values are shown as dashed lines. (a) k = 1, first mode; (b) k = 2, second mode.

1 Fig. 7 .

 17 Fig.7. Nonlinear normal mode shape at peak deflection for k = 2, K = 13,800, p = 0, and "' = 10 6 , which is close to an internal resonance. Here u~ = 0.04 and the linear mode shape for the same ui value is shown as a dashed curve.

Fig. 8 .

 8 Fig. 8. The invariant manifolds Q 31 and P 31 versus u 1 and v 1 for K = 10 4 , p = 0, and 11 = 10 6 •

Acknowledgments

This work was supported in part by grants from the National Science Foundation and the Army Research Office. The author is grateful to Christophe Pierre, Nicolas Boivin, the reviewers, and, especially, Steve Wiggins for several helpful suggestions regarding this work.

Note that this result is specific for this nonlinearity, but not to the specified boundary conditions. Also, note that to this point the fact that p and TJ are constant has not been used. In fact, the development above is valid for the case in which the nonlinear foundation stiffness depends on s. One simply uses p = p( ~) and TJ = Tt(g) in equation (41). Also, note that the subscripts on v and f.L are completely interchangeable.

For the boundary conditions at hand and constant p and TJ, the coefficients are obtained from simple integrations involving sine functions. The required terms for the first-order mode shapes and second-order dynamics are

37]

I.Lkkkk = 2.

Note that the values of vkkk and ILkkkk are irrelevant for the mode shapes since a 3 kk and a 6 kk must, from the formulation, be zero by definition [see equation (17) and the surrounding discussion]. However, they are required for the modal dynamics.

The corresponding values of a 3 ;k and a 51 k are, for i ¢ k, then

where .6. 2 ik is easily determined from equation ( 28) and the !1/s for this problem are

Vi + 3k,

'

Note that K = 9k 4 7T 4 is the resonance condition for which 11 3 k = 311k.

(51

The solution is now available in a form sufficient for plotting mode shapes of the beam as described in the previous section. The peak-amplitude nonlinear normal mode shape is given by equation (33). In this case the correction to the kth nonlinear mode arises from only the 3k th linear mode and the mode shape is easily expressed as Figure 6 depicts the first two mode shapes for two amplitudes. Note that in this symmetric case the "up" and "down" peak configurations are simply reflected images of one another, so only the "up" version is shown. Also, note that the second mode is virtually unaffected by the given . nonlinearities at the amplitudes shown. However, Figure 7 depicts the same amplitudes but with K at such a value that the second mode is near a 3:1 resonance with the first mode, and here the nonlinear distortion is significant. This is an indication that the modes are becoming highly coupled, exactly as expected near the resonance.

Note that in this case, the foundation has a lesser effect on higher modes, just as in the linear case. This can be directly seen from the coefficients in equation (53), which decrease like k-4 as k increases.

Another instructive way to view the normal modes is directly in terms of their manifolds. This can be done most simply for the present case since only two linear mode shapes combine to form the nonlinear normal mode. Consider the (infinitedimensional) phase space spanned by the state variables (q;, P; ). In terms of these variables, the normal mode manifold for the kth mode is given by the functions (Qik' Pik) [see equation (10)], which are locally represented by the series approximations (15). In order to plot the kth nonlinear normal mode in the phase space, the functions (Q 1 •k , P;k) are plotted versus (qk , pk) for all i. (The case i = k will give a trivial identity and is not of interest). While the complete representation generally involves an infinite number of such plots, only those linear modes which provide significant contributions to the nonlinear mode are of interest. In the present case, the single linear mode ¢ 3 k provides the complete correction to first nonlinear order, and thus plots of (Q 3 k ,b P 3 k, k) completely capture the first-order modal manifolds. An example of this for the first mode, k = 1, is depicted in Figure 8, which clearly indicates the amplitude-dependent nature of the nonlinear normal mode in terms of the bending of the manifolds.

Appendix

Case[

The first-order nonlinear solution for the modal oscillator (34), subject to initial conditions (uk(O), vk(O)) = (X, 0), is obtained using Lindstedt's method and is given as follows:

Note that while this starts at a displacement of X with zero velocity, one-half cycle later the velocity is again zero but the displacement is not equal to -X, but is given by Case II 2ak -X-X 2 --+•••. The first-order nonlinear solution for the modal oscillator (34) with ak = 'Yk = 0, subject to initial conditions (uk(O), vk(O)) = (X, 0), is obtained using Lindstedt's method and is given as follows:

where wk = nk + X 2 (3,Bkj80V + ... . This starts at a displacement of X with zero velocity and one-half cycle later the velocity is again zero and the displacement is -X.