Hanno Niemann
email: hanno.niemann@tu-bs.de

Joseph Morlier
email: joseph.morlier@isae.fr

Amir Shahdin
email: amir.shahdin@isae.fr

Yves Gourinat
email: yves.gourinat@isae.fr

Damage Localization using Experimental Modal Parameters and Topology Optimization

This work focuses on the developement of a damage detection and localization tool using the Topology Optimization feature of MSC.Nastran. This approach is based on the correlation of a local stiffness loss and the change in modal parameters due to damages in structures. The loss in stiffness is accounted by the Topology Optimization approach for updating undamaged numerical models towards similar models with embedded damages. Hereby, only a mass penalization and the changes in experimentally obtained modal parameters are used as objectives. The theoretical background for the implementation of this method is derived and programmed in a Nastran input file and the general feasibility of the approach is validated numerically, as well as experimentally by updating a model of an experimentally tested composite laminate specimen. The damages have been introduced to the specimen by controlled low energy impacts and high quality vibration tests have been conducted on the specimen for different levels of damage. These supervised experiments allow to test the numerical diagnosis tool by comparing the result with both NDT technics and results of previous works (concerning shifts in modal parameters due to damage). Good results have finally been archieved for the localization of the damages by the Topology Optimization.

Introduction

In recent years the use of fibre composite materials in aeronautical structures has vastly increased. Due to their superior characteristics, concidering e.g. the specific tensile strength, composites are used for a variety of lightweight structures, recently even for pressurized airplane fuselages [START_REF] Tsai | Introduction to Composite Materials[END_REF]. Unfortunately, composite structures show a very complex mechanical behavior concerning dynamic loads and a variety of damage mechanisms, that are hard to classify and to predict. Some of these damages are fibre or matrix cracking, fibre matrix debonding and ply delaminations [START_REF] Kashtalyan | Mechanisms of internal damage and their effect on the behavior and properties of cross-ply composite laminates[END_REF].

Composite laminates are susceptible to damages from a wide variety of sources which include fabrication stress, environmental cyclic loading and foreign object impact damage [START_REF] Gude | Fatigue failure criteria and degradation ruls for composites under multiaxial loadings[END_REF]- [START_REF] Gadelrab | The effect of delamination on the natural frequencies of a laminated composite beam[END_REF], and which may lead to severe degradation of the mechanical behavior due to the loss of structural integrity. Therefore, it is even more important to understand the creation and evolution of damages in composites (damage detection), to identify affected regions (damage localization) and to evaluate their influence on the structure as a whole (damage classification and quantification). These steps are the cornerstones of structural health monitoring (SHM) [START_REF] Yan | Development in vibration-based structural damage detection technique[END_REF]- [START_REF] Zou | Vibration-based model-dependent damage (delamination) identification and health monitoring for composite structures -a review[END_REF].

In recent years, structural health monitoring using vibration based methods has been rapidly expanding and has shown to be a feasible approach for detecting and locating damage. A detailed and comprehensive overview on the vibration based techniques has been presented in references [START_REF] Doebling | Damage identification and health monitoring of structural and mechanical systems from changes in their vibration characteristics: A literature review[END_REF]- [START_REF] Carden | Vibration based condition monitoring: A review[END_REF].

The basic principle of vibration based damage detection can be explained as follows. Any structure can be considered as a dynamic system with stiffness, mass and damping. Once some damages emerge in the structures, the structural parameters will change, and the frequency response functions and modal parameters of the structural system will also change. This change of modal parameters can be taken as the signal of early damage occurrence in the structural system. Vibration-based structural damage detection refers in this context to detection methods for structural damage using only the structural characteristics, such as natural frequencies, modal damping, mode shapes, etc. Structures can be excited by ambient energy, an external shaker or embedded actuators. Accelerometers and laser vibrometers can be used to monitor the structural dynamic responses. A variety of broadband excitation signals have been developed for performing shaker measurements with FFT analyzers, e.g. burst random, burst chirp, etc.. Since the FFT provides a spectrum over a band of frequencies, using a such an excitation signal makes the sprectral measurements much faster than using sine dwell or swept sine excitations [START_REF] Schwarz | Experimental modal analysis[END_REF].

Change in natural frequencies is the most common parameter used in the identification of damaged regions [START_REF] Gadelrab | The effect of delamination on the natural frequencies of a laminated composite beam[END_REF], [START_REF] Shahdin | Correlating low energy impact damage with changes in modal parameters: A preliminary study on composite beams[END_REF]- [START_REF] Richardson | Correlating minute structural faults with changes in modal parameters, Proceedings of Spie[END_REF]. Various such methods that use natural frequency information are reported by Salawu [START_REF] Salawu | Detection of structural damage through changes in frequency: a review[END_REF]. The advantage of using the change of structural natural frequencies to detect damage is its convenient measurement and high accuracy. However the measurement of natural frequencies may not provide enough information for structural damage detection to relate the changes to a correct damage location. Furthermore, natural frequencies are often not sensitive enough to initial damage in structures. Therefore, other damage indicators should also be considered in the damage localization process, like damping or anti-resonances [START_REF] Zhang | Relation of damping and fatigue damage of unidirectional fibre composites[END_REF].

As listed in [START_REF] Doebling | Damage identification and health monitoring of structural and mechanical systems from changes in their vibration characteristics: A literature review[END_REF], several authors have already worked on the idea to correlate damages with a degradation of structural stiffness and/or mass in numerical simulations, which consequently changes the dynamic properties of the numerical model. Adapting this idea, a new approach, using Topology Optimization design variables for localizing damages, has been published [START_REF] Lee | Damage detection by the topology design formulation using modal parameters[END_REF]. The key of this approach is that, due to the character of Topology Optimization, a search for defects is performed globally over the entire structure, which distinguishes this method from earlier proposed methods and offers interesting prospects.

The original goal of Topology Optimization was to find, within a defined discretized solution domain, a structure of minimal compliance (highest rigidity) by connecting and seperating trusses or insertion and removal of holes in the structure, meaningly by changing the topology. This approach had first been proposed in [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF] based on a homogenization of microstructural elements with rectangular holes whose size are defined by the design variable. The dimensions of the holes could vary from 1 (void) to 0 (solid), normalized to the element's size, and also from element to element. By homogenization of the discrete elements, the composite microstructure could be transformed in an equivalent homogeneous material. Since this method had shown to be rather cumbersome and to deliver hazy results for the optimized structure, a slightly different approach has been presented in [START_REF] Bendsoe | Optimal shape design as a material distribution problem[END_REF], which is often referred to as the Power Law approach. Here, the design variables are assumed to be an additional element property that can be understood as a relative density of the element. By reducing the stiffness and mass of certain elements with their newly assigned density fraction property, a local change of structural stiffness and mass can be obtained.

Concerning the damage localization approach, this Topology Optimization variable is now applied to the design domain of a Finite Element model that is based on the undamaged structure. By optimizing the system's stiffness and mass matrix towards those of the damaged structure by matching modal parameters, the correct location and geometry should then theoretically be found in terms of elements with a lower density. In this case, the localization of damage is done by estimating the most probable equivalent damage (local loss of rigidity) which leads to a minimization of the norm between the baseline (undamaged) FRF and experimental FRFs of different damage states. This is the principle idea behind the set up topology optimization formulation for damage localization described in the mentioned article in ref. [START_REF] Lee | Damage detection by the topology design formulation using modal parameters[END_REF].

In the here presented work, this previous method is implemented in a widely-used commercial code (MSC.Nastran), adding important preliminary tests by computing different Objective functions for fractional mass penalization and validating the results with experimental tests on composites. Thus, from the numerical methodology in [START_REF] Lee | Damage detection by the topology design formulation using modal parameters[END_REF], we propose a direct application to replace NDE with damage localization by topology optimization with codes, validations and results.

More precisely, in the cited work, the design responses have obviously been calculated using a commercial finite element solver, whereas the topology optimization is performed independently for each iteration loop. The here presented damage localization approach is translated into a Nastran input file, since MSC.Nastran is already equipped with a topology optimization routine as well, which is a supposingly more advanced and more reliant considering the optimization results.

Finally, an experimental validation of the damage localization approach is also tried out on a CFRP specimen by localizing barely visible damage due to low-energy impacts which is an important topic in SHM. At the moment, the diagnostic approach is however still limited on homogeneous isotropic materials due to restrictions in Nastran. Therefore it has been tested on homogenized quasi-isotropic composite beam models. For acquiring the necessary modal data for the damage localisation, vibration tests have been performed on such composite beams. The results offer interesting perspectives considering the damage detection in composite materials and the ultimate goal would be the development of an automatic diagnosis tool as a mean of SHM using only experimental modal data.

Numerical Estimation of Modal Parameters

Mechanical Background

In this section, the theoretical background for the numerical approach used in this work is briefly presented. Basic equations for modal analysis and the calculation of frequency response functions for a discretized system are given. Also, the optimization problems statement for the presented damage localization method is developed.

Eq. 1 shows the set of equations of motion in matrix notation for a random discretized structure. Hereby, the system is considered to be discretized by Finite Elements and M, C and K are the system's mass, damping and stiffness matrices, respectively. The vector ( f ) is a timedependent load vector.

M ü + C u + Ku = f(t) (1)
The analyses consist of the two principal steps, which are displayed in fig. 1. A modal analysis is performed to determine the principal dynamic characteristics of the system, like natural frequencies and mode shapes. From these modal parameters the response of the structure to excitation can be calculated. This is similar to solving the homogeneous and particular solution for the system of differential equations.

Figure 1: Principal steps of dynamic analyses for structural systems, where x is the displacement vector as a response to the excitation force vector f. H is the matrix of frequency response functions, where Ω is the excitation frequency

Modal Analysis

For calculating the natural frequencies, the system is considered to vibrate harmonically, damping and external forces are neglected. By replacing the displacement vector and its derivative by a sinusoidal solution function, the general eigenvalue problem can be formulated (eq. 2).

K -ω 2 M u 0 = 0 (2)
Assuming the existence of the vector u 0 , meaningly a non-trivial solution of this set of equations, the system's determinant has to vanish. This leads to the characteristical equation displayed in eq. 3.

det (K -λM) = 0 (3)
The solutions of this equation are the n eigenvalues for a system of n degrees of freedom, where each eigenvalue λ i corresponds to a natural angular frequency ω i as stated in eq. 4.

λ i = ω 2 i with i = 1, 2, ..., n (4)
The amplitude vector corresponding to an eigenvalue is the eigenvector or mode shape vector ψ i . The combinations of natural frequencies and mode shape vectors describe preferred states of vibration of the structure [START_REF] Ewins | Modal Testing: Theory, Practice and Applications[END_REF]. For the numerical estimation of the eigenvalues, several algorithms are available, e.g. the Lanczos method [START_REF] Hörnlund | Analysis and measurments of vehicle door stru dynamic response[END_REF], which is also available in MSC.Nastran.

Modal Frequency Response Analysis (Harmonic Analysis)

The cornerstone of modal frequency response analyses is the diagonalization of the system's matrices by a modal transformation. In a seperation approach the primary variables are replaced by a linear combination of the so-called modal coordinates and the mode shape vectors, as shown in eq. 5, where Ψ is the modal matrix and q is the vector of modal coordinates.

u(t) = ψ 1 q 1 (t) + ψ 2 q 2 (t) + ... + ψ n q n (t) = Ψq(t) (5)
Due to the orthogonality properties of the mode shape vectors, inserting the approach into eq. 1 and premultiplying the equation by Ψ T leads to eq. 6, which are the decoupled equations of motion for the analyzed system. For diagonalizing the damping matrix by the mode shapes vectors from the modal analysis, a proportional damping approach [START_REF] Ewins | Modal Testing: Theory, Practice and Applications[END_REF] has to be assumed.

Ψ T MΨ q(t) + Ψ T CΨ q(t) + Ψ T KΨq(t) = Ψ T f(t) (6)
The system now consists of n equations of motion for single degree of freedom oscillators in form of eq. 7, where m i is the modal mass, c i a modal damping coefficient, k i the modal stiffness and p i a modal force.

m i qi (t) + c i qi (t) + k i q i (t) = p i (t) with p i (t) = ψ T i f (7)
The equation can be normalized by the modal mass, so that the coefficients are expressed in terms of natural angular frequencies ω i and modal damping ratios ζ i , as they are defined in eq. 8. This leads to eq. 9.

ω 2 i = k i /m i and ζ i = 1 2 c i √ k i m i (8) qi (t) + 2ζ i ω i qi (t) + ω 2 i q i (t) = p i (t) m i (9)
Eq. 9 is an ordinary differential equation of second order that has the solution stated in eq. 10 which can be derived by a complex approach where Ω is the excitation frequency. Here, the variable q has already been transformed in the frequency domain which is noted by the circumflex.

qi (Ω) = 1 ω 2 i -Ω 2 + 2 jζ i ω i Ω • pi m i (10)
From here the matrix of frequency response functions (FRFs) can easily be derived, since the frequency response functions are defined as the displacement response functions multiplied by the inverse of the excitation force vector (eq. 11).

H(Ω) = û(Ω) f-1 (Ω) = n i=1 ψ i qi (Ω) f-1 (Ω) = n i=1 ψ i 1 ω 2 i -Ω 2 + 2 jζ i ω i Ω • ψ T i m i (11)
For a certain degree of freedom k, the frequency response function due to an excitation force at the l-th degree of freedom can than be calculated by eq. 12.

H kl (Ω) = n i=1 ψ ik ψ il m i ω 2 i -Ω 2 + 2 jζ i ω i Ω (12)
If the natural frequencies and damping ratios are known, this equation can be evaluated for discrete excitation frequencies.

Experimental Procedures

Material and Specimens

The goal is to provide the data which is necessary for performing and validating the damage localization approach presented later. Several specimens have been fabricated and have undergone vibration testing in the undamaged and damaged state [START_REF] Shahdin | Correlating low energy impact damage with changes in modal parameters: A preliminary study on composite beams[END_REF], [START_REF] Shahdin | Significance of low energy impact damage on modal parameters of composite beams by design of experiments[END_REF].

Resin-containing carbon-fiber/epoxy prepregs of T300/914 are used to fabricate the test specimens. The material is supplied by Hexcel composites, the physical properties are set out in tab. 1. The specimens are processed in a press. The curing cycle of the laminates is 2h at 180 • C with a warming-up cycle of 0.5h at 135 • C. The laminates are cut into beams using a diamond wheel cutter, following the ASTM D3039/D3470 standards. The specimens have a thickness of 3mm and consist of 24 plies. The lay-up is chosen as such as the delamination is said to have more profound effects on the dynamic characteristics [START_REF] Saravanos | Effects of delaminations on the damped dynamic characteristics of composites[END_REF]. The vibration tests are carried out with two steel masses attached at the ends. The aim of putting these masses at the ends is to enhance the difference in the modal parameters between the undamaged and the damaged test specimens [START_REF] Vanhoenacker | The use of multisine excitations to characterise damage in structures[END_REF].

Vibration Tests

The experimental equipment used for vibration testing is shown in fig. 2b. The experimental set-up is that of a free-free beam excited at its center, based on the Oberst beam method [START_REF] Wojtowicki | New approach for measurements of damping properties of materials using oberst beams[END_REF]. This method states that a free-free beam excited at its center shows the same dynamical behavior as that of a half length cantilever beam. The test specimen is placed centered on a B&K force sensor (type 8200), which is then assembled on a shaker supplied by Prodera, having a maximum force of 100 N. Furthermore, a fixation system is used to connect the test specimens to the force sensor. This fixation is glued to the test specimen with a HBM X60 rapid adhesive. Using a laser vibrometer of type OFV-505, provided by Polytec, the response displacements are measured with high precision and no contact. The shaker, force sensor and the laser vibrometer are manipulated by a data acquisition system running with LMS.Testlab software. Before testing, small reflecting stickers were glued to each of the measurement points on the composite specimens to ensure a good reflection of the vibrometer's laser beam (fig. 2a).

The center of the test specimens is then excited by a burst random signal. The excitation point is marked as point 17 in fig. 3. The signal is averaged 10 times for each measurement point and the chosen frequency band is 0 -1600Hz. Hanning windows are used for both the output and the input signals. The influence of excitation level is checked and a high frequency resolution (∆ f = 0.25Hz) for precise modal parameter estimation is used. While the level of the excitation signal is chosen as 1N, the response is measured at 33 points, which are symmetrically spaced in three rows along the length of the beam. The modal parameters are extracted by frequency domain parameter estimation methods (Polymax) integrated in the data acquisition system.

The composite beams are analysed in the three different states: An undamaged state (UD), a first damage state due to four impacts (D1) and a second damage state after eight impacts (D2). The locations of the impacts are also shown in fig. 3. Vibration tests are carried out at each of these three states. A simple case is studied where the impact points are chosen as such as the damage is symmetrical on both sides of the two axes of symmetry After the estimation of the set of frequency response functions, the modal parameters can be determined. This is done by decomposing the frequency response functions into the discrete modes of the system by calculating the poles of the FRF [START_REF] Brown | Parameter estimation techniques for modal analysis[END_REF]. An enhanced version of this identification is automatically performed by the LMS.Testlab software by using the PolyMAX method. Tab. 2 shows the experimentally obtained natural frequencies for a composite beam in the undamaged state and for the damage states D1 and D2 after impacts with 8J.

Furthermore, fig. 4 shows the experimentally obtained frequency response functions for one measurement point on the surface of the composite laminate that has been impacted by 8J. A remarkable shift in resonances due to impact damage damage can be noticed, especially for higher frequencies.

Optimization Problem Formulation for Damage Localization

For formulating general optimization statements suitable for computation, some general definitions have to be established. In the beginning, an objective of the optimization has to be identified as a function of the design variables that are changed in the course of the optimization. Since the principle objective of the performed optimizations is the matching of a set of modal parameters, any kind of matching function (e.g. least square formulations) can be considered. In the present case, due to the nature of the BigDOT optimization algorithm [START_REF] Vanderplaats | Very large scale optimization[END_REF], which is based on exterior point penalty functions, a pseudo objective function formulation is constructed where the parameters, which are to be matched, are included as constraints. By the automatic performance of convergence checks for each of the frequency constraints, the number of considered design responses is practically dynamically adapted in every optimization cycle, depending on the fulfillment or non-fulfillment of the constraints. In preliminary tests, this formulation also generally gave more precise results compared to alternative formulations based on finding a search direction from minimizing the least square of all responses together, for example.

Such an exemplary pseudo objective function is shown in eq. 13, where f is the original objective function, and P a penalty function consisting of a set of inequality or equality constraints, which are g and h, respectively. The parameter r p is a further penalization coefficient that can be applied to the constraints. min :

Φ(χ) = f (χ) + P(r p , h(χ), g(χ)) (13)
The vector of design variables χ consists of the Topology design variables corresponding to the Power Law approach [START_REF] Bendsoe | Optimal shape design as a material distribution problem[END_REF].These design variables are principally an additional element property that can be understood as a relative density of each element e as stated in eq. [START_REF] Shahdin | Significance of low energy impact damage on modal parameters of composite beams by design of experiments[END_REF].

χ e = ρ e ρ 0 subject to 0 ≤ χ min ≤ χ ≤ 1 ( 14
)
Since the variable is normalized by the original densities, it can only assume values between 0 and 1, as stated in the side constraints. Due to numerical reasons, 0 is replaced by a minimum threshold. This design variable is then multiplied by the corresponding element's stiffness and mass and thereby alters the properties of the element (eq. [START_REF] Tracy | Effect of delamination on the natural frequencies of composite laminates[END_REF]. E e = χ p e E 0 and m e = χ q e m 0 [START_REF] Tracy | Effect of delamination on the natural frequencies of composite laminates[END_REF] The design variable is also penalized, which is supposed to help in getting a clearer solid-void solution by making intermediate design variables more "costly". The values of the penalization exponents p and q are problem-dependent, but a penalization factor of about 3 is generally proposed in common literature [START_REF] Bendsoe | Topology Optimization -Theory, Methods and Applications[END_REF].

In the following, the constraints consisting of modal parameters are defined, where the primed value (ω) always denotes the reference data of the damaged structure, and the plain value (ω) corresponds to the current data of the optimized model. Since the goal is to minimize the difference between these values, a proximity is usually defined. The first set of constraint equations g ω requires the n ω chosen angular resonance frequencies to be within the proximity ω of corresponding resonance frequencies of the damaged structure (eq. 16).

g ω i = ω 2 i -ω 2 i 2 < ω for i = 1, 2, ..., n ω (16)
This formulation has a global character, since the resonant frequencies are global parameters for the structure. In [START_REF] Lee | Damage detection by the topology design formulation using modal parameters[END_REF], it is shown that local modal data is necessary for archieving better results, which is why the use of anti-resonances is proposed. Whereas resonances can be found at the peaks of the frequency response functions, their local minima are called anti-resonances (fig. 5). These minima result from the additional contribution of two neighboring modes with a 180 degree phase shift [START_REF] Avitabile | Experimental modal analysis -a simpel non-mathematical presentation[END_REF]. Anti-resonances are unique for any point on the structure and are therefore considered as local data. The squares of the natural angular frequencies are used for a closer correspondance with the governing equations of this subproblem, namely the undamped system's equations of motion in eq. ( 2). Due to limitations in Nastran, anti-resonances cannot be directly considered as constraints, but nodal frequency response functions, so that anti-resonances can still be included in the optimization by evaluating the FRFs at the corresponding frequency. Thus, the second used set of constraint equations (g H ) requires that the magnitude value of the nodal FRF (FRF of an FE node) at a certain excitation frequency Ω j is within the proximity of that of the damaged structure at the same frequency [START_REF] Richardson | Correlating minute structural faults with changes in modal parameters, Proceedings of Spie[END_REF]. As before, the constraints are formulated based on the governing equations in eq. [START_REF] Carden | Vibration based condition monitoring: A review[END_REF]. It has however to be noted that the frequency response functions have to be transformed in polar form and only the magnitude is yet used as a constraint.

g H j = |H kl (Ω j )| -|H kl (Ω j )| 2 < H for j = 1, 2, ..., n H (17)
Since preliminary simulations with the programmed input file delivered unsatisfying results, a further penalization function is implemented as an additional objective in the optimization. As in these first tests not only the density of the elements in the anticipated region was reduced, but the density of random elements all over the structure, which then led to a unreasonable reduction in global structural mass, the additional function had the objective to penalize the structural mass and thereby guide the optimization towards anticipated results. Furthermore, since the fibre composites do not actually lose mass due to the impacts, it is desired to keep the structural mass close to the original mass, for example. Still, the mass penalization has to be adjusted to the given damage type or material, as shown in the paragraphs dealing with the actual case studies.

Since the design variables are the elements' fractional densities, it is more practical to use also fractional masses for the objective functions, as defined in eq. 18, where m 0 is the mass when the design variables is 1 (solid).

m f ractional = n e m e m 0 (18)
Here, n e is the number of elements in the design domain. Due to the constant design volume, this fractional mass is equal to the sum of all fractional densities as shown in eq. 19. Several objective functions have been tested, which are laid out in tab. 3.

m f ractional = n e V 0 ρ e V 0 ρ 0 = n e χ e (19)
L = m f rac (1 -m f rac ) L = 1 -2m f rac (2) L = |1 -m f rac | L = -1 for m f rac < 1 (3) L = (1 -m f rac ) 2 L = -1 + 2m f rac
In fig. 6 it can then be seen that all functions show a relative minimum at the right end of the domain where all fractional densities are 1, but they possess different derivatives, which are important for the sensitivity analysis and basically characterize the penalization. The derivative of function ( 2) is constant over the whole domain ([0;1]), the derivatives of the other functions are variable with the fractional mass. Whereas the derivative value of function (1) approaches the value of the derivative of function (2) close to the right boarder, the derivative of function (3) approaches zero. Therefore, similar results have been obtained using functions (1) and (2) during the simulations. Also, they generally gave better results than function (3).

Sensitivity Analysis

Since mathematical programming methods usually involve the calculation of gradients of the objectives for the sensitivity analysis, the derivatives of the constraints with respect to the design variables are presented in this section. Eq. 20 gives the gradients of the resonance frequencies that can be derived from eq. 2 in the frequency domain form as shown in reference [START_REF] Bendsoe | Topology Optimization -Theory, Methods and Applications[END_REF], where ψ ie and ψ je are the components of the eigenvectors with respect to the nodal degrees of freedom of the element e.

∂ω 2 i ∂χ e = ψ T ie ∂E e ∂χ e -ω 2 i ∂M e ∂χ e ψ je (20)
The partial derivatives of the frequency response functions have been derived in [START_REF] Ma | Structural topology and shape optimization for a frequency response problem[END_REF]. Based on eq. 11, the frequency response function can be written as eq. 21. Then, the gradients of the frequency response functions with respect to the element density fraction χ e are calculated as in eq. 22.

H(Ω) = Ψ q(Ω) f-1 (Ω) (21)

∂H(Ω)

∂χ e = ΨS q(Ω) f-1 (Ω)

Here, S is the sensitivity matrix as defined in eq. 23, where ω i is the ith natural frequency and ζ i the corresponding damping ratio. As before, Ω is the exciting frequency.

S i j = -1 ω 2 i -Ω 2 + i2ζ i ω i Ω ψ T ie ∂E e ∂χ e -Ω 2 ∂M e ∂χ e ψ je (23)
In both cases the gradients of the stiffness and mass matrices can be calculated as in eq. 24. These calculations are principally performed in an enhanced form in Nastran.

Set-up of Topology Optimization Formulation for Damage Localization in Nastran

This section presents the Nastran input file parameters for the damage localization approach. The models that are used for a numerical validation of this method, as well as the modeling of the specimens described in the previous sections are also explained. The latter models are supposed for an experimental validation of the damage localization method.

MD Nastran seemed therefore suitable for an implementation of the localization approach by Lee et al. [START_REF] Lee | Damage detection by the topology design formulation using modal parameters[END_REF], as with the release of the version MSC Nastran 2005 r2, the embedded optimization libraries were enhanced by a topology optimization capability and the introduction of the BIGDOT optimizer. Especially this new optimization library provided the capability to deliver practical results dealing with the large amount of design variables that are common in topology optimization. Thus, only Nastran was available for an "all at once" implementation of the described method combined with low calculation time.

For the optimization of the structure, several new cards have to be added to the input file. At first, the design variable has to be defined. In the case of topology optimization, this is explicitly done by the card TOPVAR. The structure of this card is exemplary shown in tab. 4 with the corresponding symbols from section 4. Here, type of property characterizes the elements, that are included in the optimization, concerning their deformation behavior. This definition has to be in accordance with the elements of the defined design region, which is called by the Property ID that has been assigned prior to the corresponding elements.

Then, the chosen modal parameters of the structure that build up the constraints and objectives have to be included in the calculations. In Nastran, these parameters are referred to design responses and they are accessed by the DRESPi entries. For each desired structural response, such a DRESPi entry has to be defined, where the index i denotes the type of design response formulation.

The design responses can then be constrained by referencing the corresponding DRESPi entry in a DCONSTR command and by defining the boundaries there. All the responses and constraints can be called, either by themselves or as a linear combination, by the global objective function (DESOBJ) or the global constraint command DESGLB in the case control section, or, as subcase dependent constraints, by a DESSUB command for each subcase. Figure 7 shows how responses, constraint and objectives can principally be connected. More precise definitions of each of the commands are given in [START_REF]MSC Nastran 2007 r1 -Users Guide for Topology Optimization[END_REF] and [34].

In each optimization cycle, two subcases have to be evaluated: the modal analysis and the modal frequency response analysis. Each subcase thereby calls its own set of design constraints via the DESSUB command. An exemplary input file can be found in the appendix A, where the design constraints and corresponding design responses, like the model itself, are sourced out into external files to keep input file mostly model-independent.

Numerical Modeling

Several numerical models are proposed to test the feasability of the described damage localization method. These models basically consist of two types of specimens, a simple cantilever beam with isotropic material properties for numerical validation and the model of the composite beam for experimental validation.

The procedure for numerical validation is to create an undamaged Finite Element model, copy it and modify the copy by deleting elements at certain points of geometry. This way, a local loss of stiffness can be generated due to the deletion of elements. Modal frequency response analyses are then performed on this template model to obtain the modal parameters of the undamaged and the damaged case. Hereby, the damage localization can be tested on numerical models that lack an experimental base but show sufficient validity for amelioration of the method prior to the performance of experiments. This process has also been performed in the article by Lee et al. [START_REF] Lee | Damage detection by the topology design formulation using modal parameters[END_REF].

A cantilever beam test model with homogeneous isotropic material properties is shown in fig. 8a, which was meshed with iso-parametric membrane-bending shell elements. This figure also shows the applied constraints and load. On the free end of the beam a nodal frequency-dependent excitation force with constant amplitude was applied. The experimentally tested beams have also been modeled by shell elements. Due to the symmetry of the beams and the vibration test setup, it was sufficient to model only a half beam with the corresponding boundary constraints at the intersection of the the two halves [START_REF] Wojtowicki | New approach for measurements of damping properties of materials using oberst beams[END_REF]. For the nodes at this intersection the displacements in direction of the excitation force were allowed, 14 whereas all other degrees of freedom had to be constrained. Furthermore, the connecting piece and the end masses that were glued to the specimens for the vibrations tests had to be explicitely modelled by three-dimensional solid elements. A distributed load with a combined amplitude of 1N has been applied to the bottom of the connecting piece to simulate the excitation force of the shaker pot. Since, at the moment, the topology optimization feature in Nastran is still limited to isotropic materials only, the application of the code is also limited to thin quasi-isotropic laminates, for which a plane stress state assumption can be made. Mean stiffnesses can then be calculated based on the laminate theory [START_REF] Tsai | Introduction to Composite Materials[END_REF], reducing whole laminate to a two-dimensional shell. Using material properties calculated from such a homogenization of the laminate, the composite beams can be modeled by membrane-bending shell elements as shown in fig. 8b. As mentioned before, damping is not considered in the updating process, since there is no response for damping values implemented in the optmizition module of Nastran. Therefore, a low constant modal damping value of 0,3% has been chosen.

Still, the simple model suceeded to reproduce the dynamic response sufficiently accurate for the present case, since the experimental measured FRF for a randomly chosen point on the surface of the beam fits well with the numerical calculated FRF of the corresponding node (fig. 9).

Results and Discussion

Numerical Validation on Models with isotropic Material Properties

The Nastran input file has then been applied to the described models. Figure 10 shows such an obtained simulation result in terms of the optimized density distribution. Here, the first ten natural frequencies and the resonance and anti-resonance peaks of the FRF at the excitation point in a bandwith from 0 to 500Hz were used as matching data. Hereby, the optimization responses were desired to be within a randomly chosen proximity of 5% ( < 0.05) of the target value and the penalization mass was set 97.5%, which is the actual fraction of the number of elements of the damaged to the undamaged structure (ten elements of a total of 400 elements were deleted for creating the cut out in the damaged template model 10). Although, this test case is based only on numerical models, where the location of damage and the loss in mass was known in advance, the approach can be considered as reproducible with physical specimen. Since, in the actual case of a loss of mass due to impacts in physical specimen, the mass target fraction can be determined by weighing the undamaged and damaged structure. If no difference in mass can be measured, the mass target has simply to be set almost equal to 100%, as it is done with the later presented composite beam.

The density distribution resulting form the optimisation (fig. 10a) is then in quite good concordance with the goal geometry (fig. 10b), since the brighter regions are regions of lower densities that are supposed to represent the damages. However, the calculation on isotropic materials is performed only for a reproduction of similar results as in [START_REF] Lee | Damage detection by the topology design formulation using modal parameters[END_REF] and validation of the method's adaption in Nastran. Hereby, the application of this method to localize such evident damage as the cutout in the aluminium beam is not really reasonable, but it becomes all the more interesting when used for localizing non-evident damage in composite beams, as presented in the next sections.

Experimental Validation on composite beams

The optimization algorithm has then been tried out on the models of the composite specimens. Here, the data used for constraints were the natural frequencies and FRFs that had been obtained by the vibration tests of the damaged specimens (cp. tab. 2 and fig. 4). Since no actual cutout was present, that would justify a penalization mass target as in the numerical validation, the mass target was set to 99%. In the case of the damaged composite specimen, it is reasonable to keep the mass as close the original mass as possible. However, a slight reduction in mass has to be included, allowing some changes in element density and the most suitable mass target value has yet to be determined depending on the problem, which is a limitation of the method at the moment.

In the following, some Topology Optimization results are presented for the modal parameters of the beam that was impacted with an energy of 8J. Thereby, two sets of damage state modal parameter data existed for the beams (D1 and D2). Thus, topology optimization was performed with the goal to either localize the first four impact sites or all eight. It has to be reminded, that the beams were considered to be symmetrical, thus each side of the beam has undergone two and four impacts, respectively. The results for the both cases using only resonance frequencies as constraints are shown in fig. 13 and 14 with the corresponding ultra-sound scan. For the first two impacts the right locations have effectively been found with similar sizes as in ultra-sonic images that have been produced after vibration testing (fig. 14; the outer two impact damage sites have been introduced at first (D1), the two inner impacts later (D2, see fig. 3); ultra-sonic testing could only be performed after all other testing had been completed.)

For the more damaged beam with four impact sites only a rather uncertain region of damaged material with lower densities could be located (fig. 14), but not the discrete impact points as shown in the ultra-sound picture, since the damage is very closely spaced and is difficult to differenciate. However, the displayed region of lower densities is still close to a damage region geometry if the delamination regions due to the impacts, that are shown in the ultra-sound image, were merged. Still, for both cases, the FRFs and resonance frequencies show a relatively good match with the damaged cases (fig. 15 andtab. 6).

From the obtained results it can be concluded that their quality depended vastly on the choice of the optimization parameters. These optimization sensitive parameters can be divided in two principle groups. The first on consist of the parameters related to the optimization process, i.e. penalization exponents and objective functions, and the other depends on the set of chosen constraints and quality of the experimental results.

This choice of design constraints is not very intuitive. Since higher resonance frequencies show usually a higher shift due to introduced damages, it seems obvious to use modal parameters of a large spectrum. But with higher frequencies the necessary modal identification becomes more difficult as higher modes are more susceptible to noise and participate less and less in the actual dynamic response. The presented Topology Optimization approach is also limited to modal parameter shifts due to loss in stiffness and mass, but does not take into account damping, which has proven in experimental results to be more sensitive to damage than stiffness changes ( [START_REF] Shahdin | Correlating low energy impact damage with changes in modal parameters: A preliminary study on composite beams[END_REF], [START_REF] Shahdin | Significance of low energy impact damage on modal parameters of composite beams by design of experiments[END_REF], [START_REF] Montalvão | A review of vibration-based structural health monitoring with special emphasis on composite materials[END_REF]). In the current formulation it is not possible to consider changes in damping as constraints in the optimization, since no direct response for damping is available in Nastran so far.

Considering the used optimization algorithm, the implemented BigDOT optimizer in Nastran is an up-to-date algorithm with high efficiency and reliability. Since it is based on mathematical programming, it shows however common problems like possible wrong local optima and is certainly subjected to further developments in numerical optimization. Therefore, for a more detailed validation of the here presented tool, more studies on the optimization parameters and other damage cases are necessary. Still, BigDOT is one of the most capable optimization solutions readily available at the moment and it is especially capable of solving topology optimization problems [START_REF] Vanderplaats | Very large scale optimization[END_REF]. Further studies using different optimisation algorithms (stochastic or evolutionary methods, etc.) are encouraged, the principal problem in this damage localisation approach is however still the inclusion of further design responses that are sensitive to damages in the fibre composites.

Conclusions

In this paper an enhanced Topology Optimization technique for damage localization is developed based on the works of Lee et al. [START_REF] Lee | Damage detection by the topology design formulation using modal parameters[END_REF]. In the first part of this paper, the theoretical aspects of this technique are explained, and in the second part, validations of this technique are carried out on numerical models and also composite beams.

The presented damage localization approach can also be divided into two sub-problems, the optimization routine using Topology Optimization and the experimental estimation of modal parameters with the correct numerical modeling.

Concerning the results of the numerical validation, it can be concluded that the feasibility of applying the Topology Optimization approach for damage localization is proven also for an implementation in Nastran. This provides the opportunity for a broader application of the method, especially to more complex structures, that can be more efficiently analyzed in standard industrial finite element codes than in self written codes.

In the process of experimental validation, we succeed in localizing seperate damage zones, but the results become less clear and the optimization might get stuck in wrong local optima regarding the results of C-Scans of the specimens. But it can be assumed, that the poor results for the D2 damage case follow from the simplistic modeling, the lack of constraint data and not taking into account damping changes, rather than due to a conceptional error in the approach. Still, with a suitable choice of parameters, some good results with correct localization of the impacts on the composite beams can be archieved.

However, the approach has yet only been validated for one composite laminate beam. In the near future, this study shall be extended to all the other composite beams tested in reference [START_REF] Shahdin | Significance of low energy impact damage on modal parameters of composite beams by design of experiments[END_REF]. The approach shall also be applied to localizing damage in standard sandwich specimens with honeycomb and foam cores, as C-Scan testing is not very capable of detecting damage in composite sandwiches. Further prospects of this work could be an optimization using more design 20 responses, like directly addressing the mentioned anti-resonances or mode shapes by referencing the imaginary part of the FRFs, and to increase the analysed frequency bandwidth to take into account higher modes.

$-- 1 20 0 $---------------Geometry------------------------------$--Include Model file here (*.bdf) INCLUDE 'modelgeometry.bdf' $---------------BCs-----------------------------------$--referencing BCs in geometry file SPCADD 2 1 $---------------LOADs---------------------------------$--referencing Loads in geometry file RLOAD1 .001 DESMAX 50 CONVDV 1e-6 CONVPR 1e-6 $--------------End------------------------------------ENDDATA

  (a) CFRP laminate specimen (b) Experimental set-up for vibration testing

Figure 2 :

 2 Figure 2: Specimens for vibration testing (a) and schematics of the used Oberst beam test set-up (b).Masses were glued to the ends of the specimens to enhance the change in modal parameters due to damage[START_REF] Vanhoenacker | The use of multisine excitations to characterise damage in structures[END_REF]

Figure 3 :

 3 Figure3: Location of impacts and measurement points on composite beams[START_REF] Shahdin | Correlating low energy impact damage with changes in modal parameters: A preliminary study on composite beams[END_REF]

Figure 4 :

 4 Figure 4: Frequency response functions for point 2 (cp. fig. 3) for a beam impacted with 8J obtained by vibration testing. This figure qualitatively highlights the shifts in resonance frequencies from tab. 2

Figure 5 :

 5 Figure 5: Resonances and anti-resonances in frequency response functions

Figure 6 :

 6 Figure 6: Objective functions for fractional mass penalization (see. tab. 3. Each function posses a relative minimum at m f rac = 1, that helps to keep the optimization close to the original model

Figure 7 :

 7 Figure 7: Construction of objective functions and global or subcase design constraints in Nastran input files.

  (a) Finite element model of a cantilever beam with applied forces and constraints that has been used as the initial design for preliminary tests of the damage localization method. (b) Finite element model for analysed composite beams. The point used for comparison of frequency response functions is marked.

Figure 8

 8 Figure 8

Figure 9 :

 9 Figure 9: Experimental and numerical estimated frequency response functions of the composite beams for the point marked in fig. 8b.

Figure 10 :

 10 Figure 10: Results of Topology Optimization on a two-dimensional isotropic model: (a) Resulting density distribution; (b) target model with a cutout.

Figure 11 :

 11 Figure 11: Results of Topology Optimization on a two-dimensional isotropic model after processing with the imhmax-function provided in Matlab [35]. This function supresses false damage regions and leads a better estimation for damaged regions (fig. 11b)

Figure 12 :

 12 Figure 12: Frequency response functions at the excitation node for the damaged (target) model and the optimization model before and after optimization

Figure 13 :

 13 Figure 13: Comparison between C-SCAN results (top) and Damage Localization by TO (bottom) for the damage state D1 -The two middle points are blackened as they correspond to the damage state D2. The applied algorithm is able to easily detect this kind of damage (2 isolated impacts) with high accuracy.

Figure 14 :

 14 Figure 14: Comparison between C-SCAN results (top) and Damage Localization by TO (bottom)for damage state D2. The developed method highlights here a real limitation: Damage can be detected but it cannot be distinguished between each of the damage sites. The applied algorithm is only capable to show a damage zone delimiting the 4 impacts.

Figure 15 :

 15 Figure 15: Experimentally obtained FRF for the composite beams and calculated FRF of the FEM model at the same point (cp. fig. 8b), before and after optimization using only resonance frequency constraints : (a) for 4 impacts; (b) for 8 impacts.

  -------------Further Parameters--------------------$--Topoly Optimisation Variable $--(Choice of correct design domain!)

Table 1 :

 1 Properties of pre-impregnated carbon/epoxy T300/914 plies

	Property	Symbol	Value
	Young's modulus in fibre direction	E 1	122000MPa
	Young's modulus in transverse direction E 2	8500MPa
	Shear modulus	G 12	3570MPa
	Poisson ratios	ν 12 ; ν 23 ; ν 31 0.25; 0.3; 0.017
	Fibre volume ratio		60%
	Density	ρ	1550kg/m 3
	Number of plies		24
	Stacking sequence		[(0/90/45/ -45) 3 ] s

Table 2 :

 2 Estimated resonance frequencies for a composite beam before and after impacting (Impact energy: 8J / PolyMAX method)[START_REF] Shahdin | Correlating low energy impact damage with changes in modal parameters: A preliminary study on composite beams[END_REF]

	Frequencies [Hz]

Table 3 :

 3 Objective functions for fractional mass penalization

	Number	Function	Derivative
	(1)

Table 4 :

 4 Definition of Topoloy Optimization variable in Nastran

	Command ID	Label	Type of property	χ 0	χ min	∆χ p Property ID
	TOPVAR	1 DREGION	PSHELL	0.99 0.0001 0.2 3	4

Table 5 :

 5 Natural frequencies of the 2D isotropic test beam with a cutout and the model with optimized density distribution

		Natural Frequencies / Hz
	Mode Target model Optimized model
	1	16	16
	2	94	94
	3	238	238
	4	256	254
	5	257	256
	6	513	516

Table 6 :

 6

			Frequencies / Hz
	Mode		4 Impacts	8 Impacts
		Target Optimization Target Optimization
	1 (1st bending)	37	30	36	32
	2 (2nd bending)	157	148	155	156
	3 (3rd bending)	687	679	678	679
	4 (4th bending)	1363	1351	1326	1326
	(a)				(b)

Resonance frequencies of tested CFRP beams after 4 and 8 impacts and of the corresponding model with optimized density distribution. Especially the higher frequencies show a relatively good match.

Table for

 for

		modal damping ratios
	TABDMP1 1	CRIT
	0.	0.003	1600.	0.003	ENDT
	$--Eigenvalue extraction table
	EIGRL

Acknowledgements

The authors gladly thank PhD student Thibault Gouache for his contributions to this article and to the structural dynamics work group of the ISAE Toulouse.

Appendix

A. Nastran Input File for Damage Localization by Topology Optimization

This Nastran input file can generally be applied to finite element models when the correct model geometry and optimization constraints are included and the commented parameters are adjusted. ---------------------- --------------Subcases----------------------------- --------------Bulk Data Section---------------------BEGIN BULK $---------------Parameters---------------------------

$--------File Managementsection


Y Hichri
email: yosra.hichri@ifsttar.fr

V Cerezo

M T Do

H Zahouani

Effect of particles' characteristics and road surface's texture on the tire/road friction

This paper presents a laboratory study to investigate the lubrication induced by particles deposited on road surfaces which is responsible for the increase of accidents at the first rain after a long dry period. Particles are extracted from sediments collected at a catchment area and characterized by their chemical composition and size distribution. Protocol to simulate the particle deposit on the road surface and their compaction by the traffic is described. The test program includes variables such as particle concentrations, particles' size fractions and surface textures. Dry friction tests are conducted using the Skid Resistance Tester which simulates the sliding friction between a rubber slider and the test surface. Friction is lowered when the surface is covered by particles, compared with a clean surface. Repeated passages of the slider induce an increase of friction coefficient until reaching a stable value. Particles' flows during a friction run are separated into those ejected from the sample, those stored by the surface macrotexture and those that stay on the test surface; this study demonstrated that the friction coefficient is closely related to the quantity of particles available on the test surface. SEM analysis shows on the other hand that these particles are mainly trapped by the surface microtexture. Similar behavior between fine particles and powder in terms of lubrication mechanisms is observed.

Introduction

Previous studies showed that the occurrence of road accidents is statistically high at the first rain after a long dry period. [START_REF] Eisenberg | The mixed effects of precipitation on traffic crashes Accident Analys is and[END_REF] found a significantly increasing crash risk with the number of dry days since the last rain. This excess risk is generally attributed to the accumulation of fine particles on the road surface. This contamination is due to both non-exhaust emission of traffic (tires, brakes, etc.) and vegetation debris. [START_REF] Egodawatta | Translation of small plot scale pollutant build-up and wash-off measurements to urban catchment scale[END_REF] observed that particle's build-up is a function of the number of dry days and other weather conditions such as wind characteristics, traffic, etc. [START_REF] Shakely | Effects of pavement contaminants on skid resistance[END_REF] found on the other hand that the tire/road friction decreases with increasing dust concentration.

Despite the above observations, the link between particles deposited on road surfaces and the tire/road friction is not well established, mainly due to a lack of understanding of involved mechanisms. At a microscopic scale it has been highlighted, in the case of metallic wheel/rail contact [START_REF] Descartes | Experimental identification and characterization of the effects of contaminants in the wheel-rail contact[END_REF], that contaminants affect surface properties, adhesion or electrical insulation, and possibly damage or protect the contacting surfaces; thus they have a significant effect on the contact behavior. [START_REF] Heshmat | Wear reduction systems for coal-fueled diesel engines -Experimental results and hydrodynamic model of powder lubrication[END_REF] has demonstrated that solid particles can form load-carrying lubricant films which exhibit a behavior resembling that of fluid films. [START_REF] Mills | The mechanisms of pedestrian slip on flooring contaminated with solid particles[END_REF] studied friction on dry floors covered by particulate contaminants and stated that there is a critical size of the particles (around 50-60 µm) above which lubrication between a shoe's sole and the floor is due to particles rolling one on the other and below which shearing of the particles' layer (formed from particles' cohesive forces) is predominant. The floor surface roughness can trap particles and contributes to the shearing of the particle layer [START_REF] Mills | The mechanisms of pedestrian slip on flooring contaminated with solid particles[END_REF].

In a previous study (Hichri et al 2017a), a test protocol to assess dry friction on road surfaces contaminated by particles was developed. These authors showed that particles lubricate the tire/road interface and confirmed the particles trappingas a lubrication mechanismstated by [START_REF] Mills | The mechanisms of pedestrian slip on flooring contaminated with solid particles[END_REF]. In this paper, efforts have been continued to assess the effect of influential factors such as the road surface texture and the particles' size and concentration. From friction tests performed on surfaces with different textures and covered by particles of different concentrations and sizes, the particles' flows at the tire/road interface are deduced using the third body concept.

Experiments

Test samples

Two aluminum rectangular slabs (surface of 130mm  80 mm ≈ 0.01m 2 ) are used. A surface roughness similar to that of a road surface is obtained by sandblasting the surface with corundum particles (diameter superior to 2 mm) to create the microtexture (surface irregularities whose dimensions are below 0.5 mm horizontally and vertically) and grooving to create macrotexture (surface irregularities whose dimensions range between 0.1 mm and 20 mm vertically, and between 0.5 mm and 50 mm horizontally). One slab is sandblasted (Fig. 1a), named "Sandblasted" in the paper, and the other is sandblasted and grooved (Fig. 1b), named "Chocolate". Sq for Sandblasted and Chocolate is respectively 13.6µm and 13.9µm (mean of Sq values calculated from three areas of 5mm  5mm). Grooves have a triangular shape (2 mm in width at its basis and 1 mm in depth). It creates a grid with regular squares of 10 mm aside.

Particles collection and characterization

Particles are collected at the catchment area of Cheviré bridge with 90,000 vehicles daily (8.5% of trucks) (Figure 2). Runoff water from the South part of the bridge is discharged to this area. Collecting particles using a vacuum machine, as practiced in previous studies [START_REF] Vaze | Experimental study of pollutant accumulation on an urban road surface Urban[END_REF]Chiew 2002, Wilson 2006), would be an appropriate option but it would provide too small quantity of particles for laboratory tests. (Hichri et al., 2017b) The sediments extracted from the catchment area are dried at 40°C during four days and sieved to keep only particles smaller than 100 µm that affect surface properties and friction [START_REF] Descartes | Experimental identification and characterization of the effects of contaminants in the wheel-rail contact[END_REF][START_REF] Arias-Cuevas | Rolling-sliding laboratory tests of friction modifiers in dry and wet wheel-rail contacts[END_REF]. SEM (Scanning Electron Microscopy) observations of particles show that they exhibit various shapes and sizes. Moreover, particles mainly contain silica and metallic elements that can be due to the wear of braking system and tires (Fig. 3). Na, K and Mg elements could come from industrial sites located in the neighborhood of the bridge. (Hichri et al., 2017b).

Particles are spread manually on the slabs (Fig. 4a) and then compacted (simulation of vehicular traffic) by means of a wheel mounted on a metallic arm (Fig. 4b). Compaction is performed through 60 repeated passes of the wheel. This protocol allows obtaining a homogenous distribution of the particles. Figure 4c shows a slab covered by particles.

Friction tests

Friction tests are performed by means of a Skid Resistance Tester Pendulum (Fig. 5), which is well-known in civil engineering field to assess skid resistance in laboratory and on real road. The Pendulum is composed of an articulated arm equipped at its end with a rubber pad, which slides on the test surface. When the arm is released from its initial horizontal position, it moves up to a height where values on a graduated scale can be read and correspond to a friction coefficient. The test is standardized (CEN 2009).

Figure 5. Skid Resistance Tester Pendulum

Test protocol

Graduated scale

Rubber slider

Articulated arm (horizontal position)

The test protocol follows the following steps:

 On dry and clean surfaces (step numbered 0): measurement of 3D height cartographies and friction coefficient; weighing of the sample (mass m0).  Particles arrangement (see 2.3).  On surfaces covered by particles (step numbered i, i  1): friction measurement; weighing of the sample before/after friction measurement.  When the friction coefficients recorded at three successive steps do not differ more than ± 0.01, stop the series.

At each step, only one run is performed for friction measurement. This is due to the fact that the distribution of particles evolves at each friction run. Nevertheless, to assess the repeatability of the test method, some test series were repeated and it was found that values of the friction coefficient at every step are very repeatable (coefficient of variationratio standard deviation/meanless than 4% at each step) (Hichri et al 2017a).

A new rubber pad is mounted on the Skid Resistance Tester for each test series. The rubber is carefully wiped after each friction measurement to limit the contamination of the pad. A high resolution optical sensor, called Infinite Focus, based on the focus variation principle is used to capture 3D images from which surface texture profiles or cartographies are extracted. 3D cartographies are measured on areas of 5mm x 5mm in dimensions, located at the center and at the two extreme ends of the specimen (where the friction slider respectively touches and leaves the test surface). A lens enabling 20 times magnification is used for the capture. The vertical resolution of the sensor is 50 nm. Analyses are performed using the Mountains Map® software. Non-measured points are first replaced by values of neighbor points. A lowpass filter is then used to remove the form of the surface. Finally, the cartography is levelled using a leastsquare plane. Texture parameters are then calculated on the extracted profiles. SEM observations are made in the same areas.

Calculation of the mass of particles

The mass loss of the rubber pad during a series of friction measurements is negligible (few milligrams) compared with the initial mass of particles (hundreds of milligrams). Moreover, the clean sample exhibits the same weight at the beginning and end of each test series. Thus, the amount of detached particles from the two contacting bodies (slider and sample) is negligible. It is assumed in the rest of the paper that the mass of the particles matches the difference in weight between contaminated and clean states.

The mass of particles at step (i) can be then calculated as follows:

m particles,i = m specimen,i -m specimen,0 (1)
where m particles,i is the mass of particles at step (i); m specimen,i is the mass of the specimen at step (i); m specimen,0 is the mass of the clean specimen.

Results

Variations of friction and mass of particles

Variations of friction coefficient on "Sandblasted" and "Chocolate" surfaces are provided in figure 6. It can be observed that the friction coefficient on contaminated surfaces (steps i; i ≥ 1) is always lower than that of a clean surface (step 0). The friction coefficient increases with the number of friction runs and reaches a stable value. Figure 6. Friction coefficient evolution on "Chocolate" and "Sandblasted" contaminated surfaces

The final value of friction coefficient is close for the two surfaces (0.82 and 0.80 in average for respectively Sandblasted and Chocolate). However, the evolution is more rapid for Chocolate than for Sandblasted.

Figure 7 shows that the mass of particles decreases when the friction coefficient increases. According to [START_REF] Mills | The mechanisms of pedestrian slip on flooring contaminated with solid particles[END_REF], a shear effect due to the slider induces a reduction of the particles layer thickness until the surface asperities emerge and restore direct contact between the slider and the surface. A decrease of particles' mass increases then the slider/surface contact area and explains the friction increase.

Figure 7. Friction coefficient and particles mass evolution during a friction test (Sandblasted surface; 20g/m 2 ; fraction 0-40µm)

Effect of surface texture

3.2.1. Particles trapping by the surface microtexture SEM observations show that the surface coverage by particles decreases after successive passages of the friction slider. In addition, the particles that remain on the surface are those trapped by the surface valleys; it would be difficult to remove these particles, which then continue to lubricate the interface and induce lower friction compared with a clean surface.

Particles storage by the surface macrotexture

In figure 9, a sharp increase of friction on the Sandblasted surface can be observed for small quantity of particles (mass of particles around 0.03g). According to [START_REF] Heshmat | Wear reduction systems for coal-fueled diesel engines -Experimental results and hydrodynamic model of powder lubrication[END_REF] who studied the powder lubrication, the sharp increase of friction can be due to the depletion of particles until the point where direct contact is established between the two first bodies (considering powder as a third body).

Figure 9. Friction coefficient evolution as a function of mass of particles (20g/m 2 ; fraction 0-40µm)

A horizontal shift can also be seen between the two studied slabs, which means that different amount of particles generate the same friction coefficient! The explanation can be found in the presence of the grooves on the Chocolate slab, which store part of the weighed particles and, yet, do not contribute to the contact. Therefore, to be comparable, weighing results of Chocolate must be corrected to remove the mass of particles stored in the grooves. [START_REF] Fillot | Wear modeling and the third body concept[END_REF] demonstrated that the analysis of the wear-induced particles' mass can help to understand particles' flows inside the contact. Given the similarities between Fillot and coauthors' study and this study (movement of dry particles trapped between two rubbing surfaces), it is proposed to adopt the same approach to better interpret friction variations.

Mass analysis

Modeling and analysis of particles' flows

Two particles' flows can be defined in the contact area (Fig. 10): 1) Ejection flow, which is composed of particles ejected from the contact area and those falling into the grooves; 2) Circulation flow, which is composed of particles retained by the surface valleys.

For friction analysis, only the circulation flow is of interest. However, weighing includes this flow and particles stored in the grooves. Therefore, a model is developed to isolate the circulation flow in the weighing.

Figure 10. Particles' flows in the rubber pad/surface contact area (from Hichri et al., 2017c) In figure 11 right, particles layer on Chocolate slab at run (i-1) is illustrated (mass mi-1). At run i (Fig. 11 left), the extra quantity of particles falling into the grooves is supposed to represent a fraction k (0 < k < 1) of mi-1. It is assumed furthermore that k remains constant during a test series (for a given particles' concentration and size fraction).

m i-1 k.m i-1 run i-1 run i m i-1
For run 1, the amount of particles in the grooves is supposed to be: ). These quantities are then used to estimate the quantity of particles in the grooves at run 2 (m particles,1 groove ) (Eq. 3), from which the quantity of particles at the surface at run 2 (m particles,1 surface ) can be deduced (Eq. 2). The calculations at run 2 can be repeated for the following runs.

m particles,1 groove = k m particles,1 (4)

Estimation of mass of particles stored by the surface macrotexture

As the microtexture of Sandblasted and Chocolate slabs is similar (Sq are respectively 13.6µm and 13.9µmsee 2.1), it is assumed that the same mass of particles from the circulation flow provides the same friction coefficient. A criterion to determine k is then the superposition of the two sets of points in figure 8.

Results are illustrated in figure 12. Figure 12a shows results for uncorrected and corrected masses of particles for the Chocolate. Values of k mean that between 8% and 35% of particles at the surface fall into the grooves at each friction run.

Regarding the size of the grooves (44mm 2 ) and the tablet area (100mm 2 ), the ratio between them is 0.44. It means that, in maximum, 44% of particles can fall into the grooves. This calculation should be refined but it shows that the adjusted values of k are in the order of what could be expected.

Figure 13. Dimensions of the tablets and the grooves

In figure 14, if one considers the k value of combination (concentration 20 g/m 2 , fraction 80-100µm) as unusual, it could be said that k increases for increasing particles' size, meaning that more coarse particles are expelled by the friction slider from the slider/surface contact area. This trend can be explained by analysis of [START_REF] Mills | The mechanisms of pedestrian slip on flooring contaminated with solid particles[END_REF] following which, on the one hand, particles of less than 50µm in size tend to stick one to the other due to cohesive forces and, as a consequence, form a layer and, on the other hand, particles bigger than 50µm tend to roll one on the other. The consequence of these two mechanisms is that it would be more difficult to expel fine particles (< 50µm) than coarse ones (> 50µm), explaining the resulted values of k.

Relationship between friction and mass of particles

3.4.1. Effect of particles' size Figures 15a and15c show variations of the friction coefficient with the number of friction runs at different particles' sizes for the two studied surfaces. For both surfaces, it can be seen that at a given step the friction coefficient increases with increasing particles' size. For fraction 80-100µm, values of friction coefficient at stabilization are close to value at the beginning (clean state). For fraction 0-40µm, values of friction coefficient at stabilization (around 0.8) are well below value at the beginning (around 1.05). These observations are a consequence of the particles' flows analyzed in the previous section. Actually, coarse particles are easily expelled by the rubber slider and the surface is cleaned more rapidly. Inversely, fine particles are trapped by the surface microtexture or form a layer and it is thought that the surface contamination will last longer. In figures 15b and 15d, friction coefficient is plotted against mass of particles (same friction data as in figures 14a and 14c respectively). It can be observed that the effect of particles' size is no more visible using this presentation. Based on analyses made by [START_REF] Mills | The mechanisms of pedestrian slip on flooring contaminated with solid particles[END_REF] on the size-dependent lubrication induced by particles (rolling for coarse particles and shearing for fine ones), one could expect some differences between the four fractions. Further understanding of particles' interactions would help to explain the obtaining of this master curve.

Effect of particles' concentration

In figures 16a and 16c, variations of the friction coefficient with the number of friction runs at different particles' concentrations for the two studied surfaces are shown. For fraction 0-40µm, it can be seen that the friction recovering is longer for concentration 40g/m 2 compared with the two other concentrations. For fraction 80-100µm, this trend is confirmed but values of friction coefficient at stabilization are higher for concentration In figure 16b and16d, friction coefficient is plotted against mass of particles (same friction data as in figures 16a and 16c respectively). It can be observed that the effect of particles' concentration is no more visible using this presentation.

Conclusions

The lubrication of tire/road interface induced by particles deposited on the road surface was investigated in laboratory. Particles of different sizes (fractions of 0-40µm, 40-50µm, 50-80µm and 80-100µm) were spread at different concentrations (10g/m 2 , 20g/m 2 and 40g/m 2 ) on two surfaces having the same microtexture (Sq of 13.6µm and 13.9µm respectively) with and without macrotexture (simulated by grooves) respectively. The test protocol helps to evaluate the variation of friction coefficient, measured by a so-called Skid resistance tester Pendulum, as well as mass of particles with successive passes of the Pendulum's rubber slider simulating sliding of tire tread blocks on a road surface.

The effect of surface texture has been observed through two mechanisms: a trapping of particles by the roughness of surface microtexture and a storage of particles by the "reservoirs" created by surface macrotexture. The first mechanism induces continuing contamination of the road surface and a lower friction coefficient compared with a clean surface. The effect of the second mechanism is not clearly understood but it can be thought that stored particles can be entrained by the tire to the road surface and contaminate it (this mechanism can be studied in the future).

Particles smaller than 40µm can be trapped by surface microtexture or form a layer; by this fact, they are not easily expelled from the contact. Fractions higher than 40µm in size are less affected by the trapping mechanism and more easily expelled from the contact (ejected from the sample or stored by surface macrotexture). The consequence of this size effect is that friction recovering (the way friction coefficient of contaminated surfaces tends toward that of a clean surface) is slower when particles are smaller than 40µm. The effect of particles' concentration is less visible even if a slower friction recovering was observed for the highest concentration (40 g/m 2 ).

Particles' flows inside the contact were modeled and a method was proposed to correct particles' weighing to analyze only the mass of particles present on the test surfaces. Within the frame of conducted experiments, it was observed that the friction coefficient decreases when the mass of particles increases. The striking point is that when friction coefficient is plotted against mass of particles, a master curve is obtained, regardless the particles' size and concentration.

  Figure 1. "Sandblasted" (a) and "Chocolate" (b) slabs.

  Figure 2. Particles' collection site(Hichri et al., 2017b)

  Dried particles are sieved to obtain four size fractions: 0-40µm, 40-50µm, 50-80µm, 80-100µm. Values are chosen with respect to the critical size (around 40 to 50µm) separating "fine" (0-40µm and 40-50µm) and "coarse" (50-80µm and 80-100µm) particles' families (based on[START_REF] Mills | The mechanisms of pedestrian slip on flooring contaminated with solid particles[END_REF] analyses).

Figure 3 .

 3 Figure 3. Chemical composition for different size fractions of particles 2.3. Particles arrangement on the surface For each particle size fraction, three quantities of 0.1g, 0.2g and 0.4g of dry particles are laid on the test surfaces, which represent respectively 10g/m², 20g/m² and 40g/m² in concentrations. The 20g/m² concentration is representative of the amount of particles in urban areas (Vaze and Chiew 2002, Deletic and Orr 2005) and the two other represent limit values.

Figure 8 .

 8 Figure 8. SEM captions of a) the surface recovered by particles and b) the surface after 5 friction runs

Figure 11 .

 11 Figure11. Quantity of particles falling into the grooves (fromHichri et al., 2017c)

  It is therefore possible to estimate the quantity of particles in the grooves m particles,

Figure 12 .

 12 Figure12. Friction coefficient evolution as a function of corrected and uncorrected masses of particles (20g/m 2 ; fraction 0-40µm)

Figure 14 .

 14 Figure 14. Variation of k factor for different particles' concentrations

  Variation of friction and mass of particles for different particles' sizes

  with 10g/m 2 . No satisfactory argument is found to explain this difference. Variation of friction and mass of particles for different particles' concentrations

Table 1 .

 1 Adjusted k factor for different particles' concentrations and fractions.

	Fraction