C Rodda
email: rodda@b-tu.de

I D Borcia

P Le Gal

M Vincze

U Harlander

Baroclinic, Kelvin and inertia-gravity waves in the barostrat instability experiment

Keywords: Baroclinic instability, differentially heated rotating annulus, waves in rotating and stratified fluids, inertia-gravity waves, Baroclinic, Kelvin and inertia-gravity waves

The differentially heated rotating annulus is a laboratory experiment historically designed for modelling large-scale features of the mid-latitude atmosphere. In the present study, we investigate a modified version of the classic baroclinic experiment in which a juxtaposition of convective and motionless stratified layers is created by introducing a vertical salt stratification. The thermal convective motions are suppressed in a central region at mid-depth of the rotating tank, therefore double-diffusive convection rolls can develop only in thin layers located at top and bottom, where the salt stratification is weakest. For high enough rotation rates, the baroclinic instability destabilises the flow in the top and the bottom shallow convective layers, generating cyclonic and anticyclonic eddies separated by the stable stratified layer.

Thanks to this alternation of layers resembling the convective and radiative layers of stars, the planetary's atmospheric troposphere and stratosphere or turbulent layers at the sea surface above stratified waters, this new laboratory setup is of interest for both astrophysics and geophysical sciences. More specifically, it allows to study the exchange of momentum and energy between the layers, primarly by the propagation of internal gravity waves (IGW).

PIV velocity maps are used to describe the wavy flow pattern at different heights. Using a co-rotating laser and camera, the wave field is well resolved and different wave types can be found: baroclinic waves, Kelvin, and Poincaré type waves. The signature of small-scale IGW can also be observed attached to the baroclinic jet. The baroclinic waves occur at the thin convectively active layer at the surface and the bottom of the tank, though decoupled they show different manifestation of nonlinear interactions. The inertial Kelvin and Poincaré waves seem to be mechanically forced. The small-scale wave trains attached to the meandering jet point to an imbalance of the large-scale flow. For the first time, the simultaneous occurrence of different wave types is reported in detail for a differentially heated rotating annulus experiment.

Introduction

Convective flows driven by temperature gradient are ubiquitous in geophysical and in astrophysical systems. In a large variety of settings, ranging from salinity stratified water basins [START_REF] Boehrer | Double-diffusive convection in lakes[END_REF] to regions of planetary and stellar interiors [START_REF] Medrano | Double-diffusive mixing in stellar interiors in the presence of horizontal gradients[END_REF]) the direction of the temperature gradient is not co-aligned with that of the local gravity. Then the interplay between stabilising density stratification and the destabilising thermal stratification can yield the emergence of double-diffusive convective instability [START_REF] Chen | Stability of thermal convection in a salinity gradient due to lateral heating[END_REF]. The cases where the initial temperature contrast has a perpendicular component to the gravity vector are particularly interesting for mid-scale ocean dynamics [START_REF] Shibley | Spatial variability of the Arctic Ocean's double-diffusive staircase[END_REF]. If planetary (or stellar) rotation is also taken into account, the complexity of the phenomena further increases and leads to nontrivial interactions between rotational and convective instabilities. In the present study, we focus on a rotating double diffusive system in a conceptual laboratory set-up that has been surprisingly rarely investigated in the literature.

A plethora of experimental works exists dealing with purely thermally-driven rotating flows (where the density of the fluid parcels is solely determined by their temperature) utilising lateral heating and cooling since the 1950s. The widely studied thermally driven rotating annulus, introduced by [START_REF] Hide | An experimental study of thermal convection in a rotating liquid[END_REF], [START_REF] Fultz | Studies of thermal convection in a rotating cylinder with some implications for large-scale atmospheric motions[END_REF], is a classic laboratory experiment designed to understand the manner in which the atmospheric circulation transports heat from equatorial to polar latitudes. The setup of this experiment consists of a tank with three concentric cylinders filled with pure water. The innermost cylinder (representing the polar region) is cooled, whereas the outer ring (representing the equator) is heated and the tank is mounted on a turntable, so it rotates around its own vertical axis of symmetry. Therefore, the working fluid in the annular cavity is subject to a radial temperature difference and a Coriolis deflection. The combined effect leads, for high enough values of the rotation rate, to the baroclinic instability with the formation of cyclonic and anticyclonic eddies covering the full water column. The baroclinic instability is at the heart of the mesoscale motions in the atmosphere and ocean and gives rise in particular to atmospheric weather systems [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]. Besides the baroclinic waves, varying either the magnitude of the temperature difference or the rotation rate, other flow regimes can develop in the annular gap. The four possible regimes are axisymmetric flow, steady waves, vacillation, and irregular flow.

Many theoretical [START_REF] Lindzen | Vacillations due to wave interference: applications to the atmosphere and to annulus experiments[END_REF], numerical [START_REF] Ohlsen | Nonlinear interference vacillation[END_REF]Hart 1989, Randriamampianina et al. 2006), and experimental [START_REF] Buzyna | Kinematic properties of wave amplitude vacillation in a thermally driven rotating fluid[END_REF][START_REF] Von Larcher | Experiments on transitions of baroclinic waves in a differentially heated rotating annulus[END_REF][START_REF] Harlander | PIV-and LDV-measurements of baroclinic wave interactions in a thermally driven rotating annulus[END_REF] studies have focused on the so-called amplitude vacillations, i.e., the periodic fluctuations in the intensity of the aforementioned baroclinic eddies, the associated temperature patterns and sought a fundamental understanding of the manner in which these vacillations arise, which are the processes leading to them and how they equilibrate by interacting with the mean flow. The findings of these studies suggest that amplitude vacillations mainly occur because of nonlinear wave-wave interactions and wave-mean flow interaction but also feedback mechanisms between the fluid interior and the boundary layers [START_REF] Früh | Amplitude Vacillation in Baroclinic Flows[END_REF]. A special case of a flow which appears like a modulated amplitude vacillation (also known as interference vacillation) is the superposition of two waves with the same horizontal wavenumber, non-orthogonal vertical structures, and different phase speeds [START_REF] Lindzen | Vacillations due to wave interference: applications to the atmosphere and to annulus experiments[END_REF][START_REF] Vincze | An experimental study of regime transitions in a differentially heated baroclinic annulus with flat and sloping bottom topographies[END_REF].

In addition to the baroclinic waves, inertial Kelvin modes are known to be resonantly excited in a rotating cylinder in various situations as an elliptic instability [START_REF] Lacaze | Elliptical instability in a rotating spheroid[END_REF]), precessional instability [START_REF] Lagrange | Precessional instability of a fluid cylinder[END_REF][START_REF] Triana | Precessional states in a laboratory model of the Earth's core[END_REF], and parametric forcing. The name Kelvin waves here refers to inertial waves modified by stratification having a specific vertical structure and these are not to be confused with coastal Kelvin waves, which are shallow water boundary trapped waves. [START_REF] Guimbard | Elliptic instability of a stratified fluid in a rotating cylinder[END_REF] have theoretically investigated the inertial Kelvin modes in a rotating cylinder and the behaviour of the dispersion relation when vertical salinity stratification is introduced. The occurrence of different waves and their interactions can lead to instabilities. For example, the resonance between Rossby and inertial Kelvin waves induces the Rossby-Kelvin instability, which has been studied numerically in the two-layer shallow water model by [START_REF] Gula | Instabilities of two-layer shallow-water flows with vertical shear in the rotating annulus[END_REF] and experimentally in a rotating annulus for a two miscible fluid layers by Flór et al. In the present work, we propose a thermohaline version of the differentially heated rotating annulus, the 'barostrat' experiment introduced by [START_REF] Vincze | Double-diffusive convection and baroclinic instability in a differentially heated and initially stratified rotating system: the barostrat instability[END_REF], where a continuously stratified salinity profile is prepared in the annulus cavity before the experiment is run with the so-called double-bucket technique [START_REF] Oster | Density Gradient Techniques[END_REF]. This particular configuration is not only interesting for investigating the wave regimes that develop in the fluid depth as well as studying inertial Kelvin/baroclinic wave coupling and the coexistence of different baroclinic waves, but also to investigate the occurrence of small-scale waves and in particular inertia gravity waves spontaneously emitted by the baroclinic wave. Recently, several numerical models [START_REF] Jacoby | Generation of inertia-gravity waves in the rotating thermal annulus by a localised boundary layer instability[END_REF][START_REF] Borchert | Gravity wave emission in an atmosphere-like configuration of the differentially heated rotating annulus experiment[END_REF][START_REF] Randriamampianina | Inertia-gravity waves in a liquid-filled, differentially heated, rotating annulus[END_REF] and laboratory experiments [START_REF] Lovegrove | Generation of inertia-gravity waves in a baroclinically unstable fluid[END_REF][START_REF] Williams | On the generation mechanisms of short-scale unbalanced modes in rotating two-layer flows with vertical shear[END_REF] have used different configurations of the rotating annulus experiment to investigate the interaction between largescale 'balanced' flow components (quasi-geostrophic baroclinic waves) and 'fast', small-scale, ageostrophic inertia-gravity waves.

The advantage of the laboratory set up used in this paper in contrast to the classical configuration is that, by introducing vertical salinity gradient to the set up, the frequency ratios between the Brunt-Väisälä (or buoyancy) frequency N = -g/ρ 0 dρ/dz -where ρ 0 is a reference density -and the Coriolis frequency f = 2Ω, where Ω is the angular velocity of the rotation, rises from N/f < 1 (classical set up) to N/f 10 (barostrat). Because the ratio is larger than unity, IGW propagation is expected to be qualitatively similar to the atmospheric case, as it follows from the dispersion relation: At the top and the bottom of the tank (green regions) N < f and therefore N < ω < f , while in the region in the middle (light blue) f < N and so, in this case, f < ω < N .

ω 2 = N 2 (k 2 + l 2 ) + f 2 n 2 k 2 + l 2 + n 2 = N 2 cos 2 (γ) + f 2 sin 2 (γ), (1)
The criterion for baroclinic instability reads [START_REF] Hide | Sloping convection in a rotating fluid[END_REF]:

Bu = N f D L 2 < 0.581, (2)
and we see that for large N/f the aspect ratio D/L needs to be small to find baroclinic instability. In the thermohaline annulus, this is automatically fulfilled since the convective layers are thin and hence D has to be replaced by the thickness λ of the convective layers.

The paper is organised as follows: in section 2 we briefly introduce the most important features of the experiment set up and the measurement techniques used, in section 3 we described the data analysis methods applied to our measurements. Section 4 first focuses on the large scale wave regimes and interactions occurring at different fluid heights in our tank and then investigates inertia-gravity waves along the jet of the baroclinic waves. Summary and conclusions are given in section 5.

Experimental set-up and data acquisition

The experiments considered in this paper were performed at the BTU Cottbus-Senftenberg laboratories using the same apparatus introduced in previous work by [START_REF] Vincze | Double-diffusive convection and baroclinic instability in a differentially heated and initially stratified rotating system: the barostrat instability[END_REF].

For a detailed description of the experimental set up, we refer the reader to their paper. In this section, we will only report the most significant parameters.

Experimental set-up

The annular tank consists of an inner cylinder, cooled by a thermostat, and an outer cylinder, where α is the thermal expansion coefficient and β is the haline contraction coefficient.

heated
Although we do not have measurements of the thermal gradients and thus no access to the Brunt-Väisälä frequency associated to the temperature stratification, we know that temperature vertical gradients only exists in the convective cells at the top and bottom of the water column, whereas in the bulk (i.e., outside the top and bottom cells) the thermal contribution to N is zero. In the convective cells, convection and the baroclinic instability mix the temperature field and consequently decreases the vertical difference of temperature in these layers. For this reason, we speculate that in the upper/lower mixed layer the value of N is rather low. Using the numerical simulations of Borchert et al.

(2014) ad specially their Figure 3c, the value for N should be constraint in the range 0.08 to 0.4

For our thermohaline experimental setup, it is useful to introduce a local version of the Taylor and thermal Rossby numbers:

T a(z) = 4Ω 2 L 5 ν 2 λ(z) , (3)
Ro T (z) = λ(z)gα∆T Ω 2 L 2 , ( 4
)
where Ω is the rotation rate, ∆T is the horizontal difference of temperature, L is the gap width, D is the total fluid depth, ν is the kinematic viscosity, g is the acceleration of gravity, and α the volumetric thermal expansion coefficient (whose value varies according to temperature, salinity, and pressure). These two nondimensional numbers, generally used to study the flow regime in the differentially heated rotating annulus, are now depending on the vertical extent of a convective cell at height z, λ(z), which is determined by the initial buoyancy frequency profile N (z) and the lateral temperature contrast ∆T . Note that, because the flow states depend not only on the lateral temperature difference and on the rotation rate but also on the convective cell thickness, different flow states may be observed at different heights. The condition that the initial saline density difference between the top and bottom of the cell cannot outreach the horizontal thermal density difference between the cylindrical walls define the thickness of the cell λ(z), which can, therefore, be calculated using Chen's formula: (b) buoyancy frequency calculated from the density measured after the rotation. The vertical red line shows the value of the Coriolis frequency f = 2Ω, the two green areas indicate the frequencies interval in which IGWs can be found (N < ω < f ) for the regions where the salinity stratification is weak. The central light blue area indicates the frequencies interval in which IGWs can be found (f < ω < N ) for the region with strong salinity stratification. In the latter case, IGWs are expected to propagate similarly to the ones in the atmosphere. (c) plot of the azimuthal velocities averaged over time for the four measured heights. It can be noticed that the zonal flows is prograde at the surface, then retrograde at z = 75 mm, corresponding to the convective region, almost zero in the middle motionless stratified layer, and again prograde at the lowest z.

λ(z) = gα∆T N 2 (z) . ( 5
(2016) found that the formation of two separate convective layers followed the Chen relation and the layers at the top and the bottom were observed to be around 2 cm thick.

PIV and temperature data acquisition

To investigate the flow regimes at different water heights (indicated by the dashed lines in Figure 2 and by the green laser planes in Figure 1), a co-rotating particle image velocimetry (PIV) system is mounted on the cylindrical annulus. A camera (GoPro Hero 4, screen resolution 1920 × 1080, fps = 30) together with a green laser (Laser Linos Nano 250-532-100, wavelength λ = 532 nm, maximum power 300 mW) are fixed on a vertical bar mounted at the outer cylinder. The laser produces a continuous horizontal light plane of 1 mm thickness which illuminates tracer particles (a mixture of hollow glass and silver coated hollow glass spheres with diameter range distribution 2 -20 µm and mean density ρ = 1.1 g/cm 3 ) in the fluid. Applying the Stokes' law, the settling velocity at which the particles move towards the middle region of the tank ranges between V t = 0.05 -4.7 cm per hour. The particles are added to the fluid when the density profile is prepared at the beginning of the experiment.

Considering that the experiment runs for maximum 6 hours, at least the smallest particles are expected to be found 1.5 cm below the surface and above the bottom of the tank at the end of it. In addition to the settling velocity, we have to take into account that the convective motions at the top and the bottom of the tank are mixing the PIV particles rather effectively.

Therefore, at the measurement heights the particles remain in neutral buoyancy during the measurement time 12 minutes.

The field of view recorded by the co-rotating camera, fixed at close range to the water surface, covers approximatively one third of the annulus and allows to gain close-ups of the velocity fields, giving a better resolution in particular on the small-scale features. The recorded videos, having each a total duration of 12 minutes, are processed with the free Matlab toolbox UVmat (see http://servforge.legi.grenoble-inp.fr/projects/soft-uvmat) and the horizontal components of the velocities are obtained. The percentage of velocity vectors excluded during the PIV analysis and then interpolated is less than 1.5%. The rms difference between the CIV (Correlation Image Velocimetry) and the smoothed velocity field is rms = 0.5 -0.6 pixels for a typical displacement of 5 pixels. This gives an estimation of the PIV error 10% -12%, comparable with the a posteriori error estimation (see Appendix A for more details).

In addition to the PIV system, five temperature sensors have been placed at the same levels for which we also do the PIV measurements, i.e., 94 mm, 75 mm, 47 mm, and 21 mm. The sensors are placed in the center of the annulus gap and diametrically opposite to the side where we do the PIV measurements. Hence the sensors, having a diameter of 500 µm and a

Reynolds number of the order Re = 1, do not affect the PIV measurements. The temperature sensors have a sampling rate ∆t = 1 s and are able to record the data for the entire duration of the experiment (about 2 hours). Thanks to the temperature sensors, we have long time measurements simultaneously at all the chosen fluid heights. A schematic drawing of the experimental set up with the temperature sensors and the PIV apparatus is given in Figure 1.

Data analysis methods

In this section we describe the two statistical methods, i.e., harmonic analysis and empirical orthogonal functions, we use to analyse our velocity fields. In particular, we are interested in studying whether the baroclinic dynamics are different in the two unstable layers.

The harmonic analysis is a form of signal demodulation in which the user specifies the frequencies to be examined and applies least square techniques to solve for the constituents, (for an extensive description of the method we refer to the textbook of [START_REF] Thomson | Data Analysis Methods in Physical Oceanography[END_REF]). This method is very useful to investigate the spatial patterns associated with a single given frequency and in particular we use it to emphasise the waves modes at different heights. Moreover, it has the advantage to be very robust even for short time series.

Besides harmonic analysis, we apply the empirical orthogonal function (EOF) analysis to the PIV data (a complete description of the method can be found in the textbook of [START_REF] Navarra | A guide to empirical orthogonal functions for climate data analysis[END_REF]). Without specifying particular frequencies in advance this method provides a description of the spatial patterns of variability of the data series and their temporal variation, breaking the data into orthogonal functions or 'modes of variability' and thus is widely used in geosciences (e.g., [START_REF] Lorenz | Empirical orthogonal functions and statistical weather prediction[END_REF]. One mode of variability can detect more than one frequency and therefore can comprise a more complex dynamics such as interactions between waves. The advantage of using both methods independently is that we can investigate time dependent patterns (with the EOF analysis, whereas the harmonic analysis gives only the spatial patterns) and give them a physical interpretation by comparing the results obtained with the two techniques. Indeed, while the interpretation of the spatial patterns obtained by the harmonic analysis is clear, connecting the EOFs to the physical modes is not trivial.

The emission of inertia-gravity wave packets from the baroclinic jet is a phenomenon highly localised in space and time [START_REF] Viúdez | Spontaneous generation of inertia-gravity wave packets by balanced geophysical flows[END_REF]. Therefore, it is very difficult to capture these small-scale waves using the two statistical methods described in this section. A quantity often used as an indicator for IGWs is the horizontal velocity divergence 'sullivan andDunkerton 1995, Borchert et al. 2014). The horizontal divergence contains a balanced part, as defined by quasi-geostrophic balance, and an imbalanced part that is related to IGWs. We use the horizontal divergence to study the small-scale waves, their behaviour in time, their wavenumbers and frequencies. The last ones are obtained by computing two dimensional fast Fourier transforms. Rodda, Borcia, Le Gal, Vincze, Harlander

∇ h •u = ∂u/∂x+∂v/∂y (O

Harmonic analysis

We consider the velocity field V , measured by PIV at an arbitrary grid point in our measured domain V (t n ). This quantity can be expressed by a Fourier expansion

V (t n ) = V + M q=1 [A q cos(ω q t n ) + B q sin(ω q t n )] + V r (t n ) (6)
with q = 0, 1, ..., M and M is the number of distinct frequencies to be analysed, V is the temporal mean, V r the residual of the time series (it could contain other kinds of components), t n = n∆t the time, ω q a constant frequency, A q and B q are the harmonic coefficients of the Fourier series. The amplitude of the frequency component q is

C q = (A 2 q + B 2 q ) 1/2 , ( 7
)
and the phase is

φ q = tan -1 (B q /A q ) (8)
The M distinct frequencies to be analysed are chosen as the main peaks in the horizontal velocity spectra at each measured fluid height. The horizontal velocity components, measured along one radial line taken in the middle of the camera field of view, are selected. The frequency spectrum for each of these points is then calculated by using a fast Fourier transform algorithm and then the spectrum is averaged for the points of this particular chosen line.

A truncated Fourier series is thereafter fitted to the time series containing the frequencies to be analysed, ω q . The variance e 2 is computed for each point

e 2 =    V (t) -   V + M q=1 [A q cos(ω q t) + B q sin(ω q t)]      2 (9)
where V (t) are the measured PIV velocities, V is the temporal mean and A q , B q are the Fourier coefficients of the harmonics obtained by a least-square fit.

We estimate the amplitudes, C q , and phases, φ q , of the various components by minimising the variance. Once the amplitudes and the phase for a certain frequency are calculated, the corresponding velocity field can be plot on the domain recorded by the camera (approximately one third of the tank). For rather steady waves, we can graphically reconstruct the entire annulus using symmetric properties. The reconstructed plots are made by combining together 3 or 4 partial plots obtained from the harmonic analysis. These are shifted by a phase φ = 120 • , when 3 images are used or φ = 90 • when 4 images are used. The figures have then been combined together (with some overlapping) to reconstruct the full annulus, relying on the hypothesis of patterns regular in space.

Empirical orthogonal functions

To analyse data that contain oscillations in time or in space and time as a propagating signal, it is useful to use a modified version of the standard EOF analysis, the so called complex empirical orthogonal functions (CEOFs) [START_REF] Pfeffer | A study of baroclinic wave behavior over bottom topography using complex principal component analysis of experimental data[END_REF]). To calculate the CEOFs we follow the same approach described in [START_REF] Harlander | PIV-and LDV-measurements of baroclinic wave interactions in a thermally driven rotating annulus[END_REF], considering the CEOFs method for a simultaneous analysis of more than one field. Such a coupled analysis is useful in our case where one field is the u-and the other is the v-component of the velocity. We proceed in the following way to find coupled propagating patterns: complex time series of the velocity components are formed from the original time series and their Hilbert transforms:

u c (x, t) = u(x, t) + iu H (x, t), (10)
v c (x, t) = v(x, t) + iv H (x, t), (11)
where u(x, t), v(x, t) are the time series of the horizontal components of the velocity measured by PIV at each location in the recorded domain, u H (x, t) and v H (x, t) are the Hilbert transforms of u(x, t), v(x, t) and i is the imaginary unit.

The Hilbert transform of the original time series X t with Fourier decomposition

X t = ω ζ(ω) exp(-2πiωt), (12)
is defined as:

X H t = ω ζ H (ω) exp(-2πiωt), (13)
where

ζ H (ω) = iζ(ω) for ω ≤ 0 and ζ H (ω) = -iζ(ω) for ω > 0.
For the computation of the Hilbert transforms, we used the algorithm described by [START_REF] Marple | Computing the discrete-time "analytic" signal via FFT[END_REF].

Thereafter, we form extended time series by combining the v c (x, t) time series with the u c (x, t)

time series and rewriting them as a row-vector in the form of U = (u 1 , ..., u M , v 1 , ..., v M ). The data matrix D is written in the form of

D =      U 1 (t 0 ) U 2 (t 0 ) . . . U M (t 0 ) U 1 (t 1 ) U 2 (t 1 ) . . . U M (t 1 ) . . . . . . . . . U 1 (t N ) U 2 (t N ) . . . U M (t N )      (14)
where the rows represent the state vector at the spatial grid points

U (t n ) = (U 1 (t n ), ..., U M (t n )
at time t n and the columns represent the time series U m (t) = (U m (t 0 , ..., t N )) T at the spatial point m.

From the data matrix, we calculate the covariance matrix

F = D T D, (15)
The CEOFs are the eigenvectors of F and the so called Principal Components (PCs) are the corresponding time-dependent coefficients.

Results

Propagating waves at different fluid heights

A puzzling result from the experiments done by [START_REF] Vincze | Double-diffusive convection and baroclinic instability in a differentially heated and initially stratified rotating system: the barostrat instability[END_REF] was that baroclinic waves have been found only in the upper convective layer (see Figure 1). The reason for this might be damping due to bottom Ekman layer effects. In order to obtain a clear picture of the waves developing at the different fluid heights in the tank, we first consider frequency spectra of the horizontal velocity, obtained as described in Section 3.1, from which we can identify the dominant frequencies at each height, (i.e., measured from the bottom of the tank: 94 mm, 75 mm, 47 mm, and 21 mm).

We have separated the spectra into two windows: 0 < ω < f (Figure 3 left) and f < ω (Figure 3 right). The grey dashed lines in the figures indicate the value of the buoyancy frequency N at the measurements heights where known. Peaks for low frequencies appear in the upper layer, for the frequency ω = 0.03Ω, at heights z = 94 mm (stronger) and z = 75 mm (weaker) and in the bottom layer (z = 21 mm), for the frequency ω = 0.14Ω. In the middle layer (z = 47 mm), instead, there is no peak in the low frequency range. Another peak, at ω = Ω, can be seen for all the layers, with a stronger signal at the bottom of the tank. Finally, in all layers there are peaks at frequencies larger than f . Using the harmonic analysis as described in section 3.1, we can reconstruct the velocity fields corresponding to the most prominent peaks (Figure 5). As we mentioned before, besides baroclinic waves also other waves can be expected to rise in our experimental setup. The middle column of Figure 5 shows the velocity pattern corresponding to the peak at frequency ω = Ω. We identify this pattern as inertial Kelvin waves modified by stratification. There are several examples in the literature of inertial Kelvin mode, also called the spin-over inertial mode, driven by precession either due to the rotation of the laboratory by the Earth (i.e., in a spherical shell [START_REF] Triana | Precessional states in a laboratory model of the Earth's core[END_REF])) or in a precessional cylinder (i.e., [START_REF] Lagrange | Precessional instability of a fluid cylinder[END_REF]. Moreover, inertial Kelvin waves might be excited by elliptical instability, as reported by [START_REF] Lacaze | Elliptical instability in a rotating spheroid[END_REF]. In our case the inertial Kelvin mode is more likely due to imperfect alignment of the rotation axis with respect to gravity. In a Rodda, Borcia, Le Gal, Vincze, Harlander λ (cm) α (10 system with a free surface this imperfection might force an inertial Kelvin mode even if a true precession of the rotation axis is absent (personal communication with Patrice Meunier).

Due to the fact that our measurements are not simultaneous, we do not have phase information and hence we cannot investigate in detail the vertical structure of the inertial Kelvin wave and compare it with the analytical solutions found by [START_REF] Guimbard | Elliptic instability of a stratified fluid in a rotating cylinder[END_REF] for a rotating stratified cylinder. Nevertheless, our experiment shows that the amplitude of the inertial Kelvin waves varies with depth as can be noticed in Figure 3 resulting from resonances between Rossby and inertial Kelvin eigenmodes for a two layer flow in a rotating annulus. However, no clear signal of such resonance can be observed in our experiment, in particular since the mean Rossby number is too small and this resonance just happens for large Ro although some interactions are detected.

In the right column of Figure 5 the reconstructed velocity field for one of the main peaks with frequencies larger than f , for instance ω = 6.8Ω, is shown. Similar structures can be found for the other prominent peaks for frequencies ω > Ω and are present in all layers. It is also interesting to notice that most of these peaks appear at the same frequency, not changing with the height. Obviously, the horizontal structure is similar to the inertial Kelvin wave, but in this case the amplitude does not decrease in the direction of the outer wall and the wave propagates prograde. Remarkably, we found very similar high frequency waves in an experiment carried out with the same set up, but where we only rotate the cylindrical tank without any lateral temperature difference and vertical salinity stratification. We speculate that these waves with frequencies ω > f might be surface wave modes of the Poincaré type.

They might be excited by a weak sloshing at the free surface and, in the case of the barostrat experiment, also at the interface between the layers of different density. We do not further investigate these weak gravity wave modes but will focus on the signature of frontal gravity waves that are very localised in space and time. Such localised wave packets move with the baroclinic jet and Fourier analysis of local time series is hence not a proper tool to detect those waves. More details on the IGW field related to the baroclinic front shall be given in Section 4.3 after discussing the baroclinic waves in the next section.

Determination of the baroclinic wave dynamics in the two unstable layers

In this section we investigate in more detail the similarities and differences between the baroclinic waves appearing in the top and bottom layer. To do so, we performed a CEOF analysis (see section 3.2) of the measured horizontal velocity field. Figure 6 shows the first ten eigenvalues of the covariance matrix in the uppermost layer in the tank (Figure 6(a)) and in the bottom layer (Figure 6(b)). Each eigenvalue stands for the variance that can be explained by the corresponding complex function, the CEOF.

For the top layer (z = 94 mm), the first and the second eigenvalues, explaining together 85% of the total variance, are related to the baroclinic mode and its first harmonic (a power spectrum of the PCs shows one single peak for each PC at ω = 0.04Ω and ω = 0.08Ω respectively, these peaks correspond to the ones detected previously by the FFT analysis).

The third eigenvalue, explaining 3.4% of the variance, corresponds to the second harmonic of the baroclinic wave and the inertial Kelvin wave (two peaks in the PC spectrum at ω = 0.12Ω

and ω = Ω). The fourth eigenvalue, explaining 2.6% of the variability, is related to the inertial Kelvin mode.

The eigenvalue spectrum for the bottom layer (z = 21 mm) is quite different as can be seen in Figure 6(b). Also in this case, the first eigenvalue is related to the dominant baroclinic mode, but it explains only 39% of the total variance and the power spectrum of its correspondent PC shows a multitude of peaks with frequencies very close to the ones shown in Figure 9(a), and in addition a peak for ω = Ω, suggesting a complex dynamics and interactions among waves. We shall discuss these interactions in more detail in the following text. The second eigenvalue, 29% of the variance, is related to the inertial Kelvin mode interacting with the baroclinic wave (the PC spectrum shows a peak for ω = 0.14Ω and ω = Ω). The third and the fourth, explaining 4.8% and 3.4% of the total variance, show a broad spectrum of frequencies and are not easy to interpret as physical modes. The fifth and sixth eigenvalues, explaining 2.6% and 1.8% of the variability, are related to the second harmonic of the inertial Kelvin mode ω = 2Ω. With respect to the eigenvalues of the covariant matrix for the uppermost layer (Figure 7(a)) we have a significant reduction of the first eigenvalue. Moreover, an interaction between the baroclinic wave and the inertial Kelvin wave seems to play an important role in the dynamic of this layer, as both frequencies are contributing to the second CEOF.

A similar behaviour, for which the dominant component is smaller by about 20% compared to the steady wave regime has been observed by [START_REF] Hignett | Characteristics of amplitude vacillation in a differentially heated rotating fluid annulus[END_REF] in case of an amplitude vacillating regime (AV).

If we now consider the first CEOFs for the top layer (Figure 7(c)) and for the bottom layer (Figure 7(d)) and their respective PCs (Figure 7(a) and 7(b)) we can see that in the top layer the fluid is in a steady wave regime, where the baroclinic wave shows a regular behaviour and has a very low phase speed (0.0057 rad/s), while in the bottom layer the baroclinic wave shows a variation of the amplitude in time with a phase speed of 0.0157 rad/s. Usually a metastable transient AV has a very high phase speed, circa 5 times faster than the finally equilibrated flow [START_REF] Früh | Wave interactions and the transition to chaos of baroclinic waves in a thermally driven rotating annulus[END_REF]. This supports our assumption that the flow in the bottom layer is in the vacillation regime. For AV the amplitude of the wave varies periodically while the shape of the wave pattern remains constant. The strength of the vacillation is characterized by a vacillation index I v defined over one vacillation cycle as

I v = A max -A min A max + A min . (16)
Hignett (1985) indicates as critical vacillation index I v = 0.05. This value separates the flow to be in steady regimes (I v < 0.05) and amplitude vacillation regimes (I v > 0.05). For the bottom layer of our experiment the value for the vacillation index is I v = 0.19, while for the top layer I v = 0.03.

For a steady wave regime it has been observed in many baroclinic wave experiments (e.g., [START_REF] Hide | Thermal convection in a rotating fluid subject to a horizontal temperature gradient: spatial and temporal characteristics of fully developed baroclinic waves[END_REF][START_REF] Hignett | Characteristics of amplitude vacillation in a differentially heated rotating fluid annulus[END_REF][START_REF] Früh | Wave interactions and the transition to chaos of baroclinic waves in a thermally driven rotating annulus[END_REF]) that the amplitude spectrum is composed almost entirely of the dominant component and its harmonics. This is also the case for the upper layer of our barostrat experiment. The peaks corresponding to the dominant wavenumber m = 3 and its harmonics are dominating the low frequency spectrum plotted for temperature measurements at z = 94 mm in Figure 8(a). In Figure 8 the reconstructed waves, baroclinic, first and second harmonics with frequencies and wavenumbers ω t 0 = 0.04Ω, experiments the second harmonic was strongly pronounced, in our case the fundamental wave and the first harmonic are the most energetic ones. Moreover, it has been shown by [START_REF] Früh | Wave interactions and the transition to chaos of baroclinic waves in a thermally driven rotating annulus[END_REF] that for steady waves the harmonic triad shows a more persistent and pronounced phase locking than the long wave triads. [START_REF] Buzyna | Kinematic properties of wave amplitude vacillation in a thermally driven rotating fluid[END_REF] have represented a particular case of amplitude vacillation regime as due to an interference of two waves with the same azimuthal wavenumber and different phase speed. This regime is referred to as interference vacillation. The basis of this approach has been that a travelling, modulated wave can be represented as a linear superposition of two waves of the same azimuthal wavenumber, but with different phase speeds or frequencies. A second nonlinear interaction can be identified in our experiment: the two baroclinic waves wave and its first and second harmonics: ω t 0 = 0.04Ω, ω t 1 = 0.08Ω and ω t 2 = 0.12Ω; the velocity fields corresponding to these most energetic modes reconstructed using the CEOF analysis are plotted in (b),(c) and (d) respectively. The three waves of frequencies ω t 0 , ω t 1 and ω t 2 (at z = 94 mm) form a harmonic triad with azimuthal wavenumbers m t 0 = 3, m t 1 = 6, and m t 2 = 9.

m t 0 = 3, ω t 1 = 0.08Ω, m t 1 = 6,
with them, as already discussed. [START_REF] Früh | Wave interactions and the transition to chaos of baroclinic waves in a thermally driven rotating annulus[END_REF] also observed a complex interaction scenario in amplitude vacillation regimes, and found that in this case long wave triads are usually observed.

It is interesting to notice that interference vacillation has been observed by [START_REF] Harlander | PIV-and LDV-measurements of baroclinic wave interactions in a thermally driven rotating annulus[END_REF] in the classical configuration of the thermally driven annulus. However, in their case there was no indication that the two waves were coupled through nonlinear interactions but they appeared to be a linear superposition of two modes of different zonal wave number drifting at different speeds. In our thermohaline version of the experiment, on the contrary, the amplitude vacillation results from linear interactions between two waves having the same wave number leading to a mean zonal flow and nonlinear interactions between the two waves and the mean zonal flow, more in agreement to the results from [START_REF] Buzyna | Kinematic properties of wave amplitude vacillation in a thermally driven rotating fluid[END_REF].

Inertia-gravity waves

So far we focused on the large-scale modes and their interactions. However, besides the already discussed baroclinic waves and inertial Kelvin waves, it is instructive to investigate IGWs occurrence at the baroclinic wave fronts. Such waves are very localised in space and time. In the paper by [START_REF] Vincze | Double-diffusive convection and baroclinic instability in a differentially heated and initially stratified rotating system: the barostrat instability[END_REF] wave trains whose characteristics are compatible It can be noticed that for the shown layer most of the divergence signal is associated to the small-scale waves. Moreover, we can see a wave train structure clearly related to small-scale phenomena attached to the baroclinic jet. Similar structures are visible in most of our data, embedded in the baroclinic wave and showing up intermittently.

It is instructive to further investigate how these small-scale waves travel with respect to the baroclinic jet. The reason is that observations in the atmosphere by Uccellini and The phase velocity of the small-scale waves, measured from this plot and indicated by the blue line is 3 mm s -1 , the same as u 0 .

Koch (1987) identified intense low frequency IGW in the jet exit regions and hypothesised spontaneous imbalance as source mechanism. Other observations by [START_REF] Fritts | Sources of mesoscale variability of gravity waves. Part II: Frontal, convective, and jet stream excitation[END_REF] have highlighted a conspicuous enhancement of gravity wave activity in the vicinity of jets and fronts and our experimental study seems to be in line with these numerical and observational findings. To examine these regions, at z = 94 mm, we took two different cuts, one parallel to the entrance of the baroclinic wave (Figure 11) and one parallel to the exit of the baroclinic wave (Figure 12). Because we record our data in the system of reference co-rotating with the tank and the baroclinic wave moves prograde, the wavefront is crossing our line just for a short time, as it can be seen in Figures 11(a), 11(b), 12(a), and 12(b).

Hence, the most prominent travelling waves can be found in the period 10 s < t < 20 s in For a quantitative comparison, we plotted in Figure 11(c) and 12(c) the drift speed of the baroclinic wave (dashed line) and the mean value of the zonal flow in the jet region (red line).

Remarkably, the phase velocity of the small-scale waves (blue lines) is in both cases similar to the one of the zonal flow. Compared to the baroclinic wave, the wave packets move faster, with phase speed equal to 3 mm s -1 at the entrance and 5 mm s -1 at the exit region of the baroclinic wave respectively, but are attached to it. Rodda, Borcia, Le Gal, Vincze, Harlander The phase velocity of the small-scale waves, measured from this plot and indicated by the blue line is 5 mm s -1 .

Subsequently, we want to verify whether the intrinsic frequency of the measured small-scale waves is consistent with IGW dispersion relation. The intrinsic frequency ω i is the frequency of a wave relative to the flow, i.e., the wave frequency measured by an observer drifting with the mean flow:

ω i = ω m -u 0 k, (17)
ω m is the frequency measured in the co-rotating system of reference, u 0 is the wind speed and k the wave vector along the jet.

For the uppermost layer (z = 94 mm), considering a zonal flow in the jet u 3 mm s -1 , consistent with the PIV measurement, wavevector k = 2π/λ 3 cm -1 and the measured frequency ω m = 1/6 s -1 = 1.04 rad s -1 the obtained intrinsic frequency is ω i 0.1 rad s -1 .

The finding of a frequency lower than f might seem surprising, however we recall that at this measurement height the salinity stratification is weak, as we already have shown in the green areas in Figure 2(b), where we can notice that the value of N due to the salinity at z = 94 mm is zero. As we previously pointed out, the contribution of temperature stratification to the value of N close to the water surface is not measured experimentally and thus the gradients are difficult to estimte. However, keeping in mind that convection and baroclinic instability imply a strong mixing, added to the fact that heat losses at the surface weaken the temperature gradients in the top layer, it is not unreasonable to consider N < 0.1 rad/s and hence the waves inside the IGW frequency range, although we cannot prove it firmly:

| N |< ω 0.1 rad s -1 < 0.836 rad s -1 =| f |, (18)
From Figure 11(c) and 12(c) we realised that the wave packet travels with the jet. This is further in agreement with our finding for the frequency, where ω m = 1.04 rad s -1 u 0 k = 0.94 rad s -1 that suggests small-scale waves travelling with the jet.

Despite the fact that the waves observed in the experiment are not plane and the dispersion relation ignores latitudinal and vertical shear, using the intrinsic frequency we can get an estimation of the vertical wavelength using the dispersion relation solved for λ z :

λ z = λ 2 x λ 2 y (λ 2 x + λ 2 y ) (ω 2 i -f 2 ) (N 2 -ω 2 i ) . ( 19
)
The estimated vertical wavelength is therefore λ 94 z 5 cm (for N = 0) this order of magnitude being compatible with the thickness of the convective layer.

Wave trains at z = 75 mm

We continue the discussion on the small-scale waves signal in this section considering now the data at the measurement height z = 75 mm. We recall that this height is in the stratified region where N > f , for this reason IGWs are expected to show similarities with atmospheric gravity wave packets; the frequencies range for IGWs at this height is highlighted by the light blue coloured central region in Figure 2(b). Furthermore we want to recall that the baroclinic wave with m = 3 is present at z = 75 mm even though it is weaker and drifts in the opposite direction with respect to the baroclinic wave observed at the height z = 94 mm (see Figure 5).

We follow the data analysis presented in the previous subsection discussing the analogies and the differences between the waves observed at the two heights.

Figure 13 shows the plot of the horizontal divergence (a) in a snapshot at the time t = 706 s where a wave train can be seen, similarly to the one found at z = 94 mm (Figure 10(a)).

The related signal-to-noise ratio, calculated as described in Appendix A, is plotted in Figure as

Ro L = U f L , (20)
where A complete analysis of the generation mechanisms of the short scale waves observed in our experiment is beyond the purpose of this paper. Nevertheless, some possible wave excitation mechanisms can be investigated with the available data. One is excitation due to convection, Rodda, Borcia, Le Gal, Vincze, Harlander and this can be excluded since the gravity waves have not been found in our experimental setup without rotation, i.e., without the baroclinic jet. Kelvin-Helmholtz and Hölmböe instability are also mechanisms that might be responsible for small-scale waves emission. To distinguish between these two instabilities, it is necessary to know the interface thickness in shear and density. When their ratio exceeds a threshold usually considered being equal to 2, or more in general in the condition where a thick shear layer and a sharp density interface develop, Hölmböe instability can occur. A detailed analysis of these two instabilities and on other small-scale waves possible generating mechanism can be found in the paper by Scolan et al.

U = √ u 2 + v 2 is
(2014). Because we do not know precisely the thickness of the interfaces, we can not distinguish between these two shear instabilities. Nevertheless, we can estimate the Richardson number

Ri = N 2 (du/dz) 2 , ( 21
)
where u(z) is the horizontal velocity profile, and check whether the Miles criterion for stratified shear instability is satisfied. Ri ≥ 1/4 is the sufficient condition for the stability of a sheared, does point to spontaneous imbalance as generating mechanism. However, even if we are fully confident for the results at z = 75 mm, our analysis is rather coarse compared to the one done for the two-layer case by [START_REF] Scolan | Frontal Instabilities at Density-Shear Interfaces in Rotating Two-Layer Stratified Fluids[END_REF] and further investigations are needed to confirm the generation mechanism of the small-scale waves in the continuously stratified case.

Finally, we want to mention that optical deformation of the free surface due to capillary waves would have a much larger phase speed than the waves we observed. We conclude that capillary waves are not responsible for the small-scale structures at the free surface.

Conclusions

The differentially heated rotating annulus is a classical laboratory experiment to study largescale baroclinic waves and their interactions. Here, in contrast, we focused on the existence and interaction of different wave types in a modified version of the classical experiment: instead of stratification due to differentially heating the lateral walls alone we added salt stratification in the axial direction and hence considered a thermohaline version of the experiment. The novel feature in this configuration is the occurrence of thin convectively mixed layers at the top and the bottom where baroclinic instability can take place and a thicker stable stratified, dynamically less active layer in between.

For this configuration, we experimentally observed two different baroclinic waves with azimuthal wavenumber 3 for the top and 4 for the bottom layer. The former was in a steady state with a slow drift rate exciting higher harmonics by self-interaction. The latter showed amplitude vacillations and a fast drift. The vacillations result from the interaction between two m = 4 waves with different phase speeds as was described earlier by [START_REF] Buzyna | Kinematic properties of wave amplitude vacillation in a thermally driven rotating fluid[END_REF] for a classical annulus experiment. Though nonlinear triadic interaction between the waves in the individual layers can clearly be seen in the spectra it is not clear yet whether the baroclinic waves in the surface and bottom layers are coupled. What is striking is the close correspondence of some prominent peaks in the spectra: the wave peak in the top layer is very close to the difference of the frequencies of the two dominant wavenumber 4 modes in the bottom layer implying a 'frequency triad' between top and bottom layer waves. Moreover, the two harmonic frequencies of the upper layer wave nearly correspond to the frequency of the two wavenumber 4 waves in the bottom layer.

Besides the large-scale baroclinic waves, we further observe an inertial Kelvin type global mode and a higher frequency surface wave of Poincaré type. The inertial Kelvin mode shows similarity with a shallow water boundary trapped Kelvin wave, has an azimuthal wavenumber 1 but has the frequency of the tank's rotation and not f . It is hence not a shallow water mode but an inertial wave Kelvin mode with a certain vertical structure. Unfortunately, since we did not measure simultaneously at different vertical levels we do not know the vertical wave number. We just know that the inertial Kelvin mode is strongest in the non-convective zone and the bottom layer. The mode is very likely mechanically driven by a slight deviation of the rotation axis from vertical.

The origin of the Poincaré type modes with frequencies larger than f is not yet clear.

They might also be triggered mechanically by a very weak sloshing of the tank. A theoretical treatment of those modes (as done by [START_REF] Mougel | Waves in Newton's bucket[END_REF] for the homogeneous case) is hampered by the nonlinear vertical density profile. For a thorough analysis, effects at the internal interfaces have to be considered in addition to the surface wave modes. Such a study is postponed to the future when more data, in particular in the high frequency range and from vertical cross sections are available.

A long standing problem is whether in differentially rotating annulus experiments instabilities of baroclinic fronts and in particular spontaneous imbalance of the frontal flow can be observed. Short interfacial wave patterns have been observed in experiments by [START_REF] Lovegrove | Generation of inertia-gravity waves in a baroclinically unstable fluid[END_REF] and [START_REF] Williams | On the generation mechanisms of short-scale unbalanced modes in rotating two-layer flows with vertical shear[END_REF], even though the source for these patterns is not fully clear yet [START_REF] Flór | Frontal instabilities and waves in a differentially rotating fluid[END_REF]. In contrast to the global modes described above such short-wave inertia-gravity waves are hard to be found in the wave spectra since they typically occur sporadically and are spatially localized since they are tightly connected to the baroclinic front.

Although spontaneous imbalance has been documented from a number of numerical simulations also for the annulus configuration [START_REF] Borchert | Gravity wave emission in an atmosphere-like configuration of the differentially heated rotating annulus experiment[END_REF] and [START_REF] Hien | Spontaneous inertia-gravity wave emission in the differentially heated rotating annulus experiment[END_REF]) for which Reynolds number like in the experiment was rather small (Re 100) it has never been experimentally observed for the continuously stratified annulus for which in general the Rossby numbers are rather small and the flow is in a nearly balanced state. In our experiments the top layer shows wave packets in the horizontal divergence field travelling mainly with the frontal mean flow and hence much faster than the drift speed of the baroclinic waves. These structures are not unlike the wave trains described by [START_REF] O'sullivan | Generation of inertia-gravity waves in a simulated life cycle of baroclinic instability[END_REF]. We have shown that their signal is above the noise level, their frequencies are in the inertia-gravity wave range and that the frequency-wavenumber plot roughly follows one branch of the dispersion curves. We speculate that compared to the classical "Hide" experiment its thermohaline version is more suitable to find such waves. However for sound evidence, whether the waves in the mixed baroclinically unstable layer are generated by shear instability or spontaneously, we need data on the temperature related N . Then we can also determine N/f in the baroclinic layers which is larger than one for the atmosphere but smaller than one for the "Hide" experiment. For the classical annulus, inertia is dominant and the baroclinic fronts do occur over the full depth of the tank and not just in thin layers. The latter might be more favourable to generate frontal waves. Moreover, we want to mention that in contrast to the classical setup, no numerical simulation is available for the barostrat case.

We have experimentally demonstrated that the thermohaline version of the differentially heated rotating annulus is not only a testbed to study large-scale wave interactions but is also a setup suitable to study interactions of different wave types. The natural layer formation offers the possibility to investigate wave resonance but also wave propagation between the different layers. The dynamically passive layer in the annulus center might be seen as a model for the stratosphere comprised between baroclinically unstable tropospheres as was proposed by [START_REF] Vincze | Double-diffusive convection and baroclinic instability in a differentially heated and initially stratified rotating system: the barostrat instability[END_REF].

  where the intrinsic frequency in the rotating frame of reference ω of the wave is determined by the buoyancy frequency N , the Coriolis frequency f , and γ = arctan(n/ √ k 2 + l 2 ), the angle between the phase velocity and the horizontal plane, set by horizontal wavenumbers k, l and vertical wavenumber n. Due to the finite domain of values of the sine and cosine functions in (1), IGWs can only exist for specific frequencies ω. More specifically, in Figure 2(b), which represents the buoyancy frequency profile N (z) as calculated from the measured density profile in panel (a), the two distinct existence regions, according to the values of f and N , are shown.

  by a heated coil, with radii equal to a = 4.5 cm and b = 12 cm respectively. The annular cavity of width L = ba = 7.5 cm has a flat bottom, a free surface, and was filled up to the height D = 10.5 cm. The lateral temperature difference ∆T was set to 10 K, after this value was reached in the experiment, the rotation rate of the turntable about the vertical axis was slowly increased (∆Ω = 0.1 rpm every 2 minutes) until the final value Ω = 4 rpm. A sketch of experimental set-up with the double diffusive convection developing in the barostrat experiment is shown in Figure1.Applying the standard two-bucket technique[START_REF] Oster | Density Gradient Techniques[END_REF] a salinity profile was created in the annular gap. The prepared stable vertical salinity profiles have been measured with a conductivity meter at the beginning of the experiment, before starting the rotation and at the end of the experiment, after stopping the rotation. Successively the conductivity has been converted to density. The calculated vertical density profiles are plotted in Figure 2(a), from which we can distinguish three layers in the tank: two shallow layers (whose thickness, indicated by the grey shading, is approximatively two centimeters) on the top and the bottom of the tank, both characterised by constant density, and a middle layer with an almost linear density profile. The dinstinct interfaces separating the regions, visible in the plot, correspond to the frontier of the convective layers. The plot of the time averaged azimuthal velocities for the four measured heights (Figure 2(c)) shows that the zonal flow is prograde at the surface (z = 94 mm), then retrograde at z = 75 mm, almost zero in the middle motionless stratified layer (z = 47 mm), and again prograde at the lowest measured height (z = 21 mm). The plots for the density, the buoyancy frequency N , and the azimuthal velocity (Figure 2) are in eccellent agreement with the ones shown by Vincze et al.(2016) for similar experimental settings. The boundary conditions of no-flux for salinity at the surface and the bottom of the water column imposes a zero concentration gradient at these boundaries even before a temperature gradient is applied. In this configuration, only the two separate shallow fluid layers can be baroclinically destabilised. Indeed, the vertical salt stratification opposes the thermal convective motions until the ratio of the (horizontal) thermal density difference and the (vertical) salinity-induced density difference exceeds a certain critical threshold. Then double-diffusive convection rolls develop in thin layers located in regions where the salt stratification is weakest. An important clarification is needed: the curves plotted in Figure2for the density and the buoyancy frequency N are obtained measuring the temperature compensated conductivity, so they depend only on the variations of the salinity content in the water column. However, density variations depend both on the salinity (δS) and the temperature (δT ) variations in the form δρ/ρ 0 = -αδT + βδS,

Figure 1 .

 1 Figure1. Sketch of the double-diffusive convection in the barostrat experiment and of the experimental set-up with PIV levels and temperature sensors positions. The red dashed lines indicate the interfaces between the convective regions and the calm region at the heights ztop = 8.9 cm and z bottom = 2 cm.

)Figure 2 .

 2 Figure 2. Vertical density (a) buoyancy frequency (b) and azimuthal velocity (c) profiles. The dashed lines correspond to the PIV measurement heights, the grey shaded areas indicate the convective layers at the top and the bottom of the tank. (a) the density is measured before starting rotation (red line) and after stopping the rotation (black line).(b) buoyancy frequency calculated from the density measured after the rotation. The vertical red line shows the value of the Coriolis frequency f = 2Ω, the two green areas indicate the frequencies interval in which IGWs can be found (N < ω < f ) for the regions where the salinity stratification is weak. The central light blue area indicates the frequencies interval in which IGWs can be found (f < ω < N ) for the region with strong salinity stratification. In the latter case, IGWs are expected to propagate similarly to the ones in the atmosphere. (c) plot of the azimuthal velocities averaged over time for the four measured heights. It can be noticed that the zonal flows is prograde at the surface, then retrograde at z = 75 mm, corresponding to the convective region, almost zero in the middle motionless stratified layer, and again prograde at the lowest z.

  For a harmonic wave of the form V (x, t) = Re[U (x) exp(-iωt)] a peculiar phase relation that indicates propagation is a quarter wavelength shift. CEOF analysis enhances this phase relation changing the available data by adding a new data set obtained by shifting all data by one quarter wavelength by using a Hilbert transform (Navarra and Simoncini 2010). Therefore, March 1, Geophysical and Astrophysical Fluid Dynamics submission˙2˙RoddaBW2 Baroclinic, Kelvin and inertia-gravity waves a single CEOF represents a single mode split into two patterns with a phase difference of π/2 (the real and imaginary part of the CEOF).

Figure 3 .

 3 Figure3. Velocity frequency spectra for the different heights. On the left the low frequencies, 0 < ω < f , on the right heigh frequencies, f < ω. Note that the spectra amplitude on the left plot differs from the one on the right. The vertical dashed line indicates the value of N .

  The left column shows the baroclinic waves, one in the uppermost layer (z = 94 mm and z = 75 mm) close to the surface, and one in the bottom layer (z = 21 mm), with different azimuthal wave numbers, m = 3 and m = 4 respectively.It is instructive now to use the T a-Ro T regime diagram to compare the flow regimes observed for the two baroclinically unstable top and bottom layers with data from the literature for the nonsalt stratified, purely thermal configuration. Once baroclinic instability has set in, for the experiments where pure water is used, the convective cell covers the whole fluid depth. For the comparison we consider two separate data set collected at the BTU laboratory using the same experimental apparatus described in Section 2. The first data set is an investigation of the flow regimes in a broad range of the T a-Ro parameter space by Von Larcher and Egbers(2005). We reproduced the regime diagram proposed in their paper in Figure4, where the different flow regimes are indicated by the lines and the wavenumbers are given. The second data set here considered is a study by[START_REF] Vincze | Benchmarking in a rotating annulus: a comparative experimental and numerical study of baroclinic wave dynamics[END_REF] which investigate the properties of the baroclinic instability by considering different initial conditions (i.e., the lateral difference of temperature was kept constant whilst the rotation rate was increased (spin up condition ↑) or decreased (spin down condition ↓ Table1) between each measurement). The four cases chosen for comparison are plotted (black-diamond shaped marker) in Figure4and indicated by the letters A-D. Our barostrat experiment data for the top and bottom layer are plotted in red and blue dots respectively. The values of T a and Ro T calculated using (4) and (3), the thickness of the layer (λ) or the total fluid depth, and the observed wavenumbers are listed in Table1for the second data set and the barostrat top and bottom layers. One can clearly see that the thin layers of the thermohaline experiment data are located in a region where baroclinic waves with m = 3 and m = 4 have both been observed in the experiments by Von Larcher and Egbers (2005) and between m = 3 and m = 4 from the study by Vincze et al. (2015) (note that our initial conditions correspond to the spin up ones).

Figure 4 .

 4 Figure 4. Regime diagram in the T a-Ro T space. The flow regime observed in the top and bottom layers in the barostrat experiment are compared with data available from literature for the same experimental setup, but for the purely thermal configuration (see text and Table1for more details).

  T and T a numbers calculated with (4) and (3) for the data from the paper by[START_REF] Vincze | Benchmarking in a rotating annulus: a comparative experimental and numerical study of baroclinic wave dynamics[END_REF] and the two baroclinically unstable layers at the top and the bottom of the barostrat experiment. λ is the total fluid depth for the classical setup and the thickness of the convective layer for the barostrat experiment. The values of the thermal expansion coefficient α have been corrected for the salinity and the temperature. The arrows at the side of the wavenumbers mark spin-up ↑ and spin-down ↓ initial conditions.

  and Figure5(see colorbar). This suggests that a structure of the wave along the vertical axis is indeed present. Moreover, it can be noticed in Figure5that at the top (z = 94 mm and z = 75 mm) and at the bottom (z = 21 mm) of the tank the inertial Kelvin wave has a radial structure that differs from the middle region (z = 47 mm), namely it shows a higher radial wavenumber. This spatial modulation along the radius might originate from the existence in the top and bottom layers of the baroclinic instability.[START_REF] Gula | Instabilities of two-layer shallow-water flows with vertical shear in the rotating annulus[END_REF] and[START_REF] Flór | Frontal instabilities and waves in a differentially rotating fluid[END_REF] reported instabilities

Figure 5 .

 5 Figure 5. Reconstructed velocity field, using the harmonic analysis, for the principal peaks in the three layers. In the first column on the left are shown the baroclinic waves, m = 3 in the uppermost layer (first and second lines from the top) and m = 4 for the bottom layer (bottom line in figure). The central columns shows the inertial Kelvin wave, for all layers. In the column on the right is shown the field obtained for the main peak at high frequencies ω = 6.8Ω. All the frequencies are normalised by Ω.

Figure 6 .

 6 Figure 6. Variance of the first 10 eigenvalues of the coupled u and v velocities for the top layer (a) z = 94 mm and the bottom layer (b) z = 21 mm.

Figure 7 .

 7 Figure 7. Comparison between the first CEOFs and PCs for the top layer (z = 94 mm) and the bottom layer (z = 21 mm). (a) PC of the first eigenvalue z = 94 mm, (b) PC of the first eigenvalue z = 21 mm. (c) Real part of the CEOF of the first eigenvalue at z = 94 mm and (d) Real part of the CEOF of the first eigenvalue at z = 21 mm

Figure 9

 9 Figure 9(a) shows the spectrum obtained from the temperature sensor at z = 21 mm. We identify two identical shaped baroclinic waves with m = 4 at the frequencies ω b 1 = 0.089Ω and ω b 2 = 0.13Ω, indicated in Figure 9(a). Three other peaks with significant variance are recognisable, one at 2ω b 1ω b 2 = 0.048Ω, at 2ω b 2ω b 1 = 0.17Ω, and at ω b 1 + ω b 2 = 0.22Ω. The first two frequencies indicate the nonlinear interaction of the baroclinic waves ω b 1 and ω b 2 and the mean zonal flow as reported in detail in Buzyna et al. (1989). Figures 9(d), 9(e) show the velocity fields, obtained with the harmonic analysis, for the frequencies 2ω b 1ω b 2 and 2ω b 2ω b 1 . The spatial pattern of a baroclinic wave m = 4 can be seen in both figures, as one would expect from the wave-mean flow interaction.

  Figure8. (a) low frequencies spectra from temperature data at z = 94 mm. The three peaks correspond to a baroclinic wave and its first and second harmonics: ω t 0 = 0.04Ω, ω t 1 = 0.08Ω and ω t 2 = 0.12Ω; the velocity fields corresponding to these most energetic modes reconstructed using the CEOF analysis are plotted in (b),(c) and (d) respectively. The three waves of frequencies ω t 0 , ω t 1 and ω t 2 (at z = 94 mm) form a harmonic triad with azimuthal wavenumbers m t 0 = 3, m t 1 = 6, and m t 2 = 9.

  Figure 9. (a) low frequencies spectra from temperature data at z = 21 mm. The five peaks (whose velocity fields, reconstructed using the harmonic analysis, are shown in (b)-(f)) correspond respectively to two baroclinic waves (ω b 1 = 0.089Ω, ω b 2 = 0.13Ω and m = 4 ), their nonlinear interactions with the mean flow (2ω b 1 -ω b 2 = 0.048Ω, 2ω b 2 -ω b 1 = 0.17Ω and m = 4), and the nonlinear interaction between the two main waves (ω b 1 + ω2 b = 0.22Ω and m = 8).

Figure 11 .

 11 Figure 11. Baroclinic wave jet and position of the line along which the Hovmöller plot is taken (a) at t 0 = 0, i.e., the beginning of the Hovmöller plot in (c) and (b) at t 1 = t 0 + 75 s, i.e., the end of the Hovmöller plot in (c). (c) Hovmöller plot for the divergence at the entrance of the baroclinic wave at z = 94 mm. The dashed line is plotted to show the drift speed of the baroclinic wave, and the red line shows the mean value of the zonal flow in the jet region (u 0 = 3 mm s -1 ).The phase velocity of the small-scale waves, measured from this plot and indicated by the blue line is 3 mm s -1 , the same as u 0 .

Figure 11 (

 11 Figure 11(c) and 15 s < t < 30 s in Figure 12(c), when the red line in Figures 11(a), 11(b) and Figures 12(a), 12(b) is along the entrance and the exit region of the jet respectively.

Figure 12 .

 12 Figure 12. Baroclinic wave jet and position of the line along which the Hovmöller plot is taken (a) at t 0 = 0, i.e., the beginning of the Hovmöller plot in (c), and (b) at t 1 = t 0 + 75 s, i.e., the end of the Hovmöller plot in (c). (c) Hovmöller plot for the divergence at the exit of the baroclinic wave at z = 94 mm. The dashed line is plotted to show the drift speed of the baroclinic wave, and the red line shows the mean value of the zonal flow in the jet region (u 0 = 3 mm s -1 ).The phase velocity of the small-scale waves, measured from this plot and indicated by the blue line is 5 mm s -1 .

  Figure 14(a)), i.e., at the middle of the gap width, to construct a space-time diagram from which we derive a frequency-wavenumber diagram. The Hovmöller plot of the divergence along this arc is shown in Figure 14(b) and the correspondent 2D spectra in Figure 14(c). By Doppler shifting the dispersion relation using u 0 = 1 mm s -1 and the horizontal wavenumber K H = 5 cm -1 we obtain the dispersion curves plotted with red, black and yellow dashed lines in Figure 14(c) for vertical wavelengths λ 75 z = 1, 0.3, 0.2 cm respectively. Obviously, the curve with 3 λ 75 z = 0.2 cm best fits the plotted data. We hence find that the estimated values for the vertical wavelengths at the two different fluid heights (by using (19) for z = 94 mm and by plotting the dispersion curves for z = 75 mm) are very different: λ 94 z

  Figure 14. (a) circle of constant radius r = 8 cm along which the data are taken for the plots (b) and (c). (b) Hovmöller plot of the divergence at z = 75 mm, the red contour lines are Ro > 1. (c) 2D fft of the divergence space-time data plotted in (b). The dashed coloured lines show the dispersion relation for gravity waves (1) considering a horizontal wavenumber K H = (k 2 + l 2 ) ranged from -12 cm -1 to 12 cm -1 and vertical wavenumbers n 1 = 2π/ 1 λz = 6.3 cm -1 , n 2 = 2π/ 2 λz = 18.8 cm -1 , and n 3 = 2π/ 3 λz = 31.4 cm -1 . The dispersion curves are Doppler shifted according to (17) with u 0 = 1 mm s -1 and a horizontal wavenumber K H = 5 cm -1

  stratified flow. We can use the values of the azimuthal velocity at z = 94 mm and z = 75 mm to estimate the vertical shear and the values of N estimated previously to calculate the Richardson number at these two locations. Note that this is only a crude estimation because we have only two measures of the velocity. If we use the value of N ( 0.1 rad/s) that permits the existence of IGWs, the Richardson number at z = 94 mm is Ri 94 = 0.39 which is above but close to 1/4, i.e. just around the limit to rule out shear instabilities. The paradox being of course that if we imagine a higher estimated value of N (say N = 0.2 rad/s for instance), then we cannot interpret the small scales waves observed at z = 94 mm neither by IGWs (the dispersion relation will not be fulfilled) nor by shear instability as the Richardson number will strongly differ from 1/4. Fortunately, at z = 75 mm, the situation is different as the dispersion relation for IGWs is clearly validated and moreover, using the value N 75 = 3.8 rad/s, the Richardson number at this location is Ri 75 = 564. This value, being strongly larger than 1/4, confirms a stably stratified shear flow that rules out any Kelvin-Helmholtz or Hölmböe instability. Moreover, the position of small-scale waves with respect to the baroclinic jet and their occurrence related to an enhanced value of the local Rossby number(Figure 14(b))

  supported by the Spontaneous Imbalance project (HA 2932/8-1) that is part of the research group Multiscale Dynamics of Gravity Waves funded by DFG (FOR1898). This study was initiated by the European High-performance Infrastructures in Turbulence (EuHIT) program and PLG and UH also thank the CNRS LIA 1092-ISTROF for its financial support. M.V. is grateful for the support by the Hungarian National Research, Development and Innovation Office (NKFIH) under Grant Number FK125024. The authors thank Ludwig Stapelfeld, Robin Stöbel, Vilko Ruoff for technical support, Joel Sommeria for the support with the UVmat software and Chantal Staquet for the helpful discussions. Moreover we thank the Spontaneous Imbalance group of MS-GWaves, (Ulrich Achatz, Steffen Hien, Joran Rolland, Lena Schoon, Christoph Zülicke, Illia Horenko, Olga Kaiser) for support and fruitful discussions. Finally we thank the two anonymous referees for their comments that helped to improve the clarity of the paper.

  An estimation of the error coming directly from the PIV software used to get the velocities field from the images (UVmat) is given by two quantities directly calculated by the software, namely the rms difference between the CIV (Correlation Image Velocimetry) and the smoothed velocity field, and secondly the number of vectors excluded because attributed to false vectors (http://servforge.legi.grenoble-inp.fr/projects/softuvmat/wiki/Tutorial/CorrelationImageVelocimetryOptimisation). For our data the values are rms = 0.5 -0.6 pixels for a typical displacement of 5 pixels. This gives an estimation of the PIV error 10% -12%. The percentile of excluded vectors is less than 1.5%, so most vectors are preserved.

  Table 1 for more details).

Table 1 .

 1 -6 K -1 ) Ro T Comparison of the wavenumbers and Ro

	T a	wavenumber m

Baroclinic, Kelvin and inertia-gravity waves

Rodda, Borcia, Le Gal, Vincze, Harlander

The a posteriori quantification of the error on the PIV data we propose here is done by taking two consecutive PIV images (∆T = 0.03 s) and subtracting the velocity fields. Because our flow is rather slow, one can expect that the fields do not differ too much and they are completely uncorrelated from one to the other because to compute the PIV we use two different sets of images. From this we can calculate the relative error for the velocity components as

Repeating this procedure for the entire time serie, we can estimate the mean error associated to the PIV data analysis for the U and V components of the velocity. The calculated mean error is less than 15% across the whole domain, comparable with the error estimation obtained from UVmat software.

The signal-to-noise ratio, used to compare the level of the IGWs signal to the level of the background noise, is calculated as

where σ is the square root of the mean variances of the background noise, i.e., calculated in an area where no IGWs signal is visible in the divergence field and µ is the locally (over a square running filter of dimensions 10 × 10 mm) averaged divergence field.


S Caro

C Dumas

S Garnier

B Furet

Workpiece Placement Optimization in Robotic-based Manufacturing

Keywords: Robotic manipulators, Machining, Quality, Optimization, Workpiece placement

Roboticists are faced with new challenges in robotic-based manufacturing. Up to now manufacturing operations that require both high stiffness and accuracy have been mainly realized with computer numerical control machine tools. This paper aims to show that manufacturing finishing tasks can be performed with robotic cells knowing the process cutting conditions and the robot stiffness throughout its Cartesian workspace. It makes sense that the finishing task of large parts should be cheaper with robots. However, machining robots have not been adapted for such operations yet. As a consequence, this paper introduces a methodology that aims to determine the best placement of the workpiece to be machined knowing the elastostatic model of the robot and the cutting forces exerted on the tool. Therefore, a machining quality criterion is proposed and an optimization problem is formulated and solved. The KUKA KR270-2 robot is used as an illustrative example throughout the paper.

INTRODUCTION

Industrial robots were originally dedicated to "pick and place" operations. They start to be used for machining operations such as trimming, deflashing, degating, sanding and sawing. Machining robots are currently developed to help the operator realize machining operations, which have been mainly performed with Computer Numerical Control (CNC) machines so far. Industrial robots can be used, first to reduce scrap rates and production costs, and secondly to increase the volume and flexibility of production lines. However, industrial robots are not as accurate as CNC machine tools, but the total cost of a machining robot is 30% less compared with an equivalent CNC machine. Therefore, it makes sense that the larger the parts to be machined, the more interesting robot machining.

This paper introduces a methodology that aims to determine the best placement of the workpiece to be machined knowing the cutting forces exerted on the tool and the elastostatic model of the robot. A machining quality criterion is proposed and an optimization problem is formulated and solved. The KUKA KR270-2 robot is used as an illustrative example throughout the paper.

Section 2 presents the workpiece test and the cutting forces applied by the workpiece on the tool along the machined This work was partially supported by the Agence Nationale de la Recherche (ANR), France (Project ANR-2010-SEGI-003-02-COROUSSO).

path. Section 3 describes the robotic cell. Section 4 deals with the formulation of a mono-objective optimization problem to find the optimum placement of the workpiece with regard to a proposed machining quality criterion. Section 5 highlights the optimum and worst workpiece placements within the robotic cell. The best and worst redundancy planning schemes associated with those placements are also determined. Section 6 is about the conclusions of the paper and the future work.

MACHINING TASK

Figure 1 illustrates the workpiece to be machined, which is made up of aluminum alloy. F W of origin O W is the frame attached to the workpiece. The five segments AB, BC, CO W , O W D and DE, of length equal to 200 mm each, have to be milled. The tool path is offset by the tool radius from the five segments to be milled. The tool path is discretized into n points and is shown in dashed line in Fig. 1. Frame F P i is attached to the ith point of the tool path, i = 1, . . . , n. X P i is along the feed direction. Z P i is along the tool axis and points toward the robot. The machining quality is affected by the robot deviation due to the cutting forces applied on the tool [START_REF] Dumas | Workpiece Placement optimization of Six-revolute Industrial Serial Robots for Machining Operations[END_REF]. The cutting conditions are given in Tab. 1 where f z , a p and a e denote the feed rate, the depth of cut and the

F x F y F z Point Number [N]
Fig. 2. Force components applied by the tool on the workpiece expressed in frame F W as a function of the tool path point number width of cut, respectively. The cutting forces are evaluated thanks to the cutting force model described in [START_REF] Tlusty | Dynamics of cutting forces in end milling[END_REF]. As a result, Fig. 2 depicts the force components applied by the tool on the workpiece along axes X W , Y W and Z W of frame F W as a function of the tool path point number.

ROBOTIC CELL AND PARAMETERIZATION

The robotic cell shown in Fig. 3 is composed of a KUKA KR270-2 robot, a FISCHER milling spindle mounted on the robot end-effector and a rotary table.

The geometric parameters of the robot are presented in Dumas et al. (2011a). An identification procedure was developed in Dumas et al. (2011b) in order to determine all joint stiffness values of the KUKA KR270-2 robot. Accordingly, the translational and rotational displacements of the tool can be predicted knowing the wrench exerted on it. It is noteworthy that the flexibility of the links is considered through the identified joint stiffness values. The joint stiffness values of the KUKA KR270-2 robot are given in Tab. 2. k θ l is the lth joint stiffness value, l = 1, . . . , 6.

The FISCHER milling spindle is of type MFW 1709/24. Its maximum rotational speed is equal to 24 000 rpm and its power is equal to 20 kW.

The orientation of the workpiece can be fully represented with the variables (Q 2 , Q 3 , Q 4 ), a subset of the quaternions coordinates. Indeed, the quaternions represent the rotations of the workpiece with a rotation axis u = [u x u y u z ] T and an angle θ. The relation between the quaternions and the axis and angle representation can be found in [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]:

Q 1 = cos(θ/2) (1) Q 2 = u x sin(θ/2) (2) Q 3 = u y sin(θ/2) (3) Q 4 = u z sin(θ/2) (4)
where u 2 x + u 2 y + u 2 z = 1 and 0 ≤ θ ≤ π. Nevertheless, the workpiece can be only rotated about the vertical axis of the robot base frame F 0 due to the geometry of the rotary table. Therefore, . It also depicts the robot base frame F 0 , the frame F 6 attached to the robot endeffector, the frame F 7 attached to the spindle and the workpiece frame F W . From the previous closed loop chain, we can write: 0 T W W T P i = 0 T 6 6 T 7 7 T P i (7) where 0 T W is the homogeneous transformation matrix from frame F 0 to frame F W expressed as:

Q 2 = Q 3 = 0 (5) and -1 ≤ Q 4 ≤ 1 (6) Z 0 Y 0 X 0 F 0 Z 6 Y 6 X 6 F 6 Z 7 Y 7 X 7 F 7 Z W Y W X W F W
0 T W =     2Q 2 1 -1 -2Q 1 Q 4 0 0 x OW Q 1 Q 4 2Q 2 1 -1 0 0 y OW 0 0 1 0 z OW 0 0 0 1     (8)
0 x OW , 0 y OW and 0 z OW being the Cartesian coordinates of point O W expressed in frame F 0 . W T P i is the homogeneous transformation matrix from frame F W to frame F P i attached to the ith point of the tool path as shown in Fig. 1. 0 T 6 is the homogeneous transformation matrix from frame F 0 to frame F 6 . 6 T 7 is the homogeneous transformation matrix from frame F 6 to frame F 7 and depends on the geometry of the spindle and how the latter is mounted on the robot endeffector. Here,

6 T 7 =     0 - √ 2/2 √ 2/2 0 1 0 0 0 0 √ 2/2 √ 2/2 0.684 0 0 0 1     (9)
7 T P i is the homogeneous transformation matrix from frame F P i to frame F 7 . Note that the KUKA KR270-2 robot has six degrees of freedom, whereas the milling operation sets only five degrees of freedom as the rotation of the spindle about the tool axis is not fixed. Therefore, the robot is redundant with respect to the task and the kinematic redundancy is equal to one as explained in [START_REF] Conkur | Clarifying the definition of redundancy as used in robotics[END_REF]; [START_REF] Wenger | Curve-following for redundant manipulators with obstacles : feasibility analysis and solutions[END_REF]. Here, the kinematic redundancy is characterized by the angle β i , which corresponds to the rotation angle of the spindle about the tool axis at the ith point of the tool path as shown in Fig. 5. As a consequence,

7 T P i =    cos(β i ) sin(β i ) 0 0 sin(β i ) -cos(β i ) 0 0 0 0 -1 0 0 0 0 1    (10)

OPTIMUM WORKPIECE PLACEMENT

This section deals with a methodology that aims to determine the best placement of the workpiece to be machined knowing the cutting forces exerted on the tool and the elastostatic model of the robot. The proposed methodology is highlighted through the definition of a machining quality criterion and the formulation of a mono-objective optimization problem.

Machining quality criterion

It makes sense that the tool displacement c i , normal to both the feed direction along axis X P i and the tool axis Z 7 , is mainly responsible for machining errors. 0 x P i and 0 z 7 are the unit vectors of axes X P i and Z 7 expressed in frame F 0 , respectively. c i can be expressed as follows:

c i = 0 δd T i 0 z 7 × 0 x P i ( 11
)
where 0 δd i is the point-displacement of the tool at the ith point of the tool path expressed in frame F 0 and × denotes the cross product. The point-displacement of the tool is evaluated at each point of the tool path thanks to the Cartesian stiffness matrix of the KUKA KR270-2 robot obtained in Dumas et al. (2011a) and from the forces applied by the tool on the workpiece that are depicted in Fig. 2. Interval TI that the geometry of the machined workpiece should respect. Here, the tolerance interval is equal to 0.25 mm. Moreover, the three areas A 1 , A 2 and A 3 are highlighted, namely,

• A 1 amounts to the sum of the geometric errors c i upper bounded by the prescribed tolerance interval along the tool path; • A 2 is the area of the surface, which is lower bounded by the geometric errors c i and upper bounded by the prescribed tolerance interval. • A 3 is the area of the surface, which is lower bounded by the tolerance interval and upper bounded by the geometric errors c i that are higher than the tolerance interval.

The machining quality criterion f M QC , which is the objective function of the optimization problem at hand, is defined as:

f M QC = (w 1 A 1 + w 3 A 3 ) n → min (12)
where w 1 and w 3 are two weighting factors bounded between 0 and 1 and n is the number of discretized points of the tool path. The values of w 1 and w 3 are up to the manufacturer. In the scope of this optimization problem w 1 = 0.2 and w 3 = 0.8.

The smaller f M QC , the better the machining quality.

Decision variables

The decision variable vector of the optimization problem at hand contains the Cartesian coordinates of point O W expressed in frame F 0 , i.e., 0 x OW , 0 y OW and 0 z OW and the orientation parameter Q 4 of the workpiece.

The decision variable vector also contains the redundancy planning scheme. Let β i be the rotation angle of the spindle about the tool axis at the ith point of the discretized tool path as shown in Fig. 5. Accordingly, the redundancy planning scheme is defined by vector β expressed as:

β = [ β 1 β 2 • • • β n ] T (13) with -π ≤ β i ≤ π, i = 1, . . . , n.
The KUKA KR270-2 robot may have eight solutions to the inverse geometric model, namely, eight postures. The robot has to keep the same posture along a given segment in order to avoid signs of wear. However, the robot posture can change from one segment to the next one. µ j denotes the solution number to the inverse geometric model of the robot along the jth segment of the tool path, j = 1, . . . , 5. As a consequence, the posture scheme is defined by vector µ expressed as:

µ = [ µ 1 µ 2 µ 3 µ 4 µ 5 ] T (14) with µ j ∈ [1, .., 8].

Constraints

Four types of constraints arise in the optimization problem at hand:

(1) The joint angles of the KUKA KR270-2 robot should be bounded between its joint limits,

-185 deg ≤ θ 1i ≤ 185 deg (15) 0 ≤ θ 2i ≤ 146 deg (16) -245 deg ≤ θ 3i ≤ 29 deg (17) -350 deg ≤ θ 4i ≤ 350 deg (18) -125 deg ≤ θ 5i ≤ 125 deg (19) -350 deg ≤ θ 6i ≤ 350 deg (20)
θ ji being the jth joint angle of the robot, j = 1, . . . , 6, when the tool is located at the ith point of the tool path, i = 1, . . . , n. Note that the robot was parameterized by using the modified Denavit-Hartenberg convention [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF].

(2) The tool displacement c i expressed with (11) should be smaller than a given value c max , which is defined with regard to the expected machining quality, i.e.,

max (|c 1 |, |c 2 |, • • • , |c n |) ≤ c max ( 21
)
where |.| denotes the absolute value.

(3) The workpiece should be located within a work volume that depends on the rotary table size and location within the robotic cell. Here,

x min ≤ 0 x OW ≤ x max (22)
y min ≤ 0 y OW ≤ y max (23) z min ≤ 0 z OW ≤ z max ( 24
)
with x min = -0.3 m, x max = 0.3 m, y min = -2.24 m, y max = -1.64 m, z min = -0.602 m and z max = -0.214 m. (4) The robot should be far from singularities while the tool follows the path, namely,

1/κ F (J N i ) > 0.2, i = 1, . . . , n (25)
where J N i is the normalized kinematic Jacobian matrix of the KUKA KR270-2 robot defined in Dumas et al. (2011a) and evaluated at the ith point of the tool path. κ F (J N i ) is condition number of matrix J N i based on the Frobenius norm.

Formulation of the optimization problem

From ( 12) to ( 25), the optimization problem to solve in order to find the best workpiece placement can be formulated as follows:

minimize f M QC over x = 0 x OW 0 y OW 0 z OW Q 4 β T µ T T subject to max (|c 1 |, |c 2 |, • • • , |c n |) ≤ c max -185 deg ≤ θ 1i ≤ 185 deg 0 ≤ θ 2i ≤ 146 deg -245 deg ≤ θ 3i ≤ 29 deg -350 deg ≤ θ 4i ≤ 350 deg -125 deg ≤ θ 5i ≤ 125 deg -350 deg ≤ θ 6i ≤ 350 deg 1/κ F (J N i ) > 0.2 x min ≤ 0 x OW ≤ x max y min ≤ 0 y OW ≤ y max z min ≤ 0 z OW ≤ z max -1 ≤ Q 4 ≤ 1 i = 1, . . . , n (26)
Optimization problem (26) aims to find the optimum workpiece placement and the corresponding optimum redundancy planning scheme that minimize the machining quality criterion f M QC defined by ( 12) while respecting the set of constraints.

RESULT ANALYSIS

A hybrid optimization algorithm was used to solve optimization problem (26). As a matter of fact, a genetic algorithm [START_REF] Goldberg | Genetic Algorithms in Search, Optimization & Machine Learning[END_REF] was used to find the initial guess of an interior-point algorithm for large-scale nonlinear programming. The Matlab ga function combined with the Matlab fmincon function wereused. T meaning that the robot posture changes between segment 3 and segment 4 only. The green curve characterizes the optimum redundancy planning scheme, i.e., β opt vector, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement

Figure 8 illustrates the tool displacement c i with respect to the tool path point number and β angle for the optimal workpiece placement. The green curve characterizes the optimum redundancy planning scheme, i.e., β opt vector, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement. The white areas can not be reached by the robot because of its joint limits.

It is noteworthy that 0 x OW , 0 y OW , 0 z OW , Q 4 and µ are the only decision variables considered by the ga and fmincon functions in this optimization problem solving. As a matter of fact, an optimal β vector is searched at each iteration of the genetic algorithm and at each iteration of the interior-point algorithm. This vector is obtained in such a way that it minimizes the objective function f M QC and avoids discontinuities in the robot joint space along each segment. It is noteworthy that the worst workpiece placement is obtained by maximizing the objective function f M QC while respecting the constraints of optimization problem (26). Let us compare the machining quality of the workpiece obtained with the four following milling conditions:

Case 1: Optimum workpiece placement with the best redundancy planning scheme (see Fig. 8) Case 2: Optimum workpiece placement with the worst redundancy planning scheme (see Fig. 8) Case 3: Worst workpiece placement with the best redundancy planning scheme (see Fig. 10) Case 4: Worst workpiece placement with the worst redundancy planning scheme (see Fig. 10)

An video of the milling operation corresponding to Case 1 can be downloaded at 1 . An video of the milling operation corresponding to Case 2 can be downloaded at 2 . An video of the milling operation corresponding to Case 3 can be downloaded at 3 . An video of the milling operation corresponding to Case 4 can be downloaded at 4 5 . Please, be sure that the tilde "∼" symbol before "caro" is kept when the previous four links are copied and pasted in your web browser for downloading the videos.

Table 3 gives the values of the maching quality criterion f M QC defined in (12) for the previous four milling conditions. It is apparent that both the workpiece placement and the kinematic redundancy affect the machining quality of the workpiece. Table 3 shows that it is important to pay attention to the workpiece placement and kinematic redundancy. Indeed, f M QC value for Case 4 is more than 14 times higher than f M QC value for Case 1.

CONCLUSIONS AND FUTURE WORK

This paper introduced a methodology to determine the optimum placement of any workpiece to be machined knowing the cutting forces exerted on the tool and the elastostatic model of the KUKA KR270-2 robot. The theoretical results presented in this paper will be checked experimentally soon. The extension of the proposed methodology to robot machining operations with a higher kinematic redundancy of the robot with regard to the task, such as a machining robot mounted on a rail or a mobile platform, are also part of the future work.

ACKNOWLEDGMENT

  Fig. 1. Workpiece 400 300 200

Fig. 4 .Fig. 5 .

 45 Fig. 4. Closed loop chain and frames F 0 , F 6 , F 7 and F W β i

  Fig. 6. c i along all five linear paths and definition of the different areas considering the tool displacements and areas A 1 , A 2 and A 3 Figure 6 illustrates the tool displacement c i as a function of the tool path point for a given workpiece placement and redundancy scheme. Figure 6 also depicts the Tolerance

Fig. 7 .

 7 Fig. 7. Convergence of the genetic algorithm for the determination of the best workpiece placement The genetic algorithm converged after 51 generations as shown in Fig. 7, each population containing 120 individuals. Then, a local optimum decision variable vector x opt of optimization problem (26) was obtained with the Matlab fmincon function and is expressed as:

Fig. 8 .

 8 Fig. 8. Tool displacement c i [m] as a function of the tool path point number and β angle for the best workpiece placement. The green curve characterizes the optimum redundancy planning scheme, i.e., β opt vector, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement

Fig. 9 .

 9 Fig. 9. Optimum workpiece placement

Fig. 10 .

 10 Fig. 10. Tool displacement c i [m] as a function of the tool path point number and β angle for the worst workpiece placement. The green curve characterizes the optimum redundancy planning scheme, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement Figure 10 illustrates the tool displacement c i with respect to the tool path point number and β angle for the worst workpiece placement. The green curve characterizes the optimum redundancy planning scheme, i.e., β opt vector, whereas the red curve represents the worst redundancy planning scheme for this workpiece placement.

Table 1 .

 1 Cutting Conditions

	Spindle speed Feed speed	fz	ap	ae
	20 000 rpm	4 m/min	0.05 mm 5 mm varies

table

  Fig. 3. Robotic cell located at University of Nantes, France

	Spindle	Robot
	Rotary

Table 2 .

 2 Joint Stiffness Values of KUKA KR270-2 Robot in [MNm/rad]

	k θ 1	k θ 2	k θ 3	k θ 4	k θ 5	k θ 6
	0.237 3.32 2.79 0.486 0.521 0.38

Table 3 .

 3 Machining Quality Criterion for Four Milling Conditions

		Case 1 Case 2 Case 3	Case 4
	f M QC [mm]	0.083	0.35	0.82	1.2

1 http://www.irccyn.ec-nantes.fr/∼caro/IFACMIM2013/video1.avi 2 http://www.irccyn.ec-nantes.fr/∼caro/IFACMIM2013/video2.avi 3 http://www.irccyn.ec-nantes.fr/∼caro/IFACMIM2013/video3.avi 4 http://www.irccyn.ec-nantes.fr/∼caro/IFACMIM2013/video4.avi 5

The authors would like to acknowledge the financial support of the ANR, France (Project ANR-2010-SEGI-003-02-COROUSSO).