Andrés Serna

Beatriz Marcotegui

Etienne Decencière

Thérèse Baldeweck

Ana-Maria Pena

Sébastien Brizion

Segmentation of elongated objects using attribute profiles and area stability: application to melanocyte segmentation

Keywords: mathematical morphology, segmentation, connected component, attribute profile, elongation, area stability, multiphoton microscopy

In this paper, a method to segment elongated objects is proposed. It is based on attribute profiles and area stability. Images are represented as component trees using a threshold decomposition. Then, some attributes are computed on each node of the tree. Finally, the attribute profile is analyzed to identify important events useful for segmentation tasks. In this work, a new attribute, combining geodesic elongation and area stability is defined. This methodology is successfully applied to the segmentation of cells in multiphoton fluorescence microscopy images of engineered skin. Quantitative results are provided, demonstrating the performance and robustness of the new attribute. A comparison with MSER is also given.

Introduction

Filtering techniques, aiming at removing noise while preserving as much as possible the desired information, are often essential prior to segmentation. Mathematical Morphology [START_REF] Matheron | Random Sets and Integral Geometry[END_REF][START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF] is a theory of non-linear operators based on a set approach in order to study the objects morphology. Classically, it uses structuring elements. It has been shown that adaptive approaches can lead to important improvements [LDM07, MV09, PD09, [START_REF] Jos | Adaptivity and group invariance in mathematical morphology[END_REF][START_REF] Angulo | Morphological Bilateral Filtering and Spatially-Variant Adaptive Structuring Functions[END_REF].

In this paper, we propose a method to segment elongated cells (melanocytes) in multiphoton fluorescence microscopy images of engineered skin. Segmenting these images with standard methods fails since melanocytes are 1 low contrasted and noisy. Moreover, it is proven that including shape prior knowledge improves the segmentation results.

The proposed method is based on the analysis of the attribute profile over the threshold decomposition of an image. We define a new attribute, called area-stable elongation, that combines the geodesic elongation and the area stability. In our experiments, we analyze important events in the evolution of this attribute and we show its efficiency in segmenting elongated objects while filtering out noisy structures.

This paper is organized as follows. Section 2 reviews related work in the state of the art. Section 3 presents the threshold decomposition, the attribute profile and introduces the new area-stable elongation attribute. Section 4 illustrates an application segmenting elongated cells in fluorescence multiphoton microscopy images of engineered skin. Finally, Section 5 concludes the paper.

2 Related work [START_REF] Jones | Connected Filtering and Segmentation Using Component Trees[END_REF] proposed connected filters using attributes signatures, i.e. the evolution of an attribute on the component tree. He has successfully applied his method to the segmentation of wood micrographs. [START_REF] Pesaresi | A new approach for the morphological segmentation of high-resolution satellite imagery[END_REF] introduced morphological profiles using the derivative of the residues from opening/closing by reconstruction. Their method is well suited for images with low contrast and low resolution. However, the maximal residue may not be the best segmentation choice. Moreover, the computational cost increases when processing large and homogeneous images. [START_REF] Beucher | Numerical residues[END_REF] proposed the analysis of the residue evolution through successive morphological operations. This evolution over each pixel of the image leads to interesting transformations such as ultimate openings and quasi-distance functions. More recently, [START_REF] Ouzounis | Differential Area Profiles: Decomposition Properties and Efficient Computation[END_REF] proposed differential area profiles for efficient point-based multiscale feature extraction in pattern analysis and image segmentation.

Maximally Stable Extremal Regions (MSER), proposed by [START_REF] Matas | Robust widebaseline stereo from maximally stable extremal regions[END_REF], is a well-known region detector. MSER are invariant to affine transformations of both intensity and image coordinates. They have a high repeatability and can be run in linear time with respect to the number of pixels in the image (5 Mpix/s using a max-tree implementation [START_REF] Nistér | Linear Time Maximally Stable Extremal Regions[END_REF]). However, the parameter selection remains its major drawback. Even when default parameters perform well in many applications, some heuristics need to be applied in order to yield appropriate regions. Moreover, MSER favors round regions, as proved by [START_REF] Kimmel | Are MSER Features Really Interesting[END_REF], making it unsuitable to detect irregular shapes such as elongated objects.

[FL07] computed SIFT descriptors on the MSER shape in order to characterize each region of the image. This approach is proven to be robust to illumination changes and nearby occlusions. They also proposed a pyrami-dal searching to reach scale invariability. The authors also suggested the use of MSER for image segmentation. [START_REF] Forssen | Maximally Stable Colour Regions for Recognition and Matching[END_REF] extended the MSER concept to color images and [START_REF] Litman | Stable volumetric features in deformable shapes[END_REF] defined stable volumetric features in deformable shapes.

This paper is an extended version of [START_REF] Serna | Attribute Controlled Reconstruction and Adaptive Mathematical Morphology[END_REF] work. In that framework, the attribute profile is built over increasing quasi-flat zones and three application are presented: image segmentation, adaptive mathematical morphology and feature extraction. In this paper, we extend the idea of analyzing the attribute profile over the threshold decomposition of an image. Additionally, we define an attribute profile based on a new attribute, called area-stable elongation, to segment elongated objects.

Methodology

Threshold decomposition and attribute profile

A binary attribute can be extended to gray-scale through threshold decomposition [START_REF] Wendt | Stack filters[END_REF][START_REF] Maragos | Threshold superposition in morphological image analysis systems[END_REF]. Let I be a digital gray-scale image I: D → V , with D ⊂ Z 2 the image domain and V = [0, ..., R] the set of gray levels. A decomposition of I can be obtained considering successive thresholds T t (I) = {p ∈ D|I(p) > t} for t = [0, ..., R -1]. Since this decomposition satisfies the inclusion property T t (I) ⊆ T t-1 (I), ∀t ∈ [1, ..., R -1], it is possible to build a tree, called the component tree, with these level sets T t (I). Each branch of the tree represents the evolution of a single connected component (CC) X t . An attribute profile is the evolution of an attribute on a branch of the tree.

Fig. 1 illustrates the threshold decomposition in the 1D case, its component tree and attribute (width) profiles for the two maxima of the function (p A and p B ). Events on this attribute profile are useful to segment objects [START_REF] Jones | Connected Filtering and Segmentation Using Component Trees[END_REF], extract features [START_REF] Pesaresi | A new approach for the morphological segmentation of high-resolution satellite imagery[END_REF] and define adaptive structuring elements [START_REF] Serna | Attribute Controlled Reconstruction and Adaptive Mathematical Morphology[END_REF].

In general, there are two types of attributes: increasing and non-increasing [START_REF] Breen | Attribute Openings, Thinnings, and Granulometries[END_REF]. On the one hand, an attribute is increasing when its value is greater or equal to the attribute computed on any subset of the object. The most common increasing attribute is the area, which is used to compute area openings and area stability, among others. On the other hand, an attribute is non-increasing when the latter property is not reached. In general, shape attributes -such as circularity, tortuousity, elongation-are not increasing.

Geodesic Elongation

In this work, we focus on geodesic elongation, simply called henceforth elongation [START_REF] Lantuéjoul | Geodesic methods in quantitative image analysis[END_REF]. The elongation E(X t ) of an object X t is a shape descriptor useful to characterize long and thin structures. It is proportional to the ratio between the square geodesic diameter L 2 (X t ) and the area of the object A(X t ), as shown in Eq. 1. The geodesic diameter L(X t ) = sup x∈Xt {l x (X t )} is the length of the longest geodesic arc of X t , i.e. the longest internal segment l x (X t ) connecting the two end points of X t [START_REF] Lantuéjoul | On the use of the geodesic metric in image analysis[END_REF].

E(X t ) = π 4 L 2 (X t ) A(X t ) (1)
The longer and narrower the object, the higher the elongation. The lowest bound is reached with the disk, where E(X t ) = 1. An example of elongation for binary objects is presented in Fig. 2. The number on each object corresponds to its approximated elongation. An efficient implementation can be found in [START_REF] Morard | Efficient geodesic attribute thinnings based on the barycentric diameter[END_REF].

A feature image can be defined assigning to each pixel p in the maximum attribute recorded over all CC X t (p in ) containing this pixel, i.e. the maximum of the attribute profile. Fig. 4(b) exemplifies the maximal elongation image max{E(X t )} on a synthetic example while Fig. 3 presents two real cases. It has to be noted that this feature image is a partition in which objects of an elongated shape are highlighted.

Maximally Stable Extremal Regions (MSER)

A simple but interesting attribute is the area A(X t ). When used to suppress small CCs, it leads to the definition of the area opening [START_REF] Vincent | Morphological area openings and closings for greyscale images[END_REF]. Since A(X t ) is increasing, events in the area profile are analyzed instead of its global maximum. For example, great changes in area are probably related to the union of different objects while small ones are related to area stable regions, as those detected by the MSER method [START_REF] Matas | Robust widebaseline stereo from maximally stable extremal regions[END_REF].

The area stability Ψ of the region X t is defined as the ratio between its area A(X t ) and its area variation dA (Xt) dt , as shown in Eq. 2:

Ψ(X t ) = A(X t ) dA(Xt) dt (2)
A MSER is a CC X t with maximal area stability. In the original proposition, every local maximum is detected, thus making possible to have nested regions. Then, some heuristics can be applied to only select the most important peaks.

A new attribute: Area-stable elongation

Favoring regular (round) regions is one of the main limitations of MSER, as proved by [START_REF] Kimmel | Are MSER Features Really Interesting[END_REF]. Thus, it is not suitable to detect irregular shapes, as elongated objects. In order to detect elongated objects taking into account their area stability, we propose a new attribute Φ(X t ), called area-stable elongation. This attribute combines the area stability of a region Ψ(X t ) with its elongation E(X t ), as defined in Eq. 3:

Φ(X t ) = Ψ(X t )E(X t ) = A(X t ) dA(Xt) dt π 4 L 2 (X t ) A(X t ) = π 4 L 2 (X t ) dA(Xt) dt (3)
The area-stable elongation function Φ(X t ) is affine-invariant since Ψ(X t ) and E(X t ) are preserved under affine transformation of intensity (domain V ) and image coordinates (domain D), as stated by [START_REF] Matas | Robust widebaseline stereo from maximally stable extremal regions[END_REF] and [START_REF] Forssen | Shape Descriptors for Maximally Stable Extremal Regions[END_REF]. However, area-stable elongated regions are not invariant to blur. If blurring invariance is required, e.g. for a matching application, two possible solutions consist in weighting the stability function by the gradient magnitude along its boundary or pre-processing the image with a deblurring filter, as proposed by [START_REF] Kimmel | Are MSER Features Really Interesting[END_REF].

The maxima of Φ(X t ) represent area stable regions with a significant elongation. We propose to build a feature image using the maximal Φ(X t ), which implies: i) the feature image is a partition of the space, useful for segmentation; ii) each pixel contains information about the shape and the area stability of its neighborhood, which can be exploitable using prior knowledge.

Let us explain this new attribute with a toy example. Consider the 9×9 image of Fig. 4. For this example, we have approximated the euclidean distance on the 8-connectivity grid, i.e. the geodesic diameter of a pixel is equal to 1, the distance between horizontal and vertical neighbors is equal to 1, and the distance between diagonal neighbors is equal to √ 2. This toy image contains 4 gray-levels enumerated from t 0 to t 3 , and 6 CCs enumerated from A to F. Fig. 4(e) presents the component tree, where A(X), E(X), Ψ(X) and Φ(X) are the area, the elongation, the area stability and the area-stable elongation of a given component X, respectively. Note that the stability for the background (object A) is not defined since it is the root of the component tree.

The component tree contains two branches. Let us suppose that we aim at segmenting object C from the left branch, and object E+F from the right one. Let us analyze each case separately.

First, object C is an elongated object nested on a spurious elongated structure B. Analyzing the elongation profile, we can see that object B (E(X B )=6.10) is more elongated than object C (E(X C )=4.60), as shown in the maximal elongation image of Fig. 4(b). However, the stability of region C (Ψ(X C )=2.43) is higher than that of region B (Ψ(X B )=1.27), as shown in Fig. 4(c). Combining these two attributes, region C (Φ(X C )=11.18) has a higher area-stable elongation than region B (Φ(X B )=7.71), as shown in Fig. 4(d).

Second, object E is an elongated object that includes another elongated object F. Analyzing the elongation profile, we can see that object E+F (E(X E )=8.43) is more elongated than the single object F (E(X F )=3.14), as shown in the maximal elongation image of Fig. 4(b). Since their area stabilities are similar (Ψ(X E )=1.50 and Ψ(X F )=1.80, as shown in Fig. 4(c)), the highest area-stable elongation is obtained for the union of these two objects (Φ(X E )=12.65), as shown in Fig. 4(d).

It is noteworthy that applying a simple threshold (e.g. Φ(X t ) ≥ 8) in the maximal area-stable elongation image (Fig. 4(d)), the objects C and E+F are correctly segmented, which is not possible on the original image nor on the other two feature images (Figs. 4(b) and 4(c)).

Fig. 5 illustrates the behavior of our method on a real DNA image. The goal is to segment the elongated and bright fiber from the noisy background. Fig. 5(b) shows the maximal elongation image, where objects of an elongated shape are highlighted. However, spurious objects can be merged at low levels resulting in CC with high feature value, such as the porous structure in the center of the image. The maximal area stability (Fig. 5(c)) keeps also many noisy and non elongated structures in the background. Finally, Fig. 5(d) shows the area-stable elongation image, where most of noisy structures have been eliminated due to their low stability.

Application: cell segmentation

To illustrate the performance of our method, we apply it to segment elongated cells in multiphoton fluorescence microscopy images. Images correspond to reconstructed skin used in cosmetic research in applications such as screening of de-pigmenting and pro-pigmenting agents (Fig. 6(a)). This model contains two types of cells: keratinocytes and melanocytes. The latter are dendritic cells, more elongated and brighter than keratinocytes. An accurate segmentation of melanocytes becomes crucial in order to quantify the melanin in the skin. This value is used to assess the efficiency of the cosmetic ingredient. Our goal here is to segment melanocytes, which appear as bright elongated structures.

Segmenting these images with standard methods fails since melanocytes are low contrasted and noisy, as shown in Fig. 6(a). A first simple solution may consist in applying automatic thresholding, e.g. Otsu method [START_REF] Otsu | A Threshold Selection Method from Gray-Level Histograms[END_REF]. However, it does not work because foreground and background gray-distributions overlap, as shown in the histogram of Fig. 6(c). Thus, cells and background are not separable with a global threshold. In this application, we propose a segmentation method using the component tree in order to solve the problem of low contrasted cells. Besides, the use of the area-stable elongation introduces shape prior knowledge and offers robustness to noise. In such a case, each cell can be segmented if it appears in the component tree, even if its gray-level is much lower than that for other cells in the image. Moreover, thanks to prior knowledge about melanocyte shape, the result is improved, justifying the use of our proposed methodology.

In our experiments we have 8 manually annotated images of 511 × 511 pixels each. The spatial resolution is equal to 0.26 µm/pixel. The ground truth definition has been carried out by experts from L'Oréal Research and Innovation. Classical Precision (P ), Recall (R) and f mean = (2 × P × R)/(P + R) statistics are computed in order to evaluate our results. The recall (or completeness) is defined as the number of correctly segmented pixels divided by the number of pixels marked in the ground truth. The precision (or correctness) is defined as the number of correctly segmented pixels divided by the total number of segmented pixels.

To exemplify our method, let us analyze the attribute profile for a single pixel belonging to a melanocyte, called seed pixel and marked with a red x in Fig. 7(a). Fig. 7(b) presents the ground truth provided by an expert. Fig. 7(c) shows four attribute profiles: area A(X t ), elongation E(X t ), area stability Ψ(X t ) and area-stable elongation Φ(X t ). For visualization purposes, each attribute has been normalized dividing by its maximum value to be in the range [0, 1]. Additionally, the f mean is plotted in order to define the best possible segmentation for this cell. Figs. 7(d) to 7(i) show the evolution of the CC X t containing the seed pixel.

Segmentation methods using threshold decomposition are based on the hypothesis that objects of interest exist at some level of the tree. In our example, the best possible segmentation corresponds to X t=34 , for which the highest f mean is obtained. Other good segmentations are in the range X t∈ [34,30] . The whole melanocyte is not retrieved for X t>34 and it is merged with other structures for X t<30 .

Let us analyze each attribute profile, starting with A(X t ). Based on a prior knowledge about melanocytes size, attributes for t < 15 are not analyzed since they correspond to structures bigger than 75% of the whole image. Analyzing Ψ(X t ), its global maximum represents the most stable region X t=42 . This is an area-stable and round region but useless in such a case since it does not match the entire melanocyte. Another interesting attribute is E(X t ) since melanocytes are long and thin. Its global maximum corresponds to a CC merging three different objects X t=28 . From an areastability viewpoint, this region is not stable because it is generated merging three different objects in a small range t ∈ [30, 28]. Finally, the global maximum of the area-stable elongation Φ(X t ) appears at X t=34 , which is the best segmentation according to f mean .

Figs. 8 and 9 present two experimental results. Figs. 8(a) and 9(a) show the two input images with their corresponding manual annotations in Figs. 8(b) and 9(b). Figs. 8(c) and 9(c) present the max{E(X t )} images. Note that all melanocytes present a significant elongation, however some post processing is required in order to eliminate porous structures on the background. Most noisy regions are not area-stable, then the area-stable elongation Φ(X t ) appears suitable for the segmentation of this kind of objects, as shown in Figs. 8(d) and9(d). This example demonstrates the use of our area-stable elongation in order to enhance elongated objects with respect to a noisy background. Using this feature image, the melanocyte segmentation becomes an easy task.

A simple three-fold segmentation algorithm is used for this purpose: i) characterization: a feature image is computed using the maximal areastable elongation max{Φ(X t )}; ii) filtering small objects: in the feature image, small regions (smaller than 500 pixels) are eliminated using an area opening followed by an area closing. This parameter is not critical since the smallest cell in the database is approximatively 3000 pixels size; finally, iii) filtering objects with low attribute value: a simple threshold removes structures with low area-stable elongation. In our experiments we have used a threshold equal to 11 for all images. However, this parameter is not critical since several values produce similar results, as shown in the overall sensibility curve of Fig. 10. It is noteworthy that thresholds between 7 and 16 produce an overall f mean over 70%. Table 1 presents quantitative results and a comparison with respect to the classical MSER by [START_REF] Matas | Robust widebaseline stereo from maximally stable extremal regions[END_REF]. MSER regions have been computed using the algorithm directly provided by the authors [MTS + 05]. The MSER parameters have been consistently set with those used by our method, i.e. the minimal MSER area has been set to 500 pixels and the threshold decomposition has been carried out for all gray-levels (one by one) stopping when the object area is greater than 75% of the whole image. Other MSER parameters such as relative area and relative margins have been kept to their default values 0.010 and false, respectively. Setting up these two parameters is not intuitive and global improvements have not been obtained in our tests. Using our method, this kind of parameters is not required, which is an advantage with respect to MSER. MSER results may be also improved using some pre-processing step. However, it would require the setting up and the selection of the appropriate filter to do it. Another advantage of our method is that pre-processing is not used since the noise robustness is included in the area-stable elongation itself, as aforementioned in Section 3.4.

Figs. 11 and 12 present our experimental results showing the input image, the ground truth, the MSER result and our segmentation result. It is noteworthy that in most cases, melanocytes are correctly segmented by our method. Some problems are shown in Fig. 11(l) where a clearly nonelongated melanocyte in the upper left part of the image has not been segmented, and in Fig. 12(r) where a melanocyte in the lower left part has not been detected; given its low contrast, it is mixed with the background. Note that our method presents much better results than MSER for all images. As aforementioned, MSER favors round and regular regions. Thus, only a partial segmentation is possible. In fact, MSER corresponds in several cases to the cell nuclei.

Conclusions

We have presented a method to segment elongated objects using attribute profiles. Images are represented as component trees using threshold decomposition. Then, the attribute profile is analyzed and important events are recorded. In this work, the geodesic elongation and the area stability are combined to define a new attribute: the area-stable elongation. The behavior of this new attribute in relation to noise, blur and geometrical distortions is discussed. The global maximum of this attribute is computed for each pixel of the input image and a feature image is built. This image is a spatial partition where objects of interest can be easily extracted using a simple threshold.

This method can be interpreted as an extension of MSER favoring objects of a given shape. Another main difference with the classical MSER is that only the global maximum of the attribute profile is chosen, thus only the most stable and elongated region is kept. More generally, this work confirms the interest of attribute profiles for image analysis and segmentation. Our method has been successfully used in a cosmetic application aiming at segmenting melanocytes cells that appear as bright and elongated structures in multiphoton images. Standard methods fail because melanocytes are low contrasted and noisy. Better segmentations are obtained providing a prior knowledge about cells shape. One of the method limitations, common to all methods based on threshold decompositions, is that it can only segment CCs present in the component tree. The analysis based on attribute profiles can be extended to other attributes and to other hierarchical partitions. For example, [START_REF] Serna | Attribute Controlled Reconstruction and Adaptive Mathematical Morphology[END_REF] compute attribute profiles on increasing quasi-flat zones, which is useful when processing pre-filtered images. In the future, other interesting attributes such as porosity and tortuousity will be studied. Additionally, extensions to higher dimensional data (color and 3D) will be analyzed.

Figure 1 :

 1 Figure 1: 1D threshold decomposition (a), component tree (b), and attribute profile (c).

Figure 2 :

 2 Figure 2: Geodesic elongation for different binary objects. The elongation values have been approximated to get integer values.

Figure 3 :

 3 Figure 3: Maximal elongation images using threshold decomposition.

  (a) Toy image (b) Max. elongation (c) Max. area stability (d) Max. stable elongation (e) Component tree

Figure 4 :

 4 Figure 4: Toy example: maximal attributes images and component tree.

Figure 5 :

 5 Figure 5: Maximal elongation, maximal area stability and maximal areastable elongation on a real DNA image.

  (a) 511 × 511 pixels (resolution 0.26 µm/pixel). (b) Ground truth (melanocytes) (c) Histogram of melanocytes and background.

Figure 6 :

 6 Figure 6: Foreground and brackground gray distributions on a multiphoton image of engineered skin containing keratinocytes and melanocytes.

  Attribute profiles. Xt<15 is not considered because A(Xt) > 0.8 (d) Xt=42 (e) Xt=34 (f) Xt=31 (g) Xt=30 (h) Xt=28 (i) Xt=16

Figure 7 :

 7 Figure 7: Attribute profiles for pixel marked with an x in the input image (a). From (d) to (i) CC of X t containing x, for diffent values of t.

( a )Figure 8 :

 a8 Figure 8: Feature images using the global maximum in the attribute profile.

( a )Figure 9 :

 a9 Figure 9: Feature images using the global maximum in the attribute profile.

Figure 10 :

 10 Figure 10: Overall sensibility curves: threshold to eliminate objects with low area-stable elongation.

  (j) Image a (k) Image b (l) Image c

Figure 11 :

 11 Figure 11: Segmentation of melanocytes using area-stable elongation. First row: input image; second row: ground truth; third row: MSER; fourth row: our segmentation result.

Figure 12 :

 12 Figure 12: Segmentation of melanocytes using area-stable elongation. First row: input image; second row: ground truth; third row: MSER; fourth row: our segmentation result.

Table 1 :

 1 Comparison with respect to MSER. In each column, numbers on the left correspond to the proposed method, and numbers between parentheses to MSER.

	Image	Precision	Recall	f mean
	a	0.84 (0.83) 0.61 (0.36) 0.71 (0.51)
	b	0.74 (0.89) 0.84 (0.42) 0.78 (0.57)
	c	0.71 (0.52) 0.84 (0.38) 0.77 (0.44)
	d	0.78 (0.84) 0.78 (0.53) 0.78 (0.65)
	e	0.85 (0.73) 0.90 (0.40) 0.87 (0.52)
	f	0.84 (0.72) 0.68 (0.48) 0.75 (0.57)
	g	0.62 (0.52) 0.92 (0.26) 0.74 (0.34)
	h	0.83 (0.50) 0.86 (0.41) 0.84 (0.45)
	Overall: 0.78 (0.69) 0.80 (0.41) 0.78 (0.51)

Acknowledgements

Part of this work has been funded by TerraMobilita project.


Thibaut Chaise
email: thibaut.chaise@insa-lyon.fr

Jun Li

Daniel Nélias

Régis Kubler

Said Taheri b Gérard Douchet

Vincent Robin

Philippe Gilles

Modelling of multiple impacts for the prediction of distortions and residual stresses induced by ultrasonic shot peening (USP)

Keywords:

Semi-analytical method

During a manufacturing process, the ultrasonic shot peening (USP) technique can be used as the final surface treatment. The aim of this operation is to introduce surface compressive residual stresses in order to prevent crack propagation advancement. Although the numerical simulation method is able to predict the level of residual stresses in a peened part, the 3D modelling of the real USP process, in which many successive and shifted impacts take place, is very delicate to perform and costly in terms of computing time and memory space required. In this paper, a two step method based at first on the calculation of the averaged plastic strain tensor in a half-space by using a semi-analytical method and in a second time on the transfer of this plastic strain field to a finite element model is proposed in order to simulate the effects of the USP process in thin structures. The accuracy and advantages of the semianalytical method are validated by a benchmark with several finite element codes. Experiments, similar to the Almen test, are performed on thin plates of Inconel 600. Numerical results in terms of distortions and residual stresses are compared with the experimental data.

Introduction

The manufacturing process of a mechanical part is made by a succession of operations that usually involve severe mechanical and/or thermal loadings. In consequence, an element taken at the beginning of its service life usually exhibits residual stresses that can have a significant impact on its fatigue life. Surface or volume treatments are often applied at the end of the process to generate compressive residual stresses and reduce the harmful effects of the forming operations.

As an example, Inconel 600 (IN600), a nickel-based alloy, is widely used in construction of nuclear reactors due to its excellent mechanical properties at high temperature and corrosion resistance. The presence of nickel enhances the resistance against corrosion. However, when a component is placed in a corrosive environment such as the primary circuit of a nuclear reactor, it becomes sensitive to the phenomenon of stress-corrosion cracking. Therefore it must be ensured that the level of residual stresses after the manufacturing process is lower than the triggering threshold of stress-corrosion cracking. In order to reduce such risks, a final surface treatment can be applied on the part. The objective is to introduce compressive residual stresses to prevent occurrence of cracking.

Among the mechanical surface treatments, several have been studied over the years such as the classical shot peening, or more recently the ultrasonic shot peening (USP), laser peening or water jet peening. The USP process which is often applied in a local zone can be seen as a complementary technology to classical shot peening. Its main benefit compared to conventional shot peening is the low roughness of the surface due to a better quality of the shots and lower impact speeds. The USP has also the advantage that balls are contained and can be recovered after treatment. This allows in situ treatment of the surfaces when USP is used as a repairing technique.

Beyond the generation of residual stresses, peening processes can also be used to deform material, in which case it is called peen forming, see for example, [START_REF] Gariépy | Shot peening and peen forming finite element modelling-towards a quantitative method[END_REF]. For this second group of applications the prediction of the distortions generated by the process becomes critical. Applications of shot peening cover also the joining of material of dissimilar mechanical properties, see [START_REF] Harada | Joining and shaping fit of dissimilar materials by shot peening[END_REF].

Due to their various applications, peening processes are nowadays the subject of many studies. The simulation of these processes started more than a decade ago with the simulation of one single impact. [START_REF] Meguid | Finite element modelling of shot-peening residual stresses[END_REF] studied the effect of various properties, including shot radius, velocity and material hardening properties and analysed their influence on the depth of the compressive residual stress layer. They also noted the few interests at this time to the modelling of the peening processes. Rouhaud et al. (2005) have investigated the effect of temperature on the residual stress field due to one single impact, showing that it has a negligible effect for low velocity and relatively large shots. This effect was deeply studied in the work of [START_REF] Rouquette | Coupled thermo-mechanical simulations of shot impacts: effects of the temperature on the residual stress field due to shot-peening[END_REF] who showed that it reduces the level of residual stresses but only as a second order effect. It should be noted that the range of temperature was below 200 • C, which has a limited effect on the strain-stress curve for the material investigated (35NiCrMo16).

However, the results obtained for one impact cannot be considered significant of the multiple impact process, as was shown for example by [START_REF] Kang | Multiple impact modelling for shot peening and peen forming[END_REF]. They noted the importance of modelling an appropriate coverage with a sufficient number of impacts to better apprehend the residual stress field. [START_REF] Frija | Finite element modelling of shot peening process: prediction of the compressive residual stresses, the plastic deformations and the surface integrity[END_REF] used symmetries in the both surface directions to multiply the number of simulated impacts but the reference cell used in their simulations sees only one impact. They showed that the results obtained with this model show large differences with the experiments, showing the importance to account for multiple impacts.

Most recent studies thus concentrate on the prediction of residual stresses, and distortions with models involving multiple impacts. The prediction of distortions is important to determine the influence of the process and apply it to peen forming. [START_REF] Levers | Finite element analysis of shot peening[END_REF] started more than 10 years ago, predicting the distortions generated by the process on wing panels. For the validation of numerical models, the distortions are usually represented by the maximum deflection of thin peened plates that is called the Almen intensity or arc height. [START_REF] Guagliano | Relating Almen intensity to residual stresses induced by shot peening: a numerical approach[END_REF] proposed analytical formulae to link the Almen intensity to the level of residual stresses based on a finite element model involving 5 impacts. More recently, [START_REF] Miao | Experimental study of shot peening and stress peen forming[END_REF] studied experimentally the number of peening passes necessary to reach a saturation and provided empirical relations to link it to the arc height. A recent study by [START_REF] Gariépy | Shot peening and peen forming finite element modelling-towards a quantitative method[END_REF] even predicts the Almen distortion of specimen before the bolts maintaining it during the peening operation are removed.

Performing realistic simulations of multiple impacts is limited by three major difficulties. First the large number of impacts to be modelled leads to unaffordable computing time. For this reason, [START_REF] Wang | Compressive residual stress introduced by shot peening[END_REF] proposed empirical models to predict the residual stresses. [START_REF] Majzoobi | A three-dimensional simulation of shot peening process using multiple shot impacts[END_REF] found that the simulation of 25 impacts allows to obtain a uniform stress state but noted that this value can vary with the peening and material parameters. They also indicated that the optimum stress state can occur before the saturation point. As for many other authors, they used various symmetries in their multiple impacts model to reduce the computation time. [START_REF] Miao | On the potential applications of a 3D random finite element model on the modeling of shot peening[END_REF] proposed a review of these various symmetries and reference cell shapes used in the literature for multiple impact models. [START_REF] Klemenz | Application of the fem for the prediction of the surface layer characteristics after shot peening[END_REF] used similar symmetries to model 121 impacts and averaged the results over a representative area covered by several impacts, and finally confronted the numerical results to experiments. Those models suppose a deterministic distribution of the impacts that necessarily differs from the reality of the process. Nevertheless, an averaging of the simulation results over a sufficiently representative area leads to fair predictions of residual stresses. The influence of the impact distribution and especially the sequence of the impacts and the number of impact cycles was studied recently by [START_REF] Kim | A 3D FE model with plastic shot for evaluation of equi-biaxial peening residual stress due to multi-impacts[END_REF].

Several authors propose stochastic models of the impact distribution to account for the complexity of the process. For example, [START_REF] Bagherifard | A numerical model of severe shot peening (SSP) to predict the generation of a nanostructured surface layer of material[END_REF] simulated a very large number of impacts, up to 134, randomly distributed on a reference area and compared the results to residual stresses measured by X-ray diffraction to predict surface nanocrystallization. [START_REF] Mylonas | Numerical modelling of shot peening process and corresponding products: residual stress, surface roughness and cold work prediction[END_REF] also use a stochastic model of the shot distribution and consequently model a few tens of shots, randomly distributed on the impact surface with a realistic velocity distribution.

A third method proposed to account for multiple impacts is to determine the eigenstrain generated by one or a few number of impacts and to multiply them to determine the effect of large number of impacts. [START_REF] Achintha | Eigenstrain modelling of residual stresses generated by laser shock peening[END_REF] determined the eigenstrain generated by one pulse of laser shock peening, equivalent to one impact in conventional shot peening, and sum this result to represent the effect of the whole process. The resolution of the equilibrium of the impacted body containing these eigenstrains allows obtaining the final residual stresses. [START_REF] Song | An eigenstrain-based finite element model and the evolution of shot peening residual stresses during fatigue of gw103 magnesium alloy[END_REF] recently obtained satisfactory prediction of residual stresses by using a similar method. They used the plastic strain obtained at the centre of one single indentation to represent the effect of multiple impacts.

The second major difficulty is the determination of a realistic loading as input to the impact simulation. In the case of laser peening or hammering [START_REF] Hacini | Effect of impact energy on residual stresses induced by hammer peening of 304l plates[END_REF], the position and intensity of the impacts or shocks is deterministic, but the physics of the impacts remains the same. In most peening model, the loading is represented by two parameters: the impact velocity, usually taken constant for all impacts and the coverage rate, which is the ratio of the impacted surface over the initial surface. [START_REF] Majzoobi | A three-dimensional simulation of shot peening process using multiple shot impacts[END_REF] showed that the impact velocity has a significant influence on the results for multiple impacts. A correct determination of the impact velocity is thus necessary. [START_REF] Miao | On the potential applications of a 3D random finite element model on the modeling of shot peening[END_REF] performed simulations of several dozens of randomly distributed impacts at various velocities and averaged them over a representative cell to model a realistic distribution of the impacts and their velocities. The prediction of these velocities lies on a model of the process itself. [START_REF] Gariépy | Shot peening and peen forming finite element modelling-towards a quantitative method[END_REF] proposed a modelling of the air nozzle of a shot peening device to predict the impact velocity. For the case of USP, models of the peening chamber, accounting for the multiple impacts of the shots, are used. They rely on the definition of the coefficient of restitution that allows to predict the rebound velocity of a shot impacting a surface or another shot. Analytical formulae, as those proposed by [START_REF] Jackson | Predicting the coefficient of restitution of impacting elastic-perfectly plastic spheres[END_REF] allow to predict these coefficients of restitution. Recent work by Chaise et al. (2011) extended these formulae to materials with isotropic hardening to improve the prediction of the rebound velocities and therefore subsequent impact velocities within the ultrasonic chamber.

The last difficulty is the determination of the material properties of the impacted surface. Indeed at the impact velocities commonly involved in shot peening, viscous effects can occur. [START_REF] Klemenz | Application of the fem for the prediction of the surface layer characteristics after shot peening[END_REF] used for example an elasto-visco-plastic law with combined isotropic and kinematic hardening to predict the residual stresses. In the case of very large strain rates, as in laser peening for example, the determination of realistic parameters for the viscous part of the behaviour can become a challenge. [START_REF] Amarchinta | Simulation of residual stress induced by a laser peening process through inverse optimization of material models[END_REF] proposed an inverse method to determine the material behaviour at high strain rates. [START_REF] Seifried | Viscoplastic effects occuring in impact of aluminium and steel bodies and their influence on the coefficient of restitution[END_REF] showed that viscous effects can be neglected for certain materials while they are critical for others. In the present study a simple isotropic hardening law was used since it was found sufficient for the applications presented here.

In this paper, a three-dimensional approach to determine the residual stresses and distortions in a thin structure after a peening treatment is proposed. A combination of a semi-analytical method (SAM) and the finite element method (FEM) is used. The advantage of the SAM is that it allows a considerable saving of computing time compared to the FEM, however it relies on the assumption of half-space for the bodies in contact, which is not true for thin structures. A methodology is proposed here to account for thin structures, i.e. when the plastically deformed layer is of thickness significant compared to the structure thickness. The method is applied to the USP treatment of thin plates but can easily be extended to other peening processes and impacted structures.

First, an experimental study of the USP process on thin plates made of IN600 is presented. In a second step, the numerical multiple impacts model is proposed. The SAM is used to determine the plastic strain field generated by the USP process that is averaged over a local representative area of the half-space. Then, this inelastic strain field is transferred into a FE model (using the commercial package Systus/Sysweld; SYSWELD, 2008) to predict the residual stress state in the thin plate. Finally, the numerical and experimental results are compared, showing the capability of the model to predict residual stresses and distortions of a peened structure.

Experimental procedure

Experimental device and principle

The USP device used in these experiments is illustrated in Fig. 1(a). A sinusoidal electric field is delivered by a generator. By using a piezo-transmitter, this electric energy is transformed into an ultrasonic vibration and then amplified by a booster. This energy is thereafter transferred into kinetic energy stored by the balls when they are projected by the sonotrode. The balls are confined in a chamber on top of which the specimen is fixed. Similar to the "Almen test", the specimen is fixed by four bolts as illustrated in Fig. 1(b). During the process, the successive impacts will generate a compressive residual stress field at the surface and in the close sublayer of the treated plate (Fig. 1(c)). A tensile residual stress field of lower amplitude balances the compressive effect at higher depth to ensure the equilibrium. At the end of the process, a small curvature of the plate can be observed. After the removal of the bolts, the residual stress field is redistributed due to the change of boundary conditions. The convexity of the surface is thus increased in the direction of the peened surface. The maximum deflection is used industrially to characterize the peening processes.

Several parameters can influence the results of the experiments. They are mainly related to the shots (material and diameter) and to the process (shots velocities and peening coverage rate). The corresponding parameters used in our experiments are presented thereafter.

Materials and shot peening parameters

In this paper, the USP experiments are carried out on thin plates made of IN600 (initial yield strength 0 = 215 MPa, Young's modulus E = 195.2 GPa, Poisson's ratio = 0.3, density = 8250 kg/m 3 ). The bearing steel AISI52100 ( 0 = 1618 MPa, E = 213 GPa, = 0.29, = 7834 kg/m 3 ) was chosen for the shots with a diameter of 4 mm. The large difference between the materials' yield strengths allows to consider the shots as purely elastic.

For the peening device, the piezoelectric transducer emits the ultrasonic wave at 20 kHz. The amplitude of displacement of the sonotrode is 25 m. Furthermore, the dimension of the chamber is 85 mm × 35 mm × 50 mm and new balls are used for each test. A reduced number of balls, compared to traditional applications, were used to allow a fair estimation of the shots velocity as indicated further.

The treated surfaces are initially painted in order to determine the coverage rate. Experimentally, the coverage rate is defined as the percentage of the surface that has been impacted at least once. A 100% coverage is defined when the whole surface has been impacted. The coverage is then equal to 200% when the exposure time is doubled.

Residual stresses, initial state and measurements

All specimens are machined by wire electrical discharge machining and both faces have been grinded. Then, the specimens have been submitted to a stress relieving heat treatment (SRHT, 3 h at 600 • C) to reduce as much as possible the stresses caused by machining. The SHRT temperature and maintaining time have been chosen not to modify the microstructure and therefore the mechanical properties. However X-ray diffraction measurements before and after SHRT, prior the USP, showed that the level of residual stresses was not significantly affected by the SHRT. A residual stress level of the order of 100 MPa is still observed in the specimens.

For determination of residual stresses in IN600, the diffraction peaks at 114.13 • in the plane 3,3,1 are recorded using Cu Kr adiation, with a position sensitive detector (PSD). The parameters are the following: 40 kV and 30 mA. For measurement in the depth of the specimen, a small thickness is released using electropolishing on a 10 mm × 10 mm surface. A measurement error of ±50 MPa was estimated and is applied on all residual stress values. Due to the small thickness of the samples, the electro-polishing reduces the stiffness of the samples and causes a reorganization of the stresses. This effect can become significant at depths higher than 400 m. Correction techniques exist, as the one developed by [START_REF] Moore | Mathematical correction for stress in removed layers in X-ray diffraction residual stress analysis[END_REF]. However this technique is only valid for the removal of a uniform thickness over the whole specimen and was not applied here.

Impact velocity

The last important input parameter is the shot velocity. Since the movements of the balls in the chamber are random, the collisions between them or with the side walls of the chamber can change the direction and the velocity of the balls. It is assumed here that all the impacts will occur at the same velocity and an average velocity is considered. This average velocity is taken to be the maximum initial velocity of the sonotrode. Considering that the sinusoidal harmonic signal delivered by the generator is:

x (t) = Asinωt (1)
where A is the amplitude and ω represents the pulsation. It yields:

dx dt = Aωcosωt (2)
According to Eq. ( 2), the maximum initial velocity is:

V inimax = 2A f (3)
where f represents the frequency of vibration of the sonotrode. The initial theoretical shot velocity calculated with our experimental parameters is equal to 4 m/s and corresponds to the average capabilities of the USP device used for the experiments. This value is considered as the average velocity of the impacts during the process and will be used as an input in the simulations. The real impact velocity in USP is difficult to estimate due to the complex kinematics of the shots in the peening chamber. A value of 4 m/s is in the range of velocities estimated by other authors, as [START_REF] Todaka | Comparison of nanocrystalline surface layers in steels formed by air blast and ultrasonic shot peening[END_REF] that indicated that the shot velocity in USP is lower than 20 m/s. [START_REF] Pilé | Modelling of the ultrasonic shot peening process[END_REF] have established a model of the USP process and obtained a shot mean speed of the order of 3 m/s.

Due to the large number of impacts between shots and with the chamber's walls during a USP process, the impacts are largely multidirectional. The estimation of the impact velocity is then complicated. Only 10 shots were used in the experiments to reduce the number of impacts between shots. The impacts between the shots during the experiments can then be considered mainly normal to the impacted surface and the previous approximation used to estimate the impact velocity. Obviously, in the case of a real USP process, the kinematics of the shots would become more complicated and a specific model would be required to determine the impact velocity distribution, as indicated in the introduction. Such a model would allow to obtain both the average impact velocities and the velocity distribution.

Numerical model

The main problem of multiple impact simulations is the computation time they require. For this purpose, it is proposed here to use a semi-analytical method (SAM) to model the impacts that has the advantage to be about one order of magnitude faster than classical finite element method.

The semi-analytical method

When analytical solutions cannot be found for complex problems, it is possible to divide these problems in a sum of basic solutions for which analytical solutions are known. Such analytical solutions are derived from the boundary elements technique. Based on this method, a semi-analytical code for the study of elastic-plastic contact problems has been developed. The initial development of the problem was made by [START_REF] Jacq | Development of a threedimensional semi-analytical elastic-plastic contact code[END_REF] and subsequent development on the plasticity can be found in [START_REF] Nélias | Elastic-plastic contact between rough surfaces: proposal for a wear or running-in model[END_REF]. The method was then extended to solve impact problems (see Chaise et al., 2011). Note that, at the impact velocities encountered during USP process, it could be assumed that inertial effects can be neglected since the impact velocity, here equal to 4 m/s, is far lower than the speed of elastic waves in metals (3000 m/s in aluminium, 5200 m/s in steel). Under these conditions, [START_REF] Johnson | Contact Mechanics[END_REF] showed that the impact problem can be treated as quasi-static. The impact velocity observed in conventional shot peening being seldom higher than 150 m/s, a quasi-static model can be used to model the shot peening process. This is confirmed by the work of Rouhaud et al. (2005) who showed that a quasi-static model can be used, in their case for impacts at 75 m/s. In addition the bodies in contact can be considered as half-space since the contact dimensions are small compared to the dimensions of the bodies.

Stick-slip should occur within the contact area for frictional normal impact (see [START_REF] Gallego | A fast and efficient contact algorithm for fretting problems applied to fretting modes I, II and III[END_REF], as well as gross slip for oblique contact. These effects will not be considered here. In what will follow the impact will be assumed frictionless and normal to the surface. Note that for frictionless contact, it is also equivalent to an oblique impact but with the projection of the velocity vector on the normal to the surface taken as the impact velocity.

The general algorithm of the semi-analytical impact model is presented in Fig. 2. The initial state is first defined by the shot velocity, shot radius and mechanical properties of the shot and the body. At each time step, the rigid body displacement is determined from the current shot velocity. Then the displacement is used as input to a classical elastic normal contact solver. The pressure field resulting from this solver is used to calculate the so-called elastic stresses using the elementary analytical solutions derived by [START_REF] Love | A Treatise on the Mathematical Theory of Elasticity[END_REF]. Initial plastic strains are used to calculate the residual stresses (based on the solutions developed by [START_REF] Chiu | On the stress fiels and surface deformation in a half space with a cuboidal zone in which initial strains are uniform[END_REF]. The total stress is calculated as the sum of the two former stresses and used as input to a classical return-mapping algorithm to calculate the increment of plastic strain in the body.

From this increment, new residual stresses are computed along with the increment of surface residual displacement. The surface definition is so updated and the contact problem recomputed. This loop is reproduced until convergence of both the contact and plastic problems are obtained. The convergence of the contact problem is checked by the convergence of the increment of surface displacement while the value of the yield function at each point of the body is used to validate the convergence of the plastic problem. When the convergence is reached, the final contact force is used to update the velocity of the ball and the next rigid body displacement can be calculated. For more details on SAM as developed to solve inelastic contact problems or contact problem with heterogeneous materials the reader may refer to [START_REF] Boucly | Modeling of the rolling and sliding contact between two asperities[END_REF], [START_REF] Fulleringer | On the tangential displacement of a surface point due to a cuboid of uniform plastic strain in a half-space[END_REF], [START_REF] Leroux | Stick-slip analysis of a circular point contact between a rigid sphere and a flat unidirectional composite with cylindrical fibers[END_REF] or Chaise and Nélias (2011). The reader may also refer to the recent work of Wang, Keer and coauthors from Northwestern University [START_REF] Chen | Transient thermomechanical analysis of sliding electrical contacts of elastoplastic bodies, thermal softening, and melting inception[END_REF][START_REF] Chen | Thermomechanical analysis of elasto-plastic bodies in a sliding spherical contact and the effects of sliding speed, heat partition, and thermal softening[END_REF][START_REF] Zhou | A fast method for solving three-dimensional arbitrarily-shaped inclusions in a half space[END_REF], Bosman and Schipper at the University of Twente [START_REF] Bosman | Transition from mild to severe wear including running-in effects[END_REF], and Wang et al. from Tsinghua Univ. [START_REF] Wang | Partial slip contact analysis on three-dimensional elastic layered half space[END_REF].

At the end of the rebound, the final residual stress and plastic strain tensors resulting from the impact are obtained.

Benchmark validation

A benchmark was realized between the following finite element codes: Abaqus, Code Aster, Systus/Sysweld and the SAM. The aim here is to validate the results obtained from an impact simulation by the SAM and to compare the respective performances of the codes. The test case presented here corresponds to a more general shot peening case and input data is then different than the one of the USP case. Nevertheless, the amount of plastic strain obtained in both cases is similar in order to validate the chosen framework.

Input data

A single impact of an elastic ball (Young's modulus E = 207 GPa, Poisson's ratio = 0.3) on an elastic-plastic flat body (E = 111 GPa, = 0.29) is considered. The impacted body hardening law is isotropic, governed by a Ramberg-Osgood law described as follows:

= 0 + kε n p (4)
with 0 = 478 MPa, k = 3270.6 MPa and n = 0.441.

The ball diameter is 0.6 mm and the impact velocity is taken equal to 40 m/s. Those conditions correspond to the case of classical shot peening of a TiAl6V material body with AISI52100 shots. The initial velocity is perpendicular to the impacted surface and the contact is considered frictionless.

Calculation conditions

In all models, the contact area was discretized with elements of dimensions 8.5 mm × 8.5 mm × 4.25 mm. Out of this area, the mesh is defined differently for each model, each one running on a different software. Linear brick elements with reduced integration were adopted in the FE codes. For the SAM, cuboidal elements of constant size are defined. Each element has a single calculation point (equivalent to a Gauss point in FEM) located at its centre. The Abaqus and Code Aster models present one plane of symmetry, while Systus/Sysweld presents one quarter of the problem only. The FE models are presented in Fig. 3. For the SAM, no symmetry is considered and the full 3D problem is modelled. Furthermore, the element size in the SAM has to be constant to perform the computation and a fine mesh is so adopted over all the discretized area. Table 1 presents the number of nodes, type of elements and CPU frequency used for each code. The approximate computation time for each calculation is also given. The differences in computation time between the FE codes are mainly elated to the number of symmetries considered. More importantly, the SAM, despite a larger number of elements, resents a computation time about one order of magnitude lower than the FEM.

Benchmark results

The normal surface residual displacement, describing the dent after impact, obtained for each code is plotted in Fig. 4. This displacement u res z and the radial coordinate x are normalized by the radius of contact a*. The plastic train components along the depth taken at the centre of the impact are plotted in Fig. 5. The results obtained with the different codes show a good agreement both at the surface and along depth. It can be noticed that the zz component is the most sensitive and the one for which the largest differences are observed between codes. The SAM can then be considered accurate to predict plastic strains even at relatively high strain amplitudes (here about 10%).

Table 1

Numerical parameters for the single impact benchmark between the semi-analytical method, Abaqus, code ASTER and Systus/Sysweld.

Finite element codes

Semi-analytical method

Abaqus Code Aster Systus

Element size in contact area The SAM presents a significantly reduced computation time despite a higher number of elements. This increase in he number of elements is due to the constant size of the SAM elements, required for the FFT calculation. The fine mesh size in the contact area has to be the same over the whole computation domain even in the zones with lower strain and stress gradients. The Systus/Sysweld results were computed in reduced CPU time compared to the other FE codes, it should be noticed that only a quarter of the domain was modelled while Abaqus and Code Aster models describe one half of the problem. The SAM models the full 3D domain, without symmetry simplifications. Note that what is costly in terms of computing time with SAM is the number of elements in the plastic zone. Despite this drawback, the computation with the SAM is significantly faster (see Table 1). Meanwhile, the FE codes present differences in the computation time despite a similar number of nodes. The construction of the model and the FE software being different for each model, it led to significant differences on the computation time. Furthermore, a mesh sensitivity analysis was performed on the SAM. A mesh size twice coarser than the one used in this benchmark provides similar results.

x = y = 8.5 m a*/ x ≈ 10 z = x/2 (in

Multiple impacts model

Model hypotheses

Multiple impact models can be sorted in two categories: deterministic and stochastic ones, see Section 1. The model proposed here is based on a deterministic approach. The impacts' velocities and distribution are set a priori.

Each impact is simulated using the SAM. Based on the studied range of velocity, the model is quasi-static. The current model considers only impacts normal to the surface to simplify the estimation of the impact velocity. As long as the number of shots is small and frictional effects can be neglected, this hypothesis seems reasonable.

The influence of constant impact velocity is discussed in Section 3.3.5. The impact distribution is described by the coverage rate. The definition of the coverage rate estimated experimentally does not match the definition of the simulated one. It is assumed that experimental and numerical results with similar values of coverage are equivalent. Furthermore, a random sequence of impacts, typically obtained in experiments is supposed to produce results equivalent to those resulting from a defined impact sequence as long as the coverages are equal and after averaging the results over a representative area as defined in Section 3.3.3. This hypothesis is very strong but allows to compare in the most simple way experimental and numerical results.

Impact distribution and coverage rate

It is assumed that after a uniform shot peening treatment on a surface, the residual stress and strain state have specific characteristics. The strain state is supposed uniform in planes parallel to the surface and isotropic. Consequently the plastic strain varies only along the depth of the specimen (z coordinate) and is independent of he position at the surface (x and y coordinates). This is assumed valid far enough from the boundaries of the art, based on Saint-Venant's principle. Consequently, instead of the whole treated surface a local area is considered to calculate the residual strain and stress field in order to minimize the number of impacts to simulate.

In this study, the impacts are supposed to occur at the three corners of an equilateral triangle the side length of which is noted d (see Fig. 6). Taking the triangle formed by the centre of three impacts, a reference area A r is defined. The coverage rate c can be calculated as the ratio between the surface covered by the impacts over the reference area as follows:

c = 2 √ 3 a * d 2 (5)
where a* is the radius of the maximum area of contact for a single elastic plastic impact. Note that this value is only slightly altered by the presence of former impacts and is so representative for all impacts occurring during the process. The value a* is used for the coverage definition since it defines the surface on which a pressure has actually been applied which corresponds experimentally to the surface where painting would have been removed.

A significant difference appears here between the definitions of the experimental and numerical coverages. The experimental value considers that each point of the surface has been impacted at least once when the coverage equal 100%. Meanwhile, the numerical value will assume a 100% coverage while some points of the surface have not been impacted but some others have been impacted at least twice. This difference has no critical influence on he results presented here.

Averaging process

For each depth z, the plastic strain tensor components are averaged over the reference area as indicated in Eq. ( 6). This average is representative of the homogeneous strain values generated by the USP process described in the former section. This hypothesis is supported by the area-averaged approach firstly developed by [START_REF] Menig | Depth profiles of macro residual stresses in thin shot peened steel plates determined by X-ray and neutron diffraction[END_REF] and recently applied on a square pattern by [START_REF] Kim | An area-average approach to peening residual stress under multi-impacts using a three-dimensional symmetry-cell finite element model with plastic shots[END_REF] which indicated that the area-averaged solution is very close to the XRD solution.

ε p ij (z) = 1 N r (x,y)∈Ar ε p ij (x, y, z) (6)
where N r is the number of elements in the reference area A r .

A triangular repartition of the impacts was chosen to minimize the number of simulations. It allows with only three impacts to obtain a fair representation of the full shot peened state but it does present only one plane of symmetry. Different distributions, like a square one, would present two symmetries and so a diagonal averaged plastic strain tensor. It is indeed expected that for an isotropic fully shot peened state, the results taken at any point of the surface would be constant, independent of the coordinate system and without off-diagonal terms. The averaged plastic strain tensor is thus modified as follows which is equivalent to average the results over four simulations with different orientations of the triangle obtained by a rotation of /2, and 3 /2:

ε p final = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ε p xx + ε p yy 2 0 0 0 ε p xx + ε p yy 2 0 0 0 ε p zz ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (7)
The results, depending only on the depth z, are consistent with the hypothesis made on the residual strain field in a fully shot peened body. Simulations performed on a square distribution of the impacts showed that the pattern of impacts has a little influence on the results that are thus mainly governed by the coverage rate.

Representative number of impacts

In the literature, the number of impacts required to obtain stabilized results appears to be a critical parameter. The residual stresses are highly influenced by the impacts surrounding the reference area. On the contrary, the plastic strains are influenced by only a small number of impacts around the reference area, leading to eigenstrain methods as applied by [START_REF] Achintha | Eigenstrain modelling of residual stresses generated by laser shock peening[END_REF] to laser shock peening. The influence of impact beyond the three reference ones is studied here.

The case of a perfectly plastic body with a 150% coverage rate subjected to 6 impacts distributed as indicated in Fig. 6 is studied (Table 2). This case can be considered as a limit case. As no hardening is involved the evolution of plastic strain from one impact to the other is the highest. In the case of a 150% coverage rate, the contact radius a* is approximately equal to 2/3 of the distance between impacts. With very large overlapping of the impacts, the influence of each impact on the previous impact's plastic strain is the highest.

For this case (Fig. 7), only a slight increase of the averaged accumulated plastic strain is observed between the 3rd and 6th impacts. Three impacts are then considered as a very reasonable, and significantly lower compared to the former models [START_REF] Majzoobi | A three-dimensional simulation of shot peening process using multiple shot impacts[END_REF], number of impacts required to accurately describe the averaged plastic strain profile. This limit value is valid for the current model only, i.e. for impacts regularly spaced out on a triangular grid and occurring at constant velocity. For stochastic models typically, as the one proposed by [START_REF] Bagherifard | A numerical model of severe shot peening (SSP) to predict the generation of a nanostructured surface layer of material[END_REF] subsequent impacts might occur in the reference area and increase the plastic strain level but the coverage rate value would also increase.

Influence of the impact velocity distribution

The influence of a variation of the impact velocity around an average value is studied here. Three simulations of three impacts are made for two cases. The first simulation considers three impacts at a constant reference velocity v ref , corresponding to a reference kinetic energy E kinref . For the second simulation, the first impact velocity is lower than the reference velocity and such that the kinetic energy of the first impact is 20% lower than E kinref ; the second impact velocity is equal to v ref and the last impact velocity is higher than v ref with a kinetic energy 20% higher than E kinref . The third simulation considers the same impact velocities but this time the impact velocity decreases for each impact. For each simulation the total kinetic energy brought to the body through the three impacts is constant and equal to 3E kinref . The remaining impact parameters are indicated in Table 3 and the results post-processed as usual are presented in Fig. 8.

Results for these three sets of simulations are almost the same. A similar set of simulations is performed with a larger kinetic energy variation (up to 50%) between each impact, results are also plotted in Fig. 8.

For this larger variation case, the maximum (averaged value) accumulated plastic strain is lower but the profile of residual stress remains almost constant. A rebalancing of the residual stresses is observed at deep depth in this case.

From these simulations it can be concluded that a variation of the impact velocity around an average constant value has a negligible effect on the results as soon as the total kinetic energy brought to the system is constant. The use of a constant velocity for the simulations is thus validated.

3.4. From plastic strains in a semi-infinite body to residual stresses in a thin structure

The half-space hypothesis considered in the SAM being unsuitable to predict the residual state within thin structures treated by USP, a technique to transfer the plastic strain field into a thin structure is proposed here.

The averaged plastic strain tensor previously defined is only dependent on the depth and is considered representative of the final state of the plate (see Section 3.3.2). This profile is transferred to a finite element model of the plate, or any other structure, the tensors values being constant in each plane parallel to the surface. These averaged strains are used as initial plastic or inelastic strain components at the Gauss points of the FE model.

The FE analysis is done with Systus/Sysweld. In most commercial FE packages it is not possible to import an inelastic strain tensor as initial state, whereas equivalent thermal strains can be used to represent the original state as presented by [START_REF] Fulleringer | On the tangential displacement of a surface point due to a cuboid of uniform plastic strain in a half-space[END_REF]. The mesh is set up by 8-node linear bricks. No external load is applied to the plate. The material parameters of the plate in the FE model are equal to those of the SAM. The boundaries of the plate are kept free, thus representing the state when the bolts are removed. The residual stress state in the plate after the peening process is obtained by solving the static equilibrium of the plate (see Fig. 9) which accommodates the plastic (or inelastic) strains. These residual stresses are observed at the central point of the structure far from the boundaries where it is assumed that side effects will influence the results.

Results and discussion

Experimental data are compared with numerical results. For each simulation three impacts are performed. The shots velocity is assumed constant and equal to the maximum velocity of the sonotrode (4 m/s). The AISI52100 properties are considered for the shots that are assumed to behave elastically. The IN600 hardening behaviour is described by a Ramberg-Osgood hardening law deduced from monotonic tensile tests with K = 1450 MPa and n = 0.8. Simulation results are averaged following the multiple impact method. The radial averaged residual stress obtained in the half-space (i.e. by SAM prior to the transfer of inelastic strains to the thin structure) will be also indicated in Fig. 10(a)-(c) presented hereafter.

The averaged plastic strain is then exported in a FE model to determine the residual stresses in the finite structure (see Fig. 10(a)-(c)) and the deflection of the samples (Fig. 11). The residual stresses in Fig. 10(a)-(c) are plotted along the depth of the plate at its centre. Fig. 11 shows the deformed profiles of the thin plate in the plane y = 0.

The residual stresses at the centre of the plate are measured in the longitudinal ( xx ) and transverse ( yy ) directions. In addition, the residual arc height of the specimen is measured by a 3D coordinate-measuring machine in the longitudinal direction of the plate and along its centre.

The residual stress profiles are quite similar for all coverages, both from experiments and simulations. Measured and simulated stresses are plotted in Figs. 10(a)-(c). For all processed specimens the stresses are compressive both at the surface and in the top sublayer. Moreover, it can be seen that the transverse stresses ( yy ) are higher than the longitudinal ones ( xx ). This can be explained by the clamping conditions of the specimen. Indeed, the specimen can sustain higher deformation in the longitudinal direction than the transverse one. Note that the gap between xx and yy is the same for the simulations and the measurements.

A significant difference can be observed between measured and calculated residual stresses in the top surface layer (up to 100 m). This difference might be caused by the machining process, indeed EDM typically leaves strong tensile residual stresses (potentially up to a few hundreds MPa) in the subsurface layer up to 100-200 m, see for example, [START_REF] Ekmekci | A semi-empirical approach for residual stresses in electric discharge machining (EDM)[END_REF]. The samples were considered initially free of plastic strain while the stress relieving treatment was not efficient enough to reduce the initial hardening due to surface machining. The measurement points at depth ranging from 0.4 to 0.5 mm (i.e. one fourth of the plate thickness) diverge somewhat from the numerical results in Fig. 10(b) (coverage rate of 150%) and Fig. 10(c) (coverage rate of 200%). As indicated in the XRD measurement conditions, this is probably due to the loss of rigidity of the plate that was not corrected.

An excellent fit is obtained between simulation and experimental results for the deflection of the plate (see Fig. 11). The maximum measured value is nicely reproduced by the model. Results tend to diverge at the border of the plate were the actual boundary conditions logically differ from those used in the simulations. For this case, no saturation of the deflection is observed for the three values of coverage rate.

Conclusion

A multiple impacts model to predict the result of an USP treatment on a thin structure was proposed. This model is based on a semi-analytical method which has the advantage of saving considerable computing time compared to the classical finite element method, without sacrificing the quality of results. An averaging and transfer process to a FE model allows considering any kind of structures, despite the semi-infinite body hypothesis of the SAM. The averaging of the results allows simulating only three impacts to simulate a full shot peening process and to bridge the results with the actual experimental measures that are also averaged over a given, yet not necessarily representative, volume. Predicting the results of a full peening process with fine mesh description becomes affordable within a few hours on a standard desktop with a relatively fine mesh.

The model was successfully applied to thin structures. The comparison between the numerical and experimental results shows a good agreement for the distortions and correct trends for the residual stresses. Several parameters remain to be controlled for an accurate simulation of the process. In particular, the impact velocity, that has a direct influence on the level of plastic strains, was here only roughly estimated. A finer description of this velocity would be necessary for a full analysis of the process and will be of concern in future works. Nevertheless it was shown that the determination of an average impact velocity is sufficient to obtain realistic results on global structures.

Fig. 1 .

 1 Fig. 1. Experimental device: (a) overall view; (b) fixation system of the specimen by bolts; and (c) typical residual stress profile due to USP treatment.

Fig. 2 .

 2 Fig. 2. General algorithm of the semi-analytical impact model.

Fig. 3 .

 3 Fig. 3. Finite element model mesh definition and zoom around the contact area for the Abaqus model (a), Code Aster model (b) and Systus/Sysweld model (c).

Fig. 4 .

 4 Fig. 4. Benchmark results: residual surface displacement u res z as a function of the radial coordinate x, normalized by the area of contact radius a*.

Fig. 5 .

 5 Fig.5. Benchmark results: accumulated plastic strain at the centre of the impact as a function of the depth z.

Fig. 6 .

 6 Fig.6. Distribution of impacts in a triangular scheme, case of a coverage rate equal to 150%. Definition of the reference area with the centre of the three first impacts and subsequent impact positions.

Fig. 7 .

 7 Fig.7. Evolution of the accumulated plastic strain p normalized by the yield strain (0.002) for 6 impacts on a perfectly plastic body (coverage rate c = 150%).

Fig. 8 .

 8 Fig. 8. Influence of a variation of the impact's kinetic energy on averaged results for 3 impacts on a linear isotropic hardening body (ET/E = 10%): (a) averaged accumulated plastic strain and (b) Von Mises equivalent of averaged residual stress components.

Fig. 9 .Fig. 10 .

 910 Fig.9. Averaging and transfer to finite structure procedure. Averaging of results over the semi-infinite body reference area, exportation as initial inelastic strains in a finite element model of a finite size structure and computation of residual stresses by solving the equilibrium.

Fig. 11 .

 11 Fig. 11. Almen-type experimental results for a maximum sonotrode velocity v = 4 m/s. Comparison of the deformed shape by simulation and after experiments for 3 coverage rates.

Table 2

 2 Simulation parameter for the study of the effect of more than three impacts.

	Number of impacts	Impact velocity [m/s]	Yield limit y [MPa]	Tangent modulus ET/E	Coverage rate c [%]
	6	10	400	0	150

Table 3

 3 Simulation parameters for the study of the effect of a variation of impact velocity.

	Number of impacts	Impact velocity [m/s]	Yield limit y [MPa]	Tangent modulus ET/E	Coverage rate c [%]
	3	[20 ± 20% E kin ]	800	10	100
	3	[20 ± 50% E kin ]	800	10	100