J M Cadou

L Duigou

N Damil
email: n.damil@univh2m.ac.ma

M Potier-Ferry
email: michel.potierferry@univ-metz.fr

Convergence acceleration of iterative algorithms. Applications to thin shell analysis and Navier-Stokes equations

Keywords: Convergence acceleration, iterative correctors, polynomial extrapolation, Padé approximants, nonlinear problems. PACS

 and an original algorithm for vibration of viscoelastic structures. We first describe the iterative algorithms for the considered nonlinear problems. Secondly, the two accelerating techniques are presented. Finally, through several numerical tests emanating from the thin shell theory, Navier-Stokes equations and vibration of viscoelastic shells, permit to show the advantages and drawbacks of each accelerating technique.

Introduction

The computation of solutions of nonlinear problems requires the use of a prediction and a correction step. The correction step is usually done by iterative methods, such as Newton type methods. When solving problems with a large number of unknowns, this correction phase can consequently increase the CPU time due to the treatment of the matrices (triangulation, backward and forward substitution).

For this reason, modified versions of the Newton method can be used to avoid large amounts of CPU time. The most common one, that will be called "modified Newton" in this paper, needs only one matrix triangulation that is computed at the first iteration. Nevertheless, the use of such a method can increase the number of iterations and sometimes leads to divergence.

In recent works [START_REF] Damil | An iterative method based upon Padé approximants[END_REF], [START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF] and [START_REF] Cadou | Projection technique to improve high order iterative correctors[END_REF], high-order iterative algorithms have been proposed. From these studies, one can distinguish two classes of methods: on the one hand, a high-order Newton algorithm that is efficient, robust and often converges after a single iteration, and on the other hand, the so-called L * algorithm which does not require any triangulation of matrix. Like Newton's modified method, the drawback of the L * -algorithm is the number of iterations and the possibility of divergence.

In this work, we propose associating these kinds of iterative correctors with some techniques to accelerate the convergence of sequences of vectors (the iterates obtained with the corrector). Only two methods are considered here. The first method is a polynomial extrapolation method and its derivatives: minimal polynomial extrapolation (MPE [START_REF] Cabay | A polynomial extrapolation method for finding limits and antilimits of vector sequences[END_REF]), the reduced rank extrapolation (RRE [START_REF] Mesina | Convergence acceleration for the iterative solution of the equation X = AX + f[END_REF]) and the modified minimal polynomialextrapolation method (MMPE [START_REF] Brezinski | Généralisations de la transformation de Shanks, de la table de Padé et de l'ǫ-algorithme[END_REF]. Applications and numerical comparison of these methods can be found in [START_REF] Jbilou | Some results about vector extrapolation methods and related fixed point iterations[END_REF][START_REF] Jbilou | Vector extrapolation methods, Applications and numerical comparison[END_REF]). Due to the results and the conclusions presented in [START_REF] Jbilou | Vector extrapolation methods, Applications and numerical comparison[END_REF], only the MMPE method is applied in this work. Moreover, a sizeable bibliographic review related to these methods can be found in [START_REF] Graves-Morris | The epsilon algorithm and related topics[END_REF][START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF]. Nevertheless, few works deal with applications of these methods in solid or fluid mechanics. The second method uses vectorial Padé approximants [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF] and has been applied in numerical solid mechanics. Despite the extensive bibliography on convergence acceleration, there are few applications of these methods to practical problems of nonlinear mechanics. We believe that significant computational savings can be obtained by coupling some modified Newton and convergence acceleration techniques. The aim of the present paper is to evaluate these opportunities and to establish computational procedures that are reliable and suitable for problems of structural mechanics.

These acceleration techniques will be evaluated on the basis of numerical tests emanating from nonlinear elasticity, vibrations of viscoelastic shells, and Navier-Stokes equations. In the first part of this paper, the convergence accelerating techniques are coupled with the modified Newton method. Then, these techniques are associated with high-order iterative correctors and numerically evaluated in the buckling of thin shell and in a fluid mechanics framework. Finally, an original iterative algorithm for the vibrations of damped structures, presented in this work, is coupled with the MMPE algorithm. The performances of this algorithm is compared to the results presented in [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF].

Iterative algorithms and convergence acceleration

In this section, two iterative correctors are recalled. An original algorithm for the vibration of viscoelastic shells is also presented. The common point of these techniques is the use of few matrices in the iterative scheme. Therefore, a great number of iterations can be necessary to get the desired accuracy. Convergence accelerating techniques are then used to decrease the number of iterations and improve the quality of the solutions obtained.

Some iterative algorithms for nonlinear problems

The equations of the considered nonlinear problem, for the buckling of nonlinear elastic thin shell or for fluid mechanics, can be written in a quadratic framework such as:

R(U, λ) = L(U ) + Q(U, U ) -λF = 0 with U ∈ IR n , F ∈ IR n and λ ∈ IR (1)
where U is the unknown vector (U includes the displacement and the stress in nonlinear elasticity, the velocity and the pressure in fluid mechanics), F is a given vector and λ is a load parameter. L(.) and Q(., .) represent linear and quadratic operators. In this work we consider two iterative correctors to improve the quality of a trial solution U 0 . The first one is the classical modified Newton algorithm and the second one is issued from recent works ( [START_REF] Damil | An iterative method based upon Padé approximants[END_REF][START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF]).

In classical iterative algorithms, the approximated solution after k iterations is:

   U (k) = U 0 + ∆U (k) λ(k) = λ 0 + ∆λ(k) (2)
where the correction (∆U (k), ∆λ(k)) is defined by the following sequence:

   ∆U (k) = ∆U 1 + ∆U 2 + • • • + ∆U k ∆λ(k) = ∆λ 1 + ∆λ 2 + • • • + ∆λ k (3)
One objective of the present work is to accelerate sequences such as (3) by using convergence accelerating techniques.

By introducing the relations (2) into the nonlinear problem (1), we obtain the following nonlinear problem to be solve at the iteration k:

R(U (k), λ(k)) = R k-1 +L(∆U k )+2Q(U (k -1), ∆U k )+Q(∆U k , ∆U k )-∆λ k F (4
) where R k-1 is the residual vector computed at the approximated solution U (k -1) and is defined by

R k-1 = R(U (k-1), λ(k-1)) = L(U (k-1))+Q(U (k-1), U (k-1))-λ(k-1)F
(5) Within iterative methods, such as Newton-Raphson, the quadratic term of the equation ( 4) is neglected and we have to solve the following linearized problem:

L k-1 t (∆U k ) -∆λ k F = -R k-1 (6) where L k-1 t
is the tangent matrix computed at the approximated solution U (k -1):

L k-1 t (∆U k ) = L(∆U k ) + Q(U (k -1), ∆U k ) + Q(∆U k , U (k -1)) (7)
When a modified Newton algorithm is considered, the previous tangent matrix is computed once and is constant during all the iterations. The vector ∆U k and the scalar ∆λ k are the solutions to the following linear equation:

L 0 t (∆U k ) -∆λ k F = -R k-1 (8)
where L 0 t is the tangent matrix computed at the trial solution U 0 . In Newton [START_REF] Mesina | Convergence acceleration for the iterative solution of the equation X = AX + f[END_REF] or modified Newton (8) algorithm, there are n+1 unknowns (∆U k , ∆λ k ) for n equations. Thus, a supplementary equation is needed to have a wellformulated problem. In this study, we consider two cases. In the first one, the load parameter is fixed during all the iterations (i.e. ∆λ k = 0, a load control). In the second one, a sort of arc-length control is used. In the latter case, the most natural way to specify the variation of λ during the iterations is to require the correction (∆U k , ∆λ k ) to be orthogonal to the slope, denoted (U ′ , λ ′ ), of the response curve. The supplementary condition associated to the problems ( 6) and ( 8) is:

< ∆U k , U ′ > +∆λ k λ ′ = 0 ( 9
)
where < •, • > represents the euclidian scalar product. In this study the slope (U ′ , λ ′ ) is known. Indeed, as the prediction curve (the trial point (U 0 , λ 0 )) is determined by using the Asymptotic Numerical Method [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF] which gives analytical solutions, the slope (U ′ , λ ′ ) is easily computed (for details see references [START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF][START_REF] Cadou | Projection technique to improve high order iterative correctors[END_REF]).

The second iterative corrector, used in this work, has been recently presented in references [START_REF] Damil | An iterative method based upon Padé approximants[END_REF][START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF]. In these works high-order correctors have been proposed which use the homotopy technique, perturbation method and Padé approximants.

To define these correctors, we consider the nonlinear problem (4) and we modify it in the following manner:

[L * + η(L k-1 t -L * )](∆U k ) + Q(∆U k , ∆U k ) + ηR k-1 = 0. ( 10
)
where η is a real parameter (0 ≤ η ≤ 1), L k-1 t is the tangent operator at the iteration k -1.

The solution ∆U k of the modified problem [START_REF] Boer | Padé approximant applied to a non-linear finite element solution strategy[END_REF] continuously passes from 0 for η = 0 to the solution of the initial nonlinear problem for η = 1. To compute this solution, a perturbation method is used and ∆U k is sought in the form of a truncated integro-power series:

∆U k = η∆U k 1 + η 2 ∆U k 2 + ... + η n ∆U k n . (11)
This asymptotic expansion is introduced in problem (10) and by equating like power of η we obtain a set of 'n' linear problems which all have the same operator to be inverted (L * ). Finally when all the terms ∆U k n are computed, the polynomial approximation [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF] is then improved by using Padé approximants [START_REF] Brezinski | Padé approximants, Handbook of Numerical Analysis[END_REF][START_REF] Baker | Padé approximants, Encyclopedia of Mathematics and its Applications[END_REF][START_REF] Boer | Padé approximant applied to a non-linear finite element solution strategy[END_REF][START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF]. The matrix to be inverted, denoted by L * , is priori chosen. In practice, the choice will be a matrix that has been triangulated in a previous step of the calculation or a matrix which is easily inverted. In other words, L * is a preconditioner. This L * -algorithm generally requires more than one iteration to obtain the required accuracy. Indeed, sometimes the number of iterations needed for the required accuracy can be high (k =10, 15, 20 iterations). This L * algorithm will be associated with some convergence accelerating techniques. One can note that, in this case, two successive accelerating methods are used. The first one concerns the series computed with expression [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF] and the second one concerns sequence (3).

The third nonlinear iterative algorithm, studied in this paper, deals with the vibration of viscoelastic shells. The unknowns, the displacement u and the vibration frequency ω, are thus solutions of a complex eigenvalue problem [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF]:

R(U ) =      [K(ω) -ω 2 M ] u t u M .M.(u -u M )      = 0 with u ∈ I C n and ω ∈ I C ( 12
)
where M is the real mass matrix, K(ω) is the complex stiffness matrix depending on frequency. The second equation of ( 12) is a normalization condition and u M is the real mode of the undamped structure, ω M is the corresponding real eigenfrequency. In this work, we propose iteratively solving problem [START_REF] Jbilou | Some results about vector extrapolation methods and related fixed point iterations[END_REF] which is nonlinear, with respect to the frequency. If one applies the Newton iterative method to [START_REF] Jbilou | Some results about vector extrapolation methods and related fixed point iterations[END_REF], the tangent matrix is given by:

   K(ω) -ω 2 M [K ′ (ω) -2ωM ]u M t u M .M 0    (13)
This tangent matrix is complex-valued, because K(ω) and ω are complex. To avoid the triangulation of a complex matrix, we propose replacing this tangent matrix by a real one. We have chosen the preconditioner introduced in [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF] within the framework of a high-order algorithm similar to (10) [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF].

Hence, at the iterate k, the problem to be solved thus becomes:

K * ∆U k =    K(0) -ω 2 M M -2ω M M u M -2ω M M u M 0         ∆u k ω k      = R(∆U (k -1)) ( 14
)
where K(0) is the elastic matrix which is symmetric and real [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF]. Therefore, with the previous equation, we have defined an iterative algorithm based on a real operator K * . At the iterate k, a sequence identical to the one obtained with the modified Newton algorithm (eq 3) is built . The convergence of this sequence is then improved by using convergence accelerating methods. Unlike a Newton algorithm, this procedure requires the inversion of one real matrix for all the iterations instead of one complex matrix per iteration.

Two convergence acceleration techniques

In this section, a brief description of two algorithms to accelerate the convergence of sequences is given: the Modified Minimal Polynomial Extrapolation and the vectorial Padé approximants. These methods will be applied to sequences obtained either by the modified Newton algorithm, the high-order iterative algorithm or the original iterative corrector defined for the vibration of viscoelastic shells.

With vector extrapolation methods or Vectorial Padé Approximants (VPA) [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF], sequence (3) is transformed into a new one with a faster convergence than the initial sequence. The initial sequence (3) is then replaced by:

∆U new (k) = k-1 i=1 e i ∆U i (15)
where the subscript "new" designates the new vector obtained with either a vector extrapolation method or with VPA. In the first case, e i are scalar coefficients and in the second one, e i are rational functions. First, the vector extrapolation method is described.

After 'k' iterations of an iterative method (such as the ones defined in the previous section), providing 'k' vectors ∆U k , one can build an approximated vector, ∆U MMPE (k), by using the following MMPE algorithm:

Step

1 : Initialization choose 'k-1' vectors Y i Step 2 : Computations M ij = (∆U j+1 -∆U j ).Y i (16)
f i = -∆U 1 .Y i with i, j = 1, .., k -1 (17) Step 3 : Solve the linear system [M ]{a} = {f } (18)
Step 4 : Build the extrapolated vector

∆U MMPE (k) = k-1 i=1 a i ∆U i (19)
In this study, in which only consider polynomial methods [START_REF] Jbilou | Vector extrapolation methods, Applications and numerical comparison[END_REF] are considered, coefficients a i are obtained from an orthogonality relation and are then solutions of the linear system [START_REF] Giraud | The loss of orthogonality in the Gram-Schmidt orthogonalization process[END_REF]. Within this extrapolation, the chosen vectors, Y i are arbitrary linearly independent ones. Thus, in the numerical tests presented in this work, the 'n-1' vectorsY i are either the residual R(∆U (i)) or the vectors ∆U i . Finally, the transformation [START_REF] Chen | Eigensolution reanalysis of modified structures using epsilon-algorithm[END_REF] has been designed to accelerate the convergence of the sequence (3), where coefficients a i are computed from system [START_REF] Giraud | The loss of orthogonality in the Gram-Schmidt orthogonalization process[END_REF].

Concerning the Vectorial Padé approximants [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF], the sequence (3) ('k' vectors ∆U k issued from an iterative scheme) is replaced by its equivalent rational representation, ∆U Padé (k), by using the following algorithm:

Step 1: Gram-Schmidt orthonormalisation

∆U n = n j=1 α nj V * j with n = 1, .., k (20)
Step 2: Compute the coefficients d i

d i = - 1 α k-i,k-i (α k,k-i + i-1 n=1 d n α k-n,k-i ) with i = 1, .., k -1 (21)
Step 3: Compute the coefficients c n

c n = k-1-n j=0 d j k-1 j=0 d j with d 0 = 1 ( 22
)
Step 4: Build the Vectorial Padé Approximant

∆U Padé (k) = k-1 n=1 c n ∆U n (23)
Let us remark that coefficients α ij are deduced from a Gram-Schmidt orthonormalisation, the vectors V * j in (20) being orthonormalized. In this work, a modified Gram-Schmidt orthonormalisation is used because it is more stable than the classical Gram-Schmidt procedure. This algorithm is summarized below:

for i = 1, • • • , n q i = ∆U i forj = 1, • • • , i -1 α ij =< q i , V * j > q i = q i -α ij V * j enddo V * j = q i ||q i || enddo
These two acceleration techniques seem quite different, but they belong to the same family of algorithms. Futhermore it has been proved [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF] that, with a proper choice of the projection vectors Y i , they become mathematically identical. With the MMPE algorithm, the orthonormalization procedure could be avoided and the possibility of choosing vectors Y i leads to more freedom to optimize the procedure.

A few comments and remarks

In this section, we discuss how the previous algorithms are managed. Firstly, within the MMPE algorithm, we have to choose vectors Y i which are needed to compute coefficients a i in the expression [START_REF] Chen | Eigensolution reanalysis of modified structures using epsilon-algorithm[END_REF]. In the literature, only the fact that these vectors are arbitrary and linearly independent is mentionned. There has been no discussion, to our knowledge, about the choice of these vectors. Therefore in the presented iterative methods, these vectors are either the residual vectors or the solution vectors ∆U i . In the following section, we discuss this choice and the consequence of it on the final result.

Secondly, another important point, especially for MMPE, is to choose the number of vectors in the sequence to be accelerated, i.e. the number k in [START_REF] Chen | Eigensolution reanalysis of modified structures using epsilon-algorithm[END_REF]. If the solution obtained after acceleration is not sufficiently accurate, the whole process will be restarted: actualisation of the initial point, computation of at most k vectors, acceleration and so on.

Thirdly, the Vectorial Padé Approximants technique involves an orthonormalisation of the sequence. There are several algorithms to achieve this orthonormalisation and this can influence the efficiency of the procedure [START_REF] Jamai | Influence of iterated GramSchmidt orthonormalization in the asymptotic numerical method[END_REF][START_REF] Giraud | The loss of orthogonality in the Gram-Schmidt orthogonalization process[END_REF]. Hereafter, we shall use the so-called modified Gram-Schmidt algorithm. It is possible to choose orthonormalized vectors Y i within MMPE. Of course, the orthonormalisation can be expensive which is why we propose considering a modified scalar product:

V.U = nred i=1 V i U i ( 24
)
with nred lower than the dimension n of the vectors.

All these points will be analysed in the next sections.

3 Numerical applications with a one-matrix algorithm

Two benchmarks from nonlinear structural mechanics

In this section, we apply the previous methods to some traditional numerical tests in thin elastic shell analysis. The objective is to define an optimal method for using accelerating methods in the finite element framework. In the second part of this section we compare the performances of the high-order correctors when associated with accelerating methods.

The first numerical test considered in this work is the classical geometrically nonlinear cylindrical roof loaded by a single force [START_REF] Eriksson | Derivatives of tangential stiffness matrices for equilibrium path descriptions[END_REF][START_REF] Kouhia | Some aspects of efficient path-following[END_REF]. The geometric and material characteristics are given in Figure (1(a)). Only one quarter of the roof is modelled by 200 triangular DKT elements (726 d.o.f.). On the nonlinear curve obtained by ANM (method of reference [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF], Padé approximants at order 25, see Figure 1(b)) we have chosen three trial points. The corresponding displacement, load parameter and norm of the residual of these three trial points are given in Table [START_REF] Wynn | On a device for calculating the e m (S n ) transformations[END_REF]. One can note that the residual norm increases from point 1 to point 3. Only one octant of the cylinder is modelled by 900 quadrilateral elements with 8 nodes per element. The corresponding number of d.o.f. is roughly 17000. This example has been chosen because the response curve is quite difficult to obtain due to turning or bifurcation points. The correction method has to be very efficient to converge to the desired solution. In the computation, the prediction step is carried out with the help of the ANM and the correction is realized with the Newton-modified method coupled with the MMPE method. In the As the computational strategy is a key-point within nonlinear problems, due to turning or bifurcation points, all the computations use an arclength control (except for point 3 of the first example where the considered iterative methods converge only with a fixed load parameter) in the correction scheme.

A first comparison

In this part, we compare the performances of the two convergence accelerating methods presented in the previous section. This comparison is achieved by considering the Newton modified method (8) and the two previous examples.

In Figure (3) the decimal logarithm of the residual norm is plotted versus the iteration number for the three trial points (point 1 in Figure 3(a) and so on) for the example of the cylindrical thick roof. For each trial point, three iterative correctors are evaluated: the classical Newton modified algorithm, the Newton modified algorithm associated with MMPE and finally the Newton modified algorithm associated with VPA. In the case of the MMPE accelerating convergence method, vectors Y i of expression [START_REF] Giraud | The loss of orthogonality in the Gram-Schmidt orthogonalization process[END_REF] are the residual vectors computed at the previous iterates. These Figures show several interesting features. Firstly, a better convergence is obtained when an accelerating convergence method is associated to the modified Newton corrector. For instance, for trial point 1 (Figure 3(a)), if the desired accuracy is fixed and equal to 10 -8 (dashed line in Figure 3(a)), the convergence is reached with 10 iterations without acceleration and 7 iterations with MMPE and VPA. For - MMPE with Y i = R i . In parenthesis the number of iterations to get a residual lower than 10 -8 with a 'pure' Newton corrector and the logarithm of the residual vector after convergence of the Newton algorithm.

trial point 2, the number of iterations to get the required accuracy is 28, 11, 10 respectively without acceleration, with MMPE and VPA. As for the third trial point (Figure 3(c)), the results obtained with the accelerating convergence techniques are very interesting. Indeed, without acceleration, the norm of the residual obtained after 30 iterations is close to 10 -1 . Whereas, when the VPA is used, the required accuracy (10 -8 ) is reached with 23 iterations. With MMPE, the residual norm after 30 iterations is 10 -5 . Similar results are obtained when considering the second example. These results are given in Table 2 Pull-out of an open cylinder. Comparison of the number of iterations to get the desired accuracy 10 -5 . Arc-length scheme.

Figure ( 4) for point 1 (see Figure 2) and in Table (2) for the two other points.

- The second interesting feature concerns the MMPE or the VPA algorithms. One can remark in Figures [START_REF] Graves-Morris | The epsilon algorithm and related topics[END_REF] and (4) that when these two accelerating techniques are used with a low number of iterations (lower than 10), the evolution of the residual versus the number of iterations are nearly the same for the two techniques. In a recent work [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF], the authors have established the mathematical equivalence between these two techniques when a specific choice of vectors Y i is made but some differences are observed in the numerical practice, especially when a large number of vectors is used.

The third feature concerns the MMPE techniques. In Figure [START_REF] Graves-Morris | The epsilon algorithm and related topics[END_REF], one can see that when the number of iterations is high the MMPE technique becomes unstable. Indeed, in this case, an increase of the number of vectors leads to a divergence of the method, see Figure (3(b)). In the following section, we discuss this particularity and try to define the best way to use the MMPE algorithm.

What is the best way to use MMPE?

Several computational strategies have been tested to avoid these convergence problems with MMPE. The first point is the choice of the projection vectors Y i . Two strategies have been evaluated: projection on the residual (Y i = R i ) or on the increments (Y i = ∆U i ). The results in Figure 5(a) show that there are no significant differences and the process does not converge if the number of vectors is greater than 10. From Figure 5(b), it appears that the processes converge very well if these vectors Y i = R i or Y i = ∆U i are orthonormalized. Hence, orthonormalization seems to help avoid the numerical instabilities and the convergence is about the same as with Padé approximants. In a recent work [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF], the authors have proved the mathematical equivalence between VPA and MMPE, with the choice of Y i . Another efficient way of using the MMPE technique is to restart the accelerating technique after k iterations, as explained in paragraph 2.3. This technique is evaluated in Figure 6, for trial point 3, where the parameter k is equal to 5, 10 or 15. The results obtained with this technique are compared to the ones obtained with the Padé approximants (considered here as the reference) and with the MMPE algorithm without restart (with k = ∞). Let us remark that these results are obtained with non-orthonormalized vectors Y i . These curves show a better convergence when the number of iterations is fixed. Indeed, when the residual criterion is equal to 10 the best results are obtained when the parameter k is low. This emphasizes once more that MMPE is a very efficient algorithm, but that it can become unstable when the number of vectors is high. In fact, these results show that when associated with a restart procedure, the MMPE algorithm is efficient and robust, even with a high number of iterations. For instance, let us assume that 23 iterations are needed: one must first compute k = 10 iterates and apply MMPE to these ten vectors, then a second sequence of 10 vectors should be computed and accelerated in the same way and finally, a sequence of 3 vectors.

About computation time

Our last topic of discussion is the amount of CPU time needed by the MMPE algorithm. Here, two examples are considered: the pull-out of an open cylinder and the cylindrical thick roof but with a greater number of unknowns: 39.366 d.o.f.. For the pull-out of an open cylinder, we take the convergence results presented in Figure [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF]. In this case, a pure Newton method needs 4 iterations to reach the required accuracy whereas 8 iterations are needed when using a modified Newton method coupled with the MMPE technique. The corresponding CPU time for each method is given in Table [START_REF] Graves-Morris | The epsilon algorithm and related topics[END_REF]. The CPU time given for the pure Newton method corresponds to the 4 iterations and more precisely: 4 tangent matrices to build, 4 matrix triangulations and 4 forward and back- ward substitutions. For the modified Newton method, the CPU time given in Table (3) corresponds to 8 forward and backward substitutions and one tangent matrix construction and triangulation (for the first iteration). As for the MMPE method, the indicated CPU time is the total amount of CPU time for the 8 iterations of the modified Newton method. The results presented in this Table show that the CPU cost of one iteration of the modified Newton method is lower than 5% of the CPU time required for an iteration of a pure Newton scheme. Moreover, the CPU time of the MMPE step (the convergence acceleration of the 8 iterations) is nearly equal to 10% of a single modified Newton iteration and 1% of the 8 modified Newton iterations. Hence, the CPU time for the MMPE method is insignificant compared to what is needed for one iteration of a pure Newton algorithm. For problems with a large number of unknowns, the amount of CPU time for the MMPE could increase due to the orthonormalisation step in this algorithm. Therefore, to avoid the increase of the CPU time we propose to orthonormalize the sequence with respect to a smaller number of degrees of freedom (see nred in equation 24 or reference [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF]). One can remark in Table (3) that for the considered problem, the nred parameter has no influence on the total amount of CPU time. Nevertheless, one has to evaluate the influence of the nred parameter on the accuracy of the convergence acceleration method.

Thus, in Figure [START_REF] Damil | An iterative method based upon Padé approximants[END_REF], we have plotted the evolution of the residual vector versus the iteration number for several values of nred. Hence, even if nred is small (for example 40 d.o.f.), MMPE is efficient except if the number of iterations is too large. For a larger number of vectors (k > 10), it seems that the instabilities described in the previous section happen earlier, if nred = 40. Thus, a good way to apply MMPE is to restart after a few iterations (k≤10) and to consider orthonormalized vectors with a modified scalar product (nred being not too large).

A third acceleration convergence method has been tested, the so-called ǫalgorithm [START_REF] Wynn | On a device for calculating the e m (S n ) transformations[END_REF][START_REF] Wynn | On the convergence and stability of the epsilon algorithm[END_REF][START_REF] Graves-Morris | The epsilon algorithm and related topics[END_REF]. The convergence acceleration obtained with this ǫ-algorithm is nearly the same as with the other two techniques. Nevertheless, the amount of CPU time needed to build an accelerated sequence with the ǫ-algorithm is very large. Consequently, the latter is not an attractive method for a large scale problem. This very substantial amount of CPU time has already been mentioned in reference [START_REF] Jbilou | Vector extrapolation methods, Applications and numerical comparison[END_REF] and is due to the construction of a 'ǫ-table' at each iteration. The construction of this 'ǫ-table' needs the computation of a great number of vectors, the ǫ-vectors, which leads to a large amount of CPU time.

A solution to this drawback is the storage at each iteration of the already computed ǫ-vectors. Unfortunately, this solution cannot be used in a finite element framework because of the large number of unknowns.

-

Numerical applications with matrix less algorithms

The numerical results presented in the previous section show that the MMPErestart algorithm seems to be the most efficient accelerating technique. In this section, the MMPE technique will be used to accelerate the convergence of high-order iterative correctors presented in [START_REF] Damil | An iterative method based upon Padé approximants[END_REF][START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF] and also the iterative algorithm, presented in section (2.1), for viscoelastic structures. As these iterative techniques use a preconditioner matrix, the number of iterations to reach the desired accuracy can be high. The use of a convergence accelerating method, in this case, is then an interesting prospect. We will now analyse the performances of the L * -algorithm, defined by expression [START_REF] Boer | Padé approximant applied to a non-linear finite element solution strategy[END_REF] when it is coupled with an accelerating technique (MMPE).

In the following numerical tests, the operator L * is the tangent matrix at the origin (figure 1(b), U = 0). We only consider trial points 2 and 3 of the example of the cylindrical roof. In Table [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF], the number of iterations to obtain the desired accuracy, 10 -4 , is given for both trial points. The order of truncature of the series ( 11) is chosen to be equal to 5.

Firstly, Table [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF] shows that the MMPE method increases the convergence of the L * -algorithm. Indeed, for point 2, the L * -algorithm needs 10 iterations to get the required accuracy (10 -4 ), whereas when using MMPE, 8 iterations are necessary (respectively 20 and 14 with an accuracy equal to 10 -8 ).

For trial point 3, the number of needed iterations is equal to 14 without MMPE and 10 when using the latter method (respectively 26 and 14 with an accuracy equal to 10 -8 ).

So the use of the L * -algorithm with a small order of truncature decreases the number of iterations. Even if this decrease is less than with the modified Newton algorithm, it leads to a reduction of the total computational cost. Indeed, in this iterative procedure, the computational cost is more or less proportional to the number of vectors to be computed. These numbers of vectors are given in Table [START_REF] Cabay | A polynomial extrapolation method for finding limits and antilimits of vector sequences[END_REF] for several orders of truncature of the polynomial approximations [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF]. The demanded accuracy is, in this case equal, to 10 -8 . Table [START_REF] Cabay | A polynomial extrapolation method for finding limits and antilimits of vector sequences[END_REF] shows that the use of the MMPE method generally reduces the cost by about 40 %. Moreover, associating the MMPE with a small order of without MMPE with MMPE Re=5.8, R(U 0 ) = 0.24 diverge 13

Re= 230, R(U 0 ) = 0.015 > 30 17

Table 6 Flow around a cylinder. Number of iterations to get the desired accuracy 10 -4 with or without MMPE. The order of truncature of the matrix less algorithm (equation 10) is equal to 6.

truncature (order 5) is very attractive. Smaller orders (order 2,3 and 4) have also been tested. In this case, the convergence acceleration of the L * -algorithm is not very efficient.

Numerical evaluation in fluid mechanics

In this section, the L * -algorithm and the MMPE are coupled within a fluid mechanics framework. The presented numerical test is the flow around a cylinder [START_REF] Lahmam | High-order predictor-corrector algorithms[END_REF]. The mesh is composed of 496 quadrilateral elements which lead to 4120 d.o.f.. We consider two trial points: the first one at a low Reynolds number (Re = 5.8) and the second one at Re= 230. The data of the algorithm are as follows: the preconditioner L * is the prediction matrix, order of truncature of the L * algorithm equal to 6, restart parameter k = 10, vectors Y i are the residuals R i .

The numbers of iterations to get the desired accuracy (in this case 10 -4 ) for these two points are given in Table [START_REF] Mesina | Convergence acceleration for the iterative solution of the equation X = AX + f[END_REF]. These results show that the acceleration of convergence by MMPE is also very attractive for fluid mechanics. Indeed, the convergence is achieved even when the initial sequence (3) is divergent.

In the latter case, parameter k is equal to 10 and is necessary to obtain the required accuracy after 13 iterations. This value is also used for the second trial point. In this case, the convergence with the MMPE algorithm is reached after 17 iterations, whereas with the matrix less algorithm the convergence is not reached after more than 30 iterations.

Numerical applications for vibration of viscoelastic beams

In this part, the iterative algorithm presented in section (2.1) is applied to sandwich beams. The considered example is a cantilever beam with small or moderate damping (see [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF]). The beam dimensions and the material characteristics are presented in Table 7. The complex eigenvalue problem has been solved by the iterative algorithm [START_REF] Brezinski | Padé approximants, Handbook of Numerical Analysis[END_REF], with an acceleration by MMPE. The parameters for the algorithm are the following: the preconditioner K * is the real matrix ( 14), restart parameter k=6, the vectors Y i are the residuals R i .

Material properties

Elastic layers Viscoelastic layer

Young's Modulus 6.9 10 Table 8 Vibration of a cantilever viscoelastic beam. Five complex modes are computed.

The proposed iterative algorithm is the one given in section 2.1 with acceleration by MMPE (k=6). Number of computed vectors to get a final residual lower than 10 -6 .

The proposed algorithm has been compared with an algorithm, that is rather similar, and has been presented in [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF]: it is a sort of L * -algorithm using the same preconditioner and it is accelerated by Vectorial Padé approximants. Fifteen calculations, which correspond to the three values of the core loss factor (η c =0.1, 0.6, 1.5) and to the first five modes in each case, have been done. The computations have been done until the euclidian norm of the residual is lower than 10 -6 . We have checked that the two algorithms lead to the same results. We have also compared the number of vectors that have been computed (see Table 8). With the proposed algorithm, each vector corresponds to an iteration. With the one in [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF], each vector is a term in series [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF].

The proposed algorithm converges for all 15 cases, and the number of computed vectors is generally less than with the algorithm Homotopy-Perturbation Padé of [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF]. The algorithm especially converges in the case of a large damping (η c =1.5), where matrix K * is far from the consistent tangent matrix. This test establishes the robustness of the MMPE accelerator.

Conclusion

In this work, two numerical methods are used to accelerate the convergence of iterative methods. From the numerical results one can conclude that the MMPE method is attractive for computational mechanics. This method is very simple and easy to use with any iterative corrector. Moreover, this method does not increase the computational time. These accelerating convergence methods are efficient when the iterative corrector shows slow or poor convergence. The iterative methods, that have been applied in this work, require one matrix (Newton modified) or zero matrix triangulation. However, convergence acceleration is generally not necessary with a classical Newton algorithm or with a high-order Newton corrector [START_REF] Damil | An iterative method based upon Padé approximants[END_REF][START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF].

In this paper, we have tried to show that some convergence acceleration techniques lead to simple tools to build efficient and reliable algorithms that are much cheaper than Newton type algorithms. This is the main result of the paper. Various physical examples have been considered for this numerical evaluation, mainly in the field of thin shell analysis. Nevertheless, the same features have been observed for a test from fluid mechanics which suggests that our conclusions can be upheld for other nonlinear partial differential problems. For instance, the results of Assidi for applications of Padé approximants within elasto-plasticity suggest that convergence acceleration can also work well in this field [START_REF] Assidi | Méthode asymptotique numérique pour la plasticité[END_REF]. By considering the numerical tests, our advice is to use MMPE with at most ten vectors and to partially orthonormalise the projection vectors for a better stability of the algorithm. Some other points are also new, to the author's knowledge, including the association of Vectorial Padé approximants and the modified Newton corrector, the iterative algorithm for viscoelastic eigenvalue problems introduced in Part 2, the comparison between MMPE and Vectorial Padé Approximants.

Fig. 1 .

 1 Fig. 1. Characteristics of the cylindrical thick roof. Response curve of the cylindrical thick roof and the three trial points.

Fig. 2 .

 2 Fig. 2. Characteristics and response curve of the pull-out of an open cylinder. response curve, Fig. (2(b)), three trial points have been chosen to evaluate the proposed corrector.

  Trial point 2 -arc-length (4 iterations, -10.38).

  Trial point 3 -fixed load (4 iterations, -8.27).

Fig. 3 .

 3 Fig.3. Cylindrical thick roof. Logarithm of the residual versus iteration number. MMPE with Y i = R i . In parenthesis the number of iterations to get a residual lower than 10 -8 with a 'pure' Newton corrector and the logarithm of the residual vector after convergence of the Newton algorithm.

Fig. 4 .

 4 Fig. 4. Pull-out of an open cylinder. Logarithm of the residual vector versus the iteration number. Point 1, ( P R D = W R = 0.46). A pure Newton corrector requires 4 iterations to get a norm of the residual vector lower than 10 -5 . Arc-length scheme.

Fig. 5 .

 5 Fig. 5. Cylindrical thick roof. Newton modified corrector with MMPE. Evolution of the logarithm of the residual vector versus iteration number for several choices of vectors Y i , trial point 2

- 8 ,Fig. 6 .

 86 Fig. 6. Cylindrical thick roof. Logarithm of the residual versus iteration number. Trial point 3. Newton modified algorithm. Several ways to use MMPE (k=5, 10, 15, ∞).

  Pull-out of an open cylinder, convergence results from Figure[START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF]. In parenthesis, the details of the CPU time are given. The nred parameter represents the number of components used in the modified scalar product[START_REF] Cadou | Projection technique to improve high order iterative correctors[END_REF].

Fig. 7 .

 7 Fig.7. Cylindrical thick roof (39.366 d.o.f., R(U 0 ) = 0.31 ). Evolution of the logarithm of the residual vector versus the iteration number. Influence of the scalar product for the MMPE method.

Table 1

 1 Cylindrical thick roof. Characteristics of the three trial points.

		Displacement	Load parameter	Norm of the residual
	Point 1	21.8613803	0.521369466	0.27406086
	Point 2	25.7580787	1.12485855	1.3698557
	Point 3	28.3927706	1.79387442	2.8514333

Table 4

 4 Cylindrical thick roof. Number of iterations to get the desired accuracy 10 -4 (in parenthesis 10 -8 ) with or without MMPE. The order of truncature of the matrix less algorithm (eq. 10) is equal to 5.

		L * -algorithm	L * -algorithm
					with MMPE
	Order	5	10	15	5	10	15
	Point 2 100	110	165	70	80	105
	Point 3 130	160	240	70	90	135

Table 5

 5 Cylindrical thick roof. Number of vectors to get the accuracy (10 -8 ) for several orders of truncature and with or without convergence acceleration 4.1 Numerical evaluation in thin shell problems

  1 0 N/m 2 1.794 10 6 N/m 2

	Poisson's ratio	0.3	0.3
	Density	2766 kg/m 3	968.1 kg/m 3
	Thickness	1.523 mm	0.127 mm

Table 7

 7 Material properties for the cantilever beam. Length = 177.8 mm, width =12.7 mm.

	ηc	eigenvalue	Proposed algorithm	algorithm from [25]
	0.1	1	6	8
		2	6	6
		3	4	6
		4	4	6
		5	4	5
	0.6	1	14	18
		2	12	10
		3	10	10
		4	9	10
		5	9	10
	1.5	1	28	> 40
		2	22	28
		3	16	27
		4	17	18
		5	16	15


Projet Iphi

J L Coacolo

σ la ième grandeur rms suivant x ou y, dans le plan transverse. Deux méthodes de minimisation sont testées, la méthode de Powell et la méthode du Simplex. Dans un premier temps, pour réduire la durée de calcul, un code matriciel est utilisé pour la simulation du transport du faisceau. Dans un deuxième temps on utilise un code particulaire prenant en compte les effets non linéaires de la force de charge d'espace. Le résultat obtenu par la première méthode pourra être utilisé comme point de départ de la deuxième méthode.

Figure 1 2 Procédure utilisée A l'aide de N profils obtenus en modifiant le gradient dans un des quadripôles, on détermine N couples de dimensions transverses (

∑ =         - + - = N i i y i y sim i y i x i x sim i x y y y x x x f 1 2 , exp 2 , , exp 2 , exp 2 , , exp ) ( ) ( ) , , , , , ( σ σ
σ σ σ σ ε β α ε β α , avec y y y x x x ε β α ε β α , , , , ,
les paramètres de Twiss et les émittances transverses. Soient une fonction f à plusieurs variables, un point P, une direction définie à partir du vecteur u r , λ un scalaire et la fonction g définie par :

) ( ) ( u P f g λ λ + =
, une coupe de la fonction f suivant la direction u passant par le point P. La méthode de Powell consiste à minimiser f de la manière suivante :

-on part d'un point P 0 et d'une direction u 0 , -par minimisation dans l'espace à une dimension suivant la direction u 0 , on obtient un point P 1 , -on change de direction pour regarder ce qui se passe suivant u 1 , -par itération, on crée P 2 , u 2 …. jusqu'à convergence de l'itération.

Lorsque les vecteurs u i sont les vecteurs unitaires du repère, on parlera d'implémentation de base de la méthode de Powell. Selon la méthode utilisée pour le changement de direction on parlera de méthode de la plus grande pente ou du gradient conjugué. C'est cette dernière qui est prise dans la méthode de Powell et qui est utilisée pour la minimisation.

Méthode du Simplex

La méthode du simplex n'est pas la méthode la plus puissante lorsque l'on veut calculer un minimum d'une fonction à plusieurs variables, mais la méthode la plus simple à mettre en oeuvre. Un simplex dans un espace de dimension N est une figure géométrique composée de N+1 points reliés entre eux par des segments. Par exemple en dimension 2 le simplex est un triangle. Soit f la fonction à minimiser, on calcul les N+1 valeurs f(P i ) et on les ordonne. L'idée est d'encadrer un volume d'espace et de le réduire autour d'un minimum trouvé (ne se trouvant pas forcément dans le volume initial). Pour cela on essaye de projeter le plus mauvais point du simplex de l'autre coté de l'hyperplan que constituent les N autres points et éventuellement appliquer une expansion de ce coté. Ces réflexions sont faites de manière à essayer de conserver le volume du simplex. On peut aussi avoir contraction du volume le long de l'axe du plus mauvais point, ou bien même des contractions selon toutes les dimensions vers le point dont la valeur de la fonction est la plus petite. A chaque itération on recalcule un nouveau simplex en générant un nouveau point pour remplacer le point le plus mauvais. On s'arrête lorsque la dimension du simplex est suffisamment petite.

Validation de la procédure en utilisant un code matriciel pour le transport

Pour juger de la validité de la méthode utilisée et accélérer le calcul, un code matriciel a été développé. Les deux méthodes sont testées avec ce programme. Dans un premier temps, avec les caractéristiques initiales du faisceau données par le code Toutatis [START_REF] Duperrier | TOUTATIS, a RFQ code[END_REF], on simule le transport en faisant varier le gradient dans un des quadripôles et on détermine N couple ) , ( B' T/m -4.9 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.1

x σ (mm) 11.02 8.22 6.71 5.24 3.80 2.40 1.06 0.55 2.77 y σ (mm) 2.74 0.30 1.77 3.42 5.12 6.86 8.63 10.44 13.80

Tableau 2

La procédure utilisée est la suivante : On génère aléatoirement les conditions initiales dans un ellipsoïde centré en P c tel que :

1 ) ( 2 , , 0 , 6 1 ≤ - ∑ = c i c i i i x x x
, avec 0 , i

x la ième coordonnée du point initial et c i

x , la ième coordonnée du point recherché. Après minimisation, on obtient un point P f de coordonnées x i,f . Comme il est difficile de représenter les points P 0 , P c et P f dans un espace à 6 dimensions,

c'est la norme N f = 2 , , , 6 1 )
( c i c i f i i x x x - ∑ = qui est représentée en fonction de la norme N i = 2 , , 0 , 6 1 ) ( c i c i i i x x x - ∑ = .
La routine de minimisation converge vers la solution à déterminer lorsque l'on a N f = 0. La figure 4 représente

) ( i f N F N =
. On peut voir que pour la quasi totalité des points, l'objectif est atteint c'est à dire que la routine de minimisation permet de retrouver les conditions initiales. La figure 5 représente α x,f en fonction de α x,i , et montre que l'on converge vers la bonne valeur pour 104 points sur 110. Pour déterminer chaque paramètre initial, on fait une moyenne sur les N1 valeurs trouvées et l'on affine le calcul en déterminant les écarts types et en éliminant les valeurs qui sont en dehors de 3 écarts types. On peut vérifier sur la figure 6 que la méthode de Powell permet de retrouver avec une grande précision, les valeurs des dimensions rms du tableau 2 Sur cette figure sont représentées les dimensions rms obtenues pour les 110 point P 0 générés aléatoirement. -6,00E+00 -5,00E+00 -4,00E+00 -3,00E+00 -2,00E+00 -1,00E+00 0,00E+00

Prise en compte des non linéarités dans la simulation du transport

Dans les cas précédents, les mesures simulées sont obtenues à partir du transport de faisceau calculé par un code matriciel. Pour tenir compte des non linéarités de la force de charge d'espace nous allons maintenant utiliser le code particulaire, Tracewin [START_REF] Uriot | TraceWin Documentation[END_REF] Les valeurs du gradient dans le dernier quadripôle et les dimensions rms transverses du faisceau correspondantes sont présentées dans le tableau 7. On remarque des écarts maxima de +20 % sur les valeurs en y entre les tableaux 7 et 2. Les dimensions suivant x sont assez proches, l'écart est inférieur à 5 % (excepté pour le minimum).

B' T/m -4.9 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.1

x σ (mm) 10.70 8.00 6.50 5.06 3.66 2.30 1.03 0.67 2.77 y σ (mm)

  exp σ est la valeur de sigma suivant l'axe des x, obtenue expérimentalement pour la ièmemesure et i x sim ,σ la valeur du sigma équivalent obtenue par la simulation.

Figure 2 :

 2 Figure 2 : Procédure utilisée

  endroit du wire scanner. Dans un second temps, avec les couples précédents de minimisation on recherche les caractéristiques initiales. Celles-ci sont définies par le point P c dans l'espace à 6 dimensions. Pour quantifier la sensibilité des deux routines aux conditions initiales P 0 (point de départ pour la routine de minimisation) et déterminer le point P c , on recommence la minimisation pour N 1 points P 0 générés aléatoirement dans un sphéroïde de centre P c . Pour toutes ces études, les enveloppes du faisceau sont celles représentées sur la figure3, à la modification prête en fin de ligne, due à la variation du gradient dans le dernier quadripôle.

Figure 3 :

 3 Figure 3 : Enveloppes du faisceau dans la ligne d'analyses

Figure 4 :

 4 Figure 4 : Influence des conditions initiales sur la minimisation

Figure 10 :

 10 Figure 10 : b x obtenues après minimisation en fonction des valeurs de b x initiales générées aléatoirement.

  Cette étude à pour but de caractériser transversalement le faisceau de protons en sortie du RFQ du projet IPHI. On suppose connues les caractéristiques longitudinales du faisceau. En modifiant le gradient dans un des quadripôles de la ligne d'analyse du projet IPHI, différents profils transverses du faisceau sont mesurés à l'aide d'un wire scanner, voir figure 1. Les effets de la charge d'espace sur la dynamique du faisceau étant importants, il n'est plus possible d'utiliser la méthode de variation des gradients (appelée communément « méthode des 3 gradients ») pour déterminer les caractéristiques du faisceau, c'est à dire les paramètres de Twiss et les émitances. C'est pourquoi un programme de minimisation est développé à partir de travaux entrepris pour d'autres besoins[1]. Les caractéristiques initiales du faisceau sont déterminées en minimisant la distance comprise entre la courbe définie par les points de

	Table des matières
	1 Introduction 1 Introduction
	2 Procédure utilisée
	3 Description des deux méthodes de minimisation	6
	3.1 Méthode de Powell	6
	3.2 Méthode du Simplex
	4 Validation de la procédure en utilisant un code
	matriciel pour le transport mesures ) , ( j i i G σ et celle qui est simulées, où i G est le gradient modifié dans le quadripôle et 7
	4.1 Méthode de Powell j i	8
	4.1.1 Influence des conditions initiales	8
	4.1.2 Prise en compte des erreurs faites sur les mesures	10
	4.2 Méthode du simplex
	4.2.1 Influence des conditions initiales	12
	4.2.2 Prise en compte des erreurs faites sur les mesures	14
	4.3 Prise en compte des non linéarités dans la simulation
	du transport
	5 Utilisation d'un code particulaire
	6 Conclusion	17

  Figure 5: Valeurs de α x obtenues après minimisation en fonction des valeurs de α x initiales générées aléatoirement.Comme le montrent les figures 7 et 8, même si la dispersion des points pour la figure 8 est plus importante que pour le cas précédent, les résultats obtenus sont encore satisfaisants, le tableau 4 le confirme. En filtrant au plus 5 points, l'écart entre les différentes valeurs obtenues et les valeurs recherchées est de l'ordre du pourcent. Les écarts types sont inférieurs à 10 %. Sur la figure 10 est représentée β x,f en fonction de β x,i . Les deux figures montrent que cette méthode donne aussi de bons résultats. Ceci est conforté par le tableau 5.

		3,50E+01	Mesures du Wire scanner simulées et obtenues par la minimisation
		3,00E-01
		1,60E-02 3,00E+01
		1,40E-02 2,50E-01
		2,50E+01
	Dimension rms (m) Ecart relatif final Beta_x,f	6,00E-03 8,00E-03 1,00E-02 1,20E-02 1,00E+01 écart type Points utilisés 1,00E-01 1,50E+01 2,00E+01 2,00E-01 cond. initiales valeurs exactes 0.991 0.08 107 α x valeurs calculées 0.983 1,50E-01	0.01 107 β x (m/rad) ε x (m.rad) 2.03 e -7 105 0.238 3.477 e -6 -2.039 0.09 106 α y 0.240 3.503e -6 -2.066	obtenues par la simulation suivant y 2.77 e -7 107 simulées suivant x obtenues par minimisation suivant x simulées suivant y Dimensions rms à l'endroit du Wire scanner β y (m/rad) ε y (m.rad) 0.03 106 0.507 3.516 e -6 0.518 3.446 e -6
		4,00E-03 5,00E+00 5,00E-02				Tableau 4
		2,00E-03 0,00E+00
	0,00E+00 0,00E+00 0,00E+00 4.2 Méthode du simplex 2,00E-01 0,00E+00 5,00E-02 1,00E-01	4,00E-01 1,50E-01	6,00E-01 Ecart relatif initial 2,00E-01 2,50E-01	8,00E-01 3,00E-01	1,00E+00 3,50E-01 4,00E-01	1,20E+00 4,50E-01
	-6,00E+00 Figure 7 : Influence des conditions initiales sur la minimisation -5,00E+00 -4,00E+00 -3,00E+00 -2,00E+00 -1,00E+00 0,00E+00 Gradient du champ magnétique (T/m) 4.2.1 Influence des conditions initiales Beta_x,i
	2,00E-01 cond. initiales valeurs exactes 0.991 4,00E-01 α x valeurs calculées 0.991 écart type ≈ 0 Points utilisés 104 110 points P 0 sont générées aléatoirement, pour chaque paramètre il y a 6 points sur les 110 6,00E-01 8,00E-01 1,00E+00 1,20E+00 1,40E+00 1,60E+00 1,80E+00 Alpha_x,i Alpha_x,f β x (m/rad) ε x (m.rad) α y β y (m/rad) ε y (m.rad) 0.238 3.477 e -6 -2.039 0.507 3.516 e -6 0.238 3.477e -6 -2.039 0.507 3.516 e -6 ≈ 0 ≈ 0 ≈ 0 ≈ 0 ≈ 0 104 104 104 104 104 Tableau 3 hors des 3 écarts types. On peut remarquer les très bons résultats obtenus. 0,00E+00 2,00E-01 4,00E-01 6,00E-01 8,00E-01 1,00E+00 1,20E+00 4.1.2 Prise en compte des erreurs faites sur les mesures 0,00E+00 5,00E-07 1,00E-06 1,50E-06 2,00E-06 2,50E-06 3,00E-06 3,50E-06 4,00E-06 4,50E-06 0,00E+00 1,00E-06 2,00E-06 3,00E-06 4,00E-06 5,00E-06 6,00E-06 De même que pour la méthode de Powell, ) ( i f N F N = 2,00E-01 4,00E-01 6,00E-01 1,20E+00 1,40E+00 1,60E+00 1,80E+00 est représentée sur la figure 9. 0,00E+00 8,00E-01 1,00E+00 Ecart relatif final
						Ecart relatif initial
			Figure 9 : Influence des conditions initiales sur la minimisation

Figure 6 : Dimensions rms du faisceau simulées et obtenues par minimisation à l'endroit du wire scanner La figue 6 montre l'excellent accord entre les valeurs des dimensions rms simulées et obtenues par la minimisation en x et y. Les valeurs obtenues par la minimisation ont été reportées pour les les 110 cas générés aléatoirement. Dans le tableau 3, sont reportées les conditions initiales recherchées, celles déterminées statistiquement après minimisation et les différents écarts types correspondants.

Dans le cas précédent nous supposions que les mesures étaient connues exactement. Pour tenir compte des erreurs faites sur celles-ci, on génère aléatoirement une erreur comprise entre plus ou moins 5% des valeurs précédentes, ceci pour chaque jeu de conditions initiales. Figure 8 : ε y obtenues après minimisation en fonction des valeurs de ε y initiales générées aléatoirement.

  Les remarques faites pour la figure 6 sont de nouveau valable pour la figure 11. On remarquera une fois de plus l'excellent accord entre les valeurs simulées et celles obtenues par la minimisation. Prise en compte des erreurs faites sur les mesures Les résultats étant proches de ceux donnés par la méthode de Powell je ne reporterai ici que le tableau 6 donnant les valeurs moyennes des caractéristiques du faisceau déterminées et les écarts types associés.

	cond. initiales valeurs exactes 0.991 α x	β x (m/rad) ε x (m.rad) 0.238 3.477 e -6 -2.039 α y	β y (m/rad) ε y (m.rad) 0.507 3.516 e -6
	valeurs calculées 0.991	0.238	3.477e -6 -2.039	0.507	3.516 e -6
	écart type	≈ 0	≈ 0	≈ 0	≈ 0	≈ 0	≈ 0
	Points utilisés	107	107	107	107	107	107
				Tableau 5
	4.2.2 Tableau 6
							Dimensions rms à l'endroit
	obtenues par minimisation suivant x simulées suivant y obtenues par minimisation suivant x du Wire scanner Les résultats sont satisfaisants et les erreurs sont de l'ordre du pourcent et des écarts types Dimension rms (m) simulées suivant x inférieurs à 10 %.
			Gradient du champ magnétique (T/m)

Figure 11 : Dimensions rms du faisceau simulées et obtenues par minimisation à l'endroit du wire scanner

  et simuler le transport jusqu'au wire scanner L'optimisation est faite avec le code matriciel.

	cond. initiales	α x	β x (m/rad) ε x (m.rad)	α y	β y (m/rad) ε y (m.rad)
	valeurs exactes 0.991	0.238	3.477 e -6 -2.039	0.507	3.516 e -6
	valeurs calculées 0.982	0.240	3.503e -6 -2.066	0.518	3.438 e -6
	écart type	0.08	0.01	2.03 e -7	0.09	0.03	2.79 e -7
	Points utilisés	108	109	107	107	107	108

Dans un premier temps la routine PICNIR est utilisée pour déterminer la force de charge d'espace. Le calcul du champ généré par la charge d'espace se fait en considérant que la charge est distribuée dans un anneau circulaire. La section transverse elliptique du faisceau est prise en compte lorsque l'on calcul le rayon effectif de l'anneau chargé

  Pour prendre en compte la composante non linéaire de la force de charge d'espace, l'optimisation est faite avec un code particulaire, PARTRAN. Les méthodes du simplex et de Powell sont utilisées. La simulation du transport est faite avec la routine PICNIC et l'optimisation utilise la routine PICNIR. Figure 12 : Dimensions rms du faisceau simulées et obtenues par minimisation à l'endroit du wire scanner Sur la figure 12, pour chaque mesure nous avons représenté les 110 points simulés, chaque point représente une mesure avec une erreur générée aléatoirement. Naturellement les 110 points obtenus par la minimisation sont aussi représentés. On remarquera que la dispersion est plus grande pour les points simulés. Il semblerait que la minimisation corrige une partie des erreurs dues aux mesures. D'autre part il y a un bon accord entre les valeurs simulées et obtenues par la minimisation.Dans le but de caractériser le faisceau en sortie du RFQ du projet IPHI, un programme d'optimisation a été développé. La détermination de ces caractéristiques s'appuie sur les mesures faites avec le wire scanner. N'ayant pas à l'heure actuelle d'acquisitions de données, ce programme est testé avec des valeurs obtenues par la simulation. Deux méthodes de minimisation sont regardées. Les résultats obtenus avec les deux méthodes en utilisant un code matriciel pour la simulation et pour la minimisation sont excellents. Dans l'étape suivante, des erreurs de ±5% sur les dimensions transverses du faisceau sont prises en compte. Les écarts entre les valeurs recherchées et obtenues par la minimisation sont inférieurs à 3%. Pour tenir compte des effets non linéaires de la charge d'espace lors du transport du faisceau dans la ligne d'analyse, la simulation a été faite avec le code Tracewin, l'optimisation prend toujours en compte des erreurs faites sur les mesures. Les valeurs obtenues s'éloignent des valeurs recherchées et l'écart peut aller jusqu'à 50%. Pour remédier à ceci et accroitre la qualité des résultats, la simulation et l'optimisation se font avec le code particulaire. Ceci permet d'avoir des écarts inférieurs à 30% entre les caractéristiques recherchées et obtenues par la minimisation. Il n'apparaît pas nécessaire de développer d'autres méthodes de minimisation pour essayer d'accroitre la précision des résultats vu que ce qui est limitatif est l'ensemble des erreurs faites sur les mesures.

	Mesures du Wire scanner simulées et obtenues par la minimisation
	cond. initiales valeurs exactes 0.991 α x 1,60E-02 1,80E-02	β x (m/rad) ε x (m.rad) 0.238 3.478 e -6 -2.039 α y	0.507 β y (m/rad) ε y (m.rad) 3.516 e -6
				Powell
	valeurs calculées 1.007 1,40E-02	0.270	4.584 e -6 -1.178	0.302	5.936 e -6
	écart type	0.06	0.01	2.02 e -7	0.06	0.01	3.66 e -7 Dimensions rms à l'endroit
	Tableau 7 Les résultats présentés dans le tableau 8 pour les deux méthodes montrent que ce qui est dit Simplex valeurs calculées 1.003 0.270 4.600e -6 -1.179 0.303 5.933 écart type 0.06 0.009 1.85 e -7 0.06 0.01 3.68 e -7 Points utilisés 109 108 108 106 106 106 8,00E-03 1,00E-02 simulées suivant x obtenues par minimisation suivant x simulées suivant y obtenues par minimisation suivant y Dimension rms (m) 3.16 0.51 2.21 4.16 6.16 8.21 10.31 12.44 16.41 Points utilisés 110 110 110 110 110 110 1,20E-02 du Wire scanner
	ci-dessus se traduit par une certaine dissymétrie dans la qualité des résultats obtenus horizontalement et verticalement. On remarquera l'excellent accord entre les deux routines Tableau 10 6,00E-03
	d'optimisation. 4,00E-03
	2,00E-03
	cond. initiales valeurs exactes 0.991 α x 0,00E+00 5 Utilisation d'un code particulaire β x (m/rad) ε x (m.rad) 0.238 3.478 e -6 -2.039 α y -6,00E+00 -5,00E+00 -4,00E+00 -3,00E+00 -2,00E+00 -1,00E+00 β y (m/rad) ε y (m.rad) 0.507 0,00E+00 3.516 e -6 Gradient du champ magnétique (T/m)
				Powell
	valeurs calculées 0.939	0.261	4.332e -6 -0.937	0.235	6.696 e -6
	écart type	0.06	0.01	2.01 e -7	0.05	0.008	5.06 e -7
	Points utilisés	110	110	110	110	110	110
				Simplex
	valeurs calculées 0.937 écart type 0.06 cond. initiales α x Points utilisés 105 valeurs exactes 0.991	0.262 0.01 β x (m/rad) ε x (m.rad) 4.338e -6 -0.939 2.01e -7 0.05 α y 98 107 92 0.238 3.478 e -6 -2.039 Powell	0.235 0.09 β y (m/rad) ε y (m.rad) 5.209e -6 2.02e -7 94 0.507 3.516 e -6 110
	valeurs calculées 1.10 écart type 0.18	0.240 0.02	3.807e -6 -1.400 Tableau 8 6.6 e -7 0.23	0.458 0.08	3.656 e -6 .9 e -7
	Ceci peut être dû à la routine de charge d'espace prenant assez mal en compte les dimensions Le tableau 9, obtenu avec la routine PICNIC, montre que les valeurs suivant l'axe y se écart type 0.13 0.019 4. e -7 0.23 0.09 1 e -6 routine PICNIC déterminant la force de charge d'espace en 3 dimensions. valeurs calculées 1.054 0.260 4.46e -6 -1.543 0.479 3.47 e -6 transverses. Pour essayer d'accroitre la précision du calcul dans le plan vertical, on utilise la Points utilisés 94 95 99 85 98 100 Simplex 6 Conclusion
	rapprochent de celle du tableau 2. En ce qui concerne la variable y, le tableau 10 donne de Points utilisés 97 97 95 85 98 100
	meilleurs résultats que le tableau 7
				Tableau 11
	B' T/m Le tableau 11 montre les résultats obtenus avec les deux routines de minimisations et en -4.9 -4 -3.5 -3 -2.5 -2 -1.5 -1	-0.1
	σ (mm) 11.06 tenant compte d'une erreur sur les mesures de ± 5%. 8.26 6.75 5.27 3.83 x y σ (mm) 2.94 0.47 2.01 3.80 5.64 Les résultats sont meilleurs que ceux obtenus avec le code matriciel, surtout suivant l'axe des 2.43 1.12 0.66 2.81 7.53 9.46 11.43 15.08 y. Le tableau montre des résultats qui diffèrent des valeurs simulés de 30% au plus.
				Tableau 9


D Zapico Álvarez

F Bertrand

J.-M Mataigne
email: jean-michel.mataigne@arcelormittal.com

M.-L Giorgi

Mechanisms and kinetics of the inhibition layer breakdown in the case of Ti IF steel grades galvanized in GA baths

Keywords: Al-Fe-Zn ternary system, interfacial reactions, galvannealing, inhibition layer, inhibition breakdown

A reaction mechanism is proposed accounting for the inhibition layer breakdown in the case of Ti IF steels. This mechanism can be summarized in two steps: 1) the transformation of the bi-layered inhibition layer into a mono-layered one, made of the δ (FeZn7) phase only, as a result of the Zn enrichment at the steel/inhibition layer interface, and 2) the local nucleation of the Г (Fe3Zn10) phase at the steel grain boundaries, breaking the δ layer off, when the Zn concentration at these locations becomes high enough (outburst formation). The kinetics of this reaction strongly depends on the Ti IF steel chemistry. The effect of the steel chemical composition on the inhibition layer breakdown kinetics would be ruled by the rate of Zn diffusion at the steel grain boundaries and the ability of the steel to accumulate the Zn atoms within a critical volume for the nucleation of Г phase at these locations.

INTRODUCTION

The development of a conventional GalvAnnealed (GA) coating on a Ti-stabilized Interstitial Free (Ti IF) steel grade passes through different and complex reactions: the inhibition layer formation, the inhibition layer breakdown, the liquid Zn consumption and the enrichment in Fe of the coating by solid state diffusion. A fine knowledge of the kinetics of each of these reactions may considerably help process and production engineers to define the thermal cycle allowing reaching the customer targets in terms of iron content and Fe-Zn phase distribution within the coating.

A very thin layer of Al-rich intermetallic compounds nucleates very rapidly (around 0.1 -0.2 s 1) ) on the steel surface during galvanizing in Al-containing Fe-saturated Zn baths. This layer is commonly known as inhibition layer, as it inhibits, at least temporarily, the alloying reactions between Fe and Zn. It has been demonstrated in a previous study [START_REF] Zapico Álvarez | Proc. of the 9 th Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF] that the inhibition layer which forms in typical baths for galvannealing production (Fe-saturated Zn baths with Al content ranging from 0.10 to 0.135 wt.% at 460 °C) is actually composed of a very thin layer of Fe2Al5Znx (becoming discontinuous for low Al contents) and a thicker layer of Al-saturated δ (FeZn7) on its top. This microstructure allows thermodynamic equilibrium to be reached along all interfaces and the system does therefore not require any other interfacial reaction.

However, increased temperature in the induction furnace, located at the exit of the galvanizing bath, activates Fe and Zn inter-diffusion through the inhibition layer and, after a while, this layer is destabilized. This reaction, called the inhibition layer breakdown, is a key step in the galvannealing process as it leads to the sudden development of Fe-Zn phases by consuming liquid Zn. It is therefore crucial to understand the mechanisms and kinetics of this reaction.

Different mechanisms have been proposed in the literature to explain the inhibition layer breakdown. Even if a very small group of researchers [START_REF] Mcdevitt | Zinc-Based Steel Coating Systems: Production and Performance[END_REF][START_REF] Nitto | [END_REF] reported that the Fe-Zn reactions are activated by the dissolution of the inhibition layer in liquid Zn, it is generally admitted that an outburst reaction is the responsible for the inhibition layer breakdown in the case of Ti IF steels [5][6][7][START_REF] Hisamatsu | Proc. of the 1 st Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Leprêtre | Étude des mécanismes réactionnels de la galvanisation[END_REF][START_REF] Leprêtre | Zinc-Based Steel Coating Systems: Production and Performance[END_REF][START_REF] Jordan | [END_REF][START_REF] Jordan | Zinc-Based Steel Coating Systems: Production and Performance[END_REF][START_REF] Guttmann | [END_REF][14][START_REF] Guttmann | Proc. of the 3 rd Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Dionne | Proc. of Int. Symp. on Materials in the Automotive Industry[END_REF][START_REF] Saito | [END_REF][START_REF] Kim | Zinc-Based Steel Coating Systems: Production and Performance[END_REF] . Two mechanisms can be found to explain the outburst formation. A few researchers [5][6][7] submitted that the diffusion of Fe towards liquid Zn causes the growth of the inhibition layer and, as a result, an Al depletion at the solid / liquid interface, leading to the outburst formation. A much more important group of researchers [START_REF] Hisamatsu | Proc. of the 1 st Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Leprêtre | Étude des mécanismes réactionnels de la galvanisation[END_REF][START_REF] Leprêtre | Zinc-Based Steel Coating Systems: Production and Performance[END_REF][START_REF] Jordan | [END_REF][START_REF] Jordan | Zinc-Based Steel Coating Systems: Production and Performance[END_REF][START_REF] Guttmann | [END_REF][14][START_REF] Guttmann | Proc. of the 3 rd Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Dionne | Proc. of Int. Symp. on Materials in the Automotive Industry[END_REF][START_REF] Saito | [END_REF][START_REF] Kim | Zinc-Based Steel Coating Systems: Production and Performance[END_REF] reported that the triggering phenomenon in the mechanism of outburst formation is the diffusion of Zn towards the steel substrate. Nishimoto et al. 5) were the first researchers to reveal that the outbursts are formed at the steel grain boundaries for long immersion times in GA baths. Hisamatsu 8) confirmed afterwards that the steel grain boundaries are preferential locations for outburst formation. Later, other authors [START_REF] Leprêtre | Étude des mécanismes réactionnels de la galvanisation[END_REF][START_REF] Leprêtre | Zinc-Based Steel Coating Systems: Production and Performance[END_REF] provided new experimental evidence of the formation of outbursts at the emergence of the steel grain boundaries.

The kinetics of the inhibition layer breakdown is strongly influenced by the galvanizing conditions and the base steel features. It is well known that the higher the Zn bath Al content, the slower the inhibition layer breakdown kinetics [START_REF] Leprêtre | Étude des mécanismes réactionnels de la galvanisation[END_REF] . As for the steel substrate, it is generally admitted that the finer the grain size, the faster the inhibition layer breakdown occurs [START_REF] Hisamatsu | Proc. of the 1 st Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Jordan | [END_REF][START_REF] Jordan | Zinc-Based Steel Coating Systems: Production and Performance[END_REF][START_REF] Toki | The Physical Metallurgy of Zinc Coated Steel[END_REF] . The steel chemical composition also plays an important role on the kinetics of the inhibition layer breakdown. Elements such as C, N, S and P may segregate to the steel grain boundaries during recrystallization annealing, hindering the diffusion of Zn at these locations and therefore delaying the inhibition layer breakdown 5,[START_REF] Hisamatsu | Proc. of the 1 st Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Leprêtre | Étude des mécanismes réactionnels de la galvanisation[END_REF][START_REF] Guttmann | [END_REF][14][START_REF] Guttmann | Proc. of the 3 rd Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Dionne | Proc. of Int. Symp. on Materials in the Automotive Industry[END_REF][START_REF] Toki | The Physical Metallurgy of Zinc Coated Steel[END_REF] . On the contrary, Ti can precipitate with the previous elements, limiting their segregation in solid solution at the steel grain boundaries and, therefore, promoting their cleanliness and a faster inhibition layer breakdown 5,[START_REF] Hisamatsu | Proc. of the 1 st Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Guttmann | [END_REF][START_REF] Guttmann | Proc. of the 3 rd Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF][START_REF] Toki | The Physical Metallurgy of Zinc Coated Steel[END_REF][START_REF] Abe | [END_REF] . In this work, the reactivity of Ti IF steels will be assessed by calculating the excess solute Ti in mole fraction, The main objective of this work is to study the mechanisms and the kinetics of the inhibition layer breakdown in the case of Ti IF steels with different excess solute Ti. The work is divided into the following stages: 1) the inhibition layer formed on steel samples galvanized in a Zn bath with 0.127 wt.% Al was observed; the evolution of this inhibition layer was then observed when the galvannealing temperature was increased; 2) the time needed for the inhibition layer breakdown was measured as a function of galvannealing temperature and excess solute Ti of the steel; 3) key mechanisms accounting for the inhibition layer breakdown were proposed together with a simple model to predict the time needed for the inhibition layer breakdown.

EXPERIMENTAL PROCEDURE Four different commercial Ti IF full hard steels, A, B, C and D, have been selected for this study. The thicknesses of these steels are very similar and equal to 0.75, 0.70, 0.74 and 0.74 mm respectively. Their average chemical composition is shown in table 1 together with their excess solute Ti (equation 1). The analysis was performed by means of Combustion and Infrared Detection for nitrogen (HORIBA EMGA-620W) and Spark-OES (SPECTROLAB M10) for the other elements in ArcelorMittal laboratories. The commercial steel sheets were cut into samples with dimensions of 210 x 99 mm 2 . These samples were cleaned with ethanol to remove residual oil from former manipulations. Samples were then annealed at 800 °C under controlled atmosphere of N2 and 5 vol.% H2 with -60 °C dew point and galvanized in a Rhesca Hot Dip Process simulator. The galvanizing conditions are summarized in table 2.

Samples galvanized in the Rhesca simulator were then cut into samples with dimensions of 120 x 30 mm 2 and submitted to isothermal heat treatments in a Gleeble 3500 device, able to provide very high heating and cooling rates (around 200 °C/s) as well as very steady and accurate holding temperatures and times. The isothermal treatments were performed on each steel at holding temperatures of 430, 445 and 460 °C and several times in order to determine accurately the time needed to achieve the inhibition layer breakdown at each of these temperatures. In the case of steel C, supplementary holding temperatures of 425, 438 and 452 °C have been tested in order to confirm the tendencies observed. The thermal homogeneity area in the Gleeble 3500 samples was defined on the basis of a thermal homogeneity criterion of ± 5 °C and corresponded to a squared area of 10 x 10 mm 2 in the centre of the sample, where the thermocouple was placed. It is obviously in this area of the sample that the characterizations have been done. Some cross-section specimens were prepared by Ar ion beam Cross Section Polishing (CSP, JEOL SM-09010) and characterized in a Scanning Electron Microscope (SEM, JEOL JSM-7001F) in order to check the structure of the inhibition layer on the galvanized samples before heat treatment in the Gleeble 3500 pilot.

The state of the inhibition layer of the samples annealed in the Gleeble 3500 pilot was characterized through SEM (LEO DSM 982) observations of the sample surface after selective electrochemical dissolution. The electrochemical cell used for these dissolutions was composed of three electrodes: a standard calomel electrode (SCE) used as reference electrode, a Pt counter electrode used as cathode and the sample to be dissolved in the form of a 32 mm diameter disc (working electrode) used as anode. The electrolyte used is an aqueous pH 4 solution composed of 200 g/L NaCl (Fisher Chemical), 40 g/L ZnSO4•7H2O, 27 g/L CH3COONa and 11 mL/L CH3COOH (VWR BDH Prolabo). The electrochemical potential to be imposed in order to achieve the selective removal of Zn and Αl-free Fe-Zn phases is -800 mV. In these conditions, the whole inhibition layer (Al-saturated δ and Fe2Al5Znx) is kept at the sample surface because of its high Al content.

The inhibition layer breakdown is said to have occurred when holes in the inhibition layer, created by outbursts formation at the steel grain boundaries, are observed. The determination of the inhibition layer breakdown by surface SEM observations is illustrated in figure 1, where an example corresponding to steel B annealed in the Gleeble 3500 device at 430 °C for 30 (left; intact inhibition layer) and 35 s (right; broken inhibition layer) is shown. The diffusion of Zn through the inhibition layer, leading to its breakdown, mainly takes place during annealing in the Gleeble device but may begin during galvanizing and cooling down in the Rhesca simulator. This part of the diffusion process is minimized by using a short immersion time (0.7 s) followed by a forced cooling (6 °C/s). As explained later, the measured experimental time needed for the inhibition layer breakdown is corrected to take into account the diffusion of Zn in the Rhesca simulator. In order to understand the mechanisms accounting for the inhibition layer breakdown, mechanically polished cross-section specimens were also prepared from some well-chosen samples that had only known the inhibition layer breakdown after heat treatment in the Gleeble 3500 device in the thermal homogeneity area but not elsewhere (lower temperature compared to the centre of the sample). The cross-section specimens were prepared from samples taken at different positions along the longitudinal axis of the original Gleeble 3500 samples and characterized by Energy Dispersive Spectroscopy (EDS, Bruker AXS XFlash 4010 SDD) in a SEM (JEOL JSM-7001F).

EXPERIMENTAL RESULTS

The structure of the inhibition layer on the samples galvanized in the Rhesca simulator has been found to be totally consistent with the one observed in a previous study [START_REF] Zapico Álvarez | Proc. of the 9 th Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF] on a steel galvanized with the same bath temperature and aluminum content (table 2) in industrial line: the inhibition layer is in both cases composed of a very thin continuous Fe2Al5Znx layer and a thicker δ layer on its top. The presence of isolated ζ rod crystal has also been observed (figure 2).

Cross-section samples were characterized at different locations along the longitudinal axis of a Gleeble sample after galvanizing in the Rhesca simulator and annealing at 430 °C for 60 s in the Gleeble 3500 device (figure 3). A gradient of temperature existed from the left-hand side (lowest temperature) to the right-hand side (highest temperature) micrographs. At the lowest temperature (figure 3a), the Zn coating remained solid and no diffusion phenomena could be observed: the inhibition layer is still composed of a very thin continuous layer of Fe2Al5Znx and a thicker layer of δ on its top. At intermediate temperature (figure 3b), the Zn coating was liquid and the diffusion of Fe from the steel to the coating took place with two consequences: 1) the growth of the δ phase can be observed and 2) the Fe supersaturation present in the coating leads to the formation of a large amount of evanescent ζ crystals on top of the δ phase during the fast cooling. At the same time, diffusion of Zn from the coating to the steel led to the complete disappearance of the thin layer of Fe2Al5Znx, replaced by δ. As the inhibition layer is now composed of the δ phase only, the associated EDS signal for Al is thicker but less intense than in figure 3a. At the highest temperature, here equal to 430 °C (figure 3c), an outburst can clearly be observed and almost no Al is detected by EDS mapping analysis. The experimental holding time performed in the Gleeble 3500 pilot necessary for the inhibition rupture thold,exp is plotted against the experimental holding temperature Thold,exp in figure 4 for each of the Ti IF steels studied. The uncertainty in each point represents the difference between the highest holding time for which the inhibition rupture has not occurred yet and the lowest time for which the inhibition rupture has been observed. The time for inhibition breakdown has been chosen as the point in the middle of the interval defined by these two times. It is obvious that the smaller this interval, the more accurate the determination of the time for inhibition breakdown, but a compromise between accuracy and number of trials had to be found. At least in the range of the temperatures tested, the steel reactivity classification is as expected according to the values of excess Ti x (table 1): steel D is more reactive (the inhibition layer breakdown is achieved faster)

than steel C and this latter is more reactive than steel B, steel A being the less reactive one. The initial state of the inhibition layer can be described by diffusion path 1, which illustrates a microstructure with all its interfaces (steel / Fe2Al5Znx, Fe2Al5Znx / δ, δ / liquid Zn) at thermodynamic equilibrium (figures 2 and 3a) [START_REF] Zapico Álvarez | Proc. of the 9 th Int. Conf. on Zinc and Zinc Alloy Coated Steel Sheet[END_REF] .

During the galvannealing heat treatment, the inhibition layer grows by Fe diffusion towards the liquid Zn. At the same time, Zn diffuses towards the steel matrix. The steel / Fe2Al5Znx interface gets progressively enriched in Zn or, what is the same, impoverished in Al (black arrow in figure 5). The Fe2Al5Znx phase disappears and thermodynamic equilibrium between steel and δ is reached along the steel / inhibition layer interface (figure 3b). This microstructure is represented by the diffusion path 2 in figure 5.

During all this time, the steel grain boundaries are short-circuit diffusion paths for Zn and become preferentially enriched in this element. When the Zn concentration at these locations is high enough, the steel / δ interface is locally destabilized: the diffusion path crosses the ternary domain involving Fe, δ and Г. The Г phase therefore nucleates at the steel grain boundaries. Its volume expansion breaks the δ layer off locally and the Г phase gets in direct contact with liquid Zn. This microstructure is represented by path 3, where the dotted line indicates that no equilibrium is possible between the Г phase and liquid Zn. This non-equilibrium situation would lead to the dissolution of the Г phase surface and the nucleation of new Fe-Zn phases on its top. This microstructure, commonly known as outburst, is not discussed here as the present analysis focuses only on the phenomena taking place until the inhibition breakdown occurs.

A model describing the inhibition layer breakdown kinetics through outburst formation is presented below. The constants of this model were determined from the experimental results shown in figure 4. Due to the complexity of the reaction mechanisms proposed above for the inhibition layer breakdown, the following simplifications were made in order to facilitate modeling work:

• The Fe diffusion is neglected so that the thickness of the inhibition layer can be considered constant.

• The rate of Zn diffusion from the middle of the steel grains towards the grain boundaries is assumed to be infinite, so that the Zn concentration along the steel / inhibition layer interface can be assumed constant.

• The diffusion of Zn through the Fe2Al5Znx and δ layers is described by a global diffusion coefficient taking into account the Zn diffusion in the entire inhibition layer.

Based on these assumptions, a simplified diffusion model can be proposed from Fick's law in steady state. The Zn diffusion flux through the inhibition layer . .l i Zn φ (mol m -2 s -1 ) can be expressed as follows:

. . 0 . . . . . . . . . . exp l i Zn a l i Zn l i Zn l i Zn l i Zn l i Zn c T R E D c D c D ∇         - = ∇ = ∇ - = φ (2) where .
.l i Zn D (m 2 s -1 ) is the diffusion coefficient of Zn in the inhibition layer, whose variation with temperature T has been expressed through an Arrhenius equation, D0 (m 2 s -1 ) being the pre-exponential factor, Ea (J mol -1 ) the activation energy for diffusion and R (J mol -1 K -1 ) the universal gas constant. As the Zn concentrations along all interfaces of the inhibition layer and its thickness are assumed to be constant, the Zn concentration gradient in the inhibition layer . .l i Zn c ∇ is also constant.

The critical time at which the inhibition layer breakdown is achieved tc (s) is such that a critical quantity of Zn Qc (mol m -2 ) diffuses through the inhibition layer. For a given steel, this critical quantity is expected to be constant, independent of time:

          - ∇ = =
Q t +           ∇ = (4)
where Qhold,exp is the part of the critical quantity of Zn that diffuses through the inhibition layer to the steel grain boundaries during the isothermal treatment. In the model proposed here (equation 2), the diffusion through the inhibition layer only is considered. Therefore, the values of Ea and D0 obtained should not vary with the steel grade. The differences in reactivity should only be linked to different values of Qc, the quantity of Zn needed to cause the inhibition layer breakdown. However, as it can be seen in figure 6, the slopes of the straight lines are different for each steel: the value of Ea varies from one steel to another and cannot finally be considered as the activation energy of Zn diffusion through the inhibition layer. This energy must so be interpreted as the activation energy of Zn diffusion in the grain boundaries, the actual rate limiting step, which varies with excess Ti

x .

The activation energy Ea decreases while the parameter ( )

. . 0 l i Zn c c D Q ∇ increases when excess Ti x increases (table 3).
When the value of excess Ti x is higher, the grain boundaries are cleaner. In this case, Zn diffusion in the grain boundaries is faster and the corresponding Ea is lower. However, a higher amount of Zn flux is needed to trigger the outburst reaction, because it requires a certain amount of Zn to be accumulated in a critical volume to stabilize the first nuclei of Γ phase.

No experiment at higher temperature has been performed in this work but figure 6 suggests that an inversion of the steel reactivity should be observed when the temperature increases enough. Accumulating the critical amount of Zn in the required volume for the nucleation of the first Γ nuclei can indeed become less rapid in the steel with too clean grain boundaries.

CONCLUSIONS

• The inhibition layer breakdown by outburst formation occurs through Zn diffusion in the steel grain boundaries and nucleation of a Fe-Zn phase at these locations. A sequence of reaction mechanisms accounting for this phenomenon has been proposed and explained by means of the concept of diffusion paths in the Al-Fe-Zn ternary phase diagram at the galvannealing temperature. The diffusion of Zn towards the steel enriches the steel / inhibition layer interface in this element. As a result of this enrichment, the disappearance of the Fe2Al5Znx layer occurs first, stabilizing a steel / δ interface. Further enrichment of the steel grain boundaries in Zn leads to the nucleation of the Г phase at these locations, which breaks the inhibition layer off locally.

• Trying to describe the inhibition layer breakdown kinetics by a simplified Zn diffusion model through the inhibiting layer only does not hold. A diffusion model has to take into account Zn diffusion into the steel grain boundaries.

• The effect of the steel chemical composition on the inhibition layer breakdown kinetics has been investigated. The parameter excess Ti x

proposed in this study seems to successfully describe the cleanliness of the steel grain boundaries and, consequently, the reactivity of the steel with respect to the inhibition layer breakdown. The experimental results obtained at holding temperatures from 425 to 460 °C have shown that high values of excess Ti x , corresponding to the clean grain boundaries, accelerate the inhibition layer breakdown. This behaviour is predicted to be reversed for higher temperatures, when the critical amount of Zn in the grain boundaries required for Γ nucleation should be reached first for steels with lower excess Ti x .

Figure 1

 1 Figure 1 Surface SEM micrographs of steel B annealed in the Gleeble 3500 pilot at 430 °C for 30 (left; intact inhibition layer) and 35 s (right; broken inhibition layer).

Figure 2

 2 Figure 2 Cross-section SEM micrographs of the inhibition layers obtained on a sample galvanized in the Rhesca simulator (left) and on a sample from a previous study 2) galvanized with the same conditions in industrial line (right).

Figure 3

 3 Cross-section SEM micrographs and EDS mappings for Al for samples taken at different locations along the longitudinal axis of a Gleeble sample: (a) lowest temperature, (b) intermediate temperature, (c) highest temperature equal to 430 °C (steel A galvanized in the Rhesca simulator and annealed in the Gleeble 3500 device at 430 °C for 60 s).

Figure 4

 4 Figure 4 Experimental holding time needed for the inhibition layer breakdown against the experimental holding temperature performed in the Gleeble 3500 pilot.

Figure 5

 5 Figure 5 Schematic representation of diffusion paths in the Al-Fe-Zn ternary phase diagram at the holding temperature describing the mechanisms accounting for the inhibition layer breakdown.

  tgrowth (s) is the time necessary for the growth of a covering inhibition layer (around 0.1 -0.2 s after the steel strip enters the galvanizing bath1) ). Both parameters a E and treatments, the activation energy Ea can be directly extracted from the slope of straight lines expected to be obtained when representing the results of figure4in the form of experimental holding time in logarithmic scale against the reciprocal of the holding temperature:

  cannot be determined from the y-intercept values of these straight lines (equation 4)because Zn has already diffused during galvanizing in the Rhesca pilot simulator. The real time for the inhibition layer breakdown tc is therefore underestimated when represented by thold,exp. As a result, this latter needs to be corrected by taking into account the diffusion of Zn during the Rhesca thermal cycle, which can be calculated by integrating the right hand side expression of equation 3 along the Rhesca thermal cycle (immersion in the galvanizing bath and subsequent cooling) once Ea has been obtained. Finally, the real holding time necessary for the inhibition rupture tc is plotted in logarithmic scale against the reciprocal of experimental holding temperature Thold,exp for the different Ti IF steels studied (Figure6). Straight lines are obtained in all cases. This indicates that the time necessary for the inhibition layer breakdown does indeed vary with temperature according to an Arrhenius-type law. The values of Ea and 6 for the different steels are presented in table3.

Figure 6

 6 Figure 6 Real holding time needed for the inhibition layer breakdown (logarithmic scale) against the reciprocal of the experimental holding temperature.Table 3 Values of

Table 1

 1 Average chemical composition (x 10 -3 wt.%) of the Ti IF steels studied.

	Steel C	Mn	P	S	N	Si	Cu	Ni	Cr	Sn	As	Nb Mo	B	Ti	Al	Ti x	excess	*
	A	3.0 100.2 12.4 14.1 3.4 9.3 25.5 17.4 21.8 5.7 2.7 0.1 2.6 0.5 57.4 44.6	0.7
	B	2.1 89.9 13.6 15.2 3.9 4.0 24.8 13.9 22.3 3.0 1.8	-	2.2	-	83.6 29.1	1.8
	C	2.0 110.0 9.1 9.0 2.1 4.0 17.0 14.3 18.2 5.7	-	0.1 1.8 0.0 67.3 49.1	2.7
	D	2.2 99.9 6.9 12.5 3.4 4.6 36.6 20.3 27.9 3.9 2.0 0.1 1.9	-	84.1 55.8	4.2

* expressed in mole fraction

Table 2

 2 Galvanizing conditions for the Ti IF steels studied.

	wAl (wt.%) a	Tbath (°C) b	tim (s) c
	0.127 ± 0.003	460 ± 2	0.7

a bath aluminium content; b bath temperature; c immersion time.