Myriam Gharbi

Patrice Peyre
email: patrice.peyre@ensam.eu

Cyril Gorny

Muriel Carin

Simon Morville

Philippe Le Masson

Denis Carron

Rémy Fabbro

Influence of various process conditions on surface finishes induced by the direct metal deposition laser technique on a Ti-6Al-4V alloy

Keywords:

The direct metal deposition (DMD) with laser is a free-form metal deposition process for manufacturing dense pieces, which allows generating a prototype or small series of near net-shape structures. One of the most critical issues is that produced pieces have a deleterious surface finish which systematically requires post machining steps. This problem has never been fully addressed before.

The present work describes investigations on the DMD process, using an Yb-YAG disk laser, and a widely used titanium alloy (Ti-6Al-4V) to understand the influence of the main process parameters on the surface finish quality. The focus of our work was: (1) to understand the physical mechanisms responsible for deleterious surface finishes, (2) to propose different experimental solutions for improving surface finish.

In order to understand the physical mechanisms responsible for deleterious surface finishes, we have carried out: (1) a precise characterization of the laser beam and the powder stream; (2) a large number of multi-layered walls using different process parameters (P(W), V(m/min), D m (g/min), Gaussian or uniform beam distribution); (3) a real time fast camera analysis of melt pool dynamics and melt-pool -powder stream coupling; (4) a characterization of wall morphologies versus process parameters using 2D and 3D profilometry.

The results confirm that surface degradation depends on two distinct aspects: the sticking of nonmelted or partially melted particles on the free surfaces, and the formation of menisci with more or less pronounced curvature radii. Among other aspects, a reduction of layer thickness and an increase of melt-pool volumes to favor re-melting processes are shown to have a beneficial effect on roughness parameters. Last, a simple analytical model was proposed to correlate melt-pool geometries to resulting surface finishes.

Introduction

The direct metal deposition (DMD) laser technique is a recentlydeveloped manufacturing technique that allows obtaining drafts of complex metallic parts from a three dimensional CAD model as reminded by [START_REF] Pinkerton | Laser direct metal deposition: theory and applications in manufacturing and maintenance[END_REF], even if technical issues still exist (stability of the DMD, process control, geometrical limits) before a widespread industrialization. The process uses a high power laser (usually Nd:YAG, fiber or diode laser) focused onto a metallic substrate to generate a molten pool, where a coaxial powder feeding increases the material volume and contributes to the formation of a solid layer. The substrate is then scanned on (x, y, z) directions relative to the laser + powder nozzle head and the layer-by-layer additive manufacturing of complex 3D shapes becomes possible. Many materials, including graded materials, have already been more or less successfully attempted, with a specific focus on aeronautical materials. For instance, considering the specific case of Ti-6Al-4V titanium alloy, [START_REF] Maisonneuve | Profil project: direct manufacturing of aerospace components by laser cladding and laser sintering[END_REF] or [START_REF] Bontha | The effect of process variables on microstructures in laser-deposited materials[END_REF] pointed out the dependence between process parameters and grain structure (columnar or equiaxed), whereas [START_REF] Brandl | Mechanical properties of additive manufacturing titanium (Ti-6Al-4V) blocks deposited by a solid-state laser and wire[END_REF] demonstrated that the use of a wire instead of powder could lead to good-yet anisotropic mechanical properties. From the previously mentioned publications and many other ones, it appears that, with a DMD laser technique, and whatever the projected material, similar or superior mechanical properties than foundry techniques can be achieved, either in static condition or under fatigue loading.

Among the various factors involved in the DMD process, one can distinguish between 1st order factors: Laser power P 0 (W), scan speed V 0 (m/min), laser diameter d 0 (mm), and mass feed rate D 0 m (g/min), and 2nd order factors: defocussing distance, gas shielding nature and velocity, time pause between passes, spatial distributions of laser and powder etc. . . . All these factors influence the thermal history T = f(x, y, z, t) of the part, and contribute not only to the melt-pool shapes, and the resulting layer growth, but also to the final metallurgical and mechanical properties Recently (2006Recently ( -2011)), intensive numerical work has been carried out to model the DMD process, starting from the laser-powder interaction, to the thermo-mechanical calculation of residual stresses, including or-not metallurgical aspects. Authors like [START_REF] Fathi | Prediction of melt pool depth and dilution in laser powder deposition[END_REF] proposed rather simplified predictive models to predict the geometrical characteristics of the walls by assuming a solid state during the process. Whereas, more complex thermohydraulic calculations considered free moving surfaces, either on 2D multilayers configurations like [START_REF] Morville | 2D finite element modeling of heat transfer and fluid flow during multi layered DMD laser process[END_REF], or on 3D single layer as shown by [START_REF] Qi | Numerical simulation of heat transfer and fluid flow in coaxial laser cladding process for direct metal deposition[END_REF] or [START_REF] Kumar | Effect of three-dimensional melt pool convection on process characteristics during laser cladding[END_REF]. Such numerical approaches allowed calculating wall dimensions, fluid flow and realistic temperature distribution inside the clad layers.

In turn, the DMD process is now globally well understood even if limitations still remain such as:

1) the need of a robust and reliable process control to maintain constant local thermal conditions, and a stable and constant layer growth; 2) the occurrence of a deleterious surface finish (usually R a > 15 m as shown by [START_REF] Maisonneuve | Profil project: direct manufacturing of aerospace components by laser cladding and laser sintering[END_REF] on Ti-6Al-4V or [START_REF] Pinkerton | The effect of laser pulse width on multiple-layer 316L steel clad microstructure and surface finish[END_REF] on AISI 316L), that requires intensive post-machining steps, and restraints the applications of DMD to the drafts of final parts.

A common thought is that optimum surface finishes are mostly expected when using thin and stable DMD layers, without more precisions concerning process parameters, and melt-pool shapes and dynamics. Moreover, as reminded by [START_REF] Pinkerton | Laser direct metal deposition: theory and applications in manufacturing and maintenance[END_REF], one can assume that melt-pool stability depends on various internal forces (Marangoni flow, gravitational forces) or external forces (gas pressure, recoil forces near vaporization point, dynamic forces applied by the projected powder grains) applied on the molten metal. For instance, negative thermo-capillary coefficients d/dT (N/m K) are expected to provoke centrifugal Marangoni fluid flow directed away from the melt-pool center. This flow could act as a shaping contribution to the melt-pool and modify the resulting surface finish, during layer additive processes, even if such an effect has to be confirmed experimentally.

Recent work has considered surface quality as an important factor to address in itself to improve the process. For instance [START_REF] Zhu | The influence of laser and powder defocusing characteristics on the surface quality in laser direct metal deposition[END_REF] has shown the benefit versus surface quality of positioning the powder focus below the melt-pool, and the laser focus above. A melt-pool enlargement was shown to be the main contributor to optimum surface smoothness. It was also shown by [START_REF] Alimardani | Surface finish in laser solid freeform fabrication of an AISI 303L stainless steel thin wall[END_REF] that an increase of the scanning speed and a real-time control of the melt-pool dimension and temperature could significantly improve surface finish.

In this context, it seemed important to understand the specific contributions of various processing parameters to surface finish, considering simultaneously the variation of melt-pool characteristics, with the final objective of reaching optimum surface finishes, and limiting post-machining steps. Considering a widely investigated Ti-6Al-4V titanium alloy, and a large range of process parameters (laser power P 0 (W), scan speed V (m/min), mass feed rate D m (g/min), laser-powder interaction distance d i (mm), beam distribution (Gaussian or top-hat), gas shielding conditions), this paper addresses surface finish variations, and aims at correlating them with melt-pool shapes and dimensions.

Experimental conditions

Direct metal deposition conditions

DMD tests have been carried out using a HL 10,002 continuous wave disk Yb:YAG laser operating at 1.03 m, with 320-700 W cw laser powers. Two configurations were used:

(1) a 1.3 mm diameter defocused Gaussian laser spot provided by a 0.4 mm optical fiber delivery, and respectively 200 mm collimation and focusing lens in the laser head. The laser beam was used with a +5 mm defocusing condition, resulting in a 1.3 mm diameter beam on the substrate; (2) a 1.68 mm diameter top-hat laser spot obtained with a 600 m optical fiber, a 100 mm collimating lens and a 280 mm focusing lens.

For the powder distribution, a helicoïdal powder delivery nozzle was used, where the Titanium powder (45-75 m, TLS Technik) is delivered coaxially with the laser beam, resulting in a d p ≈3 mm powder focus diameter, and with average mass feed rates D m in the 1 g/min to 3 g/min regime. The relative positioning of laser beam versus powder stream is presented in Fig. 1a (d i : interaction distance between powder stream and laser beam). Argon was used as a driving and shielding gas, in order to ensure powder transportation, and to limit oxidation phenomena. All the DMD tests have been carried out with 40 mm-length/6-15 mm-height walls, starting from 2 mm-thick titanium sheets (Fig. 1b), and with scan speeds in the 0.1 m/min to 0.6 m/min range. Experimental conditions are summarized in Table 1. Last, it has to be mentioned that most of the tests were carried out with a local argon shielding (through the nozzle), whereas a few additional tests were carried out on a Lens Optomec-450 industrial set-up, with a controlled O 2 level of 20 ppm in the whole cell.

The experimental determination of mass feed rate distribution versus x-axis was carried out using a specific device where a metallic sheet with a 0.3 mm diameter laser-drilled hole is moved below the powder nozzle, and allows determining the local mass feed rate at any location of the powder stream. The laser beam profile was analyzed by an industrial beam analyzer. The 1.3 mm diameter laser spot, and powder distributions were both shown to be quasi-Gaussian in shape, as shown in Fig. 2.

In turn, coincidently, similar analytical formulations could be used for the two spatial distributions (Eqs. ( 1) and ( 2)).

D(g • s -1 • m -2 ) = 5 • D m r 2 p • exp - 5x 2 r 2 p (1)
with r p = powder stream radius (=2 mm), v p = powder velocity (m/s), d i = interaction distance (mm), D m = average mass feed rate (g/s)

∅(W • m -2 ) = 5 P 0 r 2 laser • exp - 5x 2 r 2 laser ( 2
)
with r laser = laser beam radius (=0.65 mm), P 0 = laser power (W) On the other hand, additional DMD tests were carried out with a "top-hat" uniform distribution (f = P 0 /(r 2 laser )), and a 1.68 mm diameter, in order to estimate the influence of beam distribution on the melt-pool shape, and surface finish.

Considering the laser-powder interaction during the time-offlight (interaction time t 0 = d i /v p ) of a powder grain, an estimated 320,400,500 0.1,0.2,0.4,0.6 1,2,3 calculation of powder temperature increase T could be made using a recent and simplified model by [START_REF] Qi | Numerical simulation of heat transfer and fluid flow in coaxial laser cladding process for direct metal deposition[END_REF] where the absorption of laser light by powder was adjusted by an absorption coefficient "A" (≈0.4 on titanium) (Eq. ( 3)). In this model, interactions between particles and particle/gas thermal losses were neglected. If we only consider the 1.3 mm diameter Gaussian laser spot, the two distinct interaction distances d i used during our tests (1 mm and 3 mm), correspond either to a powder mostly maintained at a solid state, or to a powder that is molten before contacting the melt-pool (Fig. 3).

T = 3A(x, y) 4r g C p . d i v p (3)
with r g = powder grain radius (mm), v p = powder velocity (≈1.5 m/s), d i = interaction distance (mm).

Real-time diagnostics

A simultaneous coaxial and lateral recording of melt-pool sizes and dynamics was carried out, using two synchronized fast cameras with C-Mos sensors (Photron), at frequency rates up to 5000 Hz, using KG3 filters to cut the laser wavelength (1.053 m), and halogen lights to improve the contrast of recorded images. This allowed us to: (1) investigate powder stream (projected powder velocity v p ), (2) investigate powder/melt-pool interaction, (3) analyze meltpool flow (powder velocity in the melt v p , fluid velocity v f when

Table 2

Thermo-physical properties of Ti-6Al-4V (A: absorptivity; S: solid; L: liquid).

(S/L) [kg/m 3 ] k (S/L) [W m -1 K -1 ] Cp (S/L) [J kg -1 K -1 ] [Pa s] L f [J kg -1 ] A Tm [K] [N/m] = d/dT [N m -1 K -1 ]
4400/4200 20/35 600/700 4.10 -3 3.10 5 0.4 1920 1.5 -2.10 -4 distinct, (4) quantify melt-pool dimensions (length L, height H, area S) and try to correlate them with surface finish data.

Surface analysis

DMD surfaces exhibit a periodic structure composed of an accumulation of lateral menisci directly due to the equilibrium shapes of the melt-pool. Manufactured walls were analyzed using: (1) low and medium magnification optical microscopy to visualize the global surface aspect and the grain microstructure (after Kroll's etching), (2) scanning electron microscopy (SEM-FEG Hitachi 4802 II) to provide a local description of DMD surfaces, (3) 2D and 3D roughness analysis (Veeco Dektak 150 stylus profiler) to quantify surface topography. The surface analysis have been carried out using a 4 mm scanning length perpendicular to the main DMD direction with a 2.5 m tip stylus, and a 0.1 mm/s scanning speed.

Experimental DMD tests

Melt-pool dynamics and melt-pool dimensions versus process parameters

When powder particles reach the molten pool, they follow two options: (1) they ricochet on the surface layers partially covered with Ti oxide, below a threshold collision angle th ; (2) they enter more or less partially the melt-pool, and are transported by the fluid flow where, after full melting, they contribute to layer growth.

Powder grains contacting the melt-pool without ricocheting are transported by a Marangoni centrifugal fluid flow at a V f velocity ≈0.45 m/s estimated with fast camera (Fig. 4). The direction of Marangoni flow is in good agreement with a negative thermocapillary coefficient d/dT (N/m K), which is mostly expected in molten titanium (Table 2).

The ratio of particles contacting the melt-pool, but not contributing to layer growth was checked analytically, by calculating the theoretical catch efficiency ratio ϑ, i.e. the integral value of the D = f(x,y) function (Eq. ( 1)) restricted by the limits of the melt-pool i.e. approximately the melt-pool width e (Eq. ( 4)). The comparison of the analytical prediction, and the experimental mass efficiency (=mass really incorporated in the wall/mass projected) indicates (Fig. 5) 15% higher analytically calculated catch efficiency data. This reveals that approximately 85% of the incident particles projected inside the contours of the melt-pool are really melted whereas, the 15% remaining are attributed to particle skipping that tend to ricochet on the melt-pool surface, to particles that are shifted when interacting with other particles, or to particles that do not fully melt and stick on the surface. This skipping phenomenon is promoted by small incidence angles at the edges of the semi-hemispherical melt-pool, and large fluid velocities in the melt-pool.

ϑ = /2 -/2 D 0 • exp(-5x 2 /r 2 p ) • dx ∞ -∞ D 0 • exp(-5x 2 /r 2 p ) • dx (4)
with r p = powder stream radius (mm), e = melt-pool width (mm), (D 0 = 5D m /r 2 p ). If we now consider the variation of wall dimensions versus (P 0 , V, D m ) process parameters, rather classical results are obtained (Fig. 6 -Wall thicknesses increase with P 0 and 1/V, due to an increase of the lineic energy P/V (J/m), but vary to a little extent with D m as the laser beam filtering is limited in the 1-3 g/min mass feeding range -Layer heights h increase with D m and 1/V, but do not vary with P 0 .

Those results can be discussed considering a recentyet-simple analytical formulation of layer height by Peyre et al.

(

) h = D * m • L/ • V (with D * m (g s -1 m -2 )) 2008
= average surface mass feed rate. However, such a formulation has to be modified to explain why h remains constant with increasing P 0 values. For this purpose, one has to take into account the notion of efficient mass feed rate: as mass feed rate distribution is Gaussian in shape (Fig. 2a), the average mass feed rate really contributing to layer growth depends on melt-pool dimensions (length L and width e). More precisely, small melt-pools will be fed by higher average feed rates (the upper part of the Gaussian distribution) that larger ones. So, an increase of layer heights h with laser power P 0 will be counterbalanced by a decrease of the efficient mass feeding D * , eff m with the melt-pool average size (L + e)/2 (Eq. ( 5)).

D * ,eff m = 2 (L + e) • (L+e)/2 -(L+e)/2 D * m • dx (5) D * m (kg s -1 m -2 ) = average mass feed rate = D m /( • r 2 p ), L (m) = melt-pool length, e(m) = melt-pool width.
The use of a uniform top-hat distribution instead of a Gaussian beam (Fig. 6) is shown to have a rather limited effect on wall dimensions, despite a reduction of thermal gradients in the melt-pool as demonstrated recently by [START_REF] Gharbi | Etude thermique des bains liquides et influence sur l'état de surface des pièces fabriquées par le procédé de fabrication directe par projection laser (FDPL)[END_REF], and a reduction of resulting centrifugal Marangoni flow. Results indicate approximately a 5% increase of wall widths "e". To understand why wall widths are nearly identical, two counterbalancing effects can be assumed: on the one hand, without considering Marangoni flow, a 1.68 mm top-hat heating is expected to induce wider fusion isotherms than a 1.3 mm Gaussian heating, but on the other hand, the reduction of centrifugal Marangoni effects tends to limit meltpool lateral expansion.

If we now consider layer heights h, a tendency to decrease is obtained with a uniform laser irradiation, especially at low scan speeds. This could be attributed to longer and wider melt-pool lengths, that tend to reduce the efficient mass feed rate D * , eff m .

Evolution of surface finishes with process parameters for a

Gaussian beam distribution 3.2.1. General approach Using SEM analysis (Fig. 7a), two distinct contributions to surface roughening were evidenced:

(1) a microscopic contribution coming from (a) particle agglomerations mostly located in inter-layers areas and (b) solidification lines, favored by small melt-pools, and rather low temperature distributions near the melt-pool;

(2) a more macroscopic contribution: the formation of periodic menisci, directly associated with the melt-pool stability (Fig. 7a).

For differentiating macroscopic waviness contribution (W a , W t ) from microscopic roughness (R a , R t ) parameters, a cut-off filter (threshold d th = 80 m) was used (Fig. 7b). The choice of a 80 m filter was justified by the maximum diameter of particles (75 m) that usually provoke micro-roughening when agglomerating on the wall edges. The R a and W a values correspond to arithmetic average values (Eq. ( 6)) whereas R t and W t correspond to peak-to-valley values. All measurements were carried out on a 4 mm-length distance, using 3 to 4 tests for each DMD condition to provide us with statistically reliable data.

R a , W a = 1 d th d th 0 Z(x) dx (6)
3.2.2. Influence of (P, V, D m ) process parameters 3D profiles of surface finish obtained with a Gaussian beam distribution are shown in Fig. 8. For constant mass feed rate and scan speed values, the benefit of using high laser powers and thin additive layers (high scan speeds) is clearly shown. Similarly, on cross-sections (Figs. 9 and 10), the periodic menisci are clearly evidenced, together with melt-pool limits, and layer heights h. The top region of the samples includes the last fusion zone, and the last heat-affected zone where the local material temperature has exceeded the ␤ transus temperature at 910 • C (Fig. 9) as already indicated by [START_REF] Maisonneuve | Profil project: direct manufacturing of aerospace components by laser cladding and laser sintering[END_REF] and [START_REF] Fachinetti | Evolution of temperature during shaped metal deposition: finite element predictions versus observations[END_REF]. This ␣/␤ metallurgical transformation induces successive transition black lines on cross-sections that allow distinguishing layer heights.

An interesting point to notice about cross-sections analysis is the evolution of the melt-pool shape with laser power: for a similar additive layer thickness (0.15 mm) a higher laser power (500 W versus 320 W) promotes a more pronounced internal concavity (quantified by 2 angle in Fig. 9) in the lower part of the melt-pool. This phenomenon is clearly due to an increase of lateral thermal gradients dT/dx and an activation of Marangoni flow (downwards at the edges of the wall, upwards in the center). Another direct effect of the Marangoni flow, is the modification of grain orientations (Fig. 9): a higher concavity of the melt-pool reduces longitudinal grain growths, and favors corolla-like but still-columnar grains structure.

For the specific case of a 0.1 m/min scan speed, thick layers are obtained (near h = 1 mm) and, for moderate laser powers (320 W on Fig. 9), the dilution factor between layers is limited. The increase in laser power (320-500 W) results in deeper and wider melt-pools, higher dilution factors (1 -h/H), and better surface finishes.

Macroscopic (W) and microscopic (R) surface finishes were reported (Fig. 11) as a function of process parameters such as: laser power P 0 (W), scanning speed V (m/s) and interaction distance d i (mm). For a constant average mass feed rate D m , results indicate that:

(1) High scanning speeds promote small waviness parameters W t (Fig. 11b), and small menisci (the mean meniscus height is approximately equal to W t /2: Fig. 10), but have a rather limited effect on micro-roughness data (Fig. 11a); (2) Long interaction distances, corresponding to powder particles that have been molten during their time-of-flight, promote low R t roughness, but do not modify W t values (Fig. 11); (3) High laser powers tend to reduce both waviness and roughness parameters.

Consequently, the best surface parameters were obtained for high scanning speeds and high laser powers. This corresponds to thin and hot liquid layers (high scanning speeds promote short melt-pools and small h values, as predicted by Eq. ( 7)).

The specific behavior of projected particles is mostly responsible for the micro-roughness topography:

(1) First, the melting of particles during their time-of-flight (with high power and long interaction distance) allows them to splat and to expand laterally when impacting the solid wall near the melt-pool, resulting in a reduced contribution to micro-roughness. On the other hand, non-melted solid particles agglomerate on the wall surface with a more or less pronounced necking effect, and contribute to the surface roughening; (2) Second, high P/V ratios (high P, low V) that promote larger meltpools, increase the catch efficiency ratio ϑ (Fig. 5), and limit the number of particles that impact the solid part of the wall and contribute to roughening. This is particularly obvious for short interaction distances where particles are kept at a solid state before reaching the melt-pool (Fig. 11a); (3) Third, a number of particles has been shown (by CCD fast camera) to be only partially melted by the melt-pool, and to form agglomerates on the wall surface. In that case, large and hot melt-pools increase the dilution rate of particles and also improve surface finish.

If we now consider the formation of periodic menisci, the meltpool behavior has to be mostly considered as a key factor instead of the particle behavior. This was made possible by using CCD camera determinations of melt-pool size and dynamics.

Correlation between surface finish and melt-pool sizes

To provide a better understanding to the variation of surface finish with process parameters, a correlation was attempted with melt-pool geometry, using melt-pool data from fast camera analysis (apparent melt-pool height H at the lateral edge of walls: Fig. 4), and final geometries of the walls (layer height h, layer width e). It was found, with a rather good fitting, that meniscus height (W t /2) decrease linearly with high H/ h ratios, and increase linearly with the shape factor of the melt-pool H/e (Fig. 12). This shows that smaller menisci, and more planar sample surfaces can be obtained with wide and deep melt-pools, and with thin additive layers h. This also indicates that large and deep melt-pools are more stable that smaller ones, even if gravity forces F g = gH are more pronounced. In other words, when considering a Ti-6Al-4 melt-pool, capillary forces (F c = 2/R) coming from surface tensions (N/m) are a dominant factor versus gravity-induced melt-pool collapse. This is confirmed by the calculation of the Bond number Bo (Eq. ( 7)) which is the ratio of gravity pressure over surface tension pressure: for a 1 mm height melt-pool, and R = 1 mm radius, Bo ≈ 0.1.

Bo = gh 2/R ≈ 0.1 (7)
The real effect from larger melt-pool volumes may be ascribed to many factors such as: hotter melt-pool surfaces and lower viscosities, provoking larger Marangoni centrifugal velocities at the lateral side, and resulting in a smoothing effect. Consequently, large temperature gradients in the melt-pool, are expected to promote smoother surfaces, due to the enhanced downward liquid flow at the lateral free surfaces of the fusion zone (Fig. 10).

Influence of a uniform laser beam distribution on surface finish

Starting from the conclusions mentioned above, a "top-hat" laser distribution is expected to result in lower temperature gradients in the melt-pools, and lower Marangoni effects. To confirm this, additional DMD trials have been carried out using a 1.68 mm diameter top-hat irradiation, and similar (P, V, D m ) conditions. As shown in Section 3.1, rather low modifications of wall dimensions were demonstrated. Concerning surface finish, the use of a 1.68 mm "top-hat" distribution indicates the following tendencies:

-The optimum surface finish is still obtained with a combination of high power and low scan speeds (high P/V values); -Surfaces are less undulated (factor 2 decrease of W t and W p values) than with a Gaussian irradiation (Fig. 13). This result seems rather contradictory with the assumption of a beneficial effect from negative Marangoni flow (Section 3.2.2) evidenced under Gaussian irradiation. A possible explanation could come from the viscosity dependence versus temperature (Pa s) = f(T(K)). If a large viscosity reduction occurs at high temperature, near the edges of the wall for a top-hat distribution, this should promote enhanced vertical flow, and reduce menisci amplitude. However, = f(T) curves are not easily available in the literature. -The micro-roughness is not modified. This is somewhat consistent with the similar wall dimensions for Gaussian and top-hat irradiations that conduce to similar powder stream/melt-pool interactions.

Analytical description of surface finish

To confirm experimental correlations between melt-pool geometries and lateral menisci, an analytical model was developed, with an analytical description for every melt-pool, considering a similar approach than [START_REF] Fathi | Prediction of melt pool depth and dilution in laser powder deposition[END_REF]. On a 2D cross section, each melt-pool was considered as the sum of two semi-ellipses x 2 /a 2 + y 2 /b 2 = 1 (one for the upper part, and one for the lower part). The analytical formulation of the upper half ellipse was validated by a comparison with the cross section of a manufactured wall (Fig. 14a). Concerning the lower half ellipse, we assumed H 2 = H 1 /n (Fig. 14b), and a concave melt-pool/solid transition due the negative Marangoni flow. The intersection of these two semiellipses, shifted by a h value, at a (x c , y c ) location, allowed us calculating the average meniscus height W p (≈W t /2). The resulting equations (Eqs. ( 9) and ( 10)) confirm that W p increases linearly with layer width e, and also depends on H/ h, as shown experimentally.

W p = e 2 -X c = e 2 1 - 1 - h H 1 (1 -(1/n)) 2 for d dt < 0 (9) W p = e 2 -X c = e 2 1 - 1 - h H 1 (1 + (1/n)) 2 for d dt > 0 (10)
One of the greatest limitations of this simple model comes from the a priori definition of the melt-pool geometry, and the unmodified melt-pool shape when process parameters change. The only adjustable factor is the (n = H 1 /H 2 ) shape factor (Fig. 14b), mostly dependent on thermo-convective flow and on the T = d/dT coefficient, which is kept constant. Considering cross-sections, it was shown that H 2 was very small compared with H 1 . For the numerical application of the model, we considered n = 4, which provides the best fitting with experimental data. The model reproduces correctly the evolution of waviness versus H/ h (Fig. 15) except for H/ h < 2 where the lower part of the melt-pool tends to be flatter and is not properly reproduced by well-established negative thermo-capillary flow. In the near future, a better analytical formulation of the melt-pool geometry, taking into account realistic angles and their evolution versus (P, V, D m ) parameters, would certainly improve the model. As a comparison, the same calculation was also carried out considering a positive Marangoni flow (Eq. ( 9)), and a convex melt-pool shape. In that case, the agreement is better with experimental data at low H/ h values because lower melt-pool curvature becomes less pronounced and n factor increases. Moreover, the elliptic shape assumption is not fully true for the lower melt-pool limit, especially near the lateral edges. This partially explains the difference observed in Fig. 15 between experiments and analytical predictions of W values.

Current 2D thermo-hydraulic calculations of the DMD process, and oncoming temperature measurements in the melt-pool using multi-wavelength pyrometer already used by [START_REF] Gharbi | Etude thermique des bains liquides et influence sur l'état de surface des pièces fabriquées par le procédé de fabrication directe par projection laser (FDPL)[END_REF] will also help us understanding the specific influence of materials' properties ( (N/m), (Pa s), (kg/m 3 )), and temperature gradients near the laser-powder-melt-pool interaction zone.

Conclusions

In this work, a specific focus has been put on surface finish issues, which are one of the main limitations of the DMD process. A distinction has been proposed between micro-roughness and macro-waviness parameters, as two specific contributors to surface modifications. It has been shown that the use of thin additive layers, and large melt-pools improve surface finish, and that increased powder/laser interaction distances, resulting in particle melting are also beneficial factors for surface finish. Last, a correlation has been established experimentally and validated analytically between periodic meniscus height, and melt-pool geometry. The confrontation of these data with recent thermohydraulic calculations will allow us, in the near future, improving the prediction of surface finish versus process parameters, and materials' properties.

Fig. 1 .

 1 Fig. 1. DMD experiments -(a) Experimental set-up and associated diagnostics; (b) detail of the laser-powder-melt-pool interaction zone (H = apparent external height of the melt-pool, h = additive layer height).

Fig. 2 .

 2 Fig. 2. Experimentally determined (a) powder distribution (for Dm = 2 g/min) and laser beam distribution (for P0 = 450 W, and 1.3 mm diameter defocused condition) (b), and associated near-Gaussian analytical formulations (f(x) = A. exp (-5x 2 /r 2 )).

Fig. 3 .

 3 Fig. 3. Temperature distribution T = f(x) inside the powder stream for 2 different interaction distances and powder diameters (based upon the analytical model by Qi et al. (2006), for P0 = 500 W.

Fig. 4 .

 4 Fig. 4. CCD analysis of laser-induced melt-pools: (a) front view with agglomerated particles, (b) side view with floating particles -a 0.45 + /-0.1 m/s centrifugal flow of particles in the melt-pool is estimated, whatever the process conditions.

Fig. 5 .

 5 Fig. 5. Comparison between the analytical calculation of the catch efficiency mass ratio ϑs and experimentally determined data: 15% of the powder projected into the melt-pool does not contribute to layer growth due to ricocheting or non-sticking phenomena.

Fig. 6 .

 6 Fig. 6. Dimensions of DMD walls for different kind of laser irradiations (Gaussian -1.2 mm diameter, Top-hat 1.68 mm diameter spot). Compared to a Gaussian irradiation, the use of a 1.68 mm top-hat irradiation tends to induce a moderate increase of layer widths e, and a tendency to decrease of layer heights h.

Fig. 7 .

 7 Fig. 7. (a) SEM picture of a DMD surface (550 m-height menisci, agglomerated particles accumulated in the inter-layers zones and solidification lines) 320 W -2 g/min -0.4 m/min, (b) 2D profile -Distinction between micro-roughness R and macro-waviness W using a cut-off filter of 80 m from a starting value z (400 W -0.2 m/min -1 g/min).

Fig. 8 .

 8 Fig. 8. 3D profiles of DMD surfaces (Dm = 1 g/min, d i = 1 mm): (a) P = 400 W, V = 0.2 m/min, (b) P = 500 W, V = 0.4 m/min: agglomerates are less pronounced when P/V increases.

Fig. 9 .

 9 Fig. 9. Cross sections of Ti-6Al-4V walls (Dm = 1 g/min, V = 0.4 m/min) after Kroll's etching -(a) P0 = 320 W, (b) P0 = 500 W. The lateral menisci are smoothed by the use of a higher laser power P0, and the melt pools lower curvature angle 2 increases with high P/V lineic energy (J/m) ratios, due to negative Marangoni flow that pushes upwards the central part of the melt-pool.

Fig. 10 .

 10 Fig. 10. Periodic meniscus formation on a Ti-6Al-4V wall (0.2 m/min, Dm = 1 g/min, 500 W), Wt = 120 m -Expected localization of Marangoni thermo-capillary flow.

Fig. 11 .

 11 Fig. 11. Influence of laser power, scanning speed, and interaction distance d i (d i = 1 mm and d i = 3 mm) on: (a) the maximum surface roughness Rt, (b) the maximum waviness amplitude Wt (Dm = constant = 2 g/min).

Fig. 12 .

 12 Fig. 12. (a) Influence of the H/ h ratio on the meniscus formation and resulting surface finish -reduced meniscus heights are obtained for high dilution ratios, (b) Influence of the H/e aspect ratio: meniscus heights are reduced by small H/e ratios.

Fig. 13 .

 13 Fig. 13. Influence of a top-hat laser irradiation on maximum waviness parameter Wt (Dm = 1 g/min) for three different scanning speeds (0.1 m/min, 0.2 m/min, 0.4 m/min).

Fig. 14 .

 14 Fig. 14. Description of the analytical model -(a) identification of the ellipse equation (a "i" exponent = 2 gives a good fitting with experimental cross section considering the melt-pool shape at the top of a DMD wall, (b) basics of the model: the intersection between two half ellipses shifted of a h value allows calculating Wp.

Fig. 15 .

 15 Fig. 15. (a) comparison between the analytical calculation of Waviness Wp (=Wt/2) dependence versus H/ h and experimental data using different thermo-capillary behaviors, (b) scheme of negative and positive thermo-capillary flows for a central and Gaussian laser irradiation, inducing severe transverse thermal gradients.

Table 1

 1 Experimental DMD conditions.

Parameters d laser (mm) P (W) V (m/min) Dm (g/min) Values 1.3

Acknowledgements

The authors would like to acknowledge the financial support of the French National Agency for Research (ANR) under the frame of the ASPECT Project (ANR-09-BLAN-0014).


Thibaud Thénint

Etienne Balmes
email: balmes@sdtools.com

Mathieu Corus

Definition of a linear equivalent model for a non-linear system with impacts

Modal characteristics of non-linear system are typically studied through response to harmonic excitation and using various definitions of non-linear modes. However, few results are available for systems under broadband excitation. The end objective sought here is to generate a linear system, in some sense equivalent to the non-linear system, whose modal characteristics evolve with a level of non-linearity. The considered application is the contact non-linearity found between the tubes of heat exchangers and their support plates. Such tubes, present in nuclear plants, participate to the nuclear safety and can be significantly excited by the fluid flow, so that their dynamic behavior is critical. The turbulent nature of the flow implies broadband excitation and the small gaps between the tubes and the support plate generate very significant non-linear behavior. The proposed equivalent linear system is based on a bilateral contact law whose stiffness and damping characteristics evolve with the amplitude of excitation. A non-linear model is first validated by correlation with experiments. It is then shown that three different indicators (bandwidth of main resonance, operational modal analysis of non-linear power spectral density and correlation of operational deflection shapes) lead to similar values of contact stiffness and damping in the equivalent linear model. This model is thus shown to be a very efficient tool to analyze the impact of the amplitude dependence of the non-linear behavior in the considered system.

Introduction

In Pressurized Water Reactors of nuclear power plants, steam generators act as heat exchangers between primary and secondary coolant fluids. They consist of a bundle of U-tubes in which flows the primary coolant fluid. Several support plates guide these tubes vertically. The secondary coolant fluid flows along the U-tubes and passes through the space between tubes and plates. This space, initially of a foliate shape, is filled with sludge deposits. Consequently, fluid flow is accelerated and the tubes are more excited. Moreover, the mechanical bonding between tube and plate is changed. The combination of these two phenomena can lead to dynamic instabilities and tube cracks. Understanding the non-linear behavior of tubes under broadband excitation but constrained in their motion by small gaps in fairly rigid plates is thus an important industrial application.

Non-linear systems with gaps have been studied for a long time [START_REF] Holmes | The dynamics of repeated impact with a sinusoidally vibrating table[END_REF][START_REF] Shaw | Periodically forced linear oscillator with impact: chaos and long-period motions[END_REF][START_REF] Gontier | Approach to the periodic and chaotic behaviour of the impact oscillator by a continuation method[END_REF][START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF][START_REF] Van De Wouw | Stochastic nonlinear dynamics of a beam system with impact[END_REF][START_REF] De Kraker | Identification of nonlinear phenomena in a stochastically excited beam system with impact[END_REF][START_REF] Sampaio | On measures of nonlinearity effects for uncertain dynamical systems -Application to a vibro-impact system[END_REF][START_REF] Collette | A combined tuned absorber and pendulum impact damper under random excitation[END_REF]. Most work is focused on harmonic excitation and the related notion of non-linear modes [START_REF] Kerschen | Past, present and future of nonlinear system identification in structural dynamics[END_REF], which are fundamentally periodic motion and are not associated with a principle of superposition. The importance of defining a level of nonlinearity is however a key ingredient of these methods that will be reused here. The few applications on broadband excitation do not provide a mechanism to analyze systematically the dependence on level.

The idea of building linear systems in some sense equivalent to the non-linear response has been considered in works on statistical linearization, Ref. [START_REF] Roberts | Random vibration and statistical linearization[END_REF] for example, and equivalent stochastic linear system [START_REF] Bouc | Analysis of multi-degree of freedom strongly non-linear mechanical systems with random input. Part I: Non-linear modes and stochastic averaging[END_REF]. The novelty of the approach proposed here is to give a physical meaning to equivalent parameters and to link their evolution with load level which is a readily accessible quantity.

Section 2 first discusses the testbed configuration and the validation of its properties with an experimental modal analysis. The retained contact model based on an unilateral stiffness distributed over the obstacle height is then introduced. Details on the numerical integration strategy are finally given.

Section 3 then shows that the proposed model correlates quite well with experiment. Features of the Power Spectral Densities (PSD) are then analyzed and shown to be usable as objectives for parametric model updating. The consistence of the updated parameters finally validates the main point of this paper that bilateral contact stiffness and damping, variable with load level, give a very relevant representation of the considered non-linear system under broadband excitation.

Model validation

Configuration of the GV-LOCAL testbed

The GV-LOCAL testbed is composed of a portion of steam-generator tube, clamped at its base. The two straight parts of the tube (legs) are inserted in a quadri-foliate obstacle that has the geometry found in true support plates. These obstacles allow a precise control of the gap between the tube and the support plates. In reality, the gap tends to fill due to oxidation and is thus an unknown and typically evolutive quantity.

Any influence of the fluid would have in true configurations is ignored here. The gap is given as fraction of the tube diameter with a nominal value of 1.3 %, an intermediate value of 0.9 % and reduced gap of 0.23%. The tube is modeled as an Euler-Bernoulli or Timoshenko beam with minor differences between the two models, that are discussed in [START_REF] Thénint | Étude d'un système non-linéaire à chocs sous excitation large bande. Application à un tube de générateur de vapeur[END_REF]. In the testbed excitation is made through a shaker whose position can be modified. This external excitation is random with a frequency content within an adjustable band. Accelerations are measured using a few accelerometers. A finer instrumentation is used around the contact zone (load measurements, accelerometers and Eddie current sensors on one of the support plate modules) allowing an estimation of contact stiffness discussed in section 2.2.

A modal analysis was performed at low amplitudes to obtain a test/analysis correlation of the free tube (without contact on the support plates). Resulting modes are shown in figure 2. Modal damping was 1.5% for the first two modes and of the order of 0.2% otherwise. Very good correlation was obtained after updating of translation and rotation stiffness at the root of the tube [START_REF] Thénint | Étude d'un système non-linéaire à chocs sous excitation large bande. Application à un tube de générateur de vapeur[END_REF].

Contact model

In this application, contact occurs between the solid support plates and the relatively thin generator tubes. When using a full 3D model, one finds visible levels of tube ovalization that cannot be reproduced with a standard beam model based on the assumption of a non-deformable section. To introduce a level of flexibility in the contact, one thus assumes a density of contact stiffness between the neutral fiber of the beams used to model tubes and the support plates, which are assumed rigid. The non-linear behavior is represented by an unilateral stiffness where contact load is null in the absence of penetration and proportional to penetration otherwise. When displaying the penetration/contact load curve in figure 3, it clearly appears that one should distinguish tangent stiffness, which corresponds to the slope of load/penetration curve, from secant stiffness, ratio of instantaneous load and penetration. For the obstacles considered here, a tangent stiffness has been experimentally measured at 4e6 N/m [START_REF] Andriambololona | Caractérisation mécanique du banc GV-LOCAL et présentation de la première campagne de mesure[END_REF] (integrated over the obstacle height). This value is of the same order than the ovalization stiffness of a pipe or than the stiffness used in Ref. [START_REF] Piteau | Vibro-impact experiments and computations of a gap supported tube subjected to single-phase fluid-elastic coupling forces[END_REF].

Under broadband excitation, the average secant stiffness is clearly something that will depend on the load level as it impacts the fraction of time spent with high secant stiffness. The key argument of this paper is thus that using a contact stiffness that depends on load level leads to a meaningful equivalent model.

For numerical applications, the physical model of a continuous elastic layer must be discretized. The number of discrete non-linear springs needed to properly represent the contact dynamics can be validated. Since the proposed equivalent model has variable contact stiffness, the model must be able to represent frequencies for a wide range stiffness values. The evolution of the first frequency with contact stiffness, shown in figure 4, first shows that various ranges can be found. For very low stiffness, one finds the free tube frequency. At intermediate values, the frequency is close to a pinned condition at the support plate center. For high values, one finds a frequency close to a mode clamped on the support plate. The figure then illustrates, that a discretization with 5 levels is a good compromise leading to good accuracy on frequencies (not just the first one shown here) for the whole stiffness range. Finally, to simplify the model, circular obstacles were considered in numerical applications. This was found sufficient for good test/analysis correlation, but one clearly expects an increase of in-plane/out-of-plane coupling due to contact forces that are not perfectly radial in foliate obstacles.

Reduced model for transient analysis

The tube non-linear dynamics are computed solving classical discretized equations of the form

[M ]{q} + [C]{ q} + [K]{q} = F ext (t) + F contact (t) (1)
where M , C and K are the mass, damping and stiffness matrices of the unconstrained tube. The excitation field has been separated between an external shaker excitation and the contact/impact forces located at the five nodes per support plate. To reduce the size of this system, a Ritz basis was used assuming

{q} = [T ]{q R } (2)
and projecting equation (1) (using M R = T T M T , ..., see [START_REF] Balmes | Methods for vibration design and validation[END_REF] among many others for details). To ensure a very accurate reduction the basis was built by creating a mass and stiffness orthonormal basis starting from the following sets of vectors

• free modes with no contact (unconstrained tube) in order to describe motion when the gap is large. While this may not be necessary, a wide band up to 7500 Hz was used based on a spectral analysis of contact loads,

• modes constrained at the support plates to describe modes when full contact is achieved,

• static responses to forces applied at each of the contact points to allow an accurate representation of transients associated with impacts.

The ability to predict modes for a model with bilateral contact and a variable contact stiffness was used [START_REF] Thénint | Étude d'un système non-linéaire à chocs sous excitation large bande. Application à un tube de générateur de vapeur[END_REF] as objective criterion used to evaluate the basis quality and very good results for this objective are obtained.

Both constrained modes and static responses are useful to guarantee proper simulations. Increasing the bandwidth for retained modes is clearly the less important aspect.

Although using a reduced model is not a necessity to achieve reasonable computation times (below the hour here), it has the advantage of easily allowing the introduction of a modal damping model [START_REF] Bianchi | Using modal damping for full model transient analysis. Application to pantograph/catenary vibration[END_REF]. In this case the mass, damping and stiffness matrices are diagonal which accelerates time integration.

The time integration was performed in SDT [START_REF]Structural Dynamics Toolbox 6.4[END_REF], using a standard non-linear Newmark average acceleration scheme with no numerical dissipation. A time-step of 2e -5s was shown to be small enough to avoid energy errors during non-linear time integration.

A characterization of non-linearity under broadband excitation

Model validation through test/analysis correlation

For experiments, a random voltage filtered between 10 and 40Hz is fed to the shaker; the injected force is modulated by the interaction between the tube and the shaker. A more realistic force is thus measured and used as excitation in the numerical simulations. Accelerations are computed or measured at the top of the tube, in the out-of-plane direction. Measurements are performed at a sampling frequency of 2048Hz. Power Spectral Densities (PSD) are obtained through a standard averaging procedure

S xx (f k ) = 2dt 2 M M m=1 |X m,k | 2 (3)
where M windowed signals are extracted using overlap and a window leading to M Discrete Fourier Transform (DFT) X m,k . The results are obviously sensitive to the overlap and windowing parameters but in the present case, total simulation/measurement time of 49s, buffer length of 10s, overlap of 65% and hanning windows gave fairly stable results and were retained. At each frequency, to take into account the variability of each windowed signal, we store the maximal and minimal value over these M DFT and compute the mean value. The thin dotted lines are then used to indicate minimum and maximum amplitude of the PSD.

Figure 5 compares predicted and computed levels for increasing RMS levels of load. This correlation is fairly good over the [10Hz, 200Hz] frequency band. Amplitudes and frequencies match for most peaks and most excitation levels. The spread of the peaks around 90Hz and 170 Hz is well represented. The anti-resonance around 120Hz is both measured and computed. The main differences in levels and frequencies occur for the peaks around 28-30Hz and 70Hz. But despite these limitations the correlation seems particularly good for this very non-linear system. This validates the use of simulations based on the model detailed in earlier sections as a guide for the analysis of non-linear characteristics.

Similarity in behavior for variable stiffness bilateral contact and non-linear behavior

Since the test was only performed at 4 levels of excitation, it is difficult to see whether peaks undergo a continuous transition. As simulations are much easier to reproduce, more detail can be obtained. Figure 6 shows the influence of increasing the load on the power spectral density of displacement at two nodes. One clearly sees a smooth transition from a low to a high level regime. As for linear modal behavior, frequencies of peaks in the PSD do not depend on the response point although point 210 (at 45 o from the top) shows a richer response with more visible peaks. To emphasize the similarity, figure 8 compares the influence of load levels on non-linear PSD to stiffness levels on linear models with bilateral contact. The trends on peak frequencies are clearly very similar. The main differences are a higher peak bandwidth for the non-linear response and a disappearance of antiresonances. In the linear model, the higher peak bandwidth can be modeled as an increased damping.

Methods for equivalence building

Since it was found that the influence of load amplitude and contact stiffness led to similar behavior of the PSD, the objective of this section is to show that similar values of contact stiffness and damping can be found using updating based on characteristics of the non-linear response. One first discusses peak position and bandwidth, operational deflection shapes and broadband PSD, which can be used as objectives for updating then shows the consistence of the resulting updated values.

Identification of the main non-linear peak

Since the response shows resonances that are typical of modal responses, the first obvious path to building an equivalent model is to use modal analysis techniques. In the present case, the excitation is not perfectly characterized so that output only methods are preferable.

The oldest and simplest approach is to seek the peak center frequency and bandwidth. The traditional approach is to use the half power bandwidth to estimate a pole with frequency f j = f max and damping

ζ j = 2
∆f peak f peak but other levels (1/4th or 1/25th are shown in [START_REF] Thénint | Étude d'un système non-linéaire à chocs sous excitation large bande. Application à un tube de générateur de vapeur[END_REF]) give similar results. From the estimated pole, updating of contact stiffness and damping is readily obtained since the dependence is very smooth for the first mode, which is isolated here.

Of course much more robust OMA methods can be used. Methods derived from the work of Verboden [START_REF] Verboven | Frequency-domain system identification for modal analysis[END_REF] were implemented and used here. Figure 9 clearly shows very similar pole estimates for the OMA and simple peak bandwidth estimates. The evolution of frequency and damping with the amplitude of excitation is easily analyzed. At low levels the number of impacts is very small and the frequency is that of free system. At high levels, the behavior is similar to bilateral contact. Between 0.1N and 10 N, a transition occurs with increasing frequency and an associated higher level of apparent damping. For linear systems, the spread of the resonance peak can only be due to damping. For the non-linear system being considered, there is no physical dissipation mechanism so this spread is not damping, although it can be consistently identified using a linear damped pole model. Rather than interpreting the spread as damping, one can realize that one can represent the system as linear with variable contact stiffness. The peak spread thus comes from averaging of a time varying resonance.

Given a pole estimate associated with the first peak, contact stiffness and damping are readily selected to obtain the same pole for the first pole of the equivalent linear system.

Shape correlation on the operational deflection shape

Modes are characterized by a pole, which was used in the previous section, and a shape which will be used here. When dealing with linear modal analysis, one of the fundamental principles is that the mode shape is invariant and generates similar spatial responses at all frequencies around the resonance. In other words the responses mostly occur within a subspace of rank one associated with the modeshape.

The singular value decomposition is the most common method to extract subspace information from an arbitrary response. Here one extracts equivalent modeshape information by computing the SVD of the PSD over a frequency band around the main peak in the form

[H(ω)] N ×N w = min(N,N w) j=1 {U j } N σ j {V j } N w (4)
where one analyzes the response at N locations (here the horizontal out of plane motion) and N w frequency points (here the frequencies in the frequency band above 1/25 t h of the main peak). Figure 10 shows the first singular vector U 1 for various load levels. At low loads, the ODS clearly corresponds to the first out of plane bending mode. As the load increases a vibration node appears and its position rises with increasing load up to the last shape which can be associated with the second free/free mode but also, with more pertinence, with a bilateral contact mode associated with a stiffness to be determined through model updating using the MAC between mode and ODS as an objective function.

The shapes are mostly real although for intermediate values of load the animation clearly shows some dephasing between various points. The analysis of singular values, in figure 11, also shows that for load values in the transition region more than one singular value contribute to the response. This emphasizes the fact that the response is really close to being linear in the extreme low and high load levels and much more complex for intermediate values.

Broadband correlation of PSD

The broadband PSD is another response feature of interest. Focusing on a wider frequency band is important if the first peak is not the only one of interest. Figure 12 illustrates how within the possible linear responses an adjustment of contact stiffness and damping can be used to obtain a response that is fairly close to the non-linear PSD.

In this search for the optimum, relating the RMS load level in the non-linear response with an appropriate input spectrum for the linear response is a difficulty. For the proposed updating of contact stiffness and damping, the correlation coefficient between the vectors of PSD at multiple frequencies and a single response point was used. This objective is clearly related with the MAC which compares responses at multiple points and a single frequency. It shares with the MAC the insentivity to vector scaling which was the desired feature.

Equivalent contact parameters obtained through updating

Each of the characteristics of the non-linear discussed in previous sections can be used as a model updating objective resulting in estimated contact stiffness and damping (in reality densities of stiffness and damping which are then distributed over the support plate height). The estimates for each of the objectives are shown in figure 13. As also true in linear modal analysis, it appears that damping cannot be estimated from the ODS shape. But otherwise the three estimates are surprisingly consistent.

Conclusion

Modeling of non-linear contact with impacts is now clearly accessible for simulation with sufficient time lengths to allow pertinent estimation of PSD. The model of this study was shown to have a very decent test/analysis correlation thus demonstrating a strong physical validity.

The impact of variable load levels was shown to be very similar to that of a bilateral contact layer with variable stiffness and damping properties. This similarity was further demonstrated by updating equivalent contact parameters using different non-linear features, which resulted in very close estimates of the parameters.

These results validate the proposition that the considered non-linear system with impact can be represented by a linear system with contact stiffness and damping that change with the load level. This equivalence is based on power spectral density and operational deflection shape information. Compared to the current

Figure 1 :

 1 Figure 1: a) Global view of the GV-LOCAL testbed. b) detail of the support plate module c) sensor configuration of the experimental modal analysis

Figure 2 :

 2 Figure 2: Experimental modeshapes of the free tube

Figure 3 :

 3 Figure 3: Left : contact model as a discretized density of contact stiffness. Right : tangent and secant stiffness.

Figure 4 :

 4 Figure 4: Evolution of the first modal frequency with respect to the global contact stiffness.

Figure 5 :

 5 Figure 5: PSD of measured (green) and computed (blue) acceleration at top in the out-of-plane direction. The vertical lines indicate modes in free free (black) and clamped (red) boundary conditions

Figure 6 :

 6 Figure 6: PSD levels from non-linear time simulation as a function of frequency and RMS load level. Levels are normalized with the load level.

Figure 7 :

 7 Figure 7: PSD of DOF 210x for amplitudes between 0.01 and 20 N RMS.

Figure 8 :

 8 Figure 8: Comparison of (left) non-linear PSD as a function of load amplitude and (right) linear PSD as a function of contact stiffness

Figure 9 :

 9 Figure 9: Left : estimation of peak center frequency and bandwidth. Right : comparison of poles obtained through OMA identification and analysis of peak bandwidth.

Figure 10 :

 10 Figure 10: Operational Deflection Shape (ODS) associated with the first singular value at various load levels

Figure 11 :

 11 Figure 11: Singular values as functions of load level

Figure 12 :

 12 Figure 12: Target non-linear PSD (black), optimal linear PSD (green) and range of linear PSD (thin lines) for two load levels.

Figure 13 :

 13 Figure 13: Equivalent stiffness and damping parameters as function of load level

practice of using linear systems with one or two boundary conditions, validating the response for a family of models with variable contact parameters seems a readily accessible and very relevant option.

Equivalence for other criteria like maximum stress is not obvious and should be analyzed in more detail. The gap is another parameter characteristic of the non-linearity. Validating that gap reduction has an effect similar to load level increase and mapping the gap/load/contact parameter relation would be an obvious extension of this work. Accounting for fluid and fluid/elastic instabilities, or applying the proposed equivalence to other non-linearities seem useful perspectives.


G Vermot Des Roches

O Chiello

E Balmes

X Lorang

Benchmarking Signorini and exponential contact laws for an industrial train brake squeal application

Contact representation of structure interactions for finite element models is nowadays of great interest in the industry. Two contact modellig strategies exist in the literature, either based on a perfect contact with no interpenetration of structures at contact points, or based on functional laws releasing the contact constraint through pressure-penetration relationships. Both strategies require very different and rarely documented numerical implementations, making difficult any objective comparison. This paper presents a benchmark between ideal contact and a functional law of the exponential type applied to squeal simulations by complex mode analysis of an industrial railway brake.

Introduction

Increase of computational power and deployment of more and more efficient solvers dramatically increased simulation capabilities in the field of structural dynamics. Industrially, simulating interactions of multiple components in mechanical assemblies is thus of particular interest nowadays -this commonly requires contact and possibly friction models. The particular application of this paper concerns brake squeal for trains during station parking operations. Noise disturbances can raise issues in the vicinity of stations over a few kilometers, limiting train exploitations. The AcouFren project piloted by SNCF (the french railway company) thus aims at providing noise indicators at early brake design stages, through simulation.

Contact modelling strategies are widely documented in the literature, from which two trends are clearly found. The oldest representation is due to Signorini from the 30's, it idealizes a perfect contact between two ideally smooth bodies, and constrains the relative displacement of contact points. The second representation strategy is functional and authorizes a level of interpenetration between bodies. First measurements were performed in the 60's, from which Greenwood and Williamson [START_REF] Greenwood | Contact of nominally flat surfaces[END_REF] developed models where contact stiffness is function of the bodies interpenetration. State-of-the-art measurements as performed by Nogueira et al. [START_REF] Nogueira | Experimental validations of elastic to plastic asperitybased models using normal indentations of rough surfaces[END_REF] still show that this model can be relevant, mostly for hard surfaces. Direct ulstrasonic measurement of contact stiffnesses, can also be performed. Biwa et al. in [START_REF] Biwa | Stiffness evaluation of contacting surfaces by bulk and interface waves[END_REF] thus showed that contact stiffness is function of the contact pressure, mostly for low interfacial loading.

Physically speaking, contact between two bodies can be interpreted as the harmony between their respective stiffness and asperity compression, as illustrated in figure 1. When one material is much softer than the other, like for brake systems, local compression of asperities becomes non-negligible at the macroscopic scale, generating local interpenetration of the nominal surfaces (with sometimes plastic asperity deformation). It can then be relevant to identify a law p(g) linking overclosure of nominal surfaces and contact pressure.

Figure 1: Representation of contact between two bodies, body stiffness and asperities Numerically speaking, implementation of functional contact laws using continuous non linear functions of displacement is easier since classical non-linear algorithms can then be used. The ideal contact law does not feature any mechanical parameter to be identified, but its numerical implementation requires specific algorithms. Their convergence is often very sensitive and fine tuning of numerical parameters is needed. Besides, numerical convergence is often obtained regarding a tolerance, such that interpenetration could still occur in the results, with levels that are difficult to control.

Choosing a contact strategy is thus difficult and the general feeling is that these different contact models yield different results. The aim of this paper is to present an objective comparison of Signorini and exponential functional contact laws in application to an industrial train brake squeal simulation. Section 2 starts by presenting numerical implementation details of functional contact, based on the work of Vermot des Roches [START_REF] Vermot | Frequency and time simulation of squeal instabilities[END_REF]. Signorini implementation is then presented, based on the work of Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF]. Section 3 then aims at demonstrating that an equivalence between both contact laws is possible, and a functional law of the exponential type is proposed. Section 4 eventually presents a validation applied to industrial railway brakes.

Contact formulations for squeal simulations

Squeal simulations in the frequency domain is nowadays widely deployed in the industry, and is based on the sliding perturbation method. This was used by Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF], Lorang [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant: Application au crissement des freins de TGV[END_REF] for train brake squeal applications, or by Vola et al. [START_REF] Vola | Friction and instability of steady sliding: squeal of a rubber/glass contact[END_REF] to study rubber/glass instabilities in sliding steady states, and by Vermot des Roches in [START_REF] Vermot | Frequency and time simulation of squeal instabilities[END_REF] for automotive brake squeal application. Contact formulations and adapted numerical schemes for squeal simulations are here presented.

Contact handling for finite element models

The representation of contact between two solids as formalized by Signorini is defined by a contact direction N along the outward normal at a contact point of one of the bodies (the chosen one, S1 is commonly called master, as contact forces will be evaluated on its surface). Figure 2 illustrates the contact configuration. The distance between the two bodies along the contact normal is called the gap, here noted g. The convention is that the gap is positive when contact is open, or not effective. When contact is closed, or effective, a contact pressure p exists between the bodies, applied at contact points, equal in amplitude on both solids. The contact law proposed by Signorini is illustrated in figure 2, it defines the condition of non-interpenetration as

     g ≥ 0 p ≥ 0 (g).(p) = 0 (1)
Alternatively, functional contact laws define a pressure-gap relationship, they are sometimes exploited as numerical regularization means, or for physical reasons, using experimental characterization. Classical laws are illustrated in figure 3. Numerically the gap is a linear observation of relative displacement between two surfaces along a normal, and can thus be written using an observation matrix [C N OR ]. Noting N the number of system DOF, and N c the number of contact points, it is written with a possible gap offset

{g 0 } {g} Nc×1 = [C N OR ] Nc×N {q(u)} N ×1 -{g 0 } Nc×1 (2)
The force resulting from the gap-pressure relationship is then defined at each contact point as

{q} T {f N } = Γ {û(q)} T N pdS ≃ e j
{û(q)} T {N } p(x j , q)ω (e) j J (e) (x j )

where f N is the global contact force, p the contact pressure, q a virtual displacement, q the displacement, x (e) j are the integration points of current element e, J (e) (x j ) the Jacobian of the shape transformation (surface associated to each integration point) and ω (e) j the weighting associated with the integration rule of element e. In practice contact forces can be recovered at nodes using the gap observation (2):

{f N } N ×1 = [C N OR ] T N ×Nc ω (e) j J (e) (x j ) p(x j , q) Nc×1 (4)
The Coulomb law is expressed as

{f T } < µ {f N } ⇔ {w} = {0} {f T } = µ {f N } ⇔ {f T } = µ {f N } {w} {w} ( 5
)
where f T is the friction force, µ the friction coefficient, and w the sliding velocity. Computation of the sliding velocity requires computation of differential velocities between slave and master surfaces in the plane orthogonal to the contact normal. A tangential displacement observation matrix [C T AN ] can be used

{w} 2Nc×1 = [C T AN ] 2Nc×N { q} N ×1 (6)
and friction forces can be recovered, from the local friction force f tj by

{f T } N ×1 = [C T AN ] T N ×2Nc ω (e) j J (e) (x j ) f tj (p(x j , q)) 2Nc×1 (7)
For 3D solids, Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF] noted that the sliding friction tangent state also features a damping term due to the possibility for the friction force to change its orientation in vibration. Physically, the higher the contact pressure, the more difficult the sliding direction variation. This is due to the priviledged sliding direction of the steady sliding state. A planar friction damping matrix [C f 2 ] derived from the expression of friction forces thus appears, written

δf T = [C f 2 ] δ q = [C T AN ] T 2 \ 1 w i µf N i \ [C T AN ] 2 δ q (8)

Functional contact

For functional contact formulations, classical non-linear methods can be applied, since contact-friction forces only depend on the system state. Given the system elastic stiffness [K el ], and a displacement vector {q}, the system must verify

[K el ] {q} = {f ext } + {f N L (q)} (9)
where the non linear force vector can be decomposed as {f N L (q)} = {f N (q)} + {f T (q)}, with {f N } the contact forces and {f T } the friction forces.

Functional contact statics

The most classical method to resolve regular non linear statics is to use the Newton method, presented in figure 4. Stakes in Newton resolution algorithms are mostly on the numerical side. Although mathematical convergence is obtained theoretically when the system is regular enough, hard laws like contact implementations are challenging.

Initial state : q 0 Residue computation The system Jacobian [J] derived from (9) can be written

r(q k ) = -[K el ] q k + f c + f N L ((q) k ) Correction [J] {∆q} = r k q k+1 = q k -∆q Convergence r k+1 fc ≤ ǫ no
[J] = [K el ] + [K c (g)] (10)
Considering a steady sliding state, the tangent friction perturbation is null. Noting ω j J(x i ) is the surface associated to each contact point. the tangent contact state defining [K c (g)] writes

∂f N (q + δq) ∂δq = [C N OR ] T \ ω j J(x i ) \ ∂p(q + δq) ∂δq [C N OR ] = [C N OR ] T \ k ci (q) \ [C N OR ] (11)
For hard laws tangent state variation can be very brutal, which is a drawback of such algorithms. Indeed, when contact opening occurs in a loaded area at iteration k, tangent contact state becomes locally null, and correction can become very large (with possibly contact overclosures). Divergence can then quickly occur as the non-linear forces returned to the large overclosures can overcome numerical conditioning and machine capabilities (infinite forces).

A full set of adaptations can be applied to the original Newton scheme, yielding numerous variations of quasi-Newton methods. For hard contact, cut-backs on corrections based on the mechanical residue norm, and controlled release on the contact part of the Jacobian for opening points are usually working and have been used here. Implementation details is key here to performance, but are outside this paper scope.

Functional contact sliding perturbation formulation

Functional contact formulations are defined as function of the gap, such that a tangent state can be derived for each contact point by deriving the gap/pressure relation, illustrated in figure 3. A key aspect to note is that for a similar contact pressure level, different functional contact laws can yield very different contact stiffnesses. Contact law calibration in statics and dynamics can then be uncorrelated.

Definition of tangent contact stiffness has been defined in [START_REF] Balmes | Vermot des Roches[END_REF]. Defining tangent friction states is less direct as it couples normal and tangential directions, which is typically non-symmetric. Friction expression ( 5) indeed shows that a variation of contact force has a direct effect on friction force, whereas friction force variation may occur without effect on contact force.

For sliding states, friction force is explicitly defined and only depends on contact force at the same point.

The tangent friction stiffness is then the tangent contact stiffness [START_REF] Balmes | Vermot des Roches[END_REF] scaled by the friction coefficient µ, and transfers normal displacement to planar force. Since planar displacements are free for a fixed contact state, this effect is not considered for real mode computation. The tangent sliding friction coupling stiffness K slide nlf is thus defined as

K slide nlf = 0 for real mode computation [C T AN ] T \ µk ci (q) \ [C N OR ] for complex mode computation (12)

Signorini contact

Ideal contact verifying the exact Signorini condition ( 9) is non regular, hence the common name of nonsmooth dynamics. {f N L (q)} has no specific link to the displacement, and the pair ({q} , {f N L }) has to be directly solved.

Signorini statics

Resolution strategies are usually based on Lagrange resolutions. The non-symmetric nature of friction requires a methodology refinement, using Augmented-Lagrange methods. De Saxcé and Feng developed the bipotential method [START_REF] Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF]. Alart et al. presented a generalization of the Newton method adapted to Lagrange formulations introducing Gauss-Siedel like algorithms [START_REF] Jourdan | A gauss-siedel like algorithm to solve frictional contact problems[END_REF].

The resolution chosen here was developped by Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF]. It is based on the status assumption of each contact points (open/closed gap). Closed contact points generate a displacement constraint under which resolution is linear. Noting with superscript (e) the observation restricted to conctact points in effectively closed contact, one resolves the problem in an orthogonal subspace T (e) defined as

{q} = T (e) {q e } / [C N OR ] (e) {q e } = 0 (13)
Equation ( 9) can thus be written

T (e) T [K el ] T (e) {q e } = T (e) T ({f ext } + {f N } + {f T }) (14)
Since contact constraints are verified, contact forces are known

{f N } = [C N OR ] (e) ([K el ] {q} -{f ext }) (15)
and are null in T (e) . Sliding friction forces (5) can thus be written

{f T } = µ [C T AN ] T [C N OR ] ([K el ] {q} -{f ext }) (16)
which is a little simplified in T (e) . Equation ( 14) becomes

T (e) T [I] -µ [C T AN ] T [C N OR ] [K el ] T (e) {q} = T (e) T [I] -µ [C T AN ] T [C N OR ] {f ext } (17) It must be noted that the factor [I] -µ [C T AN ] T [C N OR
] is not symmetric. Direct resolution can then be challenging for large problems, like the million-DOF system presented in section 4. Use of the latest PARDISO librairies [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with pardiso[END_REF], implemented in SDT [START_REF] Balmes | Vermot des Roches[END_REF] was here required.

Contact status must be predicted, and can be updated depending on the result of (17). This is formalized by considering the group of contact nodes considered in effective contact in (17) and verifying Signorini (1). This group noted C k is defined by

     ∀x j ∈ C k , f N (x j ) ≥ 0 and ∀x j / ∈ C k , g(x j ) ≥ 0 (18)
At each iteration, contact points violating Signorini conditions are classed in group D k such that

     x j ∈ C k , f N (x j ) < 0 or x j / ∈ C k , g(x j ) < 0 (19)
Status updating allows forming a new [C N OR ] (e) and consists in switching the status of points in D k ,

C k+1 = (C k -C k ∩ D k ) ∪ (D k -D k ∩ C k ) (20)
Iterative status correction can be resolved by an Uzawa algorithm, presented in figure 5. This method convergence has been proved without friction and by switching contact statuses one by one. Its performance is however very satisfying in this paper application, with very little iterations even for complex systems and high friction coefficients.

Prediction :

{f N } 0 , {f T } 0 Mechanical resolution [K el ] {q} k+1 = {f c } -{f T } k -{f N } k Correction          {f N } k+1 = [C N OR ] T [K el ] {q} k+1 {f T } k+1 = µ [C T AN ] T {f N } k+1 Contact Status update Convergence ?
Signorini verified at all contact points ? no Figure 5: Implementation of an Uzawa resolution algorithm for the Signorini-Coulomb contact-friction laws

Signorini modal analysis

Following the exact contact condition of Signorini, the system tangent state only on the contact surface, materialised by C ∞ , as presented in 2.3.1. The mechanical problem is here formulated as

λ 2 [M el ] + λ ([C el ] + [C f 2 ]) + [K el ] {q} = {f N L } [C N OR ] (e) {q} = 0 (21)
Like for static resolution, projection in the subspace orthogonal to the contact constraints T (e) yields

λ 2 M (e) µ + λ C (e) µ + [C f 2 ] + K (e) µ {q e } = 0 (22)
where matrices X (e) µ are the so-called sliding matrices, defined as

X (e) µ = I -[C T AN ] T [C N OR ] T (e) T [X el ] T (e) (23)

Definition of an exponential contact law

The two contact implementation methods presented in section 2 are numerically very different, and require different resolution algorithms. The concept of figure 1 is here exploited to show that functional laws can be calibrated to provide equivalent results to the ideal law. This section thus aims at presenting a relevant type of functional contact law and its numerical calibration.

The clamped/sliding block

The railway industry pads can be of different shapes and materials, a sample one named G35 is presented in figure 6. A friction block (made of friction material) is extracted for this section illustration. The block is here solidar to a metallic ring making the link to the backplate. The friction material has a non-negligible loss-factor, modelled by a material Rayleigh damping (easy compatibility for transient simulations). When computing steady sliding states, and sliding perturbations, the friction blocks effective contact area is supposed to be fully sliding (although some contacts can be opened). A sample demonstration of the blocks behavior is thus to consider them in clamped/sliding condition. Ideal contact condition is materialized in this case by clamping normal displacements of the friction block bottom side. Figure 7 shows some modes associated to the pad block in such conditions. To demonstrate the effect of contact stiffness (possibly derived from a functional contact law), normal spring elements of parametered stiffness are added to the blocks bottoms instead of clamping. The previous case would then correspond to a spring with infinite stiffness. Mode damping evolution as function of contact stiffness can be observed in figure 8. Damping rates increase with frequency as a consequence of the stiffness proportionality coefficient used to form the material Rayleigh damping matrix.

An interesting behavior can here be seen when contact stiffness is located of the frequency curve inflection area, seen in figure 8. For low contact stiffnesses, ground coupling is small and the friction block behaves as free, that-is-to-say the whole strain energy is located in the block. Damping ratio of the clamped/sliding block is thus linked to the damping ratio of the block itself. The same observation can be done for high contact stiffnesses, when coupling is strong and the friction block behaves as clamped. The expected damping ratio increase is obtained regarding the frequency increase, materialized bythe dotted line in figure 8.

In the curve inflection area, a non negligible strain energy is located in the coupling springs, so that a non negligible part of the total system strain energy is outside the friction block. The damping ratio thus varies less than the linear evolution expected, and can even decrease when coupling becomes non negligible while still being low. This observation is very general regarding impact of components strain energy distribution in damping of assemblies. More detailed illustrations of this effect can be found in [START_REF] Bianchi | Using modal damping for full model transient analysis[END_REF].

Exponential contact law definition

There exists a large variety of functional contact laws that could be calibrated, the exponential type only is retained here. This choice is more fonded on the expected transient behavior, where contact opening transitions and penetration saturation effects also need to be properly approached for a potentially large range of pressures. Numerical results obtained in [START_REF] Vermot | Frequency and time simulation of squeal instabilities[END_REF] with such law were satisfying.

Given {g} the gap vector (the opposite of the overclosure vector), obtained from an observation of the relative displacement of both surfaces along the local contact normal, pressure {p} is defined as

{p} = p 0 e -λ{g} (24)
The exponential law has two parameters p 0 and λ, whose effect on contact behavior is illustrated in figure 9.

To interpret these curves, values have to be observed as function of a reference force level, conceptually named F . Horizontal intersections are thus iso-values regarding system loading outside the contact area.

Parameter p 0 only drives the curve offset, thus the level of gap obtained for a given force. Identification of such parameter can be obtained if measurements of surface rougthness are accessible, providing acceptable values of interpenetrations. Coefficients are only numerically calibrated in this paper so p 0 can be arbitrarily set. Choosing a p 0 such that significant contact pressures is obtained for closed (or negative) gaps, simplifies implementation and allows easier physical interpretations.

-0.01 -0.005 0 0.005 0.01 0 Parameter λ impacts contact stiffness at a given force level, the higher, the stiffer. This parameter itself can thus be updated to choose a contact stiffness that complies with the fact that the block frequencies in clamped/sliding conditions are converged towards the behavior in perfect contact. Looking at the diagrams of figure 8, target contact stiffnesses values can be pin-pointed in figure 10. For calibration, the level of contact pressure generated by the system must be evaluated, either by using a status implementation or by choosing a linear contact law for a preliminary computation, or by exploiting the system command (the case here). Noting p 1 the obtained reference, contact pressure expression and the contact stiffness k c (g) (first order derivative of contact pressure) yield

p 1 = p 0 e λg k c (g) = λp 0 e -λg (25)
The product λg can be substituted in the first equation of ( 25) such that

λ = k c (g) p 1 ( 26
)
Choosing a target value of 10 7 M P a in figure 10, and knowing that the contact force command on the TGV brake system is here of 5kN over a surface of 21.4 10 3 mm 2 , one obtains λ = 4.28 10 7 mm -1 . Higher values could also be used for security margin regarding system variability. It must however be kept in mind that values with too high stiffness may alter numerical conditioning of the stability problem.

For comparison means, a converged value kc = 10 5.3 M P a at the limit, and a non converged value kc = 10 4 M P a will also be tested. Correspondances between λ and target stiffnesses are reported in table 3.2.

Application to industrial railway brakes

The french high speed train, TGV, features 4 disc brake systems on each axle, as presented in figure 11. A single brake system is here modelled by finite elements on which subassemblies (the disc, rig and pads) have been updated. The full model is free-meshed using second order ten nodes tetraedrons, yielding between 500,000 and 1,000,000 DOF depending on the fitted pad model.

The mesh is kinematically positioned depending on the pad thickness to have zero gap at the origin. Disc sections underlying the pads are remeshed to obtain compatible contact interfaces. Braking force is applied by an actuator at the system rear, which translates the pad holder towards the disc. The translation is kinematically realized by the rotations of levers and rods, constituting the brake rig. This paper application focuses on complex mode evaluation for station parking operations. Comparisons between calibrated exponential laws, defined in section 3.2 and the ideal contact implemenation are presented for steady state solutions in section 4.1 and complex modes in section 4.2.

Steady sliding simulations

A contact force command of 5kN is applied by the actuator. Figure 12 presents global steady sliding solution with ideal contact implementation. Levels of displacement and contact force resultants are presented. Largest displacements are seen by the lever with values aournd 140µm. Contact force distribution greatly varies over the contact surface. Friction blocks of the front end are the most loaded, unloading of the trailing end seems complicated and is function of the pad holder fixations.

Figure 12: Steady sliding state absolute displacement of a TGV brake using Signorini contact implementation, contact resultant at contact nodes of the disc top side, and strain energy density distribution Static results are compared between Signorini contact and the family of exponential contact laws defined in table 3.2. In the following, contact nodes define nodes of the mesh in the contact area, while Gauss contact points define the Gauss contact integration points used for functional contact implementation.

Figure 13 presents contact results in terms of pressure and gaps. Points with opened gaps at contact nodes are presented in figure 13a. It can be seen that no global difference rises for these displacements. Although contact pressure from an exponential law is strictly positive, exponential decay as function of opening is sufficient not to alter behavior of unloaded contact points.

Figure 13b presents gap results for closed contact. Great differences can here be observed, the softer the law, the deeper the penetration. Penetration levels seem reasonable regarding global displacements for converged λ values only. Figure 13c presents contact pressure at contact nodes. Contact pressures at the interface are here a posteriori computed from displacements to allow comparison between implementations, using a nodal resultant to nodal pressure operator. Opened contact points will thus show negative pressures at the interface, which would correspond to depressions on a membrane. Such observation shows that pressure levels seen by the structure outside the interface are very similar. Forces transmitted by the interface are naturally recovered in all cases, which is physically natural.

Global results in figure 13c can be put into perspective with the local results of figure 13d. Contact pressures are here plotted at the Gauss contact points for the functional laws only. Clear differences occur for low pressure areas, with a convergence pattern as function of λ. Such fluctuations for very low pressures are however not percieved by the structure.

Figure 14 presents absolute differences in displacement between the Signorini response and functional contact laws. It can be seen that for all cases, maximum differences scale in nanometers, a thousand times smaller than displacement scale. Larger differences can here occur, with very limited strain energy errors for converged cases. Athough static displacement seems better for the middle law, contact stiffness distribution is not as good, as discussed in section 4.2.

Sliding perturbation results

Complex modes are here directly computed based on the formulations presented in section 2. Stability diagrams are presented in figure 15 for 500 modes, providing a frequency bandwidth of 0 -10kHz. The global stability diagram of figure 15b shows the effect of material Rayleigh damping, with very large damping ratios on the first modes (mass factor) and an average damping ratio linear increase with frequency (stiffness factor).

The first four most unstable modes are plotted in figure 15a and15c. Modes 27 and 69 are pad/disc coupling modes, with effective displacement of the pad holder. Modes 81 and 59 are planar modes, instability coming from friction block modes (possibly coupled with in-plane disc modes). No unstable modes are found for high frequencies in this configuration -Rayleigh damping coefficient stills needed experimental fitting in this study's model. One of the main interest of functional contact laws is the representaton of a contact stiffness distribution depending on contact pressure, whose variation can be non negligible for low pressures [START_REF] Biwa | Stiffness evaluation of contacting surfaces by bulk and interface waves[END_REF]. This is impossible with Signorini contact, where this pattern can only be assessed by evaluating contact opening thresholds as function of complex mode amplitudes.

Globally, stability diagrams fit relatively well for all exponential calibrations. Larger differences are obvious for the soft value, highlighting a convergence as function of λ. Zooms in are provided in figures 15d, 15e and 15f. Structure modes (with clear displacement of the brake rig) are in the 2kHz range, and show correct fittings, although the soft law presents some peculiarities.

Modes in the 5.5kHz and 9kHz ranges mostly feature friction block modes. Slight differences between converged λ values and Signorini are mostly due to contact stiffness variation of low contact pressure areas.

A very good stability diagram fitting is here only obtained for the hardest value of λ. Comparison between unstable modes with MAC is eventually presented in figure 16. Due to the large amount of modes computed and model size, each mode basis weights over 3 GBytes, making direct comparisons difficult. Unstable modes are of higher interest due to their propensity to generate squeal, hence the choice to compare only these shapes, with a detection threshold set under a damping ratio of -10 -3 %.

Convergence as function of λ is well observed regarding the number of detected unstable modes and their shape correlation, figure 16c displays a squared matrix with unitary diagonal. More unstable modes with no representation in the Signorini results are found for the two softest λ values.

Conclusion

Simulation of structure-structure interaction is now at stakes for industries, requirering performant implementation of relevant contact models. For the two dominant contact modelling strategies found in the literature, efficient implementation is available, but lacks benchmarking.

Physically, transmission of contact forces between structures should naturally show correct contact pressures independently from the contact law. Ideal laws avoids identification of physical parameters, but also lacks representation of contact stiffness variation as function of contact pressure. Functional laws require identification from difficult to realize experimentation. It was however shown that numerical calibration is sufficient and can be performed by considering the system apparent stiffness seen at interfaces.

The application presented in this paper, concerning squeal simulation of industrial brake squeal with complex mode analysis in the sliding perturbation framework, allowed benchmarking ideal and exponential contact. The advantage of the exponential law over other types is the quick decay of contact forces for opened gaps, while keeping relevant saturation patterns. Satisfying results are thus obtained for frequency and transient simulations. The objective comparison obtained confirmed that pressures transmitted between structures are independent from the contact strategy.

Implementation choices can thus be performed considering experimental capability and performance, with an insurance of results convergence. A clear perspective is a benchmark extension to transient simulations, available with methods such as presented in [START_REF] Vermot | Frequency and time simulation of squeal instabilities[END_REF]. Such simulations are there so intensive that choice can only be directed by numerical performance.

Figure 2 :

 2 Figure 2: Contact normal (N) and gap (g) definition between solids, Signorini and Coulomb laws

Figure 3 :

 3 Figure 3: Classical contact laws, and their derivatives (for contact stiffness)

Figure 4 :

 4 Figure 4: Implementation of a static Newton resolution scheme

Figure 6 :

 6 Figure 6: A sample railway braking pad. The friction blocks in orange are linked to the backplate in green, using rings in blue. A dovetail joint (in red) fixes the pad in the brake rig.

Figure 7 :

 7 Figure 7: Sample block modes in clamped/sliding conditions with ideal contact for the G35 ring fitted block.

Figure 8

 8 Figure 8 presents complex mode frequency and damping results of the clamped sliding block as function of contact stiffness, an S-shaped frequency evolution can be observed. Saturation occurs around 10 4 M P a, after which frequencies are equivalent to an ideal contact implementation.

Figure 8 :

 8 Figure 8: Frequency (left) and damping (middle) evolution as function of the contact stiffness applied to the lining block of the G35 pad. Right: zoom in on the first modes.

1 -

 1 p0 = 10 4 , λ ∈ (10 2 ; 10 10 )m -λ = 1.3.10 9 , p0 ∈ (10 2 ; 10 7 )P a

Figure 9 :

 9 Figure 9: Exponential contact law evolution as function of its parameters

Figure 10 :

 10 Figure 10: Identification of target contact stiffnesses of the clamped/sliding block, and identified λ

Figure 11 :

 11 Figure 11: The TGV brake system

Figure 13 :

 13 Figure 13: Comparative gaps and contact pressures between Signorini and exponential laws resuls.

Figure 14 :

 14 Figure 14: Absolute displacement differences in nanometers between Signorini and functional static results, from top to bottom along the x, y and z axis, from left to right in increasing λ. Strain enrgy densities of differential displacements are presented for elements with highest levels only

Figure 15 :

 15 Figure 15: Stability diagrams of Signorini solutions and functional law solutions. Most unstable TGV modes

  λ = 4.28.10 7

Figure 16 :

 16 Figure 16: Unstable modes (ζ ≤ -10 -3 %) MAC between Signorini and exponential laws solutions