Pascal Bourdon
email: pascal.bourdon@univ-poitiers.fr

David Helbert

A DATA-DRIVEN APPROACH TO FEATURE SPACE SELECTION FOR ROBUST MICRO-ENDOSCOPIC IMAGE RECONSTRUCTION

Keywords: Medical imaging, microendoscopy, optical flow, image mosaicing

In the article we propose a new on-line feature space selection strategy for displacement field estimation in the context of multi-view reconstruction of biological images acquired by a multi-photon micro-endoscope. While the high variety of targets encountered in clinical endoscopy induce enough texture feature variability to prohibit the use of recent supervised learning or feature matchingbased visual tracking methods, we will show how on-line learning combined with a classical method such as Digital Image Correlation (DIC) can contribute to the improvement of convex optimization-based template matching techniques.

INTRODUCTION

Matching visual content across images or videos is a major research topic in computer vision, with many applications such as velocimetry, multi-view image reconstruction, or human motion analysis [START_REF] Thirion | Image matching as a diffusion process: an analogy with maxwell's demons[END_REF][START_REF] Baumberg | Reliable feature matching across widely separated views[END_REF][START_REF] Riethmuller | Particle image velocimetry[END_REF][START_REF] Sotiras | Deformable medical image registration: A survey[END_REF][START_REF] William R Crum | Non-rigid image registration: theory and practice[END_REF][START_REF] Arnold | Visual tracking: An experimental survey[END_REF]. While many papers focus on supervised learning techniques to address this problem in the context of reproducible scenarios such as face analysis or vehicule tracking, specific applications such as Particle Image Velocimetry (PIV) remain a challenge. A typical example is medical imaging, where strong acquisition noise combined with non-rigid distortions induce feature-point matching systems, such as the wellknown SIFT/FLANN association, into errors. As a result, classical techniques based on cross-correlation maximization (e.g. DIC) or mean square error minimization (e.g. Lucas-Kanade [START_REF] Bruce D Lucas | An iterative image registration technique with an application to stereo vision[END_REF]) are still very popular, despite certain drawbacks we will discuss later.

In recent years, multi-photon microendoscopy has become an essential tool in cell and tissue biology research

The authors would like to thank the XLIM-BIOElectroPhot team for their support, assistance, and most of all for introducing us to this fascinating and challenging subject matter. [START_REF] Bedard | Real-time video mosaicing with a high-resolution microendoscope[END_REF][START_REF] Hoover | Advances in multiphoton microscopy technology[END_REF][START_REF] Ducourthial | Development of a real-time flexible multiphoton microendoscope for label-free imaging in a live animal[END_REF]. The benefits are its ability to achieve the optical sectioning, high resolution of thick samples (2D image with a low depth of field), high contrast, low sensitivity to the diffusion of biological material, high penetration depth, minimal photo-toxicity and the possibility to dispense with exogenous carcinogen marking. It is another application where the classical DIC method is usually preferred over other methods to estimate displacements prior to the multi-view reconstruction of extracted image data [START_REF] Guizar-Sicairos | Efficient subpixel image registration algorithms[END_REF].

The displacement estimation method presented in this paper was specifically developed to deal with low field-of-view images acquired with a two-photon microendoscope [START_REF] Ducourthial | Development of a real-time flexible multiphoton microendoscope for label-free imaging in a live animal[END_REF]. Such images are obtained by scanning pixel intensity over time along a spiral path. As a result, together with noise, illumination changes, and rigid transformations induced by the hand motion of clinicians operating the endoscope, non-rigid deformations are also generated, mostly in the forms of a fish-eye effect with an additional defective whirl effect in the center of the image, making visual matching a challenging task (see Fig. 1).

Fig. 1. Consecutive micro-endoscopic acquisitions

While frequency-domain methods such as DIC work very well for rigid transformations, in the future we intend to rely on non-rigid transformation models. As a result, we decided to study model-based methods such as Baker-Matthews' Inverse Computational Image Alignment (ICIA) [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF], which is a generalization of Lucas-Kanade [START_REF] Bruce D Lucas | An iterative image registration technique with an application to stereo vision[END_REF]. Fur-thermore, our application requires sub-pixel accuracy, which is natural with ICIA, as opposed to DIC. This however does not mean that DIC is useless to our application. It will be used for ground-truth generation in an on-line learning context, as it is able to address the issue of flow estimation over simple yet large displacements.

The next section of this paper addresses, in terms of the theoretical framework proposed by Baker-Matthews, the issue of image alignment. We will show why the registration of noisy and distorted images cannot always be treated as a conventional convex optimization problem. In section 3, we will discuss the advantages of a feature space approach, and explain how we intend to use DIC results in a feature selection context. Finally, in section 4 we will present experimental results obtained on actual micro-endoscopic image data, and open a discussion about the advantages and drawbacks of this method.

VISUAL MATCHING AND CONVEX OPTIMIZATION

Theoretical framework

The problem of recovering a large field-of-view 2D image projection of a scene using low-field acquisitions can be written as follows : at time t , the endoscope can only capture a projection of said scene on a subset Ω ∈ ℜ 2 . Let function I represent this projection for any x position:

I : Ω → ℜ x → I (x, t ). ( 1
)
Given a small enough deformation or hand displacement during acquisition, there exists a geometric transform W(x, θ t ) associated with a subset Υ ∈ Ω such that:

I (x, t + ∂ t ) = I (W (x, θ t ) , t ) ∀x ∈ Υ , ( 2
)
where ∂ t is the time between two consecutive images, and where W(x, θ) denotes a parametrized set of warps that takes a pixel in the coordinate frame of I (x, t ) and maps it to a new, possibly sub-pixel, position in the coordinate frame of I (x, t +∂ t ). W(.) can be a free-form function or be modelled after rigid or non-rigid transforms. Subset Υ can either define one or several sets of connected positions, in the case of block or template matching, or a set of sparse positions, called keypoints, where function I is often replaced by its projection onto a new representation domain called feature space f(I ), in the case of feature matching. In all cases, image registration consists in estimating a set of parameters θ t which satisfy Eq. ( 2), or more likely finding the minimum of the squared residuals:

θt = arg min θ t (ξ (θ t )) ξ (θ t ) = x∈Υ I (x, t ) -I (W (x, θ t ) , t + ∂ t ) 2 2 .
(3)

Practical solvability

Objective function ξ (θ t ) is often considered to be convex and smooth enough for 1 s t or 2 nd order derivatives to exist (either natively or by filtering input images). Therefore Eq. ( 3) can be solved by taking the partial derivatives of ξ (θ t ) and setting them equal to zero in a gradient-descent fashion such as the one introduced by Lucas-Kanade:

∂ ξ (θ t ) ∂ θ t θ t = θt ≈ 0, ( 4
)
where θt → θt is the expected output.

Fig. 2 shows the result of trying to match a small textured patch between two micro-endoscopic images, using a simple translation model (θ t ∈ ℜ 2 ) such as:

W (x, θ t ) = x + θ t .
(5) Two different convex optimization algorithms have been tested and compared to global cross-correlation maximization: the Conjugate Gradient (CG) method, and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. As one can notice on Fig. 2, both methods fail to match the appropriate position in the target image. From a strictly mathematical point of view, it is easy to understand that both methods failed because assumptions made about the convexity of cost function ξ (θ t ) are inaccurate. Indeed, both CG and BFGS converged into local minima, neither of which are a global minimum. From a visual point of view, we can notice that the two matched image patterns could be good matches if it wasn't for the pattern matched globally. In other words, the current representation domain of function I is not discriminative enough to distinguish between visually acceptable matches and good, unique matches (e.g. see Zhu et al.'s representation model for unsupervised feature selection [START_REF] Zhu | Unsupervised feature selection by regularized self-representation[END_REF]).

DATA-DRIVEN FEATURE SPACE SELECTION

Feature matching

Feature-matching is often defined as establishing a correspondence between distinct points in images i.e. keypoints [START_REF] Baumberg | Reliable feature matching across widely separated views[END_REF]. One can however define a feature space where every pixel within the image has a dual representation (e.g. Gabor jets or differences of Gaussians). We define the following vector-valued function:

F (x, t ) = f (I (x, t )) =    f 1 (I (x, t )) f 2 (I (x, t )) • • • f N (I (x, t ))    , (6)
where F (x, t ) is the new representation of function I in a feature space made of a given set of N features f(.) extracted from image I (x, t ). Eq. ( 3) can be written as:

θt (β) = arg min θ t ξ β (θ t ) = θt ξ β (θ t ) = x∈Υ β T . (F (x, t ) -F (W (x, θ t ) , t + ∂ t )) 2 2 (7)
with β = β 0 β 1 . . . β N T a vector of weighting coefficients being quantitative expressions of the importance of each element of f(.), and where θt are the actual good-match warp parameters, recovered from exhaustive (brute-force) search or global minimization on the original image domain (Eq. 3).

On-line space selection

Let θt (β) be the set of warping parameters that, given a β configuration, provides us with an expected solution for the minimization of ξ (θ t ). As an example:

∂ ξ β (θ t ) ∂ θ t θ t = θt (β) ≈ 0. ( 8
)
Despite DIC not being able to estimate complex warping models, it still provides good results with simple transforms. Therefore we propose to use it for automatic feature selection during a pre-acquisition phase t ∈ [0, t 0 ], where the endoscope operator is required to capture a few images for ground truth generation, with DIC running and providing the algorithm with θt samples. While such a phase should be performed whenever substantial changes in the acquisition protocol occur (e.g. when new organs or tissues are targeted), we want to emphasize that this process does not require painful tasks such as manual annotation of ground truth data. The operator only has to capture a few images beforehand and let the algorithm self-calibrate itself. As a matter of fact, given sufficient computational resources, we expect initial DIC-based estimations to perform in real-time, thus providing reconstruction results from the very first captures (t = 0) before the optimal feature space-based estimator progressively takes over to reduce computer usage.

The estimation of β is written as:

β = arg min β α β 1 + t ∈[0,t 0 ] θt -θt (β) 2 2 , ( 9
)
where the 1 -norm constraint on β aims towards producing a sparse solution, depending on α values.

Looking for a sparse vector β has multiple purposes:

zeroing elements of β reduces the number of features to compute, hence computation time. It also acts as a safety against redundant or non-orthogonal descriptions, reducing the risk of extracting same-information about texture patterns multiple times and having to minimize a nonbijective function in Eq. ( 9).

EXPERIMENTAL RESULTS

The Photonics test data consists of 448 450 × 450 lowfield images of a mouse kidney specimen taken with a 250µm fiber endoscope. Ground truth data is generated using global cross-correlation maximization (DIC) to estimate 2D shifts on image patches centered around 205 point positions for each of the first 50 sub-images. Crosscorrelation peaks are estimated over full target sub-images in order to easily detect outliers, which are removed from the data set. Eq. ( 9) is a general 1 -minimization problem of the form arg min β β 1 + H (β), where quadratic function H (•) could be assumed as a convex and differentiable function to use an iterative shrinkage solver. However, given the preliminary and experimental nature of this work we decided to limit such assumptions and use the non-realtime stochastic Simulated Annealing (SA) technique for global optimization [START_REF] Khachaturyan | The thermodynamic approach to the structure analysis of crystals[END_REF].

The algorithm was tested with two well-known image feature spaces : scale-variant polynomials, using Gaussian (0 t h order) and Sobel (1 s t and 2 nd orders) kernels of sizes 3 × 3, 15 × 15 and 31 × 31, and Gabor wavelets [START_REF] John G Daugman | Complete discrete 2-d gabor transforms by neural networks for image analysis and compression[END_REF], using 4 equally-spaced orientations for spacial frequencies {0.05, 0.15} and scales (i.e. σ values) {3, 7}. The BFGS algorithm was used for optical flow estimation on original (raw pixel intensity) and transformed feature spaces.

Table 1 shows results in terms of Mean-Square Error (MSE) and Standard Deviation (SD) computations between estimated displacement fields and ground truth data over the full Photonics sequence. As we can see, our method outperforms intensity-based estimations, with the β-weighted polynomial feature space estimated with α = 1.0 providing the best results. It is very interesting to notice that a lower α value of 0.1, which releases the constraint on the 1 -norm of β, actually leads to worse results than a value of 1.0. We assume this is due to possible nonbijective properties of the cost function discussed earlier, especially since we already know that the initial feature spaces we chose to test are not orthogonal bases.

Feature space MSE SD

Raw pixel intensity Both Fig. 3 and Fig. 4 show additional results between pixel intensity and β-weighted, α = 1.0, polynomial feature space matching. Fig. 3 shows results for individual field estimations on a random selection of low-field images, while Fig. 4 presents an example of actual estimated displacement fields. Both figures confirm that estimations are improved using our method.

CONCLUSION

In this paper we presented a data-driven approach to feature space selection for robust micro-endoscopic image reconstruction. We first discussed how classical template matching techniques such as Digital Image Correlation are still preferred over learned trackers, that are unlikely to succeed with image defects and texture patterns typical of medical images. We also explained how, for the very same reasons, convex optimization-based methods applied in the raw pixel intensity domain can fail, converging into local minima when texture patterns are not unique enough within image data. As a result, we introduced a new feature space selection method to generate template matching-friendly patterns thanks to the on-line learning of ground truth data from non-optimal yet robust global methods such as DIC. Experiments on actual micro-endoscopic images have shown promising results and potential for going further into this research. In future work, we intend to integrate a fully rigid and non-rigid deformation model for global low-field image alignment (as opposed to multiple patch tracking) and improve the on-line learning procedure with a trained detector of data priors such as untrackable points or areas.

Fig. 2 .

 2 Fig. 2. 2D shift estimation: global, exhaustive (brute-force) search (yellow) succeeds while gradient descent (purple, green) falls into local minima of the cost function.

Fig. 3 .

 3 Fig. 3. MSE results between pixel intensity and β-weighted polynomial feature space (α = 1.0) for individual displacement fields on a random selection of low-field images.

Fig. 4 .

 4 Fig. 4. β-weighted polynomial feature space (c) has more accurate matches than pixel intensity (b).

Table 1 .

 1 Estimation

		12.67 ± 19.40
	Multi-scale polynomial, α = 0.1	11.05	±16.70
	Multi-scale polynomial, α = 1.0	10.37 ± 14.57
	Gabor wavelet basis, α = 0.1	11.21	±16.66
	Gabor wavelet basis, α = 1.0	11.05	±16.70

results pixel intensity and βweighted feature space representations for MSE and SD.


Matteo L Facchinetti

Xavier Boutillon

Andrei Constantinescu

Numerical and experimental modal analysis of the reed and pipe of a clarinet a)

A modal computation of a complete clarinet is presented by the association of finite-element models of the reed and of part of the pipe with a lumped-element model of the rest of the pipe. In the first part, we compare modal computations of the reed and the air inside the mouthpiece and barrel with measurements performed by holographic interferometry. In the second part, the complete clarinet is modeled by adjoining a series of lumped elements for the remaining part of the pipe. The parameters of the lumped-resonator model are determined from acoustic impedance measurements. Computed eigenmodes of the whole system show that modal patterns of the reed differ significantly whether it is alone or coupled to air. Some modes exhibit mostly reed motion and a small contribution of the acoustic pressure inside the pipe. Resonance frequencies measured on a clarinet with the mouthpiece replaced by the cylinder of equal volume differ significantly from the computed eigenfrequencies of the clarinet taking the actual shape of the mouthpiece into account and from those including the ͑linear͒ dynamics of the reed. This suggests revisiting the customary quality index based on the alignment of the peaks of the input acoustical impedance curve.

I. INTRODUCTION

The clarinet is usually considered as the association of a linear resonator, the pipe, and a nonlinear excitor, the reed, subject to the air flow from the mouth. Alternatively, one can consider the air column and the reed as a linear system subject to nonlinear boundary conditions. This is the approach retained in this article where the reed is considered as a linear mechanical system coupled to the pipe and where the interaction with the player is not treated. Nonlinear phenomenon such as the interaction between the reed and the jet across the reed-slit, the contact forces between the reed and the lay, and the interaction between the reed and the player's lip will be included in a subsequent piece of work as nonlinear boundary conditions to the normal modes that are described here. Humidity of the reed and the player's lip also have a damping role which is not considered in this modal analysis of a pipe coupled to a ͑dry͒ reed.

Acoustical studies of the clarinet have so far represented the mouthpiece of a wind instrument by its equivalent volume. This study goes beyond this approximation and presents the three-dimensional distribution of the pressure in the upper part of the instrument.

Studies of the pipe of the clarinet have traditionally been expressed in the frequency domain and were based on measurements or computations of input acoustic impedance. However, numerical simulations of this instrument operate in the time domain and are usually based on the reflection function of the pipe. Recent experimental studies have adopted the time domain approach with direct measurements of this reflection function. Abundant literature extensively covers these subjects: for general presentations, see Refs. 1-4. Studies of the reeds are far less extensive and the mechanical behavior of cane is still subject to discussion. The simplest reed model, a spring, is implicitly used by reedmakers when they rate them by their so-called ''strength,'' which corresponds to the mechanical compliance. Experimental studies have proposed values for the compliance of the reed. [START_REF] Nederveen | Acoustical Aspects of Woodwind Instruments ͑Illinois University Press ͑first[END_REF][START_REF] Gilbert | E ´tude des instruments a `anche simple[END_REF][START_REF] Boutillon | Evaluation of the acoustical stiffness of saxophone reeds under playing conditions by using the reactive power approach[END_REF] Associated with various models of the pipe and excitation, this model has been used in numerical simulations which were successful in describing basic features of the dynamics of clarinet-like system. [START_REF] Schumacher | Ab initio calculations of the oscillations of a clarinet[END_REF][START_REF] Mcintyre | On the oscillations of musical-instruments[END_REF][START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF] Music-oriented algorithms have also been proposed in which the values of the parameters describing the excitor and the resonator are adjusted in order to obtain realistic sounds instead of accurately describing their mechanical behavior. [START_REF] Rodet | Nonlinear dynamics in physical models: Simple feedback-loop systems and properties[END_REF][START_REF] Smith | Physical modeling using digital wave-guides[END_REF] However, this model is obviously insufficient to describe quality-based criteria: otherwise all reeds in a given commercial box ͑with a͒ Part of this work was presented in ''Application of modal analysis and synthesis of reed and pipe to numerical simulations of a clarinet,'' invited paper at the 140th meeting of the ASA, Newport Beach, CA, December 2000 ͓J. Acoust. Soc. Am. 108, 2590͑A͔͒, in ''E ´tude modale d'une clarinette,'' Proceedings of the Colloque National en Calcul de Structures, Giens, France, May 2001, and in ''Modal analysis of a complete clarinet,'' Proceedings of the International Conference on Acoustics, Rome, Italy, September 2001. b͒ Electronic mail: boutillon@lms.polytechnique.fr; present address: Laboratoire de Mecanique des Solides, E ´cole Polytechnique, 91128 Palaiseau Cedex, France. similar strength͒ would suit a given player, but this is not the case.

The next modeling step is the single-degree-of-freedom oscillator. Although some simulation algorithms [START_REF] Gazengel | Time-domain simulation of single reed wind instrument-From the measured input impedance to the synthesis signal-Where are the traps?[END_REF] have been very successful in producing realistic sounds, [START_REF] Ducasse | Modelisation et simulation dans le domaine temporel d'instruments a `vent a anche simple en situation de jeu[END_REF][START_REF] Ducasse | Models of musical-instruments for sound synthesis: Application to woodwind instruments[END_REF] this is not sufficient in itself to assert the physical validity of a model. One degree of freedom is not sufficient to account for criteria such as reed quality. Stewart and Strong [START_REF] Stewart | Functional-model of a simplified clarinet[END_REF] and Sommerfeld and Strong 16 used a refined model of the reed as a nonuniform beam. In the latter study, the pipe was only slightly simplified compared to a real clarinet and the player's air column ͑including the lungs͒ was also taken into account. There is no fundamental difference between this simulation and those based on a simple oscillator model for the reed since the interaction with the rest of the system is averaged along the beam. The beam model is needed if one wants to take into account the curved shape of the mouthpiece on which the reed beats during large amplitude motions. Gazengel [START_REF] Gazengel | Numerical simulations in time and frequency domains-Comparative-study, application to single-reed woodwind instruments[END_REF] derived a simple oscillator model from a beam equation. In his time-domain simulation, the mass of the oscillator is recalculated at each time step as a function of the position of the reed, introducing by this means the nonlinear behavior of the reed contact.

Modeling the reed as a continuous system is the current state of research. Several examples of modal analysis of clarinet reeds with holographic interferometry have been presented in conferences over the recent years, 18 -20 but never published. One example of finite-element modeling based on measurements of the mechanical properties of cane has been reported. [START_REF] Casadonte | The perfect clarinet reed? Vibrational modes of realistic clarinet reeds[END_REF] Another ͑unpublished͒ pioneering study has been done by Pinard and Laine when they were students at the E ´cole Polytechnique ͑France͒. The experimental modal analysis and the finite-element modeling of isolated reeds that are presented in the following are a development of this unpublished work. To our knowledge, no model of the reed as a continuous system in association with the air column has been proposed.

The model proposed here is aimed at overcoming several limitations of previous approaches. Besides giving a means to review the approximations of the classical model, this new approach is also a first step toward numerical simulations of the instrument based on modal projection [START_REF] Facchinetti | Etude des vibrations de l'anche de la clarinette'' and ''Analisi del comportamento dinamico di un clarinetto[END_REF][START_REF] Facchinetti | Application of modal analysis and synthesis of reed and pipe to numerical simulations of a clarinet[END_REF] rather than on propagation schemes represented by reflection functions.

The different parts of a clarinet-reed, mouthpiece, barrel, upper and lower parts of the pipe, bell-are shown in Fig. 1 together with their respective models. Fluid and solid finite-element models ͑FEM͒ for the reed and the beginning of the pipe and a lumped elements model for the main part of the pipe are used.

The work presented here begins with the modal analysis of the isolated reed. In each subsequent section, another element of the model is added, finally resulting in a complete instrument. In addition, the modes of the reed associated with the mouthpiece and barrel are compared with the results of experimental modal analysis.

II. THE REED

A. Construction of the numerical model

Establishing a finite element model requires the determination of the geometry of the reed, the choice of a constitutive law, the determination of the mechanical parameters, as well as the appropriate boundary conditions.

A series of three reeds have been measured. The thickness of each reed was measured with a coordinate measuring machine ͑Mitutoyo EURO-M 574 and Johansson Saphir 7 were used͒. Approximately 200 points, arbitrarily chosen on the reed surface, have been measured ͓Fig. 2͑a͔͒. The geometrical data for the model are interpolated from the measured values. Interpolation between measured points was done by using a fourth-order polynomial, resulting in and giving the thickness map shown in Fig. 2͑b͒. The reed is assumed to be symmetrical with regard to its longitudinal axis.

The shape of the reed was measured using a high precision optical projector ͑Macro Dynascope 5D, by Vision Engineering with Metronics Quadra-Check 4000 interpolating software͒ with the results shown in Fig. 2͑c͒. The precision of the geometrical measurements of the reed can be estimated to Ϸ2 m.

Reeds are made out of cane which is considered here as a purely elastic, transversely isotropic, homogeneous material. Viscosity and plasticity, related to energetic losses, have been neglected at this step of the analysis. The homogeneity hypothesis will be analyzed a posteriori in Sec. V. In the current state of knowledge, we have found no other plausible description that could be expressed quantitatively.

A discussion of losses in cane has been given lately by Marandas et al. [START_REF] Marandas | Mechanical characterization of woodwind reeds[END_REF] and Obataya and Norimoto. [START_REF] Obataya | Acoustic properties of a reed ͑Arundo donax L.͒ used for the vibrating plate of a clarinet[END_REF] The former found out that dry cane is viscoelastic and turns viscoplastic when wet. This implies that static tests on wet cane are not appropriate to measure Young's moduli. Obataya proposed values of the quality factor Q of the order of magnitude of FIG. 1. The clarinet: its parts and their respective models ͑not to scale͒. 100 varying with frequency, relative humidity, and internal state of cane. Since only individual modes of the reed are considered here losses can be ignored. They would need to be taken into account in modeling the actual dynamics of the instrument.

Under these assumptions, five parameters are needed to describe the material: density s , longitudinal and transverse Young's moduli E L and E T , transverse to longitudinal shear modulus G LT , and longitudinal-transverse Poisson ratio

LT . The values adopted here are given in Table I. The values for s , E L , and LT were obtained by Pinard and Laine and result from static measurements on a piece of dry cane given by a reed maker. Obataya and Norimoto give roughly the same value for the main Young's modulus E L of dry cane in the frequency that is relevant here ͑2-6 kHz͒. Their measurements show that this value decreases linearly with the relative humidity level ͑RH͒, E L decreasing by around 30% for a variation of 100% in RH. The other parameters were also obtained by Pinard and Laine. Their work has been pioneering in several respects. In particular, they were the first to match eigenfrequencies and modal patterns of reeds obtained by holographic interferometry with those obtained with a finite-element model.

As boundary condition, we consider the reed rigidly clamped on the section corresponding to the ligature and having a stress-free boundary elsewhere.

B. Computed eigenmodes

This model has been implemented on a standard PC ͑450 MHz, 250 Mbyte RAM, Linux͒ using linear Love-Kirchoff plate elements in the CAST3M finite-element code. The first modes of a reed are presented in Fig. 3. A classification of the modes is needed for referencing and an attempt is made here. Since modal patterns with closed modal lines have not been encountered, it is intuitively appealing to label the modes according to the number of intersections between the nodal lines and the edges of the reed. For the symmetric reed considered here, a mode is labeled LnTm. ''L'' stands for longitudinal and the first index n is the number of intersections of nodal lines with the edge͑s͒ parallel to the main axis. Such nodal lines include the one imposed by the boundary condition at the ligature. ''T'' stands for transverse and the index m is the number of intersections of the nodal lines with the tip edge of the reed. Modes appear in an order which can be expected (L1T0, L1T1, L2T0, L1T2, L2T1), given the larger flexibility in the direction transverse to the reed and the thickness distribution.

The generalized mass of a mode is:

mϭu T M s u, ͑1͒
where u represents the reed displacement for the mode and M is the mass matrix of the reed. For a unit value of the maximum displacement in each mode, the modal masses are 7, 0.35, 0.47, 0.063, and 0.094 mg for the L1T0, L1T1, L2T0, L1T2, and L2T1 modes, respectively. Along with modal patterns, these values establish a comparison between modes. These mass values can also be compared to the order of magnitude of the real mass of the moving reed. At the tip of the reed, the thickness is about 1/10 mm and the width 13 mm. For a density s ϭ450 kg m Ϫ3 , the mass of a moving portion of the reed of length l ͑in mm͒ is (0.59ϫl) mg.

III. MODAL COMPUTATION OF THE REED ASSOCIATED WITH MOUTHPIECE AND BARREL

This section analyzes how the dynamics of the reed is influenced by air loading and provides a comparison between results given by the model and experiments presented in Sec. II. The system considered now is composed of the reed, the mouthpiece, and the barrel and is represented using a coupled fluid-solid model.

A. Numerical model

The full model of reed, mouthpiece, and open barrel is shown in Fig. 4. The internal shape of the mouthpiece ͑a Selmer HS*) has been carefully measured by means of the coordinate measuring machine used for the reed. The barrel is considered as a cylindrical bore with a diameter of 15 mm. The air volume inside the mouthpiece and the barrel is modeled with linear tetrahaedric and prismatic finite elements of compressible elastic fluid.

The acoustic pressure at points of the open air surfaces is considered to be zero. The normal derivative of the acoustic pressure on the walls of the mouthpiece and the barrel, corresponding to air flow, is also set to zero. The boundary condition coupling the reed and the mouthpiece involves the stress in the solid and the velocity of the fluid and will be given explicitly in the following.

The eigenvalue problem for a coupled solid-fluid system is expressed in the continuous formulation by the following: [START_REF] Gibert | Vibrations des Structures-Interactions avec les Fluides ͑Eyrolles[END_REF] div CٌuϪ 2 s uϭ0, ͑2͒

div 1 f "pϩ 2 1 c 2 f pϭ0, ͑3͒
where p represents the acoustic pressure in the fluid. The densities of solid and fluid are s and f , respectively. The speed of sound is c, the angular frequency of the motion is , and C denotes the elasticity matrix of the solid. The boundary conditions coupling the fluid and the solid parts are "nϭϪpn, ͑4͒ ץp ץn ϭ f 2 u"n, ͑5͒

where n represents the unit vector normal to the solid surface and ϭCٌu the stress tensor.

In order to formulate these equations as a standard eigenvalue problem, a new variable ϭϪ(1/ 2 )p must be introduced. [START_REF] Gibert | Vibrations des Structures-Interactions avec les Fluides ͑Eyrolles[END_REF] The equations and boundary conditions become div CٌuϪ 2 s uϭ0, ͑6͒

div 1 f "Ϫ 1 c 2 f pϭ0, ͑7͒ 2 ϩpϭ0, ͑8͒ "nϭϪpn, ͑9͒ ץ ץn ϭϪ f u"n. ͑10͒
To the preceding equations, we can associate the following Lagrangian L denoting the variational formulation of the problem: Lϭ

1 2 ͵ ⍀ s ٌuCٌu dvϩ 1 2 ͵ ⍀ f 1 f c 2 p 2 dv Ϫ 2 ͩ 1 2 ͵ ⍀ s s u 2 dvϪ 1 2 ͵ ⍀ f 1 f c 2 ͑ " ͒ 2 dv Ϫ ͵ ⍀ f 1 f c 2 p dvϪ ͵ ץ⍀ "n ds ͪ ,

͑11͒

where ⍀ s and ⍀ f represent the solid and fluid volumes, respectively, and ץ⍀ represents the boundary between these volumes.

Finally, the problem is expressed in its discrete form by the following eigenvalue problem:

ͩͫ K s 0 0 0 K f 0 0 0 0 ͬ Ϫ 2 ͫ M s 0 ϪN 0 0 K f ϪN T K f T ϪM f ͬͪ ͫ u p ͬ ϭ ͫ 0 0 0 ͬ ,
where K s ͑respectively, K f ) and M s ͑respectively, M f ) are rigidity and mass matrices of the solid ͑respectively, fluid͒ part of the system and N is the operator corresponding to the coupling boundary condition ͑10͒ related to the normal vector n. Details of the derivation can be found in Ref. 26.

B. Experimental modal analysis

An experimental modal analysis on reeds by means of holographic interferometry was performed in order to check the validity of the numerical model of the reed coupled to air. Recent works have been reported in short communications. 18 -20 For various reeds mounted on a mouthpiece under dry conditions Pinard and Laine observed one mode corresponding to a longitudinal flexion at around 2200 Hz; one family of modes around 3500-3700 Hz, with patterns varying from reed to reed, some of them being in-dicative of torsion, others being closer to flexion; and one family of modes around 5800-6000 Hz, with more complex patterns.

In measurements presented here, the reed was attached to the mouthpiece exactly as on the real instrument. Since the ligature was producing strong light reflexions, it was replaced with adhesive tape placed slightly further from the tip. A sinusoidally driven loudspeaker was placed close to the reed to excite its vibration. For determining the resonance frequencies a very thin PVDF piezoelectric film ͓poly͑vinylidenefluoride͒, thickness 0.05 mm, mass 0.06 g, of which only a part was actually moving͔ was glued onto the lower thicker part of the reed, yielding the average deformation near the ligature. Resonance frequencies were determined using the maximum of the piezoelectric signal. The experiments were performed under natural humidity. A saturated atmosphere would have been preferable but was not possible with the interferometry equipment.

The eigenmodes were visualized by means of laser transmission interferometry. Complete details of the implementation of this classical method are described in Ref. 27. The images in Fig. 5 represent variations of equal normaldisplacement of the reed. The resolution of the system is half the wavelength of the laser, approximately 0.3 m.

The reed was measured either alone, associated with an open mouthpiece, or with the mouthpiece mounted on an open barrel. The first four measured modes shown in Fig. 5 correspond to the barrel configuration ͑see Fig. 4͒. They are compared with the corresponding computed modal patterns ͑see the next section for computation of the eigenmodes͒ on top. Results of the holographic measurements show that the maximum displacement of the reed is negligible compared to the distance between the mouthpiece and the reed at that level of excitation. Thus one can be confident that contact between the reed and the lay, which could possibly make the system nonlinear, does not occur. FIG. 5. Projection of four eigenmodes on the reed ͑see the text for labels͒. Top pictures: computed normalized eigenmodes of the association of a reed with mouthpiece and barrel. In this representation, a cyclic gray scale produces fringes of equal differences in normal displacement, allowing a comparison with the modal patterns observed experimentally. Bottom pictures: modal patterns measured by holographic interferometry on one good reed mounted on the mouthpiece attached to the barrel. The resonance is not very sharp owing to damping, hence the rounded eigenfrequencies. The photographed section of the reed does not have the same height between the various experiments and the simulations.

C. Results

A comparison between computed and measured modes is displayed in Fig. 5 for the situation described by Fig. 4. In this comparison with holographic measurements, the ''ligature'' of the reed had to be placed slightly beyond its normal position. This led to a slightly more flexible reed than in the normal situation. When the reed is coupled to air, one should stress that eigenmodes concern the whole system, not just the reed. Strictly speaking, expressions such as ''reed modes'' are inappropriate and refer instead to modes for which energy is mostly localized in the reed. Each mode has been labeled using the notation proposed earlier. The ''R'' prefix indicates that we regard the result just as the projection of the four first eigenmodes on the reed subspace. In order to simplify the discussion, we have not attempted to label the air configuration. One can notice that the L2T2 pattern did not appear in the isolated reed case. One can also notice that the L1T0 pattern of the reed appears in the two first modes of the coupled system.

The computed modes appear in the same order as the measured ones with eigenfrequencies deviating by 10%-20% from measured resonance frequencies. The modal patterns are globally the same despite the fact that no real reed is symmetric whereas the numerical model has been chosen symmetric. As expected, the modes are mainly localized at the tip of the reed where it becomes very thin, showing the importance of a precise measurement of the geometry. Although some of the mechanical parameters come themselves from a fit between observation and computation of modes of an isolated reed, the mixed fluid-solid model can be considered as valid within the range of approximations retained here.

Real reeds have natural asymmetries due to their geometry or to nonuniform mechanical properties. One notices that the asymmetry seems stronger for the lowest mode than for any other one.

D. Sensitivity analysis

The sensitivity analysis of the eigenfrequencies to variations in mechanical parameters describing the reed and in acoustical properties of the air is presented in Table II. The air volume is that of Fig. 4. Parameters are varied by 5% above and below their average values ͑i.e., 10% overall͒ and the corresponding overall variations of eigenfrequencies are reported. The value of the Poisson ratio appears to be irrel-evant. Eigenfrequencies 1190, 2680, and 4010 Hz vary linearly with the speed of sound. This is also almost the case for the mode at 5280 Hz. Without looking at the modal pattern of air pressure or reed displacement, one can infer that they correspond to '' 'air modes,'' with energy mostly localized in the ͑short͒ pipe. Conversely, the mode at 3700 Hz is not influenced by air characteristics; sensitivity to the shear modulus G LT indicates that the reed is subject to torsion ͑see the second mode of Fig. 3͒ and is poorly coupled to the pipe ͑Fig. 6͒. To a lesser degree, this is also the case of the mode at 6300 Hz. The mode at 4740 reveals a (E L / s ) 1/2 dependency of the eigenfrequency. It is mostly a ''reed mode'' involving primarily a longitudinal deformation. The mode at 2010 Hz is apparently a mode in which air and reed are strongly coupled. It is interesting to notice that the transverse Young's modulus does not seem to influence any frequency. The measurement of its precise value is therefore less particularly important.

E. Evolution of the eigenfrequencies

Another way of examining how the reed is coupled to the acoustic field is to follow the evolution of the eigenfrequencies when the reed is loaded by the air volume of mouthpiece and barrel. A decrease of the eigenfrequencies and a dominance of the longitudinal flexion occurs in the eigenmodes ͑Fig. 7͒.

The frequencies of the first two modes of the ͕reed, mouthpiece, barrel͖ system are mainly imposed by the resonance of the air cavity. In both modes, the reed undergoes mainly longitudinal flexion. The frequency of the torsion mode L1T1 ͑3257 Hz for the isolated reed͒ does not vary FIG. 6. Computed eigenmode at 4119 Hz in a mixed solid-air situation: acoustic pressure inside the mouthpiece and barrel. significantly, meaning that this mode is weakly coupled to the air cavity. The same phenomenon can be noticed for the mode L1T2 at 5840 Hz for the isolated reed. One can conclude from Fig. 7 and from the observation of the rather uniform pressure in the pipe at these modes ͑not shown here͒ that this mode also is weakly coupled to the pipe.

IV. MODAL COMPUTATION OF THE WHOLE CLARINET

In order to simulate the modal behavior of the complete clarinet, we have associated a finite-element model of Ϸ10 cm of pipe with lumped elements representing the rest of the pipe and matching its acoustic input impedance. This can be done since at the outlet of the barrel, the acoustic field consists essentially of plane waves. An example of acoustic pressure in the mouthpiece is represented in Fig. 8. The mode is that of a complete clarinet and corresponds to the lowest mode at 311 Hz of the medium C fingering combined with the opening of the register key ͑see the following for the complete list of modes in this configuration͒. The length of mouthpiece represented here is 32 mm and corresponds to the tapered part. One can see that the acoustic waves can already be considered as plane waves within a very good approximtation.

The lumped-element oscillators ͑shown in generic form in Fig. 1͒ are coupled to the finite-element barrel by means of a rigid plate of negligible mass, as shown in Fig. 9. The plate and the lumped-element oscillators are supposed to move only in the longitudinal axis of the instrument. The lumped elements are placed at the ͑virtual͒ junction between the barrel and the lower part of the clarinet.

It is now explained how the numerical values of the lumped elements are calculated on the basis of measurements provided by Gibiat [START_REF] Gibiat | Acoustical impedance measurements by the 2-microphone-3-calibration ͑tmtc͒ method[END_REF] on several notes of a Noblet B clarinet. Results of these measurements are supposed to represent the input acoustical impedance of the instrument. In order to measure this input impedance, a reference plane was defined by Gibiat et al. by replacing the mouthpiece with a portion of cylindrical tube of equal volume. This is the usual ''equivalent volume'' approximation which we discuss later on. Prior to matching the impedance of the lumped elements to the measured input acoustical impedance of the pipe, the latter must therefore be transported from the input plane toward the open end of the pipe. The ''transportation distance'' is equal to the length of a cylinder having the volume of the mouthpiece and the barrel. [START_REF] Facchinetti | Etude des vibrations de l'anche de la clarinette'' and ''Analisi del comportamento dinamico di un clarinetto[END_REF][START_REF] Gibiat | Mesures d'impe ´dance acoustique pour la clarinette[END_REF] An oscillator is associated with each measured impedance peak. At the angular frequency the mechanical impedance of each elementary oscillator in Fig. 1 is

Z m ͑ ͒ϭi ͩ 1 m Ϫ kϩir ͪ Ϫ1 , ͑12͒
where m, r, k are respectively the mass, damping, and stiffness of the lumped elements.

In this ''comb-like'' association, the impedances of the oscillators add. The dual association where the admittances add is ''chain-like.'' Each elementary oscillator of Fig. 1 is a mass chained with a comb of a damper and a spring, leading to Eq. ͑12͒.

The parameters m i , r i , k i of each oscillator ͑a tooth of the large comb͒ are identified by minimizing a cost functional J measuring the distance between computed and measured moduli and phase of the impedance:

Jϭ␣͉Mod͑Z comp ͒ϪMod͑ Z meas ͉͒ ϩ␤͉Arg͑Z comp ͒ϪArg͑ Z meas ͉͒. ͑13͒
The initial values of the parameters for each oscillator are obtained by identifying each single resonance peak and the final values are obtained by running a Nelder-Mead simplex search algorithm. A comparison between the measured and the identified modulus and phase of the acoustic impedance of the lowest F fingering ͑E heard͒ of the clarinet is presented in Fig. 10. The impedance represented is not the input acoustical impedance but the impedance of the lower part taken at the ͑virtual͒ junction between the barrel and the lower part of the clarinet. Therefore, the peak frequencies are not the eigenfrequencies of the instrument. The acoustical impedance represented here is the ratio of the acoustical pressure to the air velocity, normalized by c. The average modulus on a logarithmic scale would be 1 for an ideal long cylindrical pipe. According to Gibiat, it is less here due to internal losses, radiation, and presumably the complexity of the pipe. FIG. 7. Evolution of the eigenfrequencies ͑left scale, in hertz͒ when the system evolves from the isolated reed ͑left͒ to ͕reedϩmouthpiece͖ ͑middle͒ and ͕reedϩmouthpieceϩbarrel͖ ͑right͒. Black lines represent ''primary reed'' modes, dotted lines ''primary air'' modes, and dash-dot lines, ''mixed'' modes. FIG. 8. Acoustic pressure inside the tip part of the mouthpiece for a 311 Hz mode of the complete clarinet. The acoustic pressure decreases monotonically from the tip to the largest section by 14%.

Two eigenmodes of the complete instrument for different fingerings are shown in Fig. 9. One eigenmode has no amplitude per se. For each eigenmode in Fig. 9 the ͑relative͒ amplitude of the motion of the oscillators is represented by the length of a straight line extending from the plate. One notices that the pressure distribution is not uniform in the mouthpiece. Examining other similar figures reveals that the motion of the reed can differ significantly from mode to mode of a given note, even if it follows a L1T1 pattern. This means that, although the first modes of the isolated reed occur at significantly higher frequencies than those considered here, a single degree of freedom for the reed is not appropriate since it would not account for these differences. When the reed undergoes mostly longitudinal flexion, it is to be expected that the beam model used by several authors [START_REF] Stewart | Functional-model of a simplified clarinet[END_REF][START_REF] Sommerfeldt | Simulation of a player clarinet system[END_REF][START_REF] Gazengel | Caracte ´risation ... des instruments a `anche simple[END_REF] would give comparable results.

For the low F fingering ͑sounding one tone lower͒, the first eigenfrequencies are 166, 464, 743, 1147, 1436, 1620, 1950, 2058, and 2201 Hz. They are 373, 1035, 1541, 1687, 1893, 1930, and 2309 Hz for the medium G fingering and 311, 735, 1213, 1467, 1578, 1865, and 2211 Hz for the high G , played with medium C fingering and opening of the register key. These frequencies are represented in Fig. 11 in order to evaluate their harmonicity. Eigenfrequencies are normalized by their ratio to the theoretical musical frequency for the note under consideration ͑respectively, 156, 349, and 740 Hz͒, rounded to the nearest integer. For example, a 900 Hz eigenfrequency for note A4 ͑440 Hz͒ would be normalized by 2, nearest integer to 900/440. For this high note, the register key does not eliminate the first mode of the instru-ment but the sound will be locked approximately on the second mode. The lowest mode is very roughly at half the pitch of the note and is therefore normalized by the integer 2.

The sets of solid lines in Fig. 11 ment where the volume of the mouthpiece has been replaced by a cylindrical pipe having the same volume and closed at one end. The sets of dotted lines represent computed eigenfrequencies of the air column with a rigid boundary on the reed surface. Instead of the completely closed pipe of the traditional modeling, one assumes here a slight opening between the reed surface and the lay of the mouthpiece with a zero pressure condition.

V. DISCUSSION AND PERSPECTIVES

A. Alignment of resonances and low-frequency approximation

The traditional model of the mouthpiece is that of a cylinder of equivalent volume. Within this approximation there is no point in measuring the input acoustic impedance above a certain limit. This limit can be evaluated by the length scale at which the mouthpiece geometry differs from a cylinder. Taking as an order of magnitude for these geometrical differences a length of 1 cm is consistent with a 2.5 kHz frequency limit beyond which input acoustical impedances would begin to differ noticeably. In the approach followed in this paper, the equivalent volume approximation is abandoned and the acoustical input impedance of the pipe would keep full utility and validity up to the frequency of the first transverse mode of the pipe ͑13.3 kHz for the clarinet͒.

The cylinder of equivalent volume approximation for the mouthpiece is assumed to be correct for low frequencies. It appears in Fig. 11 that this approximation is not acceptable enough to be used in conjunction with an alignment of peaks criteria. One can see in Fig. 11 that variations in eigenfrequencies due to the model change are significant with regard to the alignment of resonances, even at low frequencies. In other words, the deviations from alignment in the traditional view ͑equivalent volume approximation͒ are of the same order of magnitude as the frequency shifts due to the presence of the reed and the prismatic shape of the mouthpiece.

B. Coupling of torsion modes to the air

The association of reed, mouthpiece, and a short open portion of the pipe is shown in Fig. 6. The modal acoustic pressure at an eigenfrequency of 4119 Hz is displayed in Fig. 6. In this mode, the reed undergoes torsion in a pattern very similar to the L1T1 mode of the isolated reed ͑Fig. 3͒. The characteristic distance of this modal deformation is significantly smaller than half the wavelength in air at that frequency (Ϸ10 cm); the resulting acoustical short-circuit prevents any efficient coupling of the reed to the air in the mouthpiece. This explains the fairly uniform acoustic pressure for this mode, except very near to the reed. However, there are several reasons why these modes may be important in the actual playing.

First of all, the flow entering the air channel between the reed and the lay is governed by a nonlinear equation. Therefore, antisymmetric reed modes may have an influence on the global flow entering the pipe.

It has been shown that the antisymmetric reed modes are very weakly coupled to the acoustic ͑far͒field in the clarinet. This is not to say that these modes play no role in the dynamics. Asymmetries or, better said, unevenness in the geometric or constitutive properties of reeds induce asymmetries of longitudinal reed modes and consequently an asymmetry in the local acoustical field. Due to its small relative modal mass, the torsion mode can be easily excited at a frequency different from its resonant frequency and therefore may play a significant role in the actual dynamics of the reed. The coupling factor would then be the local acoustic field. This may be an explanation for the player's experience that for different mouthpieces, the preferred reeds are also different.

This modal analysis is performed on a symmetric reed. This is not the case in reality as shown for example by the first mode in Fig. 5. The so-called torsion modes are likely to be associated in the fluid domain to a flow different from zero and therefore couple to the plane waves inside the pipe.

C. Symmetry

Experimental modal analysis shows that some reeds have strong asymmetries. Makers can be expected to be successful in controlling the symmetry of the geometry; therefore, the cause of modal asymmetries lies most probably in the lack of homogeneity of the cane used for the reed due to its natural character. Pinard, Laine, and Vach 31 examined 24 reeds, ranked by two professional players. They observed that the two reeds ranked as good and very good were symmetric whereas the poor reed had asymmetrical high modes. Based on limited sampling of reeds and players, no definite conclusion can be drawn. Intuition would suggest that asymmetry is not a desirable feature for a reed. However, we think that it might not be so.

Visualizing the lip motion in brass playing shows that lips do not move symmetrically and that this factor varies from player to player. Since brass mouthpieces are symmetric, one can conclude that the mechanical properties of lips ͑possibly coupled to dentition and the mouth cavity͒ are not symmetric for all brass players. One can hypothesize that the same is true among clarinet and saxophone players. Another observation is that different players do not always prefer the same reeds in a given box, even for common musical tasks, style, etc., and the same clarinet and mouthpiece. A good match between a player and a reed could mean that a given asymmetry in a reed would fit well the natural asymmetry of a given player and not so well with another one. It has even been observed that a few players use reeds which fit almost none of their colleagues. It would be interesting to test these players and their preferred reeds with regard to the symmetry hypothesis.

VI. CONCLUSION

The modal analyses of reeds and of a few notes of the whole clarinet were performed. Results have shown the following points.

͑1͒ A numerical model of cane based on the hypothesis of transverse isotropy is suited to describe modal patterns of reeds. Some of the numerical hypotheses ͑homogeneity, symmetry, damping͒ can be released but this would necessitate additional measurements.

͑2͒ When coupled to air, the reed is subject to deformation patterns which are not always those of its own normal modes. Therefore, the normal modes of isolated reeds cannot be taken as a source for the acoustic field in the mouthpiece. Specifically, coupling must be taken into account.

͑3͒ Torsion modes of reeds generate a strong but very localized acoustic field in the mouthpiece. It remains to be examined how this would interact with asymmetries in lower modes through the excitation process.

͑4͒ Acoustic waves are already plane within a very good approximation in the cylindrical part of the mouthpiece. Since finite-element modeling of air is interesting insofar as the waves are not plane, the air volume in the barrel and a large proportion of that in the mouthpiece can be included in the lumped-element model, reducing significantly the computational burden.

͑5͒ The shape of the mouthpiece and the dynamics of the reed influence the alignment of resonances in the same proportion as the misalignment derived from the customary observation of the input acoustical impedance. Therefore, the approximation of the equivalent volume is too coarse to be used when looking at harmonicity of resonances.

This study shows the need for input impedance measurements at higher frequencies than usually performed. It calls for simplified formulations of the acoustic field in the mouthpiece. The procedure outlined here could be used to test these formulations. Finally, the method paves the way for numerical simulations of the dynamics of the clarinet based on modal projection and taking into account the whole complexity of the reed.

FIG. 2 .

 2 FIG.2. Geometry of the reed, with dimensions in mm: ͑a͒ points actually measured, ͑b͒ interpolated thickness, ͑c͒ estimated contour.

FIG. 3 .

 3 FIG.3. First five computed modes of an isolated reed. Modes are labeled according to the number of modal lines perpendicular to the main axis ͑Ln͒ and parallel to it ͑Tm͒.

  FIG. 9. Modal representation of a complete clarinet: amplitude of the motion of the lumped-element oscillators ͑left͒, air pressure in the upper part of the pipe ͑middle͒, and deformation of the reed ͑right͒. Eigenmode 2 for note treeble F# ͑fingering of C# medium plus opening of register key͒ and eigenmode 8 for note low E ͑low F fingering͒.

FIG. 11 .

 11 FIG. 11. Normalized eigenfrequencies ͑logarithmic scale͒ of the complete clarinet, pipe with reed ͑solid symbols͒, of the pipe with a fixed reed ͑dash dot͒, and normalized resonance frequencies measured on the pipe where the mouthpiece replaced by its equivalent volume ͑dashed͒. See the text for the definition of the normalization. Fingerings are low F, medium G, and high G ͑medium C with register key͒ corresponding to notes E 3, F 4, and G 5.

TABLE I .

 I Material properties for dry cane used in reeds, as given by Pinard and Laine.

	Density	s ϭ450 kg/m 3
	Longitudinal Young modulus	E L ϭ10 000 MPa
	Transverse Young modulus	E T ϭ400 MPa
	Shear modulus	G LT ϭ1300 MPa
	Poisson ratio

LT ϭ0.

22

TABLE II .

 II Sensitivity analysis: changes in eigenfrequencies when mechanical characteristics of the reed and acoustical properties of the air vary. Changes are given in % for a 10% variation of each parameter.

	⌬ϭ10% mean values	E L 10 4 MPa	E T 400 MPa	G LT 1300 MPa	LT 0.22	s 450 kg m Ϫ3	c 340 m s Ϫ1	f 1.23 kg m Ϫ3
	1190 Hz	0	0	0	0	0	9.8	0
	2010 Hz	2.4	0	0	0	Ϫ2.2	0.7	Ϫ0.4
	2680 Hz	0.1	0	0	0	Ϫ0.2	9.6	0
	3700 Hz	1.5	0	3.1	0	Ϫ4.6	0	0
	4010 Hz	0.2	0	0	0	Ϫ0.2	9.3	0
	4740 Hz	4.9	0	0	0	Ϫ4.8	2.7	Ϫ0.1
	5280 Hz	0.6	0	0	0	Ϫ0.9	8.1	Ϫ0.1
	6300 Hz	1.7	0.9	4.9	0	Ϫ6.4	3.2	0

ACKNOWLEDGMENTS

We express our gratitude to Holger Vach for his decisive help in the experimental part of this study, to Vincent Gibiat for providing us with the measurements of the acoustic input impedances and the associated software, and to Brian Katz for many language corrections.


G Bertolino
email: bertolin@cab.cnea.gov.ar

J Perez Ipin ˜a

G Meyer

Influence of the crack-tip hydride concentration on the fracture toughness of Zircaloy-4

The influence of a hydrogen concentration gradient at the crack-tip and hydride platelet orientation on the fracture toughness, fracture mode and micromechanisms of a Zircaloy-4 commercial alloy was studied. Fracture toughness was measured on CT specimens and the analysis was performed in terms of J-integral resistance curves at temperatures ranging from 293 to 473 K. Fracture toughness results of specimens containing higher hydrides concentration near the crack-tip region, preferentially orientated in the crack plane, were compared to those obtained from specimens with a homogeneous hydrogen distribution and different platelet orientation; specimens were obtained by charging them in loaded and unloaded condition, respectively. Changes on both macroscopic and microscopic fracture behaviour were observed at temperatures ranging from 293 to 343 K, and the results show the relevance of both hydride concentration and platelet orientation. The existence of a ductile-to-brittle transition is discussed at the light of these new results.

Introduction

The combination of radiation, oxidation and hydriding degrades the mechanical behaviour of fuel cladding during nuclear reactor operation. Zirconium alloys are the sole material employed for cladding tubes of nuclear fuel because of their neutron transparency. As zirconium is a hydride forming metal, the hydrogen absorption leads to the formation of different hydrides, with the subsequent degradation of mechanical properties. Inservice, as the total amount of hydrogen in a cladding increases, hydrides precipitate preferentially close to the outer surface of the cladding, which is the coolest one.

The factors reported to affect hydride platelet orientation in zirconium alloy are crystallographic texture, prior strain and applied stress [1]. Hydride precipitation under stress shows preference for the orientation of hydrides, where platelets precipitate parallel to the planes of highest principal stresses 1 [1][2][3][4], this phenomenon being called stress-reorientation.

The influence of an uniform distribution of hydrides on the mechanical behaviour of zirconium-based alloys has been extensively studied, usually through uniaxial tensile tests [5][6][7][8]. The mode I crack growth behaviour of Zircaloy-4 has been also examined in a previous work [9]. The occurrence of a ductile-to-brittle transition with temperature and/or H content was reported in zirconium alloys with homogeneous hydrogen distribution, also under uniaxial tensile loading. This transition seems to be strongly dependent of the material microstructure, thermal history, hydrogen content, morphology and orientation of hydride precipitates [10]. However, fracture toughness testing did not reveal such transition for hydrogen levels under 2000 ppm at temperatures between room temperature and 470 K [9].

When the above-mentioned results are intended to be transferred to industrial service conditions, a question immediately arises: is the homogeneous hydrogen distribution representative of the actual distribution present in a material submitted to service conditions, especially in cases where cracks could be present? When a crack is present in the in-service component, a non-homogeneous hydride precipitates distribution can be expected, with a larger presence of H in the crack-tip region, leading to a different behaviour of the crack and, very important, to a possible non-conservative condition when results from tests with homogeneous H distribution are employed.

The influence of a non-homogeneous distribution of hydrides has been examined in two recent works, when hydrides are present in the form of a layer or rim in cladding tubes [START_REF] Daum | 13th ASTM International Symposium on Zr in the Nuclear Industry, ASTM STP 1423[END_REF], and also in the case of blisters in sheets [START_REF] Pierron | [END_REF], showing a significant loss of ductility with increasing hydrogen content and hydride rim thickness, although in the case of blister the fracture process was controlled apparently by the crack growth resistance of the substrate.

The aim of this study was to examine the influence of a concentration gradient of hydride precipi-tates in the near crack-tip region on the fracture behaviour of unirradiated Zircaloy-4 alloy under mode I loading. The objective of this hydrogen distribution has been to obtain specimens with a hydride distribution representative of a loaded component containing a crack or a defect. The study was carried out, by testing CT specimens, on recrystallized Zircaloy-4 with a crystallographic texture similar to that of as-fabricated Zircaloy-4 cladding tubes. Some aspects related to the ductile-to-brittle transition in this material are discussed.

Materials and methods

Material and specimens

The material of this study was a Zircaloy-4 in plate shape, hot and cold rolled to a final thickness of 6 mm and recrystallized leading to an equiaxed grain structure, with an average grain size of 8 lm. A summary of the principal alloying elements is given in Table 1. As a result of the previous plastic deformation processing, the sheet exhibited a strong texture with its basal pole tilted around ±30°away from the normal direction toward the transverse direction [9].

Standard fatigue pre-cracked compact tension (CT) specimens with dimensions in accordance with ASTM E 1820-99 [START_REF] Astm | Standard Test Method for Measurement of Fracture Toughness[END_REF] were employed to the characterization (W = 25.4 mm), Fig. 1. The specimens were machined with the same thickness of the plate (6 mm), implying that the rolling plane was normal to the crack front, and the notch parallel to the rolling direction, Fig. 1.

Hydrogen charging method

The hydriding process was performed on specimens (pre) mechanically loaded in mode I. The mechanical loading process consisted of loading the specimen in a test machine until a desired K (stress intensity factor), then by using a simple wedge the imposed displacement remains after the load is removed. At this point some load relaxation occurs, so the effective K is not well known. All the specimens were loaded at similar K values, around 18 Mpa p m. In this way, a stress distribution characterized by high stress gradient, Westergaard-Irwin type, was introduced in the specimen, making more hydrogen diffuse into the areas of highest stress levels, leading to a high hydrogen concentration gradient in the near crack-tip region. Two consecutive techniques were employed for the hydrogen charging; an electrochemical technique first, and then an annealing in a hydrogen rich atmosphere. Electrochemical charge was carried out at room temperature in a solution of 0.1 M KOH with a current density of 10 mA/cm 2 during a period of time determined by the desired final concentration. The annealing process was performed at 450 °C at hydrogen pressure between 1.5 and 3 bar. Then, the specimens were slowly cooled inside the furnace, and keeping the mechanical load. The specimens were unloaded after the cooling process was complete.

Fracture tests

ASTM E 1820-99 standard was followed, evaluating the J-R curves by means of the single-specimen unloading compliance technique. The stable crack growth initiation value (J IC ) was calculated when possible, while J C was evaluated for the specimens exhibiting brittle fracture. J C corresponds to the J value at the fracture point. After the test, those specimens that exhibited only ductile behaviour were fatigue post-cracked to delineate the stable crack growth. Because the compliances were not measured on the load line, data corrections had to be made [14]. The tests were performed at temperatures between 293 and 473 K. Specimens were heated by using two strip heaters and temperature measurement and control were performed by means of a K-type thermocouple placed close to the crack-tip. The temperature was controlled within ±1 K.

Some J initiation values, mostly for tests at high temperature, were higher than the maximum valid J (133 kJ/m 2 ) in accordance with standards for this material and with specimen thickness of 6 mm. In these cases the J initiation value is not J IC but J Q in accordance with the standard.

Metallographic and fractographic analysis

A standard procedure of mechanical polishing and chemical etching was carried out for the metallographies, with a solution of 45 pct glycerol, 45 pct HNO 3 and 10 pct HF to reveal the microstructure and then a solution of 92 pct HNO 3 and 8 pct HF to reveal the Zr-hydrides distribution.

After the tests, one-half of each specimen was observed through scanning electron microscopy (SEM) to study the fracture surface and hydride distribution; the remaining halves were prepared for hydrogen content measurements by inert gas fusion technique using a LECO TM hydrogen analyser. The accuracy of the hydrogen content measurement is ±10 ppm. The H content was measured at least at five different locations, four of them along the remaining ligament as indicated in Fig. 2.

Results and discussion

The employed hydrogen charge produced hydride precipitates of d phase type, platelet shaped, with a concentration along the crack path exponentially diminishing from its highest value at the crack-tip -Fig. 2 shows this behaviour for three specimens. This was the expected result due to the variation of hydrogen solubility with stress. Fig. 3 shows two micrographs of a specimen with 1634 H ppm at the crack-tip and 173 H ppm in the bulk (face LT), the variation on the hydride concentration as a function of the distance to the crack-tip can be clearly observed. Hydrides generally precipitated at grain boundaries, following their curvature; however, the formation of a continuous network of hydride along the grain boundaries was often observed for regions with medium hydrogen content, i.e. 173 ppm, Fig. 3(b). Transgranular precipitation was also found in the near crack-tip region, which corresponded to the highest concentration zone. Micrographs showed that the plane orientation in which platelets precipitated tended to be parallel (in average) to the LS plane and normal to the T direction, at least in the regions near the crack path, Fig. 3. Instead, in specimens charged without pre-stress, i.e. homogeneous distribution of hydrogen, the platelet orientation was parallel to the LT plane, Figs. 3(c) and 4. It had several times been reported that hydride precipitation under stress shows preferential orientation, where the normal to the platelet is parallel to the tensile stress direction [2]. This phenomenon is related to preferred nucleation and growth orientation and has been reported to occur in pressure tubes where the precipitation direction changes from circumferential to radial under certain specific stress state [1]. The platelets orientation in the present specimens was consistent with this process, Fig. 4. According to Hardie [15], the hydride orientation due to stress is extremely dependent of the temperature and the stress level, with a threshold stress (r th ) depending on the annealing temperature and the material yield strength. Bell and Duncan [16], suggested, in a work on experimental determination of r th , that r th should reduce to a null value at a temperature of 704 K. Regarding our experiments, all the tested specimens with homogeneous hydrogen distribution throughout the volume were heat treated at 723 K and no change in platelet orientation was observed; only those specimens with a mechanical load imposed during the hydrogen charge showed platelet preferentially orientated to the crack plane. These results are not in accordance with BellÕs statement, but with SinghÕs who found no preferred orientation for tests carried out at 723 K with no load.

All specimens containing a non-homogeneous H distribution and tested at room temperature failed by brittle mode without previous stable crack growth in a range of J C values from 15 to 25 kJ/ m 2 (corresponding to K values from 38 to 50 MPa p m). Fig. 5 shows the corresponding J C values plotted as a function of the H content far away from the crack-tip, i.e. in the region of the bulk hydrogen content; and, for comparison, the J IC values obtained by testing similar specimens with homogeneous hydrogen distribution in the specimen volume [9]. Comparing these results, it seems that, instead of the bulk H level, the hydride distribution and the platelet orientation near the crack-tip control the fracture behaviour, rendering it brittle due to a high local hydride concentration even when the bulk H level is low, or the high concentration zone is small. Results from SEM in situ test [17] at room temperature showed that when hydrides were oriented parallel to the crack propagation plane, the crack grows following a path of cracked hydrides either without or with minimum matrix plastic deformation.

The observation of the fracture surfaces of specimens charged under stress showed that the dominant fracture micromechanism at room temperature was cleavage, Fig. 6; in contrast to the one observed for the homogeneous H distribution where all specimens fractured by nucleation and coalescence of microvoids. As can be observed in Fig. 7, some specimens showed a fracture Ôtransi-tionÕ region between the fatigue pre-crack zone and that of fast crack propagation during the test, where microcracks normal to the crack growth direction are present. This area could be interpreted as produced by a sub-critical crack growth mechanism developed during the hydrogen charging stage. The applied pre-stress in these specimens The fracture behaviour for the tests performed at 343 K was also brittle, as it was at room temperature, although the records showed a clear non-linearity, and some stable crack growth was observed before unstable fracture. At higher temperatures, instead, the influence of hydrogen diminished and ductile behaviour occurred. When results are plotted as a function of temperature, Fig. 9, a fracture mode transition can be observed between 343 and 393 K. Analysis of data as a function of hydrogen content is more difficult because at high temperatures no tested specimen kept a high bulk H level. The observation of the fracture surfaces showed that the fracture micromechanism at high temperature was nucleation and growth of microvoids, Fig. 10, as was observed on specimens with homogeneous hydrogen distribution.

When the behaviour of specimens with preferentially orientated to the crack plane hydrides are compared to others with no preferentially oriented hydrides, although having similar hydrogen content, fracture mode and fracture toughness prove quite different. The following examples depict this difference in tests at room temperature. For homogeneous hydrogen distribution up to 2000 ppm, ductile fracture was observed, while for hydrogen content around 1000 ppm homogeneously distributed, the J initiation values were around 90 ± 5 kJ/m 2 . In these cases, the hydrides were not preferentially oriented. In contrast, for a non-homogeneous distribution with 1200 ppm near the cracktip where the hydrides were orientated to the crack plane, brittle fracture occurred at J C values around 15 kJ/m 2 , Fig. 5. Therefore, both high H concentration and preferentially oriented to the crack plane platelets influence the fracture mode, enhancing brittle behaviour and displacing the ductile to brittle transition region. It seems that the transition from ductile to brittle behaviour could be associated with both the hydride platelet orientation with respect to array on the crack path and the crack-tip region H content. A tensile stress field around a crack makes the hydride precipitate preferentially in an orientation favourable to crack propagation. An array of hydrides oriented preferentially to the crack plane near the crack-tip leads to high probability of finding a crack path mostly composed by hydrides, resulting in a sudden crack instability at a low fracture toughness value. When the matrix ductility is high enough (e.g. at high temperature), the matrix will deform around the fractured hydrides, thus avoiding brittle fracture and resulting in a stable crack growth and higher fracture toughness.

Summary

The influence of a hydrogen concentration gradient at the crack-tip on the fracture toughness, fracture mode, and micromechanisms of a Zircaloy-4 alloy was studied within a temperature range from 293 to 473 K. Hydride stress-orientation due to the mechanical load imposed during hydrogen charging was observed. Platelets preferential precipitation plane was LS rather than LT, the precipitation plane for unstressed condition.

In the heterogeneous H distribution case, brittle fracture without previous stable crack growth was observed at room temperature, brittle fracture with some previous stable crack growth occurred at 323 K, and the fracture mode was ductile at the highest testing temperature. In contrast, fracture behaviour was always ductile for homogeneous H distribution even at 2000 ppm H. The main fracture micromechanism observed was cleavage for specimens that fractured in a brittle manner at low temperatures, and nucleation, growth and coalescence of microvoids for specimens that failed in a ductile way, especially at high temperature.

A change in fracture mode from ductile to brittle occurred when the hydrides were concentrated ahead the crack-tip and oriented parallel to the crack plane at the lowest temperatures. The ductile-to-brittle transition is dependent not only on the crack-tip hydrogen content but, also on the hydride platelet orientation.

Fig. 1 .

 1 Fig. 1. Geometry and dimensions of employed CT specimens, W = 25.4 mm, B = 6 mm. The reference system corresponds to the specimen orientation related to the plate: L -long direction (rolling direction), T -transverse direction, and S -short direction.

Fig. 2 .

 2 Fig. 2. Variation of H content along the remaining ligament.

Fig. 3 .

 3 Fig. 3. Micrograph showing the hydride distribution in the same specimen (charged under stress), (a) near the crack-tip zone (1634 H ppm) and (b) far away from the crack-tip (>8 mm) and not far from the crack path (173 H ppm). Image (b) shows hydrides ÔagglomerationÕ forming arrays on the L direction, LT plane. (c) Micrograph from a specimen charged without stress.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Orientation of hydride platelets with (LS) and without (LT) mechanical loading.

Fig. 6 .

 6 Fig. 6. Fracture surfaces of specimens tested at room temperature.

Fig. 7 .

 7 Fig. 7. Fracture surface showing the appearance of a noncommon area between the fatigue pre-crack region and the fracture itself.

Fig. 8 .

 8 Fig. 8. J values at different temperatures as a function of the H level: (a) near the crack-tip, and (b) far away from the crack-tip (at more than 8 mm from the crack-tip, corresponding to the bulk H content, see Fig. 2).

Fig. 9 .

 9 Fig. 9. J values as a function of temperature.

Fig. 10 .

 10 Fig. 10. Fracture surfaces of specimens tested at 473 K.

Table 1

 1

	Chemical composition of material used
	Alloying elements (wt%)				Main impurities (ppm)
	Fe	Cr	FeCr	S	O	Al	C	Hf	Si
	0.21	0.1	0.31	1.42	0.125	25	110	52	39