Wissam Benmalek

Manuel Collet

Emmanuel Foltête

Mathieu Morvan Ouisse

Emmanuel Foltete

Mathieu

NUMERICAL SIMULATION AND EXPERIMENTAL VALIDATION OF GAP SUPPORTED TUBE SUBJECTED TO FLUID-ELASTIC COUPLING FORCES FOR HYBRID CHARACTERIZATION TESTS

come

Numerical simulation and experimental validation of gap supported tube subjected to fluid-elastic coupling forces for hybrid characterization tests

INTRODUCTION

Steam generators are heat exchangers used to convert water into steam from heat produced in a nuclear reactor core. They are used in pressurized water reactors between the primary and secondary coolant loops.

Figure 1. Steam generator

When the steam generator is operating, water in the secondary loop partially changes to steam. When rising up, this fluid interacts with the U-tubes (see Fig. 1), which are therefore subjected to flow excitation. This excitation can be split into two kind of forces: turbulence forces which are independent of the movement of the tube and the so-called fluid-elastic coupling forces, depending on acceleration, velocity, displacement and fluid reduced velocity [START_REF] Tanaka | Fluid elastic vibration of tube array in cross flow[END_REF][2]& [START_REF] Hassan | Time Domain Models for Damping-Controlled Fluidelastic Instability Forces in Tubes With Loose Supports[END_REF] . The total flow excitation can be finally expressed as: 𝐹 𝑇 + 𝐹 𝑓 (𝑦, 𝑦, 𝑦, 𝑉𝑟) = 𝐹 𝑇 -𝑀 𝑓 . 𝑦̈-𝐶 𝑓 (𝑉𝑟). 𝑦̇-𝐾 𝑓 (𝑉𝑟). 𝑦 𝑉 = 𝑉/(𝑓 * 𝐷)

Where 𝑦, 𝑦̇ and 𝑦 are acceleration, velocity and displacement vectors. Under some specific conditions of fluid reduced velocity, 𝐶 𝑓 is negative enough to make the structure instable. This phenomenon is called fluid-elastic instability and can damage the structure. The tubes are however supported by plates (see Fig. 1) which guide them and limit their vibration amplitude. In fact, the impacts between the tubes and the plate tend to stabilize the tubes. Thus, we can finally represent the whole problem as below:

ℳ • 𝑦̈+ 𝒞 • 𝑦̇+ 𝒦 • 𝑦 = 𝐹 + 𝐹 𝑓 (𝑦, 𝑦, 𝑦, 𝑉 ) + 𝐹 (2)
Premultiplying the Eq.( 2) by modal base Φ we get a set of independent equations of motion in the modal coordinates

ℳ • 𝑞̈+ 𝒞 • 𝑞̇+ 𝒦 • 𝑞 = 𝑓 + 𝑓 𝑓 (𝑞̈, 𝑞̇, 𝑞, 𝑉 ) + 𝑓 (3)
Where 𝑞̈, 𝑞̇, 𝑞 are generalized acceleration, velocity and displacement vectors. Because these mechanisms are complex and difficult to realize in an experimental set up, the main aim of our study is to develop a hybrid control loop to simulate this coupling effect in the frame of an experimental characterization test bench.

STRUCTURE & MODAL UPDATING Figure 2. Gap supported tube

The studied structure Fig. 2 is composed of a tube attached to a slender plate clamped in rigid block. At the middle height of the tube, two gap stops located at 0.5mm create punctual impacts depending on the vibration amplitude of the tube.

A finite element shell model was developed to generate the mass, damping and stiffness matrices. This model was updated in order to match the numerical behavior with experimental one. Two criterions were used to compare the numerical and experimental model: Modal Assurance Criterion (MAC) and frequency error criterion. Fig. 3 & Tab.1 summarizes the results obtained for the 6 first modes.

EXPERIMENTAL IDENTIFICATION

Impact stiffness:

The only source of nonlinearity in our problem comes from the impact forces located at x c (middle height of tube), which are computed in an explicit manner as the Eq.( 4) shows.

Where K c is the impact stiffness. The value of the parameter K c is identified through experimental measurements. Several impact tests were performed in which the impact forces where measured by using force sensor. Knowing the mass sensor "M c " and measuring the duration of the impact "Tc", we can deduce an estimation of the impact stiffness using the following approximation:

𝑇𝑐 = 𝜋 𝜔 = 𝜋√𝑀 /𝐾 (5)

Figure 4. Impact stiffness identification

Figure 4 presents a histogram of the obtained results for 25 tests. K c was estimated at 10 6 N/m on average (for T c = 0.3ms on average), which matches with the results found by P.Piteau [START_REF] Piteau | Experiments and computations of a loosely supported tube in a rigid bundle subjected to single-phase flow[END_REF] & T.Thenint [START_REF] Thénint | Etude d"un système non-linéaire à chocs sous excitation large bande. Apllication à un tube de générateur de vapeur[END_REF].

Fluid-elastic forces parameters:

In order to model fluid-elastic forces, a significant research effort has been conducted over the last four decades. These researches leaded to several theoretical models. We can quote M.Paidoussis [START_REF] Païdoussis | Fluid-Structure Interactions: Cross-Flow-Induced Instabilities[END_REF], S.Price [START_REF] Price | A Review of Theoretical Models for Fluidelastic Instability of Cylinder Arrays in Cross-Flow[END_REF] or S.Granger [START_REF] Granger | Motiondependent excitation mechanisms in a square in-line tube bundle subject to water cross-flow: an experimental modal analysis[END_REF] models etc. In our study we have chosen the CEA one [START_REF] Piteau | Vibro-Impact Experiments and Computations of a Gap-Supported Tube Subjected to Single-Phase Fluid-Elastic Coupling Forces[END_REF] which is semianalytical model. In this approach, the parameters M f , C f and K f are identified experimentally.

Referring to works carried out by CEA [START_REF] Piteau | Experiments and computations of a loosely supported tube in a rigid bundle subjected to single-phase flow[END_REF], the first mode is predominant in the tube response and the effects of higher modes are almost negligible, thus the fluid elastic force is projected only on the first mode and the other modes are not influenced.

𝑓 𝑓. (𝑞, ̈𝑞, 𝑞) = -(𝑚 𝑓 𝑞̈+ 𝑐 𝑓 (𝑉 )𝑞̇+ 𝑘 𝑓 (𝑉 )𝑞) [START_REF] Païdoussis | Fluid-Structure Interactions: Cross-Flow-Induced Instabilities[END_REF] Therefore, the fluid-elastic effect can be finally modeled as fluid added mass, fluid added damping and fluid added stiffness. These last two coefficients are assumed depending on reduced fluid velocity V r [START_REF] Piteau | Experiments and computations of a loosely supported tube in a rigid bundle subjected to single-phase flow[END_REF].

[ 𝑚 + 𝑚 𝑓 0 0 0 𝑚 0 0 0 ⋱ ] 𝑞̈+ [ 𝑐 + 𝑐 𝑓 (𝑉 ) 0 0 0 𝑐 0 0 0 ⋱ ] 𝑞+ [ 𝑘 + 𝑘 𝑓 (𝑉 ) 0 0 0 𝑘 0 0 0 ⋱ ] 𝑞 = 𝑓 + 𝑓 (7)
At this stage, in order to completely identify fluid-elastic forces we should determine the three parameters m f , c f (V r ) and k f (V r ). The modal characteristics of the structure are assumed to be known in the air and the fluid-added mass m f doesn"t depend on fluid velocity. This parameter is then identified through the measurement of the frequency of the first mode in stagnant water and by comparing it to the equivalent measurement in air m f . Concerning fluid added stiffness k f and damping c f , we infer these coefficients by measuring for each fluid velocity first mode frequency and damping coefficient (see Eq.( 6) & Eq.( 7)).

𝜔 = √ 𝑘 + 𝑘 𝑓 𝑚 +𝑚 𝑓 (8)
𝜁 = 𝑐 + 𝑐 𝑓 2𝜔(𝑚 + 𝑚 𝑓 ) (9)
In the Fig. 5 & Fig. 6 we recall the results of measurements given by [START_REF] Piteau | Experiments and computations of a loosely supported tube in a rigid bundle subjected to single-phase flow[END_REF] for each value of fluid velocity. As we can notice, the damping become negative from 2.1 m/s and the fluid elastic instability phenomenon appears. Nevertheless in non-linear case (with impacts), the structure remains stable because the gap stop limit the vibration amplitude of tube. In order to model this phenomenon, we have to introduce another dimensionless concepts, reduced velocity, coupling damping coefficient C d and coupling stiffness coefficient C k . These dimensionless coefficients will help to explain the stabilizing effect of the impacts.

In his works, S.Caillaud [START_REF] Caillaud | Excitation forcée et contrôle actif pour la mesure des forces fluide-élastiques[END_REF] introduced the dimensionless coupling coefficient as follow:

𝐶 𝑑 = 8𝜋(𝑚 + 𝑚 𝑓 ) 𝜁 • 𝜔 𝜔 -𝜁 𝜌 𝑓 • 𝐷²𝐿 𝑉 (
Where L e is the modal equivalent length calculated from the first mode shape, ζ o and ω o are the modal damping coefficient and circular frequency in stagnant water.

When including impacts in the problem, the expression of the reduced fluid velocity in Eq.( 1) becomes a little more complicated, since it depends on vibration apparent frequency of tube and not anymore on the first modal frequency [START_REF] Price | A Review of Theoretical Models for Fluidelastic Instability of Cylinder Arrays in Cross-Flow[END_REF].

There are different way to estimate this vibration tube frequency, we can quote for example zero-crossing method [START_REF] Boashash | Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals[END_REF] & [START_REF] Au-Yang | THE CROSSING FREQUENCY AS A MEASURE OF HEAT EXCHANGER SUPPORT-PLATE EFFECTIVENESS[END_REF] and Rice frequency method f R [START_REF] Rice | Mathematical analysis of random noise[END_REF]. This last method was chosen for our problem.

f R is proportional to velocity RMS over displacement RMS of the tube free end within sliding size windows τ. f R (t, τ) = σ ẏ( t, τ) 2πσ y (t, τ) (

The dimensionless coupling coefficients C d and C k can be inferred from ζ and ω already measured. In addition, as shown in Eq.( 10) & Eq.( 11), they depend on reduced velocity. Figure 7 presents the variation of dimensionless coupling coefficient depending on reduced velocity. It can clearly be seen that C d become completely negative from about 3.8. After determining experimentally all parameters, the next step is to create a numerical solver including all concepts we introduced to solve the nonlinear problem. The Newmark time integration solver has been implemented to get the tube response. Numerous studies have been conducted in linear as well as in nonlinear problems [START_REF] Géradin | Théorie des vibrations: application à la dynamique des structures[END_REF].

Instable

In our approach, a mechanical conservation energy criterion has been used in order to verify trustworthiness of the numerical results. 𝐸𝑚 is constituted by sum of three terms: kinetic energy, potential energy and shocks potential energy. At each computation step, theses energy are evaluated as follows and the sum is divided by the initial mechanical energy in order to get the conservation in percentage.

𝐸𝑐 ( ) = 1 2 . (𝑦( ) ) 𝑇 . ℳ. (𝑦( ) ) (15)
𝐸𝑝 ( ) = 1 2 . (𝑦 ( ) ) 𝑇 . 𝒦. (𝑦 ( ) ) ( 16)

𝐸𝑝𝑐 ( ) = 1 2. 𝐾𝑐 . (𝐹𝑐 ( ) ) 𝑇 . (𝐹𝑐 ( ) ) (17)
First we applied this criterion on Newmark method alone without any correction algorithm taking in consideration the impacts, but it does not insure the conservation energy especially for long time simulation. To avoid this behavior, a Newton-Raphson algorithm [START_REF] Paultre | Dynamique des structures : Application aux ouvrages de génie civil[END_REF] slightly modified has been implemented to in order to fit our tolerance requirement concerning the impact time. Figure 10 presents the energy conservation criterion after implementing Newton-Raphson algorithm for a computation step of 2.10 -5 s, we reached around 100% of energy conservation with less than 1% of error Fig. 11.

IMPACTS STABILIZING EFFECT

At this stage, the parameters describing the fluid elastic force are depending on reduced velocity. The later includes the effect of the impacts through apparent frequency f R .

The impact stabilizing effect can be explained by the conjunction of two phenomena: The first one is the dissipation due to the higher orders modes of the structure. In fact, when impact occurs, many modes not excited by the turbulent flow are excited, and will then dissipate a part of energy increasing the stability of the tube. The second phenomenon is the increase of the stiffness of the structure by being in contact with a stiff obstacle. During contact, since the apparent stiffness is increased, apparent frequency of the structure should also increase, hence decreasing the reduced velocity (Eq.( 12)), which makes the fluid elastic load dissipating. The numerical results in Fig. 9 illustrate well this phenomenon.

NUMERICAL VS EXPERIMENTAL RESULTS

In this part we compare our first numerical results to the experimental one [START_REF] Piteau | Vibro-Impact Experiments and Computations of a Gap-Supported Tube Subjected to Single-Phase Fluid-Elastic Coupling Forces[END_REF] for gap equal to 0.5mm and flow velocity of 2.1m/s in order to validate our numerical approach. On the left side, a detailed plot of the free end displacement for numerical simulation and experimental measurement are shown. On the right, there associated histograms are presented.

Figure 12. Computed response & Experimental response respectively

The direct comparison of the detailed plots is not appropriate and results should be compared in a statistical sense, since turbulence forces are modeled by random excitation with a bandwidth of 40Hz. But through this plot we can highlight the displacement which is bounded at about double distance of the gap. This can be apprehended from the first mode shape. Concerning the histogram, the distribution is centered on zero, namely the initial equilibrium position. Also notice the good agreement of the estimated time between the two results.

ACTIVE CONTROL

In order to design the active part of the bench which simulates the fluid-elastic force effect, the characterization of each control loop component (shaker, accelerometer etc) must be performed. This stage is crucial because it allows complete mastering of the system. In this paper we do not give details about this modeling stage but we can refer to Dorf and Bishop works [START_REF] Dorf | Modern Control Systems[END_REF] for more details.

Figure 13. Control loop

After this phase of modeling, our efforts were conducted on plant designing [START_REF] Preumont | Vibration Control of Active Structures: An Introduction[END_REF]. The aim is to impose a root for first pole (first mode) in such a way that the whole system fit as much as possible the evolution of both first modal frequency and modal damping coefficient previously measured for different value of the control gain G. In this way, by controlling the gain the corresponding reduced velocity configuration can be reproduced. Figure 14 gives an example of a root locus of the system for a plant that we are still working on. As it can be noticed, the system becomes instable around the modal frequency of the tube in stagnant water for G = 0.25 and fluid elastic instability occurs. We are currently working on the experimental setup to validate these results before adding gap stops and studying the stabilizing effect of the impact.

Figure 3 .

 3 Figure 3. Modal updating (MAC)

Figure 5 .

 5 Figure 5. Identified first modal frequency depending on fluid velocity (m/s)

{ 4 )Figure 6 .

 46 Figure 6. Identified first modal damping coefficient depending on fluid velocity (m/s)

10 )

 10

Figure 7 .Figure 8 .

 78 Figure 7. Identified dimensionless fluid-elastic damping coupling coefficient depending on dimensionless reduced fluid velocity

Figure 10 .Figure 11 .

 1011 Figure 10. Newmark with Newton Raphson method conservation energy

Figure 9 .

 9 Figure 9. (a) Tube response & estimated impact force (b) Apparent frequency (c) Reduced velocity

Figure 14 .

 14 Figure 14. Root locus

  It can be observed that there is a good agreement between the two models. The lowest values of MAC is 85% and the mean frequency error is lower than 4%.

		Secondary
		steam exit
		Secondary
		water entrance
	Anti-vibration
	bars
	Tube support plate
	Tube of primary
	loop
	Primary water exit	Primary water entrance
		Modes	Numerical modal	Experimental modal	Error
			frequency	frequency
		1	23.1	23	0%
		2	103	100	3%
		3	305	318	4%
		4	426	407	4.6%
		5	774	798	3%
		6	1300	1404	8%
				Mean	3.7%


Christian Lexcellent

M Collet

M Ouisse

E Foltete

Investigation on time response to impact test of a 2D shape memory structure

Keywords:

A numerical implementation of the RL model for shape memory alloys (SMA) is presented in the paper, adapted to high strain rate loading. The objective is to predict the time response of a 2D SMA structure subjected to an impulse force (duration 50 ms) and following free vibration with a decreasing amplitude. The choice of the material mechanical properties has been done in order to have phases transformations during the oscillations. The damping due to these phase changes is clearly identified when the results are compared with a linear model without induced martensite. The thermomechanical constitutive relation of the SMA has been fullfiled to be able to take into account time reaction when the strain rate is very high. The full model has been implemented in a finite element code and tested on a 2D sample.

The SMA description which has been implemented is based on the RL model [START_REF] Raniecki | Thermodynamic models of pseudoelastic behavior of shape memory alloys[END_REF], built from thermodynamics of irreversible process. A Helmholtz free energy is defined for the two-phase reference volume element containing the volume fraction (1ξ) of austenite phase and ξ of martensite phase: are the internal energy and entropy of austenite phase, T is the temperature, π f 0 (T ) is the driving force of the stress free state, ε is the strain tensor, L is the Hooke tensor, φ it is a coefficient of internal interaction, and ε tr is the the phase transformation strain tensor:

φ(ε, T, ξ) = u (A) 0 -T s (A) 0 -ξπ f 0 (T ) + 1 2 ε -ε tr : L : ε -ε tr + ξ(1 -ξ)φ it (1) u
ε tr = 3 2 dev(σ) σ vm γξ ( 2
)
σ vm is the Von Mises equivalent stress and dev(σ) the deviator stress tensor. The driving force of phase transformation is:

π f (σ, ξ, T ) = γσ vm ρ + π f 0 (T ) -(1 -2ξ)φ it (3)
The Clausius Duhem inequality leads to the kinetic laws of the phase transformations:

ξ = f (σ vm , π f , ξ) ( 4 )
The function f is zero except for austenite to martensite transformation (0 < ξ < 1 and σvm > 0 and π f (T ) > 0):

f (σ vm , π f , ξ) = γ σvm ρ A 1 1-ξ -2φ it (5)
and for martensite to austenite transformation (0 < ξ < 1 and σvm ≤ 0 and π f (T ) ≤ 0):

f (σ vm , π f , ξ) = γ σvm ρ A 2 ξ -2φ it (6)
The model has been implemented in COMSOL Multiphysics in a finite elements context, in the same way as the one described in details in ref. [START_REF] Thiebaud | Implementation of a model taking into account the asymmetry between tension and compression, the temperature effects in a finite element code for shape memory alloys structures calculations[END_REF] in the context of classical structural dynamics analysis. An implicit timestepping scheme is used for the time-dependent solver algorithm, based on DASPK solver [START_REF] Brown | Using Krylov methods in the solution of large-scale differential-algebraic systems[END_REF] (variable-order and variable-stepsize backward differentiation). The nonlinear system of equations at each time step is solved using a Newton-Raphson scheme. The tangent unsymmetric sparse linear system is finally solved using direct UMFPACK solver [START_REF] Davis | Algorithm 832: UMFPACK V4.3, an unsymmetric-pattern multifrontal method[END_REF]. In the context of transient dynamics, such as responses to impacts, some numerical tests have shown that the convergence of the solvers was quite bad, inducing very large computation times or even not converged results. One of the main reasons for this is undoubtely the strong nonlinearity of the kinetic law 4 whose expression is conducted by the value of the thermodynamic force π f and of the equivalent Von Mises Stress σ vm . The way which has been chosen in this case for smoothing the numerical difficulties induced by the kinetic law, is to add an inertial term in it. In a physical point of view, the austenite to martensite transformation is not instantaneous, and it has been observed that the speed of transformation is linked to the speed of sound in the alloy. This term could have a significant impact when structural solicitations induce speeds which are comparable to the sound velocity, while being negligible for quasi-static or "low-frequency range" applications. The kinetic law 4 has then been modified to :

t m (ξ) ξ + ξ = f (σ vm , π f , ξ) ( 7 )
in which t m (ξ) is a relaxation function which can be chosen depending on the speed of sound in the alloy. This term acts as a time delay in the response function and regularizes the phase transformation behavior. Its own dependance on the martensite volume fraction ξ insures Lyapunov stability of this equation and also prevents overbounded values of ξ superior to 1 or inferior to 0. In the numerical tests, t m (ξ) is a function which is close to 10 -3 when ξ belongs to [0, 1] and zero otherwise. As a first observation from the numerical tests, one can observe that the phase transformation has a significant effect on the global damping. On figure 1b, two results are shown : the first one is related to the SMA calculation corresponding to the above description, and the second one corresponding to a linear equivalent material (same resolution, with an imposed zero martensite rate). The curve shows the vertical displacement of right observing point, and one can clearly observe that the SMA decrease is faster than the linear one. A very simple 1-dof equivalent model can be used to evaluate the added damping, using a log-decrement approach. The time response of the linear model has a frequency of 360 Hz associated to a damping ratio of 0.11, while the time response of the SMA model has a frequency of 243 Hz and a damping ratio of 0.24. In both cases, the first mode is the main component of the response. This 1-dof equivalent linear model is clearly not useful for system modeling since there is a strong nonlinearity of the response versus amplitude of excitation force, but it can be used to obtain qualitative information about added damping. The other physical parameters (martensite rate, thermodynamic force and equivalent Von Mises Stress) have some coherent physical behaviors, that will be confirmed using comparisons with experiments. Moreover, the added term in the kinetic laws has a great positive effect on convergence, without affecting the results given by the original model.

Figure 1 .

 1 Figure 1. a) 2D-structure; b) Vertical displacement of observing point


Morvan Ouisse
email: morvan.ouisse@univ-fcomte.fr

Emmanuel Foltete
email: emmanuel.foltete@univ-fcomte.fr

Une méthode de modification des modes complexes couplés en identification modale pour la vérification de la propreté vibroacoustique

Dans le cadre de l'analyse modale expérimentale en dynamique des structures, l'identification passe par la détermination des modes complexes. Si l'on s'intéresse à la représentation du comportement de la structure sur une plage de fréquence réduite, les déformées modales obtenues doivent vérifier la condition de "propreté" [1], qui traduit la possibilité de représenter ce comportement par un système réduit équivalent, du second ordre. L'extension de cette propriété au cas du couplage vibroacoustique se fait facilement. Les déformées modales complexes permettant de représenter correctement la dynamique couplée du système sont obtenues à partir des déformées identifiées, qui sont modifiées si elles ne vérifient pas les conditions de propreté. Le problème à résoudre est alors un problème de minimisation sous contrainte, qui a une solution analytique dans le cas de la structure dans le vide [1]. Dans le cas du couplage, on obtient un système matriciel non linéaire qui est explicité dans les travaux présentés ici. Sa résolution globale n'a pas encore été effectuée, mais on peut néanmoins en traiter la partie "diagonale" par blocs, par des techniques classiques faisant notamment intervenir la résolution d'équations de Riccati. Cette procédure permet d'améliorer sensiblement la propreté du système couplé, comme le montrent les simulations numériques et les résultats expérimentaux présentés dans ce papier. Une seconde procédure, appelée "sur-nettoyage", est présentée, et permet de vérifier toutes les conditions requises, mais en imposant en plus des conditions théoriquement non nécessaires. Ces deux procédures sont appliquées et comparées sur un cas-test numérique, et sur deux jeux de données de mesures vibro-acoustiques.

I. INTRODUCTION

Les comportements de systèmes vibro-acoustiques peuvent être à l'origine de situations où de forts couplages sont présents, qu'ils soient désirés comme dans le cas des instruments de musique, ou non pour un couplage entre le domaine fluide d'une cavité automobile et les vibrations de la caisse par exemple. Nous nous intéressons ici au cas de l'analyse modale vibro-acoustique de systèmes couplés, constitués d'une structure vibrante fermant un domaine acoustique fluide (vibroacoustique interne).

L'analyse modale vibro-acoustique [START_REF] Wyckaert | Vibro-acoustical modal analysis : Reciprocity, model symmetry and model validity[END_REF] permet d'identifier expérimentalement les modes couplés du système, de la même façon que l'analyse modale classique en dynamique des structures permet l'identification des modes d'une structure dont les vibrations sont supposées découplées du comportement du fluide. La synthèse modale consiste notamment à rechercher une représentation du comportement du système sur une bande de fréquence par un modèle expérimental construit sur une base réduite constituée d'un certain nombre des modes identifiés. Cette identification conduit généralement à la détermination des modes complexes de la structure, ces grandeurs étant souvent entachées d'erreurs dont les origines peuvent être multiples, du fait de l'approche expérimentale utilisée. Ces erreurs ont des ordres de grandeur variables suivant les techniques utilisées et l'expertise de l'opérateur, mais même de petites erreurs peuvent conduire à des reconstructions de système différentes.

Il existe une condition mathématique permettant d'assurer que le système étudié peut être représenté par un système discret, appelée condition de propreté des vecteurs propres complexes. Cette condition est très bien explicitée en référence [START_REF] Balmes | New results on the identification of normal modes from experimental complex ones[END_REF]. Lorsque les vecteurs ne vérifient pas la condition de propreté, l'auteur présente une méthodologie de modification de ces vecteurs afin que la condition soit vérifiée. Nous proposons ici d'étendre cette approche au cas de la vibroacoustique.

II. DESCRIPTION DU PROBLÈME ET DÉCOMPOSITION MODALE

A. Équations de mouvement

La discrétisation d'un problème de vibro-acoustique interne, basée sur une description en déplacement de la structure et en pression de la cavité, conduit au système matriciel suivant [START_REF] Morand | Fluid Structure Interactions[END_REF] :

M s 0 C T M f M ẍ p Ẍ + B s 0 0 B f B ẋ ṗ Ẋ + K s -C 0 K f K x p X = F s q F (1)
où :

x est le vecteur déplacement généralisé de la structure p est le vecteur des pressions acoustiques -M s est la matrice de masse de la structure, obtenue par discrétisation de la forme bilinéaire Ωs ρ s uδudΩ, ρ s étant la masse volumique de la structure. Cette matrice est symétrique, définie, positive.

-M f est la matrice de "masse" fluide, obtenue par discrétisation de la forme bilinéaire 1 ρ f c 2 Ωf p∇δpdΩ, c étant la célérité du son dans le milieu fluide, et ρ f la masse volumique du fluide. Cette matrice est symétrique, définie, positive. Ses composantes n'ont pas la dimension d'une masse physique, elle est appelée matrice de masse par analogie avec les vibrations de structure.

-K s est la matrice de raideur de la structure, obtenue par discrétisation de la forme bilinéaire Ωs σ ij (u)ǫ ij (δu)dΩ.

Cette matrice est symétrique, semi-définie, positive. -K f est la matrice de "raideur" fluide, obtenue par discrétisation de la forme bilinéaire q est le vecteur représentant les sources acoustiques (débits volumiques) placées dans la cavité. Dans la suite, on notera n s le nombre de degrés de liberté de la structure, n f le nombre de degrés de liberté du fluide, et n = n s + n f .

B. Réponse du système en régime harmonique par décomposition modale

Afin de résoudre le problème couplé, il est possible de procéder par décomposition modale. La nature non symétrique du problème fait qu'il est nécessaire de calculer les modes propres gauches et droits du problème. Cette résolution passe par la représentation d'état du système :

U Ż -AZ = F (2)
où :

U = B M M 0 , A = -K 0 0 M Z = X Ẋ , F = F 0 (3)
La résolution de ce problème aux valeurs propres est souvent peu efficace d'un point de vue numérique car il met en jeu des matrices non symétriques pour lesquelles les algorithmes d'extraction de valeurs propres ne sont généralement pas optimisés. Néanmoins, cette extraction est possible, et conduit à l'obtention des valeurs propres complexes s ν du système (à chaque valeur propre s ν du système est associée s ν qui est également valeur propre). La matrice des valeurs propres est notée :

S = diag(s ν ) 0 0 diag(s ν ) (4)
A ces valeurs propres sont associés : -Les vecteurs propres de droite dy ν qui sont tels que U T gyS = A T gy ou gy T U S = gy T A (6)

(s ν U -A)dy ν = 0,
Les relations d'orthogonalité s'écrivent alors, en notant ξ la matrice de "masse modale", diagonale :

gy T U dy = ξ ou gy T Ady = ξS (7)
La réponse harmonique du système projetée sur la base modale s'écrit alors, pour une pulsation d'excitation ω :

Z(t) = dy (ξ(jωE 2n -S)) -1 gy T F e jωt (8)
où E 2n est la matrice d'identité de taille 2n. Cette relation peut s'écrire sur les vecteurs de taille n dans le domaine fréquentiel sous la forme :

X(ω) = dY ΞgY T F (9)
avec :

Ξ = diag( 1 ξν (jω-sν ) ) 0 0 diag( 1 ξ ν (jω-sν ) ) (10)

C. Particularités du couplage vibro-acoustique

Le couplage vibroacoustique caractérisé par l'équation (1) n'est pas auto-adjoint, mais les termes extra-diagonaux des matrices de masse et de raideur sont reliés entre eux. On montre facilement que les vecteurs propres de gauche sont reliés aux vecteurs propres de droite par une relation particulière [START_REF] Wyckaert | Vibro-acoustical modal analysis : Reciprocity, model symmetry and model validity[END_REF] :

Si dY ν = X ν P ν alors gY ν = X ν -Pν s 2 ν (11)
Ainsi, la résolution du problème aux valeurs propres de gauche n'est pas requise pour déterminer la base modale complète.

Cette relation peut aussi être écrite de la façon suivante :

Si dY = X P alors gY = X P ∆ (12) où ∆ = -S -2 .
Dans ces conditions la masse modale peut s'écrire plus simplement :

ξ ν = 1 s ν X T ν (M s s 2 ν -K s )X ν - 1 ρ f s 2 ν P T ν (M s s 2 ν -K s )P ν (13)

III. CONDITION DE PROPRETÉ DES VECTEURS PROPRES COMPLEXES

A. Cas des vibrations de structure

La notion de propreté des vecteurs propres complexes est très bien détaillée en référence [START_REF] Balmes | New results on the identification of normal modes from experimental complex ones[END_REF]. A partir d'un certain nombre de modes complexes identifiés, cette relation traduit l'aptitude du système à être complètement représenté (sur la bande de fréquence d'interêt) par un système réduit équivalent limité construit à partir des seuls modes identifiés. Si un ensemble de vecteurs propres ne vérifie pas cette relation, c'est que le système ne peut être représenté de cette manière, ou que l'identification est entachée d'erreurs ayant conduit à une mauvaise identification des vecteurs complexes. Dans ce dernier cas, Balmès propose une méthodologie de modification des vecteurs afin qu'ils vérifient la relation de propreté, tout en étant aussi proches que possible des vecteurs initiaux. Dans le cas de la structure pure, la relation de propreté s'écrit, si les modes sont normalisés par ξ ν = 1 :

dY dY T = 0 (14)
Dans le cas où les vecteurs propres ne vérifient pas cette relation, on cherche dY , approximation de dY , tel que :

1) dY dY T = 0 2) dY -dY minimale
Il a été démontré qu'une solution de ce problème était :

dY = E n -δδ -1 dY -δdY (15)
où δ est une matrice de multiplicateurs de Lagrange, qui est obtenue par résolution de l'équation de Riccati suivante :

dY dY T -δdY dY T -dY dY T δ + δdY dY T δ = 0 (16)

B. Cas d'un problème non auto-adjoint

L'obtention de la condition de propreté dans le cas d'un problème non auto-adjoint est immédiate, à partir de l'expression de la réponse forcée (8), ou de l'inversion des relations d'orthogonalité (7), dans le cas où la normalisation est telle que ξ = E 2n :

U -1 = dygy T (17)
Soit :

B M M 0 -1 = 0 M -1 M -1 -M -1 BM -1 = dY gY T dY SgY T dY SgY T dY S 2 gY T (18) Et : A -1 = dySgy T (19)
Soit :

-K 0 0 M -1 = -K -1 0 0 M -1 = dY S -1 gY T dY gY T dY gY T dY SgY T (20)
Il est ainsi clair que la relation de propreté dans le cas d'un système non auto-adjoint s'écrit :

dY gY T = 0 (21)
Si cette condition est vérifiée, les matrices caractéristiques du système peuvent être directement déterminées par les relations suivantes :

M = dY SgY T -1 (22) K = -dY S -1 gY T -1 (23) B = -M dY S 2 gY T M (24)
Dans le cas particulier du couplage vibro-acoustique, les vecteurs propres gauches sont déduits des vecteurs propres droits, par la relation (12). La condition de propreté (21) peut ainsi s'écrire en fonction des seuls modes de droite : Dans cette partie, on suppose connue la matrice modale de droite dY contenant les vecteurs propres de droite identifiés mais ne respectant pas la condition de propreté :

XX T X∆P T P X T P ∆P T = 0 (
dY = X P (26)
Le vecteur propre de gauche est alors obtenu par la relation (12). On cherche alors le vecteur dY , approximation de dY , vérifiant au mieux la condition de propreté vibroacoustique (25), et minimisant la norme dY -dY . On note W et Q les composantes associées au déplacement de structure et à la pression acoustique du vecteur recherché :

dY = W Q (27)

A. Mise en équation du problème de minimisation complet

On cherche à résoudre le problème de minimisation sous contrainte suivant :

Trouver W et Q minimisant W -X et Q -P sous la contrainte W W T W ∆Q T QW T Q∆Q T = 0 (28) où X et P sont deux matrices complexes n s × 2n et n f × 2n, connues, et ∆ une matrice diagonale complexe 2n × 2n connue.
La résolution de ce problème passe par la définition du Hamiltonien H, faisant intervenir 4 matrices de multiplicateurs de Lagrange réelles δ i (i=1 à 4), que l'on peut écrire sous forme symbolique :

H = W -X + α Q -P +δ 1 ⊗ W W T + δ 2 ⊗ QW T +δ 3 ⊗ W ∆Q T + δ 4 ⊗ Q∆Q T (29)
L'écriture des conditions de minimisation du Hamiltonien, conduisent après développement au système d'équations matricielles suivant :

               0 = W Q - X P + 1 2 δ 1 + δ T 1 δ T 2 δ 2 0 W Q + 1 2 0 δ 3 δ T 3 δ 4 + δ T 4 W ∆ Q∆ 0 = W W T W ∆Q T QW T Q∆Q T (
30) Notons que l'on retrouve bien l'expression fournie par Balmès dans le cas de la structure seule :

W -X + δW = 0 (31)
où δ est réelle et symétrique. De cette équation est déduite l'expression de X en fonction de δ :

W = (E n -δδ) -1 (X + δX) (32)
δ est alors déterminé en écrivant la condition de propreté qui fournit l'équation de Riccati à résoudre (16).

Le système (30) est un système d'équations matricielles, non linéaires, qui nécessite pour sa résolution l'utilisation d'algorithmes numériques lourds à mettre en oeuvre. Nous n'avons pas tenté de résoudre numériquement ce problème, mais avons utilisé des stratégies permettant d'obtenir des solutions non optimales au sens du problème (28), mais néanmoins plus proches de la solution optimale que ne le sont les vecteurs complexes de départ. Ces stratégies sont présentées dans la suite.

B. Nettoyage de la partie diagonale de la condition de propreté

Une première approche permettant le nettoyage partiel de la solution est de procéder uniquement au nettoyage de la partie diagonale de la condition de propreté (25). Ainsi, on obtient deux problèmes (structure et fluide) qui sont découplés. Une première stratégie permettant d'obtenir des vecteurs complexes approchant la relation de propreté consiste donc à résoudre les deux problèmes indépendamment, c'est-à-dire à chercher W et Q, approximations de X et P , vérifiant :

1) W W T = 0 2) X -W minimale 3) QQ T = 0 4) P -Q minimale
La résolution de ce problème de minimisation sous contrainte ne pose pas de difficulté pour la partie structurale, une des procédures possibles de résolution étant celle décrite par Balmès. En ce qui concerne la partie fluide, la présence du terme ∆ empêche l'exploitation directe de ces résultats, mais en remarquant que Q∆Q T = (Q∆ 1/2 )(Q∆ 1/2 ) T , on peut appliquer la procédure proposée au vecteur

Q = Q∆ 1/2 = jQS -1 .
Cette première stratégie de modification des vecteurs propres complexes conduit à une amélioration sensible des résultats, mais n'est pas complètement satisfaisante dans la mesure où la partie "extra-diagonale" de la relation de propreté n'est pas prise en compte.

C. Sur-nettoyage des vecteurs propres complexes

La solution explicite du problème d'optimisation sous contrainte obtenue par Balmès étant valable pour tout problème mettant en jeu une condition du type xx T = 0, on peut utiliser pour vecteur d'optimisation :

x =   X P P ∆   (33)
Alors :

xx T =   XX T XP T X∆P T P X T P P T P ∆P T P ∆X T P ∆P T P ∆ 2 P T   (34)
Ainsi, on retrouve bien les 4 termes intervenant dans l'écriture de la condition de propreté vibro-acoustique (25), mais l'optimisation impose alors des contraintes supplémentaires au problème, qui sont P P T = 0 et P ∆ 2 P T = 0, conditions non requises a priori par la solution du problème. De plus, le processus d'optimisation va supposer P et P ∆ indépendants, ce qui revient à modifier les valeurs propres complexes en plus des vecteurs propres complexes.

Néanmoins, nous verrons dans la suite que cette procédure conduit à de bons résultats sur les cas-tests traités.

V. ILLUSTRATION : NETTOYAGE DE VECTEURS PROPRES COMPLEXES VIBRO-ACOUSTIQUES PAR LES PROCÉDURES PROPOSÉES

Dans les trois cas-tests proposés ici, la procédure utilisée est identique :

1) A partir d'une données de fonctions de réponses en fréquences, des modes complexes sont identifiés par la technique de lissage de fonction de transfert [START_REF] Fillod | Research method of the eigenmodes and generalized elements of a linear mechanical structure[END_REF] Ainsi toute la matrice a été nettoyée, mais l'algorithme a en plus imposé des conditions supplémentaires non requises en théorie.

A partir des vecteurs propres nettoyés, la reconstruction des matrices est alors possible. Les matrices reconstruites sont indicées D pour la reconstruction directe, G pour la reconstruction à partir du nettoyage diagonal, et S pour le sur-nettoyage.

On obtient finalement pour les matrices de masse :

M D =
3.2954 -0.0005 -1.6777 0.0138

M G = 3.2947 -0.0005 -1.6781 0.0138 M S = 3.2987 -0.0005 -1.6825 0.0138 (43)
Les matrices de raideurs reconstruites sont :

K D = 992.
Et les matrices d'amortissement valent :

B D = -1.5495 -0.0056 -0.6357 -0.0072 B G = 2.1676 -0.0031 -1.4051 -0.0004 B S = 3.2946 -0.0080 -0.8763 -0.0015 (45)
On s'aperçoit que les trois reconstructions fournissent quasiment les mêmes matrices de masse et de raideur, en particulier pour les valeurs diagonales des matrices. Par contre, elles ne vérifient pas la topologie du système vibro-acoustique décrit par l'équation [START_REF] Balmes | New results on the identification of normal modes from experimental complex ones[END_REF]. En ce qui concerne la reconstruction de la matrice d'amortissement, son principal défaut concerne la valeur négative de l'amortissement B f identifiée dans les trois cas. Par contre, le nettoyage a permis de faire passer le terme B S d'une valeur négative dans le cas de la reconstruction directe à une valeur positive si les vecteurs complexes sont nettoyés. Notons cependant que le point de départ est constitué de vecteurs complexes bruités, et donc non forcément représentables par un système physique du second ordre. L'objectif du nettoyage est en fait de trouver le système physique ayant le comportement le plus proche du système caractérisé par la base modale complexe identifiée.

Une remarque importante est que la procédure de reconstruction des matrices (en particulier pour la matrice B) ainsi que la relation de propreté sont très sensibles au bruit sur les vecteurs complexes, puisque de faibles écarts sur ces vecteurs propres conduisent à de fortes différences sur les matrices.

Une fois ces matrices identifiées, et afin de comparer l'efficacité de la reconstruction, il est possible de calculer les fonctions de réponse en fréquence de chacun des systèmes reconstruits, et de les comparer au point de départ, c'est-à-dire à la recomposition modale des deux modes complexes bruités initiaux. Le tracé de la fonction de transfert structure-structure (mais prenant en compte le comportement couplé du système) est présenté sur la figure [START_REF] Balmes | New results on the identification of normal modes from experimental complex ones[END_REF], sur laquelle on peut observer que la technique de reconstruction directe ne permet pas de retrouver une évolution fréquentielle convenable (notamment en phase), tandis que les deux techniques de nettoyage donnent satisfaction, surtout la seconde technique, dite "sur-nettoyage" qui permet de reconstruire une FRF très proche de l'originale, ce qui est le but de l'opération. Rappelons cependant que les matrices ainsi déterminées n'ont pas une forme complètement en accord avec la topologie du système vibro-acoustique. Une campagne de mesure a été menée sur une guitare, en utilisant la technique d'excitation par marteau de choc. Un accéléromètre a été placé sur la table d'harmonie et un microphone mesure la pression acoustique au niveau de l'évent de l'instrument. Ce couplage diffère de la formulation proposée ici puisqu'il met en relation directe l'air de la cavité et le milieu extérieur. Néanmoins, en première approximation, on peut assimiler l'existence du milieu ouvert par une condition limite de pression nulle, qui est compatible avec la formulation utilisée. Dans la référence [START_REF] Foltête | Analyse modale expérimentale de systèmes vibro-acoustiques : application aux modes A0 et T1 de la harpe[END_REF], une technique de reconstruction des matrices utilisant la technique des moindres carrés a été proposée, fonctionnant très bien sur le cas de la guitare proposé ici, mais ne donnant pas satisfaction dans le cas de la harpe, qui est présenté plus loin dans cet article.

Les modes complexes ont donc été identifiés à partir des mesures, et les techniques de nettoyage ont été appliquées, pour comparer ensuite les FRF obtenues dans chaque cas considéré. Les résultats sont présentés sur la figure [START_REF] Wyckaert | Vibro-acoustical modal analysis : Reciprocity, model symmetry and model validity[END_REF], où l'on s'aperçoit que l'inversion directe du système conduit à des matrices erronées (notamment pour l'amortissement), tandis que le nettoyage permet d'améliorer de façon très efficace la reconstruction des fonctions de transfert.

Ce jeu de matrices permet de représenter correctement les FRF obtenues à partir des modes identifiés, en respectant la relation de propreté, mais pas exactement la topologie des matrices vibro-acoustiques. La prochaine étape dans ce travail consistera donc à déterminer une méthodologie permettant de respecter à la fois la relation de propreté vibro-acoustique et la topologie des matrices reconstruites.

C. Cas-test 2 : exploitation de données expérimentales issues d'une campagne de mesure sur une harpe

Ce dernier cas-test porte sur un jeu de mesures portant sur une harpe, pour laquelle la table d'harmonie était excitée avec un pot vibrant, et les mesures effectuées par vélocimétrie laser sur la table d'harmonie pour la partie structure, et par microphone situé dans un évent pour la partie acoustique. Ce cas avait posé problème lors de la reconstruction par la technique des moindres carrés, comme indiqué dans la référence [START_REF] Foltête | Analyse modale expérimentale de systèmes vibro-acoustiques : application aux modes A0 et T1 de la harpe[END_REF].

La mise en oeuvre des procédures de nettoyage proposées ici permet d'obtenir de bons résultats en ce qui concerne la reconstruction des FRF obtenues à partir des modes complexes identifiés, comme illustré sur la figure [START_REF] Morand | Fluid Structure Interactions[END_REF]. (47) Le système, s'il donne satisfaction, possède le même inconvénient que celui obtenu avec les mesures de la guitare : il ne respecte pas la topologie de la formulation utilisée pour le couplage vibro-acoustique.

VI. CONCLUSION

Dans ce papier, nous avons étendu la notion de propreté des modes complexes au cas du couplage vibro-acoustique. La condition de propreté a été écrite dans ce cas, son expression différant de celle existant en dynamique des structures car le système est ici non auto-adjoint. Une technique efficace de "nettoyage" des vecteurs complexes existant pour la structure, son extension au cas du couplage vibro-acoustique a été tentée, conduisant à un problème de minimisation sous contrainte. Bien que ce problème n'ait pas été résolu à ce jour, nous avons proposé ici deux approches permettant un nettoyage partiel (ou trop contraignant) des vecteurs complexes. Ces deux approches ont été appliquées sur un cas-test numérique et sur deux cas-tests expérimentaux, et dans tous les cas ont permis de retrouver les bonnes fonctions de réponse en fréquence après nettoyage des vecteurs, tandis que dans tous les cas testés la reconstruction directe des matrices caractéristiques mène à un système dont le comportement est très éloigné du système initial. En ce sens, les procédures proposées sont efficaces. Malheureusement, les matrices reconstruites ne vérifient pas toujours la topologie des matrices constitutives de la formulation vibro-acoustique, ce qui nous fournit le prochain objectif de l'étude : améliorer les procédures de modification présentées ici afin d'obtenir un système matriciel vérifiant la topologie spécifique du couplage vibro-acoustique.

REMERCIEMENTS

Les auteurs souhaitent remercier chaleureusement Jean-Loïc Le Carrou et François Gautier du Laboratoire d'Acoustique de l'Université du Maine, qui ont effectué l'ensemble des mesures utilisées pour l'identification des modes complexes présentés ici.

  25) IV. PROCÉDURES DE MODIFICATIONS DES VECTEURS PROPRES COMPLEXES AFIN QU'ILS VÉRIFIENT LA CONDITION DE PROPRETÉ

Fig. 1 .

 1 Fig. 1. Comparaison des FRF calculées à partir des vecteurs complexes nettoyés suivant les procédures proposées : cas-test simulé bruité

Fig. 2 . 1 B = 3 .

 213 Fig. 2. Comparaison des FRF calculées à partir des vecteurs complexes nettoyés suivant les procédures proposées : cas des mesures sur une guitare

Fig. 3 .

 3 Fig. 3. Comparaison des FRF calculées à partir des vecteurs complexes nettoyés suivant les procédures proposées : cas des mesures sur une harpe

  1 ρ f Ωf ∇p∇δpdΩ. Cette matrice est symétrique, définie, semi-positive. -C est la matrice de couplage vibro-acoustique, associée à la forme bilinéaire Sc pδu n dS. Cette matrice est rectangulaire si le nombre de degrés de liberté du fluide est différent du nombre de degrés de liberté structuraux.

-B s et B f représentent respectivement les pertes structurales et acoustiques. Cette formulation suppose qu'il n'y a pas de perte au niveau de la surface de couplage entre les deux milieux et que les dissipations peuvent être modélisées par des effets visqueux. -F s est le vecteur représentant les forces généralisées exercées sur la structure.

  Pour ces valeurs numériques, les valeurs propres complexes sont s 1 = -0.1065 + 10.1630j et s 2 = -0.1920 + 19.6811j, et la matrice modale de droite est :Remarquons que pour être cohérent avec les notations utilisées auparavant, la matrice dY devrait avoir 2 lignes et 4 colonnes, où les deux dernières colonnes sont constituées des vecteurs conjugués. Nous avons volontairement omis cette partie de la matrice afin d'alléger les notations et de simplifier la lecture du document.Ces vecteurs sont alors bruités numériquement, avec un bruit aléatoire de ±10% sur les parties réelle et imaginaire de chaque composante, de densité de probabilité uniforme. Dans le cas traité ensuite, les vecteurs initiaux sont après bruitage (on indice D les grandeurs associées à la reconstruction directe) :

							. Ces
	modes complexes, composés de valeurs propres com-
	plexes et de vecteurs propres complexes, sont utilisés
	comme point de départ des calculs effectués ensuite.
	2) Les modes complexes identifiés sont utilisés pour re-
	construire directement les matrices caractéristiques du système à l'aide des relations (22) à (24). Ces matrices sont alors employées pour calculer la réponse fréquen-	dY D =	-0.0243 + 0.0278j 0.0582 -0.0614j 11.4021 -11.4469j 9.9273 -11.2118j	(37)
	tielle du système reconstruit, qui est comparée à la	Ce niveau de bruit est important, mais il permet de tester
	réponse fréquentielle initiale. Ce cas de reconstruction	la robustesse de la procédure présentée ici. La propreté de ce
	est appelé "Direct", les différences entre les fonctions	jeu de vecteurs est alors :
	de réponse en fréquence (FRF) étant dues aux termes associés à la non-vérification de la relation de propreté : dY gY T = 0 dans les équations (18) et (20).		dY D gY T D =	-0.0011 0.0003 -0.1394 -0.2883	(38)
	3) Les vecteurs complexes de départ sont nettoyés par la	La mise en oeuvre de la technique de nettoyage n˚1 permet
	procédure de nettoyage "diagonal", puis les matrices ca-	d'obtenir les vecteurs "nettoyés" suivants (on indice G les
	ractéristiques sont reconstruites, et utilisées calculer	grandeurs relatives à cette technique de nettoyage) :
	les FRF du système, qui sont comparées aux réponses
	fréquentielles initiales. Ce cas de reconstruction est appelé "Diagonal". 4) Les vecteurs complexes de départ sont nettoyés par la procédure de "sur-nettoyage", puis les matrices carac-téristiques sont reconstruites, et utilisées pour calculer	-0.0246 + 0.0275j 11.3037 -11.5944j -10.1147 + 11.0071j -0.0605 + 0.0592j et dans ce cas la propreté devient : dY G =	(39)
	les FRF du système, qui sont comparées aux réponses fréquentielles initiales. Ce cas de reconstruction est		dY G gY T G =	0 -0.0642 5.0 × 10 -16 5.3 × 10 -5	(40)
	appelé "Sur-nettoyage".				On observe bien que seule la partie diagonale de la matrice
	A. Cas-test 1 : simulation numérique			a été nettoyée. La mise en oeuvre de la technique de nettoyage n˚2 permet
	Pour ce premier cas-test, nous utilisons des données simu-	d'obtenir les vecteurs "nettoyés" suivants (on indice S les
	lées, en partant d'un système du second ordre à deux degrés de	grandeurs associées à cette technique de nettoyage) :
	liberté dont la topologie respecte celle d'un problème vibro-
	acoustique couplé : 3.234 -1.465 1.275 × 10 -2 0	ẍ p		0.0242 -0.0277j -11.7433 + 11.1237j 10.2244 -10.8952j 0.0606 -0.0591j et dans ce cas la propreté devient : dY S =	(41)
		+	1.123 0	0 3.185 × 10 -3	ẋ ṗ	(35)	dY S gY T S =	-9.8 × 10 -17 -3.3 × 10 -16 -1.3 × 10 -14 -2.3 × 10 -14	(42)
		+	1000 1.465 0 1.65	x p	= 0
	dY =	-0.0266 + 0.0263j 12.1474 -11.8894j 10.3920 -10.2175j 0.0597 -0.0598j	(36)


Thanh Phuong Nguyen
email: nguyentp@loria.fr

Isabelle Debled-Rennesson

Curvature estimation in noisy curves

come

Introduction

A lot of applications in image processing requires the geometrical measuring of represented discrete objects.

In the framework of the discrete geometry, estimators of geometrical parameters have been proposed but they rely on the recognition of discrete line segments which is very sensitive to the noise existing in the studied curves [START_REF] Feschet | Optimal time computation of the tangent of a discrete curve: Application to the curvature[END_REF][START_REF] Coeurjolly | Discrete curvature based on osculating circle estimation[END_REF][START_REF] Lachaud | Analysis and comparative evaluation of discrete tangent estimators[END_REF]. The boundary of such discrete objects is often noisy due to acquisition process. Therefore the concept of blurred segment was introduced [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF], it allows the flexible segmentation of discrete curves, taking into account noise. Relying on an arithmetic definition of discrete lines [START_REF] Reveillès | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF], it generalizes such lines, admitting that some points are missing. A curvature estimator was proposed in [START_REF] Debled-Rennesson | Estimation of tangents to a noisy discrete curve[END_REF] and used in an application to the arc detection in technical documents [START_REF] Salmon | A new method to detect arcs and segments from curvature profiles[END_REF], the complexity of this curvature estimator algorithm is in O(n 2 ) (n is the number of points of the studied curve) and it can only be applied to 8-connected simple curves. We proposed in this paper an extension to general curves of this algorithm. First the recognition algorithm of blurred segments proposed in [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] is extended to general curves [START_REF] Buzer | An elementary algorithm for digital line recognition in the general case[END_REF]; the problem of adding or removing a point is studied. Then thanks to the decomposition of curve into maximal blurred segments, we proposed an algorithm of calculation of the curvature at each point of a general discrete curve in O(nlog 2 n).

The paper is organized as follows. In Section 2, after recalling definitions related to blurred segments, we study the problem of adding (or removing) a point to (from) a blurred segment of width ν in the case of a general discrete curve. Then we propose an extension to the noisy curves of the notion of maximal segment of a discrete curve. An algorithm to determine all maximal blurred segments of a discrete curve is given in Section 3. In Section 4, after recalling the definition of the curvature estimator adapted to noisy curves, we proposed a new algoritm for the determination of the curvature at each point of a discrete curve. Examples are also given.

2

Blurred segment of width ν

Definitions

The notion of blurred segments relies on the arithmetical definition of discrete lines [START_REF] Reveillès | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF] where a line, whose slope is a b , lower bound µ and thickness ω (with a, b, µ and ω being integer such that gcd(a, b) = 1) is the set of integer points (x, y) verifying µ ≤ axby < µ + ω. Such a line is denoted by D(a, b, µ, ω). Let us recall definitions [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] that we use in this paper (see Figure 1 In this section we study the problem of adding (or removing) a point to (from) a blurred segment of width ν. The algorithm of recognition of blurred segments of width ν presented in [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] is executed in linear time however it only considers the incremental addition of a point to a blurred segment in the first octant. We present here the general case which requires the incremental calculation of both height and width of the convex hull after adding or removing a point. To do that, we use the results given in [START_REF] Overmars | Maintenance of configurations in the plane[END_REF] and [START_REF] Buzer | An elementary algorithm for digital line recognition in the general case[END_REF] that we briefly recall below.

Dynamic estimation of the convex hull

The problem of dynamic estimation of the convex hull of a set of points when adding (or removing) a point to (from) this set was proposed by M.H. Overmars, J. van Leeuwen [START_REF] Overmars | Maintenance of configurations in the plane[END_REF]. In this work, a convex hull is considered as the union of two parts: the upper convex hull (U hull ) and the lower convex hull (L hull ). U hull and L hull are updated after each operation of addition or removal of a point, the cost of these operations are estimated by the following theorem [START_REF] Overmars | Maintenance of configurations in the plane[END_REF].

Theorem 1 The convex hulls U hull and L hull of the set S of n points may be dynamically kept, in the worst case, in O(log 2 n) by an operation of addition or removal.

Determination of height and width of the convex hull

We use the double technique of binary search [START_REF] Buzer | An elementary algorithm for digital line recognition in the general case[END_REF] to determine the height and width of the convex hull. In [START_REF] Buzer | An elementary algorithm for digital line recognition in the general case[END_REF], the convex hull is also considered as the union of two parts U hull and L hull . The double technique of binary search permits to find the vertical width of the convex hull by using the concavity property of the function height

(x) = U hull (x) -L hull (x) in O(log 2 n).
3 Maximal blurred segments of width ν

Definitions and first property

The notion of maximal segment of a discrete curve was proposed in [START_REF] Feschet | Optimal time computation of the tangent of a discrete curve: Application to the curvature[END_REF][START_REF] Lachaud | Analysis and comparative evaluation of discrete tangent estimators[END_REF] and relies on the discrete line segments. This structure enables a global understanding of the discrete curve to be analyzed. We propose here an extension of that notion to the blurred segments, adapted to noisy curves, by using the same notations as in [START_REF] Lachaud | Analysis and comparative evaluation of discrete tangent estimators[END_REF].

Let us consider a discrete curve called C, the points of C are indexed from 0 to n -1. C is a general curve and the points of C can be disconnected . We note C i,j a set of successive points of C increasingly ordered from index i to j. Definition 4. The predicate "C i,j is a blurred segment of width ν" is denoted by BS(i, j, ν). The first index j, i ≤ j, such that BS(i, j, ν) and ¬BS(i, j + 1, ν) is called the front of i and noted F (i). Symmetrically, the first index i such that BS(i, j, ν) and ¬BS(i -1, j, ν) is called the back of j and noted B(j). Definition 5. C i,j is called a maximal blurred segment of width ν and noted M BS(i, j, ν) iff BS(i, j, ν) and ¬BS(i, j + 1, ν) and ¬BS(i -1, j, ν).

It is obvious that an equivalent characterization for a maximal blurred segment of width ν, M BS(i, j, ν), is to show that F (i) = j and B(j) = i. In this work, we use the notion of blurred segment of width ν which is maximal on the right or left sides: Definition 6. C i,j is called a maximal blurred segment of width ν on the right side (resp. on the left side) and noted M BS R (i, j, ν)

(resp. M BS L (i, j, ν)) if F (i) = j (resp. B(j) = i). Property 1 Let C be a discrete curve, M BS ν (C) the sequence of maximal blurred segments of width ν of the curve C. M BS ν (C) = {M BS(B 1 , E 1 , ν), M BS(B 2 , E 2 , ν), ..., M BS(B m , E m , ν)}, it sat- isfies B 1 < B 2 < ... < B m . So we have: E 1 < E 2 < ... < E m . Proof: We consider 2 consecutive maximal blurred segments M BS(B i , E i , ν) and M BS(B i+1 , E i+1 , ν). By hypothesis, B i < B i+1 , let us suppose that E i ≤ E i+1 , then M BS(B i+1 , E i+1 , ν) becomes a part of M BS(B i , E i , ν). Therefore M BS(B i+1 , E i+1 , ν
) is not a maximal blurred segment, that is contradictory.

Algorithm for the segmentation of a curve C into maximal blurred segments

We propose an algorithm (see Algorithm 1) which determines all maximal blurred segments of width ν of a discrete curve C according to the conditions given in section 3.1. To do that, property 1 is used.

Complexity

Each point of the curve is scanned at most twice in this algorithm. The cost of determining a new optimal bounding discrete line when we add (or remove) a point to (from) a blurred segment is in O(log 2 n). Hence the complexity of this algorithm is in O(nlog 2 n).

4 Discrete curvature of width ν

Definition

We recall in this section the curvature estimator which is adapted to noisy curves [START_REF] Debled-Rennesson | Estimation of tangents to a noisy discrete curve[END_REF]. It is directly deduced from the estimator proposed by D. Coeurjolly [START_REF] Coeurjolly | Discrete curvature based on osculating circle estimation[END_REF] for the 2D curves without noise. This technique can be seen as a generalization of the classical order m normalized curvature [START_REF] Rosenfeld | Angle detection on digital curves[END_REF].

Let C be a discrete curve, C k is a point of the curve. Let us consider the points C l and C r of C such that : l < k < r, BS(l, k, ν) and ¬BS(l -1, k, ν), BS(k, r, ν) and ¬BS(k, r + 1, ν).

The estimation of the curvature of width ν at the point C k shall be determined thanks to the radius of the circle passing through the points C l , C k and C r . To determine the radius R ν (C k ) of the circumcircle of the triangle [C l , C k , C r ], we use the formula given in [START_REF] Harris | Handbook of mathematics and computational science[END_REF] as follows (see Figure 2.a):

Let s 1 = || ---→ C k C r ||, s 2 = || ---→ C k C l || and s 3 = || ---→ C l C r ||, then R ν (C k ) = s 1 s 2 s 3 (s 1 + s 2 + s 3 )(s 1 -s 2 + s 3 )(s 1 + s 2 -s 3 )(s 2 + s 3 -s 1 )
Then, the curvature of width ν at the point [START_REF] Coeurjolly | Discrete curvature based on osculating circle estimation[END_REF], the degenerated cases, which correspond for example to colinear half-tangents, may be independently tested and, thus, a null curvature is affected to the considered point.

C k is C ν (C k ) = s Rν (C k ) with s = sign(det( ---→ C k C r , ---→ C k C l )) (it

indicates concavities and convexities of curve). As indicated in

B i E i B i+1 E i+1 (b) Ei (Bi+1) is front (back)
of points in first (second) bold edge Let us then consider the decomposition of C into maximal blurred segments:

M BS ν (C) = {M BS(B 1 , E 1 , ν), M BS(B 2 , E 2 , ν), ..., M BS(B m , E m , ν)} with B 1 < B 2 < ... < B m and E 1 < E 2 < ... < E m .
We look for the indices i and j such that i is the first index such that E i ≥ k and j is the last index such that B j ≤ k. So it is obvious that l = B i , r = E j and that the curvature of width ν at the point C k is the inverse of the radius of the circumcircle of the triangle [C l , C k , C r ]. More generally, we have the following simple result: Property 2 Let L(k), R(k) be the functions which respectively present the indices of the left and right extremities of the maximal blurred segments on the left and right sides of the point C k .

-∀k such that E i-1 < k ≤ E i , then L(k) = B i -∀k such that B i ≤ k < B i+1 , then R(k) = E i
This method is used in the algorithm 2 (see Figure 2.b).

Complexity

Both steps of labelling and estimation of the curvature at each point are executed in linear time. However, the determination of the maximal blurred segments are executed in O(nlog 2 n). Thus the complexity of our method is in O(nlog 2 n). It is more efficient that O(n 2 log 2 n) when we work with general curves as well as O(n 2 ) for the simple curves with the existing methods.

m = |M BSν |; E-1 = -1; Bm = n; for i = 0 to m -1 do for k = Ei-1 + 1 to Ei do L(k) = Bi; for k = Bi to Bi+1 -1 do R(k) = Ei; end for i = 0 ( * ) to n -1 ( * ) do Rν (Ci) = Radius of the circumcircle to [C L(i) , Ci, C R(i) ]; Cν(Ci) = s Rν (C i ) ; end end
(*) The bounds mentioned in the algorithm are correct for a closed curve. In case of an open curve, the instruction becomes: For i = l to n -1 -l with l fixed to a constant value. Indeed it is not possible to calculate a maximal blurred segment on the left side (resp. on the right side) at the first point (resp. at the last point) of the curve. Thus the calculation of the curvature begins (resp. stops) at the l th (resp. (n -1l) th ) point of the curve.

Results

Two discrete curves

Conclusion

We have proposed the notion of maximal blurred segments of a discrete curve for a given width as an extension to noisy curves of the definitions proposed in [START_REF] Feschet | Optimal time computation of the tangent of a discrete curve: Application to the curvature[END_REF][START_REF] Lachaud | Analysis and comparative evaluation of discrete tangent estimators[END_REF]. This decomposition of a discrete curve enabled us to propose an optimized version of the algorithm of calculation of the curvature of width ν at each point of a general discrete curve. However, we expect to obtain a specific algorithm with a better complexity for simple curves where points are added in one direction. The extension of the results of this paper to 3D discrete curves will be subject to a next publication.

): Definition 1 .Fig. 1 .

 11 Fig. 1. D(5, 8, -8, 11), optimal bounding line (vertical distance = 10 8 = 1.25) of the sequence of gray points

Algorithm 1 :

 1 Algorithm for the segmentation of a curve C into maximal blurred segments of width ν Data: C -discrete curve with n points, ν -width of the segmentation Result: M BSν -the sequence of maximal blurred segments of width ν begin k=0; S b = {C0}; M BSν = ∅; a = 0; b = 1; ω = b, µ = 0; while ω-1 max(|a|,|b|) ≤ ν do k++; S b = S b ∪ C k ; Determine D(a, b, µ, ω) of S b ; end bSegment=0; eSegment=k-1 ; M BSν = M BSν ∪ C bSegment,eSegment ; while k < n -1 do while ω-1 max(|a|,|b|) > ν do bSegment++ ; S b = S b \ C bSegment ; Determine D(a, b, µ, ω) of S b ; end while ω-1 max(|a|,|b|) ≤ ν do k++ ; S b = S b ∪ C k ; Determine D(a, b, µ, ω) of S b ; end eSegment=k-1; M BSν = M BSν ∪ C bSegment,eSegment ; end end

4. 2

 2 Algorithm for the estimation of the curvature of width ν at each point of CWe propose in this section a new algorithm for the determination of the curvature of width ν at each of the n points of a curve C. The complexity of this algorithm is better than the one of the naive algorithm, in O(n 2 ), which consists in calculating at each point C k , the maximal blurred segment on the right side, M BS R , the maximal blurred segment on the left side, M BS L , then the circle Estimation of the curvature at the point T with width 2

Fig. 2 .

 2 Fig. 2.

Algorithm 2 :

 2 Width ν curvature estimation at each point of C Data: C discrete curve of n points, ν width of the segmentation Result: {Cν (C k )} k=0..n-1 -Curvature of width ν at each point of C begin Build M BSν = {M BS(Bi, Ei, ν)} ;

  (3.a and 3.c) are represented on the Figure 3 with the graph of their curvature values calculated at each point of the curves with width 2 (3.b and 3.d). The points of the curve 3.c corresponding to the peaks of the associated curvature graph 3.d are indicated by black pixels.