Gwenaëlle Girault

Mathilde Gay

Sebastien Lobo

Laurent Bramerie

Michel Joindot

Jean-Claude Simon

Alexandre Shen

F Blache

H Gariah

F Mallecot

Quantum Dash Actively Mode-locked Fabry-Perot Laser Module demonstrated as part of a Wavelength Tunable RZ Transmitter

  L'archive ouverte pluridisciplinaire

Quantum Dash Actively Mode-locked Fabry-Perot Laser Module demonstrated as part of a Wavelength Tunable RZ Transmitter G. Girault, M. Gay, S. Lobo, L. Bramerie, M. Joindot, J.C. Simon A. Shen, F. Blache, H. Gariah, F. Mallécot, O. Le Gouezigou, F. Poingt, L.

Le Gouezigou, F. Pommereau, B. Rousseau, F. Lelarge, G-H. Duan

A quantum dash Fabry-Perot actively mode locked laser module is tested as part of a 42.7 Gbit/s transmitter with more than 10 nm wavelength tunability. Its low chirp level is also assessed through chromatic dispersion tolerance measurements.

Introduction: Mode-locked laser diodes (MLLDs) have a great potential for many applications in optical communications such as optical time-division multiplexing [START_REF] Arahira | 160-gb/s OTDM signal source with 3R function utilizing ultrafast mode-locked laser diodes and modified NOLM[END_REF] or optical clock recovery [START_REF] Renaudier | First demonstration of alloptical clock recovery at 40 GHz with standard-compliant jitter characteristics based on a quantum-dots self-pulsating semiconductor laser[END_REF][START_REF] Roncin | Multi-Data-Rate System Performance of a 40 GHz All-Optical Clock Recovery Based on a Quantum-Dot Fabry-Perot Laser[END_REF]. Furthermore, multi-wavelength lasers have become attractive as they enable to obtain a large number of wavelength-division-multiplexing (WDM) channels from a single source [START_REF] Nishikawa | Optical highfrequency pulse generation with wavelength tunability by using a DBR laser with vernier-type gratings[END_REF][START_REF] Chung | All-optical multiwavelength conversion of 10 Gbit/s NRZ/RZ signals based on SOA-MZI for WDM multicasting[END_REF]. Future WDM light sources must meet a number of criteria such as simple and stable operation, cost effectiveness and low energy consumption. Quantum Dash (QD) Fabry-Perot (FP) MLLDs are promising candidates as they meet all of the above criteria. In addition these components have a small footprint, hold possibilities for monolithic integration and have recently shown the possibility to generate short pulses over a wide spectral bandwidth through passive filtering [START_REF] Shen | Nearly Fourier-transform limited tunnel injection quantum dash mode-locked Fabry-Perot laser module for tunable pulse generation[END_REF]. This last reference deals with spectral and temporal characterizations performed thanks to the linear spectrogram method [START_REF] Dorrer | Investigation of the spectrogram technique for the characterization of picosecond optical pulses[END_REF].

In this Letter, a 42.7 Gbit/s RZ transmitter using the QD FP MLLD module is assessed thanks to Bit Error Rate (BER) measurements, and its chirp is evaluated through chromatic dispersion tolerance investigations.

Source description and characterisation: The chip and the module were fabricated at III-V lab [START_REF] Shen | Nearly Fourier-transform limited tunnel injection quantum dash mode-locked Fabry-Perot laser module for tunable pulse generation[END_REF]. The QD-based hetero-structures were grown by GSMBE on a S-doped (100) InP wafer. The active core consists in 6 InAs QDs layers and a GaInAsP injector quantum well which is separated by a thin InP barrier designed for phonon-assisted tunnel injection (TI) lasers emitting at 1.55 µm. The TI layers are enclosed within barriers and two separate confinement hetero-structure (SCH) layers. Along with the FP chip, a temperature probe, a Peltier cooler and a micro-wave V-type connector have been integrated into a butterfly module. The optical spectrum from the passively MLLD module is presented on Fig 1.

As shown in [START_REF] Shen | Nearly Fourier-transform limited tunnel injection quantum dash mode-locked Fabry-Perot laser module for tunable pulse generation[END_REF], the spectrum is quite flat and centred at 1555 nm, with a 3 dB width of 19 nm. This wide spectrum allows to achieve a tunable source using an external optical filter (represented by a dashed line on the left hand side of Fig. 1) Thus when using the module in a 40 GHz actively mode-locking configuration with a Gaussian optical filter with a 3 dB bandwidth of 0.6 nm, we obtained an optical clock showing a full width at half maximum (FWHM) of 6.3 ps and a timing jitter of 500 fs. The module was then driven by a 20 dBm, 42.7 GHz electrical clock and a bias current of 300 mA was used. Temporal characterizations, shown on the right hand side of Fig Finally the signal was analyzed after detection in an optically preamplified receiver followed by an electrical demultiplexer (ETDM).

Firstly we assess the tunability of the QD FP MLLD module through receiver sensitivity measurements. For such a purpose we measured the required receiver input power for a BER of 10 -9 (called sensitivity) at different wavelengths (from 1550 nm to 1560 nm) by tuning the optical filter (OF#1). We compared three kinds of transmitters: the reference source and the QD FP MLLD with two different Gaussian optical filters (OF#1) with 3 dB bandwidths of respectively 0.25 and 0.6 nm.

Results are presented on Fig. 3 where the sensitivity is depicted versus the wavelength. In the right hand side of the figure, the corresponding eye diagrams are displayed. The pulse FWHM of the reference source is 8.25 ps; those of the QD FP MMDL with respectively the 0.6 nm and 0.25 nm wide filters are 6.3 ps and 14.4 ps.

In the case of the reference source, the sensitivity is equal to -31 dBm over the entire C-band. In the case of the 0.6 nm filter, the sensitivity was measured at -32 ± 0.25 dBm over the studied spectral bandwidth (from 1550 nm to 1560 nm), showing a quite stable transmitter. The sensitivity improvement compared to the reference result is due to the shorter pulse width, a well known effect [START_REF] Idler | System performance and tolerances of 43 Gbit/s ASK and DPSK modulation formats[END_REF]. In the case of the 0.25 nm filter, the mean sensitivity is equal to -31.5 ± 1 dBm. The larger variation of the sensitivity might be linked to the higher wavelength selectivity of the filter. This BER assessment demonstrates for the first time the high quality of the QD FP MLLD module over 10 nm when used as part of a 40 Gbit/s transmitter.

Secondly, we evaluated the chirp of the QD FP MLLD module by investigating the influence of chromatic dispersion (CD) on the receiver sensitivity. Indeed in the case of an unchirped signal, the CD leads to a degradation of the receiver sensitivity (pulse broadening) and thus the penalty versus CD curve is symmetrical with respect to the normal or anomalous dispersion. On the contrary for a chirped signal, the CD can lead to pulse narrowing and thus improving the penalty and shift the minimum penalty from the zero dispersion point [START_REF] Agrawal | Nonlinear Fiber Optics[END_REF].

For this experiment we included a tunable CD module (centred at 1550 nm and with a 1 ps/nm step) in front of the receiver. Then we measured the penalty (for a BER of 10 -9 ) induced by the CD at 1550 nm for the three same configurations as previously. The penalty is defined here as the receiver sensitivity variation when chromatic dispersion is added.

Results are plotted on Fig. 4. As the chirp results into a shift of the minimum penalty from the zero CD point, we looked at this parameter. As the minimum penalty point is not easy to identify accurately, we introduced <CD> defined as the average value of the two dispersion values associated with a given penalty, here 3 dB. It equals +6.5 ps/nm in the reference case, showing a really low chirped transmitter. In the case of the 0.25 nm filter,  <CD>, equal to +5 ps/nm, and the CD tolerance are equivalent to the case of the reference signal. Now in the case of the 0.6 nm filter, the CD tolerance is weaker as the pulses are shorter, however <CD> is still close to zero (-5 ps/nm). To sum up, in both configurations of the QD FP MLLD module, offsets of the penalty curve are comparable to that of the reference and close to zero indicating a low chirped source module.

Conclusion:

In this paper we demonstrated the high quality of a QD FP actively MLLD module used in a 42.7 Gbit/s transmitter. The pulse width of the transmitter was adjusted using optical filtering. The receiver sensitivity, typically -32 dBm for 6.3 ps width pulses and -31.5 dBm for 14.4 ps width pulses, is comparable to that obtained with a classical RZ 33% optical source. Moreover we established that this module was tunable through a passive optical filtering over 10 nm without any consistent impact on the receiver sensitivity (less than 1 dB). Furthermore, we showed the low chirp of this simple transmitter through chromatic dispersion tolerance investigation demonstrating the high potential of this module for transmission applications. QD MLLD, 0,6 nm filter RZ 33% QD MLLD, 0,25 nm filter

Figure captions:Fig. 1

 1 Figure captions:

Fig. 2 Fig. 3

 23 Fig. 2 Experimental setup.

Fig. 4

 4 Fig. 4 Penalty at a BER of 10 -9 versus the chromatic dispersion.

Figure 1 Wavelength

 1 Figure 1

  Figure 3

  assessment: The experimental setup is presented on Fig 2.Either a reference source or the QD FP MLLD module was used in the transmitter in order to make a comparison. The reference source, delivering a classical RZ 33% signal, consisted in a tunable continuous laser followed by an x-cut LiNbO 3 modulator, driven by a 21.3 GHz sinusoidal clock. The QD FP MLLD module was driven by a 42.7 GHz electrical clock and a bias current (Ibias) of 300 mA. Both optical clocks (reference and QD FP MLLD) were coded via an x-cut LiNbO 3 modulator driven by a 2

1, were obtained with an optical sampling oscilloscope (with 1 ps time resolution, and 150 fs jitter resolution).

BER

31

-1 PRBS at 42.6 Gbit/s, generated by multiplexing electrically four trains at 10.66 Gbit/s. In the case of the QD FP MLLD, a tunable optical filter (OF#1) shaped the signal, according to the principle displayed on the left hand side of Fig 1.

Acknowledgments: This work is partly supported by the French ANR projects TONICS and PERSYST II.


Abir Zendagui
email: abirzendagui@etu.u-bourgogne.fr

Gaëtan Le Goïc

Hermine Chatoux

Jean-Baptiste Thomas

Jochum P c , Maniglier S d Yuly Castro

Marvin Nurit

Alamin Mansouri

Quality assessment of dynamic virtual relighting from RTI data: application to the inspection of engineering surfaces

Keywords: Material appearance, Reflectance Transformation Imaging, Psychometric evaluation, Visual assessment

This paper aims to evaluate the visual quality of the dynamic relighting of manufactured surfaces from Reflectance Transformation Imaging acquisitions. The first part of the study aimed to define the optimum parameters of acquisition using the RTI system: Exposure time, Gain, Sampling density. The second part is the psychometric experiment using the Design of Experiments approach. The results of this study help us to determine the influence of the parameters associated with the acquisition of Reflectance Transformation Imaging data, the models associated with relighting, and the dynamic perception of the resulting videos.

INTRODUCTION

Controlling the appearance of engineered surfaces is of primary importance for many high added-value industrial products. A major issue is that surfaces are often characterized through geometric specifications, [START_REF] Lemesle | Surface reflectance: An optical method for multiscale curvature characterization of wear on ceramic-metal composites[END_REF] and the functional relation between these geometric descriptors. [START_REF] Goic | Geometric quality and Appearance of surfaces: Local and global approach[END_REF] The optical properties of surfaces, and more generally the visual appearance (how they are perceived) are not explicitly considered. [START_REF] Pernkopf | 3d surface acquisition and reconstruction for inspection of raw steel products[END_REF][START_REF] Coules | Reflectance transformation imaging as a tool for engineering failure analysis[END_REF] Indeed, this interaction and its perception by people is related to a complex multi-physical and multi-scale process. In an industrial context, the Reflectance Transformation Imaging technique can be seen as simplification of BRDF 5 allowing to estimate the angular reflectance at each point/pixel of the inspected surface. Angular reflectance plays an essential role in the perception of surfaces, [START_REF] Malzbender | Enhancement of shape perception by surface reflectance transformation[END_REF] and its capture [START_REF] Zendagui | Quality assessment of reconstruction and relighting from rti images: Application to manufactured surfaces[END_REF][START_REF] Pintus | Objective and Subjective Evaluation of Virtual Relighting from Reflectance Transformation Imaging Data[END_REF][START_REF] Giachetti | Multispectral RTI Analysis of Heterogeneous Artworks[END_REF][START_REF] Ponchio | Relight: a compact and accurate rti representation for the web[END_REF] fills the gap between sensory evaluation and geometric measurement as illustrated in Figure 1. This process is bidirectional: either reproducing a design, i.e. a visual emotion, on the manufactured surfaces and/or, inversely, checking the conformity of sensory visual descriptors during or at the end of the production process.

Reflectance Transformation Imaging (RTI) acquisition consists in acquiring a series of images by a fixed camera positioned orthogonal to the studied surface while varying the light direction. Each image of the series corresponds to an illumination direction. The images are used afterwards in a modelling process to reconstruct the angular reflectance and relight the surface under any arbitrary illumination direction.

The objective of this study is to investigate the relevance of this information for the appearance inspection 11 of engineered surfaces. A previous study in this direction made it possible to understand the influence of the acquisition and the modeling parameters on the accuracy of the relighting 7 using the RTI data. This preliminary study highlighted an important limitation in the proposed experiment. This limitation is due to the fact that the psychometric evaluation was carried out on fixed images. Each image is associated with a fixed angle and the different study parameters. The assessment is therefore based on static images, which represents a significant and critical difference with the way in which humans carry out this type of evaluation. For example, in the industrial context of a visual sensory assessment, controllers are required to vary the incidence of light by changing the orientation of the surface to increase the robustness of the process. To address this limitation, this study implements the analysis of dynamic relighted surfaces to quantify the impact of the acquisitions parameters closely to the functioning of human perception. To this end, we have implemented a dynamic relighting psychometric experiment, which means the assessment is made from reconstructed videos of the surface appearance in which the light direction varies according to a defined illumination path, consistent with industrial visual control practices, where the controllers are required to follow an exploration protocol, which includes the variation of the azimuth and theta illumination angles (θ, φ).

In the present study, the first step aims to determine the acquisition parameters (Exposure time, Gain, Sampling density) per experiment surface (three engineered surface were used with different reflectance behaviour) using the Design of experiment Approach (DoE). The second step aims to evaluate the accuracy of dynamic relighting using RTI data where we noticed a preference of DMD as an approximation model to relight virtually a surface.

METHODS

The quality of RTI reconstruction depends on several parameters related to both the acquisition and the modeling process. In a non-exhaustive way, significant factors that we have identified are, for example, the acquisition scale, the spatial-angular resolution, the exposure parameters, the illumination wavelength, the implemented angular reflectance model, the number of descriptors implemented in the chosen model, or even the reconstruction light path. The quality of RTI reconstruction also depends on the intrinsic complexity of the surface reflectance behaviour, and interactions are potentially existing and significant between many of these factors. The surface samples, the RTI acquisition device and parameters, and the Design Of Experiment (DOE) proposed method and implementation are presented below.

Experiment Samples. For this study, three engineered surface samples have been selected because they show different reflectance behaviour. These samples are presented in figure 2. They are respectively an industrial paper sample, noted S 1 in the following. S 2 is a metallic sand blasted sample (Spectralon® SRS-05-020 Diffuse Reflectance Standard), and S 3 is a glossy metallic sample (brushed-polished surface). S 1 and S 2 are diffuse material contrary to S 3 .

RTI acquisition setup and parameters. The acquisitions were carried out with an in-house developed RTI system (figure 3a). The vision system is based on a monochromatic 2/3 active pixel-type CMOS sensor, allowing low-noise acquisitions with a resolution of 12.4 Mp (4112 × 3008). The camera is equipped with a precision micro-imaging modular lens system equipped with a motorized magnification and focus, which makes it possible to adapt the field of view and to quickly focus on the area of interest. Regarding the light source, we use a white uniform light (high power LED) equipped with a lens. In order to carry out RTI acquisitions, the light source is mounted on a motorized hoop allowing to control the angular positioning of the source (φ and θ angles). All the RTI acquisition parameters can be controlled with a user interface (figure 3b), which also allows the visualization of acquired data and post-processed features mappings. This software has been developed for the specific needs of this study to allow the automation of RTI acquisition cycles on a surface with different acquisition parameters, such as the zoom factor or the number/density of acquisition positions.

Evaluation and results

Experimental protocol 4.1.1. RTI setup

The acquisitions were carried out with a custom device (figure 9a) designed and developed to allow the implementation of the HD-RTI method. The image sensor is a monochromatic 2/3 active pixel-type CMOS sensor, allowing lownoise acquisitions with a resolution of 12.4 Mp (4112 ⇥ 3008). The camera is equipped with a precision micro-imaging modular lens system equipped with a motorized magnification and focus, which makes it possible to adapt the field of view and to quickly focus on the area of interest. Regarding the light source, we use a white uniform light (high power LED) equipped with a lens to limit the solid angle of the light to 15 degrees. In order to carry out RTI acquisitions, the light source is mounted on a motorized hoop of radius 35 cm capable of positioning the source in any position in the hemisphere of the same radius around the object by varying ( and ✓) angles. All the parameters of the HD-RTI acquisition process can be controlled with a user interface (figure 9b), which also allows the visualization of acquired data and post-processed features mappings.

Data set description

The method is applied on 3 samples with distinct properties in terms of isotropy and reflectance (figure 10). Paper is an industrial matte paper sample. Its roughness is isotropic. Plastic is a plastic polymer color-checker chart. This surface is an array of patches of di↵erent colors used for color calibration in photography. Each patch has a di↵erent roughness. Metal is a brushed metal sample (hot rolling). The acquisitions on these three samples were carried out with a field of view of 6.89 ⇥ 5.11mm (pixel size of approx. 1.7 µm 2 , with 149 lighting angles distributed homogeneously (see figure 5). A first step (in this study) consisted in defining the optimal setting of the RTI acquisition parameters not retained in the experimental design. Indeed, certain parameters, and particularly those related to illumination, can strongly influence the quality of the data obtained and thus bias the results sought through the psychometric experiment. The evaluation metric chosen to estimate the quality of an RTI acquisition is based on the number of not-measured points in the RTI acquisition, i.e. the sum of the number of under exposed pixels (grey level value of 0) and over-exposed ones (255). After various experiments on the surfaces presented above (figure 2), we chose to study in particular the choice of two parameters for this preliminary step: the angular density of acquisition and the exposure time.

Angular density of acquisition. Indeed, for certain surfaces and particularly anisotropic ones with oriented textures, the acquisition density can influence the chosen metric. For example, for a brushed surface, saturated values are observed when the illumination is orthogonal to the direction of the brushing texture (theta angle), and the number of saturated points will strongly depend on the chosen acquisition density and the width of the solid angle which generates the saturated values when the illumination is close to being orthogonal to the surface texture. The density of acquisition experiment values are 50, 100, 200 and 400, which corresponds to [50 -400] positions homogeneously distributed in the (θ, φ) angular space. The results obtained show that an asymptote is reached from a density corresponding to 100 acquisition positions, even for anisotropic surfaces. The tests carried out to determine the optimal adjustment of the acquisition time for each of the surfaces of the experimental plan were therefore carried out with this angular density value.

Exposure time. This parameter has a strong influence on the quality of the data obtained, and its adjustment is generally carried out arbitrarily by the operator in charge of the acquisition. The method then generally consists in choosing an exposure value compromise between the raking illumination angles and higher angles, close to 90 degrees, that induces more over-exposed pixels. For this experiment, we choose 10 exposure values linearly distributed between 2 and 125 milliseconds, which correspond to the upper limit of the lighting time of our device (the led works at high intensity mode, which needs to meet discontinuous operating mode). Moreover, the exposure time value must also be adjusted according to the magnification setting which influences the amount of light received by the sensor. Thus, several images were taken in different level of magnification (from 30% to 100%) to define the optimum exposure time value per magnification. The plot obtained (figure 4) shows the exposure time adjustment values which experimentally makes it possible to obtain a constant luminosity of the images captured at different magnifications with our RTI acquisition setup. As expected, it is observed that the sensitivity of the vision system decreases for large magnifications values. The evolution follows a quadratic law whose least squared approximation is given below in equation 1 (red curve in figure 4). Once the optimal exposure time has been estimated for a surface at an arbitrarily chosen zoom value, this law is then used to determine the optimal exposure time values associated with the different magnification values of the experiment. The results of this preliminary step are presented in figure 5 for samples S 1-3 . The exposure time used corresponds to the crossing of the curves of the evolution of the percentage of under-saturated points and oversaturated ones, which appears to be a robust value always close to the overall minimum of unmeasured points. The obtained results are synthesized in table 1 and table 2.

Psychometric experiment. The evaluation is performed through a psychometric experiment by evaluating the quality of reconstruction with the different settings of the Design Of Experiment (DoE) study parameters. Indeed, the quality of RTI reconstruction and its perception depends on several parameters, related to both data acquisition and information modeling and reconstruction. There are, in a non-exhaustive way: the acquisition scale, spatial and angular resolution, the lighting exposure time, the wavelength of light, the reconstruction model, the number of decomposition descriptors, the scale of decomposition, and the reconstruction light path, etc. The quality of RTI reconstruction also depends on the complexity of the surface reflectance behaviour. In this study, we have chosen, based on previous experiments, to focus on 4 parameters that seem to be particularly significant and difficult to control in the case of the RTI process. These parameters, noted P i , are: • P 1 , the reconstruction model: three different approximation models of angular reflectance are considered (Polynomial Texture Mappings (PTM), [START_REF] Malzbender | Polynomial texture maps[END_REF] Hemispherical Harmonics (HSH) [START_REF] Gautron | A novel hemispherical basis for accurate and efficient rendering[END_REF] and Discrete Modal Decomposition (DMD). [START_REF] Pitard | Reflectance-based surface saliency[END_REF] • P 2 , the angular density of acquisition: acquisitions are performed from 50, 100, 200 and 400 positions homogeneously distributed in the (θ, φ) angular space (see figure 6).

• P 3 , the scale of acquisition: the acquisition system controlling the magnification through a precision motorized lens. Two modalities are defined in the DOE (40 % and 80%).

• P 4 , the surface reflectance behaviour: three engineered surfaces, already presented above, noted S 1-3 have been chosen for this experiment.

The influence of the different light path on the reconstruction is not investigated in this study. The light path was defined and kept constant over all experiments, as shown in Figure 6a. The assessment is performed through a psychometric experiment where the participants were asked to choose the most similar video to the reference one. The reference video is associated to a light path in Figure 6.a, varying the azimuth and elevation angles (θ: from 0 • to 359 • , φ: from 6 • to 30 • ) and constructed from a dense raw RTI acquisition. The relighted video frame's are obtained from the virtual re-lightening of the surface in each light position of the light path used in the reference video. The original and the relighted videos were simultaneously displayed to the experiment screen divided into four parts (see Figure 7). The top-left tagged for the original video through all the experiment. The three other parts are randomly positioned relighted video from parameter P 3 . All the other parameters are fixed per sequence. Each parameter combination from the DoE is projected twice in the experiment. The experiment is divided into two parts: training and main session. The training session sequences (5 sequences) were chosen based on their temporal and spatial information, so users span the full range of temporal and spatial information of the test sequences in the experiment conditions. The video sequences duration is 13 s and randomly projected for each user. The psychometric experiment was designed using MATLAB® Design app and started with a data collection of participants information using a text form with: name, last name, age, research field and gender for statistics reasons. The form is followed by a text window to explain instructions and the task. Participants were accompanied in the first step (training session) and allowed to ask questions for a better comprehension of the instructions. The participant's choice is expressed using the keyboard's arrows (up, down, left arrows) and validate their choice using the return button to pass to the next sequence.

RESULTS AND DISCUSSION

The experiment was held in a dark room, nine volunteers (from the ImViA laboratory team) have participated (7 males, 2 females) with a mean age of 31 years old. All participants have a good visual acuity with correction if needed. All of them have an expertise in image processing and computer vision. The psychometric experiment took an average time of 20.86 minutes per participant, divided to 3.74 minutes for the training session and 17.11 minutes for the main session. The results are presented in the form of heat maps (see figure 8, where the grey cells represents the unmeasured data or outliers), for each of the samples of the study separately, and globally (by averaging the results obtained with all the samples).

The overall results show a significant difference between the appearance reconstruction methods used for RTI relighting. We observe that the PTM reconstruction technique is very rarely chosen (less than 5% of choices in average) as being the most similar to the reference video. The HSH technique is located in an intermediate way, whereas the Discrete Modal Decomposition is chosen, on average, much more often (approximately 60% of choices in average) as being the one which offers a perception closest to the reference video. These data confirm Concerning the effect of the angular density of acquisition, it is observed that the performance of the PTM relighting is globally stable (at a very low level) whatever the angular density of sampling of the RTI acquisition. However, even though outliers were removed from the experiment for the angular density of 50 (the lowest), the results associated with this value seem to show an increase in the performance of the PTM method for this type of acquisition (i.e. with a limited number of acquisition points). For the HSH method, there is a tendency for the performance of the method to decrease when the acquisition sampling is increased, and vice versa for the DMD method. This can be explained by the lower number of descriptors, i.e. basis vectors, used with the HSH technique (16) than with the DMD technique (50). Indeed, this higher number of descriptors allows, for dense acquisitions, a better approximation of the local variations (higher frequencies) of the angular reflectance in each pixel.

We also observe that the average performance of the HSH technique decreases when the magnification increases, and vice versa for the DMD technique. This tendency is significant in the overall results heat map. For example, for the acquisition density of 100, the HSH method is chosen on average at 40.74% with a magnification of 40% , and this percentage is 33.33% with a magnification of 80% , while the DMD method is chosen. at 55.55% at 40% magnification and 62.81% at 80% magnification. In our opinion, this is explained by the fact that the sensor integration effect is directly related with the size of the acquired pixel. This integration effect produces a smoothing of the information. Thus, for small acquisition scales, the complexity of the measured angular reflectance tends to increase, which explains the observed increase in performance of the DMD method in comparison with HSH relighting.

Finally, the analysis of the results surface sample by surface surface sample also highlights significant elements in terms of similarity of perception with the reference video. In particular, it is observed that the sample S 2 leads to a similar evaluation between HSH and DMD relighting techniques, with a slight advantage to HSH. Conversely, the samples S 3 and S 1 show a very significant difference between the two techniques, particularly the sample S 3 where the DMD technique is preferred up to four times more often (for example, for the angular density of 100, the HSH obtains 19.45% of preferences, while DMD obtains 77.77%). These results show that the complexity of angular reflectance can come from both variations in intensity (local specularities, case of sample S 3 ) and directionality (anisotropic texture, case of samples S 1 and S 3 ). In this case, advanced methods, such as DMD, make it possible to improve the perceived quality of RTI reconstructions. Conversely, for surfaces with simpler behavior in terms of reflectance, i.e more Lambertian, without directionality (case of the standardized reflectance surface, manufactured by very homogeneous sandblasting), methods such as the HSH technique provide good modeling performance.

CONCLUSION

This article firstly presented a method for setting the optimal acquisition parameters in RTI. In a second time it presented a psychometric evaluation of the perceived quality of dynamic appearance reconstruction (relighting) from RTI acquisitions on manufactured surfaces. The studied parameters are the angular density of sampling, the measurement scale and the reconstruction model. Three models are studied (PTM, HSH and DMD). In addition, the psychometric experiment is carried out on three surface samples with distinct reflectance behavior (paper, sandblasted metal, and brushed / polished metal). The results show the importance of the choice of the reconstruction model for the perception of the RTI relighting. The DMD reconstruction model shows the best perceptual performance of the proposed RTI video reconstructions. In addition, a strong interaction is observed between the angular acquisition sampling, the measurement scale, and the choice of the reconstruction model. This study shows that the use of the RTI approach in an industrial context to aid the evaluation of the visual quality of surfaces is relevant and promising, but it is necessary to investigate and control the influence of the parameters of acquisition and modeling in terms of perception. The influence of the light path should also be considered.

Figure 1 .

 1 Figure 1. The production line of a high quality manufactured surfaces: from the artistic creation to the production of the surface.

Figure 2 .

 2 Figure 2. Experiment Samples: (a) S1 -Industrial paper (b) S2 -Spectralon® Diffuse Reflectance Standard (SRS-05-020) (c) S3 -Metallic brushed/polished surface.

Figure 9 :

 9 Figure 9: Custom acquisition system for HD-RTI

Figure 3 .

 3 Figure3. Custom RTI acquisition system (ImViA Laboratory), from .[START_REF] Nurit | Hd-rti: an adaptive multi-light approach for the visual quality assessment of manufactured surfaces[END_REF]

12

 12 Figure3. Custom RTI acquisition system (ImViA Laboratory), from .[START_REF] Nurit | Hd-rti: an adaptive multi-light approach for the visual quality assessment of manufactured surfaces[END_REF]

Figure 4 .

 4 Figure 4. Exposure time (ET ) vs Magnification level (z), estimated from acquisitions on sample S3.

S

S 1 S 2 S 3 Figure 5 .

 235 Figure 5. Mean percentage of unmeasured pixels Vs Exposure time for RTI acquisitions of Samples S1, S2 and S3 (100% Magnification).

Figure 6 .

 6 Figure 6. a) Experiment light path: 397 positions (b-d) position of the RTI acquisitions for respectively 50, 100 and 200 positions respectively.

Figure 7 .

 7 Figure 7. Psychometric experiment: User Interface

Figure 8 .

 8 Figure 8. Psychometric experiment: Heat Maps of the percentage of choices for samples S1-3

Table 1 .

 1 Optimum Exposure time results for samples S1-3, with our RTI acquisition system, at magnification value of 100%.

			1	S 2	S 3
	Gain		0	10	18
	Magnification (z)	100	100	100
	Optimal Exposure Time (E T , in ms) 60.65 108.21 31.07
	Mean percentage of NM pixels (%)	0.013	0.41	11.52
			S 1	S 2	S 3
	Optimal Exposure Time	Magnification (40 %) 32.69 61.68 60.83
	(E T , in ms) (E T , in ms)	Magnification (80 %) 41.78 88.25 84.55

Table 2 .

 2 Optimum Exposure time results for samples S1-3, with our RTI acquisition system, for the used magnification levels: 40% and 80%.

ACKNOWLEDGMENTS

This work benefited of the funding of French National Research Agency (ANR) through the project NAPS (https://anr.fr/Projet-ANR-17-CE10-0005).


Abir Zendagui

Jean-Baptiste Thomas

Gaetan Le Goic

Yuly Castro

Marvin Nurit
email: marvin_nurit@etu.u-bourgogne.frgaetan.le-goic

Alamin Mansouri
email: alamin.mansouri@u-bourgogne.frjean.b.thomas

Marius Pedersen
email: marius.pedersen@ntnu.no

Gaëtan Le Goïc

Quality assessment of reconstruction and relighting from RTI images: application to manufactured surfaces

Keywords: Reflectance Transformation Imaging, psychometric evaluation, relighting

I. INTRODUCTION

The inspection of manufactured surfaces is generally a visuo-tactile process realized by human experts. It typically consists in changing the product surface orientation to change the light incidence angle, which seems to be the most influential parameter. To ensure the reproducibility and to reduce the variability of results and decisions to make [START_REF] Debrosse | Sensory perception of surfaces quality[END_REF] [START_REF] Guerra | Sensory metrology in the environment of visual qulity control[END_REF], the inspection process has to be automatised or at least semiautomatised (at least light source could be changed automatically). In this sense, the Reflectance Transforming Imaging (RTI) approach seems to be a very adequate technique to automatise and make this Non Destructive Inspection (NDI). It simplifies the inspection conditions since only light position changing is to be considered. This technique is widely used for image relighting, i.e. for reconstructing a scene under an arbitrary (virtual) direction of illumination with applications in Video-games, Cultural heritage and recently in manufacturing.

Agence Nationale de la Recherche through SUMUM and NAPS projects, The Norwegian Colour and Visual Computing Laboratory -NTNU RTI acquisition process consists in taking a set of images with a fixed camera, orthogonal to the inspected surface, while varying the light direction at each image capture. From this set of images, each pixel is described by a set of discrete values (measured gray-levels, which are considered to be proportional to the luminance [START_REF] Durou | Shape from shading -Eclairages, réflexions et perspectives[END_REF]). To model the surface visual appearance continuously (corresponds also to relighting), these discrete values can be approximated (for each pixel) by projecting the vector containing the discrete measured gray-levels on a set of basis functions. The main methods used to model this information are Polynomial Texture Mappings (PTM, based on 2nd order polynomial functions) [START_REF] Malzbender | Polynomial texture maps[END_REF], Hemispherical Harmonics (HSH) [START_REF] Gautron | A Novel Hemispherical Basis for Accurate and Efficient Rendering[END_REF] and a more recent method named Discrete Modal Decomposition (DMD) [START_REF] Pitard | Discrete Modal Decomposition: a new approach for the reflectance modeling and rendering of real surfaces[END_REF].

In this study, we investigates the accuracy of the reconstruction/relighting from RTI-base images using the three models: PTM, HSH, DMD objectively using the Peak Signal to Noise Ratio (PSNR) and Structural Similarity error metric (SSIM) [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] and subjectively using a psychometric experiment for comparing models with different density of sampling in order to choose the best configuration of parameters to reproduce a perceptual aspect of physical manufactured surfaces.

II. BACKGROUND

RTI consists in acquiring photometric stereo images corresponding to varying light positions in two directions (elevation and azimuth) with a sensor positioned orthogonal to the surface. According to Durou et al. to [START_REF] Durou | Shape from shading -eclairages, ŕeflexions perspectives[END_REF], the luminance of a pixel is proportional to the grey intensity of the pixel and the formula is written in equation 1

G = k * L(U c ) (1)
Where k is the proportionality factor. L(U c ) is the luminance in sensor's direction. The result of a RTI acquisition is, for each pixel, a vector containing the luminance values corresponding to the different light positions used during the acquisition stage. To be able to relight a scene under any arbitrary direction an approximation (reconstruction) operation is necessary to build a continuous surface from the sparse point cloud (lunimance values). In the literature, three approximation models are mainly used to achieve this reconstruction: the Polynomial Texture Mapping (PTM), the Hemispherical Harmonics (HSH) based decompostion and the Discrete Modal Decomposition (DMD).

A. Polynomial Texture Mapping (PTM)

The Polynomial Texture Mapping was developed initially to improve the photo-realistic rendering in 3D imaging [START_REF] Malzbender | Polynomial texture maps[END_REF]. Widely used in the cultural heritage [START_REF] Macdonald | Polynomial texture mapping and 3d representations[END_REF], PTM approximate the the luminance surface in each pixel by a polynomial regression of 6 vectors (a 0 -a 5 ). These 6 vectors are sought for each pixel from the quadratic surface representing the reflectance defined by the second degree polynomial equation below:

L(l u , l v ) = a 0 + a 1 l u + a 2 lv + a 3 l u l v + a 4 l 2 u + a 5 l v (2)
The equation ( 2) is solved by a regression using the least squares formulation. The PTM's description of the surface is thus still simple and very smooth (no sharp frequency), but it demonstrates its efficiency with a small stack of images to correctly relight the surfaces, especially lambertian ones.

B. Hemispherical Harmonics (HSH)

The Hemispherical Harmonics is an improvement of the PTM technique [START_REF] Westin | Predicting reflectance functions from complex surfaces[END_REF] that implements more complex basis functions [START_REF] Gautron | A novel hemispherical basis for accurate and efficient rendering[END_REF], [START_REF] Kautz | Fast, arbitrary brdf shading for low-frequency lighting using spherical harmonics[END_REF], [START_REF] Ramamoorthi | An efficient representation for irradiance environment maps[END_REF], [START_REF] Ludwig | Environment Map Based Lighting for Reflectance Transformation Images[END_REF] P m l (cosθ) = P m l (2cosθ -1) and θ ∈ 0,

π 2 (3)

C. Discrete Modal Decomposition (DMD)

The Discrete Modal Decomposition is a technique initially used in dynamic mechanics to decompose a vibration, for the geometrical tolerancing [START_REF] Adragna | Analysis of shape deviations of measured geometries with a modal basis[END_REF] [16], multi-scale analysis of topographic surface roughness [START_REF] Julien | Analyse multi-échelle des composants articulaires d'une prothèse de hanche par la méthode modale[END_REF] [18] [START_REF] Goic | Multi scale modal decomposition of primary form, waviness and roughness of surfaces[END_REF] and for heat source evaluation [20] [21]. This method was then applied in the field of surface appearance, to model the local behaviour of the reflectance angular component, at each point/pixel. The DMD technique consists in decomposing the signal which represent the measured surface into a family of functions (eigen-modes) resulting from the vibration of basis forms (planar, cylindrical, spherical. . . ) and forming the eigen-basis. The eigen-basis is defined by eigen-vectors called Modal vectors. The modal bases derives from the resolution of the classical dynamic equation (4). K.q + M.q " = 0 (4) With q = q(x, y, t)

= +∞ i=1 Q i (x, y).cos(ω i t) (5)
Where M is Mass matrix, K is the stiffness matrix and q(x, y, t) is the displacements that characterizes the modal shapes. The projection of the measured surface onto the modal basis is then the sum of a linear combination of modal vectors and the decomposition residual (6)

V = Nq i=1 λ i Q i + (N q ) (6)
With Q i is the modal vectors composing the modal basis. λ i is the modal coordinates and N q is the number of modes of decomposition.

In Figure 1, reconstruction results for the same pixel using the three approximations PTM, HSH and DMD are presented. The RTI acquisition device we used in this study is a home-made one built in a manner to be able to position the light source precisely for any couple θ and φ since the light source is fixed on a rotated hoop. This configuration allows to construct a virtual lightening dome around the surface with flexible positioning of the light source. The camera and the lenses were chosen to ensure high quality and high resolution images. Fig. 2. Home-made x-RTI system for surface appearance acquisition.

2) Samples: The manufactured surfaces used in this study are ten metallic sides of manufactured rings. These surfaces, called LGR, were polished and finished with a tribofinition process, composed of three main successive treatments: Softening, polishing, brightening. The ten rings are grouped two by two in five process groups. The only difference between groups is the time that surfaces passed in the softening step and it is varying from 40 minutes and 80 minutes which makes the rugosity varying from 0.22µm to 0.17µm.

3) Data preparation: To build our dataset for the experiment, a dense (650 measurements: 130 in θ and 5 in φ) photometric stereo acquisition was taken with the system described previously for three plane sides of the ten rings presented in section III-A2. For each model the contributions are calculated with 6 terms for PTM, 9 harmonies for HSH and 50 modes for DMD. To estimate the visual quality of the RTI-based reconstructions, the direction of the virtual lighting used is the same for the measured image, see example in

B. Objective evaluation

The first metric we used for the objective evaluation is the Peak Signal to Noise Ratio (PSNR) since it has already proven its relevance in RTI data. For instance, it has already been used to evaluate the efficiency of the fitting models reproductions [START_REF] Pitard | Discrete Modal Decomposition for surface appearance modelling and rendering[END_REF] [23], to evaluate the RTI compressed images quality to be visualised on a mobile device [START_REF] Schuster | Rti compression for mobile devices[END_REF] or on a web pages [START_REF] Ponchio | A compact representation of relightable images for the web[END_REF]. This pixel-to-pixel quality metric is used in this study in perspective to correlate it with subjective experiment and to continue the investigation on the perceptual quality of the fitting models in [START_REF] Pitard | Discrete Modal Decomposition for surface appearance modelling and rendering[END_REF]. In Table II, the peak signal-to-noise ratio for the reproduced image using the three approximation models with the reference (captured photos). For each sampling density (All, Half, Fifth, Seventh), the images were reproduced in 5 light directions to study the impact of the exposure of the surface to the light on the reproduced image. In Table III the PSNR is computed over the 5 surfaces corresponding to the 5 fabrication process, using the DMD with the chosen sampling densities.

In a second step for the objective evaluation, we added a new metric namely Structural Similarity (SSIM) index, which seems to take more in consideration the spatial aspect and consequently the human visual perception. This metric uses in its computation three terms: luminance, the contrast, and structural term [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF]. Same for PSNR, we computed the structural similarity in Table VII of the reproduced images using the three approximation models in different sampling density (All, Half, Fifth, Seventh) in 5 light directions. In Table VIII the SSIM is computed only using DMD with different density sampling for the 5 used surfaces in the study. Each surface were token from a different group process to investigate on future studies the impact of the fabrication process on the evaluation of the surfaces. In order to see the mean and the median ranking of models, the Table V were computed to see the behaviour of the quality of the reproductions over the approximation models, and over sampling densities in the Table VI.

C. results and discussion

We notice in Table I that the global mean and medium of PSNR can rank easily the reproduced images quality as PTM, HSH, DMD respectively from the worst to the best. This rank is less obvious with 1/7 sampling density between HSH and DMD using PSNR in Table II where HSH and DMD are having respectively 30.78 and 30.37 in (θ : 180 • , φ : 14.5 • ) and 25.83, 25.89 in (θ : 180 • , φ : 28 • ) and 23.29, 23.20 in (θ : 88.61 • , φ : 41.5 • ). The same results for SSIM in Table VII were the difference between the models is reduced using 1/7 of the measurements (93 image) to at least 0.03. Regarding the quality of reproductions using DMD over the 5 manufactured surfaces in Table III using PSNR and in Table VIII using SSIM, a higher scores where for 1/5 of measurements in two angles (θ : 180 • , φ : 28 • ) and (θ : 180 • , φ : 55 • ) for all the surfaces which mean a higher quality of reproduction with this sampling density in this two angles comparing to the captured image. For the angle (θ : 180 • , φ : 41.5 Despite the results of PSNR and SSIM that DMD outperforms the other models and some other conclusions that we can make considering these results, it's important to mention that these error metrics are a pixel to pixel comparisons and no consideration to the biological human vision system, for that reason a further psychometric experiment were done to evaluate subjectively the perceived quality of the reproductions using RTI protocol.

D. Visual assessment

To evaluate the performance of the reconstruction (fitting) models (PTM, HSH, DMD) to reproduce the visual aspect of the manufactured surfaces used in the study, a psychometric experiment was designed in two parts to answer two main goals:

• The impact of the fitting models on the reproduction of the aspect of manufactured surfaces. • The impact of the sampling density to approximate perceptual aspect of surfaces. 1) Controlled psychometric experiment: The psychometric experiment were carried out on a dedicated web-based platform for psychometric scaling experiments: QuickEval [START_REF] Ngo | QuickEval: a web application for psychometric scaling experiments[END_REF]. In order to increase the perceived contrast of images with the environment, a black background was picked-out for the webpage and the experiment held in a dark room (no lighting with matt grey wall). The experiment is a paired comparison [START_REF] Engeldrum | Psychometric scaling: a toolkit for imaging systems development[END_REF] where participants' task is to choose the most similar representation (the reconstructed images using two of the three fitting models each time, positioned on the left and on the right) to the original image (the acquired image positioned in the middle). The experiment was setup as a forced-choice experiment with flipping of the pairs. The experiment divided on two parts took in average 30 minutes (10 minutes for the 1 st part and 20 minutes for the 2 nd part) for 240 comparison in total. The two parts answer respectively to the main goals presented in section III-D.

2) Participants: The participants involved in this experiment were 19 volunteer researchers working at the Norwegian Colour and Visual Computing Laboratory. They can be considered as experienced observers in psychometric experiments.

3) Results and discussion: In Figure 6, a comparison of users' preferences comparing the three approximation mod-Fig. 5. Screen shot of one example of the proposed experiment in QuiqEval platform, in the center the original image, in left and right two reproduced images using two of the three approximation models els using all the measurements (All the acquired Images to approximate the visual aspect) of the acquisition. The five groups of the graph are according to the chosen five angles for the experiment and as shown, more than 60% of users picked out DMD as the more accurate reproduction to the reference shown (raw image). Roughly 30% of them chose HSH which is a good evaluation of the quality of images using this model, but we believe that this amount of votes is due to the cases were the user had to choose between HSH and PTM, this latter is a model rarely chosen for the quality of reconstructed images. Only 6% of user's preferences were for PTM in (θ : 180 • , φ : 14.5 • ) using all the measurements and 1% using one seventh of measurement in the same angle and 0% in (θ : 180 • , φ : 55 • ) and this is due to poorness of PTM to approximate the luminance of the pixel defined by a quadratic surface (simple geometry) and which average the pixel's luminance under exposed (θ : 180 • , φ : 14.5 • ) and overexposed (θ : 180 • , φ : 55 • ) to the light. In (θ : 88.6154 • , φ : 41.5 • ) PTM do a better in Table 7, comparing to its performance in the other cases, the reason is that PTM's regression formulation (Least squares) is more efficient with less information : minimum residual. Fig. 6. Number of times a reproduction has been selected using the three fitting models : DMD, HSH, PTM using all the measurements These results are detailed and extended in Table IX for all the angles and sampling density used in the study, where the Fig. 7. Number of times a reproduction has been selected using the three fitting models : DMD, HSH, PTM using 1/7 of the measurements standard deviation to the mean is ∼ 2 for DMD in 15/20 of cases which makes it two standard deviation better than the average of the votes. And ∼ -2 for PTM in 15/20 of cases. In the other cases, and as explained before the gap between the performance of models decreases with less measurements and increases in some angles (phi = 14.5 • and 55 • ).

The Table XIII shows that the global mean and median of all the tested cases in the experiment, rank the accuracy of the reproduced images as PTM for the worst similarity to the original image, HSH as more accurate reproduction and DMD as the most accurate reproduction to the aspect of the surface. This conclusion match with results using the objective quality metrics (PSNR and SSIM). The second part of the experiment was devoted to study the impact of the sampling density as explained in III-D1. Considering the number of times a reconstruction has been selected in Figure 8 computed using one surface and roughly the same result for the other four surfaces, we notice that more measurements does not mean necessarily a better visual aspect. In the grazing angle (θ : 180 • , φ : 14.5 • ) a higher sampling density where chosen: ALL= 30% of votes, HALF= 30%, FIFTH=30% as a higher quality representation of the original image. This ranking is changed for other angles (higher elevation) and less density of measurements is more accurate to the acquired image, 45% of votes were for the lower density (1/7) for the three angles and less then 15% for the highest sampling density and this is due to the extra perceived information that RTI provides to the image by integrating in computing the pixels' luminance for one angle, the luminance of the pixel in all lighting directions. And since the task of users is to choose the closer representation to the reference, the extra information disturbs their perceptual evaluation of the reconstructed images.

The results of mean z-score for the five angles and the five surfaces chosen for the experiment is presented in Table XII. For the different sampling density tested, the standard deviation to the mean is higher for a 1/7 of the measurements for the three angles (θ : 180 • , φ : 41.5 • ),(θ : 180 • , φ : 55 • ),(θ : 88.61 • , φ : 41.5 • ) for all the surfaces. A full measurements density were preferred for a lower elevation (θ : 180 • , φ : 28 • ),(θ : 180 • , φ : 14.5 • ) because in these angles, images are dark, so the perceptual evaluation is biased due to the limitation of the sensor which represent the values of pixel in only 1 byte. A future improvement of this investigation is to extend our acquisition to be able to build high dynamic range images. Computing the the global mean of the z-score in Table XIII for all the surfaces and the angles for each density ranks the quality of the reconstructed images using different sampling density respectively as: All, Half, Fifth and the seventh from the worst to the best perceived quality.

IV. CONCLUSION

The study presented in this paper aims to evaluate the quality of the reconstruction and relighted images using RTI stereo photo-metric images. The reconstruction models used: PTM, HSH, DMD with different sampling densities. The evaluation of the relighted images is done first with objective error metrics: PSNR and SSIM and second with psychometric experiment to involve the human subjective evaluation of surfaces as a reference for the perceptual evaluation of manufactured surfaces aspect. The subjective and objective evaluation of the surfaces' aspect shows that DMD outperforms PTM and HSH as a the most accurate model to approximate the perceptual aspect of a surface. The results shows also that the quality of the reconstruction using RTI-based images decreases with a high sampling density due to disturbance of the perceived quality with the amount of the not perceptual information to the human vision system.

  derived from Hemispherical harmonics. Gautron et al. defined the hemispheric harmonic bases H m l (θ, φ) as an adapted version of spherical harmonies S m l (θ, φ) that rely on derivation of Legendre polynomials (3).

Fig. 1 .

 1 Fig. 1. Reconstruction of the same pixel using, from left to right, using PTM, HSH and DMD modeling

Figure 4 .

 4 Fig. 3. Manufactured rings

Fig. 4 .

 4 Fig. 4. Example of a measured image with RTI-based reconstructions in same lighting direction (θ : 180 • , φ : 41.5 • )

Fig. 8 .

 8 Fig. 8. Number of times a reproduction has been selected using four sampling density : All, Half, Fifth and a Seventh of the measurement using one of the five manufactured surfaces: A01

  • ) in Table VIII the 1/2 sampling density where more accurate to the captured image using SSIM and PSNR for the five surfaces.

TABLE I GLOBAL

 I MEAN AND MEDIAN OF PSNR COMPUTED OVER 5 REPRODUCED SURFACES' ASPECT THE 5 LIGHT DIRECTIONS

	Method	Mean Median
	PTM	17.46	19.46
	HSH	24.52	23.47
	DMD	26.76	25.67

TABLE II PEAK

 II SIGNAL-TO-NOISE RATIO FOR THE REPRODUCED IMAGES USING PTM, HSH, DMD AND 4 SAMPLING DENSITY (ALL, HALF, FIFTH AND SEVENTH) WITH THE ACQUIRED IMAGES IN DIFFERENT ANGLES.

	Sampling	Angle	PTM	HSH	DMD
		θ : 180 • , φ : 14.5 •	12,13	31,06	34,87
	All (1):	θ : 180 • , φ : 28 •	20,65	26,17	29,11
	650	θ : 180 • , φ : 41.5 •	19,46	22,64	25,43
	position	θ : 180 • , φ : 55 •	14,90	19,52	22,37
		θ : 88.61 • , φ : 41.5 •	20,26	23,60	26,06
	Half (1/2): 324 position	θ : 180 • , φ : 14.5 • θ : 180 • , φ : 28 • θ : 180 • , φ : 41.5 • θ : 180 • , φ : 55 • θ : 88.61 • , φ : 41.5 •	12,10 20,66 19,47 14,87 20,26	31,03 26,17 22,66 19,51 23,57	34,84 28,83 25,53 22,23 25,82
	Fifth (1/5): 130 position	θ : 180 • , φ : 14.5 • θ : 180 • , φ : 28 • θ : 180 • , φ : 41.5 • θ : 180 • , φ : 55 • θ : 88.61 • , φ : 41.5 •	12,09 20,66 19,47 14,89 20,25	30,91 26,28 22,63 19,60 23,37	33,62 30,14 24,69 22,99 24,21
		θ : 180 • , φ : 14.5 •	12,10	30,78	30,37
	Seventh	θ : 180 • , φ : 28 •	20,59	25,83	25,89
	(1/7): 93	θ : 180 • , φ : 41.5 •	19,42	22,50	24,42
	position	θ : 180 • , φ : 55 •	14,80	19,31	20,60
		θ : 88.61 • , φ : 41.5 •	20,25	23,29	23,20

TABLE III PEAK

 III SIGNAL-TO-NOISE RATIO FOR THE REPRODUCED IMAGES USING DMD AND 4 SAMPLING DENSITY (ALL, HALF, FIFTH AND SEVENTH) WITH THE ACQUIRED IMAGES IN DIFFERENT ANGLES FOR THE FIVE MANUFACTURED SURFACES

	Surfaces	Angles	ALL	HALF FIFTH	SEVEN
		(180 • , 14.5 • )	34,87	34,84	33,62	30,37
	Surface 1:A01	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	29,11 25,43 22,37	28,83 25,53 22,23	30,14 24,69 22,99	25,89 24,42 20,60
		(88.61 • , 41.5 • ) 26,06	25,82	24,21	23,20
		(180 • , 14.5 • )	36,10	36,05	34,80	31,52
	Surface 2: A04	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	31,14 27,29 23,75	30,84 27,38 23,61	32,26 26,76 24,55	27,92 26,30 21,69
		(88.61 • , 41.5 • ) 24,71	24,48	22,81	21,86
		(180 • , 14.5 • )	35,21	35,15	34,09	30,90
	Surface 3: A06	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	30,64 27,96 24,74	30,39 28,04 24,59	31,45 27,14 25,35	27,26 26,93 22,63
		(88.61 • , 41.5 • ) 23,56	23,30	21,73	20,65
		(180 • , 14.5 • )	36,37	36,31	35,15	31,55
	Surface 4: A07	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	30,32 26,76 23,59	29,98 26,88 23,41	31,61 26,17 24,48	27,40 25,85 21,60
		(88.61 • , 41.5 • ) 26,01	25,68	24,17	23,17
		(180 • , 14.5 • )	37,20	37,16	35,83	32,63
	Surface 5: A10	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	31,88 27,49 23,35	31,47 27,60 23,14	32,58 26,50 23,98	28,51 26,00 21,53
		(88.61 • , 41.5 • ) 26,58	26,28	24,86	23,94

TABLE VIII STRUCTURAL

 VIII SIMILARITY FOR THE REPRODUCED IMAGES USING DMD AND 4 SAMPLING DENSITY (ALL, HALF, FIFTH AND SEVENTH) WITH THE ACQUIRED IMAGES IN DIFFERENT ANGLES FOR THE FIVE SURFACES.

	Surfaces	Angles	ALL HALF FIFTH	SEVEN
		(180 • , 14.5 • )	0,88	0,88	0,84	0,74
	Surface 1:A01	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,83 0,76 0,66	0,82 0,77 0,64	0,84 0,70 0,72	0,69 0,69 0,51
		(88.61 • , 41.5 • )	0,75	0,73	0,64	0,59
		(180 • , 14.5 • )	0,88	0,88	0,84	0,74
	Surface 2: A04	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,84 0,76 0,65	0,82 0,77 0,63	0,86 0,71 0,72	0,70 0,69 0,49
		(88.61 • , 41.5 • )	0,69	0,66	0,56	0,50
		θ : 180 • , 14.5 • )	0,84	0,84	0,80	0,72
	Surface 3: A06	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,84 0,79 0,68	0,82 0,80 0,66	0,84 0,73 0,73	0,70 0,71 0,53
		(88.61 • , 41.5 • )	0,68	0,65	0,56	0,50
		(180 • , 14.5 • )	0,90	0,89	0,85	0,76
	Surface 4: A07	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,85 0,75 0,64	0,83 0,76 0,61	0,86 0,69 0,72	0,71 0,67 0,48
		(88.61 • , 41.5 • )	0,71	0,67	0,58	0,53
		(180 • , 14.5 • )	0,91	0,90	0,87	0,79
	Surface 5: A10	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,86 0,76 0,63	0,84 0,77 0,59	0,87 0,69 0,71	0,73 0,66 0,47
		(88.61 • , 41.5 • )	0,76	0,73	0,66	0,62

TABLE X GLOBAL

 X MEAN AND MEDIAN OF Z-SCORE COMPUTED OVER 5 REPRODUCED SURFACES' ASPECT IN 5 LIGHT DIRECTIONS

	Method	Mean Median
	PTM	-1.87	-2.12
	HSH	-0.02	-0.05
	DMD	1.90	2.13

TABLE XI GLOBAL

 XI MEAN AND MEDIAN OF Z-SCORE COMPUTED OVER 5 REPRODUCED SURFACES' ASPECT IN 5 LIGHT DIRECTIONS USING ALL, HALF, FIFTH AND SEVENTH OF MEASUREMENTS

	Method	ALL	HALF FIFTH SEVENTH
	Mean	-0.65	-0.43	0.50	0.58
	Median	-0.77	-0.56	0.48	0.89

TABLE XII Z

 XII -SCORE FOR THE REPRODUCED IMAGES USING DMD AND 4 SAMPLING DENSITY (ALL, HALF, FIFTH AND SEVENTH) WITH THE ACQUIRED IMAGES IN DIFFERENT ANGLES FOR THE FIVE SURFACES.

	Surfaces	Angles	ALL	HALF FIFTH SEVEN
		(180 • , 14.5 • )	0,19	0,41	0,38	0,97
	Surface 1:A01	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,27 0,75 0,63	0,16 0,73 -1,04	0,32 0,30 0,22	0,11 1,19 1,44
		(88.61 • , 41.5 • )	0,95	-1,01	0,57	1,38
		(180 • , 14.5 • )	0,67	0,17	0,35	-1,19
	Surface 2: A04	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,68 0,91 -1,06	0,16 -1,04 0,52	0,51 0,30 0,67	0,34 1,66 0,90
		(88.61 • , 41.5 • )	0,68	0,89	0,67	0,90
		(180 • , 14.5 • )	0,49	0,11	0,46	-1,06
	Surface 3: A06	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,35 -1,00 -1,12	0,51 0,92 0,84	0,75 0,73 0,38	0,11 1,19 1,57
		(88.61 • , 41.5 • )	-1,00	0,79	0,51	1,28
		(180 • , 14.5 • )	0,38	0,67	0,06	-1,12
	Surface 4: A07	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,78 0,97 -1,19	0,44 0,38 0,69	0,95 0,52 0,38	0,27 0,84 1,50
		(88.61 • , 41.5 • )	0,87	0,66	0,66	0,87
		(180 • , 14.5 • )	0,05	0,49	0,34	0,77
	Surface 5: A10	(180 • , 28 • ) (180 • , 41.5 • ) (180 • , 55 • )	0,83 -1,28 -1,37	0,16 0,89 0,58	0,71 0,73 0,38	0,28 1,44 1,57
		(88.61 • , 41.5 • )	-2,01	0,66	1,09	1,70


B Bou-Saïd

G Grau

Keywords: foil bearings, rotor dynamics, stability, gas lubrication

On Nonlinear Rotor Dynamic Effects of Aerodynamic Bearings With Simple Flexible Rotors

The last decades have experienced a growing enhancement of aeronautical oil free turbomachinery. The classical linear approach of rotor dynamics commonly uses stiffness and damping coefficients to model journal bearings. In the present study, a nonlinear time dependant calculation is used for the dynamic simulation of a rotor mounted with aerodynamic (gas) bearings. A comparison between the two approaches indicates that the dynamic behavior of such bearings can be nonlinear in operating ranges where the rotor eccentricity reaches high values. In that case, the linear approach may lead to incorrect results and the nonlinear approach should be performed for better rotor dynamic prediction. A numerical procedure which analyzes the dynamic behavior of simple flexible rotors taking into account the nonlinear (transient regime) characteristics of aerody-namic bearings is presented. A simple example highlights the needs of nonlinear simu-lations in order to predict dynamic performance in oil-free turbomachinery.

Introduction

The designers of turbine engines have always had to deal with very severe restrictions concerning weight. Consequently, these designers focused on an increase in the specific power involving a reduction in the size of the machines, and an increase in rotating speed. In ranges of 30,000 to more than 100,000 rpm, an aerodynamic ͑gas bearing͒ approach seems well adapted compared to the traditional technological elements which cannot meet certain requirements, e.g., those concerning adequate lifetime. Thus air bearings have received considerable attention in these last decades and find increasing applications whenever elevated speeds, low loads, and a high degree of accuracy are necessary.

There are two operating conditions for aerodynamic bearings which are quite distinct: a dry contact-mixed lubrication mode, and an aerodynamic lubrication mode. The first condition corresponds to the transient mode and specifically the start-stop phases. This study deals with the latter behavior, when a complete film of air separating the rotating member from the static member is created. The determination of the pressure field in the bearings from the solution of the Reynolds equation in the case of compressible fluids, in particular with perfect gases, has been intensively investigated ͓1-3͔. Air bearing static behavior has been previously treated by several authors ͓3-11͔.

Once static properties are known, designers must investigate the dynamic properties of the system consisting of the shaft and the bearings. In this paper, we will first treat the dynamic behavior of a rigid shaft mounted symmetrically on two aerodynamic bearings. Two methods of calculation, linear and nonlinear, are presented and compared in this simple case. This comparison reveals differences in results according to the loading, particularly in terms of speed instability and response to an unbalance excitation force. This portion of the study will lead us to suspect that the linear approach may be inadequate in some range of operating conditions compared to the nonlinear approach. Next, an improve-ment taking into account damping within the deformable structure of the bearings is highlighted. The last part of this study enables us to improve our understanding of the dynamic behavior of a deformable shaft by including nonlinear aspects of the problem.

A first model is based on a simplified approach to the problem. The step by step calculation is validated by comparison with the traditional linear rotor dynamics analysis. An improvement of this model is then made possible using the Rayleigh-Ritz method based on three elementary vibration modes taken for the shaft. Thus, we can evaluate the influence of the flexibility of the rotor on its critical engine failure speeds, and also on the level of stability of the system. The encouraging results of this modeling lead us to a more complicated finite element method approach. In our conclusions we will point out the principal results and place the present contribution in context. Finally, some directions of future research will be proposed.

Bump Foil Bearing

Many numerical and experimental studies have been concerned with the static and dynamic characteristics of so-called bump foil bearings. These works are generally concerned with two separate issues:

• Tribological phenomena related to the stop/start phases. Dry friction occurring during these phases is the main cause of deterioration of the bearings. The texture and shape of the surfaces and the materials in contact play a major role in these phenomena ͓12-15͔; • Development of models to predict aerodynamic bearing performance in terms of carrying capacity, power loss and dynamic behavior ͓3,7,16,17͔.

This work is a contribution to the second research area. The bump foil technology allows us to take a very simple structural approach to the analysis. The upper foil ensures the continuity of the profile. The geometry of the bearing and the structural behavior are governed by the springs on which the upper foil is pressed, see Fig. 1.

Architecture and Geometry.

Bump foil bearings present several alternatives. Indeed, the number of lobes and foils ͑for which 1 Corresponding author.

1 explain the higher damping of the foil bearings is based on the presence of internal friction inside the deformable structure. Consequently, to add damping, some researchers assumed it was sufficient to increase the friction coefficient between the foils by using various types of coatings. Peng and Carpino ͓47͔ presented an analytical study in which they showed that damping in air bearings increases when the air film stiffens. As the eccentricity or the load increases the air film becomes stiffer and thus activate the dry friction between the bump and the housing. This phenomenon induces additional damping. They noticed that when the air film is flexible, i.e., when the load is light or the rotational speed is relatively low, ball bearings offer more damping than foil bearings. On the other hand, for higher load and speed, when the film of air is stiff, foil bearings provide more damping than their rigid counterparts due to the damping generated by friction between the foils of the deformable structure ͓48-50͔. Damping can be considered using a viscous model when the load is small and the temperature high. Apparently, the addition of internal friction does not have any notable influence ͓28͔.

Rotor Model. Classical theoretical studies have been interested in the dynamic characteristics considering the simplified case of a rigid rotor supported symmetrically by two identical bearings. The modeling of the dynamic action of the bearings is performed using the linear approach. Few authors considered the nonsymmetrical problem with a deformable rotor, and no existing study presents the direct coupling ͑with nonlinear analysis of the bearings͒ between the dynamic analysis of the rotor and its interaction with the foil bearings.

Comparison Between Linear and Nonlinear Methods

It is interesting to compare the linear and nonlinear methods to determine the range of validity. To this end, results in terms of stability level and unbalance response are given. For simplification, a rigid rotor symmetrically mounted on two identical three lobes rigid bearing is considered. The geometry of the lobes is circular and no geometrical preload is applied. The main characteristics are given in Table 1. The load on the bearing is only that of the rotor mass. For the stability study ͑Fig. 2͒, the critical speed is investigated for different rotor mass values. For the unbalance study ͑Fig. 3͒, the rotor mass is fixed at 2.25 kg and the mass unbalance is fixed at 1 g mm.

The results presented in Figs. 2 and3, thus confirm the nonlinear dynamic behavior of the air bearings according to the operating range. Concerning stability, at high loads ͑corresponding to high eccentricity͒ the stability threshold given by the linear approach differs from that given by a nonlinear approach. Concerning response to a mass unbalance, there is a transition between the ranges where the two methods are in agreement and where they present strong differences. At high velocity ͑low values of eccen-tricity͒, the two methods give identical results, but they present dramatically different results at low velocity ͑high values of ec-centricity͒.

However, the differences in phenomena observed concerning the domains of validity of the two approaches remain to be clarified. The definition of the limits of validity of the linear theory is thus a relevant direction of study for aerodynamic bearings. The comparison of the results given by linear and nonlinear approaches, Figs. 2 and3, highlights the differences between these two methods. The curves obtained are very distinct. These differ- the stiffness varies͒ can be adapted according to the desired performance. The first bearing developed with this technology consisted of a rigid sleeve onto which a corrugated foil and an upper foil are welded extending along all the bearing circumference. This foil assembly determines the general profile of the bearing film and its deformability ͑Fig. 1͒.

A correct implementation of the fixed deformable structure of the bearing foil must consider a variation of its structural stiffness with the angular position of the considered point ͓18͔. The characteristics of this type of bearing are very close to its rigid counterpart ͑a smooth journal bearing͒ and this allows the bearing to support high loads,͓9,19͔.

Aerodynamic Performance. The literature is replete with numerical and experimental work presenting the static and dynamic performance of foil bearings ͓19-30͔. The description of their behavior follows the same format as for hydrodynamic bearings. However, in addition, one has to take into account the effect of compressibility of the fluid coupled with the structural deformation, i.e., an aero-elastic coupling.

The analysis of the foil bearing dynamic behavior can be carried out according to two different methods. In the linear method, the dynamic coefficients ͑stiffness and damping matrices͒ of the bearing are first calculated. Reynolds equation is solved for small increments of relative journal and bearing displacement and velocity, enabling evaluation of these coefficients. These matrices are then used to study the rotor dynamics. ͓31-34͔. With the nonlinear method, the rotor dynamic movements are calculated in a step by step fashion based on the coupled resolution of Reynolds equation. This process gives the forces due to the bearings, and the equations of rotor dynamics ͓35-38͔.

The assumption of small displacements of the shaft center inside the bearing causes the linear analysis to have certain limitations. The method assumes that the variation of aerodynamic force with respect to its equilibrium position is nearly linear. A nonlinear analysis calculates the bearing forces at each time step to study the rotor dynamic behavior. No assumptions are needed for the bearing force evolution ͓30,39,40͔. This type of study is rarely presented in the literature of foil bearings because of large computer demands in terms of time and memory ͓41,42͔. On the other hand, similar work with incompressible fluids ͓35,31,43͔, confirms the importance of considering this type of analysis in certain operating conditions. Indeed, some isolated works highlight interesting phenomena on the specifics of the linear and nonlinear approaches ͓44-46͔. According to these works, the linear approach presents a restricted field of validity. It is thus important to investigate this point in order to avoid erroneous results.

Bearing Damping. The phenomena related to damping within the bearing structure are the subject of many numerical and experimental studies. One of the assumptions used in the past to ences occur both in term of stability threshold and response to unbalance. When aerodynamic bearings are highly loaded, it seems to be more accurate to use the nonlinear analysis.

Accounting for Structural Damping

Following the conclusion of the study where we compared the linear and nonlinear approaches, we set up a model taking into account a structural damping similar to that presented for the linear analysis in ͓51͔ and ͓33͔, but adapted to nonlinear calculation ͓26͔.

Modeling. This model has been detailed in ͓52͔. The main idea here consists in rewriting the balance equations for deformable bearings taking into account the speed of deformation. The compressive forces generated in the air film are balanced by forces due to the bearing stiffness ͑or the compliance S c ͒ and due to its damping ͑viscous damping coefficient per unit of area noted C e ͒, Fig. 4.

Starting from the balance forces on the upper foil the expression of the variation of the pressure field between t and t + dt can be written as

͑P t+dt -P t ͒ = 1 S c ͑H t+dt -H t ͒ + C e ͩͯ dH dt ͯ t→t+dt -ͯ dH dt ͯ t-dt→t ͪ ͑1͒
Expressions of the terms dH / dt described previously are as follows:

ͯ dH dt ͯ t→t+dt = H t+dt -H save t+dt ⌬t and ͯ dH dt ͯ t-dt→t = H t -H save t ⌬t ͑2͒
where H save t is the film thickness at time t before deformation, H t is the film thickness at time t after deformation, H save t+dt = H t + displacement of the shaft between t and t + dt ͑the film thickness at time t + dt before deformation͒, and H t+dt is the film thickness at time t + dt after deformation. From the above, ͑1͒ can be expressed in the following form:

͑P t+dt -P t ͒ = 1 S c ͑H t+dt -H t ͒ + C e ͩ H t+dt -H save t+dt ⌬t - H t -H save t ⌬t ͪ ͑3͒
Therefore, the relation giving the film thickness can be obtained as

H t+dt = H t + S . ⌬t ⌬t + S c . C e ͩP t+dt -P t + C e ⌬t ͑H save t+dt -H save t ͒ͪ ͑4͒
This film thickness is then introduced in the Reynolds equation in order to obtain the pressure field at each time step,

ץ ץ ͫ PH 3 ץP ץ ͬ t+⌬t + ץ ץZ ͫ PH 3 ץP ץZ ͬ t+⌬t = ⌳ͩ ץPH ץ ͪ t+⌬t + 2⌳ͩ P t+⌬t -P t ⌬t H t+⌬t + P t+⌬t H t+⌬t -H t ⌬t ͪ ͑5͒
Influence of Damping on the Threshold of Stability. The structural damping coefficient influence on the stability threshold is then studied using the set up model. The viscous damping coefficient used in the model can be replaced, in an equivalent way, by the damping coefficient used classically in rotor dynamics ͓53,54͔,

␣ = C e 2 ͱ KM
where: C replaces the viscous damping coefficient C e , K replaces the reverse of the compliance coefficient S c , and M represents the mass of the rotor assigned to the bearing. All the studied profiles of the bearings are with three deformable smooth lobes ͑Fig. 5͒ and all the grooves are identical and have an extent of 10 deg ͑Table 2͒. We thus have represented the critical engine failure speed diagram for various values of ␣ associated with this profile, Fig. 6. For comparison, we include the curve corresponding to the rigid profile.

It appears clearly that the addition of damping within the deformable structure, brings about a striking increase of the stability level for a compliant bearing ͑as for the low values of ␣͒. As could be expected, the higher the damping coefficient, the higher the critical speed. Furthermore, the critical speed does not decrease further with load increase. The relation between the critical Comments. In standard approaches the dynamic behavior of a deformable shaft supported by hydrodynamic bearings or ball bearings, is evaluated with the assumption that the bearing action is linear. The stiffness and damping matrix of the bearings are calculated and then introduced into a rotor dynamic computer code. This type of coupling is described as "weak." In the case of the aerodynamic bearings, we have been able to demonstrate that a linear dynamic analysis can lead to erroneous results according to the considered operation range.

The nonlinear character of a journal aerodynamic bearing is thus a significant fact to be considered in order to carry out an analysis of the rotor dynamic behavior. Therefore, we propose a "strong" coupling whereby the interaction of the bearings with the shaft is performed directly, step by step in time, in only one mode of the total calculation. In this first part the rotor was rigid and symmetrically supported by two identical aerodynamic bearings.

It is now interesting to study the behavior of a more realistic case where the rotor is deformable and the bearings are not symmetrically mounted on the rotor. Because the approach presented is new, we have gradually proceeded towards the more complete problem of a deformable shaft with a nonlinear action of the journal aerodynamic bearings. Thus, the first part of the analysis consists of a dynamic model of a simple deformable rotor resulting from ͓14͔. The Rayleigh-Ritz method with a three mode modal base is employed.

The dynamic behavior of a rotor mounted on two ball bearings was then compared to existing results for validation. Subsequently, the model was used to study the behavior of a deformable shaft running on two air bearings. The comparison with a linear "traditional" dynamic modeling confirmed once again the unique nonlinear dynamic behavior of the system. Finally, a comparison between a ball bearing and an air bearing points out the importance of the bearings on the overall dynamic behavior.

Rotor Dynamic Study With Foil Bearings

Short Description of the Rayleigh-Ritz Model. The Rayleigh-Ritz model uses a modal base restricted to three simple modes: pumping ͑or translation͒, pitching ͑or swinging͒, and first bending, Fig. 7. This modeling constitutes a first approach taking into account the two rigid body modes and the first bending mode of the shaft ͓54͔.

The symbol R 0 ͑X , Y , Z͒ indicates the initial frame ͑Fig. 8͒, in which the rotor axis at rest lies along the Y direction. The rotational speed ⍀ is considered constant. The rotor is considered as a beam and the neutral fiber displacement is studied. Each point on this fiber has two degrees of freedom that are displacements in the X and Z directions.

The expressions of the shaft neutral fiber displacements in the X and Z directions are, respectively, written in the following form:

u͑y,t͒ = f a ͑y͒q 1 ͑t͒ + f b ͑y͒q 3 ͑t͒ + f c ͑y͒q 5 ͑t͒ = f a ͑y͒q 1 + f b ͑y͒q 3 + f c ͑y͒q 5 ͑6͒ w͑y,t͒ = f a ͑y͒q 2 ͑t͒ + f b ͑y͒q 4 ͑t͒ + f c ͑y͒q 6 ͑t͒ = f a ͑y͒q 2 + f b ͑y͒q 4 + f c ͑y͒q 6 ͑7͒
where q i are the generalized independent coordinates. The displacement functions, f a , f b , and f c are, respectively, selected to represent the geometry of the vibration modes: pumping, pitching, and bending of a beam of constant section supported at its ends, = -g a ͑y͒q 1 -g b ͑y͒q 3 -g c ͑y͒q 5 ͑12͒

The second order derivatives for u and w displacements are necessary to express the elastic energy of the shaft,

ץ 2 u ץy 2 = d 2 f a ͑y͒ dy 2 q 1 + d 2 f b ͑y͒ dy 2 q 3 + d 2 f c ͑y͒ dy 2 q 5 = h a ͑y͒q 1 + h b ͑y͒q 3 + h c ͑y͒q 5 ͑13͒ ץ 2 w ץy 2 = d 2 f a ͑y͒ dy 2 q 2 + d 2 f b ͑y͒ dy 2 q 4 + d 2 f c ͑y͒ dy 2 q 6 = h a ͑y͒q 2 + h b ͑y͒q 4 + h c ͑y͒q 6 ͑14͒
One may observe from the above expressions that g a ͑y͒ = 0 and h a ͑y͒ = 0 because f a is of zero order. In the same way, h b ͑y͒ =0 because f b is of order one.

In the calculations below, we will not entirely develop the equations because manipulations are extensive but straightforward. We point out the expressions of various energies which should be calculated to write the Lagrange equations system ͓53͔.

Disk Element. The disk is considered as symmetrical. Its kinetic energy T D can be written in the following way:

T D = 1 2 M D ͑u ˙2 + w ˙2͒ + 1 2 I Dx ͑ ˙2 + ˙2͒ + 1 2 I Dy ͑⍀ 2 + 2⍀ ˙͒ ͑15͒
The energy of deformation of the disk is nil since it is considered as being rigid.

Rotor Element. The rotor elements are considered to be axisymmetric. The expression of kinetic energy T S is as follows:

T S = S 2 ͵ 0 L ͑u ˙2 + w ˙2͒dy + I 2 ͵ 0 L ͑ ˙2 + ˙2͒dy + IL⍀ 2 + 2I⍀ ͵ 0 L ˙dy ͑16͒
As the term IL⍀ 2 is constant, we do not have to introduce this term in the writing of the Lagrange equations. The energy of deformation of the shaft, U S is written as

U S = EI 2 ͵ 0 L ͫͩ ץ 2 u ץy 2 ͪ 2 + ͩ ץ 2 v ץy 2 ͪ 2 ͬ dy ͑17͒ U S = EI 2 ͵ 0 L h c 2 ͑y͒dy͑q 5 2 + q 6 2 ͒ ͑ 18͒
Unbalance. The shaft can be unbalanced. This unbalance is defined by a mass m u at a distance d to the geometric center of the shaft. Its kinetic energy is written as

T u = m u d⍀͑u ˙cos ⍀t -w ˙sin ⍀t͒ ͑ 19͒
Having presented the energies, it is now appropriate to look at the actions of the ball bearings and/or the air bearings. In addition, to represent the operation of a shaft on air bearings, the influence of gravity must be added to the equations.

Ball Bearings. The shaft can possibly be supported by ball bearings. Their stiffness and damping coefficients are assumed to be known. The influence of the deflection is neglected, Fig. 9. The expression of the forces' virtual work exerted on the shaft is written in the following way:

␦W rlt = -k xx u␦u -k xz w␦u -k zz w␦w -k zx u␦w -c xx u ˙␦u -c xz w ˙␦u -c zz w ˙␦w -c zx u ˙␦w ͑20͒ or ␦W rlt = F rltx ␦u + F rltz ␦w ͑21͒
where F rltx and F rltz are the components of the generalized forces of the ball bearing. In matrix form it can be written as

ͫ F rltx F rltz ͬ=-ͫ k xx k xz k zx k zz ͬͫ u w ͬ-ͫ c xx c xz c zx c zz ͬͫ u ẇ ˙ͬ ͑22͒
Air Bearings. Taking the air bearings into account in the behavior of a shaft requires the expression of the virtual work due to these bearings. As presented previously, the forces generated by the air bearings on the shaft depend at the same time on the local deflection of the rotor on the level of each bearing as well as the rotational speed. The expression of the virtual work of the forces generated by an air bearing is as follows:

␦W pal = F palx ␦u + F palz ␦w ͑23͒ with F palx , F palz , the components of the forces generated by the air bearing on the rotor. The components of the forces generated by the bearing are calculated at each time step, using Reynolds equation expressed in dynamic regime ͑Eq. ͑5͒͒.

The Reynolds equation provides the field pressure which leads, upon integration, to the fluid film forces acting in the bearing. This kind of coupling is detailed in ͓26͔ in the case of a rigid rotor As the functions h a and h b are nil, the expression of the energy of deformation of the shaft is written according to the generalized coordinates, symmetrically supported by two air bearings. In order to simplify calculations, we do not consider misalignment of the bearing generated by pitching, bending of the shaft or a difference in value of the magnitude of the static loads applied to the bearings.

Gravity. The action of gravity creates an offset of the rotor in the air bearings. We consider that the total mass of the shaft exerts a force located at the shaft center of gravity. The expression of the virtual work due to gravity is as follows:

␦W gravity = -M total g␦w͑y cdg ͒ ͑ 24͒ with M total as the total mass of the shaft, g is the acceleration of gravity, and y cdg is the Y-coordinate of the shaft's center of gravity.

We obtain the system of equations by application of the Lagrange formalism,

d dt ͩ ץT c ץq ˙i ͪ -ץT c ץq i + ץU ץq i = ␦W i ͑25͒
with i =1-6, T c being the kinetic energy, U the energy deformation of the whole shaft components, ␦W the total virtual work, and q i the generalized coordinates of the system. Thus we can obtain the system of equations for the rotor dynamic which can be written in the following matrix form,

͓M͔͕q ¨͖ + ͓C͔͕q ˙͖ + ͓K͔͕q͖ = ͕F͖ ͑ 26͒
with ͓M͔ the matrix of mass, ͓C͔ the matrix governing the gyroscopic effect and possibly the ball bearing damping, ͓K͔ the matrix of the rotor stiffness and possibly the ball bearing stiffness, ͕F͖ the vector of the external actions ͑the unbalance, the air bearings forces and gravity͒.

Dimensionless Equations. The system of equations can be written in the following form:

͓M ¯͔͕Q ¨͖ + 1 ⍀ ͓C ¯͔͕Q ˙͖ + 1 L 3 ⍀ 2 ͓K ¯͔͕Q ¯͖ = ͕F ¯͖ ͑ 27͒
with the following dimensionless variable definitions:

Q ¯i = q i C 0 , Q ˙i = q ˙i C 0 ⍀ , Q ¨i = q ¨i C 0 ⍀ 2 , ͓M͔, M ij = M ij L 3 , C ij = C ij L 3 , F ¯i = F i L 3 C 0 ⍀ 2 ͑28͒
simple. Generally, the rotors used in oil-free turbomachinery are quite stiff, thus considering only rigid and first bending modes is a reasonable assumption.

Results and Discussion

The Rayleigh-Ritz model will now be used to study the dynamic behavior of a deformable shaft supported by air bearings. Thus we will be able to show the usefulness of this type of analysis while observing, in particular, the influence of the rotor flexibility on the response to an unbalance and the stability threshold of the system. We consider a solid or hollow shaft of external length 300 mm and 50 mm of diameter, as shown in Fig. 12. In the speed range of 2000-18,000 rpm, the dynamic behavior is close to the symmetrical rigid configuration ͑Fig. 13͒. This rotor seems too short and thus too stiff to observe deflection phenomena in the considered speed range. Once again this validates the Rayleigh-Ritz model. The influence of the bending mode has not been shown by this example.

To increase the bending mode, a rotor with 900 mm length is chosen for the rest of the study. The air bearings, as in the preceding case, are located at 50 mm from the ends of the rotor. The dynamic calculation of this shaft's behavior with nonlinear action of the air bearings gives an interesting response curve to unbalance. We always observe a peak of amplitude towards 4000 rpm. It is similar to the case of a rigid rotor but a second peak appears towards 9500 rpm. We compare this response curve in Fig. 14 with that obtained from the symmetrical rigid rotor simple model in order to highlight the effect of the rotor's flexibility.

It is interesting to represent the deformation of the rotor for speeds of 4000 and 9500 rpm, Figs. 15 and 16. These figures give the representation of the trajectory of the rotor's center for various Y-coordinates. We took the precaution to remove the static deformation of the shaft due to gravity, so that the observation of the deformation in rotation is easier.

It appears clearly that the mode of vibration at the speed of 4000 rpm is of the "rigid body" type because we obtain a very small twist of the shaft. The amplitude of displacements along the Z-axis is one third of those along the X-axis, Fig. 15.

For the speed of 9500 rpm, the associated mode of vibration introduces the inflection of the rotor whose form is complex. It should be noted that the amplitude of the deformation along the Z-axis ͑direction of the weight͒ is approximately four times higher than that in the X-direction, Fig. 16.

One can also plot on the same graph the response curves to unbalance in configuration ball bearings and air bearings, Fig. 17. It appears that the two curves present notable differences. First of all, we note that the air damping of the bearings leads to a decrease of the amplitude of the rotor center trajectory compared that of a "ball bearing" configuration ͑with identical unbalance͒. At high speed, the eccentricity is low, film thickness is high, and the air bearing stiffness is less than that of a classical ball bearing. Hence, a part of the second mode energy is in the bearing and not in the rotor, as with ball bearings. In the air bearing configuration, we thus obtain displacements which are due to both the rigid and with time T=t⍀, ⍀ as the rotational speed. The parameter C 0 denotes an arbitrary value of the order of the bearing radial clearance. The symbols␦and L, respectively, denote the density and the length of the rotor.

This dynamic system must be solved step by step in time to take into account the nonlinear variations of aerodynamic bearing forces in the right hand side member. The numerical integration scheme is the same as that presented in ͓26͔.

Validation

The results from the model are compared to those given in ͓14͔ for a rotor described in Fig. 10, supported by two ball bearings ͑Table 3͒. In the reference, the rotor is modeled by finite element method with a pseudo-modal reduction.

Four critical speeds are given by the Rayleight Ritz model ͑Table 4͒. The first two critical speeds are very close to the reference value, Fig. 11. The last two critical speeds are quite different. In this case, modes three and four involve bending modes that are more complex than the first bending mode involved in the Rayleigh-Ritz model. The first two modes are the first bending mode. This study confirms that the simplified Rayleigh-Ritz model gives good results if the rotor bending mode remains bending modes. At low speed, a rigid vibration mode is present with the air bearing configuration and not with that of the ball bearings. This simple case outlines how bearings influence the whole dynamic response of the system.

General Conclusions

This study of aerodynamic journal bearing dynamic characteristics has been conducted because it is crucial to understanding rotor dynamic behavior. Initially, the rotor is considered to be rigid and supported symmetrically by two identical radial bearings. We have shown that the linear modeling, i.e., representing the bearings with stiffness and damping coefficients, has a re-stricted field of validity. For instance, in the case of a three lobe rigid bearing without geometrical preload, this range is limited to relative eccentricities lower than 0.65.

Thus, the prediction of the bearing stability threshold and the response to dynamic excitations can be erroneous if one uses a linear approach beyond the appropriate limit for this type of aerodynamic bearing. In addition, we illustrated the fact that for high loads the deformable bearing profiles have a lower stability than their rigid counterparts. Since we were fully aware that this observation goes against some physical phenomena observed in experiments, we set up a model which takes into account a structural damping of viscous type in the nonlinear approach. Even with very low coefficient values of the bearing structural damping, the level of stability could be increased, in accordance with the phenomena presented in the literature.

Having clearly highlighted the nonlinear character of aerodynamic bearings dynamic behavior, we considered direct coupling with a deformable shaft. At this point we developed a calculation based on the Rayleigh-Ritz method using a base of three modes ͑pumping, pitching, and bending͒. Thus, we were able to confirm the coherence of the obtained results compared to those from the rigid rotor model. The deformation of the rotor reveals modes of vibrations which cannot be predicted by the rigid model. At low speed, a rigid mode of vibration is present with the air bearing configuration. At high speed we obtain displacements which are due to both rigid and bending modes. This simple case brought into light the influence of the bearing on the dynamic response of the complete system. It appears clearly in the first part of this paper that the addition of damping within the deformable structure leads to the increase of the deformable bearing's stability level.

As could be expected, the higher the damping coefficient, the higher the critical speed. Furthermore, the critical speed does not decrease further with load increase. The relation between the critical engine failure speed and the static load seems to deviate from linear behavior and approach a higher order polynomial form. It is clear that future work will relate on the one hand to the modeling of dry friction between the foils and the bumps to properly quantify the damping brought by this process, and on the other hand to the development of a more realistic rotor dynamic model so as to properly characterize the dynamic behavior of the whole system.

Fig. 2

 2 Fig. 2 Evolution of critical speed value against the dimensionless static load

Fig

  Fig. 5 Schematic of three lobed deformable bearing

Fig. 9

 9 Fig. 9 Stiffness and damping coefficients of a ball bearing

Fig. 10

 10 Fig. 10 Configuration of the reference case

Fig. 11

 11 Fig. 11 Dynamic response at y = l 2

NomenclatureC

  ϭ radial clearance C ij ϭ damping coefficients D ϭ bearing diameter e ϭ bearing eccentricity E ϭ Young modulus f a,b,c ϭ displacement functions F rltx , F rltz ϭ roller bearing forces F palx , F palz ϭ bearing forces g ϭ gravity acceleration h ϭ film thickness H ϭ dimensionless film thickness, H = h / C I ϭ rotor inertia in the X direction I Dx ϭ disk inertia in the X direction I Dy ϭ disk inertia in the Y direction K ij ϭ bearing stiffness coefficients K s ϭ stiffness coefficient of the structure L ϭ bearing length M ϭ rotor mass M D ϭ disk mass O ϭ bearing center O a ϭ shaft center p ϭ pressure p a ϭ ambient pressure P ϭ dimensionless pressure, p / p a q i ϭ generalized coordinates R ϭ bearing radius R a ϭ shaft radius S c ϭ compliance S ϭ bearing section S ¯c ϭ dimensionless compliance S ¯= S͑P a / C͒ t ϭ time T ϭ dimensionless time, T = t T S ϭ rotor kinetic energy T D ϭ disk kinetic energy T C ϭ total kinetic energy T u ϭ unbalance kinetic energy U 1 ϭ speed component in X direction U S ϭ rotor energy deformation u͑y , t͒, v, w ϭ shaft displacement in the X, Y, and Z directions X, Y, Z ϭ initial frame x, y, z ϭ local frame W o ϭ static load W ϭ carrying capacity W ¯ϭ dimensionless carrying capacity, W ¯= W / ͑p a LD͒ W dx,y ϭ external dynamic forces ␦W ϭ virtual work ␣ ϭ damping coefficient ⍀, ϭ rotational speed ⌰, ⌿ ϭ shaft angular displacements ⌳ ϭ bearing number

Fig. 15

 15 Fig. 15 Deformation of the rotor at 4000 rpm

Table 1 Data input for the lobed bearing

 1

	L / D	0.75, 1, 1.25
	D	50 mm
	Number of pad	3
	Pad angle, deg	110
	Radial clearance,

m 50 Fig. 1 Bump foil bearing description

Table 3 Data for the reference case

 3 / m, K zz =7ϫ 10 7 N / m, K xz = K zx =0 C xx =5ϫ 10 2 N / m / s, C zz =7ϫ 10 2 N / m / s, C xz = C zx =0

		l 1 =0.2 m, l 2 = 0.5 m, l 3 =1 m, L = 1.3 m
		Plain section, = 0.05 m
	Rotor	E =2ϫ 10 11 N / m 2 , = 7800 kg/ m 3 , = 0.3
			Inner diam,	Outer diam,
	Disk	Thickness, m	m	m
	0.05 0.05 0.06 K xx =5ϫ 10 7 N Unbalance Disk 1 ͑l 1 ͒ Disk 2 ͑l 2 ͒ Disk 3 ͑l 3 ͒ Bearing stiffness and damping	0.05 0.05 0.05 200 g mm on disk 2	0.12 0.2 0.2
	Velocity	from 0 to 30000 rpm

Table 4 Critical speeds

 4

		Critical speed,	Critical speed,
		Rayleigh-Ritz	Reference
		rpm	rpm
	1	3610	3620
	2	3840	3798
	3	12850	10018
	4	15080	11279
	5 6 7	¯16785 ¯24408 ¯26615