B Claude

L Duigou

G Girault

Y Guevel

J M Cadou

Numerical comparison of eigenvalue algorithms for vibroacoustic problems

Keywords: Linear vibroacoustic problem, eigenvalues, pertubation method

This paper concerns the study of vibroacoustic problems. By considering a displacement-pressure formulation, a non-symmetric eigenvalue problem is obtained. In order to solve it, three numerical schemes are compared: the classical ARPACK solver, an indicator method (initially proposed in B. Claude et al. Comptes Rendus Mécaniques, 2017, 345(2)) which has the property to be null at the eigenvalues, and an original method based on the analysis of Taylor series expansions near a singularity. Numerical results show all the evaluated numerical methods give accurate results but the indicator method requires the lowest computational times. Nevertheless, the original method based on the behavior of the perturbation method close to eigenvalues seems to be a very promising technique.

Introduction

Vibroacoustic problems are encountered in many industrial applications (automotive, aerospace, ...). One of the difficulty to deal with this kind of problems is to well represent the coupling between fluid and solid. Several formulations can then be considered [START_REF] Morand | Fluid Strucutre Interaction[END_REF]. Among them, one can cite a pure displacementformulation, a three fields formulation (displacement, pressure and potential for example) or a displacement-pressure formulation. In this study, this latter is chosen. One of the difficulties for this displacement (solid) -pressure (fluid) formulation is that it leads to an unsymmetric eigenvalue problem. In order to solve this problem, some numerical methods can be used : they are listed and described in Ref. [START_REF] Voss | Structural-Acoustic Vibration Problems in the Presence of Strong Coupling[END_REF] for different cases of structures and fluid.

This paper concerns the free vibrations of an elastic structure which is completely filled with a homogeneous, inviscid and non-weighting compressible fluid. The displacementpressure equations are spatially discretized with the finite element method. The discrete linear non-symmetric eigenvalue problem is then solved with three numerical methods. One objective of this paper is to compare these numerical methods. The first algorithm is based on the common Arnoldi method with ARPACK solution [START_REF] Lehoucq | Arpack Users Guide: Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. For the two others methods, a perturbation method is used. The second method has been recently proposed in [START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF]. This method consists in introducing a right-hand side (r.h.s.) in the eigenvalue problem. This r.h.s. is a scalar multiplied by a random vector. This modified eigenvalue problem is then solved with the help of a perturbation method where the angular frequency is the perturbation parameter. Finding the couple of unknowns (eigenvalues and modes) consists then in finding the value of the angular * J.M. Cadou, jean-marc.cadou@univ-ubs.fr frequency for which the introduced scalar is null. In the last numerical method, the forced vibration problem is considered. A perturbation method is also used to compute the response curve of the considered problem. This third numerical method is based on the behavior of the series in the vicinity of a singularity. Then according to the work presented in Ref. [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF], close to a singularity (in this study the eigenvalue), a geometric progression appears in the asymptotic expansion. So determining the eigenvalues consists in finding the emergence of a geometric progression during the computation of the response curve of the forced vibration problem.

Governing equations

In this study, the displacement (u in the solid) and pressure (p in the fluid) formulation is chosen. The governing equations of an elastic solid (Ω s ) entirely filled with a homogeneous, inviscid and non-weighting compressible fluid (Ω f ) are the followings:

                 σ i j, j (u) + ω 2 ρ S u i = 0 in Ω s p + ω 2 c 2 p = 0 in Ω f σ i j (u)n f j = pn f j on ∂Ω sf ∂p ∂n f = p. n f = ω 2 ρ f u. n f on ∂Ω sf (1)
where subscripts s and f stand for solid and fluid domain respectively. The symbol ∂Ω sf is the solid-fluid interface where the couplings are. The scalars ω, ρ f , ρ s and c are the angular frequency, the density of the fluid and of the solid and the velocity of sound respectively. By using the finite element method, the discrete form of the previous equations is obtained:

K s -C 0 K f -λ M s 0 ρ f C t M f u p = 0 0 (2)
In the previous equations, the symbols K f or s , M f or s and C represent the stiffness matrix (for the fluid or the solid), the mass matrix (for the fluid or the solid) and the matrix due to the fluidsolid coupling. The scalar λ represents the eigenvalue (square of the angular frequency). For the sake of simplicity, the previous system of equations ( 2) is written as the following generalized eigenvalue problem:

[K -λM] U = 0 (3)
where K and M are unsymmetric matrices and U is a mixed unknown vector (the displacement and the pressure). This problem is now solved with the help of three different numerical methods.

Numerical methods

Classical numerical method

In this study, all the numerical results obtained with the proposed methods are compared to the results obtained with the Arnoldi method by using ARPACK solution [START_REF] Lehoucq | Arpack Users Guide: Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. As it has been pointed out in Ref. [START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF], the system (2) is poorly conditioned and then ARPACK is unable to compute the corresponding eigenvalues. So, according to the analysis reported in [START_REF] Maess | Substructuring and model reduction of pipe components interacting with acoustic fluids[END_REF], the system (2) needs to be preconditioned and becomes:

K s -C k 0 Kf -λ M s 0 ρ f C t m Mf u s p = 0 0 (4)
with the following definitions:

                     Kf = abK f and Mf = abM f C k = aC and C m = bC p = 1 a p a = K s F M f F and b = M s F K f F (5)
where symbol K F stands for the Frobenius norm of the matrix K defined by K F = n i=1 n j=1 |K i j |. With these modifications of the initial problem, ARPACK is able to compute the eigenvalues and eigenmodes of a vibroacoustic problem [START_REF] Maess | Substructuring and model reduction of pipe components interacting with acoustic fluids[END_REF]. These results are considered in the following as reference results.

The indicator method

The method proposed in Ref. [START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF] to compute the eigenvalues and eigenmodes of the system (3) is reminded. The key point of this method is to modify the Eq. ( 3) by introducing a right hand side µF:

[K -λM] U = µF ( 6
)
where F is a random load vector and µ is a scalar. Hence, finding the solutions (λ, U) verifying the eigenvalue problem (3) consists in determining when the scalar µ is null. In order to have a well-posed problem, the following additional equation is chosen:

U -U 0 , U 0 = 0 (7)
where •, • represents the Euclidian scalar product. The vector U 0 is the solution of the problem (6) when λ is equal to zero (i.e.

the static problem). In this case, the scalar µ is chosen equal to one. In order to solve the problem ( 6) with a perturbation method, the eigenvalue λ is defined by:

λ = λ 0 + λ (8)
where the value λ 0 is supposed to be known, equal to zero for the beginning of the calculation. The scalar λ is the new parameter of the problem. To determine accurate values of λ for which the scalar µ is null, the unknown X = (U, µ) are searched as an integro-power series with respect to the parameter λ:

X = N i=0 λi X i (9)
where N stands for the truncation order of the asymptotic expansions. By introducing these latter in the equations ( 6) and ( 7) and balancing terms with identical powers of λ, a set of linear problems is obtained:

                                     Order 0 : (K -λ 0 M) U 0 = µ 0 F U 0 , U 0 = U 0 , U 0 Order 1 ≤ i ≤ N : (K -λ 0 M) U i = µ i F + MU i-1 U i , U 0 = 0 (10)
Once, all these linear problems have been solved, the polynomial expansion ( 9) is replaced by a rational equivalent one, namely Padé approximants [START_REF] Baker | Padé approximants, Encyclopedia of Mathematics and its applications[END_REF][START_REF] Najah | A critical review of asymptotic numerical methods[END_REF]:

X Padé, N ( λ) -X 0 = N-1 k=1 R (N-1-k) ( λ) Q (N-1) ( λ) λk X k (11)
where R k , Q k are polynoms of degree k. Finally, finding the values for which the scalar µ is equal to zero consists in determining the roots, λr , of the polynoms R (N-1-k) ( λ). The accuracy of a given root is checked by the following criteriom:

t U r KU r -λt r U r MU r λr ≤ (12)
where is a chosen user tolerance. A continuation technique [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF] is used to determine all the eigenvalues between the initial value (λ 0 = 0) and the given maximum value of λ. Let us remark that this method is not sensitive at all to the poorly conditioned matrices [START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF].

It is finally pointed out, that this indicator method is quite similar to the numerical algorithm proposed by Avery et al. in Ref. [START_REF] Avery | A Padé-based factorization-free algorithm for identifying the eigenvalues missed by a generalized symmetric eigensolver[END_REF] to identify the eigenvalues missed by a symmetric eigensolver. Indeed, these authors proposed to define a scalar transfer function which is a rational function whose poles are exactly the eigenvalues of the considered problem. This transfer transfer is replaced by an equivalent Padé approximant, by computing Taylor expansions around chosen initial points. Analysing the poles of the Padé approximant permits to detect and identify the eigenvalues of the considered problem.

Analysis of the series

The last proposed numerical algorithm to solve the eigenproblem (3) considers the forced vibration problem:

[K -λM] U = F (13)
where F represents the amplitude of a time periodic load vector applied to the solid part of the considered problem : F(t) = F cos(ωt). The previous problem is also solved by using a perturbation method where the unknown X, with X = (U, λ), is defined by:

X = N i=0 a i X i ( 14
)
where the path parameter a can be written as:

a =< U -U 0 , U 1 > +(λ -λ 0 )λ 1 ( 15
)
By introducing the asymptotic expansions ( 14) into equations ( 13) and ( 15) and by equating like powers of a, a set of linear equations is obtained:

                                                                           Order 0 : (K -λ 0 M)U 0 = F with λ 0 known Order 1 : (K -λ 0 M)U = MU 0 λ 1 = 1 1+<U ,U > U 1 = λ 1 U Order 2 ≤ i ≤ N : (K -λ 0 M)Uq(i) = i-1 ir=1 λ ir MU (i-ir) λ i = -<Uq(i),U 1 > λ 1 +<U ,U 1 > U i = λ i U + Uq(i) (16)
So by solving the (N + 1) equations previously defined, one can compute a part of the solution. Next, the polynomial approximation ( 14) is replace by the Padé approximants [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. By applying the continuation technique as previously [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF], the response curve of a vibroacoustic problem submitted to a harmonic load vector is then computed in a range λ ∈ [0, λ max ] where λ max is chosen by the user.

The question is now how can be determined the eigenvalue and eigenmode solutions of Eq. 3 ? The proposed algorithm is based on the behavior of the asymptotic expansions close to a singularity. Indeed, in a recent paper, Cochelin and Medale [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF] have shown the emergence of a geometric power series in the asymptotic expansion near a singularity. In this latter reference, the singularity is a steady bifurcation. This technique has been used to detect steady bifurcation in Newtonian [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF][START_REF] Guevel | Numerical bifurcation analysis for three-dimensional sudden expansion fluid dynamic problem[END_REF] or non-Newtonian [START_REF] Jawadi | Asymptotic Numerical Method for steady flow of Power-law fluids[END_REF] fluid mechanics and Hopf bifurcation in rotating beam framework [START_REF] Bekhoucha | Eigensensitivity analysis of rotating clamped uniform beams with the asymptotic numerical method[END_REF]. From the work presented in Ref. [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF], the unknown X can be splitted into two parts:

X = X + X error ( 17
)
where X is a flawless solution and X error represents the solution due to the defaults or approximation errors in the computation. In the reference [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF], the authors have demonstrated that the term X error in the initial series X is a geometric progression series. So determining a singularity consists in finding a geometric progression into the polynomial representation [START_REF] Jawadi | Asymptotic Numerical Method for steady flow of Power-law fluids[END_REF]. Hence, this is carried out numerically by verifying a collinearity condition and a relative error test applied to the last terms of the series as follows:

For N -3 ≤ n ≤ N -1 compute :

β n = < X n , X N > < X N , X N > and X ⊥ n = X n -β n X N (18)
If

N-2 n=N-3 |β n | 1/(N-n) |β N-1 | -1 2 ≤ δ 1 ( 19
)
and

N-1 n=N-3 X ⊥ n X h ≤ δ 2 ( 20
)
where the two parameters δ 1 and δ 2 are chosen equal to respectively 10 -3 and 10 -6 according to the numerical results presented in Ref. [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF].

If both previous criteria are checked, then the common ratio of the geometric progression (β N-1 ) is exactly the distance to the singularity, denoted by a 0 , and leads to the computation of all the terms Xn of the flawless solution:

For 1 ≤ n ≤ N -1 : Xn = X n -a (N-n) 0 X N (21)
These terms are then used to compute a "cleaned" solution. Finally, once a geometric progression has been numerically detected, the eigenvalue and the eigenmode can be computed by using the following expression:

X c = X 0 + N-1 i=1 a i 0 Xi (22)
Moreover, it is known from [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF] that the vector X N is related to tangents at singular solutions. So, when the continuation is performed on the fundamental branch and in the case of a pitchfork bifurcation, X N is exactly the direct mode at the singluar solution, see Ref. [START_REF] Guevel | Numerical bifurcation analysis for three-dimensional sudden expansion fluid dynamic problem[END_REF]. Using the Lyapunov-Schmidt reduction [START_REF] Guevel | Numerical bifurcation analysis for three-dimensional sudden expansion fluid dynamic problem[END_REF], one shows that this direct mode is the orthogonal tangent to the fundamental branch.

In the considered vibroacoustic problem, the flawless solution (i.e. X) is the solution of the generalized eigenvalue problem (3) and the "error" solution (solution with default) is due to the load applied to the structure (Eq. 13). So, close to the eigenvalue, a geometric progression should appear in the asymptotic expansions [START_REF] Jawadi | Asymptotic Numerical Method for steady flow of Power-law fluids[END_REF].

Numerical results

The three previous numerical methods are now applied to determine the eigenvalues and eigenmodes of a steel cavity (linear elastic structure) which is filled with an inviscid compressible fluid (air or water) (see Figure 1). This problem has been solved in references [START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF] and [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. Linear quadrilateral elements are used both for solid and fluid domains. air and water obtained with ARPACK and the indicator method are given in Table 1 and compared to those from Ref. [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. The Air Water PMs

Reference [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF] PMs

Table 1: Comparison of the first ten angular frequencies obtained with the proposed methods (PMs) and those given in Ref. [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF] for the elastic cavity presented in Fig. 1.

solutions obtained with the indicator method (Section 3.2) and those obtained with ARPACK are exactly the same. Some of the eigenmodes (Mode 1, 2 and 4) are represented in Fig. 2 for a cavity filled with air. One considers now the forced vibration problem of an elastic cavity filled with water. The response curve is given in Fig. 3. This curve is obtained by choosing a truncation order N of the asymptotic expansion ( 14) equal to 20. More than one hundred steps of the continuation technique is required to compute the solution up to an angular frequency equal to 3000 rad/s. The interesting feature is the step accumulation near the eigenfrequencies, ω 1 =654 rad/s and ω 2 = 2159 rad/s in Fig. 3. Each step of the continuation technique is represented by a cross symbol in Fig. 3. This behavior of the perturbation method close to singular values is classical and has been used in some applications as a visual indicator to detect instabilities. The algorithm presented in Section 3.3 and based on the detection of a geometric progression is now applied to detect eigenvalues by analysing the series. Then, when the latter method is used, only 55 steps of the continuation technique are needed to computed the solution for an angular frequency in the range [0, 8000]. This permits to compute the first ten eigenvalues (for the water case) given in Table 1. These eigenvalues are exactly the same as those computed with ARPACK or with the indicator method. With the analysis technique of series, the eigenmodes are easily carried out by introducing the distance to the closest singularity (a 0 ) into the expression (22). These modes are plotted in Fig. 4 for the air-filled cavity. In this figure, one can clearly see that these modes are not the same that those obtained with ARPACK or with the indicator method (see Fig 2). For a more precise idea of the accuracy of the modes computed by using expression (22), the logarithm of the residual of the mode 2 has been calculated by using Eq. ( 12) and is equal to 7.83. This confirms that the modes obtained with expression (22) are not solution of the initial eigenvalues problem (3). This is quite surprising compared to the previous applications of this technique for computing steady bifurcations for example. Nevertheless, there exists another way to compute these eigenmodes. One considers for example a computation (the forced vibration case) carried out near the second eigenvalue of the air-filled cavity (close to ω 0 = 1000 rad/s with λ 0 =ω 2 0 in Eq. ( 16)). For this value of the angular frequency, a geometric progression has been detected and the next eigenvalue is then found. So, for this computation, some of the terms of the asymptotic expansion (the component U i of the mixed vector X i of Eq. 14) are plotted in Figures 5 and6 for several truncation orders. These figures show that the component U i (displacement and pressure) becomes similar to the mode 2 (see Figures 2(c) and 2(d)) when the truncation order increases. To underline this behavior, the evolution of the residual (12) of the vector U i is plotted versus the evolution of the truncation order in Fig. 7. This figure shows that the vector U i is exactly the eigenmode when the truncation order increases. The accuracy of the vector U i seems to linearly depend on the truncation order. With this property, accurate eigenmodes can be then computed instead of using expression (22). In a way, Dirac delta functions in the response curve (Fig. 3) act as second famillies of solution crossing the fundamental branch at the singular solutions. As those functions are orthogonal to the curves, it is consistent to Figure 7: Evolution of the logarithm of the residual vector (Eq. 12) versus the truncation order for the vector U i close to the second eigenvalue for an air-filled cavity.

see the eigenmode in the high order terms of the series. The last point concerns the computational times required for each presented technique. So the computation times needed to 1) for the elastic cavity filled with air [START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF][START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. 1) for the elastic cavity filled with water [START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF][START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF].

compute the first ten eigenvalues are plotted versus the number of unknowns considered in Figures 8 and9 for air and water respectively. These figures show that the algorithm based on the analysis of the series is less time consuming than the ARPACK method. The indicator technique is the method requiring the lowest computational times. The numerical tests (not presented in this work) have shown that, on the contrary of the indicator technique, the analysis of series method is sensitive to the fact that matrices are ill-conditioned. So for this latter, as with ARPACK in section (3.1), all the matrices have been modified according to the expressions (5).

Conclusion

In this paper, three numerical techniques to compute the eigenvalues of a vibroacoustic interior coupled problem have been compared. A classical technique using the Arnoldi method with ARPACK solution is compared to two techniques using a perturbation method. A first method which has been initially proposed in [START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF] consists in determining the null values of an indicator. A second technique, proposed in a forced vibration framework, is based on the analysis of the behavior of the series close to singularities. Thanks to the detection of a geometric progression in the series, accurate eigenvalues can be determined. When a geometric progression is found, the vibration modes are the terms of the Taylor expansions at a high order (around 15 according to the presented results). The computa-tional times are lowest with the indicator technique. Nevertheless, the method using the analysis of the series close to the eigenvalues seems to be a promising technique. Some works are in progress to improve this technique, for example by extracting singularities in the vicinity of the eigenvalues [START_REF] Soliman | Extended perturbation series analysis and improvement applied to a porous catalyst slab[END_REF] or by using Euler transform [START_REF] Bromwich | An introduction to the theory of infinite series[END_REF]. All these techniques are based on Domb-Sykes plot [START_REF] Van Dyke | Analysis and improvement of perturbation series[END_REF][START_REF] Domb | On the susceptibility of a ferromagnetic above the Curie point[END_REF] which easily leads to the closest singularity of the initial point of the perturbation method.

3 FFigure 1 :

 31 Figure 1: Geometry and mechanical properties for the steel cavity [11].

Figure 2 :

 2 Figure 2: Eigenmodes for the elastic cavity filled with air (from reference [4]).

Figure 3 :Figure 4 :

 34 Figure 3: Response curve of the forced vibration problem for the elastic cavity filled with water. The displacement corresponds to the point where the load is applied (see Fig. 1).

Figure 5 :

 5 Figure 5: Terms of the asymptotic expansions close to the second eigenvalue for an air-filled cavity according to the truncation order.

Figure 6 :

 6 Figure 6: Terms of the asymptotic expansions close to the second eigenvalue for an air-filled cavity according to the truncation order.

Figure 8 :

 8 Figure 8: Comparison of CPU times obtained with ARPACK and the proposed methods to get the ten first eigenvalues (Table1) for the elastic cavity filled with air[START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF][START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF].

Figure 9 :

 9 Figure 9: Comparison of CPU times obtained with ARPACK and the proposed methods to get the ten first eigenvalues (Table1) for the elastic cavity filled with water[START_REF] Claude | Eigensolutions of a vibroacoustic interior coupled problem with a perturbation method[END_REF][START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF].


J M Cadou
email: jean-marc.cadou@univ-ubs.fr

Keywords: Linear vibroacoustic problem, eigenvalues, pertubation method

In this paper, an efficient and robust numerical method is proposed to solve non-symmetric eigenvalue problems resulting from the spatial dicretization with the finite element method of a vibroacoustic interior problem. The proposed method relies on a perturbation method. Find the eigenvalues consists in determining zero values of a scalar which depends on the angular frequency. Numerical tests show that the proposed method is not sensitive to poorly conditioned matrices resulting from the displacement-pressure formulation. Moreover, computational times required with this method are lower than those needed with a classical technique such as for example the Arnoldi method.

Introduction

The study of linear vibrations of fluid-solid interaction (FSI) problems is encountered in many industrial applications. In this paper, the considered FSI problem consists on an elastic solid filled with an inviscid and compressible non-weighting fluid. A variational principle based on a displacement(for the solid)-pressure (for the fluid) formulation is adopted. In this approach, spatially discretized equations lead to a non-symmetric eigenvalue problem. Moreover, this latter can be poorly conditioned mainly due to large difference of magnitude between the fluid and the solid stiffness and mass matrices [START_REF] Maess | Substructuring and model reduction of pipe components interacting with acoustic fluids[END_REF]. This poor conditioning can harm the convergence of classical eigensolvers and sometimes gives inaccurate results. So, in this paper, a very basic numerical method is proposed to solve the non-symmetric and poorly conditioned eigenvalue problem associated with an interior vibroacoustic problem. This numerical method is based on a pertubation method and has been initially proposed in Ref. [START_REF]Etude des vibrations linéaires de plaques par une méthode asymptotique numérique et les approximants de Padé[END_REF] to solve symmetric eigenvalue problems arising from elastic linear solid framework. The method consists in introducing a right-hand side (rhs) in the initial problem. This rhs is a scalar multiplied by a random vector. The objective is then to find the angular frequency value for which this scalar becomes null. Indeed, this null value indicates a solution of the initial problem (i.e an eigenvalue). To find the null values of this scalar, unknowns of the problem (solid displacement, fluid pressure and the scalar) are searched as asymptotic expansions of the angular frequency. Resolution of a set of linear problems permits to build an analytical expression whose roots can be numerically computed. From these roots, the eigenvalues of the vibroacoustic problem are determined.

The present paper is organized as follows, the discretized non-symmetric eigenvalue problem to be solved is introduced in Section 2. The section 3 is devoted to the presentation of the proposed numerical algorithm. In section 4, numerical tests, a deformable cavity filled with air or water, permit to show efficiency and accuracy of the proposed method.

Governing equations

In this study, the problem of an inviscid compressible fluid contained in an elastic solid is considered. Using the finite element method to discretize the variational <u-p> formulation, it is obtained the following classical matrix system:

K s -C 0 K f -λ M s 0 ρ f C t M f u s p = 0 0 (1)
In this equation u s and p are the solid displacement and the fluid pressure respectively. The subscripts s and f stand for the solid and the fluid. The matrices K, M and C represent the stiffness matrix, the mass matrix and the coupling matrix due to the fluid-solid interaction. The scalar ρ f is the fluid density. Eigenvalue λ is the square of the angular frequency (i.e. λ = ω 2 ). Formally, the previous equation ( 1) can be written as the following generalized eigenvalue problem:

(K -λM) U = 0 ( 2
)
where U is a mixed unknown vector containing the displacement of the solid, u s , and the pressure, p, in the fluid. The previous system is non-symmetric due to the presence of the coupling matrix, C, in the mass and stiffness matrices, respectively M and K. Moreover, as the previous system is poorly conditioned (see [START_REF] Maess | Substructuring and model reduction of pipe components interacting with acoustic fluids[END_REF] for more details), numerical preconditioning is strongly encouraged to ensure convergence and accuracy of the solutions. In this study, with a classical eigensolver, ARPACK [START_REF] Lehoucq | Arpack Users Guide: Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF], a preconditioning technique should be added to find eigensolutions. On the contrary, with the proposed method based on perturbation technique, no preconditioning will be required to find accurate solutions.

The proposed method

The numerical method designed to solve linear problem (2) has been initially proposed to solve eigenvalue problems arising from vibrations of elastic solids [START_REF]Etude des vibrations linéaires de plaques par une méthode asymptotique numérique et les approximants de Padé[END_REF]. So to solve the problem (2), this latter is modified by introducing a right-hand side µF where µ is an unknown scalar and F is a random load vector. The problem (2) then becomes:

(K -λM) U = µF (3)
As the number of unknowns is greater than the number of equations, the following orthogonality condition is chosen to be the additional equation:

U -V, V = 0 (4)
where •, • stands for the Euclidian scalar product of two vectors. The vector V is a known vector and will be precised below. Hence, according to the Eq. ( 3), determine solutions (λ, V ) which satisfy the initial problem (2) consists to determine the values λ for which the scalar µ is equal to zero. A way to find these null values of µ is to compute the scalar µ for several discrete values of λ. This 'direct' method requires a lot of computations (i.e. a lot of matrix triangulations K) and then is not efficient in term of computational times. Instead of this 'direct' method, it is proposed to solve the problem (3) by using a perturbation method, λ being the perturbation parameter. At the end of the computations, an analytical expression of the scalar µ is then obtained. The determination of the value of λ for which the scalar µ is equal to zero should be then straightforward. In order to use the perturbation method, it is suggested to define the eigenvalue λ under the following form:

λ = λ 0 + λ (5)
where λ 0 is the initial value. For the first step of the presented method, λ 0 is chosen equal to zero. Problem (3) can be rewritten:

K -(λ 0 + λ)M U = µF (6)
The unknowns (U, µ) are then sought as an integro-power series with respect to the parameter λ:

U µ = N i=0 λi U i µ i ( 7
)
where N is the truncation order of the asymptotic expansions. Inserting Eq. ( 7) in equations ( 4) and ( 6) and balancing terms with identical powers of λ a set of linear problems is obtained: Order 0:

(K -λ 0 M) U 0 = µ 0 F U 0 , V = V, V (8)
Order 1 ≤ i ≤ N:

(K -λ 0 M) U i = µ i F + MU i-1 U i , V = 0 (9)
All linear problems (8) and ( 9) have the same operator and differ from their right-hand sides. So one single matrix triangulation and (N+1) backward and forward substitutions are required to computed all the unknowns (U i , µ i ) of the polynomial approximation [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. At the end of the computation, asymptotic expansions are replaced by equivalent rational approximations, called Padé approximants [START_REF] Baker | Padé approximants, Encyclopedia of Mathematics and its applications[END_REF][START_REF] Najah | A critical review of asymptotic numerical methods[END_REF] X Padé, N ( λ) -

X 0 = N -1 k=1 R (N -1-k) ( λ) Q (N -1) ( λ) λk X k ( 10
)
where R k , Q k are polynoms of degree k and X is a mixed vector containing the vector V and the scalar µ. This rational representation permits to increase the validity range of the polynomial approximation. At the end of the computations, finding the couples (U, λ) which are solutions of the generalized eigenproblem (3), consists in determining the roots of numerators R (N -1-k) in rational approximations [START_REF] Boumediene | Nonlinear forced vibration of damped plates coupling asymptotic numerical method and reduction models[END_REF]. These numerators are defined by the following expression:

R (J) = 1 + J j=1 λj d j for 0 ≤ J ≤ N -2 (11)
where scalar coefficients d j are expressed as:

d 1 = - α (N,N -1) α (N -1,N -1)
and

d j = - α (N,N -j) α (N -j,N -j) - j-1 l=1 α (N -l,N -j) α (N -j,N -j) d l for 2 ≤ j ≤ N -1 (12)
The coefficients α (I,J) are determined by performing a modified Gram-Schmidt orthonormalization of the initial set of vectors X N . The denominator Q (N -1) ( λ), which is the same for the rational fraction (Eq. 10) is computed in the same way as the numerators R (N -1-k) by using expressions [START_REF] Bobillot | Iterative techniques for eigenvalue solutions of damped structures coupled with fluids[END_REF] and [START_REF] Deü | Vibration and transient response of structural-acoustic interior coupled systems with dissipative interface[END_REF].

The roots of the numerators, R (N -1-k) , are denoted in the following by (λ r ). Among all these roots, some of them do not lead to a solution of the initial problem (3). To check if a root, λ r , gives a true and accurate eigenvalue, the following criterion is defined:

( t U r KU r -λ t r U r M U r ) λ r ≤ ǫ (13)
where U r is computed by introducing the root λ r into Padé approximant [START_REF] Boumediene | Nonlinear forced vibration of damped plates coupling asymptotic numerical method and reduction models[END_REF]. The parameter ǫ is an user parameter which specifies the quality of the numerical solutions. In this study, the value of ǫ is chosen equal to 10 -8 . Once all the roots have been carried out for one step (i.e. for a given value of λ 0 ), a new step of the method is done by defining a new value of the scalar λ 0 . In this work, this latter is computed in two ways. If during the previous step, a solution of the eigenvalue problem has been found, the root λ r verifying the criterion (13), a new starting value λ 0 is defined according to the following expression:

λ 0 = 1.1 * λ r (14)
Otherwise, if no eigenvalue has been found, then the range of validity of the Padé approximants is estimated by using a criterion introduced in Ref. [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF] and leading to a new initial value λ 0 . Lastly, the method requires an additional condition (4) in which a vector V must be introduced. The first step of the current algorithm permits do define this vector. To get it, λ 0 is considered to be null and the initial value of the scalar µ is arbitrary chosen equal to one (µ = 1). Then, from equation (3), vector V is obtained as solution of:

KV = F (15)

Numerical results

To prove the efficiency of the proposed method, a steel cavity (linear elastic structure) is filled with an inviscid compressible fluid (air or water). This example is issued from Ref. [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. Geometric and material characteristics are given in

Ω f Ω s Γ s Γ f
Figure 1: Geometry and mechanical properties for the steel cavity [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF] elements both for the fluid and the structure. Solid element has eight degrees of freedom (dof) corresponding to two displacements and four nodes. Fluid element has four dof (pressure at each node of the fluid element). Finally, interface element has twelve dof (eight for the solid displacement and four for the fluid pressure). In this study, Fig. 1 interface). Then the number of dof for this mesh is close to 4000. The first ten eigenvalues obtained with the proposed method are reported in Table 1 and compared to values given in the work of Bermúdez and Rodríguez [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF].

Results in this table show an excellent agreement. In Fig. 2, fluid and solid modes corresponding to the first four eigenfrequencies given in Table 1 are presented. For the proposed method, the results given in Table 1 are obtained with a truncation order, parameter N in Eq. ( 7), equal to 15 and an accuracy parameter, ǫ in Eq. ( 13), chosen equal to 10 -8 . For these parameters, 37 continuation steps are required to compute the first ten eigenfrequencies for cavity filled with air and 43 steps for cavity filled with water. In Fig. 3, it is plotted the evolution of the number of steps versus the truncation order to get the first ten eigenvalues (Table 1). Truncation order varies between 10 and 35 whereas accuracy parameter is fixed to 10 -8 . This plot shows that the number of continuation steps is relatively constant regardless of the truncation order chosen. In fact the key point of the proposed algorithm is the roots computation of numerator R k of the Padé approximants [START_REF] Boumediene | Nonlinear forced vibration of damped plates coupling asymptotic numerical method and reduction models[END_REF]. So, to get accurate values of these roots, two numerical methods are used [START_REF] Press | Numerical Recipes in Fortran 77, The Art of Scientific Computing[END_REF]. The first one is the Bairstow's method and the second one consists in finding the eigenvalues of a matrix whose characteristic polynomial is the numerator of the Padé approximants [START_REF] Boumediene | Nonlinear forced vibration of damped plates coupling asymptotic numerical method and reduction models[END_REF]. With these two techniques, accurate roots are found whatever truncation order varying between 10 and 35. Neverthe-less, truncation order has an influence on the computational time required to obtain eigenvalues. Indeed, a truncation order equal to N requires a single matrix triangulation but (N+1) backward and forward substitutions to compute the (N+1) couples of unknowns (U i , µ i ). So in Fig. 4, it is plotted the evolution of the computational time to get the ten first eigenvalues (given in Table 1) when truncation order is varying between 10 and 35. This plot shows that the lowest computational times is obtained for a truncation order close to 15. Performance of the present algorithm is now compared to a classical numerical method to compute eigenvalues of a non-symmetric problem. In this study the Arnoldi method is chosen by using the ARPACK solution [START_REF] Lehoucq | Arpack Users Guide: Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. Unfortunately, ARPACK is unable to compute eigenvalues of the problem (2). Indeed, this problem is poorly conditioned due to the fact that components of matrices M s and ρ f C t are very large compared to entries in matrix M f [START_REF] Maess | Substructuring and model reduction of pipe components interacting with acoustic fluids[END_REF]. So, the generalized eigenvalue problem (Eq. 1) needs to be preconditioned to give fast and accurate eigenvalues. According to the analysis reported in Ref. [START_REF] Maess | Substructuring and model reduction of pipe components interacting with acoustic fluids[END_REF], the preconditioned generalized eigenvalue problem to be solved with ARPACK is rewritten:

K s -C k 0 Kf -λ M s 0 ρ f C t m Mf u s p = 0 0 (16)
with the following definitions:

           Kf = abK f and Mf = abM f C k = aC and C m = bC p = 1 a p a = Ks F M f F and b = Ms F K f F (17)
where symbol K F stands for the Frobenius norm of the matrix K defined by

K F = n i=1 n j=1 |K ij |.
A key feature of this study is that the proposed method gives accurate eigenvalues with or without preconditioned matrices. So in Fig. 5, computational times required to get the first ten eigenvalues of the numerical tests (Table 1) are plotted versus the number of unknowns (from 500 to 4000). These times are given for both numerical methods (ARPACK and the proposed numerical method). This figure shows that the proposed algorithm needs computational times lower than those required with ARPACK whatever the size of the problem. Hence, with the finest mesh and in the case of a cavity filled with air, the proposed algorithm demands a computational time four times lower than the ARPACK computational time.

Conclusion

In this paper, it is proposed to use a perturbation method to compute the eigenvalues of a vibroacoustic interior coupled problem. Eigenvalues are determined in analyzing roots of the numerator of a fraction (Padé approximants). The proposed method is easy to implement in a computational software because it only requires a linear solver and a subroutine to realize matrix-vector products. This algorithm is efficient even if two eigenvalues are very close each other which is the case when considering an elastic cavity filled with air. The proposed numerical method does not require use of a preconditioning technique contrary to classical eigensolver such as Arnoldi method. The optimal truncation order for the polynomial approximation is close to 15 according to the analysis carried out in this study. Computational times needed with the proposed numerical method are lower than the ones required with the classical ARPACK solver. Nevertheless, the number of steps to get ten eigenvalues is relatively important (approximately 40 for the considered examples) leading to a great number of matrix triangulations. So, for large scale problem, involving several thousand of dof, the presented algorithm can require too large computational times. Future works concern the use of reduced order models based on the perturbation method (recently proposed to analyze fluid bifurcation problems [START_REF] Heyman | Computation of Hopf bifurcations coupling reduced order models and the asymptotic numerical method[END_REF] and vibrations of plates [START_REF] Boumediene | Nonlinear forced vibration of damped plates coupling asymptotic numerical method and reduction models[END_REF]) to decrease computational times. Application of this kind of method to analyze vibrations of damped structures [START_REF] Bobillot | Iterative techniques for eigenvalue solutions of damped structures coupled with fluids[END_REF][START_REF] Deü | Vibration and transient response of structural-acoustic interior coupled systems with dissipative interface[END_REF] are also in progress. 1). Elastic cavity filled with air or water [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. Parameters used for the proposed method are N=15 and ǫ = 10 -8 .

Fig. 1 .

 1 Finite element discretization uses linear quadrilateral Steel: E=1.44 10 11 Pa ν= 0.35 ρ s = 7.7 10 3 kg.m -3 Air's properties: c f = 340 m.s -1 ρ f = 1 kg.m -3 Water's properties c f = 1430 m.s -1 ρ f = 1000 kg.m -3

Figure 2 :

 2 Figure 2: First four modes of the coupled system, fluid (air) pressure mode (a, c, e and g) and solid displacement mode (b, d, f and h).11

Figure 3 :

 3 Figure 3: Number of continuation steps of the proposed method versus truncation order to get the first ten eigenvalues shown in Table1. Elastic cavity filled with air or water, 4000 dof.

Figure 4 :

 4 Figure4: Evolution of the CPU times (s) versus truncation order to get the first ten eigenvalues shown in Table1. Elastic cavity filled with air or water, 4000 dof. Accuracy parameter ǫ is equal to 10 -8 .

Figure 5 :

 5 Figure5: Comparison of CPU times obtained with ARPACK and the proposed method to get the ten first eigenvalues (Table1). Elastic cavity filled with air or water[START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. Parameters used for the proposed method are N=15 and ǫ = 10 -8 .

Table 1 :

 1 Comparison of the first ten eigenvalues obtained with the proposed method (PM) and those given in Ref.[START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF].

	, the finest mesh is composed of 2340


Pierre Charrier

Elisabeth Ostoja-Kuczynski

Erwan Verron

Laurent Gornet

Grégory Chagnon

INFLUENCE OF LOADING CONDITIONS ON FATIGUE PROPERTIES FOR FILLED ELASTOMERS

Fatigue characterisation of elastomers is an essential prerequisite of new car project schedules imposed by carmakers. Consequently, durability prediction softwares are needed to reduce the time allowed to design elastomeric parts. Until now, the only theoretical variable used to predict duration life is the tearing energy proposed by Rivlin and Thomas. Classically, fatigue properties are determined using plane stress experiments (trousers or pure shear samples for example) under simple loading conditions.

In the present paper, a brief bibliographical review presents the limitations of this approach and the need of more complex loading conditions is emphasised. Then, a new simple shear sample is proposed to simplify the treatment of experimental results, said the calculation of the tearing energy, under complex loading conditions. The test sample is initially notched and is submitted to controlled displacement loading conditions. The procedure adopted to analyse the results is presented. Assuming that the tearing energy and the crack growth rate remain constant during the test, they can be easily calculated using both the decreasing stiffness of the sample during the test and a FE model. Note that no complex system, such as video recording, is necessary to measure the crack length. Then, experimental results are reported. Both filled crystallising and non-crystallising elastomers are studied in order to highlight the influence of the testing temperature and of the imposed displacement ratio (i.e.: minimum/maximum displacement) on the fatigue crack growth rate law. This fatigue crack growth rate law obtained agrees well with previously published papers. Furthermore, other complex tests are in progress to determine the influence of the Mullins stress-softening on fatigue properties. Finally, in regards with these results, it is demonstrated that the tearing energy is not a sufficient variable to completely characterise the fatigue life of elastomers.

INTRODUCTION

Tires manufacturers have studied the durability of elastomers for more than 60 years. Numerous scientific papers induced the definition of a well-adapted criterion for this type of industrial applications: the tearing energy T. Moreover, a large experimental data bank is available, but most experiments are restricted to plane stress test samples under simple loading conditions. Automotive anti-vibration system (AVS) companies have now a great interest in elastomers durability. However, the corresponding industrial components and loading conditions are greatly different than those involved in tires industry. Indeed, the loading conditions are close to the state of plane strain and loading histories are more complex: they include variable pre-loading, non-relaxing conditions… In order to investigate the durability of elastomers in AVS applications, it is necessary to determine the limitations of the tearing energy criterion to ensure that its use for AVS products is well-adapted. Second, new test sample geometries have to be defined to reproduce the local stress state encountered in these applications.

In the first part of the present paper, limitations of the tearing energy T are investigated. More precisely, the assumptions that were made by Rivlin and Thomas to define T, the effect of pre-loadings (maintained or not) and the limitations due to the elastic nature of T are thoroughly examined. In a second part, a new test sample geometry is presented. This new sample provides plane strain conditions and could be used to study the influence of loading history, and especially of non-relaxing conditions. The corresponding experimental methodology, said the determination of T and the measurement of the crack propagation, is described in details. More precisely, although the crack growth measurement is easily performed for thin rubber sheets [GEN64], in the present case of a massive rubber block, it is more difficult to determine the crack length. So, a new method for the crack growth measurement is proposed. It is based on the determination of the global stiffness of the sample. Note that the present experimental approach does not need any video equipment to measure the crack propagation. Finally, first experimental results under different temperatures and pre-loading conditions are presented. Two filled materials are considered: a NR that exhibits strain crystallisation and a PU that does not

THEORETICAL FRAMEWORK

DEFINITION OF THE TEARING ENERGY

In order to study the fracture mechanics of rubber-like materials, the first step consists in defining a failure criterion.

Opposite to metallic materials, present geometric and material non-linearities do not permit the use of simple criteria based on the stresses, because the stress state can not be simply calculated analytically. Rivlin and Thomas [RIV53] proposed an energetic failure criterion: the tearing energy T. This approach can be seen as the extension of Griffith's theory [GRI20] to rubber-like materials. The corresponding crack propagation criterion can be written as:

( ) c ext T A W W T ≥ ∂ - ∂ =
, where T c is the critical value of the tearing energy at the failure point, W is the stored elastic energy, W ext represents the work of external forces and A stands for the crack area. Under constant displacement loading conditions, the work of external forces is equal to zero and the previous equation becomes:

c l T c W t T ≥       ∂ ∂ - = 1
in which t is the thickness of the sample and c is the crack length. This definition of the tearing energy is widely used in the bibliography. It is to note that the Griffith's criterion, that is the basis of the definition of T, considers that the decrease of the stored elastic energy during the crack growth is only balanced by the energy dissipated in the crack front during propagation. Nevertheless, in the case of elastomeric materials, the decrease of the stored elastic energy is balanced by both the energy dissipated in the new crack area and by other dissipated energies that are due to the material behaviour (e.g. Mullins effect, viscoelasticity, strain-crystallisation).

Consider now the tearing of thin sheets and let us examine the assumptions of Rivlin and Thomas in order to determine if it is possible to consider the previous definition of T relevant in that special case. According to the authors, the energy dissipated in the volume of the sample, far from the crack front, can be neglected in the case of low dissipative materials. This is not the case of dissipation phenomena that occur in the neighbourhood of the crack front. In this area, the strain magnitude is higher and the stress state is more complex (anisotropic, viscoelastic, compressible…) than in the rest of the sample, said far from the crack. Rivlin and Thomas assume that these locally dissipated energies remain proportional to the crack propagation dc and to geometric characteristics of the crack front. They concluded that the energy required to propagate the crack of a length dc does not substantially depend on the sample geometry and on the way loading conditions are applied. Consequently, this energy, classically measured during experiments, only characterises the tearing of thin sheets and includes both the tearing energy T and the energy dissipated in the crack front.

CLASSICAL TEST SAMPLES

In order to study the crack propagation in elastomeric materials, Rivlin and Thomas proposed three different geometries for the test samples. Every samples are thin sheets in plane stress conditions close to the crack front. These three samples are the tensile test sample, the trouser sample and the pure shear sample; they are presented in Figure 1. They have been widely used in the papers dealing with the durability of elastomers.

LIMITS TO THE USE OF T

Before using the tearing energy for samples that are different than thin sheets, we have to review the experimental and theoretical limits of the tearing energy defined above, as already highlighted in literature.

First, strictly regarding to the definition of T, the experimental procedures proposed in literature are not well-adapted to its measurement. Indeed, during experiments, it is not possible to separate the energy used for the propagation of the crack area, said the tearing energy, from other energies that are dissipated both in the volume and in the neighbourhood of the crack front. Only a global value of the dissipated energy can be obtained. Moreover, the tearing energy as defined by Rivlin and Thomas is a purely elastic parameter. It does not take into account the complex behaviour of elastomers that includes strain-softening, viscoelasticity ... Gent and Kim proposed a new experimental definition of T that partially overcomes this difficulty [GEN78]. To take into consideration the hysteresis of the material under cyclic loading, they determine the tearing energy using the elastic energy restored during the unloading part of the cycle and not the elastic energy that is stored during the loading part. Finally, the bibliography exhibits different experimental results which prove that T is not an intrinsic parameter, said it can not be used as a criterion. This is the case of Thomas who demonstrates that the critical tearing energy T c is modified by cyclic loadings of the test sample [THO58]: the material is strengthened during this cyclic loading. The works of Gent and Kim under static conditions [GEN78] and of Busfield et al. under cyclic conditions [BUS96] highlight the effect of pre-loading on the tearing energy. In the case of pure shear test samples, the authors show that a pre-loading parallel to the crack weakens the material, said the corresponding critical tearing energy is divided by 25. These results demonstrate that the tearing energy is not an appropriate parameter to characterise the durability properties of elastomeric materials.

DEFINITION OF A NEW TEST SAMPLE

GEOMETRY OF THE SAMPLE

As mentioned above, classical test samples used to study the crack propagation in elastomeric materials are considered under plane stress conditions. So, they can not be used to study the local stress state encountered in AVS components. In order to reproduce the plane strain conditions, we should consider test samples under plane strain conditions. This consideration has led us to a new design of samples.

Lindley et al. ([LIN77,LIN79,LIN82]

) provide the methods needed to calculate T for elastomeric blocks in which a crack propagates at the interface between the material and a rigid insert. Using their results, we consider a "2 blocks -simple shear sample" with a pre-cut at the interface elastomer/insert. Figure 2 presents this sample. As proposed in the previous cited papers, both the block width and the crack length are at least twice larger than the block thickness. Moreover, using a well-defined pre-cut length, our new sample can be used for crack propagation experiments in which the tearing energy and the crack growth rate remain constant when the crack propagates.

DETERMINATION OF THE TEARING ENERGY

The determination of the tearing energy that corresponds to an elastomeric block under enforced displacement is directly based on the definition of Rivlin and Thomas:

l elas c W b T       ∂ ∂ - = 1
, where b is the block width and W elast is the elastic energy stored in the sample. This elastic energy is calculated assuming that the material behaviour is linear under simple shear loading. Denoting L max and the maximum enforced displacement and resulting force, respectively , the maximum tearing energy obtained during the test is denoted T max and can be written as: ). In the present work, this global stiffness is determined by a finite element analysis. The rubber block is modelled under plane strain conditions and simple shear experiments are simulated using the software ABAQUS. For a given value of the crack length increment, said dc, the crack propagation effect on the global stiffness is determined by performing computations on samples in which the crack length value is set to c, c+dc, c+2dc, …. For the simulations, the material behaviour is modelled by the nearly incompressible hyperelastic Mooney-Rivlin constitutive equation.

The results of these simulations can be summarised by the curve of the resulting force versus the crack length. More precisely, we consider non-dimensional variables, said the non-dimensional resulting force F/F 0 , where F 0 is the resulting force obtained for the block without pre-cut, and the non-dimensional crack length c/ l c where l c is the sample length. The corresponding curve is presented in Figure 3. As shown in this figure, the curve can be fitted by a linear equation for sufficiently large values of the crack length. The fitting results are shown in Figure 3. The slope of the line is denoted A and given by: Figure 4 : linear fit.

c F F l A c ∂ ∂ - = max 0 c/lc α = F(c)/Fo c/lc α = F(c)/Fo ( ) B l c A F c F c +         - = 0
Using this new parameter, the maximum tearing energy can be written as:

c l b F L A T ⋅ ⋅ ⋅ = 2 0 max max
This equation expresses that T max is a function of experimental parameters and of the resulting force of a block without pre-cut F 0 . As this force can be related to the resulting force F C0 (corresponding to a pre-cut block) by a constant parameter 1/α as shown in Figure 3, the maximum tearing energy only depends on experimental parameters:

c c l b F L A T ⋅ ⋅ ⋅ = 2 0 max max α
where b and l c describe the geometry of the sample, L max is the enforced displacement and F c0 is the resulting force corresponding to the accommodated sample without crack propagation.

The method adopted for the determination of T max presents a major limitation. The pre-cut length and the crack length considered during the test should be chosen in order to ensure that the parameters A and α remain constant during the cyclic experiments.

MEASUREMENT OF THE CRACK GROWTH RATE

In the bibliography, the crack growth rate is often measured using a video recording system that gives ∆c every ∆n cycles. This method is well-adapted to experiments under plane stress conditions, but it is no more suitable for our plane strain samples, because of edge effects. Consider a massive elastomeric block, the crack front shape is not straight because it is deviated by edge effects: the crack growth is not constant along the sample width, as shown in the left-hand side of Figure 5. Consequently, the measure of the crack length close to the edge of the block may be highly different than its measure in the middle of the block, as described in the right hand-side of Figure 5. In order to overcome this difficulty, a new method based on measurement of the resulting force is proposed in order to estimate the crack growth rate. It is to note that the maximum force highly varies during the test, whereas the minimum force that corresponds to zero displacement remain almost constant. This could be simply explained. On one hand, the zero displacement configuration represents the closure of the crack; it can be concluded that the crack length has no effect on F mini . On the other hand, the curve of F maxi can be divided in three different parts, that correspond to three different physical phenomena:

• The beginning of the curve (indexed 1 in Figure 6) shows a high decrease of the resulting force during approximately 1000 cycles. This decrease is a consequence of the dissipative nature of elastomers. It includes the accommodation to the Mullins effect, the viscoelasticity of the material and the stabilisation of the block temperature. • The main part of the curve (indexed 2) corresponds to the crack propagation. The chosen geometry and the test conditions imposes that the tearing energy remain constant. So, the crack growth rate is also constant and the measured curve can be fitted by a linear equation. Its slope is equal to

N F ∆ ∆ max
and its intersection with the force-axis is F c0 (see Figure 6).

• The last part of the curve (indexed 3) corresponds to the catastrophic crack propagation that takes place before the breakage of the sample.

Considering the previous remarks, the determination of the crack growth rate is quite easy. Indeed, the crack growth rate only depends on the slope of the experimental curve and on some parameters previously obtained by FE calculations:

n F AF l n F F c dn dc c ∆ ∆ = ∆ ∆ ∆ ∆ = max 0 max max

EXPERIMENTAL RESULTS

CRACK CLOSURE EFFECT

Experiments under constant enforced maximal displacement were conducted to determine the influence of the minimal strain on the crack growth rate. Two materials were tested: a natural rubber (NR) which exhibits strain crystallisation effects and a polyurethane (PU) which does not. Experimental results are presented in Figure 7 for NR and in Figure 8 for PU.

1.E-08

1.E-07

1.E-06 -1.0 -0.5 0.0 0.5 1.0

Rl dc/dN

Maximal strain = 50%

1.E-09

1.E-08

1.E-07

1.E-06

1.E-05 -0.1 0.0 0.1 0.2 0.3 0.4 0.5

Rl dc/dN

Maximal strain = 125%

Maximal strain = 187.5%

Maximal strain = 250%

Figure 7 : influence of the minimal strain on the crack growth rate for PU Figure 8 : influence of the minimal strain on the crack growth rate for NR

These figures present the influence of the displacement factor, said the ratio of the minimum on the maximum displacements of the loading cycle, on the crack growth rate. Moreover, as each curve stands for a constant value of the maximum displacement (or strain), the figures highlight the influence of the minimum enforced displacement on the crack growth rate. They demonstrate that the magnitude of the minimal enforced displacement has no influence on the crack growth rate as long as it remains lower than a threshold value. This threshold value corresponds to the special case of a loading cycle without compressive part. Above this threshold value, the crack growth rate significantly decreases as R l increases. In regards with metallic material state of the art, this effect can be named the "crack closure effect". In order to account of this effect, experimental results are not analyse considering the maximum tearing energy T max but the maximum amplitude of the tearing energy during the cycle denoted ∆T and defined by: ∆T=T max -max(T min ,0)

Finally, it is to note that the "crack closure effect" is not sensitive to the strain crystallisation.

EXPERIMENTS UNDER NON-RELAXING CONDITIONS

Assuming now that the crack closure effect should be considered, we examine more precisely the influence of the pre-load on the crack growth rate law. As mentioned above, this influence is analysed by performing experiments with different values of the displacement ratio R l =L min /L max . Experimental results are presented in Figure 9 for PU and in Figure 10 for NR and lead to the following conclusions:

• If the test sample is compressed during the cycle, said 0 ≤ L R , then the crack growth rate law is not modified by the magnitude of compression (see Figure 9). Consequently, there is a unique crack growth rate law for the cases 0 ≤ L R .

• If the displacement ratio is higher than a positive threshold value, the effective stress cycle is tensile / tensile. In that case the crack growth rate law is shifted but the slope remains constant for all values of R l (see Figure 9 and Figure 10). It can be said that the material is reinforced whatever it crystallises or not.

-10,0 In conclusion of this study, we propose a first attempt to model the influence of non-relaxing conditions on the crack growth rate law. In order to develop our theory, the curves presented in Figure 8 are redrawn in a different way: for each value of the enforced displacement ratio, the corresponding crack growth rate is non-dimensionalised by the crack growth rate corresponding to R l =0. Here, this method is only applied to NR and the corresponding master curve is presented in Figure 11.

) ( )    > + × + × = < = ⋅ ⋅ = d Lthreshool L L L L d Lthreshool L L L R R c R R a R f R R R f with T B R f dn dc for b for 1 ) ( 2 0 0 α
The corresponding fitted curve is shown in Figure 11.

INFLUENCE OF THE TEMPERATURE

Automotive AVS components are classically used at different temperatures in service conditions. Indeed, it is important to be able to quantify the temperature effect on the crack growth rate law. That is the reason why cyclic experiments were performed on both materials at different temperatures. The crack growth rate laws are presented in Figure 12 for PU and in Figure 13 for NR. Examining the results and in accordance with previously published results, a major difference between the two materials can be observed:

• PU is very sensitive to the temperature: the crack growth rate law is shifted to higher values of dc/dn as the temperature increases (the material becomes stiffer as the temperature decreases). • NR results are more scattered but the material seems to be relatively insensitive to the temperature in comparison with PU, whatever the value of the enforced displacement. This observation is confirmed by Figure 14 in which the crack growth rate versus temperature curves are horizontal lines for each values of the enforced strain. The authors suppose that this sensitivity to the temperature can be explained by the strain crystallisation properties of the material considered. Some works are in progress to prove this proposition.

CONCLUSION

The present work is the preliminary step of a wider experimental study on fatigue life of elastomers. The aim of this study consists partly in adapting the classical approach developed previously to investigate the crack propagation under simple loading conditions to the case of anti-vibration system parts.

In this paper, the limitations of the tearing energy as defined by Rivlin and Thomas were highlighted and a new experimental simple shear sample was presented. This sample reproduces plane strain conditions in the neighbourhood of the crack front and can be submitted to complex loading histories. Moreover, the methodology used to determine the crack growth rate was detailed. First experimental results on both PU and NR were thoroughly analysed and the corresponding crack growth rate laws were established. Moreover, the influence of the pre-loading conditions and of the temperature on this law was also considered. All of these results are in good agreement with previously published results and validate our new sample.

Further works are in progress and new experiments are planned in order to investigate the influence of the loading frequency and of ageing conditions on the fatigue life of elastomers. The next step of the study will be concerned with the effect of complex loading conditions (enforced force or enforced displacement, form of the cycle, sample accommodation…) on the crack growth rate law.

  Figure 1: thin sheet test samples

Figure 2 :

 2 Figure 2 : Design of our new sample

  In regards with this equation, the maximum tearing energy only depends on experimental parameters ( the global stiffness of the sample that is a function of the crack length (

Figure 3 :

 3 Figure 3 : resulting force versus crack length.Figure4: linear fit.

Figure 5 :Figure 6 :

 56 Figure 5 : shape of the crack front for a sample under plane strain conditions Figure 6 presents the evolution of the maximum and minimum forces obtained during a fatigue test under controlled displacement conditions.

Figure 9 :

 9 Figure 9 : effect of non-relaxing conditions on the crack growth rate law for PU Figure 10 : effect of non-relaxing conditions on the crack growth rate law for NR

Figure 11 :

 11 Figure 11 : master curve for the crack growth rate law under non-relaxing conditions. Case of NRFinally, the crack growth rate law is approached by the following function:

(

Figure 12 :

 12 Figure 12 : crack growth rate law at different temperatures for PU

Figure 14 :

 14 Figure 14 : influence of the temperature on crack growth rate NR

ACKNOWLEDGEMENTS

The authors acknowledge the CETIM Nantes for its participation in the study and more especially for the fractographic analysis of the broken samples.