Ahmet Iscen

Yannis Avrithis

Giorgos Tolias

Teddy Furon

Ondřej Chum

Hybrid Diffusion: Spectral-Temporal Graph Filtering for Manifold Ranking

State of the art image retrieval performance is achieved with CNN features and manifold ranking using a k-NN similarity graph that is pre-computed off-line. The two most successful existing approaches are temporal filtering, where manifold ranking amounts to solving a sparse linear system online, and spectral filtering, where eigen-decomposition of the adjacency matrix is performed off-line and then manifold ranking amounts to dot-product search online. The former suffers from expensive queries and the latter from significant space overhead. Here we introduce a novel, theoretically well-founded hybrid filtering approach allowing full control of the space-time trade-off between these two extremes. Experimentally, we verify that our hybrid method delivers results on par with the state of the art, with lower memory demands compared to spectral filtering approaches and faster compared to temporal filtering.

Introduction

Most image retrieval methods obtain their initial ranking of the database images by computing similarity between the query descriptor and descriptors of the database images. Descriptors based on local features [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF][START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] have been largely replaced by more efficient CNN-based image descriptors [START_REF] Gordo | End-to-end learning of deep visual representations for image retrieval[END_REF][START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF]. Regardless of the initial ranking, the retrieval performance is commonly boosted by considering the manifold structure of the database descriptors, rather than just independent distances of query to database images. Examples are query expansion [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF][START_REF] Arandjelovic | Three things everyone should know to improve object retrieval[END_REF] and diffusion [START_REF] Donoser | Diffusion processes for retrieval revisited[END_REF][START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF][START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Query expansion uses the results of initial ranking to issue a novel, enriched, query [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF] on-line only. Diffusion on the other hand, is based on the k-NN graph of the dataset that is constructed off-line, so that, assuming novel queries are part of the dataset, their results are essentially pre-computed. Diffusion can then be seen as infinite-order query expansion [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF].

The significance of the performance boost achieved by diffusion has been recently demonstrated at the "Large-Scale Landmark Recognition"1 challenge in conjunction with CVPR 2018. The vast majority of top-ranked teams have used query expansion or diffusion as the last step of their method.

Recently, efficient diffusion methods have been introduced to the image retrieval community. Iscen et al . [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] apply diffusion to obtain the final ranking, in particular by solving a large and sparse system of linear equations. Even though an efficient conjugate gradient (CG) [START_REF] Hackbusch | Iterative solution of large sparse systems of equations[END_REF] solver is used, query times on large-scale datasets are in a range of several seconds. A significant speed-up is achieved by truncating the system of linear equations. Such an approximation, however, brings a slight degradation in the retrieval performance. Their method can be interpreted as graph filtering in the temporal domain.

In the recent work of Iscen et al . [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF], more computation is shifted to the offline phase to accelerate the query. The solution of the linear system is estimated by low-rank approximation of the k-NN graph Laplacian. Since the eigenvectors of the Laplacian represent a Fourier basis of the graph, this is interpreted as graph filtering in the spectral domain. The price to pay is increased space complexity to store the embeddings of the database descriptors. For comparable performance, a 5k-10k dimensional vector is needed per image.

In this paper, we introduce a hybrid method that combines spectral filtering [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] and temporal filtering [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. This hybrid method offers a trade-off between speed (i.e., the number of iterations of CG) and the additional memory required (i.e., the dimensionality of the embedding). The two approaches [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF][START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] are extreme cases of our hybrid method. We show that the proposed method pairs or outperforms the previous methods while either requiring less memory or being significantly faster -only three to five iterations of CG are necessary for embeddings of 100 to 500 dimensions.

While both temporal and spectral filtering approaches were known in other scientific fields before being successfully applied to image retrieval, to our knowledge the proposed method is novel and can be applied to other domains.

The rest of the paper is organized as follows. Related work is reviewed in Section 2. Previous work on temporal and spectral filtering is detailed in Sections 3.1 and 3.2 respectively, since the paper builds on this work. The proposed method is described in Section 4 and its behavior is analyzed in Section 5. Experimental results are provided in Section 6. Conclusions are drawn in Section 7.

Related work

Query expansion (QE) has been a standard way to improve recall of image retrieval since the work of Chum et al . [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF]. A variety of approaches exploit local feature matching and perform various types of verification. Such matching ranges from selective kernel matching [START_REF] Tolias | Visual query expansion with or without geometry: refining local descriptors by feature aggregation[END_REF] to geometric consensus [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF][START_REF] Chum | Total recall II: Query expansion revisited[END_REF][START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF] with RANSAC-like techniques. The verified images are then used to refine the global or local image representation of a novel query.

Another family of QE methods are more generic and simply assume a global image descriptor [14-18, 7, 6]. A simple and popular one is average-QE [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF], recently extended to α-QE [START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF]. At some small extra cost, recall is significantly boosted. This additional cost is restricted to the on-line query phase. This is in contrast to another family of approaches that considers an off-line pre-processing of the database. Given the nearest neighbors list for database images, QE is performed by adapting the local similarity measure [START_REF] Jégou | A contextual dissimilarity measure for accurate and efficient image search[END_REF], using reciprocity con-straints [START_REF] Delvinioti | Image retrieval with reciprocal and shared nearest neighbors[END_REF][START_REF] Qin | Hello neighbor: Accurate object retrieval with k-reciprocal nearest neighbors[END_REF] or graph-based similarity propagation [START_REF] Zhang | Query specific fusion for image retrieval[END_REF][START_REF] Donoser | Diffusion processes for retrieval revisited[END_REF][START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. The graph-based approaches, also known as diffusion, are shown to achieve great performance [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] and to be a good way for feature fusion [START_REF] Zhang | Query specific fusion for image retrieval[END_REF]. Such on-line re-ranking is typically orders of magnitude more costly than simple average-QE.

The advent of CNN-based features, especially global image descriptors, made QE even more attractive. Average-QE or α-QE are easily applicable and very effective with a variety of CNN-based descriptors [START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF][START_REF] Gordo | End-to-end learning of deep visual representations for image retrieval[END_REF][START_REF] Tolias | Particular object retrieval with integral maxpooling of cnn activations[END_REF][START_REF] Kalantidis | Cross-dimensional weighting for aggregated deep convolutional features[END_REF]. State-of-the-art performance is achieved with diffusion on global or regional descriptors [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. The latter is possible due to the small number of regions that are adequate to represent small objects, in contrast to thousands in the case of local features.

Diffusion based on tensor products can be attractive in terms of performance [START_REF] Bai | Regularized diffusion process on bidirectional context for object retrieval[END_REF][START_REF] Bai | Ensemble diffusion for retrieval[END_REF]. However, in this work, we focus on the page-rank like diffusion [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF][START_REF] Zhou | Ranking on data manifolds[END_REF] due to its reasonable query times. An iterative on-line solution was commonly preferred [START_REF] Donoser | Diffusion processes for retrieval revisited[END_REF] until the work of Iscen et al . [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF], who solve a linear system to speed up the process. Additional off-line pre-processing and the construction and storage of additional embeddings reduce diffusion to inner product search in the spectral ranking of Iscen et al . [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. This work lies exactly in between these two worlds and offers a trade-off exchanging memory for speed and vice versa.

Fast random walk with restart [START_REF] Tong | Fast random walk with restart and its applications[END_REF] is very relevant in the sense that it follows the same diffusion model as [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF][START_REF] Zhou | Ranking on data manifolds[END_REF] and is a hybrid method like ours. It first disconnects the graph into distinct components through clustering and then obtains a low-rank spectral approximation of the residual error. Apart from the additional complexity, parameters etc. of the off-line clustering process and the storage of both eigenvectors and a large inverse matrix, its online phase is also complex, involving the Woodbury matrix identity and several dense matrixvector multiplications. Compared to that, we first obtain a very low-rank spectral approximation of the original graph, and then solve a sparse linear system of the residual error. Thanks to orthogonality properties, the online phase is nearly as simple as the original one and significantly faster.

Problem formulation and background

The methods we consider are based on a nearest neighbor graph of a dataset of n items, represented by n × n adjacency matrix W . The graph is undirected and weighted according to similarity: W is sparse, symmetric, nonnegative and zero-diagonal. We symmetrically normalize W as W : = D -1/2 W D -1/2 , where D : = diag(W 1) is the diagonal degree matrix, containing the row-wise sum of W on its diagonal. The eigenvalues of W lie in [-1, 1] [START_REF] Chung | Spectral graph theory[END_REF].

At query time, we are given a sparse n × 1 observation vector y, which is constructed by searching for the nearest neighbors of a query item in the dataset and setting its nonzero entries to the corresponding similarities. The problem is to obtain an n × 1 ranking vector x such that retrieved items of the dataset are ranked by decreasing order of the elements of x. Vector x should be close to y but at the same time similar items are encouraged to have similar ranks in x, essentially by exploring the graph to retrieve more items.

Temporal filtering

Given a parameter α ∈ [0, 1), define the n × n regularized Laplacian function by

L α (A) : = (I n -αA)/(1 -α) (1)
for n × n real symmetric matrix A, where I n is the n × n identity matrix. Iscen et al . [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] then define x as the unique solution of the linear system

L α (W)x = y, (2)
which they obtain approximately in practice by a few iterations of the conjugate gradient (CG) method, since L α (W) is positive-definite. At large scale, they truncate W by keeping only the rows and columns corresponding to a fixed number of neighbors of the query, and re-normalize. Then, (2) only performs re-ranking within this neighbor set. We call this method temporal2 filtering because if x, y are seen as signals, then x is the result of applying a linear graph filter on y, and CG iteratively applies a set of recurrence relations that determine the filter. While W is computed and stored off-line, ( 2) is solved online (at query time), and this is expensive.

Spectral filtering

Linear system (2) can be written as x = h α (W)y, where the transfer function h α is defined by

h α (A) : = L α (A) -1 = (1 -α)(I n -αA) -1 (3)
for n × n real symmetric matrix A. Given the eigenvalue decomposition W = U ΛU of the symmetric matrix W, Iscen et al . [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] observe that h α (W) = U h α (Λ)U , so that (2) can be written as

x = U h α (Λ)U y, (4)
which they approximate by keeping only the largest r eigenvalues and the corresponding eigenvectors of W. This defines a low-rank approximation h α (W) ≈ U 1 h α (Λ 1 )U 1 instead. This method is referred to as fast spectral ranking (FSR) in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Crucially, Λ is a diagonal matrix, hence h α applies element-wise, as a scalar function

h α (x) : = (1 -α)/(1 -αx) for x ∈ [-1, 1].
At the expense of off-line computing and storing the n × r matrix U 1 and the r eigenvalues in Λ 1 , filtering is now computed as a sequence of low-rank matrix-vector multiplications, and the query is accelerated by orders of magnitude compared to [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. However, the space overhead is significant. To deal with approximation errors when r is small, a heuristic is introduced that gradually falls back to the similarity in the original space for items that are poorly represented, referred to as FSRw [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF].

We call this method spectral filtering because U represents the Fourier basis of the graph and the sequence of matrix-vector multiplications from right to left in the right-hand side of (4) represents the Fourier transform of y by U , filtering in the frequency domain by h α (Λ), and finally the inverse Fourier transform to obtain x in the time domain by U .

Hybrid spectral-temporal filtering

Temporal filtering ( 2) is performed once for every new query represented by y, but W represents the dataset and is fixed. Could CG be accelerated if we had some very limited additional information on W?

On the other extreme, spectral filtering (4) needs a large number of eigenvectors and eigenvalues of W to provide a high quality approximation, but always leaves some error. Could we reduce this space requirement by allocating some additional query time to recover the approximation error?

The answer is positive to both questions and in fact these are the two extreme cases of hybrid spectral-temporal filtering, which we formulate next.

Derivation

We begin with the eigenvalue decomposition W = U ΛU , which we partition as Λ = diag(Λ 1 , Λ 2 ) and U = (U 1 U 2 ). Matrices Λ 1 and Λ 2 are diagonal r × r and (n -r) × (n -r), respectively. Matrices U 1 and U 2 are n × r and n × (nr), respectively, and have the following orthogonality properties, all due to the orthogonality of U itself:

U 1 U 1 = I r , U 2 U 2 = I n-r , U 1 U 2 = O, U 1 U 1 + U 2 U 2 = I n . (5)
Then, W is decomposed as

W = U 1 Λ 1 U 1 + U 2 Λ 2 U 2 . (6)
Similarly, h α (W) is decomposed as

h α (W) = U h α (Λ)U (7) = U 1 h α (Λ 1 )U 1 + U 2 h α (Λ 2 )U 2 , (8)
which is due to the fact that diagonal matrix h α (Λ) is obtained element-wise, hence decomposed as

h α (Λ) = diag(h α (Λ 1 ), h α (Λ 2 )).
Here the first term is exactly the low-rank approximation that is used by spectral filtering, and the second is the approximation error

e α (W) : = U 2 h α (Λ 2 )U 2 (9) = (1 -α) U 2 (I n-r -αΛ 2 ) -1 U 2 + U 1 U 1 -U 1 U 1 (10) = (1 -α) U 2 (I n-r -αΛ 2 )U 2 + U 1 U 1 -1 -U 1 U 1 (11) = (1 -α) I n -αU 2 Λ 2 U 2 -1 -U 1 U 1 (12) = h α (U 2 Λ 2 U 2 ) -(1 -α)U 1 U 1 . (13)
We have used the definition (3) of h α in [START_REF] Hackbusch | Iterative solution of large sparse systems of equations[END_REF] and [START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF]. Equation ( 12) is due to the orthogonality properties [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF]. Equation [START_REF] Tolias | Visual query expansion with or without geometry: refining local descriptors by feature aggregation[END_REF] follows from the fact that for any invertible matrices A, B of conformable sizes,

(U 1 AU 1 + U 2 BU 2 ) -1 = U 1 A -1 U 1 + U 2 B -1 U 2 , (14)
which can be verified by direct multiplication, and is also due to orthogonality. Now, combining (8), ( 13) and ( 6), we have proved the following.

Theorem 1. Assuming the definition (3) of transfer function h α and the eigenvalue decomposition (6) of the symmetrically normalized adjacency matrix W, h α (W) is decomposed as

h α (W) = U 1 g α (Λ 1 )U 1 + h α (W -U 1 Λ 1 U 1 ), ( 15
)
where g α (A) : = h α (A) -h α (O) = (1 -α) (I n -αA) -1 -I n ( 16
)
for n × n real symmetric matrix A.

For x ∈ [-1, 1] in particular, g α (x) : = h α (x) -h α (0) = (1 -α)αx/(1 -αx).
Observe that Λ 2 , U 2 do not appear in [START_REF] Jégou | A contextual dissimilarity measure for accurate and efficient image search[END_REF] and indeed it is only the largest r eigenvalues Λ 1 and corresponding eigenvectors U 1 of W that we need to compute. The above derivation is generalized from h α to a much larger class of functions in appendix A.

Algorithm

Why is decomposition (15) of h α (W) important? Because given an observation vector y at query time, we can express the ranking vector x as

x = x s + x t , (17)
where the first, spectral, term x s is obtained by spectral filtering

x s = U 1 g α (Λ 1 )U 1 y, (18)
as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF], where g α applies element-wise, while the second, temporal, term x t is obtained by temporal filtering, that is, solving the linear system

L α (W -U 1 Λ 1 U 1 )x t = y, (19)
which we do by a few iterations of CG as in [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. The latter is possible because

L α (W -U 1 Λ 1 U 1 ) is still positive-definite, like L α (W).
It's also possible without an explicit dense representation of U 1 Λ 1 U 1 because CG, like all Krylov subspace methods, only needs black-box access to the matrix A of the linear system, that is, a mapping z → Az for z ∈ R n . For system [START_REF] Tolias | Particular object retrieval with integral maxpooling of cnn activations[END_REF] in particular, according to the definition (1) of L α , we use the mapping

z → z -α Wz -U 1 Λ 1 U 1 z /(1 -α), (20)
where product Wz is efficient because W is sparse as in [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF], while U 1 Λ 1 U 1 z is efficient if computed right-to-left because U 1 is an n × r matrix with r n and Λ 1 is diagonal as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF].

Discussion

What is there to gain from spectral-temporal decomposition (17) of x?

First, since the temporal term [START_REF] Tolias | Particular object retrieval with integral maxpooling of cnn activations[END_REF] can recover the spectral approximation error, the rank r of U 1 , Λ 1 in the spectral term [START_REF] Zhang | Query specific fusion for image retrieval[END_REF] can be chosen as small as we like. In the extreme case r = 0, the spectral term vanishes and we recover temporal filtering [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. This allows efficient computation of only a few eigenvectors/values, even with the Lanczos algorithm rather than the approximate method of [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Most importantly, it reduces significantly the space complexity, at the expense of query time. Like spectral filtering, it is possible to sparsify U 1 to compress the dataset embeddings and accelerate the queries online. In fact, we show in section 6 that sparsification is much more efficient in our case.

Second, the matrix W -U 1 Λ 1 U 1 is effectively like W with the r largest eigenvalues removed. This improves significantly the condition number of matrix L α (W -U 1 Λ 1 U 1 ) in the temporal term [START_REF] Tolias | Particular object retrieval with integral maxpooling of cnn activations[END_REF] compared to L α (W) in the linear system (2) of temporal filtering [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF], on which the convergence rate depends. In the extreme case r = n, the temporal term vanishes and we recover spectral filtering [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. In turn, even with small r, this reduces significantly the number of iterations needed for a given accuracy, at the expense of computing and storing U 1 , Λ 1 off-line as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. The improvement is a function of α and the spectrum of W, and is quantified in section 5.

In summary, for a given desired accuracy, we can choose the rank r of the spectral term and a corresponding number of iterations of the temporal term, determining a trade-off between the space needed for the eigenvectors (and the off-line cost to obtain them) and the (online) query time. Such choice is not possible with either spectral or temporal filtering alone: at large scale, the former may need too much space and the latter may be too slow.

Analysis

How "easier" for CG is W -

U 1 Λ 1 U 1 in (19) compared to W in (2)?
In solving a linear system Ax = b where A is an n × n symmetric positivedefinite matrix, CG generates a unique sequence of iterates x i , i = 0, 1, . . . , such that the A-norm e i A of the error e i : = x * -x i is minimized over the Krylov subspace K i : = b, Ab, . . . , A i b at each iteration i, where x * : = A -1 b is the exact solution and the A-norm is defined by x A : = √

x Ax for x ∈ R n . A well-known result on the rate of convergence of CG that assumes minimal knowledge of the eigenvalues of A states that the A-norm of the error at iteration i, relative to the A-norm of the initial error e 0 : = x * , is upper-bounded by [START_REF] Trefethen | Numerical linear algebra[END_REF]

e i A e 0 A ≤ φ i (A) : = 2 κ(A) -1 κ(A) + 1 i , (21)
where κ(A) : = A A -1 = λ 1 (A)/λ n (A) is the 2-norm condition number of A, and λ j (A) for j = 1, . . . , n are the eigenvalues of A in descending order.

In our case, matrix L α (W) of linear system (2) has condition number

κ(L α (W)) = 1 -αλ n (W) 1 -αλ 1 (W) = 1 -αλ n (W) 1 -α . ( 22
)
The first equality holds because for each eigenvalue λ of W there is a corresponding eigenvalue (1 -αλ)/(1 -α) of L α (W), which is a decreasing function.

The second holds because λ 1 (W) = 1 [START_REF] Chung | Spectral graph theory[END_REF]. Now, let W r : = W -U 1 Λ 1 U 1 for r = 0, 1, . . . , n -1, where Λ 1 , U 1 represent the largest r eigenvalues and the corresponding eigenvectors of W respectively. Clearly, W r has the same eigenvalues as W except for the largest r, which are replaced by zero. That is, λ 1 (W r ) = λ r+1 (W) and λ n (W r ) = λ n (W). The latter is due to the fact that λ n (W) ≤ -1/(n -1) ≤ 0 [START_REF] Chung | Spectral graph theory[END_REF], so the new zero eigenvalues do not affect the smallest one. Then,

κ(L α (W r )) = 1 -αλ n (W) 1 -αλ r+1 (W) ≤ κ(L α (W)). ( 23
)
This last expression generalizes [START_REF] Bai | Ensemble diffusion for retrieval[END_REF]. Indeed, W = W 0 . Then, our hybrid spectraltemporal filtering involves CG on L α (W r ) for r ≥ 0, compared to the baseline temporal filtering for r = 0. The inequality in ( 23) is due to the fact that |λ j (W)| ≤ 1 for j = 1, . . . , n [START_REF] Chung | Spectral graph theory[END_REF]. Removing the largest r eigenvalues of W clearly improves (decreases) the condition number of L α (W r ) relative to L α (W).

The improvement is dramatic given that α is close to 1 in practice. For α = 0.99 and λ r+1 (W) = 0.7 for instance, κ(L α (W r ))/κ(L α (W)) = 0.0326. More generally, given the eigenvalues λ r+1 (W) and λ n (W), the improvement can be estimated by measuring the upper bound φ i (L α (W r )) for different i and r. A concrete example is shown in Figure 1, where we measure the eigenvalues of the adjacency matrix W of a real dataset, remove the largest r for 0 ≤ r ≤ 300 and plot the upper bound φ i (L α (W r )) of the relative error as a function of rank r and iteration i given by ( 21) and [START_REF] Zhou | Ranking on data manifolds[END_REF]. Clearly, as more eigenvalues are removed, less CG iterations are needed to achieve the same relative error; the approximation error represented by the temporal term decreases and at the same time the linear system becomes easier to solve. Of course, iterations become more expensive as r increases; precise timings are given in section 6.

Experiments

In this section we evaluate our hybrid method on popular image retrieval benchmarks. We provide comparisons to baseline methods, analyze the trade-off between runtime complexity, memory footprint and search accuracy, and compare with the state of the art.

Experimental setup

Datasets. We use the revisited retrieval benchmark [START_REF] Radenović | Revisiting oxford and paris: Large-scale image retrieval benchmarking[END_REF] of the popular Oxford buildings [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] and Paris [START_REF] Philbin | Lost in quantization: Improving particular object retrieval in large scale image databases[END_REF] datasets, referred to as ROxford and RParis, respectively. Unless otherwise specified, we evaluate using the Medium setup and always report mean Average Precision (mAP). Large-scale experiments are conducted on ROxford +R1M and RParis +R1M by adding the new 1M challenging distractor set [START_REF] Radenović | Revisiting oxford and paris: Large-scale image retrieval benchmarking[END_REF].

Image representation. We use GeM descriptors [START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF] to represent images. We extract GeM at 3 different image scales, aggregate the 3 descriptors, and perform whitening, exactly as in [START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF]. Finally, each image is represented by a single vector with d = 2048 dimensions, since ResNet-101 architecture is used. Baseline methods. We consider the two baseline methods described in Section 3, namely temporal and spectral filtering. Temporal filtering corresponds to solving a linear system with CG [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] and is evaluated for different numbers of CG iterations. It is used with truncation at large scale to speed up the search [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] and is denoted by Temporal †. Spectral filtering corresponds to FSR and its FSRw variant [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Both FSR variants are parametrized by the rank r of the approximation, which is equal to the dimensionality of the spectral embedding. Implementation details. Temporal ranking is performed with the implementation3 provided by Iscen et al . [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. The adjacency matrix is constructed by using top k = 50 reciprocal neighbors. Pairwise similarity between descriptors v and z is estimated by (v z) 3 + . Parameter α is set to 0.99, while the observation vector y includes the top 5 neighbors. The eigendecomposition is performed on the largest connected component, as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Its size is 933,412 and 934,809 for ROxford +R1M and RParis +R1M, respectively. Timings are measured with Matlab implementation on a 4-core Intel Xeon 2.60GHz CPU with 200 GB of RAM. We only report timings for the diffusion part of the ranking and exclude the initial nearest neighbor search used to construct the observation vector.

Comparisons

Comparison with baseline methods. We first compare performance, query time and required memory for temporal, spectral, and hybrid ranking. With respect to the memory, all methods store the initial descriptors, i.e. one 2048dimensional vector per image. Temporal ranking additionally stores the sparse regularized Laplacian. Spectral ranking stores for each vector an additional embedding of dimensionality equal to rank r, which is a parameter of the method. Our hybrid method stores both the Laplacian and the embedding, but with significantly lower rank r.

We evaluate on ROxford +R1M and RParis +R1M with global image descriptors, which is a challenging large scale problem, i.e. large adjacency matrix, where prior methods fail or have serious drawbacks. Results are shown in Figure 2. Temporal ranking with CG is reaching saturation near 20 iterations as in [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. Spectral ranking is evaluated for a decomposition of rank r whose computation and required memory are reasonable and feasible on a single machine. Finally, the proposed hybrid method is evaluated for rank r = 400, which is a good compromise of speed and memory, as shown below.

Temporal [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] Temporal † [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] Spectral (FSRw) [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] Spectral ranking variants (FSR, FSRw) are not performing well despite requiring about 250% additional memory compared to nearest neighbor search in the original space. Compared to hybrid, more than one order of magnitude higher rank r is required for problems of this scale. Temporal ranking achieves good performance but at much more iterations and higher query times. Our hybrid solution provides a very reasonable space-time trade-off.

Runtime and memory trade-off. We report the trade-off between number of iterations and rank r, representing additional memory, in more detail in Figure 3. It is shown that the number of iterations to achieve the top mAP decreases as the rank increases. We achieve the optimal trade-off at r = 400 where we only need 5 or less iterations. Note that, the rank not only affects the memory and the query time of the spectral part in a linear manner, but the query time of the temporal part too [START_REF] Kalantidis | Cross-dimensional weighting for aggregated deep convolutional features[END_REF].

Sparsification of spectral embeddings is exploited in prior work [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. We sparsify the embeddings of our hybrid method by setting the smallest values of U 1 to zero until a desired level of sparsity is reached. We denote this method by SHybrid. This sparse variant provides memory savings and an additional speedup due to the computations with sparse matrices. Figure 4 shows that performance loss remains small even for extreme sparsity e.g. 99%, while the results in Figure 2 show that it offers a significant speedup in the global descriptor setup. Performance-memory-speed comparison with the baselines is shown in Table 1. Our hybrid approach enjoys query times lower than those of temporal with truncation or spectral with FSRw, while at the same time achieves higher performance and requires less memory than the spectral-only approach. We summarize our achievement in terms of mAP, required memory, and query time in Figure 5. Temporal ranking achieves high performance at the cost of high query time and its truncated counterpart saves query time but sacrifices performance. Spectral ranking is not effective at this scale, while our hybrid solution achieves high performance at low query times.

Comparison with the state of the art. We present an extensive comparison with existing methods in the literature for global descriptors at small and large scale (1M distractors). We choose r = 400 and 5 iterations for our hybrid method, 20 iterations for temporal ranking, r = 2k and r = 5k for spectral ranking at small and large scale, respectively. Temporal ranking is also performed with truncation on a shortlist size of 75k images at large scale. The comparison is presented in Table 2. Our hybrid approach performs the best or second best right after the temporal one, while providing much smaller query times at a small amount of additional required memory.

Conclusions

In this work we have tested the two most successful manifold ranking methods of temporal filtering [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] and spectral filtering [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] on the very challenging new benchmark of Oxford and Paris datasets [START_REF] Radenović | Revisiting oxford and paris: Large-scale image retrieval benchmarking[END_REF]. It is the first time that such methods are evaluated at the scale of one million images. Spectral filtering, with both its FSR and FSRw variants, fails at this scale, despite the significant space required for additional vector embeddings. It is possible that a higher rank would work, but it wouldn't be practical in terms of space. In terms of query time, temporal filtering is only practical with its truncated variant at this scale. It works pretty well in terms of performance, but the query time is still high. Our new hybrid filtering method allows for the first time to strike a reasonable balance between the two extremes. Without truncation, it outperforms temporal filtering while being significantly faster, and its memory overhead is one order of magnitude less than that of spectral filtering. Unlike spectral filtering, it is possible to extremely sparsify the dataset embeddings with only negligible drop in performance. This, together with its very low rank, makes our hybrid method even faster than spectral, despite being iterative. More importantly, while previous methods were long known in other fields before being applied to image retrieval, to our knowledge our hybrid method is novel and can apply e.g. to any field where graph signal processing applies and beyond. Our theoretical analysis shows exactly why our method works and quantifies its space-timeaccuracy trade-off using simple ideas from numerical linear algebra.

Fig. 1 .

 1 Fig. 1. (a) In descending order, eigenvalues of adjacency matrix W of Oxford5k dataset of n = 5, 063 images with global GeM features by ResNet101 and k = 50 neighbors per point (see section 6). Eigenvalues on the left of vertical red line at j = 300 are the largest 300 ones, candidate for removal. (b) Contour plot of upper bound φi(Lα(Wr)) of CG's relative error as a function of rank r and iteration i for α = 0.99, illustrating the space(r)-time(i) trade-off for constant relative error.

Fig.Fig.

  Fig. mAP vs. CG iteration i and mAP vs. time for temporal, spectral, and hybrid filtering. Sparsified hybrid is used with sparsity 99%.

Fig.Fig. 5 .

 5 Fig. mAP vs. level of sparsification on our hybrid method for r = 400. Dashed horizontal line indicates no sparsification.

Table 1 .

 1 Performance, memory and query time comparison on ROxford +R1M with GeM descriptors for temporal (20 iterations), truncated temporal (20 iterations, 75k images in shortlist), spectral (r = 5k), and hybrid ranking (r = 400, 5 iterations). Hybrid ranking is sparsified by setting the 99% smallest values to 0. Reported memory excludes the initial descriptors requiring 8.2 GB. U1 is stored with double precision.

	SHybrid

Table 2 .

 2 mAP comparison with existing methods in the literature on small and large scale datasets, using Medium and Hard setup of the revisited benchmark.

		ROxford	ROxf +R1M	RParis	RParis +R1M
		Medium Hard Medium Hard Medium Hard Medium Hard
	NN-search	64.7	38.5	45.2	19.9	77.2	56.3	52.3	24.7
	α-QE [27]	67.2	40.8	49.0	24.2	80.7	61.8	58.0	31.0
	Temporal [8]	69.9	40.4	61.6	33.2	88.9	78.5	85.0	71.6
	Temporal † [8]	-	-	59.0	31.4	-	-	79.7	65.2
	FSR [9]	70.4	42.0	22.6	5.8	88.6	77.9	66.6	40.2
	FSRw [9]	70.7	42.2	42.1	18.8	88.7	78.0	77.4	59.7
	SHybrid (ours)	70.5	40.3	62.6	34.4	88.7	78.1	82.0	66.8

https://landmarkscvprw18.github.io/

"Temporal" stems from conventional signal processing where signals are functions of "time"; while "spectral" is standardized also in graph signal processing.

https://github.com/ahmetius/diffusion-retrieval/

Acknowledgements: This work was supported by MSMT LL1303 ERC-CZ grant and the OP VVV funded project CZ.02.1.01/0.0/0.0/16 019/0000765 "Research Center for Informatics".

A General derivation

The derivation of our algorithm in section 4.1 applies only to the particular function (filter) hα [START_REF] Gordo | End-to-end learning of deep visual representations for image retrieval[END_REF]. Here, as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF], we generalize to a much larger class of functions, that is, any function h that has a series expansion

We begin with the same eigenvalue decomposition (6) of W and, assuming that h(W) converges absolutely, its corresponding decomposition

similar to [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF], where U1, U2 have the same orthogonality properties [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF]. Again, the first term is exactly the low-rank approximation that is used by spectral filtering, and the second is the approximation error

Again, we have used the series expansion [START_REF] Tong | Fast random walk with restart and its applications[END_REF] of h in ( 27) and (29). Now, equation ( 28) is due to the fact that (U2Λ2U

for i ≥ 1, as can be verified by induction, while for i = 0,

In both (30) and (31) we have used the orthogonality properties [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF]. Finally, combining ( 25), ( 29) and ( 6), we have proved the following.

Theorem 2. Assuming the series expansion (24) of transfer function h and the eigenvalue decomposition (6) of the symmetrically normalized adjacency matrix W, and given that h(W) converges absolutely, it is decomposed as

for n × n real symmetric matrix A. For h = hα and for

This general derivation explains where the general definition of function g (33) is coming from in [START_REF] Delvinioti | Image retrieval with reciprocal and shared nearest neighbors[END_REF] corresponding to our treatment of hα in section 4.1.


Abbas Krayani

Frédéric Dufour

Gilles Pijaudier-Cabot

Elastoplastic nonlocal damage model for concrete and size effect analysis

come

Elastoplastic nonlocal damage model for concrete and size effect analysis

INTRODUCTION

One of the specific behaviours of cementitious materials is that in uniaxial loading the traction strength is much smaller than compression strength. From that consideration, [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure[END_REF] has developed a damage model which takes only the positive principal strains into account. This criterion is widely used for its simplicity and relative accuracy to represent mode I failures of concrete structures under monotone loading especially.

For some specific structures (thick and/or prestressed) the stress state may be locally bi-or tri-axial. As for any geomaterials the confinement has a positive effect on the structural strength and depending on the level of confinement the failure mode may change to mode II. For this reason, the modified von Mises criterion [START_REF] De Vree | Comparison of nonlocal approaches in continuum damage mechanics[END_REF] seems to be more attracive as the second invariant of the stress tensor is introduced into the yield criterion. However, this criterion is symetric (traction-compression) and does not represent the reality.

The idea followed by [START_REF] Jason | An elastic plastic damage formulation for concrete: Application to elementary tests and comparison with an isotropic damage model[END_REF] is to couple plasticity, based on the effective stress, and local damage formulation in one single constitutive relation in order to represent a correct volumetric response of concrete, i.e. the confinement effect and the non symetric response of the material. The damage process, controlled by elastic strain , is described by the isotropic model developed by [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure[END_REF]. The plastic process is described by means of a yield sur-face, inspired from [START_REF] Etse | Fracture energy formulation for inelastic behavior of plain concrete[END_REF] and modified by [START_REF] Crouch | Application of a stress return algorithm for elasto-plastic hardening-softening models with high yeild surface curvature[END_REF]. In this approach, damage controls softening, while plasticity controls harding, in compression especially.

This model generates ill-posed mathematical problems due to the loss of ellipticity of the governing differential equations. Numerically, the results suffer from pathological sensitivity to the size and the orientation of the finite element mesh and the total energy dissipated by the fracture process tends to zero upon mesh refinement. Several nonlocal constitutive laws, which can be regarded as a remedy of the pathalogical mesh dependence, can be found in the literature. Whether they are in an integral form [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] or in a gradient form [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF], a salient characteristic of both types is the presence of a characteristic length in the constitutive relation.

The aim of this paper is to extend the previous approach to a non local model in order to simulate properly concrete structure behaviour with mixed mode failure. The regularization technique, based on the implicit gradient definitions of the nonlocal strain tensor, is adopted to overcome the mathematical problems due to the softening constitutive law (damage part), while the plastic part remains local.

In this contribution, the constitutive law of the model is first presented. Then, the validation and the performance of the model is evaluated in some structural cases.

MODEL FORMULATION

The effective stress approach has been chosen for coupling plasticity and damage effects. The total stress σ is then a function of the damage variable D and of an effective stress σ ′ that solves the equation of a plastic yield surface

σ = (1 -D)σ ′ (1)
2.1 Plasticity Elastoplastic constitutive law is defined by the following set of equations in which the effective stress σ ′ has been substituted to the applied stress σ :

ε σ ′ εp κ = ε e + ε p = Eε e = λm (σ ′ , κ) = λh (σ ′ , κ) (2)
where ε, ε e and ε p are the total, elastic and plastic strains respectively, E is the elastic stiffness tensor, m is the flow vector, κ is the set of internal variables and h is the plastic modulus. The plastic multiplier λ is given by the loading-unloading criterion (Kuhn-Tucker form)

F (σ ′ , κ) ≤ 0, λ ≥ 0, F (σ ′ , κ) λ = 0, (3)
where F (σ ′ , κ) is the yield function that defines the effective stress state. The plastic flow vector m is defined as (associated law)

m = ∂F ∂σ ′ (4)
The plastic yield surface, inspired from [START_REF] Etse | Fracture energy formulation for inelastic behavior of plain concrete[END_REF] and modified by [START_REF] Crouch | Application of a stress return algorithm for elasto-plastic hardening-softening models with high yeild surface curvature[END_REF], depends on four main functions ρ (effective stress invariant), k (hardening function), ρc (deviatoric parameter) and r (deviatoric shape function) (for more details, see [START_REF] Jason | An elastic plastic damage formulation for concrete: Application to elementary tests and comparison with an isotropic damage model[END_REF]))

F = ρ2 (σ ′ ) - k (σ ′ , k h ) ρc (σ ′ ) r 2 (σ ′ ) (5)
where k h is the hardening variable given as a function of confinement:

kh = √ 2 3 ε p :ε p ζ(σ ′ ) if k h < 1, kh = 0 if k h = 1, (6)
where ζ depends on the first normalized invariant

ζ = -A h + A 2 h -B h ξ + C h if ξ ≤ 0, ζ = -A h + A 2 h + C h if ξ > 0, (7)
A h ,B h ,C h are three model parameters. Note that equation ( 7) assumes that k h ranges between 0 and 1. When k h = 1, the yield surface becomes a limit surface with no hardening in pure elasto-plastic model. Combining equations ( 2) and (7) yields the following expression for the plastic modulus h:

h = √ 2 3 ∂F ∂σ ′ : ∂F ∂σ ′ ζ(σ ′ ) if k h ≤ 1, h = 0 if k h = 1, (8)
The plastic part of the constitutive relation contains 10 parameters in which four parameters remains at a fixed value: α = 0.5, γ = 0.99, A h = 7 × 10 -5 and k 0 = 0.1 (for more details, see [START_REF] Jason | An elastic plastic damage formulation for concrete: Application to elementary tests and comparison with an isotropic damage model[END_REF]).

Only two parameters r c and r t will be given in the applications where the others are taken as p = 0.4,

B h = 2 × 10 -2 , C h = 2 × 10 -6 and A = -0.5.

Damage model

Damage growth, controlled by elastic strain tensor, is governed by the following loading function:

g(ε e , k d ) = ε eq (ε e ) -k d (9)
where ε eq is the equivalent strain. k d initially equals the damage threshold ε D 0 and during the damage process equals the largest ever reached value of ε eq . The evolution of damage is governed by the standard loading/unloading condition:

g(ε e , k d ) ≤ 0, kd ≥ 0, kd g(ε e , k d ) = 0 (10)
For quasi-brittle materials such as concrete, the damage propagates much more easily under tension than under compression. Taking that into account, one possible choice is the so-called modified von Mises definition according to de Vree et al. ( 1995)

ε eq = k -1 2k (1 -2ν) I 1 (11) + 1 2k (k -1) 2 (1 -2ν) 2 I 2 1 + 6k (1 + ν) 2 J 2
where I 1 is the first invariant of strain tensor, J 2 is the second deviatoric invariant of the strain tensor and k is the ratio between uniaxial compressive and tensile strengths. This model uses only one single damage function

D = 1 - 1 -A ε eq + A exp (B (ε eq -ε D 0 )) (12)
where A andB are two constant of the model. In this contribution, Mazars' definition of equivalent strain is used

ε eq = 3 i=1 ε e i + 2 (13)
where ε e i + are the positive principal elastic strains. The damage variable D is determined as a linear combination of two damage variables D t and D c , that represent tensile damage and compressive damage respectively, by the help of two coefficients α t and α c computed from the elastic strain tensor and depending on the type of effective stress state:

D = α t D t + α c D c (14) D t,c = 1 - 1 -A t,c ε eq + A t,c exp (B t,c (ε eq -ε D 0 )) (15) α t,c = 3 i=1 ε t,c i ε e i + ε eq (16)
where A t , B t , A c , B c are four model parameters. ε t and ε c are the principal values of strains associated to the positive and negative effective stresses respectively.

The combination of plasticity and damage exhibits strain-softening and all the inherent difficulties attached to this specific material behavior, i.e., spurious strain localization and dependence of the energy dissipation on mesh refinement. To avoid this undesirable behaviour, the model may be regularized considering a nonlocal approach. In the present paper, a regularization technique, similar to that presented by [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF], is chosen. It is based on implicit gradient definition of the nonlocal strain tensor which is calculated for each component of the elastic strain tensor as:

ε e ij = εe ij -c∇ 2 εe ij (17)
where ∇ 2 is the Laplacian operator. The parameter c is of the dimension length squared and characterizes the non local interaction. This equation must be completed by the homogeneous Neumann boundary conditions

∇ε e .n = 0 (18)
with n the outward unit normal to the physical boundary. Now the local equivalent strain in equation ( 13) is replaced by its nonlocal counterpart:

ε eq = 3 i=1 ε e i + 2 (19) 3 STRUCTURAL APPLICATION 3.

Bar under uniaxial tension

The validation of the regularization technique introduced to the damage part of the model is illustrated by means of a tension bar with a defect in the middle. The bar length is 1 m and the cross section is 2 × 2 cm 2 . In simple tension tests, when the yield surface is reached, the effective response of the plastic damage constitutive law presented in this paper is perfectly plastic and no further damage can evolve due to a constant elastic strain. For this reason, the plastic process in this example is described by the Von Mises' model with linear hardening and no limit surface. The following parameters are used in this analysis: E = 3.3 × 10 10 Pa, ν = 0.2, A t = 0.9, B t = 10000, ε D 0 = 10.5 × 10 -5 , σ y = 3 × 10 6 Pa and E t = 1.1 × 10 10 Pa, where σ y the yield stress and E t the hardening modulus are the marerial properties for plastic part. Fig. 1 shows the damage distribution over the bar close to failure for a non local parameter c equal to 0.002 m 2 . The properties of the plastic gradient damage model are assessed by carrying out the analysis with several values of the gradient parameter c and a fixed mesh of 80 pentagon elements. Fig. 2 and Fig. 3 show the load-displacement response and the profile of damage for several values of the internal length parameter c. As desired, both the ductility in the loaddisplacement response and the width of the final damage profile increase with the internal length while for c= 0 the response of the regularized version coincide with that of the local one where the damage is localized in the center.

Shear failure mode: local approach

To demonstrate the abality and the performance of our constitutive law to predict correctly the failure mode in structural applications, the experimental behavior of reinforced concrete Continuous Deep Beam (reported by [START_REF] Asin | Behavior of statically interminate deep beams[END_REF]) has been simulated. The specimen, whith a minimum amount of stirups, has been selected for numerical simulation. However, the vertical shear reinforcement was not taken into account to simplify the study. The material parameters for the concrete are: E = 3.47 × 10 10 Pa, ν = 0.2, A t = 1, A c = 2.75, B t = 13000, B c = 2470, ε D 0 = 1 × 10 -4 , r c = 120 × 10 6 Pa and r t = 11 × 10 6 Pa. Moreover, a comparative numerical investigation has also been performed with two isotropic damage models in order to show advantages achieved by the introduction of plasticity into the damage formulation. The first one using the classical damage formulation defined by [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure[END_REF] with the following parameters: E = 3.47 × 10 10 Pa, ν = 0.2, The geometry and the loading setup are presented in Fig. 4 where the thickness is 0.15 m. Because of symmetry, only (the right) half of the specimen has been meshed with 3048 eight node cubic elements. With these geometry and parameters, the simulations are performed with the local approaches, since, using the nonlocal model presented in this paper would require a mesh small enough to take into account non local interaction and would reach computational lim-its due to the memory size of the calculation. A satisfactory agreement of the experimental and numerical results is observed. obtained by the damage models, where only flexural cracks are formed. Fig. 7 compares the load deflection curves of the plastic damage and damage models. As expected, the peaks for the damage models are underestimated as the failure is in mode I. This strong difference is due to the fact that for Mazars' criteria is written only as a function of positive elastic strain (Eq. ( 13)) which means that it is only due to microcracks in mode I, while for the modified von Mises definition, the strength under biaxial tensile-compressive stress for concrete is overestimated which leads to overestimation of the shear resistance of the material in this model (see e.g. [START_REF] Patzák | Adaptivity resolution of the localized damage in quasi-brittle materials[END_REF]).

A t = 1, A c = 2.75, B t = 8000, B c = 1768, ε D 0 = 1 × 10 -4 ,

Size effect in mixed mode tests

Quasibrittle materials like concrete are charactarized by gradual softening in a fracture process zone (FPZ) that is independent and not negligeable compared to the structure size D. Hence, the response of geometrically similar specimens is not proportional and yield a material size effect [START_REF] Bažant | Size effect in blunt fracture: concrete, rock, metal[END_REF]. It can be explained by stress redistribution and the associated energy release due to the developement of large cracks or a large FPZ prior to failure. A complete explanation of the size effect can be found in the textbook by [START_REF] Bažant | Fracture and Size Effect in Concrete and Other Quasi-Brittle Materials[END_REF]. From the modelling viewpoint, the size effect can be described only with failure models that contain internal length which is related to the characteristic length of the material, i.e. the width of the FPZ [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure[END_REF]. The aim of this section is to highlight the existence of the size effect on the mixed mode failure as it is predicted by the plastic gradient enhanced damage formulation presented in this paper. Asymmetric four-point bending tests proposed by [START_REF] Gálvez | Mixed mode fracture of concrete under proportional and nonproportional loading[END_REF] In order to avoid some possible mesh bias and ensure that the global energy dissipation in failure process is captured correctly, the size of an element in the process zone should be kept constant and small enough compared to the characteristic length. c n is a constant depending on geometry which plays no role in size effect analysis. The ultimate loads P u for the three sizes and its corresponding nominal strength are listed in the Table 1. The numerical results are interpreted with the help of Bazant's size effect law [START_REF] Bažant | Size effect in blunt fracture: concrete, rock, metal[END_REF].

σ N = Bf ′ t (1 + D/D 0 ) -1/2 (21)
where B is a dimentionless geometry-dependent parameter, D 0 is a characteristic size and f ′ t is the tensile strength of the material. Bf slope -1/2 where these two lines intersect at the abscissa D = D 0 . Obviously, it can be seen that there is a pronounced size effect that is in a good agreement with Bazant's size effect relation. The larger the specimen, the lower the nominal strength. It is also observed that the plasticity plays more important role for the smaller specimen (the hardening variable k h is more important as the specimen is smaller). On the other hand, for the larger specimen the developement of damage is dominant and its behavior is closer to linear elastic fracture mechanics and more different from a limit analysis solution.

CONCLUSION

In this work, the hardening plasticity is combined with a gradient enhanced damage formulation for modelling the structure behaviour with mixed mode failure. The example of tensile bar shows the effectivity of the regularization technique to overcome the mathematical problems due to the strain softening behaviour (damage part) where the damage distribution zone is controlled by a material heterogeneity parameter. Furthermore, the simulations of the Continuous Deep Beam show the improvements achieved by the coupled model compared to the damage models. The model predicts rather accurately the failure mode (e.g. mode II) and the ultimate load in complex structures.

Finally, the numerical simulations of three asymmetric four-point bending tests on notched specimens of different sizes show a dependence of the nominal strength on the specimen size. The larger the specimen, the lower the nominal strength. Thus, the nonlocal approach exhibits a size effect in mixed mode failure which is in a good agreement with Bazant's size effect relation.

Figure 1 :Figure 2 :

 12 Figure 1: Development of damage in a tensile bar with a defect in the middle.

Figure 3 :

 3 Figure 3: Damage profiles along the beam for different values of the internal length scale.

  and the other one employing the modified von Mises definition of the equivalent strain according to de Vree et al. (1995) with the following parameters: E = 3.47 × 10 10 Pa, ν = 0.2, A t = 0.9, B t = 8000 and ε D 0 = 1 × 10 -4 . The reinforcement bars are modelled with a von Mises plasticity model with E = 2.1 × 10 11 Pa, ν = 0.3, σ Y = 5 × 10 8 Pa (yield stress), Et = 4.5 × 10 10 Pa (tangent modulus).

Figure 4 :

 4 Figure 4: Geometry and loading setup of CDB.

  Fig. 5 provides a comparison of the experimental crack pattern with the numerical damage profile of the plastic damage model. As seen in the experiments, flexural damage band (mode I) starts first at midspan (stage A). Further load increase leads to the development of the flexural damage band and a propagation of damage (mode I & II) continues over the support (stage B), then a shear damage band (pure mode II) formed suddenly (stage C). Finally, failure occurs at the load introduction zone near the loading column (stage D).The ultimate failure load obtained experimentally was 1180 kN , the corresponding numerical value is 1285 kN , which is overestimated by 11%. Thus, the model can reproduce correctly the mode of failure and predict rather well the ultimate load.

Figure 5 :Figure 6 :

 56 Figure 5: Developement of crack pattern in CDB; Experiment versus simulation with plastic damage model.

Figure 7 :

 7 Figure 7: Load-deflection curves. Comparison between plastic damage model, Mazars' model and modified von Mises' model.

Figure 8 :

 8 Figure 8: Geometry and loading setup for asymmetic fourpoint bending tests.

Figure 9 :

 9 Figure 9: Profile of damage for small size specimen.

Fig. 10

 10 Fig. 10 shows the load deflection responses for the three different sizes. The nominal strength is obtained with the formula σ N = c n P u bD (20)

Figure 10 :

 10 Figure 10: Load-displacement curves for small size specimen (D=75 mm), medium size specimen (D=150 mm) and large size specimen (D=300 mm).

′tFigure 11 :

 11 Figure 11: Size effect on nominal strength.

Table 1 :

 1 Numerical results for three different sizes.

	Size	Depth Ultimate load Nominal strength
		D	P u	σ N
		(mm)	(kN)	(MPa)
	small	75	8.213	2.19
	medium	150	14.67	1.956
	large	300	24.57	1.638


Hussam El Cheikh
email: hussam-elcheikh@hotmail.fr

Bruno Courant

Jean-Yves Hascoët

Ronald Guillén

H El Cheikh

J.-Y Hascoët

Prediction and analytical description of the single laser track geometry in direct laser fabrication from process parameters and energy balance reasoning

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Since the invention of LASER by [START_REF] Maiman | Stimulated optical radiation in ruby[END_REF], laser applications growth rapidly and widely in many domains. A variety of industries (automotive, aerospace, navy, defense and many other sectors) use the laser technology for welding, cutting and hardening [START_REF] Courant | Surface treatment of titanium by laser irradiation to improve resistance to dry-sliding friction[END_REF]. Direct laser fabrication is a strategic application of laser technology based on laser cladding. A powder is injected in a laser beam and fused onto a substrate. The result of the deposition is a clad and the accumulation of those clads allows to fabricate a 3D object from a CAD file. This technique has several advantages over conventional fabrication techniques and can also be combined with subtractive processes as shown by [START_REF] Kerbrat | Manufacturability analysis to combine additive and subtractive processes[END_REF]. These advantages include faster processing speed, no requirement of tooling, ability to fabricate complex shapes and retention of a metastable microstructure/composition [START_REF] Majumdar | Laser processing of materials[END_REF]. Meanwhile, one of the disadvantages is the need of more mastery of the deposit laser track geometry [START_REF] Williams | Structural-Property Relationships in Materials -ME 6796[END_REF].

Several works and publications have been developed in the literature to improve and understand the DLF process and the theory of deposition. [START_REF] Wen | Modeling of coaxial powder flow for the laser direct deposition process[END_REF] present a numerical model which predicts the coaxial powder flow for the laser direct deposition process, including the particle stream flow in and after the nozzle and laser-particle interaction process. By solving the coupled momentum transfer equations between the particle and gas phase while incorporating particle temperature evolution, the dynamic and thermal behavior of multi-particles in the stream is completely modeled. [START_REF] Partes | Analytical model of the catchment efficiency in high speed laser cladding[END_REF] study a model which evaluates the catchment efficiency with respect to the molten pool geometry. Its model takes into account the melt pool geometry for low speeds, laser powers and feed rates, and by taking into account the particle melting during the flight for higher speed, laser powers and feed rates. [START_REF] Picasso | Finite element simulation of laser surface treatments including convection in the melt pool[END_REF] developed a 2D, stationary, finite element model to determine the shape of the melt pool in a known clad height by solving a stationary Stefan equation. Picasso et al. (1994) developed a simple geometrical model which allows the determination of the laser beam velocity and the powder feed rate from beam diameter, clad height and powder jet geometry for a specified laser power. [START_REF] Bamberger | Calculation of process parameters for laser alloying and cladding[END_REF] used the Mie theory in a simplified theoretical model to show the influence of the injected particles on laser beam energy transfer to treated surface. [START_REF] Toyserkani | 3-D finite element modeling of laser cladding by powder injection: effects of laser pulse shaping on the process[END_REF] use a transient three-dimensional finite element modeling, to show the effect of laser pulse shaping on the process. This model can predict the clad geometry as a function of time and process parameters including laser pulse shaping, travel velocity, laser pulse energy, powder jet geometry, and material properties. The decrease of the laser beam velocity or the increase of powder feed rate and laser pulse energy increases the clad height. [START_REF] Lalas | An analytical model of the laser clad geometry[END_REF] take into account the process speed and feed rate of the powder to predict the clad characteristics by an analytical approach based on the phenomenon of the surface tension considered as dominant. [START_REF] Onwubolu | Prediction of clad angle in laser cladding by powder using response surface methodology and scatter search[END_REF] investigate the prediction of the geometric form of the clad in laser cladding by powder using response surface methodology and scatter search optimization technique. [START_REF] De Oliveira | Analysis of coaxial laser cladding processing conditions[END_REF] present a theoretical and experimental study linking the process parameters to geometrical characteristics of an individual laser track using a numerical resolution. [START_REF] Lin | Clad profiles in edge welding using a coaxial powder filler nozzle[END_REF] suggest a simplified mode function to estimate the clad profiles based on the Gaussian distribution in the jet flow. [START_REF] Liu | Effects of powder concentration distribution on fabrication of thinwall parts in coaxial laser cladding[END_REF] presented a model to study the effects of process variables on wall thickness, powder primary efficiency and speed of forming a thin metallic wall in single-pass coaxial laser cladding. The static model of powder mass concentration distribution was defined as a Gaussian function. Either these models give incomplete information or they use the finite element method to calculate the laser track geometry and the induced thermal field and so are very time consuming because of the very small dimensions of the molten bath compare to the building entire part and the very short time of the interaction between one powder grain with the laser beam compare to the time to build the entire part. A way to improve the calculation quickness is to use an analytical description of the surface laser track. In this work, our goal is to propose a pragmatic way to optimize the process linking the laser track geometry to the process parameters, and also to obtain an analytical description of the laser track geometry to improve the thermal calculation quickness.

It is shown that it is possible to simulate the shape and the geometrical characteristics of the deposit clad, with an analytical formula, taking into account the process parameters, the width of the clad and the powder efficiency. First it is observed that the cross-section of the clad is always a part of a virtual disk. Then an expression is obtained linking the radius and the center of that disk to the clad width and the powder efficiency.

Materials and methods

A coaxial nozzle is used with an IPG Photonics © 700 W fiber laser (150 m). A carrier gas (Ar) flow, fixed at 3 L/min, imports the powder while a secondary gas (Ar) flow, fixed at 5 L/min, is shaping the powder stream. The laser beam is focused on the substrate surface while the distance of work between the nozzle and the focus plan is 5 mm. The beam analysis was carried out with the "FOCUS-MONITOR" from PRIMES Gmbh. The laser beam is scanned with a rotation pinhole (about 20 m in diameter). This pinhole takes a small part of radiation and 2 mirrors direct the signal to a detector. By scanning the laser beam layer by layer, the laser beam geometry and the spatial distribution of irradiance can be determined. This nozzle is also equipped with a cooling channel to dissipate heat. The focus plan of the powder is beneath the laser at 3.5 mm from the nozzle, the diameter of this focus plan is 0.6 mm.

In this study, powder of 316L stainless steel is used to fabricate the experimental clad. Particles size is between 45 and 75 m.

A low carbon steel substrate is used with the dimension mentioned in Fig. 1.

On each substrate, 6 laser tracks are made. In order to avoid any board effect, the 6 laser tracks are carried out at a distance of 10 mm from the substrate boards. The process duration is adjusted in order to obtain always the same length track of 9 cm.

To study effects of the main laser cladding parameters on the clad geometry, tracks are produced in a wide variety of power P (180, 280 and 360 W), velocity V (300, 600 and 900 mm/min) and mass feed rate Q m (0.025, 0.05 and 0.075 g/s). Samples are polished and etched with Nital (2%). Every combinations of the 3 values of the 3 parameters are explored which represents 27 single tracks. In order to evaluate the accuracy of the experimental results, two cross-sections are made and studied for each of the 27 laser tracks. And for one parameters combination four laser tracks have been made and analyzed. Cross-section characteristics change with the process parameters. Fig. 2 shows for each clad the height (H), the width (W), the depth penetration (H f ), the area (S) and the penetration area (S f ).

Theory

The powder is injected from a coaxial nozzle into the laser, melt and solidify on the substrate.

For a single laser track, the 3 important process parameters are the laser power, the powder feed rate, and the velocity of the laser advancement. These three parameters have a big influence on the clad dimensions.

Actually, we have experimentally observed that if dimensions of the clad depend on the process parameters, the cross-section always remain a part of a virtual disk to which the center can be below or above the work piece surface, so the model uses a circular form to describe the clad cross-section. The disk center is found above the work piece surface when the velocity is not great enough to avoid the accumulation of the deposited material. The usefulness of the proposed model is in one hand to predict the height, the width and the general form and surface of the cross-section in any cases and in the other hand to obtain an analytical description of the clad surface, which then, can be useful to simulate the thermal field induced by the laser heating.

This circular geometry suggests that surface tension forces are dominant. The calculation of the capillary length gives a result of 0.004 m. The liquid bath has a characteristic length smaller than this capillary length and so the gravity effect is negligible and the liquid shape is almost spherical (circular in the cross-section plane and slightly lengthened along the laser displacement axis because of the laser movement).

A similar approach exist, using the Simpson's rule, but is limited to cases where the clad is flattened. It allows describing the clad height and width by [START_REF] Davis | Optimum deposition parameters for the direct laser fabrication (DLF) of quasi-hollow structures[END_REF], and to optimize the clad parameters by [START_REF] Pinkerton | Modelling the geometry of a moving laser melt pool and deposition track via energy and mass balances[END_REF] and to avoid excessive dilution and overlap porosity by [START_REF] Colaco | A simple correlation between the geometry of laser cladding tracks and the process parameters. Laser processing: surface treatment and film deposition[END_REF].

Our approach allows the prediction of the form, the height, the width and the area of the clad from the process parameters.

Relationships between the geometrical characteristics R, W and the process parameters

Fig. 3 shows the shape of the virtual disk supposed to describe the laser track shape.

The disk is defined by the location of its center and its radius (R). They both depend on the process parameters. W is the clad width and also the chord of the disk at z = 0. Taking into account only the powder which is incorporated in the clad, the actual surface of the clad cross-section is given by P e Q m / V.

They are distnguished by the position of the circle center which can be located below or above the substrate surface. The powder mass deposited is related to the volume of the laser clad above the substrate surface. In the limit case the circle center is located on the substrate surface and one can show that P e Q m / V = R 2 /2 and then that the radius of the disk is given by:

R = 2P e Q m V
This case is a very interesting one to build a structure because of its 1/2 circular geometry.

Outside this limit case:

• Type 1: If P e Q m / V < R 2 /2, the center of the disk is below the substrate surface and at a distance R 2 -(W/2) 2 from it. The radius R is obtained by solving the equation:

P e Q m V = 1 6W R 1 -1 - W 2R 2 × 3R 2 1 -1 - W 2R 2 + 4W 2 (1)
This case is also interesting and corresponds to the flattened one.

• Type 2: If P e Q m / V > R 2 /2, the center of the disk is above the surface substrate and at a distance R 2 -(W/2) 2 from it. The radius R is obtained by solving the equation:

P e Q m V = R 2 - 1 6W R 1 -1 - W 2R 2 3R 2 1 -1 - W 2R 2 + 4W 2 (2)
This case cannot be used to build a full density structure but can be interesting for a porous one.

Finally the clad geometry can be expressed as an analytical expression which is the one of the circle of radius R and centered at the distance R 2 -(W/2) 2 below or above the substrate surface.

The radius of the disk can be calculated knowing the width (W) and the powder efficiency (P e ). It is to be noticed that the laser power is hidden behind P e . We suggest two different methods to predict P e and W.

Determining P e and W with parametric relationships

Each clad is defined by the corresponding process parameters (P, Q m , V). A parametric study is searched linking W and P e to the process parameters. A relationship as y = a(P ˛Q mV ) + b is thus established by El [START_REF] El Cheikh | Analysis and prediction of single laser tracks geometrical characteristics in coaxial laser cladding process[END_REF], where "y" is one of the laser track geometry parameters while a and b are constants. These relationships are very useful to predict the clad width and height but do not help to understand the physical phenomenon involved during the process. They are used in complement to our theoretical study (Table 1).

Determining W using a thermal study

The relationship between W and the experimental parameters does not involve the parameter Q m . This observation shows that a thermal process governs the parameter W. This is also in agreement with the usually accepted fact that only the powder which falls in the melting bath stays in the clad [START_REF] Lin | A simple model of powder catchment in coaxial laser cladding[END_REF]. The width of the clad Table 1 Relationships between the geometrical characteristics of the cross-sections and the process parameters [START_REF] El Cheikh | Analysis and prediction of single laser tracks geometrical characteristics in coaxial laser cladding process[END_REF].

Combined parameter R 2 H P 1/4 Q 3/4 m /V 0.944 W P 3/4 V 1/4 0.922 S PQm/V √ V 0.957 m PQm/V √ V 0.934 Pe P 3/4 Q -1/3 m V -1/2 0.910 S f Ln(P 4/5 Q -1/4 m ) 0.649 H f Ln(P 2 V 1/4 Q -1/4 m ) 0.765
is given by the molten pool dimension. The laser beam moves in the x direction at a constant velocity, many assumptions can be made to simplify the interaction between the laser and the substrate as follow:

1. The substrate is considered as a semi-infinite solid initially at the room temperature T 0 . 2. The molten zone is supposed to be a sphere whose surface is always at the molten temperature T m . 3. Thermo-physical properties are temperature independent.

In the steady state the pool geometry is constant. An energy balance is established between the absorbed laser energy in unit time, Q l , and the energy transported through the boundary of the molten pool in contact with the substrate and the deposit clad in unit time Q c by the heat conduction.

Q l = Q c (3)
The laser energy in unit time absorbed by the work piece is expressed as

Q l = Abs P (4)
In these conditions and in case of transient conduction, an analytical solution is given by [START_REF] Holman | Heat Transfer[END_REF]:

Q c = kA(T m -T 0 ) √ ˛t ( 5
)
where A is the surface at the temperature T m crossed by the heat flux, i.e. the area of the molten pool in contact with the substrate and the deposit clad.

In laser cladding, the melting zone is almost spherical and the surface contact A represents a part of its surface as shown in Fig. 4,

A = f d 2 .
The characteristic time for the passage of the molten zone above a point is given by [START_REF] Liu | Effects of powder concentration distribution on fabrication of thinwall parts in coaxial laser cladding[END_REF] as From Eqs. ( 5) and ( 6) the diameter of the molten zone is given by

d = Abs f 2/3 P √ ˛ k √ V (T m -T 0 ) 2/3 (7)
The width of the molten zone W is approximately equal to this diameter.

The width is given versus the molten temperature and the thermo physical proprieties of the clad material, and depends on the power and the velocity of the laser beam.

Results

Using the parametric relationships

The mathematical model presented above, allows the prediction of the cross-section geometry of the clad, knowing the width and the powder efficiency. The geometrical characteristics of the clad are measured with the aid of a software measuring tools integrated in the optical microscopy system. The width and the powder efficiency are linked to the laser power (P), the powder feed rate (Q m ) and the laser velocity (V) with the following equations:

W = 0.0030 P 3/4 V 3/4 -0.0108 (8)
P e = 0.312 Ln

P 3/4 V 1/2 Q 1/3 m -1.995 (9)
Figs. 5 and6 show the relationships between the combined process parameters and the experimental results.

We can now predict the geometrical form of the clad with the proposed mathematical model. This is presented in Fig. 7. The continuous line represents the prediction from the experimental measurements of the width and the powder efficiency. The dashed line is the prediction from the parameters (P, Q m , V) and the analytical relationships previously determined. It clearly shows that in the limits determined by our experiments we are able to predict the clad geometry from the analytical parametric relationships presented above.

In Fig. 8 the cross-section simulations obtained with the parametric relationships (dashed lines) are compared with those obtained with the experimental measured width (continuous line) and superimposed with the experimental photographs. These overlays are convincing.

Using the thermal analytical relationships

The thermal analytical model is now used to simulate the crosssection. Results obtained with the analytical model are given in Fig. 9.

The continuous line is the clad representation when W is calculated with the parametric relationships while the dashed line is the clad representation when the width W is calculated from the thermal analytical model.

Discussion

We showed in a previous study [START_REF] El Cheikh | Analysis and prediction of single laser tracks geometrical characteristics in coaxial laser cladding process[END_REF] that the laser clad geometrical characteristics are not due to the powder distribution in the powder jet stream but induced by the surface tension forces applied on the melting bath. It is experimentally established that the clad width depends only on the laser power and the laser velocity (Fig. 5). This shows the thermal dependency of the clad width and justifies the analytical thermal modeling of the laser clad width. The powder efficiency depends on the three process parameters through a logarithm function (Fig. 6).

Fig. 8 shows variation of the clad cross-sections with process parameters. Power remaining constant, the laser track height (H) and the cross-section area (S) increase if the velocity decreases or if the powder feed rate increases. For low velocities and high feed rate powder, the powder accumulation leads to the formation of a well defined cylindrical laser tracks. At the opposite for high velocity and low feed rate powder the clad profile is flattened and matches with a small part of a cylinder. When the laser power increases, the powder efficiency (Fig. 6) and the mass of incorporated powder per meter (g/m) increase.

In our case, the difference between calculated and experimental values always remains under 100 m (results obtained in cases 280 W, 0.025 g/s, whatever is the chosen speed and 180 W, 0.025 g/s, 600 mm/min) and in most cases under 50 m concerning the clad height and the width.

The two suggested ways linking the width to the process parameters are the parametric study or the balance energy reasoning and they lead to very narrow results (Fig. 9). This confirms the validity of the thermal model. Only a parametric study seems to be able to link the powder efficiency and the process parameters. However, the powder efficiency depends on the particular nozzle used and seems difficult to generalize with an analytical relationship. [START_REF] Partes | Analytical model of the catchment efficiency in high speed laser cladding[END_REF] proposed an analytical model of the catchment efficiency in which powder grains are incorporated if they fall in the melting bath and/or if they melt during the flight. The model depends on the melting bath geometry measured with a good thermal camera or difficult to simulate without the time consuming finite element method. And as authors say "toward higher melt pool lengths the difference between experiment and calculation of the catchment efficiency was strongly increasing". An analytical formula proposed by [START_REF] Liu | Effects of powder concentration distribution on fabrication of thinwall parts in coaxial laser cladding[END_REF] has been tried with our experimental results but this has not been successful. A simple and predictive analytical model for the powder efficiency catchment remains to be done. So P e is kept as an experimental entry in this model. The model proposed in this work, can predict the general geometry of the clad cross-section. This model is particularly interesting because it allows the width prediction of the clad layer analytically. The width depends on laser power, velocity of advancement, material proprieties as the conductivity, the diffusivity and the molten temperature. The constant K depends on the interaction between the laser beam and the material with the absorption coefficient and from the surface fraction of the melt pool through which the heat flows in the material. In our case, these good results are obtained with a constant K equal to 0.26. Independently of the used model the absorption coefficient has to be experimentally measured. This is particularly true concerning the direct laser fabrication process with coaxial powder injection because of the complex interaction process between the material and the laser beam. Few experimental trials have to be done in order to adjust the model whatever it is. These few experiments can be used to the K constant and the relationship between P e and the experimental parameters. This model can then be used to optimize the process parameters for one particular application. Choosing the width W and the radius R, one can find the process parameters P, Q m and V. The validity domain of this model is limited by the chosen extreme values of the experimental parameters P, Q m and V. One can choose a larger domain than our but must find a good relationship between P e and the process parameters P, Q m and V. The thermal model assumptions are not restrictive.

W = K P √ ˛ k √ V (T m -T 0 ) 2/3

Conclusion

Laser cladding process success depends on the ability to predict the laser clad geometry (form and dimensions) from the process parameters. Acting the fact that the laser clad cross-section looks like a disk we have related the circle center position and the circle diameter to the clad width W and the powder efficiency P e . Two ways are explored to calculate W from the process parameters. Finally the cross-section simulations are compared to the experimental photographs. These overlays are convincing and so this modeling is able to predict the complete clad geometry with an analytical expression.

Assuming the fact that the clad geometry is governed by the surface tension on the melted bath and so that the clad geometry is a part of a disk it is shown that it is possible to predict, with a good uncertainty, the complete geometrical characteristics of laser clad. The other interest point of this model is that this analytical formula can help in a thermal study to determine the shape of the molten pool, and the distribution of temperature in the piece.

Fig. 1 .

 1 Fig. 1. Substrate dimensions and positions of single laser tracks on the steel substrate.

Fig. 2 .

 2 Fig. 2. Cross-section measures for each single laser track.

Fig. 3 .

 3 Fig. 3. Schematic of the cross-section.

Fig. 4 .

 4 Fig. 4. Schematic view of the spherical molten pool.

Fig. 5 .

 5 Fig. 5. Width via the combined parameters.

Fig. 6 .

 6 Fig. 6. Power efficiency via the combined process parameters.

Fig. 7 .

 7 Fig. 7. Confrontation between the geometrical representation of the clad geometry obtained with experimental measures of W and Pe (continuous line) and analytical determination of W and Pe from the parametric relationships (dash line) at P = 180 W.

Fig. 8 .

 8 Fig. 8. Confrontation between experimental results and simulated ones obtained with the analytical relationships at different power: (a) 180 W, (b) 280 W and (c) 360 W.

Fig. 9 .

 9 Fig. 9. Confrontation between the geometrical representation of the clad geometry obtained with the analytical relationships (continuous line) and with the thermal analytical model (dash line) at different power: (a) 180 W, (b) 280 W and (c) 360 W.