Modeling of composite multilayer rubber-steel: vibro-acoustic insulation of vehicle brake system

– The study of nonlinear dynamic behavior of laminate composed of steel and rubber layers also referred as “Shim”, used for vibro-acoustic insulation in brake system, is investigated. The simulation of the vibro-acoustic nonlinear behavior of Shim depending on frequency, taking into account the large deformations and various nonlinear hyper-viscoelastic laws of rubber are considered. This paper presents a solution to contribute in the identiﬁcation of the best design of Shim in terms of damping vibration of brake systems, using analytical and numerical method. The choice of the best structure depends essentially on the nature of rubber, on the stacking sequence of materials, on their thickness, on the number of layers and on volume fraction of rubber. An analytical study, with the use of the transfer matrix method is presented. A model on the ﬁnite element software ANSYS is constructed. The results lead to conclusions about the best structure and design of Shim in term of vibro-acoustic insulation.


Introduction
The properties of rubber (large deformations, hyperelasticity and damping) make their use very interesting.These materials are increasingly used in industries such as automotive and aerospace industries.They are introduced generally in pieces submitted to high mechanical loads (static and dynamic).The behavior model is strongly nonlinear.Rubber constitutes a solution of several major problems in the industrial field.Their specific mechanical properties in terms of hyper elasticity and visco-elasticity assuring the capacity of damping, justify their use in several industrial applications.For instance, rubber plays a fundamental role in the development of the vibrations in order to reduce noise of the braking action.These viscohyperelastic multilayered materials works mainly in compression and in shear modes and automotive.Among the several applications, rubber is used with steel to form a multilayered structures, called Shim, placed in the brake system.In this structure, thin rubber layers are introduced between steel laminates.The role of the Shim is to absorb vibrations in order to reduce noise of the breaking a Corresponding author: moez.chafra@ept.rnu.tnaction.The viscoelastic multilayered work mainly in compression and shear mode and their effectiveness depends essentially on the nature of rubber, disposition of materials, thickness and number of layers.Today, automotive industry use empirical methods to define the best structure to employ.Each test of these methods is repeated for each brake system and for each product (car, train, truck . . .).These trials have an important cost and the structure found is not necessarily the best structure.This work presents a solution to contribute in the identification of the best design of Shim in term of vibro-acoustic insulation with less cost and better efficiency.
Several studies have been employed to model the behavior of viscoelastic multilayer structures.For example, many previous works have focused on the design formulations for thin multilayer sandwich plates [2,3].Furthermore, several theories have been formulated and numerous approximations have been developed to simplify the modeling of this multilayered structure [6,7].Recent publications describe the analysis of the effect of nonlinear material on the global multilayer dynamic behavior [9].In this paper, the objective is to simulate analytically and numerically the non linear vibro-acoustic behavior of these Shims depending on the vibration frequency.The Article published by EDP Sciences model is based on a linearization procedure involving a finite static deformation which is a compression preload submitted to shear vibrations under this preload [15,16].The model is a mixed model considering linear viscoelastic vibrations over nonlinear large static deformations.Finally we have to analyse a linear viscoelastic problem in which creep and relaxation functions depend on the large static deformations and of visco-hyperelasticity laws.
The complex nature of the rubbers such as large strains and the various nonlinear hyper-viscoelastic laws can be taken into account.In this present study, an analytical approach using Mooney-Rivlin law is detailed and solved by means of MAPLE.Analytical results are compared with numerical results coming from homogenization methods using finite elements software (ANSYS).Finally, this model will permit to take into account a parametric study using specified criteria of potential industrial customer to get the best structure of the Shim in terms of damping and vibro-acoustic insulation.
In our approach we first use classical linear elasticity to derive the model frequency dependence and later the visco-hyperelasticity to introduce the hyper-viscoelastic effect of elastomer layers in terms of preload.It is the purpose of this study to analyse a non linear material's effect on the global multilayer behaviour in various geometric configurations.

Linear field
We are considering a multilayer cuboid composed [10] of N layers of total thickness h.In the present analysis, the origin of the orthonormal coordinate system (0, x, y, z) is chosen at the last-surface of the laminated plate [11], see Figure 1.The displacements of the plate in the x, y and z directions are denoted by u, v and ω, respectively.
We suppose that the top face of the shim is submitted to a vertical compression preload F z = F 0 in the z direction and to a transverse harmonic shear F y in the second direction y (Fig. 1).The bottom face of shim is clamped.
At each point of the shim, equations of motion in terms of Cauchy stress σ = (σ ij ) and displacement U = (u i ), are fulfilled: where b represents the body force vector and ρ the mass density.
We define a field displacement as follows: where U is a shear displacement in plane (x, z) defined by 2 unknown functions f (z) and g(z) (g(z) characterizes the preload).The determination of the constitutive equation in linear viscoelasticity can be further simplified if we used a Carson-Laplace transform [12].Carson-Laplace transform is defined as follows: and f (p) is defined as The constitutive equation for an isotropic non ageing material is as follows [13]: where λ and μ are Lamé parameters and δ = (δ ij ) is the Kronecker tensor.X denotes the tensorial nature of X = δ, ε or σ.

Complex modulus of nonlinear hyper-viscoelastic material
The rheological model of the rubber treated follows Zener rheological model with 3 parameters (k, τ and E ∞ ) and modeled by Mooney-Rivlin hyperelastic law defined by this expression of energy density w defined as follows [14]: In which I 1 and I 2 are the first and the second invariants of the left-hand Cauchy-Green dilatation tensor and A and B are Mooney-Rivlin parameters.In the present analysis, we treat an incompressible material where the expressions of I 1 and I 2 are defined as follows: where λ 0 designate the main compression elongation.
The tangent Young's modulus at elongation λ 0 , in static prompting, is given by the following relationship, Nashif et al. [16]: with We suggest, in dynamic prompting, a mixed model considering linear vibrations over a nonlinear static compression deformation.Beda and Chevalier [15] show that the total Piola-Kirchhoff stress is the sum of static stress (independent on time for long time) and dynamic stress (small vibrations).Considering viscoelastic point of view, the dynamic tensional stress at frequency ω is proportional to the complex tangent Young's modulus E * (λ 0 , ω) at static elongation λ 0 which is defined through relationship ( 8) by the following equation: where In which, A, B are the Mooney-Rivlin parameters and h * A (ω), h * B (ω), two complex functions representing the material characteristics.At low frequencies (ω ∼ 0) and at very high frequencies (ω ∼ infinite) the previous assumptions (9), (10) must lead to static tangent Young's moduli (8), that means the following conditions must be verified: Conditions (11) are satisfied in this particular case: Providing that Ferry [17] specifies that h(ω) can be a real function and then ( 9) becomes and through relationship ( 12): Functions h(ω) and ξ(ω) are determined by the knowledge of viscoelastic behavior of the material E 0 (ω) and ξ 0 (ω) without static load (λ 0 = 1) and then Equations ( 14) lead to Thus the real and the imaginary parts of the complex tangent Young's modulus is complex are defined as follow: Furthermore, h(ω) and ξ(ω) are two functions deduced from the viscoelastic characteristic without preload.In the case of shear vibration over a nonlinear static elongation we assume that the Poisson ratio ν of rubber is constant, generally about 0.47, so the shear tangent modulus G of the material is proportional to the tangent Young's modulus E, even in viscoelastic behavior (complex moduli).That means and thus the previous modeling can be used.Note that Gacem et al. [9] consider coupling between elongation preload and shear preload under Gent-Thomas hyperelastic law.

Young modulus without preload
The rheological model of the rubber treated follows the rheological model of Zener with 3 parameters.This model is presented in Figure 2.
The Young modulus, without preload (λ 0 = 1), is defined by the following equations:

Transfer matrix
Taking into account Equation ( 2) and equation of motion (1), we obtain, in the case of elastic behavior, the following solutions for the unknown functions f (z) and g(z) in each nth layer n = 1, 2, . . ., N: f n (z) = p n e iqnz + r n e −iqnz (18) In which q n = ρnω 2 μn is a wave number, ρ n and μ n = G n are respectively the mass density and the shear modulus for the nth layer and a n , d n , p n and r n are constants which are determined using the continuity conditions at the interface of each layer (continuity of stress and displacement).From Equations ( 2) and ( 8), we can deduce that the following expression of the shear displacement field in the y direction is as follows: In the nth layer the strain tensor is in the following form: Notice that shear strain ε 23 and compression strain ε 33 are not coupled (Eqs.( 17) and ( 18)).So we deduce the following stress tensor: In this work, we have used the transfer matrix method to link the mechanical properties of different layers of Shim [9].This method is needed to determine the global characterization of our structure, and in order to simplify equations, it requires a change of variable as stated below: hj e iqnς , affects only unknown constants p n and r n .So we obtain the displacement and the stress, Equations ( 19) and ( 23), fields for each layer in terms of ς: To proceed further, it is expedient to define matrix transfer equations relating magnitude of stress y , (25) (26) as follows in layer n [18]: −μ n q n sin(q n ς) iμ n q n cos(q n ς) Equation ( 27) can be rewritten as follows: where and Taking into account continuity of fields between the (n−1)th layer and the nth layer, we deduce that (Eq.( 28)): then relationship (28), at ς = h n , becomes: (30) So we write Equation (30) as following: In which With recursion, we obtain a relation between the displacement and the stress fields in the top layer (number N ) and in the bottom layer (number 1): ) We get this matrix form: In the case of homogeneous medium, in static, the shear stress τ is constant, thus the shear modulus G of material is directly proportional to the thickness h and inversely proportional to the shear displacement U 0 .Using Equation (28), the global shear modulus of the shim can be obtained as follows: where h = N j=1 h j is the thickness of the global structure.
Notice that Equation ( 35) is the result of static case and can be used, in dynamic case, only at low frequency range or for large wavelength in comparison with the Shim's thickness, see Gacem et al. [9].
The same procedure can be used for linear viscoelastic materials.We have just to replace elastic moduli by complex moduli at circular frequency ω.In the case hyperviscoelastic mixed model considering linear shear vibrations of the Shim over nonlinear static elongation, the complex shear moduli of each layer depends on the circular frequency ω and on elongation λ 0 .The wave-number q n , in Equation ( 18) becomes: Then matrices [α (n) ] and [α] become complex: and the effective complex shear modulus (35) becomes: The results of this analytic work were a subject of a parametric study.To avoid repeating the calculations for each modified parameters (stacking sequence of materials, their thickness and number of layers), we have used the analytical part in MAPLE software.

Parametric study
The main purpose of this study is to implement a numerical code to identify the best design of Shim in its role as vibro-acoustic insulation.To achieve this vision, we must go through a parametric study.We can change 3 parameters in this study: -Thickness of layers steel/rubber connected to the volume fraction v r of rubber.-Number of layers.
-Main compression elongation λ.This parameter is introduced due to the introduction of rubber nonlinear hyper-viscoelastic behavior.

Interpretation of the results
The figures below represent the damping factor as a function of the frequency for various cases of structure.The damping factor is a dimensionless quantity that expresses the ease with which a material dissipates energy in its vibration.The characteristics of rubber are presented in the previous paragraphs.For steel, we used classical characteristics Young's modulus equal to E = 220 GPa, a mass density equal to ρ = 7800 kg.m −3 and a Poisson's ration ν equal to 0.28.We should note that the thicknesses of layers presented in figures below are in cm.
The characteristics of the rubber are given in Section 2.3 for Young's modulus, and assuming that the Poisson ratio is constant (n = 0.47), the shear modulus without preload (λ 0 = 1) is given by Zener's modeling in which: -G ∞ = lim t→∞ G(t) is the asymptotic shear modulus at t = ∞; k is a dimensionless positive parameter; τ is a relaxation time; ω = 2π f is the applied circular frequency, t denotes the time.
Numerical values of parameters of the treated model are as follows: G ∞ = 1.15 GPa, k = 0.185, τ = 32.2 × 10 −6 s.We represent in Figure 3 the damping factor of the rubber material used.
When analyzing the Zener's viscoelastic modeling it can be shown, Chevalier and Vinh [19] Section 2.8.3, that the damping factor is maximum at frequency For the rubber material used f m = 4540 Hz as shown in Figure 3.
In Figure 4, we plot the variation of G's damping factor (tan δ) as function of the thickness of the layers of a 3-layer Shim.As well evident from Figure 4, for a fixed Shim thickness h, damping increases with the thickness of the rubber layer.It is interesting to analyze the variation of the effective shear damping factor versus the volume fraction of rubber v r = total thickness of rubber/h (in %). Figure 4 shows that the damping factor increases with v r : large variations when 10% < v r < 30%, small variations when 30% < v r < 80%.We can get a good damping factor with thin rubber layers.Figure 5 shows the variation of damping factor (tan δ) as a function of the number of layers of Shim structure.We note that the damping factor (tan δ) increases with the volume fraction of rubber v r .In this figure tg δ is quite when 30% < v r < 40%, so curves show that the variation of effective damping factor of the shim, versus number of rubber layer and versus the thickness h, is very small: 1.3%.
Figure 6 shows the variation of damping factor (tan δ) as function of the main elongation λ for the case of one rubber layer.
We note that as the elongation decreases as the damping factor increases. Figure 7 shows the variation of the shear storage (elastic portion) and loss moduli (viscous portion) of the shear modulus.In this figure, G decreases and G increases when the elongation increases.The variation of G and G , for different numbers of layers, depends on the thickness of rubber layer and specially on the volume fraction of rubber v r because the thickness of the shim is constant (10 −2 m).The storage shear modulus and the maximum of loss shear modulus decrease quite linearly with the volume fraction, when 20% < v r < 40%.

Finite-elements calculations
In the analytical approach, the solution has been presented for cuboid geometries.However, this approach cannot be developed when complex geometries, arbitrary boundary conditions or other nonlinearities are involved.Therefore, a numerical method is required by Miller [20] and Bhashyam [21] to solve such problems.
Numerical simulations were carried out using the finite element program ANSYS.In order to determine the effective shear modulus of the shim, we should apply a numerical shear test.The first step in the numerical method  is to introduce structure, mesh (Fig. 8) and constitutive laws of materials.We can then compute all components of the stiffness matrix.
To introduce the real behavior of rubber, we have used the shear modulus represented in terms of Prony's series, as assumed by Tzikang [22]: Here, G 0 and K 0 are, respectively, the shear and bulk moduli at the fast load limit (the instantaneous modulus), and G ∞ and K ∞ are the moduli at the slow limit.On ANSYS numerous hyperelastic models are presented, so we can define the rubber as a material modeled by a Mooney-Rivlin hyperelastic law.

Homogenization method
When choosing the displacement field given by Equations (2) we assume that the Shim is infinite in x and y directions.The analytical problem is then a "one dimensional problem" in z direction.The numerical shear test, using ANSYS software, needs 3-dimensional boundary conditions on the R.V.E (Representative Elementary Volume representing the Shim) as close as possible of the analytical shear test in plane (2,3), equivalent as shear test in plane (1,3).For a numerical shear test we suggest the following boundary conditions.The plane y 3 = 1 is subjected to a displacement equal to 1 in the direction of y 1 .The rotation in the plane y 3 = 1 is not allowed and the plane y 3 = 0 is clamped.The displacement of two planes y 2 = 0 and y 2 = 1 is blocked in the direction of y 2 .The rotation in the plane y 1 = 1 in the direction of y 2 is not allowed (Fig. 9).The homogenization procedure consists in expressing the mean value of stress on the R.V.E as a linear function of strains on the R.V.E.In the case of shear test, expressed by the previous boundary conditions, the procedure leads to the effective shear modulus of the shim.
For a homogeneous medium (Fig. 9) shear displacement in x direction is proportional to z, then the shear stress σ 13 is constant, thus the effective shear modulus G is directly proportional to the thickness h and inversely proportional the shear displacement U 0 in x direction at z = h.Finally, because of the linearity of the problem, there is no restriction to choose the shear displacement U 0 equal to 1 at z = h = 1 and no unities are necessary.Notice that the previous boundary conditions are not true when the shim is in operation (planes y 1 = 0 and y 2 = 1 are free).
To determine the function of the shear modulus G(t), we have used the curve fitting method.This method is a technique for identification curve by the identification of the function by the constants of the Prony's series.After identification of the function G(t), we apply the Carson-Laplace transformation (p = iω), or Fourier tansform, to identify the complex shear function G * (f ) for each frequency f = ω/2π.V r= 20% (5 layers)   Fig. 11.Variations of shear modulus damping factor of shims versus frequency: comparison between numerical and analytical cases.

Interpretation of the results
In Figure 10 we plot both our numerical results as well as the results from the analytical approach for one rubber layer.
This figure shows excellent agreement between the numerical and analytical results.The error between these two methods is about 1.7%.Figure 11 shows the variation of damping factor (tan δ) as a function of frequency for a fixed total thickness h = 10 −2 m and by varying the number of layers and the volume fraction of rubber v r .It can be notice that the accuracy between numerical and analytical results increases lightly with the number of layers.The rising of layers can explain this phenomenon it is the fact the stiffness of the schim increases with steel layers and considering the boundary conditions on the schim's thickness, the stress and strain fields in numerical approach are close to stress and strain fields in analytical approach.
We note a good correlation between numerical and analytical models.Figure 12   numerical results of the evolution of G for 2 values of the elongation.We remark an excellent agreement between the both methods.
In Figure 13, we plot both our analytical results as well as the results from the numerical analysis.This figure shows the variation of damping factor (tan δ) as a function of the frequency for different structures of Shim.We considered 2 Shims including one composed of a single rubber layer and the second is composed of 3 layers (2 steel and 1 rubber).We remark a good correlation between the numerical and the analytical models.
Unlike analytical analysis (infinite medium in x and y direction) numerical finite elements methods need three dimensional boundary conditions.So we conclude the validation of the numeric method (choice of boundary conditions) in comparison of the analytical method in the evaluation of effective shear modulus of the schim.Notice that the numerical approach is not true when the schim is in operation.These 2 methods may be useful for the automobile industry.In this project, industrial data about the nature of the rubber are welcome to use, we can provide the best design in terms of vibro-acoustic damping.

Conclusion
In this paper, analytical and numerical methods have been developed to contribute in the determination of the best design of Shim in terms of damping vibrations of the brake systems.The choice of the best structure depends essentially on the nature of rubber, and on the volume fraction of rubber, stacking sequence of materials, their thickness, and the number of layers.The vibro-acoustic nonlinear behavior of Shim, the large deformations and various non linear hyperviscoelastic laws of rubber were also considered.In the present paper, we have attempted to expose and explain the different steps followed to achieve this work.
Our results indicate a good accuracy between analytical and numerical models.Thus, the analytical method has the advantage to obtain accurate results in simple computation process while the numerical approach is closer to reality and can take into account more complex systems and situations.

Fig. 1 .
Fig. 1.Example of Shim structure submitted to compression and shear load.
24)where 0 ≤ ς ≤ h n and h n is the thickness of layer n.The property of exponential function, e iqnz =

Fig. 3 .
Fig. 3. Variation of shear modulus damping factor of rubber material as a function of frequency (one rubber layer).

Fig. 4 .
Fig. 4. Variation of shear modulus damping factor of a three-layered shim (S/R/S) as a function of frequency with changing the thickness of layers (or the volume fraction of rubber vr in %).The thickness h of the shim is constant (h = 1 cm).

Fig. 10 .
Fig. 10.Analytical and numerical variations of shear modulus damping factor of rubber versus frequency (the case of a single rubber layer).

Fig. 12 .
Fig. 12.Comparison between numerical and analytical results for shear modulus G (loss modulus) versus frequency for various elongations λ.

Fig. 13 .
Fig. 13.Variation of shear modulus damping factor of rubber and of 3-layered shim versus frequency for different values of main elongation λ: comparison between numerical and analytical cases.

-
The shear loss shear modulus of the shim G decreases and the storage modulus G increases when the elongation increases.-The damping factor of the shim increases with the volume fraction of rubber v r : large variations for small values of v r , small variations for large values of v r .That we can get a good damping factor with thin layers of rubber.-The effects of the number of rubber layers and of the global thickness of the shim are very small.