Martin Poncelet
email: poncelet@lmt.ens-cachan.fr

Hugo Leclerc
email: leclerc@lmt.ens-cachan.fr

A Digital Image Correlation algorithm with light reflection compensation

Keywords: Full field measurement, Digital Image Correlation, Light reflection, Speckle, Uncertainty assessment

This work focuses on the influence of light reflection on Digital Images Correlation results at the macroscopic scale, and on a way to circumvent this problem. It shows that the local displacement uncertainty rises up to 5 times the usual one when a reflection occurs. An observation of the topography of the speckle reveals an important sub-pixel roughness that explains the grey level fluctuations at the pixel scale, spoiling the calculation the gradient of the texture. To circumvent this problem, a new DIC algorithm is proposed, based on a single minimization with several pairs of images where the reflections are located in different regions. For each image, weighting functions are used to 'exclude' the reflection regions from the calculation, while the necessary information is obtained from the other images.

Introduction

The increasing performances of optical full-field kinematic measurement techniques has lead to their widespread use. Among the many possible ones [START_REF] Grédiac | Full-Field Measurements and Identification in Solid Mechanics[END_REF], such as photoelasticity, moiré, holographic and speckle interferometry, grid method and Digital Image Correlation (DIC), this last one [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF] is very often chosen for its 'simple' use conditions in comparison with the other techniques. Consequently, DIC measurements are now performed out of the research laboratories [START_REF] Yoneyama | Bridge deflection measurement using digital image correlation[END_REF], where the testing conditions are no longer optimal, or performed on industrial products which surfaces are no longer flat [START_REF] Berring | Torsional performance of large wind turbine blades experimental and numerical analysis[END_REF][START_REF] Beaubier | CAD-Based Calibration and Shape Measurement with StereoDIC -Principle and Application on Test and Industrial Parts[END_REF]. With these new configurations, lightning problems appear: non-constant illumination, cast, coreshadows and reflections.

Non-constant global illumination issues may be prevented by cunning protocols: Yoneyama et al. [START_REF] Yoneyama | Bridge deflection measurement using digital image correlation[END_REF] performs tests on a bridge at night with artificial illumination, Grant et al. [START_REF] Grant | High-temperature strain field measurement using digital image correlation[END_REF] uses a combination of illumination and optical filter at the same wavelength to prevent disturbances.

Large scale heterogeneity of illumination issues due to reflection can be partially solved during measure by use of polarized filters [START_REF] Kim | Characterization of the postnecking strain hardening behavior using the virtual fields method[END_REF]. However every reflection can not be withdraw with this technique, so that post-processing techniques are necessary.

One may want to 'enhance' images with several existing techniques, virtually lighting up shadows or correcting burned areas. Based on empirical or scientific models, some techniques aim aesthetic or visibility purposes [START_REF] Sedgewick | Scientific Imaging with Photoshop: Methods, Measurement, and Output[END_REF][START_REF] Babaloukas | Evaluation of three methods for retrospective correction of vignetting on medical microscopy images utilizing two open source software tools[END_REF], while others aim quantitative purposes, like automatic detection [START_REF] Gómez-Sanchis | Correction of the effects of the light source on spherical objects. An application to the analysis of hyperspectral images of citrus fruits[END_REF]. In any case, these visual enhancements are not of interest from a DIC point of view, potentially introducing spurious grey level gradients that spoil the results. Consequently, some full-field measurement algorithms consider an optical flow conservation relaxation. Local DIC algorithms use normalized cross-correlation calculation for each sub-image for example [START_REF] Bornert | Assessment of digital image correlation measurement errors: methodol-ogy and results[END_REF][START_REF] Sutton | Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[END_REF], while global DIC algorithms use a local grey level correction at the scale of the analyzed elements [START_REF] Besnard | Analysis of Image Series through Digital Image Correlation[END_REF].

To go further with this idea, it is possible to use the variation of the optical flow as a pertinent information in itself, for example to measure both kinematic and thermal fields from IR images [START_REF] Maynadier | One-Shot measurement of thermal and Kinematic fields: InfraRed Image Correlation (IRIC)[END_REF].

Dealing with shadows is more difficult. The best theoretical solution to light up shadowy areas is the use a homogeneous light source like dome lights, but it is experimentally impossible in many cases. Several light sources are then often preferred, but it multiplies the number of unwanted reflections. Reflections are thus nearly unavoidable in the case of complex shapes.

Because of the strong variations due to reflections, images are often characterized by very low grey levels, near saturation or saturation situations. The saturated areas do not carry information, so that the measurement is strictly impossible. A simple way to prevent this is to use a slightly higher dynamic range (for example a semi-professionnal Digital Single Lens Reflex camera with up to 14 bit depth), or High Dynamic Range devices (scientific grade cameras with more than 16 bit depth and very low noise) during the experiment.

The remaining important issue is that during the translation of the object, the texture of the surface follows it while the shadows and reflections do not exactly (Fig. 1). During a rotation the texture may nearly not move from the camera point of view near the rotation axis, while reflections change tremendously. Shadows may also evolve, usually only slightly but non-convex or partially flat objects can present sudden form and cast shadows. If an adequate lightning setup is used, the difficulties linked to shadows may be suppressed as stated before, but reflections still exists and are generally characterized by high grey level gradients. Consequently their gradients may be 'misinterpreted' as a texture displacement. The simplest way to circumvent this problem afterward is to exclude the problematic areas from the measurement zone. The information lost in the disturbed areas may be important, all the more if stereovision techniques (like stereo-DIC [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[END_REF][START_REF] Besnard | Analysis of necking in high speed experiments by stereocorrelation[END_REF]) are used, because a reflection may affect a different zone on each picture, contrary to shadows.

A second problem is related to the fluctuations of the apparent texture, due to changes of the angles between the lights, the cameras and the surfaces. They may lead to local additional uncertainty which should be taken into account by a local weighting.

Up to now, most full-field kinematic measurements [START_REF] Sutton | Advances in two-dimensional and three-dimensional computer vision[END_REF] use an homogeneous weight, assuming a constant apparent noise.

The present paper proposes to study the influence of reflections on DIC measurements and then to develop a 'reflection-proof' DIC technique. One focuses on the case of a DIC measurement at the 'macroscopic scale' with a regular, white-light, as in many industrial and academic applications. The study uses regular speckles obtained with mat sprays and no satured areas to quantify the measurement quality loss. The results must therefore be considered as obtained 'in the best measurement conditions' (no gloss paint, no saturated areas). It will be shown that it is not possible to correct these reflections by numerically taking into account the geometry of the specimen because of the inner characteristics of the coating, which induce unpredictable reflections at the pixel scale.

The proposed 'reflection-proof' technique requires a simple experimental setup and is adapted to non-satured cases as well as satured ones. It is based on the use of a quick succession of different light sources: several different light sources sequentially illuminate the surface during a negligible short time in comparison with the object velocity. This technique enables to exploit any area of the studied object since the reflections appear at different locations on the images while the mechanical configuration is the same. A dedicated algorithm solves a global problem, i.e. takes into account all the images corresponding to the different light sources, with an appropriate local weighting of each image. A validation of this technique will be last shown in the case of a sphere illuminated by 2 light sources, and observed by a single camera.

The paper is divided as follow. In a first part, the influence of the reflection on DIC measures is studied: a setup is proposed to study the influence, the results are presented for a standard algorithm and one with a local grey level correction, and then the reasons of this influence are searched with a measurement of the evolution of the grey level with the angle of the incident light, and a sub-pixel scale roughness study of the coating. In a second part, the reflection-proof technique, called 'Multiple Illumination Digital Image Correction' is explained and mathematically formulated. A third part presents the validation of the technique, and a discussion follows.

2 Light reflection on a typical paint speckle

Experimental setup and protocol

To characterize the reflection of light on a painted speckle, the optical bench presented Fig. 2 and 3 is used. All the devices and paints used thereafter are similar to the ones generally used for DIC measurement. The only special feature of the bench is the possibility to accurately position every devices. Its principe is to observe an object with a camera, while the angle of the object αc (and so the angle of the incident light α i used to illuminate the object) is changed by use of rotation stages. The camera is a Canon EOS 60D Digital Single Lens Reflex equipped with 60 mm focal lens. The image has a resolution equal to 5184×3456 pixels, and is saved in a lossless format (.raw). The white light source (EFFISHARP NearField 000-02 LED lamp, Effilux, Orsay (France)) is considered as point-wise and homogeneous. The observed objects are respectively a 50x75-mm polished aluminum plate and a 82.5-mm steel sphere covered with a white layer of mate paint and then a spray of black mate paint (RAL 9010 'Noir pur mat', RAL 9010 'Blanc pur mat', Ront Production, France). The optical axis of the camera is in the same horizontal plane as the centre of the sphere, while lightning b 1 (reps. b 2 ) is slightly above (resp. below).

Light reflection influence on kinematic measurement

To study the reflection influence on the DIC measure, the sphere is placed on the rotation stage with αc = 0 at the beginning of the test. In this configuration, the rotation axis is considered on the point of the sphere closest to the camera (Fig. 3).

Light source b 1 is on, and αc + α i is constant during the test (i.e. the light is fixed with respect to the camera). The aperture is so that the observed surface is in the depth of field. The exposure time is so that no pixel is saturated (the difference between the maximum grey level value and the saturation value (65535) is around 60, i.e. around a thousandth of the dynamic range). Images are taken for different αc, ranging from 0 to 14 • . Five 50 × 100 pixels areas noted A0 to A4 are chosen along the horizontal axis X (Fig. 4). The reflection displacement is clearly visible between the initial position in area A4 (Fig. 4 (a)) and the last one in area A2 (Fig. 4 (b)). The value of the grey level (in percent of the total dynamic range) along the horizontal line crossing the disturbed area f line (x) is plotted to better show the displacement of the reflection.

DIC measurement is first done with a TRI3 'finite element' like algorithm [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF][START_REF] Leclerc | Integrated Digital Image Correlation for the Identification of Mechanical Properties[END_REF]. The reference image is the one at αc = 0, while deformed images correspond to the different other αc values. The mesh size, as for the next algorithm results, is 25 pixels. The residual is calculated for each deformed image as the standard deviation of the difference between the reference image and the deformed image corrected by the found displacement using a bicubic interpolation for the sub-pixel displacement.

The average residual for each image is around 2% of the effective dynamic range (i.e.

maximum minus minimum grey level values of the images), which proves that the calculation has converged (usual mean residual without reflection is about 2%). the displacement is far more important, so that the previous noisy area is no longer directly visible. Since one may be interested in displacement heterogeneities with small magnitude embedded in a large quasi-uniform displacement (for example due to a rigid body motion), it is proposed thereafter to subtract the long wave components of the displacement field to better understand the reflection disturbance.

To quantitatively assess the perturbation of the displacement field due to the reflection, the 5 areas will be specifically studied. A high-pass filtered field Uy is calculated by subtracting to Uy a field obtained by applying a radial gaussian filter to Uy. The standard deviation of the gaussian filter is equal to 40 pixels, so that Uy does not contain the displacement corresponding to the rotation of the sphere and the integer locking artifacts ( [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF]), but only keeps the high frequency heterogeneity of the field in the studied areas. For each of the 5 areas, the mean value of Uy is plotted Fig. 6, showing a linear evolution with αc. The perturbation P is calculated as the standard deviation of Uy on each area. P is plotted as a function of the rotation angle on Fig. 7 (b). Three distinct behaviours are notable :

-For area A0, the reflection is far away for any image, so that it is not affected.

The value of the perturbation is very low (P ≈ 4 × 10 -3 pix) for small angle values (≤ 2 • ) and increases slowly with the rotation (up to P ≈ 4 × 10 -2 pix).

The perturbation is directly linked to the measurement uncertainty of the DIC technique.

-For area A4, the perturbation increases up to 4 times higher because the reflection on the reference image ( 0 • ) is over this area. This 'additional' perturbation exists even for very small rotation angles (< 1 • ) and is nearly constant once the reflexion is away from the area on the current image.

-For areas A1 to A3, the perturbation is nearly equal to A0 one for small rotation angles because the reflection is then not over these areas. Nonetheless for higher rotation angles, the reflection 'crosses' the different areas -For area A0, the perturbation is slightly lower than with the standard algorithm.

-For areas A4 and A3, the perturbation is also lower, but the important increase occurring for small αc values (≤ 2 • ) is nearly unchanged.

-For areas A1 to A2, the maximum value is clearly less than with the standard algorithm, but it is still higher than for area A0.

One concludes that the grey level correction is able to compensate for the global variation of illumination taking place during the large rotation of the sphere. However this linear grey level correction at the scale of the element is not adequate to delete the important increase of the perturbation for small rotation angles (i.e. the most common case during an experiment).

Grey level evolution with light orientation

The high effect of reflection on full field measurement results for large rotation angle was expected, since the grey levels evolve importantly. The smaller, but non negligible, effect for αc ≤ 2 • in the case of both DIC and GLC-DIC is less obvious because the reflection apparently does not change. To better understand the evolution of the grey level value in the disturbed area, the flat specimen is now used to study the evolution of the grey level with the angle of the incident light α i .

The specimen is placed on the rotation stage instead of the sphere (see dashed lines on Fig. 3), with the speckle plane on the rotation axis O. The normal axis On of the specimen is tilted with αc = 15 • . The light source b 1 is now in the same horizontal XZ plane as the optical axis. The angle α i is then increased step by step with its own rotation stage with a picture taken for every position. The grey level value for 40 pixels located along the rotation axis O is plotted as a function of α i (Fig. 9). For high values of α i (90 to 30 • ) the evolution is smooth and corresponds to diffuse reflection, with total extinction at 90 • . Two bundles of curves start from 90 • , the upper one corresponding to the white background of the speckle (high reflectivity) and the lower one to the black dots (low reflectivity). For lower α i values, the specular reflection contribution is visible. One expects a smooth increase with a maximum at α i = 15 • but the measure shows a very noisy evolution spreaded from roughly 30 to 0 • . Consequently, one understands that it is no use to try to correct the reflection on the images with the sole knowledge of the optical configuration (theoretical geometry of the objet, relative positions of the camera and light source). The speckle itself induces a non trivial grey level evolution with the incident light angle.

Sub-pixel surface study

To understand the reason of this erratic evolution of the grey level value with the reflection, the roughness of the speckle of the flat specimen is studied. To measure it, a confocal chromatic sensor (130 µm range, 0.016 µm uncertainty (standard deviation), CHR-150-L, STIL (France)) with a high precision translation stage (0.1 µm resolution, 1 µm flatness defect over 100×100 mm 2 , Micromesure 2, STIL (France)) is used. A 10×10 mm 2 area is measured on the flat specimen with a 5 µm pitch and a mean over 5 consecutive measures. Both the intensity I and the height h can be analyzed afterward.

As a first approximation, the variation of height h(x, y) (Fig. 10 (a)) is related to the presence of black paint droplets as seen on the intensity map I(x, y) (Fig. 10 (b)).This variation is less than 10 µm. The variation of the local slope arctan(|∇h|) on a 1×1 mm 2 region of Fig. 10 (a) is presented Fig. 11 (a). The uncertainty of the slope is evaluated to 0.23 • if one considers a 0.1 µm uncertainty of h measurement. The map is nonetheless very clear. The roughness is different on the white and black surfaces, probably due to the difference of filler or drying time (the black paint spray is dry in a few seconds while the white thick layer takes several minutes). In both case it reveals that the roughness is not random, but is characterized by a 'wrinkle' pattern with an average distance about 25 µm, i.e. at a lower (but not negligible) scale than the pixel physical size (55 µm, represented by black lines on Fig. 11). The roughness is rather sharp (up to 15 • slope). This explains the high intensity variations in the reflection region for a small rotation in the previous experiment. One deduces that this effect can not be taken into account with an a priori knowledge of the theoretical surface of the object, as done for example in [START_REF] Gómez-Sanchis | Correction of the effects of the light source on spherical objects. An application to the analysis of hyperspectral images of citrus fruits[END_REF], where the high frequency fluctuations are not of interest.

It can also not be withdrawn with techniques used for fringe projections -where the reflection has very low frequencies in comparison with the carrier frequency [START_REF] Takeda | Fourier transform profilometry for the automatic measurement of 3-D object shapes[END_REF] because its spectrum is close to the one of the speckle. In other cases -object and speckle characteristic length several orders of magnitude over the roughness of the paint speckle, other speckle creation method, different light wavelength or image acquisition method -the conclusion may be different. 3 Multiple Illuminations Digital Image Correlation

Principle

The principle of the proposed method consists in calculating the data needed for the solving steps of each DIC procedure for each lighting condition, with weights allowing to exclude the regions where the reflection lies. By a choice of multiple illumination (Fig. 12), the different weights allow some kind of complementary contribution. A single global minimization procedure is then launched for each displacement field update iteration.

Illumination 1 Illumination 2

Initial state images

Current state images

f 1 f 2 g 2 g 1 (Reflection)

Mathematical formulation

In a DIC algorithm, an initial image f is related to a deformed image g, the displacement field u, and the camera noise n by g(x + u) = f (x) + n(x) .

(

) 1
Solving the usual mono-illumination DIC problem, i.e. measuring u, consists in minimizing a functional Φ over a set of possible displacements v

Φ(v) = (g(x + v) -f (x)) 2 dx . (2)
In the present case, to take into account the fact that the same displacement is observed with different light sources (here 2), the functionnal is rewritten with two pairs of initial (f 1 , f 2 ) and deformed images (g 1 , g 2 ) as follow

Φ(v) = (g 1 (x + v) -f 1 (x)) 2 p 1 (x) dx + g 2 (x + v + i β i .ϕ i (x)) -f 2 (x) 2 p 2 (x) dx (3)
where β i are additional discrete unknowns, ϕ i are predefined rigid body motion functions, and p 1 and p 2 are two weighting functions.

If the p j functions are not used (or all of them set to a uniform value equal to 1), several pairs of images are used to assess the displacement field without respect of the existence of disturbed areas. Such a calculation leads to the results presented Fig. 13: the perturbation is nearly equal to the one obtained with one pair of images, only slightly lower. These functions are thus mandatory to delete the influence of the disturbed areas on both the reference image and the deformed image by taking into account the fact that under non-uniform lightning conditions (i.e. realistic conditions), some regions might not check the optical flow conservation. A first method to calculate the p j functions needs the knowledge of the geometry of the setup (specimen as well as lightnings and camera, and of course the motion of the specimen) and the bidirectional reflectance distribution function (BRDF) corresponding to the speckle. It is then easy to define the theoretical disturbed areas (roughly corresponding to the region between +/-15 centered on the specular reflection on Fig. 9) prior to the MIDIC calculation. Unfortunately the geometry is not always known (no CAD of the setup, no stereo DIC to measure the displacement of the surface of the specimen).

The BRDF is also rarely available. One may use the informations given in Fig. 9 to model the BRDF as a first approximation, but the measurement was restricted to one direction and only the incident angle has been probed. That would probably leading to a poorly relevant BRDF model.

A second method, more arbitrary but relevant, is to define the area that is disturbed (for example by using a criterion like the standard deviation of the displacement being over a given threshold of the usual standard deviation found on an unaffected area) after a first standard DIC calculation on the first pair of deformed images. Once the area has been assessed, one measures the corresponding average grey level on the perimeter of this area on the corresponding blurred images and uses it as a threshold to define the weighting functions for every images. This method has to be tuned for every different setup, but it is simple, fast and does not require to know the BDRF and geometry, so that it is presented thereafter.

We will only consider binary weighting functions, computed with the simple combination of a blur filter and a threshold. These weighting functions will merely suppress the regions where the light reflections may induce difficulties. For every angle value, one has first to blur the images taken with each light source to compute the weighting functions. A 2D isotropic filter is chosen, with a radius set to 256 pixels in the Fourier space, i.e. 1/4 of the image size. It thus only keeps the average grey level evolution due to the lightning, and not due to the texture itself. An example is given for the angle equal to 0 • 20 and the b1 lightning on Fig. 16, where (a) represents the image and (b) the blurred image.

A threshold has then to be chosen, over (resp. below) which the weighting function is set to 0 (resp. 1), so that it excludes the disturbed area of each image. To find the value of this threshold, one studies the perturbation obtained with a standard DIC algorithm for the smallest angle value presented Fig. 5(a). One defines the disturbed area as the circular region where the perturbation (the standard deviation of the displacement field) is equal to 150% of the perturbation on an unaffected region. In the case presented Fig. 5(a), the disturbed area radius equals 165 pixels. Considering now the corresponding blurred image Fig. 16(b), the average grey level value on the perimeter of the disturbed area is 69% of the total dynamic range. To have a conservative threshold, 66% (2/3) of the dynamic range is chosen and applied thereafter to every weighting function calculation. Fig. 16(c) presents the p j function obtained for the image Fig. 16(a) with this criterion. ϕ i functions are optional. They may be useful if the shooting with the different lighting conditions are not done exactly with the same point of view, for example because the rapid mirror or shutter motions on a digital single lens reflex induce camera displacements. Without the ϕ i functions, a band of high displacement uncertainty appears in the present case at the frontier defined by the weighting function p j . The affected elements are the ones crossed by the step of p j function. The ϕ i functions correspond to displacement fields defined on all the image, such as rigid body motions.

For this reason, and because the disturbed areas are excluded by the p j functions during the minimization of Eq.3, the influence of the long wave variation of the grey level due to the lightning change is negligible in comparison with the short wave grey level heterogeneity due to the texture, so that the β i coefficients are well evaluated.

For what we have tested so far, ϕ i functions for uniform camera plane motions were enough to circumvent this problem.

Validation

The same optical setup as in section 2 is used for validation. The observed object is the same 82.5mm-sphere with the same speckle. For each angle, two light sources (b 1 (as in section 2.2) and b 2 in Fig. 2) are successively illuminating the sphere, with one image taken for each light source. areas where the masks are used, the perturbation is similar to the one obtained with a standard algorithm using a single image for each angle and an area not affected by a reflection (i.e. A0 on Fig. 7). Area A4, where the reflection is centered in the reference image has also exactly the same trend. One thus considers the algorithm as satisfactory for the goals presented previously.

Discussion

Apart from the main result which is the reflection-proofness of the algorithm, three other facts are worth noting.

First, the perturbation slowly increases for every area. The same phenomenon was observed with the DIC and GLC-DIC algorithms for the A0 area, but to a lesser extent in the latter case: GLC-DIC algorithm is able to partially compensate reflection problems for large angles. This proves that for such important rotation amplitudes, the texture evolution does not strictly satisfy the optical flow conservation assumption.

A grey level correction improves the results, but even for GLC-DIC algorithm, the perturbation slightly rises in areas unaffected by the reflection. Based on the speckle study presented in section 2.4 -emphasizing important normal variations -a reasonable hypothesis is that the texture evolution can not be accurately modeled by a strain and a variation of illumination at the element scale.

Second, the use of several pairs of images to assess a displacement field instead of one will necessarily lead to a lower uncertainty. The camera noise n(x) was assessed to be a white gaussian one. Since the displacement field is obtained through the summation of the contribution of k pairs of images (here k = 2), the central limit theorem implies that the measurement uncertainty (the standard deviation of the displacement field) is divided by √ k (here √ 2) in comparison with a standard algorithm using only 1 pair of images. If the perturbation of A0 area is compared on Fig. 7 and Fig. 17, one remarks a ratio slightly higher than √ 2. This small difference may be explained by the size of the areas, chosen small enough to have several of them along the reflection trajectory, but not large enough to be statistically representative.

Third, the choice of the weighting function is here arbitrary. A binary function, excluding a region around the reflection, is a first proposition and may be improved.

Theses functions were not specifically chosen to be optimal in term of displacement noise reduction. However one notes that even such a simple function leads to a tremendous improvement of the perturbation without any noticeable disturbance at the sharp frontier of the exclusion region.

Conclusion

The goal of this work was first to assess the displacement uncertainty due to the presence of a reflection when a DIC algorithm is used at the macroscopic (millimeter, centimeter) scale. The proposed protocol has revealed that the local displacement uncertainty rises up to 5 times the usual one when a reflection occurs. In the case of a DIC algorithm with grey level correction, this increase is less important for large rotation angle values (> 2 • ) but nearly unchanged for small ones. The point source lightnings that were used are of course more critical than more homogeneous ones, but the observed object has a rather soft curvature, so that the obtained results are considered representative of usual experiments.

An observation of the topography of the speckle (a black and white paint spray)

has shown an important sub-droplet (and thus sub-pixel) roughness that explains the grey level fluctuations at the pixel scale, spoiling the calculation the gradient of the texture.

To circumvent this reflection issue, a new DIC has been proposed. It performs a single minimization with several images where the reflections are located in different regions. For each image, weighting functions are used to 'exclude' the reflection regions from the calculation, while the necessary information is obtained from the other images.

Even with non-optimized weighting functions, very good results are obtained: the displacement is no longer disturbed by the reflection, especially for small angle values (< 2 • ) where the DIC algorithm with grey level correction was unable to improve the results.

Eventually, this new technique with multiple images opens the way to new experimental protocols: more practical ones, like using a single, lightweight, mobile, lightning equipment to illuminate large objects that are submitted to slow evolution (building, ...), or more efficient ones, like using several constant lightnings with different colors and an RGB sensor to capture fast phenomena without reflection issues.

Fig. 1

 1 Fig. 1 Artificial example of evolution of the reflection and shadow positions in the case of a rotation and a translation of a cylinder. First, the apparent displacement near the rotation axis is nearly zero, while the reflection get closer to it. Second, the reflection apparent displacement is smaller than the real translation of the object.

Fig. 2 Fig. 3

 23 Fig. 2 Optical setup: (a) object, (b 1 ) light source, (b 2 ) additional light source used in section 3.3, (c) DSLR, (d 1 ) rotation stage for objet tilt, (d 2 ) used rotation stage for light tilt, (e) translation stage.

Fig. 5 (

 5 a) shows the measured displacement Uy along Y for a rotation equal to 0 • 20 . This displacement along Y (parallel to the rotation axis) is due to the perspective effect. It is presented here instead of the displacement along X because it has a smaller magnitude, so that the perturbations are comparatively more visible. Even though the general trend is satisfactory, a slight noisy area is noticeable. It corresponds to the region A4 of the 0 • image depicted Fig. 4 (a) where the reflection occurs. The displacement along Y corresponding to a 14 • rotation is presented Fig. 5 (b). The magnitude of

Fig. 4

 4 Fig. 4 Images of the sphere with the reflections for an angle equal to 0 • (a) and 14 • (b). The areas used for the perturbation assessment are represented by the yellow rectangles. The value of the grey level (in percent of the total dynamic range) along the white line crossing the disturbed area f line (x) is also plotted to better show the displacement of the reflection.

Fig. 5

 5 Fig. 5 Vertical displacement Uy measured for a 0 • 20 (a) and 14 • (b) rotation.

Fig. 6 Fig. 7

 67 Fig. 6 Mean displacement for the different areas, as a function of the rotation angle of the object αc in the case of a standard DIC algorithm.

Fig. 8

 8 Fig. 8 Perturbation P for the different areas, as a function of the rotation angle of the object αc in the case of a DIC algorithm with grey level correction.

Fig. 9

 9 Fig.9Evolution of the grey level for several pixels located along the rotation axis O as a function of the rotation angle α i .

Fig. 10

 10 Fig. 10 Confocal chromatic measure of the speckle surface: (b) roughness (a 4-order bipolynomial fit have been subtracted to the raw measure) and (a) intensity map.

Fig. 11

 11 Fig. 11 Surface analysis: local slope (a) and Intensity map (b) on the same 1 mm 2 area. Black lines stand for the pixel size.

Fig. 12

 12 Fig. 12 Principle of multiple illuminations Digital Images Correlation.

Fig. 13

 13 Fig. 13 Perturbation P for the different areas, as a function of the rotation angle, with the proposed algorithm and no weighting.

Fig. 14

 14 Fig. 14 Example of p j weighting function calculation: (a) image, (b) blurred image, (c) p j obtained with subsequent values inversion and binarization (threshold = 66% of the dynamic range).

  The obtained displacement fields Uy are presented Fig.15for the first and last angle increments (0 • 20 and 14 • ). No disturbed areas are visible. The p j weighting functions are presented Fig.16: P 1 (resp. P 2 ) corresponds to the images taken with b 1 (resp. b 2 ) lightning source in the reference position and with αc = 0 • 20 or 14 • 00 . The excluded areas of each pair are not overlapping, which would result in a displacement field defined by any optical informations on the intersection area. For the large angle value, one clearly sees the two exclusion areas of the reference image and the deformed image for each p j , since the disturbed areas have 'moved' about 200 pix. between these two images. The perturbation results are presented on Fig.17

  . The perturbation of the different areas is now following the same trend. For area A0 which is never affected by the reflection, the perturbation is lower than with the standard algorithm since two pairs of images are used to assess the displacement field for each angle. For other

Fig. 15

 15 Fig. 15 Vertical displacement Uy measured for a 0 • 20 (a) and 14 • (b) rotation with the proposed algorithm.

Fig. 16 PFig. 17

 1617 Fig. 16 P j functions for the first and the last rotation angle values generated and used by the algorithm. The red rectangles correspond to the areas A0 to A4.

Acknowledgments

The authors thank the research team Eikology of the LMT-Cachan laboratory for the helpful discussion on the subject, Y. Quinsat (LURPA, ENS Cachan / CNRS EA 1385) for his help on the confocal chromatic measurement and C. Galetta and S. Mottola (Master 1 students) for their preliminary work with the authors on the reflection influence on DIC.


Y B Kim

Y T S Noah

Choit

PERIODIC RESPONSE OF MULTI-DISK ROTORS WITH BEARING CLEARANCES

The forced steady state response of a multi-disk rotor system involving a clearance at one of the bearings is determined by using a harmonic balance approach. The impedance method is applied to each of the harmonic steady state components in order to reduce the system to its displacement at the non-linear bearing support. The results reveal the interrelated roles of the bearing clearance, mass eccentricity and side force in producing dangerous subharmonics. The significant effects of the strong non-linearity of a bearing clearance are studied as related to the various system parameters. The results show that the approach developed in this study is computationally superior to numerical integration methods in analyzing multi-disk rotor systems with strong non-linearity.

INTRODUCTION

Several experimental studies [1][2][START_REF] Ehrich | ASME Paper No. 66-MD-1 Design Engineering Conference and Show[END_REF] have demonstrated the occurrence of undesirable subsynchronous whirling of rotors interacting through a clearance with their housing or supports. The present work extends the incremental harmonic balance method developed earlier [START_REF] Choi | [END_REF] for a modified Jeffcott rotor model with bearing clearances to the determination of the periodic response of non-linear multi-disk rotor systems. The significant effects of the strong non-linearity of a bearing clearance are studied as related to the various system parameters.

Several methods have recently been advanced for determining the periodic solution of single-disk rotor systems with non-linear supports [START_REF] Choi | [END_REF][START_REF] Yamamoto | Memoirs of the Faculty of Engineering[END_REF][START_REF] Childs | [END_REF][7]. For application to larger, multi-disk rotor systems, Nataraj and Nelson [START_REF] Nataraj | American Society of Mechanical Engineers Design Conference[END_REF] developed a periodic solution method based on a collocation approach for the response of the rotor. They utilized a subsystem approach [START_REF] Cipra | [END_REF] to reduce the size of the resulting system of algebraic equations.

In general, there exist three different response types for a rigid rotor with a bearing clearance: (i) periodic response of the same period as that of the rotational speed of the rotor, (ii) subor supersynchronous response, and (iii) aperiodic response. To determine the response of the first or the second type, a harmonic balance method (HBM) with OFT (Discrete Fourier Transform)/IDFT (Inverse Discrete Fourier Transform) can be utilized for a multi-disk rotor with a bearing clearance, as presented in this paper. However, the third response type, which is usually more complicated than the first two [START_REF] Neilson | response of a rigid rotor mounted on discontinuously nonlinear supports[END_REF], cannot be obtained by using this approach. A specially modified DFT/IDFT algorithm is needed to obtain this response since the frequency components are incommensurate. Such a modified algorithm is being developed separately and is not included in the present paper.

In the present work, it is shown that the harmonic balance method, as applied to a large, non-linear rotor support system, can be made to be highly efficient. This is achieved by using a version of an impedance, or dynamic stiffness [START_REF] Noah | [END_REF][START_REF] Fan | [END_REF], in which the steady state imbalance response of a spinning rotor is described by displacements at the bearings. This is done by using a partition of the impedance matrix to relate these displacements to the imbalance forces. A similar approach was independently developed in reference [START_REF] Jean | American Society of Mechanical Engineers, Vibrations Conference[END_REF] to reduce a set of n algebraic equations, resulting from a collocation method, applied to a rotor system, to only those co-ordinates which are directly coupled to the system's non-linear components. The impedance method is applied to each of the harmonic steady state relations corresponding to each of the harmonic components of the assumed periodic solution. For illustration, the method is applied to a two-disk cantilevered rotor system with one of its bearings possessing a deadband clearance.

ANALYSIS

EQUATIONS OF MOTION

To demonstrate the application and efficiency of the harmonic balance method as applied to multi-degree-of-freedom finite element rotor models, a typical rotor with non-linear bearing is analyzed to determine its non-linear response to imbalance. The rotor model considered is shown in Figure 1. A bearing with deadband clearance between the outer race and the housing is located at node 2. The equation of motion for the rotor-support system in terms of finite elements can be written as

[M] {R} + [ C] {R} + [K] {R} = {f} + { g( t)}, (1)
where [M] is a mass matrix, including diametral moment of inertia for the disks, [ C] includes damping terms and gyroscopic terms, and [K] is a stiffness matrix (a list of nomenclature is given in the appendix). The displacement vector {R} is defined as { R} = { Rtx, R2x, R3x, Rly, R2y, R3y, {34y, /3sy, {36y, f34x, f3sx, f36x} T' where R; and {3; denote displacement and rotation, respectively. A vector R has the dimension of 12 x 1. The nodal displacements 7x, 7y, 8x and 8y of Figure 1

I

When contact/non-contact occurs between the bearing outer race and the bearing support at node 2, the non-linear force terms can be expressed as {f} = {O,fx, 0, OJn 0, 0, 0, 0, 0, 0, O}T, in which [START_REF] Choi | [END_REF], with F, as a constant side force, (

) =KR 1- +F J, B 2x J 2 2 .< R2x+ R2,. 0
where 6 denotes the bearing radial clearance and K 8 stands for the bearing support stiffness, andf, u;.) represents the X (y) direction non-linear force at node 2. The sinusoidal excitation force vector which occurs at nodes 1 and 3, { g( t)}, can be described by {g(t)} = {gs} sin wt+{gc} cos wt.

(

) 3
In order to facilitate the formulation, the co-ordinates are rearranged so that the coordinates at the gap are listed first, {f} = Ux..f.., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} 1 , in which the matrices [M], [C] and [K] and the force vector {g(t)} are also changed accordingly.

STEADY STATE SOLUTION

The steady state, periodic solution can be assumed as follows. For the ith node, J;.=Co,.+ L C""cos-+D""sin-.

n ~I 11 11 (4)
(5)

A selection for the integer m to assign to 11 allows the determination of the 1/ m order subharmonic responses. Substituting in equation (I) from equations (3), ( 4) and ( 5) and applying the harmonic balance method result in

[K]{A} = { C}, (6)
where [K] is a stiffness matrix and {A}= {A~x, A~y, A6x, ... , A~,.} T' {C}={Co., Co,,O, ... ,O}T.

For matrix reduction, the vectors in equation ( 6) can be written in the following partitioned form: for the constant term, where

K ({A%.}) _ ({ C~}) [ ] {A~} -{ C~} ' {A%.}= {A~x. A~yV, {A~}= {A6x
, A6v, A~x. A~"' ... , Agx, Agv}T, { C%.} = {Cox, Co,}T, { C~} = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

Here the subscripts m and s denote master and slave, respectively. The [K] matrix can be partitioned as

( [KmmJ(2x21 [KmsJ(2xiO)) [K] = [K,m](IOx2l [K,,]ooxiOl •
The aim of this partition is to express the motion in terms of a reduced model in the master co-ordinates at the gap. From equation ( 7), [KmmHA%.}+ [KmsHA~} = { C%.}, Combining equations ( 7) and ( 8) gives It is noted that the constant coefficients Cox, Cov in equation ( 5) are now expressed as functions of the unknowns A~x, A~n ... , B;,x, B~Y. Equation ( 10) is an implicitly defined non-linear algebraic matrix to which the Newton-Raphson second order iterative method can be applied. To apply the Newton method, equation ( 10) is cast in an incremental form,

[[Kmm]-[Kms][K.sr 1 [K,m]]{A%.} = { C%.}.
Similarly, the coefficients associated with the cosine terms lead to

-(n;r [M]{A}+C;) [C]{B}+[K]{A}={Cn}+{~~cl, :~:::}. ( 13
)
where

(n = 1, 2, ... , N),
For the sine terms, Combining equations ( 13) and ( 14) leads to

-(nw) 2 [M]{B}-(nw) [ C] {A}+ [K] {B} = {Dn}+ {{gs} ~ff n = 11 }, ( 14
)
[ S] { q} = { w} + { u}, where {q} ={A~" A~v• B~x, B~" A~" A:,,., B~" B~,, ... , A~., A~"' B~" B~,r.

{u} = {g;, 0, 0, g;, g:., 0, ... , O}r.

To reduce the size of equation ( 15), it is partitioned as

( [SmmJ(4x4) [SmsJ(4x20)) ({Qm}) =({Wm})+({Um}) [SsmJ(20x4) [S,](20x20) { Q,} 0 { U,} '
where

{ Qm} == {A~x, A~v, B~x, B~y}r, From equation (16), { { 2 2 2 2 }T Wm} = C nx• C nv• Dnx• Dny • (15) fl6) [Smm]{ Qm}+ [Sms]{ Os}= { Wm}+{ Urn},
Combining equations (17a) and (17b) gives [SN) where n denotes the nth harmonic mode. A similar expression can be written for .JC"Y and .JDnv for they direction. From equation ( 21 Here the superscripts denote the degree of freedom number and a subscript denotes the harmonic term in x or y direction. The total dimension of the unknown vector {.Jr} is 4N + 2. This [P] matrix has to be found by using numerical differentiation with discrete and inverse discrete FFT, and it has to be updated at each iteration until the Euclidean norm of the unknown vector {.Jr} converges. More details on the calculation procedure are given in the paper by Choi and Noah [START_REF] Choi | [END_REF].

[Sm,HOm}+[S,]{Q.} = { U,}. (17a, b) [[SmmJ-[SmsUSssr 1 [S,m]]{ Om} = { Wm} + {Urn}-[SmsJ[Ss,r 1 { U,}. (18) With [S] = [SmmJ-[Sms][Sssl-1 [S,m] (4 X 4), equation (

NUMERICAL RESULTS AND DISCUSSION

The methods as formulated in the previous section have been applied to the rotor shown in Figure 1. A density of 7835 kg/ m 3 and an elastic modulus of 2•07 x 10 11 N/ m 2 were used for the shaft. The geometric data of its three elements is listed in Table 1.

The first three critical speeds of the linear rotor in absence of gyroscopic terms were computed by using the harmonic balance method. All the critical speeds are compared with results for the eigenvalues as computed by the commercial package (IMSL) subroutine EIGZF. The results are identical and a Bode plot is shown in Figure 2. A comparison of the accuracy of the direct numerical integration (fourth order Runge-Kutta) and the harmonic balance method (HBM) is shown in Figure 3. The steady state response obtained by using numerical integration was calculated by using equation (1 ). For the --, HBM; e numerical integration. application of the HBM method in a numerical procedure, the first four harmonic coefficients were used, which is the most efficient in calculating subharmonics [START_REF] Choi | [END_REF]. As shown in Figure 3, a dominant 1/2 order subharmonic response obtained by using the HBM method demonstrates good agreement with the direct numerical integration results, not only for the waveforms for displacement but also for phase. The slight difference (within 5% error) is due to ignoring the higher harmonic terms in the HBM method. As more harmonics are added, the HBM method results converge to those obtained by using numerical integration. The calculation time of the HBM method is superior to that of the numerical integration method. For most simulation cases, the numerical integration method takes several hundred times longer than the HBM method to yield the steady state response for a typical multi-degree-of-freedom rotor model such as the one utilized herein. However, better computational performance of the numerical integration can be achieved if (i) the integration step size is taken larger, (ii) the system damping is large, and (iii) the system has a small number of degrees of freedom.

When the rotor is subjected to imbalance forces, damaging 1/2 order or the lower subharmonics near twice or n times the natural frequency appear in the presence of the deadband clearance and a side (radial) force due to gravity and misalignment. To demonstrate the occurrence of significant subharmonic response, the damping was set to a very small value. The results presented in Figure 4 and subsequent figures were obtained with account taken of the gyroscopic moments due to rotation of the two disks. As shown in Figure 4, the 1/2 order subharmonic response occurs near twice the second natural frequency and is dominant. This can cause a damaging effect over the operating speed. The shaft center may trace either circular or elliptical paths under the harmonic forces. However, in cases of sub harmonic response, the paths will be of more complicated forms. In Figure 5, a typical 1/2 order subharmonic response with harmonic shaft trajectory is shown. The existence of the 1/3 order and the lower order subharmonic components .:: I 020 rad/ s; radial were examined for various values of the side force, gap size and eccentricity, near the natural frequencies. As well as the 1/2 order subharmonic response, the 1/3 order subharmonic response can be as dangerous since the peak to peak precession amplitude of a 1/3 order subharmonic can be very large in some extreme cases, as shown in Figure 6. The figure shows that the dominant response is the 1/3 order subharmonic component. A 1/4 order subharmonic response can also exist under certain conditions, as shown in Figure 7. These results confirm that high-speed rotating machinery with asymmetry in

0•5 X 0•000006
"' I E Q "' -c Q) E Q) u 3•2 2•2 1•2 % -0•8 (/) -o "' -2•8 \ ' ', -I 0 I 2
x displacement x 10-5 (m) the bearing supports might have a lower order subharmonic vibration response [START_REF] Ehrich | [END_REF]. In a recent paper submitted by the present authors [START_REF] Kim | Nonlinear Dynamics[END_REF], the same whirling shapes for a model presented in reference [START_REF] Ehrich | [END_REF] were obtained by using the harmonic balance method. The effects of the side force, eccentricity and gap size play a significant role in the occurrence of dangerous subharmonics as mentioned above. In earlier work, Choi and Noah [START_REF] Choi | [END_REF] showed that dangerous sub harmonic responses could develop for a rotor system for a certain range of the side force magnitude. The effects of the side force and mass eccentricity were investigated, and the results are presented in Figures 8, 9 and 10. In Figure 8, with gap sizes smaller than 2 •54 x 10 -s m, the eccentricity has the dominant effect on the subharmonic response. The figure shows that for a certain range of side force magnitudes (between 50 N and 250 N), the subharmonic response can exist. In Figure 9 are shown the side forces and eccentricities at which the maximum 1/2 order subharmonic component occurs. The figure shows that for a given bearing clearance higher side force is required to produce significant subharmonic motion when the eccentricity is small (lower than 5 x 10-5 m). However, for side forces higher than some value, the eccentricity has no effect on the occurrence of subharmonics, since the bearing outer race will be already in contact with the housing. This trend is the same for all gap sizes. The effects of the side force and gap size were studied and the results are shown in Figure 10. This figure shows values of the side forces and gap sizes at which the maximum 1/2 order subharmonic response occurs. The figure also shows that the gap size is significant in determining whether or not large subharmonics occur. This figure shows that if the gap size becomes larger, more side force is required for the occurrence of the subharmonics. This is because for larger gap sizes more side force is needed to keep intermittent contacting between bearing and housing. Therefore, a well balanced rotor will not have the damaging subharrnonic response (or subsynchronous motions) even if there exists a deadband between bearing outer race and housing. If the rotor is poorly balanced, a possible procedure to avoid a subharmonic response is to decrease the side force or the gap size. The effect of damping at the bearing was investigated and the result is shown in Figure 11. The figure shows that the increase in damping results in a smaller subharmonic response, with a more significant effect on reducing or eliminating the subharmonics being observed for values of the viscous damping coefficient larger than 175 N s/ m. The side force does not affect this tendency.

CONCLUSIONS

An effective numerical algorithm for determining the steady state forced vibration of a strongly non-linear multi-degree-of-freedom rotor system has been developed. The method provides an accurate (with ~5% error) and efficient means for determining the non-linear effects on rotor dynamic behavior. In addition, the complicated non-linear steady state periodic motions of multi-degree-of-freedom rotor systems can readily be studied by using the present method. Even though the method can accommodate any number of retained harmonics, there is a trade-off between level of accuracy and computational cost. In general, four harmonics were found to be adequate for good accuracy at a reasonable cost. Another advantage of the harmonic balance method is that it is capable of providing the unstable as well as the stable solutions. This is necessary if the global dynamic behavior of the system is to be determined. This includes bifurcation analysis and possible occurrence of chaos.

The newly developed algorithm was applied to a typical multi-disk rotor/bearing system. Dangerous subharmonic resonances may occur in the presence of bearing clearances, mass eccentricities and side forces due to gravity or misalignments. The results show that the ranges of the above three parameters are mutually related in determining the sub harmonic occurrence. The side force and size of radial clearance within certain ranges can induce significant subharmonic responses. These results can be explained by the contacting motion between bearing outer race and housing which might render the system highly non-linear. Once the stiffness starts changing due to the contact motion, the stiffness variation is not affected by increasing the level of the side forces. These characteristics of stiffness variation due to the bearing-clearance effects well agree with other previous experimental results [2] and analytical results [START_REF] Childs | [END_REF][START_REF] Jean | American Society of Mechanical Engineers, Vibrations Conference[END_REF]. APPENDIX: NOMENCLATURE master degree of cosine coefficient of steady state solution (2 x 1) incremental vector of {A?.,} nth cosine coefficient of the steady state solution in x, y direction at ith node slave degree of cosine coefficient of steady state solution (2 x I) nth sine coefficient of the steady state solution in x, y direction at ith node damping matrix (12 x 12) constant cosine coefficient vector of non-linear force (12 x 1) nth cosine coefficient of non-linear force in x, y direction master degree of non-linear force vector with cosine coefficient (2 x 1) slave degree of non-linear force vector with cosine coefficient ( 10 x 1) constant sine coefficient vector of non-linear force ( 12 x 1) nth sine coefficient of non-linear force in x,y direction stiffness matrix (12 x 12) bearing stiffness reduced impedance matrix with constant terms (2 x 2) mass matrix (12 x 12) master degree of trigonometric coefficient vector of steady state solution ( 4 x 1) master degree of trigonometric coefficient vector of steady state solution with nth harmonic mode ( 4 x 1) incremental vector of { Q?.,,} slave degree of trigonometric coefficient vector (20 x 1) displacement vector (12 x 1) x,y rectlinear displacement at ith node total trigonometric coefficient matrix (24 x 24) reduced impedance matrix with trigonometric terms ( 4 x 4) reduced impedance matrix with trigonometric terms of nth harmonic mode ( 4 x 4 l overall coefficient matrix (4N + 2, 4N + 2) equivalent trigonometric force vector ( 4 x 1) master degree of trigonometric coefficient vector of non-linear force (20 x 1) incremental vector of {worn,} master degree of trigonometric vector of non-linear force with nth harmonic mode ( 4 x I) Non-linear restoring vector (12 x 1) non-linear restoring force vector in x,y direction imbalance force vector (12 xI) imbalance force vector of cosine, sine terms ( 12 x 1) constant trigonometric coefficient vector of steady state solution (24 x 1) overall incremental coefficient vector ( 4N + 2, 1) total imbalance force vector (24 x 1) overall incremental coefficient vector of non-linear force (24 x 1) constant trigonometric coefficient vector of non-linear force (24 x 1) overall force vector (4N + 2, I) rotational displacement in x,y direction at ith node radial clearance

Figure

  Figure I. Finite element rotor model.

  A~x = L A~x cos-+ B~x sin-' ,.+ I A~vcos-+B~ .. sin-,If the system reaches the steady state, the non-linear force can be expressed as[START_REF] Choi | [END_REF]

  where [K] is a reduced impedance matrix with constant terms since { C%.} is a function of the displacement of node 2 at the bearing clearance as shown in equation (7), or Cox= Cox(A~x. A~"' Ai" Ai, ... , A;_,x, A~Y' B;_,x, B;_,y), (11a)Coy= C0y(A~., A~Y' Ai., Aiy, ... , A~., A~Y' B~., B~y).(llb)

  18) can be expressed as [S]{ Om} = { Wm} +{ 0}, 09) where [S] is a reduced impedance matrix with trigonometric terms. In an incremental form, equation (19) can be written as (20) where the subscript n indicates the nth harmonic term, n = 1, 2, ... , N. Now, upon combining equations (12) and (20), the overall incremental equation of the system takes the form where and [T]{.1r} -{.1v} = {z}, {.1r} = {.1A~x' .1A~Y' .1Afx, .1Afv, .1Bix, .1Bi" .1A~" .1A~.., .1B~x' .1B~ .. , ... , .1A7-,x, .1A7-,v, .1B7.,., .1B7-,,r, (21) { .1v} = {.1Cox, .1Cov, .1C~." .1Ctv, .1D." .1Dil , .:lC2x, .:lC2_.., .:lD2., .:lD2_,., ... , .1CNx, .1CNy• .1CNx• .1CNy}T ( [K] [T]~ : [S.] ~ ) (4N +2, 4N +2).

  ), {.Jv} is written with partial derivatives of {.Jr} terms since {.Jv} is a function of {.Jr}, as shown by equation (22). The partial derivative terms yield the Jacobian matrix, so one can finally obtain the equation where ( aCox/aA~x [P]= aCo<aA~x aDNvfaA~x ([ T] + [P]){.Jr} = {z}, aCoxfaA~v aCov/aA~v aCo,/aAix aCovfaAix a Cox! aB~,.) aC0 .. /aB~}'

Figure 2 .

 2 Figure 2. Bode plot of the radial response at the bearing with clearance (no clearance, no gyroscopic effect).

• 0 "Figure 3 .

 03 Figure 3. Comparison of non-linear rotor response at the bearing with numerical integration at rotational speed = 1020 rad/s; radial clearance = 1•27 x 10-5 m, mass eccentricity = 2•54 x 10-5 m, side force = 10 N.

5 __

 5 ~~~o~i=-____________ -'---"""-~-~-~-=-=-=--=-~-0•000004 IX 0•000002 0 Whirling frequency (rad/s) x 100

Figure 4 .Figure 5 .

 45 Figure 4. Comparison between subharmonic and harmonic responses at the bearing; radial clearance 1•27 x 10-5 m, mass eccentricity = 1•27 x 10-5 m, side force = 44•5 N.

Figure 6 .

 6 Figure 6. Orbit of the subharmonic response at the bearing with rotational speed clearance = 1•27 x w-s m, mass eccentricity = 5•08 x w-s m, side force = 44.5 N.

Figure 7 .

 7 Figure 7. Orbit of the subharmonic response at the bearing with rotational speed clearance = 1•27 x w-s m, mass eccentricity = 0•89 x 10-5 m, side force = 22•25 N.

Figure 8 .

 8 Figure 8. Effect of side force and eccentricity at the bearing with frequency = 1180 rad/s; radial clearance= 2•54x 10~5 m. Eccentricity x 10~5 (m): -e-. 2•54; -+-, 5•08; -0-, 7•82; -o-, 10•18; --x-. 12•70.

Figure 9 .Figure 10 .

 910 Figure 9. Side force and mass eccentricity for peak subharmonic component at 1180 rad/s.

Figure 11 .

 11 Figure 11. Effect of damping at bearing on the sub harmonic response at 1180 rad/ s; radial clearance = 2•54x 10-5 m and mass eccentricity= 12•7 x 10-5 m. Side force (N): -----• 10; -+-, 150; -0-, 200.

TABLE 1

 1 Rotor configuration and physical parameters

	Item
	Young's modulus
	Mass of first disk
	Mass of second disk
	Shaft outer diameter
	Shaft inner diameter
	1st element length
	2nd element length
	3rd element length
	Non-linear spring stiffness
	1st disk diametrical moment of inertia
	2nd disk diametrical moment of inertia

ACKNOWLEDGMENTS

This work was carried out as part of a research project supported by NASA, Marshall Flight Center under contract No. NASS-36293. The authors are grateful to Thomas Fox, the technical monitor, for his enthusiastic support and interest.