Marc-Antoine A Campana

Morvan Ouisse

Emeline Sadoulet

Massimo Ruzzene

Neil Smith

Fabrizio Scarpa

E Sadoulet-Reboul

S Neild

Impact of non-linear resonators in periodic structures using a perturbation approach

Keywords: Periodic structures, Resonators, Non-linearities

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

Introduction

The Floquet-Bloch theorem is widely used to investigate the wave propagation in periodic structures, initially introduced in mathematics by Floquet for [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF], before being extended in structures [START_REF] Rayleigh | on the maintenance of vibrations by forces of double frequency, and on the propagation of waves through a medium endowed with a periodic structure[END_REF][START_REF] Brillouin | Wave Propagation in Periodic Structures[END_REF] and damped systems [START_REF] Collet | Floquet-bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems[END_REF] by analysing the periodic unit cell only. One of the main interests for analysing the wave propagation in these periodic media is to identify "stop bands" or "bandgaps", i.e. frequency areas in which waves are no longer free to propagate. The knowledge of these areas in the wave number/frequency maps can be used to build frequency filters using periodic systems and to create therefore frequency bands in which wave propagation is strongly attenuated. Two types of stop bands can be identified in periodic systems. Bragg bandgaps, related to the material properties or the geometry of the unit cell, and resonant bandgaps, caused by a resonating system included in the periodic structure. In a one dimensional spring-mass resonant system, the combination of a resonating mass with the mass of the principal chain results to an apparent frequency-dependent mass, which has the specificity to become negative for a certain range of frequencies close to the eigenfrequency of the resonator. This phenomena is called negative mass effect [START_REF] Yao | Experimental study on negative effective mass in a 1d mass-spring system[END_REF], and it is related to the presence of a bandgap in those systems.

Several studies have been performed by using resonators systems in periodic structures. Examples are multi-resonators in a single unit cell [START_REF] Zhou | Effects of relevant parameters on the bandgaps of acoustic metamaterials with multi-resonators[END_REF], or the effect of a graded metamaterial [START_REF] Banerjee | Frequency graded 1d metamaterials : A study on the attenuation bands[END_REF] to generate pseudo-periodic systems. A large portion of studies related to wave propagation in periodic arrays (also including damping effects) have been described and reviewed in [START_REF] Banerjee | Waves in structured mediums or metamaterials : A review[END_REF]. Some studies have also focused on beam sandwich structures composed by homogenised continuum media with local resonators [START_REF] Chen | Dynamic behaviour of sandwich structure containing spring-mass resonators[END_REF]. Qian [START_REF] Qian | Bandgap properties in locally resonant phononic crystal double panel structures with periodically attached spring-mass resonators[END_REF] has investigated two-dimensional periodic sandwich plate continuum systems using resonators between the two skin panels.

Studies about periodic non-linear systems have been performed to observe the effect of the amplitude of the response during wave propagation, and more specifically on the boundaries of the bandgap compared to a linear periodic system.

Non-linear resonators in a linear periodic chain have been investigated in previous studies. Lazarov and Jensen [START_REF] Lazarov | Low-frequency band gaps in chains with attached non-linear oscillators[END_REF] have used an harmonic balance method to analyse the asymptotic behaviour of the relative displacement between the principal mass and the resonant mass. Georgiou and Valakis [START_REF] Georgiou | An invariant manifold approach for studying waves in a one-dimensional array of non-linear oscillators[END_REF] have used a different type of oscillator using an geometric non-linearity originating from the angle of a pendulum attached to the principal mass of the chain.

The work presented in this paper focuses on developing a formulation to identify the dispersion curve of a non-linear resonator using the perturbation method developed by Nasiretti et al [START_REF] Narisetti | A perturbation approach for predicting wave propagation in one-dimensional nonlinear periodic structures[END_REF] for discrete non-linear spring-mass systems, and then extended to finite element structures by Manktelow et al [START_REF] Manktelow | Finite-element based perturbation analysis of wave propagation in nonlinear periodic structures[END_REF]. The method has been so far applied to a full nonlinear periodic system, and here the behaviour of the nonlinear resonator only is analysed to observe the impact of the eigenfrequency of the resonant system in the corrected term of the dispersion relation. Another novelty presented in this work is the development of a numerical inverse method used on a finite periodic chain structure to obtain the dispersion relation. The method consists in imposing the wave number to the whole chain and then obtain the corresponding frequency analyzing the temporal signal obtained after releasing the system. This method is different from usually applied approaches that consist on imposing an excitation and identifying back the wave number [START_REF] Manktelow | Analysis and experimental estimation of nonlinear dispersion in a periodic string[END_REF], and has the particularity to avoid space aliasing as well as control the amplitude we want to impose to be compared with the perturbation technique. Also, the precise observation of the boundaries of the bandgap can be performed since the wavenumber is imposed.

Dispersion analysis for infinite periodic structures

This section details strategies to determine dispersion relations for infinite periodic structures : the case of a linear spring resonator is first considered, then more attention is dedicated to the case of a nonlinear spring-mass resonator.

Linear spring-mass resonator

Let consider a one-dimensional principal periodic chain made of spring-mass units (mass m and stiffness k). A resonator represented by another springmass system (mass m R and stiffness k R ) is attached to each mass. This chain is considered as infinite, and an appropriate unit cell is defined to apply the Floquet-Bloch Theorem (Fig. 1).

   m d 2 uj dt 2 + k(2u j -u j-1 -u j+1 ) -k R (q j -u j ) = 0 m R d 2 qj dt 2 + k R (q j -u j ) = 0 (1)
where u j and q j respectively denote the displacements of the j th principal mass and the j th resonator. In harmonic regime the displacements can be expressed as x j = X j e iωt , with i being a complex number such as i 2 = -1. Eq. 1 becomes :

70    -ω 2 mu j + k(2u j -u j-1 -u j+1 ) -k R (q j -u j ) = 0 -ω 2 m R q j + k R (q j -u j ) = 0 (2)
The term q j depends on u j from the second line of Eq. 2. By applying the Floquet-Bloch Theorem to the first line (u j+1 = e jµ u j ) (µ being the reduced wave number) one can obtain the following dispersion equation :

-mω 2 + 2k(1 -cos(µ)) + (k R - k 2 R k R -m R ω 2 ) = 0 (3)
The roots of Eq. 3 are ω 1 and ω 2 , representing the two branches of the dispersion curve :

ω 1 = (m+m R )Ω 2 R +4k sin 2 (µ/2)- √ ((m+m R )Ω 2 R +4k sin 2 (µ/2)) 2 -16km sin 2 (µ/2)Ω 2 R 2m ω 2 = (m+m R )Ω 2 R +4k sin 2 (µ/2)+ √ ((m+m R )Ω 2 R +4k sin 2 (µ/2)) 2 -16km sin 2 (µ/2)Ω 2 R 2m (4)
with

Ω R = k R
m R being the eigenfrequency of the resonator. The dispersion curve is shown Fig. 2a, with the propagative part in red and the evanescent part in blue for m R = m = 1 kg and k R = k = 1 N.m -1 . To validate numerically the result, a finite spring-mass chain of 20 unit cells has been simulated to obtain the Frequency Response Function (FRF) of the system. The boundary conditions in the numerical model are representative of a free-free state at the two extremities of the chain, while the mass located in one of the ends is excited by imposing a longitudinal propagating wave generated by the force F = F 0 cos(ωt) with F 0 = 1N . An harmonic analysis is performed by varying the value of ω, and the amplitude of the displacement is captured on the mass at the other end of the structure. The results are shown in Fig. 2b.

The FRF exhibits the presence of a frequency region in which the amplitude is considerably attenuated. This area represents the bandgap (resonant in this case). A resonant bandgap can be differentiated from a Bragg one by observing the boundaries of the propagating part [START_REF] Sharma | Local resonance and bragg bandgaps in sandwich beams containing periodically inserted resonators[END_REF]. In the resonant case, a wave will stop propagating for a value of µ = π and then propagates again for µ = 0.

The boundaries of the bandgap can be expressed by calculating the values of ω 1 (µ = π) and ω 2 (µ = 0), leading to the following expression :

ω 2 g = (m+m R )Ω 2 R +4k- √ ((m+m R )Ω 2 R +4k) 2 -16kmΩ 2 R 2m ω 2 d = Ω 2 R (1 + m R m ) (5)
with ω g and ω d being the pulsation of the left and the right boundary, respectively. It is also important to notice that the value of Ω R is bounded by the values of ω g and ω d , and this has consequences on the value of the amplitude of the resonant mass. The rewriting of the second line of Eq. 2 leads to

q j = Ω 2 R Ω 2
R -ω 2 u j . This expression indicates that the value of the amplitude of q j The term q j will however never be able to reach an infinite value in the propagative zone, since ω -< ω g < Ω R < ω d < ω + , with ω -and ω + representing respectively the value of the pulsation before ω g and after ω d . This property of a resonating system is different from a representation of spring-mass systems with Bragg bandgaps, like the diatomic mass system of Hussein et al. [START_REF] Hussein | Dynamics of phononic materials and structures : Historical origins, recent progress, and future outlook[END_REF]. This particular periodic configuration shows amplitudes of the two masses that depend directly to the boundaries of the bandgap. This may have repercussions for the representation on the nonlinear system.

Nonlinear spring-mass resonator system

Perturbation method for identification of non-linear dispersion diagram

A Duffing spring mass resonator is added between the masses m and m R . In this case the non-linear force is expressed as :

f nl = k R (q j -u j ) + εΓ(q j -u j ) 3 ( 6
)
where Γ is the nonlinear stiffness constant and ε a small perturbation parameter. Replacing f nl in Eq. 1 leads to :

   m d 2 uj dt 2 + k(2u j -u j-1 -u j+1 ) -k R (q j -u j ) -εΓ(q j -u j ) 3 = 0 m R d 2 qj dt 2 + k R (q j -u j ) + εΓ(q j -u j ) 3 = 0 (7)
The objective here is to idenify the correction term in the dispersion relation using the Lindstedt-Poincaré adapted method developed in [START_REF] Narisetti | A perturbation approach for predicting wave propagation in one-dimensional nonlinear periodic structures[END_REF]. The main goal is to generate a first order asymptotic development on the displacements of the masses and on the frequency, corresponding to the following expressions :

u j = u (0) j + εu (1) j + O(ε 2 ), q j = q (0) j + εq (1) j + O(ε 2 ), ω = ω 0 + εω 1 + O(ε 2 ). (8)
Replacing expressions 8 in Eq. 7 gives, after development and separation of the equations in the ε 0 and ε 1 orders :

ε 0 :    ω 2 0 d 2 u (0) j dτ 2 + (2u (0) j -u (0) j-1 -u (0) j+1 ) -β(q (0) j -u (0) j ) = 0, κ 2 ω 2 0 d 2 q (0) j dτ 2 + (q (0) j -u (0) j ) = 0, (9)
ε 1 :    ω 2 0 d 2 u (1) j dτ 2 + (2u (1) j -u (1) j-1 -u (1) j+1 ) -β(q (1) j -u (1) j ) = -2ω 0 ω 1 d 2 u (0) j dτ 2 + Γ(q (0) j -u (0) j ) 3 κ 2 ω 2 0 d 2 q (1) j dτ 2 + (q (1) j -u (1) j ) = -2κ 2 ω 0 ω 1 d 2 q (0) j dτ 2 -Γ β (q (0) j -u (0) j ) 3 (10)
with

ω 2 n = ω 2 n Ω 2 0 , Ω 2 0 = k m , α = m R m , β = k R k , κ = α β , Γ = Γ k , τ = ωt. 120
Assuming an harmonic regime is reached, the displacements u (0) j and q (0) j can be expressed in the following form :

u (0) j = Au 2 e jiµ e iτ + Āu 2 e -jiµ e -iτ q (0) j = Aq 2 e jiµ e iτ + Āq 2 e -jiµ e -iτ (11)
The second order derivatives are then :

d 2 u (0) j dτ 2 = -u (0) j d 2 q (0) j dτ 2 = -q (0) j (12)
Replacing Eq. 11 and 12 into Eq. 9 yields

ε 0 :    -ω 2 0 u (0) j + (2u (0) j -u (0) j-1 -u (0) j+1 ) -β(q (0) j -u (0) j ) = 0 -κ 2 ω 2 0 q (0) j + (q (0) j -u (0) j ) = 0 (13)
Eq. 13 is equivalent to Eq. 2. Following the same process as in the linear case, roots are now expressed as

(ω (1) 0 ) 2 = 1+α+4κ 2 sin 2 ( µ 2 )- √ (1+α+4κ 2 sin 2 ( µ 2 )) 2 -16κ 2 sin 2 ( µ 2 ) 2κ 2 , (ω (2) 0 ) 2 = 1+α+4κ 2 sin 2 ( µ 2 )+ √ (1+α+4κ 2 sin 2 ( µ 2 )) 2 -16κ 2 sin 2 ( µ 2 ) 2κ 2 . ( 14
)
When adopting the same approach used during the linear analysis, the purpose is to isolate the corrected displacement q

(1) j in the second line of Eq. 10 and inject it in the first line. In this way it is possible to obtain all the nonlinear terms in the same equation. To do so, it is assumed that the displacement q (1) j can be expressed under the form :

q (1) j = Be iτ + Be -iτ (15)
Where B the amplitude of q

(1) j , whose the value does not have importance for the following calculations. Replacing Eq. 15 in the second line of Eq. 10 we obtain :

ε 1 :    ω 2 0 d 2 u (1) j dτ 2 + (2u (1) j -u (1) j-1 -u (1) j+1 ) -β(q (1) j -u (1) j ) = -2ω 0 ω 1 d 2 u (0) j dτ 2 + Γ(q (0) j -u (0) j ) 3 q (1) j = 1 1-κ 2 ω 2 0 (u (1) j -2κ 2 ω 0 ω 1 d 2 q (0) j dτ 2 -Γ β (q (0) j -u (0) j ) 3 ) (16)
By injecting expression of q

(1) j in the first line of Eq. 10 we have :

ω 2 0 d 2 u (1) j dτ 2 + (2u (1)
j -u (1)
j-1 -u (1) j+1 ) -β κ 2 ω 2 0 1 -κ 2 ω 2 0 u (1)
j = F (τ ) (17)
with :

F (τ ) = 2ω 0 ω 1 α + (1 -κ 2 ω 2 0 ) 2 (1 -κ 2 ω 2 0 ) 2 u (0) j -Γ κ 2 ω 2 0 1 -κ 2 ω 2 0 4 (u (0) j ) 3 (18)
Replacing u (0) j by its expression 11, the equation becomes :

F (τ ) = ω 0 ω 1 α + (1 -κ 2 ω 2 0 ) 2 (1 -κ 2 ω 2 0 ) 2 A u - 3 Γ 8 κ 2 ω 2 0 1 -κ 2 ω 2 0 4 A 2 u Āu e iτ e ijµ +d 1 e 3iτ e 3ijµ (19)
The terms d 1 are associated to the 3 rd order of the nonlinearity. The linear kernel of Eq. 17 is similar to the one in Eq. 3 ; this implies that one needs to have all the coefficients in e ijµ must be equal to 0 not to obtain a secular term in the temporal expression of u

j . By imposing this condition and rearranging the equation, we obtain the following expression for ω 1 :

ω 1 = 3 Γ|A u | 2 8 κ 8 ω 7 0 (1 -κ 2 ω 2 0 ) 2 (α + (1 -κ 2 ω 2 0 ) 2 ) (20)
Equation 20 leads to the establishment of the final expression describing the corrected dispersion relation for the periodic structure :

ω (p) = ω (p) 0 + ε 3 Γ|A u | 2 8 κ 8 (ω (p) 0 ) 7 (1 -κ 2 (ω (p) 0 ) 2 ) 2 (α + (1 -κ 2 (ω (p) 0 ) 2 ) 2 ) (21)
with p = 1 representing the left branch, and p = 2 being the right branch of the dispersion curve. Eq. 21 illustrates the importance of the pulsation of the resonator to determine the value of ω 1 . If κ 2 ω 2 0 approaches 1 (i.e., equivalent to say that ω approaches Ω R ), the value of ω 1 will increase and might reach a point where it becomes greater than ω 0 . This would however contradict the hypothesis underlying the current perturbation method. That condition shows that the value of the parameters for the resonators plays a crucial role in the distance that Ω R must possess compared to the boundaries of the bandgap.

During the following analyses we will consider the parameters α,κ, Γ are considered to be equal to 1 (unless Γ = 0, and this to obtain the linear case).

The graphic representation of the corrected dispersion curve obtained is shown Fig. 4.

Domain of amplitude validity

The perturbation analysis assumes an asymptotic development of the term ω to obtain the nonlinear corrected result. However, the asymptotic development of the term ω 2 = (ω 0 + εω 1 ) 2 = ω 2 0 + 2εω 0 ω 1 leading to Eq. 10 is true only if 2εω 0 ω 1 ω 2 0 , which can be rewritten as 2ε leading to the following condition :

3|A u | 2 ε 4 ω 6 0 (1 -ω 2 0 ) 2 (1 + (1 -ω 2 0 ) 2 ) 1 (22)
By assuming that the order of magnitude of is the same as ε, the following condition for the amplitude can be written :

|A u |h(ω 0 ) < 1 (23)
where h(ω 0 ) = 3 This means that the hypothesis of linearisation will be true, even for very high values of amplitude. However, the amplitude results to be much lower when the wave number increases, approaching a maximum theoretical limit of 0.48. The opposite behaviour however happens in the right branch (Fig. 5b) ; in this case it is possible to observe a limit for A u , with a minimum value around 0.6 in the low wave numbers ranger and reaching 1.8 for high values of µ.

Identification of dispersion curves from analysis of finite nonlinear periodic structures

This section presents a numerical method to obtain the apparent wave number of the structure, followed then by results related to this approach and a comparison with the method described in 2.

Inverse method

In general the techniques used to identify the apparent wave number of a periodic structure are based on applying an harmonic excitation to one of the masses of the periodic structure and then evaluating the response. This paper introduces here other numerical method to compare the results obtained from the perturbation theory. The idea underpinning this methodology is to impose the wave number of the whole structure as an initial condition, releasing the system and observing how it evolves in the time, and then obtain the frequency of the corresponding imposed wavenumber as output. One of the advantages of this novel methodology is that a time-domain analysis is required, rather than the space-domain analysis one to obtain the wave number in direct methods [START_REF] Berthaut | K-space identification of apparent structural behaviour[END_REF]. The use of a time-domain analysis helps to reduce the error made due to the spatial discretisation of the periodic structure. In particular for spring-mass systems, the distance between two consecutive masses does not need to be represented in the model, and this allows to fix the element length to unity. According to the Shannon theorem, a minimum of 2 elements per wavelength is necessary to avoid space aliasing. However, to represent a sinusoidal excitation in a correct way, it is common to use at least 6 elements per wavelength, and that condition can not be reached when the wavenumber exceeds 1. Fig. 6 illustrates this problem by showing the space representation for two different wave numbers. It is however important to note that this method fits mainly numerical simulations, since imposing a displacement as an initial conditions could be much less practical in an experiment than imposing an harmonic excitation to the whole structure.

The finite structure here is represented by an assembly of 300 unit cells in a linear chain. The equation of motion of the finite structure with a free-free boundary condition is written and solved using ODE45 in Matlab. A perfectly matched layer (PML) is also used at the boundaries of the structure to avoid wave reflection. The PML consists in a viscous damping force c(x p ). ẋp applied to all the masses of the principal chain under the form :

c(x p ) = Ce 1- 1 (1-xp ) n (24)
with C a constant positive value, x p ∈ [0; 2] a value interpolating the mass 210 at the p th position such as c(0) = c(2) = C and c(1) = 0, and n an even number representing the order of the PML. The higher the value of n is, the less the masses close to the boundaries will be affected by the damping. Fig. 7 shows the trend of the PML for two different values of n and C = 2.

The resolution of the linear system (i.e. with Γ = 0) is has to be done in 215 order to obtain the values of ω 0 (µ) corresponding to the dispersion curve in the linear case. This step is necessary to impose the initial condition of the resonating system, written as

A q = A u 1 -κ 2 ω 2 0 . (25)
Those initial conditions are then written under the following form :

                     U j (t = 0) = A u cos(µ(x j -x 0 )) Uj (t = 0) = 0 Q j (t = 0) = A q cos(µ(x j -x 0 )) Qj (t = 0) = 0 j ∈ [[1; 300]] (26)
The terms U j and Q j are the values of the displacements of the j th principal and resonating mass of the system, respectively. The location of the mass considered to obtain the time-displacement response is x 0 , and µ is the reduced wave number. The value of x 0 has to be chosen carefully in order to obtain a coherent representation of the results : if the chosen mass is too close to the boundaries, the PML and eventually residual wave reflections will provide incoherent results.

For that reason the observed mass should be located in the middle of the structure. Also, even if the displacement is imposed, the spatial discretisation will still be inaccurate for wave numbers greater than one ; this implies however the risk of not always obtaining the maximum value of the amplitude if a random mass is selected as an output. The writing of the initial spatial displacement as in Eq. 26 ensures that the amplitude A u will be observable for the mass located in x 0 for any value of µ, when x j = x 0 .

After the masses being released for an imposed wave number µ 0 and after a time equivalent to t = 25T 0 such as T 0 = 2π ω0(µ0) , the frequency of the time signal coming from the mass located in x 0 is measured. One example of input and output signal obtained for a µ = 1 on the right branch of the dispersion curve is shown Fig. 8. Fig. 8b shows that the amplitude does not remain constant and equal to A u as long as the time goes by. Consequently, a Fast Fourier Transform would not give accurate results here, since the obtained signal have variations of amplitude and therefore variations of frequency due to the nonlinear behaviour. Alternatively, an average of the periods of the signal is done to get the global value of the frequency. Depending on the imposed wave number, the value of the output average amplitude will become different, which makes not A u by comparing it with the output amplitude of the system (noted A out ). After an iteration, the value of A out is estimated performing a Hilbert transform of the signal :

A out = |H(U 0 (t))| (27)
If A out < A u , the value of A u0 is slightly increased and the process is repeated until A out = A u with a certain tolerance. Once this condition is done, one can capture the average frequency of the signal to obtain the combination (µ, ω).

Fig. 8 shows the results using initial conditions shown in Fig. 9, with a value of A u0 = 0.53 to get an average value A out = 0.5. The dashed line represents the value of amplitude equal 0.5.

Numerical results

Before obtaining the numerical nonlinear dispersion curves, a verification of the condition estimated in subsection 2.2.2 has to be performed to verify that the domain of validity of the amplitudes matches with the prediction when using

f err = |f num -f per | f per (28)
with f num the frequency obtained with the numerical simulations and f per the frequency obtained with the perturbation method. On the other hand, the value of 2 ω 0 ω 1 representing the condition leading to Eq. is calculated 23 is calculated. The highest the value of this function, the higher the error of the perturbation hypothesis. Fig. 10a shows that for high values of µ and amplitudes, the error increases rapidly. This curve actually generalises figure 5a, in which the relative error was fixed to 1, to give an overview of how fast this error is evolving for higher values of amplitude. Fig. 10b shows that the error tends to follow patterns similar to the previous result, and it also shows the same behaviour for high values of amplitude. Similar comments can be provided for the results presented in Fig. 11. Also this case, the error rapidly increases after A u 0.3. Those estimations are close to the ones found with the asymptotic assumption shown in 2.2.2 ; they are however not exactly the same because of the numerical precision of the solver used for the finite structure, and the estimation of taking 280 ε as "small".

The corrected dispersion curve for a amplitude value of A u = 0.3 is represented Fig. 12, and zooms on the branches are shown Fig. 13 and14. From the close looks it is possible to observe the difference between the numerical method and the perturbation approach. One can observe that for this particular value 285 of amplitude, numerical and theoretical results provide a close match, and the dispersion curve is shifted to higher frequencies compared to the linear case. demonstrates that the results will be inaccurate up to a certain value of wavenumber (around 0.8 in that case) but will still remain correct after this value.

Conclusion

Nonlinear periodic resonating systems made of baseline linear chains and distributed Duffing oscillators have been evaluated in this work using a perturbation approach. A quantitative study of the maximum admissible amplitude has been performed to provide an overview of the limits of the proposed method.

A new inverse approach imposing the wavenumber has also been developed and applied to a finite structure with a large number of unit cells to confirm the results obtained with the perturbation method. This numerical inverse identification approach provides a different way of analysing the wave propagation of nonlinear periodic structures. The approaches proposed in this work could be also extended to multi dimensional structures, as well as finite element structures to observe the influence of the only resonant part of a system and perform the inverse method in more complex cases.

Acknowledgements

This work has been supported by the MSCA ITN VIPER (Vibroacoustics of Periodic Media) programme.

  Floquet-Bloch provides the dispersion relation of an assumed infinite periodic structure Preprint submitted to Mechanical Systems and Signal Processing 11 juin 2019

Figure 1 :

 1 Figure 1: 1D resonant spring mass system

Figure 2 :Figure 3 :

 23 Figure 2: Results for one dimensional linear resonant spring-mass system

Figure 4 :

 4 Figure 4: Dispersion curve using the perturbation approach in the nonlinear case for Au = 1 and ε = 0.05

2 0 ) 2 )

 22 . Plotting the function 1 h will give the 165 maximum admissible value that |A u | can assume to satisfy the equation 23.The results are shown figure5. Fig.5ashows that for values of µ up to 1.0 almost any value of A u would not affect the dispersion in the nonlinear regime.

Figure 5 :

 5 Figure 5: Maximum of amplitude versus wave number analysis

Figure 6 :

 6 Figure 6: Spatial representation of masses for two different values of µ.

Figure 7 :

 7 Figure 7: Representation of the PML

  Time signal obtained for the mass x0

Figure 8 :

 8 Figure 8: Exemple of input and output signal for Au = 0.5 and µ = 0.1

Figure 9 :

 9 Figure 9: Input and output signals obtained after correction

Figure 10 :

 10 Figure 10: Model error estimation for the left branch of the dispersion curve

Figure 11 :Figure 12 :

 1112 Figure 11: Model error estimation for the right branch of the dispersion curve

Fig.15b confirmFigure 13 :

 13 Fig.15bconfirm that the high level of amplitude does not satisfy the hypothesis of the Lindstet-Poincaré approach. 295

Figure 14 :

 14 Figure 14: Zoom on the right branch of the dispersion curve (Au = 0.3, ε = 0.05)

Figure 15 :

 15 Figure 15: Results for an amplitude overpassing the nonlinear assumptions limit

This project has received funding from the European Union's Horizon 2020 research and innovation programme under Marie Curie grant agreement No 675441. This work was performed in cooperation with the EUR EIPHI program (ANR 17-EURE-0002).


Anne-Clémence Corbin

Damien Soulat

Manuela Ferreira

Ahmad-Rashed Labanieh

Xavier Gabrion

Pierrick Malécot

Vincent Placet

Towards hemp fabrics for high-performance composites: Influence of weave pattern and features

Keywords: A. Polymer-matrix composites (PMCs), B. Mechanical properties, C. Numerical analysis, E. Weaving

Recent developments in the field of bio-based composite materials are mainly focused on the use of unidirectional reinforcements. The production of woven fabrics and required yarns or rovings is still complex for composite applications due to the finite length of plant fibers and to the high number of process parameters which can be tuned. This study focused on the influence of weave pattern and process parameters on the resulting material properties at different scales. Results from mechanical characterizations and X-ray nanotomography show that very competitive tensile properties can be obtained for woven hemp fabric composites made from low-twisted rovings, in particular when compared to the front-runner flax cross-ply laminate.

Introduction

The increase in environmental awareness has stimulated the idea of reducing synthetic-based products in the production of composite materials and more specifically of reinforcements [START_REF] Dhar Malingam | Mechanical Characterisation of Woven Kenaf Fabric as Reinforcement for Composite Materials[END_REF]. To achieve this objective of reducing the use of synthetic fibers, the incorporation of cellulosic fibers such as flax, hemp or kenaf has received the attention [START_REF] Gourier | Influence of PA11 and PP thermoplastic polymers on recycling stability of unidirectional flax fibre reinforced biocomposites[END_REF][START_REF] Bourmaud | Towards the design of high-performance plant fibre composites[END_REF][START_REF] Masto | Investigation of the possible origins of the differences in mechanical properties of hemp and flax fibres: A numerical study based on sensitivity analysis[END_REF]. Indeed, their specific mechanical properties are competitive when compared to glass fibers and they have greater environmental benefits [START_REF] Faruk | Biocomposites reinforced with natural fibers: 2000-2010[END_REF][START_REF] Gurunathan | A review of the recent developments in biocomposites based on natural fibres and their application perspectives[END_REF][START_REF] Martínez | LCA of biocomposites versus conventional products[END_REF][START_REF] Bachmann | Environmental analysis of innovative sustainable composites with potential use in aviation sector-A life cycle assessment review[END_REF]. Since the mechanical properties of reinforced composites depend mainly on the nature of the reinforcing fibers, their quantity (fiber volume fraction), their placement and distribution, textile preforms (woven, braided, knitted fabrics, etc.) are widely studied to optimize reinforcement architectures [START_REF] Katnam | Towards balancing in-plane mechanical properties and impact damage tolerance of composite laminates using quasi-UD woven fabrics with hybrid warp yarns[END_REF][START_REF] Mouritz | Review of applications for advanced threedimensional fibre textile composites[END_REF]. In the textile process, there is direct control over fiber placements and ease of handling of fibers. Textile preforming operations play a key role in most of the composite manufacturing processes. The basic textile yarn and fabric forming processes have been modified and developed to a significant extent to meet the increasing demand from the composite manufacturing sector.

At the scale of the composite material, the best uniaxial mechanical properties are generally obtained when the fibers are aligned in the direction of the applied load [START_REF] Swolfs | Recent advances in fibre-hybrid composites: materials selection, opportunities and applications[END_REF]. Unfortunately, due to their finite length, it is more difficult to obtain an alignment with plant fibers than with continuous synthetic fibers [START_REF] Pickering | A review of recent developments in natural fibre composites and their mechanical performance[END_REF]. Consequently, many studies on plant fiber composites are related to the development and use of different types of reinforcements (unidirectional (UD), non-woven and woven preforms) [START_REF] Gourier | Influence of PA11 and PP thermoplastic polymers on recycling stability of unidirectional flax fibre reinforced biocomposites[END_REF][START_REF] Cadu | What are the key parameters to produce a high-grade bio-based composite? Application to flax/epoxy UD laminates produced by thermocompression[END_REF][START_REF] Moudood | Environmental effects on the durability and the mechanical performance of flax fiber/bio-epoxy composites[END_REF][START_REF] Bar | Influence of flax/polypropylene distribution in twistless thermally bonded rovings on their composite properties[END_REF][START_REF] Pisupati | Influence of resin curing cycle on the characterization of the tensile properties of flax fibers by impregnated fiber bundle test[END_REF][START_REF] Khalfallah | Innovative flax tapes reinforced Acrodur biocomposites: A new alternative for automotive applications[END_REF][START_REF] Bensadoun | Fatigue behaviour assessment of flax-epoxy composites[END_REF]. In these studies, the mechanical properties (bending, impact, tensile) are identified at the scale of composite samples, but also at the fibers' scale using back-calculation method from the mixture law, and/or by the IFBT (Impregnated Fiber Bundle Text) method described in [START_REF] Bensadoun | Impregnated fibre bundle test for natural fibres used in composites[END_REF]. The use of UD reinforcements is justified by the fact that it is challenging to manage the alignment of reinforcement, at the different scales, in interlacing fabrics when using plant fibers. The main difficulties are associated with the yarn twist, which is necessary in textile technologies to manufacture interlaced fabrics. The twist induces a reduction of the ultimate fiber strength in the final composite by a factor cos²(2α) where α is the angle of twist [START_REF] Bar | Influence of flax/polypropylene distribution in twistless thermally bonded rovings on their composite properties[END_REF][START_REF] Moothoo | A study of the tensile behaviour of flax tows and their potential for composite processing[END_REF][START_REF] Shah | Modelling the effect of yarn twist on the tensile strength of unidirectional plant fibre yarn composites[END_REF][START_REF] Madsen | Plant fibre composites -porosity and stiffness[END_REF]. It can also restrict the polymer penetration into the reinforcing fiber bundles [START_REF] Liotier | Role of interface formation versus fibres properties in the mechanical behaviour of bio-based composites manufactured by Liquid Composite Molding processes[END_REF]. For fabric architectures with high crimp, the in-plane mechanical properties (e.g. stiffness and strength in tension and compression) can also be reduced [START_REF] Katnam | Towards balancing in-plane mechanical properties and impact damage tolerance of composite laminates using quasi-UD woven fabrics with hybrid warp yarns[END_REF]. Another difficulty when using interlaced fabrics is related to the prediction of their properties. The use of modified mixture laws, as proposed by different authors [START_REF] Virk | Modulus and strength prediction for natural fibre composites[END_REF][START_REF] Summerscales | Enhanced rules-of-mixture for natural fibre reinforced polymer matrix (NFRP) composites (comment on Lau et al. in volume 136)[END_REF] and taking into account the degree of alignment, requires the knowledge and measurement of the torsion. However, during the composite manufacturing stages, such as the preforming step [START_REF] Capelle | Complex shape forming of flax woven fabrics: Design of specific blank-holder shapes to prevent defects[END_REF], the injection molding step [START_REF] Pisupati | Influence of resin curing cycle on the characterization of the tensile properties of flax fibers by impregnated fiber bundle test[END_REF] or the impregnation according to the viscosity of the matrix [START_REF] Liotier | Role of interface formation versus fibres properties in the mechanical behaviour of bio-based composites manufactured by Liquid Composite Molding processes[END_REF][START_REF] Mbakop | Effect of compaction parameters on preform permeability and mechanical properties of unidirectional flax fiber composites[END_REF], particularly on complex shapes, the alignment of reinforcement is significantly modified. In addition, the drapability of fabrics on double-curvature molds cannot be achieved with UD preforms, to achieve this goal at least two directions of interlaced (or stitched) reinforcements must be present. So, one of the main advantage of the textile fabrics when compared to the UD tapes is effectively their very good drapability and their capability for complex shape formation with no gap. These elements allow the reduction of the manufacturing costs and times.

The development of interlaced reinforcement architectures based on natural fibers with controlled characteristics is therefore a requirement. This work is first dedicated to the manufacturing of woven fabrics made from hemp rovings for composite applications. A roving consists of a continuous alignment of fibers, parallel to each other and slightly twisted. In this work, rovings were preferred to yarns to manufacture woven fabrics. They have higher linear mass density, a higher number of fibers in cross section, and a lower twist level when compared to yarn. This increases the impregnation capability and they are therefore more suited for composite applications. As explained in previous studies [START_REF] Goutianos | The optimisation of flax fibre yarns for the development of high-performance natural fibre composites[END_REF][START_REF] Corbin | Multi-scale analysis of flax fibres woven fabrics for composite applications[END_REF][START_REF] Omrani | Mechanical properties of flax-fibre-reinforced preforms and composites: Influence of the type of yarns on multi-scale characterisations[END_REF][START_REF] Madsen | Hemp yarn reinforced composites -I. Yarn characteristics[END_REF], a sufficient tenacity is required for manufacturing preforms, which is not necessarily achieved with this type of yarn. In the textile industry, chemical treatments, mainly mercerization processes, are classically used to improve yarns properties and their weavability [START_REF] Ku | A review on the tensile properties of natural fiber reinforced polymer composites[END_REF][START_REF] Xie | Silane coupling agents used for natural fiber/polymer composites: A review[END_REF][START_REF] Baley | Specific features of flax fibres used to manufacture composite materials[END_REF]. Flax fibers are commonly used to manufacture reinforcement with low twisted rovings, which are available in the market and well-studied in the literature [START_REF] Capelle | Complex shape forming of flax woven fabrics: Design of specific blank-holder shapes to prevent defects[END_REF][START_REF] Corbin | Multi-scale analysis of flax fibres woven fabrics for composite applications[END_REF][START_REF] Omrani | Mechanical properties of flax-fibre-reinforced preforms and composites: Influence of the type of yarns on multi-scale characterisations[END_REF][START_REF] Baets | Determination of the optimal flax fibre preparation for use in unidirectional flax-epoxy composites[END_REF][START_REF] Goutianos | Development of Flax Fibre based Textile Reinforcements for Composite Applications[END_REF][START_REF] Yan | Flax fibre and its composites -A review[END_REF]. However, only few studies focused on composites reinforced by low twisted yarns manufactured with hemp rovings, and in these studies, hemp fibers are comingled or cowrapped with thermoplastic polymers [START_REF] Zhang | Commingled natural fibre/polypropylene wrap spun yarns for structured thermoplastic composites[END_REF][START_REF] Baghaei | Novel aligned hemp fibre reinforcement for structural biocomposites: Porosity, water absorption, mechanical performances and viscoelastic behaviour[END_REF][START_REF] Baghaei | Characterisation of polylactic acid biocomposites made from prepregs composed of woven polylactic acid/hemp-Lyocell hybrid yarn fabrics[END_REF][START_REF] Baghaei | Characterization of thermoplastic natural fibre composites made from woven hybrid yarn prepregs with different weave pattern[END_REF][START_REF] Baghaei | Manufacture and characterisation of thermoplastic composites made from PLA/hemp co-wrapped hybrid yarn prepregs[END_REF]. In this paper, woven reinforcements are developed with low twisted rovings made of 100% hemp fibers.

In this study, several characterization steps are performed at different scales (fiber, roving, fabric and composites). The performances of the composites made from these woven fabrics are compared to the ones obtained with a cross-ply flax laminate, considered to be the best plant fiber composite currently available on the market. The influence of the weaving step on the fiber features is also discussed.

Materials

Hemp rovings

The mono-dimensional reinforcement used in this study are hemp rovings supplied by Linificio e Canapificio Nazionale, an Italian Company. The rovings are 100% made of hemp fibers. They have a linear density of 259 Tex, measured according to the standard NF G07-316 [START_REF]NF G07-316 -Textiles -Tests of yarns -Determination of linear density[END_REF], and a twist level of 36 turns/m, measured according to the standard NF G07-079 [START_REF]NF G07-079 -Textiles -Testing threads -Determining the twisting of threads by untwisting/retwisting with a double re-test[END_REF]. These rovings were submitted to a chemical treatment classically used in textile industry, called mercerization, to improve their tenacity in view of the weaving [START_REF] Madsen | Plant fibre composites -porosity and volumetric interaction[END_REF].

Woven fabric manufacturing

Woven fabrics were manufactured with the above-mentioned hemp rovings on a manual dobby loom, a

Leclerc Weavebird weaving machine. The weaving process used is shuttle weaving. Before the weaving process, several preparations steps are required. The first one, called warping, consists in preparing the warp beam: warp yarns are wound parallel to each other, with the length of the future fabric, around the drum of the warping machine and then they are transferred on the warp beam. The further step is the drawing in: once the warp beam is fixed on the dobby loom, the warp yarns are inserted one by one through the heddles and the weaving reed. The weaving can then begin: the weave pattern is created with the movement of the dobby frames and between to movements, the shuttle is passed inside the shed to insert the weft yarn and the weaving reed is refolded to compress and form the fabric. Three fabrics with different weaving diagrams were produced, i.e. plain weave (H1), satin 6 weft effect (H2) and twill 6 weft effect, see Figure 1. On schematic views, made with WiseTex software, warp yarns are represented in blue color and weft yarns in red color. They have the same warp density of 6 yarns/cm and different weft densities, namely 5 yarns/cm for H1 and 9.5 yarns/cm for H2 and H3. Warp and weft densities have been measured according to the standard NF ISO 4602 [START_REF] Afnor | NF ISO 4602 -Reinforements -Woven fabrics -Determination of number of yarns per unit length of warp and weft[END_REF]. Composites plates were manufactured by hot pressing with the produced hemp fabrics and an epoxy resin, the GreenPoxy 56 from Sicomin Company. Their dimensions are 200 x 300 mm². The fabrics were conditioned at a temperature of 23°C and a relative humidity of 50% during at least 24 h prior to the composite manufacturing. After this pre-conditioning, the plies were hand layed-up, impregnated manually with the epoxy resin and then cured at 130°C under a pressure of 3 bars during 1 h. A hot-press Agila Press 100 kN was used. As there are three types of fabrics, three types of composites were manufactured (noted as C1, C2 and C3 with the H1, H2 and H3 reinforcements, respectively). Four plies were stacked for C1 and C2 composites and two plies for C3. All the plies were stacked with the same orientation. During composite manufacturing, there is 1.6 times the mass of resin compared to the mass of fibers.

The mold is open at the ends, which allows the excess resin to be removed. Following the manufacturing, the composite plates were conditioned at 23°C and 50% RH during at least four weeks to reach the moisture content equilibrium. In order to compare the properties of the composites reinforced by the woven fabrics, a "reference" composite was also manufactured with a pure unidirectional flax reinforcement. The FlaxTape TM 110 by Lineo Company was used to manufacture a cross-ply laminate with the same protocol.

IFBT specimens were also prepared by aligning the hemp rovings and impregnating them with GreenPoxy 56 resin. After impregnation, the specimens were cured at 60°C during 24 under a pressure of 2 bars. The dimension of the manufactured specimens are approximately 200 mm x17 mm x 1 mm. Prior to the manufacturing of the IFBT specimens, the hemp rovings were conditioned at a temperature of 23°C and a relative humidity of 50% during few hours to reach the equilibrium. After manufacturing, the IFBT specimens were also stored in a climatic chamber (23°C -50% RH) during a minimum of three weeks to reach the moisture content equilibrium.

For both composite plates and IFBT specimens, the quantity of resin was calculated to reach a fiber volume fraction in composites of approximately 40%. The real fraction of each specimen and each plate was measured after manufacturing.

Characterization methods

Fiber properties -Impregnated Fiber Bundle Test (IFBT)

The IFBT test was used to determine by inverse method the effective tensile properties of the hemp fibers constituting the rovings. In this method, fully described in [START_REF] Bensadoun | Impregnated fibre bundle test for natural fibres used in composites[END_REF], rovings are aligned in the same direction and impregnated with the resin to manufacture composite samples. Tensile tests are then done on these samples, and the effective modulus of fibers is obtained by back-calculation using a rule of mixtures. In this work, we used the following law [START_REF] Shah | Modelling the effect of yarn twist on the tensile strength of unidirectional plant fibre yarn composites[END_REF][START_REF] Madsen | Plant fibre composites -porosity and stiffness[END_REF][START_REF] Virk | Modulus and strength prediction for natural fibre composites[END_REF][START_REF] Summerscales | Enhanced rules-of-mixture for natural fibre reinforced polymer matrix (NFRP) composites (comment on Lau et al. in volume 136)[END_REF][START_REF] Madsen | Plant fibre composites -porosity and volumetric interaction[END_REF][START_REF] Summerscales | A review of bast fibres and their composites. Part 2 -Composites[END_REF][START_REF] Aslan | Effect of consolidation pressure on volumetric composition and stiffness of unidirectional flax fibre composites[END_REF], fully described by Eqs. 1-3:

= + 1 - (1) = cos (2 ) (2) = tan (2 ) (3)
where E C is the composite modulus, η 0 a fiber orientation factor, η 1 the fiber length factor, V f the fiber volume fraction, V m the matrix volume fraction, V p the porosity volume fraction, E f is the fiber modulus (in the longitudinal direction), E m the matrix modulus, n the porosity factor, α the surface twist angle, r the radius of roving and T the twist level of roving.

Taking into account the non-linear tensile behavior generally observed for plant fiber composites, and for comparative purpose, the composite modulus was measured on two different strain ranges: E 1C between 0 and 0.1% of longitudinal strain and E 2C between 0.3 and 0.5%. The corresponding modulus determined by back-calculation at the scale of the fibers were noted E f1 and E f2 . The fiber length factor (η 1 ) was considered equal to 1 (which is generally the case when the length to diameter ratio of the fibers is higher than 50). Based on the empirical results reported in literature [START_REF] Madsen | Plant fibre composites -porosity and stiffness[END_REF][START_REF] Aslan | Effect of consolidation pressure on volumetric composition and stiffness of unidirectional flax fibre composites[END_REF], the porosity factor (n) was considered equal to 2.

Tensile tests were conducted on five IFBT specimens using a MTS Criterion 45 universal machine, with a crosshead displacement rate of 1 mm/min. The longitudinal strain was measured with an Instron 2620-601 extensometer with a gauge length of 50 mm.

Roving properties

At the roving scale, the twist level was studied according to the standard NF G07-079 [START_REF]NF G07-079 -Textiles -Testing threads -Determining the twisting of threads by untwisting/retwisting with a double re-test[END_REF] on 10 samples per type of yarn. The tenacity was measured using tensile tests. For each yarns, 10 specimens were tested using a MTS Criterion 45 tensile machine, with a load cell of 10 kN, a gauge length of 200 mm, a crosshead displacement rate of 200 mm/min and a preload of 0.5 cN/Tex, according to the standard NF EN ISO 2062 [START_REF] Afnor | EN ISO 2062 -Textiles -Yarns from packages -Determination of single-end breaking force and elongation at break using constant rate of extension[END_REF].

The same properties were also determined after weaving to evaluate the influence of the weaving process on these features. The rovings were manually extracted from the fabrics.

Fabric properties

The areal density and the thickness of the different fabrics were measured on 10 specimens according to the standards NF EN 12127 [START_REF]NF EN 12127 -Textiles -Fabrics -Determination of mass per unit area using small samples[END_REF] and NF EN ISO 5084 [START_REF]-Textiles -Determination of thickness of textiles and textile products[END_REF], respectively.

Compaction tests were also performed on woven fabrics, on three specimens for each fabric, to characterize and compare the compressibility behavior of the different weave patterns and to determine the fiber volume fraction as a function of applied pressure. This is an important parameter for the composite manufacturing processes. The compaction tests were performed on a MTS Criterion 45

universal machine with a full-range load cell of 100 kN at a crosshead displacement rate of 1 mm/min.

The pressure was applied with a compaction cylinder with a diameter of 62.5 mm until a pressure of 20

MPa. The displacement was measured with a LVDT sensor, mounted between the compaction cylinder and the bottom plate, with a full range of 2.5 mm and an accuracy of 2.5 µm. The compressive frame stiffness of the test set-up was also measured and taken into account to correct the measured displacement, as suggested in [START_REF] Ivanov | Local compressibility of draped woven fabrics | Request PDF[END_REF][START_REF] Danzi | A model for the time-dependent compaction response of woven fiber textiles[END_REF][START_REF] Mbakop | Experimental analysis of the planar compaction and preforming of unidirectional flax reinforcements using a thin paper or flax mat as binder for the UD fibers[END_REF]. Square samples of the fabrics with dimensions of approximately 80x80 mm² were laid-up in the same material direction (see Figure 3). For all the tests, four plies were used. The initial thickness of the lay-ups was measured on the machine with a pre-load of 10 N. Three successive loading/unloading steps were then applied. Only the first loading is analyzed in this study. The pressure is calculated by dividing the applied force by the area of the compaction cylinder, and the strain by dividing the measured displacement by the initial thickness.

Composite properties

The microstructure of the manufactured composites was characterized using X-ray nanotomography.

Image acquisition was realized on an Easytom micro-nano tomograph (RX Solutions, France). The system is equipped with an X-ray source Hamamatsu Open Type Microfocus L10711 having maximum voltage 160 keV. The images reported herein were collected using X-rays of voltage 40 keV and a current of 22 µA. The X-ray transmission images were acquired using a Fluoroscopic High Speed imaging subsystem PaxScan 2520DX. The entire volume is reconstructed at a full resolution with a voxel size of 5.77 µm, using filtered back-projection. The data were processed using VG Studio software.

Tensile specimens were cut in the manufactured composite plates using a Trotec laser cutter. The mechanical characterization was conducted according the standard ASTM D3039-00 [START_REF]ASTM D3039-00 -Test Method for Tensile Properties of Polymer Matrix Composite Materials[END_REF], on a MTS Criterion 45 universal machine, with a crosshead displacement rate of 1 mm/min, a specimen width of 15 mm and a gauge length of 150 mm. For samples per direction (warp and weft) are tested for all the composites materials. The longitudinal strain was measured with an Instron 2620-601 extensometer with a gauge length of 50 mm.

Results and discussions

Influence of weave pattern on the fabric properties

The textile properties measured on the woven fabrics are presented in Figure 4.a.for the areal density and in Figure 4.b. for the thickness and weft density. H2 and H3 fabrics have a similar areal density of approximately 400 g.m -2 , which is approximately 1.4 times higher than the one of H1. This can be explained by the interlacement and the weft density. Indeed, the H2 and H3 fabrics, which are satin 6 and twill 6 weave patterns, have one binding point for 6 rovings, whereas the H1 fabric (plain weave pattern) has one binding point for two rovings. A binding point is the interlacement between warp and weft rovings, represented by black squares in Figure 1. Thus, as rovings are less interlaced in H2 and H3 fabrics, they are more compacted on top of each other and the weft density is higher. As a result, there are more rovings per unit area in H2 and H3 fabrics, and the areal density is higher than for H1 fabric. The compressibility of the different fabrics was also characterized. The compaction curves are shown in Figure 5. The three fabrics exhibit a very similar behavior with a maximum fiber volume fraction of approximately 70%. In this way, it can be expected that the composite made with these three fabrics will have similar fiber volume fraction after manufacturing steps.

The compaction response is strongly non-linear (see Figure 5). As explained by different authors [START_REF] Nguyen | Analyse expérimentale et numérique de la compaction des renforts fibreux: Application pour la perméabilité[END_REF][START_REF] Grieser | Investigation of the compaction behavior of carbon fiber NCF for continuous preforming processes[END_REF] this results from the expression of different phenomena. In the first stage, for low pressure, the

Thickness

Weft density response is quasi-linear. It results from the combination of the deformation of the roving cross-section, the roving flattening and nesting. In the second stage, the curve becomes nonlinear and in addition to the first mechanisms, the bending of the roving takes place, which leads to the reduction of void and slippage between fiber and/or roving. At last, a third phase, linear again, is observed. The compaction is then attributed to the transverse compression of fibers

Influence of the weaving process on the roving features and fiber orientation

During the weaving process, the roving properties are impacted. The values of the twist level and of the breaking tenacity of rovings before and after weaving for each fabric are shown in Figure 6.a and Figure 6.b, respectively. In the warp direction, the values of twist level remains in the standard deviation whereas in the weft direction, the twist level decreases significantly after weaving. This untwisting is more important for the weft rovings inside H2 than in H3 and H1 fabrics. A decrease of approximately 94% in the twist level is measured for the roving from the H2 fabrics against a decrease of 61% for the roving from H3 fabric and 51% for the roving from H1 fabric. This trend is similar for the decrease of breaking tenacity, always in weft direction: there is a decrease of approximately 46% for H2 fabric, 34% for H3 fabric and 16% for H1 fabric. This is consistent with the decrease in twist level. When the rovings are less twisted there is less cohesion between fibers and then the strength is reduced. The decrease in the twist level of weft rovings is attributed to the technique used for the weft insertion in shuttle weaving. Whereby, the weft rovings need to be wound into a small bobbin, called a quill, which is inserted in the shuttle to pass the weft roving through the shed.

In the warp direction, even if the standard deviations are quite high, an increase of tenacity at break after weaving can be observed. This increase is more important for the H2 fabric. For the satin weave pattern, the movement of frames moves rovings which are distant from each other, whereas for plain weave and twill weave patterns, rovings are raised side by side. Indeed, there are less fibers hanging between the warp rovings during the satin weaving process, and they are less "damaged" in the process, which is why they have the best tenacity at break. The warping process and the tension applied on the rovings during the weaving process can stiffen the warp rovings, and confer to them better mechanical properties, especially a better tenacity at break, which is more or less damaged by hanging due to the movement of the frames.

Figure 6 : Influence of weaving process on yarns properties a. Twist level b. Tenacity at break

The values of the twist level and tenacity at break before and after weaving are summarized, in Figure 7, to highlight the influence of twist level on tenacity at break. As mentioned in the literature [START_REF] Goutianos | The optimisation of flax fibre yarns for the development of high-performance natural fibre composites[END_REF][START_REF] Omrani | Mechanical properties of flax-fibre-reinforced preforms and composites: Influence of the type of yarns on multi-scale characterisations[END_REF], the increase in the twist level leads to an increase in the breaking tenacity. The twist level gives more cohesion between the fibers and thus more strength to the rovings.

Warp yarns before weaving Warp yarns after weaving

Weft yarns before weaving Weft yarns after weaving The tensile properties of the fibers were also determined before and after weaving using IFBT tests (Figure 8). Before weaving, the rovings were extracted from the bobbins, and after weaving, the rovings are extracted in warp and weft direction from the fabrics. Figure 8 shows the values of the fiber moduli (E f1 and E f2 ) before and after weaving for the three fabrics. These results show the high rigidity of the hemp fibers used for this work, with an initial apparent rigidity (E f1 ) greater than 60 GPa.

In the weft direction, the fiber modulus E f1 after weaving is slightly lower than the one before weaving for the H1 and H3 fabrics. This apparent stiffness reduction could be both attributed to the remaining crimp in the rovings (and the resulting fiber disorientation, neglected when back-calculating the fiber properties) as well as fiber damage induced by the weaving process. Interestingly, the mean E f1 value for the H2 fabric is higher after weaving. This improvement of the fiber stiffness can be attributed to a stiffening effect [START_REF] Placet | Nonlinear tensile behaviour of elementary hemp fibres. Part I: Investigation of the possible origins using repeated progressive loading with in situ microscopic observations[END_REF] induced by the mechanical loading applied during the weaving process. This effect could counterbalance the potential decrease in fiber stiffness due to their misorientation (remaining crimp effect) and hypothetic damage induced by the weaving. All these underlying mechanisms and microstructural features deserve to be thoroughly studied in forthcoming studies. In the warp direction, the fiber stiffness is more negatively impacted by the weaving process. A decrease of approximately 30% is observed after weaving for the H1 and H3 fabrics. Again, this effect is less pronounced for the H2 fabric. As a conclusion, it can be underlined that the fiber stiffness is well preserved during the processing

Influence of the weave pattern on the composite properties

Using these three fabrics (H1, H2 and H3), composite materials were manufactured, respectively noted CH1, CH2 and CH3. The thickness and the microstructural features are shown in Table 1. Two methods were used to determine the constituent volume fractions. The first one is by weighting, as described in [START_REF] Berges | Influence of moisture uptake on the static, cyclic and dynamic behaviour of unidirectional flax fibre-reinforced epoxy laminates[END_REF]. The second one is by using X-ray nanotomography images analysis and reconstruction. In order to conduct observations inside the composite materials, nanotomography imaging is used. Some images are presented in Figure 9, Figure 10 and Figure 11, respectively for CH1, CH2 and CH3 composites. For the three composites, some porosity can be observed, which is mainly located between the plies but also inside the rovings. However, there is not a lot of porosity, and none of the composites have major impregnation defects. The process used is thus suitable for the manufacture of the composite plates. With nanotomography analysis, it is also possible to study the distribution of porosities throughout the thickness of the composite samples (see Figure 12). For CH1 and CH2, porosities are evenly distributed over all plies whereas for CH3, porosities are more located in the upper ply (as a reminder, CH3 is made of only two plies).

Table 1: Composites composition and thickness

Composites

Figure 12: Distribution of porosities through the thickness of the composite samples

The mechanical properties of the three composites are presented in Figure 13.a. for the modulus and in Figure 13.b for the strength. Since the fiber volume fraction was slightly different for each composite type, and to facilitate the comparison, the properties were normalized with a fiber volume fraction of 45%, using a rule of mixture. It can be observed that the satin weave pattern leads, at the scale of the composite material, to the best mechanical properties in the weft direction. In the weft direction, the modulus is 10% higher than in the warp direction for CH1, 44% higher for CH2 and 31% higher for CH3. For the strength at failure, the value in the weft direction is 20% higher than in the warp direction for CH1, 49% for CH2 and 41% for CH3. Indeed properties are more similar between the warp and weft direction in CH1 composite than in the CH2 and CH3 composites. As expected, the plain weave pattern 0 0.2 0.4

-1.5 -1 -0.5 0 0.5 1 1.5 Porosity cumulate (%)

Thickness (mm) CH1 (1.8 mm) CH2 (2.4 mm) CH3 (1.2 mm)
gives the better balanced properties to its composite. In the H2 and H3 fabrics, the weft density is higher than in the H1 fabric due to the weave pattern, there are also more fibers in this direction, and the resulting mechanical properties are higher in the weft direction. Moreover, in the weft direction, the rovings are untwisted during the weaving process, as shown previously in Figure 6.a., and thus the fibers are more aligned with the loading direction. In this way, the mechanical properties are higher in the weft direction than in the warp direction, thanks to a higher fiber volume fraction and a better alignment of the fibers in the direction of the mechanical stress.

The mechanical properties of these composites (CH1, CH2 and CH3) made from a woven fabric were also compared to the one obtained for a cross-ply laminate stacking sequence, consisting of 0° and 90° purely unidirectional flax fiber plies. In terms of mechanical properties, the composite reinforced by unidirectional plies has a tensile strength of 182 MPa and a modulus of 18.7GPa in both directions. In comparison with the reference composite, the modulus in the warp direction of each woven composite is lower (-23% for CH1, -40% for CH2 and CH3) but quite similar or higher in the weft direction (-14% for CH1, -7% for CH2 and -12% for CH3). This trend is similar for the strength at failure: the values in the warp direction are lower than for the reference composite (-30% for CH1, -47% for CH2 and -46% for CH3) and similar higher in the weft direction (-3% for CH1, 17% for CH2 and 5% for CH3). By using a plain weave pattern for woven fabric reinforcement, the properties obtained are closer than those which can be achieved with a symmetrical cross-ply laminate. The use of a more complex weave pattern, such as satin or twill, leads to more unbalanced properties, but which can be higher in one direction than those of a composite manufactured with unidirectional plies. These results are consistent with what had already been demonstrated in another study with hemp/PLA comingled rovings, in which composite reinforced by satin woven fabric has higher properties than the one reinforced by basket woven fabric [START_REF] Baghaei | Characterization of thermoplastic natural fibre composites made from woven hybrid yarn prepregs with different weave pattern[END_REF].Thus, it can be assumed than these results are promising, and high performance composites can also be manufactured with woven fabrics made of 100% hemp low twisted rovings as reinforcements.

Influence of rovings features on composite properties

As seen previously, the roving properties are influenced by the weaving process. The main parameter altered by the weaving process is the twist level. This factor is really influential on the fiber alignment in composite materials and thus on the mechanical properties which result at the scale of the composite material. The percentage of fibers in the two main directions was determined from the tomography analysis. The values are determined from a single sample for each composite and shown in Table 2. The total percentage of fibers in warp and weft directions is lower than 100% due to misalignments of fibers mainly due to twisting of rovings. By twisting the rovings, fibers are positioned in other directions than the main direction of the rovings, and thus, after weaving, they are oriented in other directions than the warp or weft direction of the fabrics. In this way, the twist level of rovings after weaving (Figure 6.a.)

was compared with the percentage of fibers in the main directions of the composite (Table 2). It can be conclude that the more the twist level of rovings after weaving is decreased during the weaving process, the higher is the fiber content in the composite direction. By decreasing the twist level of the roving during weaving process, the fibers are more aligned in the direction of the mechanical stress.

By plotting the strength of the composites measured in both material directions for the different woven fabrics according to the twist level of the rovings (measured after weaving), we confirm a direct relationship already mentioned in the literature: the more the fibers are aligned in the roving main axis and the load direction, the higher the strength of the composite [START_REF] Shah | Modelling the effect of yarn twist on the tensile strength of unidirectional plant fibre yarn composites[END_REF][START_REF] Goutianos | The optimisation of flax fibre yarns for the development of high-performance natural fibre composites[END_REF][START_REF] Omrani | Mechanical properties of flax-fibre-reinforced preforms and composites: Influence of the type of yarns on multi-scale characterisations[END_REF].

H3 Weft

The properties of the composites were also predicted using a V f equal to 45%, the rule of mixtures, and the back-calculated properties of the fibers, their twist level in the rovings and the properties of the matrix. The results are presented in Figure 15. The rule of mixture used is presented in Eqs. 1-3. In ROM1, the values of the fiber modulus and the twist level correspond to the ones measured before weaving, whereas in ROM2, the values after weaving are considered. The moduli of fibers E f2 before and after weaving used are the ones presented in Figure 8. The composite properties obtained using ROM1 are clearly under-or overestimated for all materials: approximately 10% of difference between experimental and predicted values in both directions. However, by using the fiber and roving properties after weaving, the composite properties are predicted with a better agreement when compared to experimental data (the relative difference between experimental and predicted values varies from 1% to 6%) excepted for warp direction in CH1 and CH3 composites (respectively 36% and 24% of difference between experimental and predicted values). To obtain a good prediction, the fiber and roving properties after processing have to be considered. One of the remaining questions is the influence of the composite manufacturing (in particular pressure and temperature) on the fiber and roving features. This will be studied in a further step using the X-ray nanotomography results.

Table of figures

Table of tables

Figure 1 :

 1 Figure 1: Schematic representation and pictures of the woven fabrics manufactured and studied in this work

Figure 2 :Figure 3 :

 23 Figure 2: Compaction test setup a. schematic diagram with preforms b. configuration for machine stiffness measurement

Figure 4 :

 4 Figure 4: Properties of the woven fabrics a. Areal density b. Thickness and weft density

Figure 5 :

 5 Figure 5: Compaction curves for the three hemp fabrics

Figure 7 :

 7 Figure 7: Breaking tenacity of the rovings as a function of their twist level. Evidence of a direct relationship

  fabric, and even improved for the fibers oriented in the weft direction. All the results abovementioned are also valid for E f2 .

Figure 8 :

 8 Figure 8: Comparison of the fiber stiffness before and after weaving for the three considered fabrics

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Nanotomography views of CH1 composite

Figure 13 :

 13 Figure 13: Composite mechanical properties estimated at V f = 45% (Reference = flax cross-ply laminate) a. Modulus b. Strength

Figure 14 :

 14 Figure 14: Relationship between the twist level of the rovings in the fabrics and the strength of composites

Figure 15 :

 15 Figure 15: Comparison of experimental modulus and predicted modulus of composites 5 Conclusion

Figure 1 : 4 Figure 2 : 7 Figure 3 : 8 Figure 4 : 9 Figure 5 : 10 Figure 6 : 12 Figure 7 : 13 Figure 8 : 14 Figure 9 : 15 Figure 10 : 15 Figure 11 : 15 Figure 12 : 16 Figure 13 : 18 Figure 14 : 19 Figure 15 :

 142738495106127138149151015111512161318141915 Figure 1: Schematic representation and pictures of the woven fabrics manufactured and studied in this work .............................................................................................................................................................. 4 Figure 2: Compaction test setup a. schematic diagram with preforms b. configuration for machine stiffness measurement .................................................................................................................................. 7 Figure 3: Schematic view of the stacking sequence in compaction tests (weft direction in red, warp direction in blue) .......................................................................................................................................... 8 Figure 4: Properties of the woven fabrics a. Areal density b. Thickness and weft density .......................... 9 Figure 5: Compaction curves for the three hemp fabrics ............................................................................ 10 Figure 6 : Influence of weaving process on yarns properties a. Twist level b. Tenacity at break .............. 12 Figure 7: Breaking tenacity of the rovings as a function of their twist level. Evidence of a direct relationship ................................................................................................................................................. 13 Figure 8: Comparison of the fiber stiffness before and after weaving for the three considered fabrics ..... 14 Figure 9: Nanotomography views of CH1 composite ................................................................................ 15 Figure 10: Nanotomography views of CH2 composite .............................................................................. 15 Figure 11: Nanotomography views of CH3 composite .............................................................................. 15 Figure 12: Distribution of porosities through the thickness of the composite samples .............................. 16 Figure 13: Composite mechanical properties estimated at V f = 45% (Reference = flax cross-ply laminate) a. Modulus b. Strength ............................................................................................................................... 18 Figure 14: Relationship between the twist level of the rovings in the fabrics and the strength of composites .................................................................................................................................................. 19 Figure 15: Comparison of experimental modulus and predicted modulus of composites .......................... 20

Table 2 :

 2 Percentage of fibers in composite in two main material directions

	Composite Warp direction (%)	Weft direction (%) Total (%)

Table 1 :

 1 Composites composition and thickness ........................................................................................ 14

Table 2 :

 2 Percentage of fibers in composite in two main material directions .............................................. 18

Acknowledgements

The authors wish to thank the Italian Company "Linificio e Canapificio Nazionale" for providing the hemp rovings used in this study. This project has received funding from the Bio-Based Industries Joint Undertaking under the European Union's Horizon 2020 research and innovation program under grant agreement No 744349.

they can endure higher in-plane shear stress, in comparison with unidirectional reinforcements. In this study, the influence of weaving process parameters, and especially the weave pattern, has been studied at different scales. Compaction tests have also been conducted on woven fabrics and do not highlight significant differences between the weave patterns, thus the fiber volume in the final composite can be expected to be similar for the three materials. Afterwards, the influence of weaving process on rovings and fibers properties has been studied. Due to the untwisting of rovings used in weft direction during weaving process, cohesion between fibers in the cross section becomes less important, which reduces the tenacity at break of the rovings. Despite this untwisting of the rovings, after weaving, the effective stiffness of fiber (studied through IFBT method) is decreased in most cases. This effect is lower when weaving a satin pattern. A fiber stiffening is even observed in the weft direction. Whatever the weave pattern used, only few pores can be observed in composites through tomography analysis. Regarding the mechanical properties of the composites, the use of several weave diagrams leads to balanced or unbalanced properties between the two main directions, depending on the yarn density and the fiber fraction in the studied direction. The percentage of fibers aligned in the direction of the mechanical stress was directly linked to the untwisting of the rovings in this direction during the weaving process. The determination of the effective fiber properties after weaving allows also better composite properties prediction through the rule of mixtures. Results at the scale of composite highlights that mechanical properties similar to the best flax cross-ply laminates can be achieved with hemp woven fabrics made from low-twisted rovings.

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Authors' names and affiliations


Safiullah Timorian
email: safiullah.timorian@gmail.com

Morvan Ouisse

Noureddine Bouhaddi

Sergio De Rosa

Numerical investigations and experimental measurements on the structural dynamic behaviour of quasi-periodic meta-materials

Keywords: Quasi-periodic structures, band diagram, frequency stop-band, forced response, metamaterials

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

Numerical investigations and experimental measurements on the structural dynamic behaviour of quasi-periodic meta-materials S. Timorian 1,2 , M. Ouisse 2 , N. Bouhaddi 2 , S. De Rosa 1 and F. Franco 1

Introduction

Recent literature review shows the design of metamaterials as one of the central topics in the vibroacoustic analysis of periodic media. There are interesting subjects in the literature linked to the current investigation which is intended mainly to the design of new architectured metamaterials (i.e. lightweight, easy-to-manufacture, and low frequency stop band properties, noise and vibration reduction) [START_REF] Matlack | Composite 3dprinted metastructures for low-frequency and broadband vibration absorption[END_REF], [START_REF] Claeys | Design and validation of metamaterials for multiple structural stop bands in waveguides[END_REF].

Periodic structures create stop bands effect due to the geometrical and/or material impedance mismatches that can result in great vibroacoustic performances. In terms of material properties, there are various cases that show in-plane and out-of-plane elastic properties. For example the dynamical behaviour of a 2D periodic waveguide, which exhibits in-plane elastic properties (Young's and shear modulus) compared to out-ofplane ones, are described in terms of elastic wave propagation in [START_REF] Ouisse | A piezo-shunted kirigami auxetic lattice for adaptive elastic wave filtering[END_REF][START_REF] Collet | Floquet-bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems[END_REF]. The possibility of designing such smart materials or structures, that can partially reduce mechanical waves on certain frequency ranges, is addressed in [START_REF] Billon | Mechanics and band gaps in hierarchical auxetic rectangular perforated composite metamaterials[END_REF][START_REF] Madeo | Modeling phononic crystals via the weighted relaxed micromorphic model with free and gradient micro-inertia[END_REF][START_REF] Martinsson | Vibrations of lattice structures and phononic band gaps[END_REF]. Geometrical discontinuity also plays a central role in creating stop band effects; as an example, sandwich beams with periodic auxetic core, exhibiting impedance mismatch generated by varying elastic and mass properties of the core, are able to produce stop band effects [START_REF] Ruzzene | Wave propagation in sandwich plates with periodic auxetic core[END_REF][START_REF] Qing-Tian | Wave propagation in sandwich panel with auxetic core[END_REF].More complex shapes like star shaped honey-combs can be designed to provide frequency stop band behaviour in low frequency regimes [START_REF] Meng | Band gap analysis of star-shaped honeycombs with varied poisson's ratio[END_REF].

The quasi-periodicity concept which is investigated in this paper can be interpreted as a certain degree of irregularity introduced in a periodic pattern. According to the literature, there are numerous examples of disordered periodicities and/or uncertainties in real structures, like bridges with column spans, and array of fuel tanks interconnected with each other by flexural links. The issues of the non-perfect periodicity in real structures can also be attributed to errors in manufacturing processes [START_REF] Carta | Transmission and localisation in ordered and randomly-perturbed structured flexural systems[END_REF]. The presence of defects and imperfections in geometric and constitutive properties of the structures is generally referred as disorder [START_REF] Chen | Study on band gaps of elastic waves propagating in onedimensional disordered phononic crystals[END_REF]. Numerical solutions of these types of disordered systems can lead to the need of full stochastic approaches [START_REF] Ichchou | Stochastic wave finite element for random periodic media through first-order perturbation[END_REF]. When the irregularity is localized in space (a 'defect in the periodic arrangement), phenomenon of Anderson localization may occur, evidencing that the vibration propagation in a structure is not entirely regular and that could be impeded by the irregularities, giving rise on the average to an exponential decay of vibration level [START_REF] Pierre | Localization of vibrations by structural irregularity[END_REF][START_REF] Hodges | Vibration isolation from irregularity in a nearly periodic structure: theory and measurements[END_REF]. They also demonstrates an example of a string with regular and irregular spacing of added masses in order to use it as a sort of passive vibration control.

On the other hand, the imperfections can be engineered and used as design parameters to tailor the dynamic behaviour. Among others, sequences of impedance mismatches built on numerical series like Fibonacci, Thu-Morse or Rudin Shapiro can be considered as design templates for the engineered irregularity [START_REF] Velasco | Elastic waves in quasiperiodic structures[END_REF].

Hou et al., [START_REF] Hou | Acoustic wave propagating in one-dimensional fibonacci binary composite systems[END_REF], stated that robust quasi-periodic design may offer new vibroacoustic properties to the structures. This robust design refers to the investigation of the transmission properties and the frequency spectra of a Fibonacci composite material with different thickness ratio of two layers. Similarly Aynaou et al. investigated theory of acoustic wave propagation on 1D phononic band gap structures made of slender tube loops pasted together with slender tubes of finite length following Fibonacci sequence [START_REF] Aynaou | Propagation and localization of acoustic waves in fibonacci phononic circuits[END_REF]. In another example similar topological configuration of Fibonacci and Thue-Morse sequences are investigated with an experimental observation of the phononic scattering band structure in 1D quasi-periodic systems [START_REF] King | Acoustic band gaps in periodically and quasiperiodically modulated waveguides[END_REF]. A 1D bar is designed using Fibonacci sequence for stop band distribution and prestress effects by [START_REF] Gei | Wave propagation in quasiperiodic structures: stop/pass band distribution and prestress effects[END_REF][START_REF] Pisano | Fibonacci's liber abaci: a translation into modern english of the book of calculation[END_REF].

The current work is based on the analysis of quasi-periodic 2D structures. A deterministic approach is taken into account to introduce a two directional Thue-Morse sequence for creating a quasi-periodic meta-material to improve the vibroacoustic performances.

In the literature there are various examples of 1D quasi-periodic systems for the dynamic analysis; on the contrary, 2D quasi-periodic systems are very rare. Some sequences are intrinsically ready to be used for 2D cases as the Thue-Morse, while others like Fibonacci are well adapted to the 1D cases.

In this paper, two dimensional lattices are built with the conventional finite element method (FEM) to comply with the Thue-Morse sequence in order to explore the opportunities offered in terms of reduction of the forced response in some frequency bands [START_REF] Timorian | Investigation for the analysis of the vibrations of quasiperiodic structures[END_REF][START_REF] Timorian | Band diagram and forced response analysis of periodic and quasi-periodic panels[END_REF]. Specifically, the work presents modified versions of star-shaped concave unit cells.

The tailored quasi-periodicity is defined by invoking a bi-directional Thue-Morse Morphism sequence on the meta-material. The geometrical impedance mismatch results in asymmetry and follows a combination of two different star-shape elements by variation of different corner angles. The main target of the study is the design of structural stop band/filters in order to isolate maximum vibration level in some frequency ranges.

The contents of this paper is structured as follows. In section 2, quasi-periodic structure is presented. Section 3 introduces the Wave Finite Element Models: dispersion analysis. On the contrary, in section 4, experimental results are validated. Section 5 describes the structural dynamic analysis using Finite Element models. A closing section summarise the conclusion of this work.

Quasi-periodic structures

In this paper the focus is on the application to 2D domains of the Thue-Morse sequence for introducing irregularity in both directions. In this section a possible approach is presented to analyse the vibration performances of quasi-periodic 2D structures.

Quasi-periodic structure section contains three sub-sections. The first one is dedicated to the cells description, followed by the second sub-section that introduces the Thue-Morse Morphism sequence, while the third sub-section concludes with the periodic sequence.

Cells

The waveguide is a star-shaped lattice. The meta-material is made of two different unit cells (A and B). Fig. 1 shows the geometry of the two cells inspired from [START_REF] Chang | Elastic wave propagation of twodimensional metamaterials composed of auxetic star-shaped honeycomb structures[END_REF].

The difference between the two unit cells is in the angles α, which defines the corner opening in the concave shape.

Unit cells A and cell B have a side of 22 mm, and the internal core with star shape is 20 × 20mm 2 . The thickness of the base plate is 3 mm and the wall thickness of the cells is 1mm. Cell A has an angle α A = 30 • and cell cell B has an angle α B = 20 • . The angles are variables of the sequences according to the generating order of the developed

Thue-Morse Morphism sequence

The Thue-Morse sequence is based on binary arrays. Starting from a 1D example, the sequence can be derived by a function s from the set of binary sequences to itself by replacing every A in a 1D sequence with AB and every B with BA [START_REF] Berstel | Combinatorics on words: Christoffel words and repetitions in words[END_REF]:

s(x, y) = s(x)s(y). (1)
Eq.1 defines a map s for all strings x, y. Defining the Thue-Morse morphism s(A) = AB, leads to the relations for increasing the quasi-periodicity pattern in one direction. The associated sequence starts with A then AB, ABBA, ABBABAAB, and so on, [START_REF] Allouche | The ubiquitous prouhet-thue-morse sequence[END_REF].

Considering the sequence in 2D, the first order starts from a 2 × 2 matrix containing [A, B; B, A] which is equivalent to the second order s 2 (A). In the next steps, second order matrix is translated to the right, left, and diagonal directions in order to built s 3 (A) and so on. A loop is generated to translate the same order to the right and left with the opposite entries and similar entries to the diagonal direction to fulfil the pattern of the Thue-Morse sequence. The orders are detailed as follows:

s(A) = AB s 2 (A) = s(s(A)) = AB | BA s 3 (A) = s(s(s(A)) = ABBA | BAAB

Wave and finite element models: dispersion analysis

This section is divided into two subparts. The first one is related to the analysis of the behaviour of infinite periodic 2D lattices by computing the dispersion curves of unit cells, for full cell-A and full cell-B configurations, depicted in Fig. 4. The second is dedicated to other infinite periodic 2D lattices by computing the dispersion curve of alternated (mixed periodic) A-B depicted in Fig. 6, that consists of periodic lattice.

Eigenfrequency analysis is used to compute the band diagram of periodic structures with first 25 Bloch modes for single cells and first 120 Bloch modes for the second case with combined (mixed periodic) A-B [START_REF] Brillouin | Wave propagation in periodic structures: electric filters and crystal lattices[END_REF][START_REF] Mace | Modelling wave propagation in two dimensional structures using finite element analysis[END_REF]. Floquet Bloch (FB) method is used; the periodic boundary conditions are written as

u Rx = e -jk x D u Lx , (2)
u Ry = e -jk y D u Ly ,

where u Rx and u Lx are the right and left side displacement along the x axis, and u Ry and u Ly are the displacement along the y axis, respectively. D is the length of reciprocal or unitary repeated element, k x and k y are the wavenumbers in the x and y directions. The first complete stop bands is observed in the ranges of 3000 -4500 Hz, and the second one is between 5000 -5800 Hz. There is a slight narrow band-gap from 5800 Hz to 6000 Hz which disappears rapidly. The third one starts from 7000 Hz to 8000 Hz which is quite wider and the fourth one with ∆ f = 2 kHz. There is also a very narrow one after 11000 Hz just before the next wider stop band with the ranges of 12000 -13800 Hz.

The overall results of the combined (mixed periodic) A-B show some kind of mix between the dispersion diagram of single cells shown in Fig. 8. Cell-A provides interesting dynamics properties with a few wide band gaps. Cell-B has more band gaps, but with narrow frequency width. The combination of the two in a periodic formation provides a large number of gaps, the largest ones being almost as wide as cell-A, while keeping a number of gaps comparable to cell-B. In the next section, FRFs will be computed on finite structures, in order to first, validate the applicability of the gaps for real structures, and second to check whether or not a quasi-periodic arrangement

Experimental measurements

In this section, the dynamical response of the lattice is validated with experimental tests. The quasi-periodic panel is manufactured by laser cutting as shown in Fig. 10.

The material properties of this test-article are already responded in Tab. 1.

White noise excitation point

Clamped boundary Observed response (output) Observed response (input) Shaker

Finite element analysis

In this section, two periodic and one quasi-periodic lattices are considered. The first two periodic lattices are modelled as finite meta-materials made of cell A and cell B in two separate lattices. The aim is to verify whether those stop bands predicted in the previous section are maintained for this specific finite element configuration i.e. with a finite domain. Secondly the quasi-periodic meta-material in Fig. 15 is compared with a combined cell A and cell B (periodic lattice), and in the last subsection an FRF of a meta-structure (filter junction) with a pure quasi-periodic lattice is described in details.

Point load on z direction

Fixed boundary conditions

Observed response point

Point load on transverse direction

Frequency response of meta-materials

The numerical models are clamped on the right edges and a harmonic point force is applied in the region close to the fixed boundary condition in Fig. 15 with a frequency ranging from 0 -35kHz. The model is meshed by 1061 free tetrahedral and 840 biquadratic triangles element 6 node. Frequency range of 35 kHz is selected for all the FRFs to check how the various cases behave in the higher frequencies.

The response is observed in dark blue spot location in Fig. 15. Fig. 16 shows an FRF for the three types of lattices. Two of them are fully periodic with cell A and cell B and the third one is quasi-periodic lattice. It can be seen from the plots that the predicted band gaps in the dispersion diagram of periodic lattices (Fig. 8) have similarities with the FRFs of the finite counter parts. There are basically three full band gaps, which are associated with the resonant and Bragg type band gaps. As there is no theoretical definition too generic to classify the types of band gaps in this dispersion diagram, however at least the first are resonant band gaps. Fig. 16 provides four band gaps marked as (BG1, BG2, BG3, and BG4). In BG1 the response are almost the same for all the cases, A, B and quasi-periodic arrangement (maybe even random arrangement of A and B) provides similar results.

BG2 has almost the same width for the 3 cases. It is the largest one but possesses a localised mode in the middle for the case A. The quasi-periodic is as large as case B but starts lower in frequency.

In BG3 the band gap is visible only in case A, very narrow on quasi-periodic but not visible as in case B. BG4 almost same width for case A and quasi-periodic, the quasiperiodic fall starts before than case A, no band gap for case B. In conclusion quasiperiodic can be seen as a good compromise between case A and case B.

I.

Operative deformed mode at frequency 2.54 kHz II.

Operative deformed mode at frequency 6 kHz III.

Operative deformed mode at frequency 7.3 kHz IV.

Operative deformed mode at frequency 19 kHz A comparison between the numerical models and the measurements is made in Fig. 22. Two band gaps can be seen in the plot for both numerical and experimental cases.

The first band gap between 2.5 and 5 kHz is almost similar in terms of width for both cases but with a low fall for the experimental part. The second band gap which starts from 8.1 kHz and continues to higher frequencies over 10.9 is similar too, for both cases. As it can be seen the second band gap in the experimental measurement has again low fall compared to the numerical one. The different depth of numerical and experimental band gaps can be function of many concurrent parameters. For example, for very response amplitude in the experimental measurements, the measured data are strongly influenced by the fundamental noise. On the contrary, the computed reduction of the structural response is simply a function of the numerical approximations. In the last section the structural dynamic behaviour of quasi-periodic lattice is computed by FE analysis and compared with the experimental results. The results show an acceptable agreement. In another case study, the FE model of quasi-periodic lattice is embedded between two bare panels which is used as a junction filter. The induced vibration energy transfers through the junction, and it acts as a meta-material filtering property to the elastic waves due to the impedance mismatch in the star-shaped unit cells. The results of FRF in terms of RMS of velocity amplitude gives some elements of novelty to reduce maximum unwanted vibration from the host structure.

The star-shaped unit cells sequences can thus be considered as a viable starting point to the optimisation of the final configurations for designing structures with desired frequency stop-bands.

Figure 1 :

 1 Figure 1: Meso-scale unit cell A and B preliminary design stage

Figure 2 :

 2 Figure 2: 2D mesh of cell-A and cell-B

s 4 (Figure 3 :

 43 Figure 3: Quasi-periodic lattice combining unit cells A and B: following Thue-Morse sequence in two directions

Fig. 5 Figure 4 :Figure 5 :

 545 Fig. 5 shows another periodic lattice with alternated unit cells A and B. As for the previous cases, the whole geometry includes 64 unit cells, in which 32 of them are cells A and 32 cells B.

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: FE model of the combined (mixed periodic) A-B with periodic conditions

Figure 9 :

 9 Figure 9: Band diagram of combined A-B lattice

  178

Figure 10 :

 10 Figure 10: Meta-material prototype under experimental test

Figure 11 :

 11 Figure 11: Manufactured meta-material with laser cutting

Fig. 12 Figure 12 :

 1212 Fig. 12 shows a measured FRF which is measured between 0 -100 Hz for the low frequency ranges. On this figure, the eigenfrequencies of the numerical model are shown using vertical lines. The first 3 modes are very well correlated in terms of frequency, which means that the low frequency behaviour of the structure is well captured by the FE model. In a second step, measurements are performed up to 10.9 kHz. The corresponding FRFs are shown in Fig. 13.

Figure 13 :

 13 Figure 13: Velocity response measured in the input point

Figure 14 :

 14 Figure 14: Velocity response measured in the output point

Figure 15 :

 15 Figure 15: Quasi-periodic lattice, numerical model setup for FRF analysis

Figure 16 :

 16 Figure 16: FRF of the two periodic and a quasi-periodic lattices

  . The two lattices are similar in terms of in-plane dimension and volume fraction.The responses provided in Fig.17in terms of velocity amplitude show similar dynamic effects with a slight difference. Only small difference can be observed as an higher depth around 15-dB for quasi-periodic one in blue line compared to strictly periodic in red line. The second difference is due to the wider frequency stop bands around 11 -19 kHz. The quasi-periodic curve has lower peaks starting from 16 kHz and continues up to higher frequencies around 25 kHz. In conclusion both designs have almost similar dynamic effects with a slight difference. The generated combination of quasi-periodicity in the lattice could give reduced response in the attenuation level of frequency range 3.5 -5 kHz compared with the regular periodic combinations. It also has a similar band gap wide in the medium frequency range 11 -19 kHz, with lower attenuation level and peaks after 25 kHz. lattice combined A,B Thue-Morse sequence Strictly periodic lattice combined A,B alternative sequence

Figure 17 :

 17 Figure 17: FRF of the periodic and quasi-periodic lattices

Figure 18 :Figure 19 :

 1819 Figure 18: Geometrical model of the meta-structure

Figure 20 :

 20 Figure 20: Operative modes of the meta-structure in four frequencies[START_REF] Billon | Mechanics and band gaps in hierarchical auxetic rectangular perforated composite metamaterials[END_REF].

Figure 21 :

 21 Figure 21: Iso-surface of three selected response frequencies

Figure 22 :

 22 Figure 22: Comparison of numerical and experimental FRFs

Table 1 :

 1 Material properties

	Material	Polymer
	Modulus of elasticity(Pa) 1.68 × 10 9
	Density (kg • m -3 )	818
	Poisson's ratio	0.38
	Structural loss factor	2.5 %

Acknowledgements

This investigation is carried out in the framework of the VIPER project (VIbroacoustic of PERiodic media). The project has received funding from the European Union's Horizon 2020 research and innovation program under Marie Curie grant agreement No 675441 and EUR EIPHI (ANR 17-EURE-0002) project and special thanks to the S.mart platform for the manufacturing process .