Aurelien Grolet
email: aurelien.grolet@ensam.eu

On a model reduction method for computing forced response using non-linear normal modes

This paper presents a method for constructing reduced order models using non-linear normal modes (NNM) in the context of non-linear vibrations. Starting from a discretised version of the non-linear problem, the non-linear normal modes (NNM) of the structure are computed using the Harmonic Balance Method (HBM). A two parameters (amplitude and phase) parametrization of the NNM is introduces and they are then used for the forced response construction, assuming that the solution involves only a single (non-linear) resonant mode. The solution is then eventually corrected by linear terms which helps decreasing the error outside resonance, in particular around anti-resonances. The procedure results in two algebraic equations containing only two variables, one controlling the amplitude of vibration and the other controlling the phase, thus leading to a drastic reduction in the number of degrees of freedom. The procedure is illustrated on a simple, but representative, example. It is shown that a single mode approximation is sufficient for computing a good approximation of the forced response around a particular mode.

Introduction

This work deals with the vibration study of non-linear structures with concervative non-linearities, such as geometric nonlinearity. In order to compute the response of nonlinear systems, the structure is generally discretized using Finite Element Model (FEM) [START_REF] Zienkiewicz | Finite element method its basis and fundamental[END_REF], leading to a set of non-linear differential equations describing the time evolution of the degrees of freedom (dof). The resolution of such non-linear equations relies most of the time on iterative procedures based on the Newton-like algorithms. For systems with large number of dof, the resolution procedure can become very time consuming, and may even not be realized in a reasonable time for non-linear industrial structures. In order to reduce the computation time, it has been proposed to build reduced models which allow to represent the dynamics of the system with a reduced number of degrees of freedom.

One of the most used reduction method, is the Galerkin projection which consist in constructing a reduction basis, and then to express the movement of the full system as a linear combination of the reduction basis vectors. The motion equations are then projected on a projection basis, which leads to a set of reduced non-linear differential equations giving the evolution of the reduced dof. In general, the projection basis can be different from the reduction basis (Petrov-Galerkin Projection), but, in practice, most of the studies choses a projection basis equals to the reduction basis (Galerkin projection). The construction of the reduction basis can be realized either a priori or a posteriori. In the first case, the reduction basis is constructed using only the available information on the structural system (eg. modal synthesis [START_REF] Thomas | Calcul de vibrations nonlinéaires de micro/nano structures piézoélectriques stratifiées par modèles réduits avec correction quasi-statique[END_REF], modal derivative [START_REF] Slaats | Model reduction tools for nonlinear structural dynamics[END_REF][START_REF] Jongh | A reduction method for nonlinear dynamic systems[END_REF], Krylov subspace [START_REF] Bai | Krylov subspace techniques for reduced order modeling of large scale dynamical systems[END_REF][START_REF] Bai | Structure preserving model reduction using a Krylov subspace projection formulation[END_REF], ...). In the second case, the reduction basis is extracted from data collected form numerical simulations, or experiments (POD [START_REF] Chatterjee | An introduction to the pod[END_REF], compagnon/dual modes [START_REF] Chang | Reduced order modelling for the nonlinear geometric response of some curved structures[END_REF],...).

Even if the Galerkin projection method can gives good results, it is not always clear how to select the vectors to be retained in the reduction basis. Indeed, the non-linearities can induce a coupling between several modes of vibration, and thus all those "coupled" modes have to be selected in the reduction basis.

In order to directly take into account the modal coupling induced by the nonlinearities, some authors have proposed to reduce the system using Nonlinear Normal Modes (NNM). NNM can be view as an extension of the linear normal modes for nonlinear dynamical systems. Several definitions are available depending on the type of considered solution. ( [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Rosenberg | On nonlinear vibration of systems with many degrees of freedom[END_REF][START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF] ). Using the idea of the center manifold theorems, the Shaw and Pierre definition [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF] has been used, for example, to propose reduced order models of nonlinear dynamic systems (see e.g. [START_REF] Apiwattanalunggarn | Component mode synthesis using nonlinear normal modes[END_REF][START_REF] Pesheck | Modal reduction of a nonlinear rotating beam through nonlinear normal modes[END_REF]). Another closely related way of reducing nonlinear system is through the use of the normal form theory, either in the first [START_REF] Jezequel | Analysis of non-linear dynamical systems by the normal form theory[END_REF][START_REF] Blanc | On the numerical computation of nonlinear normal modes for reducedorder modelling of conservative vibratory systems[END_REF] or second [START_REF] Neild | Applying the method of normal forms to second-order nonlinear vibration problems[END_REF] order formalism. More recently, a method based on the use of complex normal mode for constructiong reduced order models of systems with friction has been presented in [START_REF] Krack | Reduced order modeling based on complex nonlinear modal analysis and its application to bladed disks with shroud contact[END_REF].

In this paper, we propose a new model reduction approach based on the use of NNM computed with HBM for reducing structural systems with conservative non-linearities. The originality of the method rely on considering the spatio-temporal evolution of the non-linear mode, whereas most of the reduction technique presented in the introduction only retains the the spatial distribution.

In a first section we will briefly recall the use of the harmonic balance method for computing solutions of non-linear vibration problems. Next, we will expose the definition of NNM used in this paper and define the parametrization. Finally, the model reduction procedure will be exposed and we will demonstrate its efficiency on numerical examples.

1 Non-linear vibration and Harmonic Balance Method (HBM)

In this paper, the HBM will be used to compute the nonlinear modes and reference solutions. In this section, one briefly recall the main steps of the HBM for computing periodic solutions of non-linear differential equations.

Let's consider a structural system with conservative non-linearities. After discretization (n dof), the motion equations are given by a set of n non-linear differential equations: Forced dynamic equations

M ü + C u + K u + F nl (u) = F (1)
where u(t) is a vector containing the time evolution of the degrees of freedom, M, C and K are respectively the mass, damping and stiffness matrices, F(t) is the excitation force vector which is assumed to be periodic with period T = ω 2π , and finally F nl (u) is the vector of non-linear forces (supposed to depends only on the position).

The HBM consists in searching for an approximated solution u in the form of a truncated Fourier series up to the H-th harmonic, as follows:

u(t) = a 0 + H k=1 a k cos(kωt) + b k sin(kωt) (2)
where a k and b k are the vectors of Fourier coefficients associated with harmonic k. Substituting Eq.(2)into Eq.( 1) and projecting over the truncated Fourier basis results is a set of non-linear algebraic equations which can be written as:

Z(ω) u + F nl ( u) -F = 0 (3) where u = [a T 0 , a T 1 , b T 1 , . . . , a T H , b T H ]
T is a vector of unknown of size n h = n(2H + 1), and Z is the matrix of dynamic stiffness defined block-diagonally as Z = blkdiag(K, (Z k ) 1≤k≤H ) where Z k is given by:

Z k = K -(kω) 2 M kωC -kωC K -(kω) 2 M (4)
F and F nl ( u) correspond respectively to the excitation force and non-linear force in the frequency domain.

For an imposed frequency ω, the solutions of the algebraic system in Eq.( 3) can be computed using Newton like algorithms. For different forcing frequencies, the solution can be computed by applying continuation procedures, see e.g. [START_REF] Groll | The harmonic balance method with arc-length continuation in rotor stator contact problems[END_REF][START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF].

Non-linear Normal Modes

Here we recall the definition of Non-linear Normal Modes (NNM) that will be used through this paper, and the numerical procedure used to computed them.

Definition

In this paper, we consider the computation of NNM of structural system with concervative nonlinearities such as Eq. [START_REF] Afaneh | Nonlinear response of a clamped-clamped beam with internal resonance under sinusoidal excitation[END_REF]. Several definition of NNM are available, here we will only use the Rosenberg's definition of NNM [START_REF] Rosenberg | On nonlinear vibration of systems with many degrees of freedom[END_REF], and one can refer to [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Peeters | Nonlinear normal modes, part ii: toward a practical computation using continuation technique[END_REF] for further detail about NNM definition and computation.

In Rosenberg's definition, the NNM are the "vibration in unisson" of the undamped, unforced system of equation. This definition correspond to the fact that for a motion on a NNM:

• (i) all dof goes to zero (resp. to their extremum) at the same time,

• (ii) all dof vibrate at the same frequency,

• (iii) the displacement of each dof can be parametrized by the displacement of a unique reference dof.

The Rosenberg definition of NNM is very simple and puts some restrictions on the kind of computed NNM. The first restriction of this definition is the fact that NNM are solutions of undamped systems. So that using this definition, one needs to drop all damping terms (frictions, structural damping, ...) before computing NNM. Another restriction of the Rosenberg definition is the fact that all dof should vibrate in unisson, implying that there can be no phase shift between dof [point (i)]. In other words, all dofs have to be perfectly tuned, and solutions such as travelling waves (solutions with a constant space shift between two contiguous dof) do not fall into this definition of NNM. Travelling wave will not be considered in this study. However, the procedure presented in this paper could be extended to travelling waves in a quite straight forward manner (provided one has computed the travelling NNM).

Computation

NNM will be sought as the solution of the undamped unforced equation:

M ü + K u + f nl (u) = 0 (5)
Here, the HBM will be used to compute the NNM of Eq.( 5). In this case the solution is again sought under the form of the truncated Fourier series in Eq.( 2), but the angular frequency ω is also an unknown now.

Applying the HBM, leads to an under-determinated set of n h = (2H + 1)n algebraic equations for n h + 1 unknowns. This indetermination is due to the absence of initial conditions and from invariance by time translation of the unforced, undamped equation. That is if u is a solution of Eq.( 5) then the time shifted signal u τ is also a solution of Eq.( 5) for all τ ∈ R.

For each dof n, the HBM approximation can be re -written as

u n (t) = a n 0 + H k=1 a n k cos(kωt) + b n k sin(kωt) = a n 0 + H k=1 c n k cos(kωt + φ n k )) (6)
where c n k and φ n k are the amplitude and phase depending in the dof n and the harmonic k.

A new equation, termed phase equation [START_REF] Peeters | Theoretical and experimental modal analysis of nonlinear vibrating structure using non linear normal modes[END_REF][START_REF] Arquier | Two method for the computation of nonlinear modes of vibrating systems at large amplitudes[END_REF], is thus added to the n h algebraic equations in order to set the phase of the solution. Most of the time, the phase equation is constructed such that the velocities of one of the dof is set to 0 at t = 0 : ẋi0 (t = 0) = 0. In the case of the HBM, this phase condition turns to the following algebraic equation:

H k=1 kb i0 k = 0 (7)
As we are using Rosenberg's definition of NNM, all dof must go to their maximum/zeros at the same time. A simple choice to ensure this property is to set each phase φ n k to zero for each dof and each harmonics: φ n k = 0, ∀k, n. This case coresponds to the the case of the so called monophase NMM wich was also studied by M.Peeters in [START_REF] Peeters | Theoretical and experimental modal analysis of nonlinear vibrating structure using non linear normal modes[END_REF].

The HBM approximation, for the mono-phase NNM, then reduces to the following:

u(t) = a 0 + H k=1 a k cos(kωt) (8)
In other word, only the cosine terms will be retained in the computation of the NNM, i.e we search for solutions that are even in time (u(-t) = u(t)) Using the approximation in Eq.( 8), and applying the HBM projection, one gets a set of n nnm = (H + 1)n algebraic equations for n nnm + 1 unknowns a k , k = 0, .., H and ω, which can be put under the following form:

Z(ω) u + F nl ( x) = 0 (9) with u = [a T 0 , a T 1 , . . . , a T H ] T , and Z = blkdiag(K, (K -ω 2 M), . . . , K - (Hω) 2 M),
where "blkdiag" is an operator that constructs a block diagonal matrix from its inputs.

Equation ( 9) is solved using continuation techniques such as the arc-length continuation. At the first iteration, the initial frequency for the Newton-Raphson procedure is taken from the frequency of the linearised system, and the form of the first harmonic a 1 is initialized with the corresponding linear mode shape (the NNM are tangent to the LNM at low amplitude). Then the arc-length continuation procedure is applied in order to compute the NNM for further amplitude.

For the forced response computation two parameters will be used to "tune" the non-linear mode with respect to the excitation: the amplitude of a control dof and the phase. Here we will simply choose the control parameter as the amplitude q ∈ R of the first harmonic of a reference dof m 0 : q = a m0

1 . The NNM will be represented as follows:

u N N M = q[ a 0 q + H k=1 a k q cos(kωt)] = q[Ψ 0 (q) + H k=1 Ψ k (q) cos(kωt)] (10)
Where, by construction, we have Ψ m0 1 (q) = 1. Note that other normation condition could be used, such as norming the first harmonic shape w.r.t the mass matrix (see e.g. [START_REF] Krack | Reduced order modeling based on complex nonlinear modal analysis and its application to bladed disks with shroud contact[END_REF]). The choice of the control dof m 0 is important, in particular its amplitude q should not go to zero when varying the continuation parameter ω.

The parametrisation w.r.t. the phase φ is simply re-introduced by considering the effect of a (global) time translation on the HBM solution:

u N N M (t + τ ) = a 0 + H k=1 a k cos(kωt + k ωτ φ ) (11)
The j-th NNM can therefore be represented as:

Ψ j (q, φ)(t) = Ψ 0 (q) + H k=1 Ψ k (q) cos(kωt + kφ) (12)
This representation has the advantage that the position on the NNM backbone is controlled by only two parameters q and φ. Also, note that not only the spacial shape but also the time evolution of the NNM is kept into account using this parametrisation, this will be a key point of the reduction procedure that follows.

3 Model order reduction using Non-linear Normal Modes

Single-mode reduction

We now turn to the resolution of the damped and forced equation ( 1) by means of NNM reduction. A single mode approximations is considered around the main resonance of the j-th NNM, so that an approximated solution u of Eq.( 1) is sought under the following form:

u(t) = Ψ j (q, φ)(t) (13)
where Ψ j (q, φ)(t) is the resonant NNM. Note that this approximation allows to retain the spatial and temporal feature of the NNM at the same time. The amplitude and phase of the resonant NNM are obtained using a variational method. First the matrix problem in Eq.( 1) is transformed into the following (space-time) weak formulation:

T 0 δu T M ü + C u + K u + F nl (u) -F dt = 0, ∀ δu (14)
Then the single NNM approximation and its variation are computed as follows:

u(t) = Ψ j (q, φ)(t), δ u = ∂Ψ j ∂q δq + ∂Ψ j ∂φ δφ (15)
and they are substituted into the weak formulation, resulting in two algebraic equations for the amplitude and phase (q, φ) of the resonant NNM:

δq T 0 ∂Ψ j ∂q T M Ψj + C Ψj + K Ψ j + F nl (Ψ j ) -F dt = 0 ( 16
)
δφ T 0 ∂Ψ j ∂φ T M Ψj + C Ψj + K Ψ j + F nl (Ψ j ) -F dt = 0 ( 17
)
The two equations in ( 16)-( 17) can simply be written as G(q, φ, ω) = 0. Their resolution for a fixed excitation frequency ω then gives the values for the variable q and φ. Continuation procedure, can be applied in order to extend the solution for further frequencies as described in the next section.

Numeric considerations

Let's consider the system in Eq.( 16), [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF] we have to evaluate the derivative of the NNM with relation to phase and amplitude. From the parametrisation of the NNM in Eq.( 12), the computation of the derivative with respect to q gives:

∂Ψ j ∂q = Ψ 0 + q ∂Ψ 0 ∂q + N h k=1 Ψ k + q ∂Ψ k ∂q cos(kωt + kφ) (18)
and the derivative with respect to the phase is given by the following:

∂Ψ j ∂φ = -q N h k=1 kΨ k sin(kωt + kφ) (19)
Plugging equations ( 12),( 18) and ( 19) into the variational equations ( 16)-( 17) and computing the integrals lead to the following form for the algebraic equation of the resonant NNM:

Ψ 0 + q ∂Ψ 0 ∂q T qKΨ 0 + F 0 nl -F 0 ex + k Ψ k + q ∂Ψ k ∂q T q[K -(kω) 2 M]Ψ k + F k,c nl cos(kφ) -F k,s nl sin(kφ) -( F k,c ex cos(kφ) -F k,s ex sin(kφ)) = 0 (20)
-q

N h k=1 kΨ T k -qkωCΨ k + F k,c nl sin(kφ) -F k,s nl cos(kφ) -( F k,c ex sin(kφ) -F k,s ex cos(kφ)) = 0 (21)
The two algebraic equations ( 20), [START_REF] Pesheck | Modal reduction of a nonlinear rotating beam through nonlinear normal modes[END_REF] resulting from the variational procedure are defined in a continuous manner over q and φ. However, the continuation procedure used to compute the NNM only gives outputs Ψ j (q, φ) at a discrete number of point q s . Therefore, for a given q a linear interpolation of the harmonic coefficients Ψ k (and their derivatives ∂Ψ k ∂q ) is computed between the two closest points of the (discrete) NNM results.

Forced response computation

For the forced response computation the two algebraic equations G(q, φ, ω) in Eqs.( 16),( 17) (or equivalently Eqs.( 20),( 21)) must be solve for different values of the excitation frequency. In order to account for possible turning points, continuation method are used in order to "follow" the response, for example using arc-length like parameters.

In the case of the proposed method, one can also adopt a different point of view and use the phase of the resonant NNM as the continuation parameter. This is motivated by the fact that the phase is always increasing during the non-linear resonance, and therefore does not undergo turning points (at least for the principal resonance). The phase then appears as a good parameter for the continuation, and a simple sequential continuation algorithm is sufficient to compute the whole resonance curve even in the presence of turning points (w.r.t. frequency) (see Fig. 4 and5).

In this case of simple sequential continuation the phase φ will be discretized from 0 to -2π and the amplitude and frequency are considered variable. In other word, we will search for q(φ), ω(φ) at discrete points φ = φ s . Note that here only a linear discretization is used for the phase, but grid adaptation could also be considered.

Linear corrections

In order to improve the solution given by the single resonant NNM approximation, a linear correction can be added as suggested in [START_REF] Krack | Reduced order modeling based on complex nonlinear modal analysis and its application to bladed disks with shroud contact[END_REF]. The corrected solution can be sought under the following form:

u(t) = Ψ j (q, φ)(t) + i =j Φ i η i (t) (22)
where the first term is the resonant NNM obtained by the procedure exposed in the previous sections and the second term is the participation of the linear modes. The second term, added a posteriori, is added mostly to take into account for eventual anti-resonances and also allows to decrease the error of the single NNM approximation away from resonance. In this study, the participations of the (non-resonant) linear modes will be assumed small so that it is reasonable to compute them using only the linear (linearised) system. The time evolution of the non resonant linear modes are then simply obtained from:

ηi + 2ξ i ω i ηi + ω 2 i η i -Φ T i F = 0 ( 23
)
where ξ i and ω i are respectively the modal damping and natural angular frequency of the i-th linear mode. Note that we have assumed that the amplitude and the phase of the resonant NNM could be computed independently of the amplitudes and phases of the linear modes, which seems reasonable as long as the modes a well separated. Consequently, the response of the linear modes can be computed in parallel (or before, or after) of the computation of the resonant NNM response.

Also, if the structure has a large number of degree of freedom, not all linear modes need to be retrained in the correction, and the sum in Eq.( 13) can be truncated as it is usually done in linear modal synthesis (the effect of the discarded modes could also be taken into account through a usual static correction).

Sum up of the reduction procedure

Finally, the reduction procedure in summarized in the following bullet points

• Compute the j-th NNM with HBM and parametrize it with relation to phase and amplitude to obtain Ψ j (q, φ) as in Eq.( 12).

• Use the following approximation for the forced solution

u(t) = Ψ j (q, φ)(t) + i =j Φ i η i (t)
• Non-linear part of the resolution -For a given excitation frequency ω find the parameters q and φ of the resonant NNM by solving the two algebraic equations in Eq.( 20)- [START_REF] Pesheck | Modal reduction of a nonlinear rotating beam through nonlinear normal modes[END_REF]. Generate the frequency response using arc-length continuation.

-Or, for a given φ find the parameters q and ω by solving the algebraic system in Eq.( 20)- [START_REF] Pesheck | Modal reduction of a nonlinear rotating beam through nonlinear normal modes[END_REF]. Generate the FRF using sequential continuation on φ ∈]0, -2π[ (no turning points with this parametrization)

• Linear part of the resolution (optional) -Compute linear corrections η i by solving Eq.( 23)

-Optionally, compute a static correction (if not all modes are included in the previous linear correction) [START_REF] Bai | Krylov subspace techniques for reduced order modeling of large scale dynamical systems[END_REF]

Numerical examples

In order to demonstrate the validity and efficiency of the previously exposed method, we propose an application to a simple (but representative) example: a linear beam with non-linear attachment.

Description of the system

In this section we consider the case of a linear beam clamped at one end and attached to a non-linear spring at the other end. The beam is made out of steel (E = 210GPa, ρ = 7800 kg.m -3 ), has a length L = 0.6 m and a square section of width a = 14 mm. The beam is modelled by standard Euler Bernoulli elements with 3 dof per node (axial disp. u, transverse disp. v, section rotation θ). A total number of 20 elements is used, leading to a problem containing n = 60 degrees of freedom (the first node is clamped). The effect of the attachment is modelled by a non-linear spring with the following expression for the restoring force:

F nl (v L ) = k 1 v L + k 3 v 3 L + k 5 v 5 L ( 24
)
where v L = v(L) represents the transverse tip displacement, and k 1,3,5 are some stiffness coefficients. Two sets of coefficient will be studied. First only a cubic restoring force will be considered (k 5 = 0). Then a quintic term will be introduced in the restoring force in order to generated more complicated behaviour for validating the proposed reduction method. The solutions computed with the HBM, in the free and/or forced cases, will be considered as reference solutions for comparison purpose.

Cubic restoring force

In this first part of the example, only a cubic restoring force will be considered, corresponding to the following set of coefficients for the formula in Eq.( 24): k 1 = 7.3e3, k 3 = 1e8 and k 5 = 0.

Non-linear mode computation

We compute the first three NNMs of the system by using the HBM with H = 5 harmonics. The backbone curve of the computed modes, in a Frequency-Energy Plot (FEP), are depicted on Fig. 1.

In this example, we will study the response of the beam in a frequency range around its first mode of vibration. Therefore, we will only analyse the first NNM (similar results could be obtained for the other modes).

It can be seen on Fig. 1. that the non-linear attachment has a stiffening effect on the first mode, and that at f = 70 Hz an 3:1 internal resonance occurs in the system between mode 1 and 2.

The information contained in the NNM can be analysed by projecting each harmonic of the NNM (i.e. each vector Ψ k (q) of Eq.( 12)) onto the (sub)space generated by (a subset of) the linear mode shapes. The results of such projection onto the first five Linear Mode (LM) shape is given in Fig. 2. Note that the abscissa is the point index, which can be seen as a curvilinear abscissa on the backbone curve.

One can see that, at first (for low point index), only LM-1 has a significant contribution onto the first harmonic (i.e. at low amplitude NNM-1 is equivalent to LM-1). When the point index increases, the contribution of LM-1 is decreasing and the contribution of mode 2 is increasing (with its highest contribution on the third harmonic): this confirms the presence of a 1:3 internal resonance between mode 1 and 2 (around point index = 125). Following the same lines, one can see that there also exist a 1:7 internal resonance between mode 1 and 3 (around point index = 155). An illustration of the evolution of NNM-1 is given on Fig. 3 where we plotted the evolution in the "phase space" and in the "configuration space" for four different positions on the NNM depicted by small circles on Fig. 2. In particular, one can see that when the internal resonance occurs, the configuration line cannot be parametrised only with one coordinates (the graph in the configuration space is not a bijection anymore).

This last remark is a direct consequence of the fact that the 2D invariant manifold intersects itself when experiencing internal resonances.

In order to avoid problems with internal resonance and parametrization, the maximum amplitude for the tip displacement in the numerical simulations will be limited to 0.024m.

As we restrict our study to a single non-linear mode approximation, NNM-1 is considered over a frequency range that does not include internal resonances (ie: [45, 70] Hz, corresponding to an amplitude range [0, 0.024]m) and its control parameter q is chosen to be the amplitude of the first harmonic of the

Forced response synthesis

We now use the previously computed NNM to synthesise an approximation for the forced response of the system. Several forcing amplitude will be considered in order to gradually increase the level of non-linearity in the responses. We recall the HBM solution will be considered as a reference solution for comparison purpose.

Figure 6 shows the FRF of the system (i.e.: the evolution of the transverse amplitude of the tip displacement as a function of the excitation frequency) for three forcing amplitudes. The NNM-synthetised solution are plotted along with the HBM solution (reference).

First one can plot the evolution of the NNM amplitude and phase q and φ as a function of the frequency (Fig. 4). Those curves exhibit the standard turning points associated with such nonlinear system. One can also take a different point of view and plot the amplitude q and the frequency as a function of the phase φ (Fig. 5). those curve don't exhibit turning point, making this representation better for continuation purpose, since a simple sequential continuation is sufficient. In the three forced cases the reduced model is able to give an excellent representation of the beam displacement (as compared to the reference solution) using only two variables q and φ, even for amplitudes close to the internal resonance.

A zoom of the resonance area is given on Fig. 7. We can observe that, close to the resonance, the NNM-reduced solution is closer to the backbone curve than the HBM solution. However, the difference remains small, and the maximum amplitude (at resonance) is well approximated by the reduction method. Finally, one can compare the time evolution between the reference solution and the NNM-reduced solution at the resonance for the highest force case. For comparison purpose we can also plot the time evolution obtained by a Galerkin projection keeping only the first linear mode, or the first four linear modes. This comparison is done in Fig. 8 where one can see that the time evolution of the solution obtained by two different methods (HBM and red NNM) are very close to each other. This also shows that keeping only the first linear mode is not sufficient at all to accurately represent the time evolution of the system at resonance. To obtain a level of accuraty close to the one obtained with the NNM reduction method, one would need to include at least the first four linear mode shapes in the reduction basis.

Quintic restoring force

In the remainder of the example a quintic term is added in the expression of the non-linear forces Eq.( 24), leading to the following set of coefficients: k 1 = 7.3e3, k 3 = 1e8 and k 5 = -7e10. This set of (somehow arbitrary) values has been chosen so that the system have a more complex behaviour compared to the previous case. In particular, the present choice of coefficients allows for the second NNM to exhibit an hardening behaviour followed by a softening one. In the remainder of the example, only this mode (i.e. the second NNM) will be considered. This will show that the proposed reduction method performs well in this non trivial case.

NNM computation

We computed the second NNM by using the HBM with H = 5 harmonics. The backbone curve of the computed mode is depicted in a Frequency-Energy Plot (FEP) on Fig. 9.

It can be seen that the non-linear attachment has, at first, a stiffening effect on the second mode (up to 215 Hz), followed by a softening effect. Also note that a modal interaction (internal resonance) occurs in the softening region (presence of a loop). In order to avoid (parametrisation) problems for the NNM, the amplitude of the tip displacement will be considered smaller than 0.05 (i.e. below the amplitude corresponding to internal resonance). The control parameter q is kept the same, i.e. the amplitude of the first harmonic of the tip displacement. We now use the previously computed NNM to synthesise an approximation for the forced response of the system. Again, several forcing amplitude will be considered in order to gradually increase the level of non-linearity.

Figure 10 shows the FRF of the system for four forcing amplitudes. The NNM-synthetised solutions are plotted along with the HBM solution (reference).

In the four cases there is an excellent agreement between the two methods. This shows that the forced response can be well approximated using only two algebraic variables q and φ, even in this non trivial case mixing hardening and softening behaviours.

Conclusion

In this paper we presented a reduction model procedure based on the use of nonlinear normal modes (NNM). First the NNM are computed using the harmonic balance method (HBM) and the Rosenberg definition of NNM. Therefore only a particular case of NNM is considered, namely mono-phase NNM.

The NNM are then used to compute single mode approximation of non-linear damped and forced equation. The solution of the equation is sought under the form of an NNM with an added phase in order to take into account the damping. This leads to only two algebraic unknowns which are the amplitude and the phase of the solution. The two equations controlling the amplitude and phase, are derived using a variational procedure. The procedure can be applied to localized or distributed non-linearity as long as they are conservatives.

The procedure is applied on a numerical example showing that this single mode approximation gives promising results. In the case of large industrial models, the computation of the nonlinear normal modes may be very costly. To reduce this computation cost, it might be judicious to first reduced the mechanical system using linear modes and a stiffness evaluation procedure [START_REF] Chang | Reduced order modelling for the nonlinear geometric response of some curved structures[END_REF] Further work should extends the procedure to travelling waves solutions and should also address the problem of multi (non-linear) modes approximation, in particular for the reduction of system with internal resonances.

annexe: Coupled version of the variational problem

The variation of the approximation is

u(t) = ∂Ψ j ∂q δq + ∂Ψ j ∂φ δφ + i =j Φ i δη i (t) (25)
Substituting into the variational formulation leads to a coupled version of the equations for the resonant NNM and the linear modes. The algebraic equations for the phase and amplitude are given by: T 0 ∂Ψ j ∂(q, φ)

T M Ψj + C Ψj + K Ψ j + F nl (Ψ j + Φ) + ∂Ψ j ∂(q, φ) T M Φ + C Φ + K Φ - ∂Ψ j ∂(q, φ) T F dt = 0 (26)
and the linear modes equations (n -1 differential equations, variable η i (t)) are given by:

Φ T i M Ψj + C Ψj + K Ψ j + F nl (Ψ j + Φ) +m i ηi + c i ηi + k i η i -Φ T i F = 0
In this paper, we made the hypothesis that we can neglect the coupling terms between linear and nonlinear modes. This leads to the set of formulas in Eq.( 16) and [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF].

Figure 1 :

 1 Figure 1: Backbone curve in Frequency-energy Plot for the first three non-linear normal modes of the cantilever beam with cubic non-linear attachment

Figure 2 :

 2 Figure 2: Projection of NNM-1 harmonic's onto the first five linear modes as a function of the point index (equivalent to a curvilinear absissa on the NNM curve)

Figure 3 :

 3 Figure 3: Evolution of NNM-1 along the backbone curve (phase space plot for the tip displacement v L /velocity dvL dt , and configuration space plot for the middle v L/2 and tip displacement v L ).

Figure 4 :

 4 Figure 4: Amplitude q and phase φ of the NNM as a function of the excitation frequency for different force level F = 20, 50, 120 N

Figure 5 :

 5 Figure 5: Amplitude q of the NNM and excitation frequency ω as a function of the phase for different force level F = 20, 50, 120 N. Note the absence of turning points as compared to Fig.4

1 Figure 6 :

 16 Figure 6: Comparison between the HBM solution (reference) and the NNMreduced solution for different force level F = 20, 50, 120 N for the case of a cubic restoring force

Figure 7 :Figure 8 :

 78 Figure 7: Detail of the previous figure around the resonance peak

Figure 9 :

 9 Figure 9: Backbone curve in Frequency-energy Plot for the second non-linear normal mode of the cantilever beam with quintic non-linear attachment

2 Figure 10 :

 210 Figure 10: Comparison between the HBM solution (reference) and the NNMreduced solution for different force level F = 50, 100, 150, 200 N for the case of a quintic restoring force


Naoki Haruta

Paulo Filho

Marques De Oliveira

Tohru Sato

Kazuyoshi Tanaka

Michel Baron

Force-Induced Dissolution of Imaginary Mode in Mechanochemical Reaction: Dibenzophenazine Synthesis

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

Mechanochemistry 1,2 is the study of chemical phenomena using mechanical actions such as grinding. In particular, mechanochemistry related to friction is well known as tribochemistry. [START_REF] Hsu | The Nature and Origin of Tribochemistry[END_REF][START_REF] Boldyreva | Mechanochemistry of Inorganic and Organic Systems: What Is Similar, What Is Different?[END_REF][START_REF] Tysoe | On Stress-Induced Tribochemical Reaction Rates[END_REF] Friction sometimes causes mechanical excitation on the order of eV, leading to luminescence called triboluminescence. 6 Such mechanical excitation can also drive chemical reactions. Mechanochemical reactions have attracted much attention recently because of their eco-friendliness as no solvent is needed. In conventional reactions, thermal excitation is a driving force that overcomes reaction barriers. However, mechanochemical reactions are not always conventional thermal reactions. Mechanochemical conditions sometimes cause unexpected chemical reactions. For example, Hickenboth et al. have reported that the Woodward-Hoffmann rule [START_REF] Woodward | The Conservation of Orbital Symmetry[END_REF] can be broken in mechanochemistry. [START_REF] Hickenboth | Biasing Reaction Pathways with Mechanical Force[END_REF] To explain these unique features of mechanochemistry, Ribas-Arino et al. have proposed that mechanical force modifies the potential energy surface and changes the reaction barriers. [START_REF] Ribas-Arino | Understanding Covalent Mechanochemistry[END_REF][START_REF] Ribas-Arino | Covalent Mechanochemistry: Theoretical Concepts and Computational Tools with Applications to Molecular Nanomechanics[END_REF] One typical example of mechanochemical reactions is the synthesis of dibenzophenazine (DBPZ), a well-known quinoxaline-containing molecule (see Figure 1). The quinoxaline structure is an important building block of organic semiconductors and pharmacologically active compounds. [START_REF] Martins | Solvent-Free Heterocyclic Synthesis[END_REF] DBPZ has been synthesized in various ways. Conventional syntheses require both a solvent and a catalyst. More specifically, Brønsted acids, [START_REF] Minsky | Novel Polycyclic Dianions: Metal Reduction of Nitrogen Heterocycles[END_REF][START_REF] Niknam | Metal Hydrogen Sulfates M(HSO 4 ) n : As Efficient Catalysts for the Synthesis of Quinoxalines in EtOH at Room Temperature[END_REF][START_REF] Niknam | Silica Bonded S-Sulfonic Acid: A Recyclable Catalyst for the Synthesis of Quinoxalines at Room Temperature[END_REF][START_REF] Rashedian | Silica-Bonded N-Propyl Sulfamic Acid: A Recyclable Catalyst for the Synthesis of 1,8-Dioxodecahydroacridines, 1,8-Dioxo-octahydroxanthenes and Quinoxalines[END_REF][START_REF] Karami | Synthesis of Aza-Polycyclic Compounds: Novel Phenazines and Quinoxalines Using Molybdate Sulfuric Acid (MSA). Polycyclic Aromat[END_REF][START_REF] Kumbhar | Brønsted Acid Hydrotrope Combined Catalyst for Environmentally Benign Synthesis of Quinoxalines and Pyrido[2,3-b]pyrazines in Aqueous Medium[END_REF][START_REF] Karami | A Modified Synthesis of Some Novel Polycyclic Aromatic Phenazines and Quinoxalines by Using the Tungstate Sulfuric Acid (TSA) as a Reusable Catalyst Under Solvent-free Conditions[END_REF][START_REF] Khodabakhshi | A Rapid and Eco-friendly Synthesis of Novel and Known Benzopyrazines Using Silica Tungstic Acid (STA) as a New and Recyclable Catalyst[END_REF][START_REF] Sajjadifar | Facile Method of Quinoxaline Synthesis Using Phenol as a New, Efficient and Cheap Catalyst at Room Temperature[END_REF][START_REF] Zhang | Synthesis and Photophysical Properties of Three Ladder-type Chromophores with Large and Rigid Conjugation Structures[END_REF] metal oxide nanoparticles, 1 6 , 1 8 , 1 9 , 2 2 , 2 5 , 2 7 , 2 8 , 3 0 , 3 6 , 4 2 , 4 3 and ionic liquids [START_REF] Zare | Ionic Liquid 1-Butyl-3-methylimidazolium Bromide ([bmim]Br): A Green and Neutral Reaction Media for the Efficient, Catalyst-free Synthesis of Quinoxaline Derivatives[END_REF][START_REF] Mulik | Proficient Synthesis of Quinoxaline and Phthalazinetrione Derivatives Using [C 8 dabco]Br Ionic Liquid as Catalyst in Aqueous Media[END_REF][START_REF] Vahdat | Sulfonated Organic Salts: Recyclable Green Catalysts for the Facile and Rapid Route Synthesis of 2,3-Dissubstituted Quinoxaline Derivatives in Water[END_REF][START_REF] Vahdat | A Green and Efficient Protocol for the Synthesis of Quinoxaline, Benzoxazole and Benzimidazole Derivatives Using Heteropolyanion-Based Ionic Liquids: As a Recyclable Solid Catalyst[END_REF][START_REF] Sajjadifar | A Practical Efficient and Rapid Synthesis of New Quinoxalines Catalyzed by Citric Acid as a Trifunctional Brønsted Acid at Room Temperature Under Green Condition[END_REF][START_REF] Sajjadifar | Synthesis of Quinoxaline Derivatives Using Sulfonic Acid Functionalized Imidazolium Salts as Highly Efficient and Reusable Brønsted Acidic Ionic Liquids Catalysts Under Solvent-free Conditions[END_REF] have been used. Recently, however, Carlier et al. have synthesized DBPZ using only mechanical ball-milling [START_REF] Carlier | Use of Cogrinding as a Solvent-free Solid State Method to Synthesize Dibenzophenazines[END_REF] without any solvent or catalyst and have obtained a very high yield (>99%).

To investigate the efficient mechanochemical synthesis of DBPZ, experimental and theoretical studies have been reported. [START_REF] Carlier | Greener Pharmacy Using Solvent-free Synthesis: Investigation of the Mechanism in the Case of Dibenzophenazine[END_REF][START_REF] De Oliveira | Comprehensive Experimental Investigation of Mechanically Induced 1,4-Diazines Synthesis in Solid State[END_REF][START_REF] Haruta | Reaction Mechanism in the Mechanochemical Synthesis of Dibenzophenazine: Application of Vibronic Coupling Density Analysis[END_REF] The yield obtained 24 h after mechanical grinding is greater than that obtained immediately after grinding. [START_REF] Carlier | Greener Pharmacy Using Solvent-free Synthesis: Investigation of the Mechanism in the Case of Dibenzophenazine[END_REF] That is, the product is generated even after grinding is complete. The yield obtained immediately after grinding is proportional to the grinding time, whereas the yield obtained 24 h after grinding is not. In addition, the reaction heat flow has two peaks, that is, early and delayed ones. [START_REF] De Oliveira | Comprehensive Experimental Investigation of Mechanically Induced 1,4-Diazines Synthesis in Solid State[END_REF] These experimental facts suggest the possibility of different reaction mechanisms.

In our previous work, [START_REF] Haruta | Reaction Mechanism in the Mechanochemical Synthesis of Dibenzophenazine: Application of Vibronic Coupling Density Analysis[END_REF] DBPZ synthesis was theoretically investigated using vibronic coupling density (VCD) analysis. [START_REF] Sato | Vibronic Coupling in Naphthalene Anion: Vibronic Coupling Density Analysis for Totally Symmetric Vibrational Modes[END_REF][START_REF] Sato | The Jahn-Teller Effect: Fundamentals and Implications for Physics and Chemistry[END_REF] One of the reactants, o-phenylenediamine (PDA), has a high-lying highest occupied molecular orbital (HOMO). Another reactant, 9,10-phenanthrenequinone (PQ), has a lowlying lowest unoccupied molecular orbital (LUMO). According to the frontier orbital theory, [START_REF] Fukui | A Molecular Orbital Theory of Reactivity in Aromatic Hydrocarbons[END_REF][START_REF] Fukui | Recognition of Stereochemical Paths by Orbital Interaction[END_REF] the regioselectivity should be indicated by these frontier orbitals. However, the HOMO of PDA is strongly delocalized over the molecule. Thus, the frontier orbital theory fails to predict its regioselectivity.

The frontier orbital theory takes into account only stabilization due to charge-transfer (orbital) interaction and does not consider any molecular deformation. In contrast, the VCD theory takes into account stabilization due to molecular deformation. The VCD can be regarded as a reactivity index in the framework of the conceptual density functional theory (DFT). [START_REF] Parr | Density Functional Approach to the Frontier-electron Theory of Chemical Reactivity[END_REF][START_REF] Geerlings | Conceptual Density Functional Theory[END_REF] VCD analysis has identified reactive regions in fullerenes, nanographenes, and so on, for which the frontier orbital theory is not effective. [START_REF] Haruta | Reaction Mechanism in the Mechanochemical Synthesis of Dibenzophenazine: Application of Vibronic Coupling Density Analysis[END_REF][START_REF] Sato | Vibronic Coupling in Naphthalene Anion: Vibronic Coupling Density Analysis for Totally Symmetric Vibrational Modes[END_REF][START_REF] Sato | The Jahn-Teller Effect: Fundamentals and Implications for Physics and Chemistry[END_REF][START_REF] Sato | C 60 Bearing Ethylene Moieties[END_REF][START_REF] Haruta | Chemical Reactivity in Nucleophilic Cycloaddition to C 70 : Vibronic Coupling Density and Vibronic Coupling Constants as Reactivity Indices[END_REF][START_REF] Haruta | Vibronic Couplings in Cycloadditions to Fullerenes[END_REF][START_REF] Sato | Vibronic Coupling Density and Related Concepts[END_REF][START_REF] Haruta | Regioselectivity in Multiple Cycloadditions to Fullerene C 60 : Vibronic Coupling Density Analysis[END_REF][START_REF] Haruta | Reactivity of Endohedral Metallofullerene La 2 @C 80 in Nucleophilic and Electrophilic Attacks: Vibronic Coupling Density Approach[END_REF][START_REF] Haruta | Reactivity Index for Diels-Alder Cycloadditions to Large Polycyclic Aromatic Hydrocarbons Using Vibronic Coupling Density[END_REF] We have applied VCD analysis to PDA and PQ and successfully identified their reactive sites.

PDA has two conformers with C 2 and C s symmetries depending on the conformation of its two pyramidal amino groups. C 2 PDA is the most stable one, whereas C s PDA is unstable with one imaginary frequency. According to VCD analysis, C 2 PDA favors stepwise formation of C-N bonds between the two reactants, whereas C s PDA favors concerted formation. Because the mechanical force must increase the amount of C s PDA during grinding, the concerted reaction is considered to be possible under mechanochemical conditions.

Although VCD analysis predicts both stepwise and concerted mechanisms, these reaction pathways and their energy barriers, which are indispensable for quantitative and even qualitative discussion, have not been determined to date. Our goal in this study is to identify all the transition states and their reaction barriers and then propose reaction pathways under mechanochemical conditions from a theoretical viewpoint. In addition, more importantly, a novel concept of mechanical modification of a reaction pathway is also proposed: force-induced dissolution of imaginary frequencies.

COMPUTATIONAL DETAILS

We performed geometry optimizations and vibrational analyses of all the reactants, intermediates, transition states, and products using DFT calculations. We calculated the intrinsic reaction coordinates to confirm the validity of the obtained transition states. We also calculated the internal energies, enthalpies, and Gibbs free energies of all the species at 298.15 K and 1.0000 atm to examine the thermochemistry. All the energies presented in this paper include zero-point vibrational energy corrections. These DFT calculations were performed at the B3LYP/6-311G(d,p) level of theory using the Gaussian 09 program. [START_REF] Frisch | Revision C.01[END_REF]

RESULTS AND DISCUSSION

3.1. Reaction Pathways. DBPZ synthesis consists of two C-N bond formations and two dehydration processes. Figure 2 shows the obtained energy diagram of the reaction pathways and their activation energies. The structures of the reactants and products, intermediates, and transition states are shown in Figures 345, respectively. The lowest-energy reaction path was found to be the stepwise path, which involves successive formation of two C-N bonds. In the stepwise path, both the first and second C-N bond formations have almost the same high reaction barriers: 161.8 and 141.4 kJ mol -1 , respectively. Furthermore, the subsequent dehydration process requires higher energy than these C-N bond formations. Thus, the first dehydration process seems to be the rate-determining step. Note, however, that a water molecule produced in the same or a neighboring system catalyzes both the first and second dehydration processes, as discussed later. The catalytic water greatly reduces the activation energies. In addition to the stepwise path, a concerted one was also identified; it involves simultaneous formation of two C-N bonds. In the concerted path, the formation of two C-N bonds has the highest energy barrier. Because the subsequent dehydration processes have much lower barriers, they are considered to occur rapidly. Although the concerted path has a higher barrier than the stepwise path, the difference is only 184.9 kJ mol -1 , which can be overcome by mechanical excitation (on the order of eV) 6 but not by thermal excitation.

Regarding the concerted path, we should note that PDA, TS1′, and I2′ have C s symmetry with the imaginary frequencies of a″ modes, which are symmetry-lowering modes. Thus, the concerted path is not generally accepted as a reaction path (see Figure 6) because a transition state should have only one imaginary frequency by definition. Later, we show that the instability of this path can be resolved by applying an external mechanical force.

Table 1 shows the calculated heat of reaction. This reaction is exothermic and produces a large amount of heat because of aromatization and hydrogen-bond formation between DBPZ and two water molecules (see Figure 3d). These hydrogen bonds stabilize the final product by 51.4 kJ mol -1 . The energy of hydrogen-bond formation was found to account for ∼40% of the overall heat of reaction. The rest can be attributed mainly to aromatization.

3.2. Catalytic Effect of Water. Hydrogen bonds can be broken under mechanochemical conditions because the bond energies are smaller than the mechanical energy. The free water molecules interact with reaction centers such as -OH and -NHgroups and catalyze the dehydration steps. An energy diagram including the catalytic effect of water is shown in Figure 7. The structures of the related species are illustrated in Figure 8. The catalytic water lowers the activation energies of the dehydration steps by 60-90 kJ mol -1 . Therefore, neither dehydration step is the rate-determining step. In DBPZ synthesis, the formation of the two C-N bonds is rather important for its kinetics.

On the basis of the obtained results, we discuss plausible reaction mechanisms under mechanochemical conditions. All the reaction barriers are too high to be overcome by thermal energy. This reaction is therefore difficult under conventional thermal condition without any catalyst. This is consistent with many experimental reports. However, under mechano- Although water catalysis requires breaking of a hydrogen bond, the bond energy is sufficiently small (∼25.7 kJ mol -1 ), as shown in Table 1. Values adjacent to transition states indicate activation energies in kJ mol -1 . chemical conditions, mechanical action can generate sufficient energy (on the order of eV) locally. 6 Therefore, both the stepwise and concerted paths are possible during grinding. After grinding, the energy difference between the two paths cannot be ignored. The stepwise path, with its lower barriers, must be favored. However, even these barriers are too high. In solid-state reactions, energy transfer can occur among neighboring reactants, intermediates, and products. Therefore, the supplied mechanical energy and large heat of reaction can be reused to conduct the reaction successively. This is in contrast to reactions in a solution, in which solvents separate reactant molecules from each other. We recently reported a detailed analysis of the experimentally observed reaction heat flow and kinetics that supports the findings of the present study. [START_REF] De Oliveira | Comprehensive Experimental Investigation of Mechanically Induced 1,4-Diazines Synthesis in Solid State[END_REF] 3.3. Removal of Instability by External Force. The obtained concerted path has the imaginary frequencies of the a″ modes, which are symmetry-lowering modes. They originate from the instability of C s PDA (see Figure 9a). Consequently, the concerted path is not generally accepted as a reaction path. We are interested in whether the instability is conserved under mechanical conditions.

During mechanical grinding, a random external force is applied to each reactant inside microcrystalline grains or at the interface of the grains. We should herein pay attention to whether an external force removing the instability of the concerted path is possible or not. If possible, such an external force can be applied at a certain moment, thereby opening the concerted path. At other moments, an external force might be useless for opening the concerted path, and hence the stepwise path is preferred.

Application of an external mechanical force is theoretically equivalent to deformation of the equilibrium structure. A molecular deformation corresponding to an external force reportedly modifies the potential energy surface and changes the reaction barriers. [START_REF] Ribas-Arino | Understanding Covalent Mechanochemistry[END_REF][START_REF] Ribas-Arino | Covalent Mechanochemistry: Theoretical Concepts and Computational Tools with Applications to Molecular Nanomechanics[END_REF] Mathematically, this can be understood as the Legendre transformation of the potential energy surface as position r → force f.

To simulate application of a mechanical force to PDA, we considered a deformation of C s PDA defined as a change in the planarization of the pyramidal amino groups to the opposite direction (see Figure 9b). This deformation is described in detail in the Supporting Information. This mode is orthogonal to the a″ mode with the imaginary frequency because it has a′ symmetry. As shown in Table 2, we found that the imaginary frequency of the a″ mode becomes real when C s PDA is deformed by this a′ mode. Note that even if the a′ deformation is applied, the force along the a″ mode is still zero because of its symmetry. That is, the deformed C s PDA is always fully optimized with respect to the a″ mode.

The strength of the applied force can be estimated using a potential gradient ∂E/∂Q d , where E denotes the total electronic energy with the nuclear-nuclear potentials, and Q d is a massweighted coordinate of the a′ deformation mode. As shown in Figure 10a, the strength of the force is proportional to the deformation. Figure 10b shows that increasing the force makes the curvature of the a″ mode positive. The positive curvature is equivalent to a real frequency. A conceptual illustration is shown in Figure 10c.

This result shows that if C s PDA is mechanically deformed by the a′ deformation mode, then the a″ instability in the concerted path is resolved. That is, the mechanical action leads to a new reaction path. In practical terms, this effect could be attributed to the interfacial stacking of microcrystalline grains of PDA and PQ under mechanical pressure, as illustrated in Figure 10d. This effect is not possible under conventional thermal conditions. The mechanical force necessary to remove the instability is small (∼10 -4 a.u.). Accordingly, it does not significantly change the obtained energy barriers of the concerted path.

3.4. Origin of Instability Removal. To elucidate the origin of the force-induced dissolution of the imaginary mode, a simple model potential was constructed, which depends only on two mass-weighted normal coordinates, Q α and Q β . In this study, the two coordinates can be regarded as the normal coordinate of the a″ imaginary mode and a′ deformation coordinate, respectively. The model potential E is given by a The optimized geometry of C s PDA is the reference geometry.

E E Q E Q E Q Q ( ) ( ) ( , ) = + + α α β β αβ α β (1)
The frequency of mode α is imaginary, and thus, E α is a double minimum potential

E Q Q Q ( ) 1 2 1 4 ( 0, 0) 2 2 4 2 ω δ ω δ = + ! < > α α α α α α α α (2)
where ω α denotes the frequency of mode α, and δ α is a quartic force constant. By contrast, the frequency of mode β is real, and thus, E β is a single minimum potential

E Q Q ( ) 1 2 ( 0) 2 2 2 ω ω = > β β β β β (3)
where ω β is the frequency of mode β. Furthermore, a cubic force constant γ was considered as the intermode coupling

E Q Q Q Q ( , ) 1 2 2 γ = αβ α β α β (4)
Application of an external mechanical force shifts the equilibrium point. Applying a deformation Q β = Δ is equivalent to applying the force given by

F E Q Q Q 0, 2 i k j j j j j j y { z z z z z z ω = ∂ ∂ = Δ β β β = =Δ α β (5)
The potential curve under the condition Q β = Δ is obtained as

E E Q E E Q Q Q ( ) ( ) ( , ) 1 2 ( ) 1 4 1 2 2 2 4 22 ω γ δ ω = + Δ + Δ = + Δ + ! + Δ α α β α β α α α α α β (6)
Consequently, the curvature of mode α depends on F β owing to the cubic force constant

E Q F Q Q 2 2 0, 2 2 2 i k j j j j j j y { z z z z z z ω γ ω γ ω ∂ ∂ = + Δ = + α α α β β = =Δ α β (7)
The DFT-based potential surface of PDA was numerically fitted by the model potential, and the fitting parameters were determined as follows: ω α = 1.14i × 10 -3 a.u., δ α = 3.47 × 10 -8 a.u., ω β = 7.74 × 10 -3 a.u., and γ = 6.73 × 10 -7 a.u. Figure 11 shows the model potentials obtained using these parameters. The values of the cubic force constant γ were 0, 1, and 2 times the DFT-based value, which correspond to weak, intermediate, and strong coupling, respectively. The obtained result clearly illustrates that the cubic force constant is essential for the force-induced dissolution of the imaginary mode.

CONCLUSIONS

The reaction pathways in mechanochemical synthesis of DBPZ were investigated using DFT. According to our computational results, the minimum-energy path involves stepwise C-N bond formations between two reactants. In addition, a concerted path with a higher energy barrier was also found. The latter is generally not accepted as a reaction path because the highest-energy point along the path is not a transition state but a higher-order saddle point. However, even a small mechanical force dissolves the imaginary frequencies so that the higher-order saddle point is changed into a transition state. The mechanical force identified here is one example of how to make the imaginary frequencies real. To find such forces in a mechanochemical reaction may become a central problem in theoretical mechanochemistry.

■ AUTHOR INFORMATION

Corresponding Author *E-mail: tsato@scl.kyoto-u.ac.jp. Phone: +81-75-711-7849. Fax: +81-75-711-7849.

ORCID

Naoki Haruta: 0000-0002-6635-1570 Tohru Sato: 0000-0002-4554-8411

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

This research was supported by the JSPS-CNRS bilateral program. The Embassy of France (Exploration Japon 2017) is also gratefully acknowledged. Parts of the numerical calculations were performed at the Supercomputer Laboratory of Kyoto University and the Research Center for Computational Science, Okazaki, Japan. We would like to thank Editage (www.editage.jp) for English language editing.

■ REFERENCES

Figure 1 .

 1 Figure 1. Synthetic scheme of dibenzophenazine (DBPZ) from ophenylenediamine (PDA) and 9,10-phenanthrenequinone (PQ). DBPZ has a quinoxaline structure.

Figure 2 .

 2 Figure 2. Energy diagram with two reaction pathways. Values in blue adjacent to transition states indicate activation energies in kJ mol -1 . Initial species have different energies because PDA has different conformers.

Figure 3 .

 3 Figure 3. Optimized structures of (a) C 2 PDA, (b) C s PDA, (c) PQ, (d) DBPZ coordinated by two water molecules, (e) DBPZ, and (f) a water molecule. The red dashed lines indicate hydrogen bonds.

Figure 4 .

 4 Figure 4. Optimized structures of (a) I1, (b) I2, (c) I1′, and (d) I3. I1′ has a structure similar to that of I2 and C s symmetry with one imaginary frequency. The red dashed line indicates a hydrogen bond.

Figure 5 .

 5 Figure 5. Optimized structures of (a) TS1, (b) TS2, (c) TS1′, (d) TS3, and (e) TS4. The red dashed line indicates a hydrogen bond. Note that TS1′ has C s symmetry with the imaginary frequencies of a″ modes as well as the a′ reaction mode.

Figure 6 .

 6 Figure 6. Schematic representation of stepwise and concerted pathways on potential energy surface. A concerted path with multiple imaginary frequencies is not generally accepted as a reaction path. However, an external force can dissolve imaginary frequencies.

a

  ΔE, ΔU, ΔH, and ΔG denote the change in electronic and nuclearnuclear repulsive energies, change in internal energy, change in enthalpy, and change in Gibbs free energy, respectively, at 298.15 K and 1.0000 atm.

Figure 7 .

 7 Figure 7. Energy diagram showing the catalytic effect of water. The energy of an additional water molecule was properly subtracted.Although water catalysis requires breaking of a hydrogen bond, the bond energy is sufficiently small (∼25.7 kJ mol -1 ), as shown in Table1. Values adjacent to transition states indicate activation energies in kJ mol -1 .

Figure 8 .

 8 Figure 8. Optimized structures of (a) I4, (b) TS3′, (c) I5, (d) I6, (e) TS4′, and (f) DBPZ with two water molecules produced by water catalysis. The red dashed lines indicate hydrogen bonds.

Figure 9 .

 9 Figure 9. (a) a″(1) normal mode of C s PDA with an imaginary frequency of 252i cm -1 . (b) Totally symmetric deformation mode with the planarization of the pyramidal amino groups in the opposite direction.

Figure 10 .

 10 Figure 10. (a) Relation between the totally symmetric mechanical force acting on C s PDA and its molecular deformation. (b) Curvature ∂ 2 E/∂Q 2 of the a″(1) mode as a function of the mechanical force. (c) Conceptual illustration of potential energy surface as a function of a normal coordinate of the a″(1) mode and the mechanical force. (d) Illustration of interfacial stacking of microcrystalline grains of PDA and PQ under mechanical pressure.

( 1 )

 1 Baron, M. Towards a Greener Pharmacy by More Eco Design. Waste Biomass Valorization 2012, 3, 395-407.

Figure 11 .

 11 Figure 11. Model potentials determined using the obtained fitting parameters. The values of the cubic force constant γ were (a) 0, (b) 1, and (c) 2 times the DFT-based value.

Table 1 .

 1 Obtained Heat of Reaction a

	factor	ΔE (kJ mol -1 )	ΔU (kJ mol -1 )	ΔH (kJ mol -1 )	ΔG (kJ mol -1 )
	two hydrogen bonds	-51.4	-49.1	-54.1	18.7
	aromatization	-70.3	-65.3	-62.8	-92.9
	overall	-121.7	-114.4	-116.9	-74.2

Table 2 .

 2 PDAs Deformed by the Totally Symmetric Deformation Mode and the Frequencies of the a″(1) Normal Mode a

	totally symmetric deformation (10 -2 Å)	0	1	2	3	4
	frequency of a″(1) mode (cm -1 )	250i	167i	50	148	166


Chaima Soussi

Mathieu Aucejo

Walid Larbi

Jean-François Deü

Sound transmission through double-glazed window: numerical and experimental analyses

Keywords: Double-glazing, Experimental modal analysis, Finite element

Double-glazing is widely used in constructions due of their sound insulation performance. According to the standards, the prediction of their acoustic performances is carried out in laboratory using two reverberant rooms. It's well known that, in the low frequency range below the Schroeder frequency, the presence of these rooms can affect the measurements due to their modal behavior. Thus, the prediction of their acoustic and vibratory behavior is very complicated and depends on several factors such as the used materials, the boundary conditions and the diffuseness of the acoustic field. In this work, experimental and numerical studies on the vibroacoustic behavior in low frequency range of double-glazing are presented. First, a finite element model for the proposed fluid-structure coupled problem is proposed. Experimental modal analysis of the structure is then detailed. The results of this analysis are used to validate a numerical model of double-glazing. Finally, various models for numerical sound reduction prediction of double glazing are presented and compared.

INTRODUCTION

Double glazing is often used in noise reduction since a relatively high transmission loss is achieved by introducing of an acoustic gap between the panels. Due to various effects of fluid -structure interaction, the prediction of the acoustic performance of the double wall remains very complicated in the low frequencies whatever the method used. We propose in this paper a reliable method based on experimental measurements and numerical analyzes able to predict the vibroacoustic behavior of double-glazing. In the literature, the article (1), ( 2) and (3) studied experimentally the transmission loss through single, double and triple glazing. In [START_REF] Xin | Analytical modeling of sound transmission through clamped triple-panel partition separated by enclosed air cavities[END_REF], the effect of boundary conditions on this acoustic indicator of double panels was investigated. A comparison of different type of glass m odels can be found in [START_REF] Miskinis | Comparison of sound insulation of windows with double glass units[END_REF]. For the numerical prediction approaches, several methods are available. The choice of the numerical method depends on the computational cost and the treated frequency band. The Finite Element Method (FEM) is suitable to treat problems in low frequencies. In [START_REF] Larbi | Vibroacoustic analysis of double-wall sandwich panels with viscoelastic core[END_REF], this method is applied to develop a model of a sandwich plate with viscoelastic core. The FEM is applied in [START_REF] Sgard | A numerical model for the low frequency diffuse field sound transmission loss of double-wall sound barriers with elastic porous linings[END_REF] by the authors for the different layers of the sound barrier coupled to a variational Boundary Element Method (BEM) to account for fluid loading.

The present paper is organized as follows: In the first part, a finite element formulation of the vibro-acoustic problem is presented. The second section is devoted to the modal experimental analysis of the double-glazing. The impact hammer test is used to evaluate the modal parameters (natural frequencies, mode shape and damping) which were used to recalibrate a numerical model of the studied system. In the last part, numerical evaluation of the acoustic performances of double -glazing is presented. The difficulties of the numerical models are related to the representati on of (i) the acoustic excitation and (ii) the effects of the emitting and receiving rooms. Thus, four configurations are proposed from which two are selected. We study the difference between the case where an acoustic room is modeled and where the Rayleigh integral is used. The effect of the damping and the acoustic impedance is also analyzed.

Finite element formulation for the vibro-acoustic problem

Let us consider the linear vibrations of an elastic structure filled with a homogeneous, inviscid and compressible fluid, neglecting gravity effects. This configuration is chosen to model the doubleglazing system. We establish in this section the local equations of the coupled problem and the corresponding finite element matrix equations.

Local equations of the problem

We consider an elastic structure occupying the domain Ω 𝑆 𝑖 (𝑖=1, 2) at the equilibrium. The structure is subjected to a prescribed displacement 𝐮 on a part Γ 𝑢 and to surface force density 𝐅 𝑑 on the complementary part of its external boundary. The interior fluid domain is denoted by Ω 𝐹 and the fluidstructure interface by Σ (see Figure 1). We note 𝝈, the stress tensor and 𝜌 𝑠 , the density of the structure. In addition, since the compressible fluid is assumed to be inviscid, instead of describing its motion by a fluid displacement vector field, which requires an appropriate discretization of the fluid irrotationality constraint, we use the pressure scalar field 𝑝. Let us note 𝑘 the wave number in the fluid (𝑘 = 𝜔 𝑐 𝐹 ⁄ ) where 𝑐 𝐹 is the speed of wave propagation in the fluid, 𝜌 𝐹 the mass density of the fluid, and 𝐧 𝑠 the normal of the structure. The local equations of the fluid-structure coupled problem described previously are:

Eq. ( 1) corresponds to the elasto-dynamic equation in absence of body force; Eqs. ( 2), (3) are the prescribed mechanical boundary conditions represent the surface forces applied on Γ 𝜎 and the displacement imposed on Γ 𝑢 ,respectively; Eq. ( 4) results from the action of pressure forces exerted by the fluid on the structure; Eq. ( 5) is the Helmholtz equation; and Eq. ( 6) is the kinematic interface fluid-structure condition on Σ.

Finite element discretization

In the finite element context, the previous vibro-acoustic problem is described in terms of structure displacement 𝑢 and acoustic pressure 𝑝. The discretization of the variational formulation associated to the previous equations ( 1) to ( 6) leads to the following matrix system (see [START_REF] Larbi | Vibroacoustic analysis of double-wall sandwich panels with viscoelastic core[END_REF] for details):

( [ 𝐊 𝑢 -𝐂 0 𝐊 p ] + 𝑖𝜔 [ 𝐃 𝑢 0 0 𝐃 𝑝 ] -𝜔 2 [ 𝐌 𝑢 0 𝐂 T 𝐌 𝑝 ] ) [ 𝐔 𝐏 ] = [ 𝐅 0 ] (7)
with 𝐊 𝑢 , 𝐌 𝑢 and 𝐃 𝑢 are the structural stiffness, mass and damping matrices, 𝐊 p , 𝐌 𝑝 and 𝐃 𝑝 are the associated acoustic matrices, 𝐂 is the fluid-structure coupling matrix and 𝐅 is the structural load.

Experimental analysis of the double-glazing

The modal analysis of the structure is essential for understanding its dynamic performance. In this section, we are interested in experimental modal characterization (resonant frequencies and damping coefficients) of double glazing. The experimental results will be used for the validation of a developed numerical model.

Description of the structure

The complete structure presented in Figure 2.a is a domestic wooden window with two symmetrical opening mounted on a frame. The model is 1.45 m wide by 1.5 m high. It weights approximately 76 kg. Each opening is composed of double-glazed surrounded by a rubber seal and incorporated in a sash which is composed of two stiles and two rails (Figure 2.b). The complete, the window is composed of 92 pieces of different sizes and different materials such are wood, glass, rubber, etc. The double-glazing, presented in Figure 3, combines two glass panes separated by a 18 mm spacer and an argon gap. It is 0.581 m wide by 1.309 m high. The polyamide spacer is filled at least half by the desiccant, which is a silicon material used to absorb the moisture within the cavity through the perforations of the spacer. The gluing of the spacer on the window is don e with a butyl type sealant or resistant silicone. All the parts are connected and sealed with a visible silicone located around the outside edge.

Experimental Modal Analysis

The Experimental Modal Analysis (EMA) is a tool to characterize the modal properties of a complex structure through the identification of its natural frequencies, mode shapes and modal damping ratios. Basically, the implementation of an EMA requires the measurement of the Frequency Response Function (FRF) of the considered structure at different points of the latter (Figure 4). In this work, the "roving hamme" test is used with an impact hammer and three reference accelerometers. This method consists in fixing accelerometers and then impact the structure at several points.

Impact testing considerations

The basic principle of the roving hammer test is to hit the structure with an impulse and then measure its response. This needs to be done at several points on the structure to obtain a fair representation of the mode shapes. The impact hammer is equipped with a force sensor and a nylon tip.

Accelerometer

Impact hammer Figure 4 -The device of modal analysis The suspended double-glazing (see Figure 5) is meshed into 153 nodes. Excitation is done in all nodes with an average of three impacts and measurements are performed up to 400 Hz, with 0.5 Hz frequency resolution.

Numerical Model

In this section, the double glazing experimentally tested is modeled numerically with the finite element method. The objective is to find a numerical representation of the structure able of reproducing its dynamic modal behavior. Thus, the model found will be used for acoustic studies.

The numerical model consists of two 2D shells representing the two glass panels separated by a three-dimensional cavity filled with argon (sound speed 317 m/s and density 1.6 kg/m 3 ). The structure is considered free in the space. The properties of the glass defined for the two planes are 𝐸 = 60 GPa; 𝜌 = 2450 (kg/m 3 ) and  =0.23. An equivalent material, whose properties are 𝐸 = 0.1 GPa ; 𝜌 = 1000 (kg/m 3 ) and  =0.49, is defined for the group «spacer + butyl + desiccant». Quadrilateral shell elements for the solid and hexahedral elements for the fluid are used for the finite element discretization. The mesh sizes for the structure and the acoustic domain are controlled by the wavelength 𝜆 (m) which depends on the frequency range of interest. For the acoustic domain, the acoustic wavelength is 𝜆 𝑎 = 𝑐 𝑓 ⁄ where 𝑐 (m/s) is the speed of sound. For the structure, the flexural wavelength is 𝜆 𝑓 = √2𝜋 𝑓 ⁄ (𝐷 𝑀 ⁄ ) 1/4 where 𝐷 (Pa.m 3 ) is the flexural rigidity and 𝑀 (kg/m)) the surface mass density. Since it is recommended to use 6 elements per wavelength, we obtain a model containing 35604 degrees of freedom.

Results and discussion

Firstly, in order to validate the experimental measurements, the principle of reciproc ity must be verified. It means that the frequency response function measured between two points must be the same regardless of which of them is input or output. Two cases are tested; the first one where the excitation point P1 and the measurement P2 are in the same glass pane. The second one where the excitation point is in the first glass pane and the measurement point is in the second one. The results presented in Figure 6 show that the reciprocity is insured since the curves are superposed following the permutation between the excitation and measurement points. 8) has been implemented. Figure 7 shows the comparison of the first ten natural frequencies of the structure obtained from the experimental and numerical analysis. For a such complex structure, results are considered satisfying since the maximum difference of natural frequency is 6.8 %, observed for Mode 8. Table 1 presents a comparison of modes shapes 1, 2 3 and 6 experimentally measured and found by finite element method. Once again, a good agreement between the two methods is observed. [START_REF] El-Kafafy | Polymax Plus estimator: better estimation of the modal parameters and their confidence bounds[END_REF]. Results of the first ten modal damping presented in the Figure 8 show that this parameter variates from 0.9 % at the 6 th mode to reaches the maximum at the 2 nd mode with a value of 11.3 %. These data will be used for vibroacoustic numerical simulations of the problem.

Measurements of sound transmission 4.1 Experimental set-up

The airborne sound insulation of building elements such as walls, doors, windows can be evaluated in laboratory according to standards norms [START_REF]Acoustics -laboratory measurement of sound insulation of building elements: Part 1: Application rules for specific products; Part 2: Measurement of airborne sound insulation; Part 4: Measurement procedures and requirements; Part 5: Requirements for test facilities and equipment[END_REF]. These standards allow the determination of the acoustic performances of the structures by the measurement of the sound reduction index.

Measurement takes into consideration two adjacent reverberation rooms, sepa rated by a common wall containing an opening in which the test element is mounted (see Figure 9). During the laboratory measurements, 𝑅 is defined as the difference between the average sound pressure levels in the source

Figure 9 -Experimental set-up of acoustic experimental test

Other configurations

Standards [START_REF]Acoustics -laboratory measurement of sound insulation of building elements: Part 1: Application rules for specific products; Part 2: Measurement of airborne sound insulation; Part 4: Measurement procedures and requirements; Part 5: Requirements for test facilities and equipment[END_REF] assume that acoustic excitation is a diffuse field for the entire frequency band of measurement [100 Hz -5000 Hz]. This condition is respected only beyond the "Schroeder frequency" of the room which marks the limit between the modal behavior and the diffuse field. This frequency is generally high and inversely proportional to the dimensions of the chamber.

At low frequencies, to treat the problem numerically by respecting the condition of the diffuse field, it is necessary to model a room whose dimensions are larger than the wavelength which is penalizing for the cost of computation. To overcome this problem, the plane wave superposition method is an effective solution. The technique consists in exciting the structure by a diffuse sound field (Figure 10.b). This solution therefore makes it possible to avoid modeling the emitting chamber. On the reception side, it can be assumed that the window radiates in a semi -infinite field without considering the receiving acoustic cavity (Figure 10.c). The combination of these two solutions, presented by the last configuration (Figure 10.d), makes it possible to reduce the size of the model and to avoid the modal effects of the rooms.

For numerical methods, the FEM method is used to model the window. For the modeling of the transmission and reception sides, the method used will depend on the chosen configuration. For example, for cases where both chambers are modeled, the FEM will be used since it is a bounded problem. For the cases of radiation in a free field corresponding to an unbounded problem, several solutions are presented (Rayleigh, BEM or also FEM combined with a perfectly matched layer (PML) for the far field.

Numerical analysis

In this section, we focus our attention on the numerical analysis of configurations (b) and (d) to calculate the sound transmission of the double-glazing model previously established. For the two cases considered, the excitation of the structure is a diffuse sound field which is a superposition of random plane waves in phase and direction making a uniform sound pressure at every point . The difference between the two analyses is the modelling of the receiving side: In the first configuration, the structure radiates in a 3.25×4×5 m 3 rigid rectangular cavity. The structure is clamped to a wall of the room. In the second case, the Rayleigh integral is used in the receiving side. This condition represents the acoustic radiation of a plane structure mounted in a rigid baffle in a half-space. Since direct frequency analyses are carried out, only materials damping is defined in the models.

Results and discussion

The carried numerical models allow to determinate the Transmission Loss of the structure (TL). This parameter is calculated in decibels (dB) as the ratio of the sound power incident 𝑤 𝑖 by the sound power transmitted 𝑤 𝑡 by the structure: TL = 10 log 10 (𝑤 𝑖 𝑤 𝑡 ⁄ ) (9) A comparison of the results in narrow band and 1/3 octave is established. The frequency range for narrow band figure goes from 1Hz to 600 Hz and from 31 Hz to 500 Hz for the 1/3 octave band presentation. A damping of 1% is used for all the material of the models. Firstly, a comparison of the TL calculated with the two configurations (b) and (d) is illustrated in the Figure 11. We can see in the narrow frequency band, a fluctuation of the TL at the low frequencies with the presence of the receiving room. This is explained by the modal behavior of the acoustic cavity. In fact, the coincidence of the excitation with the acoustic natural frequencies, for example Mode 2 at 43 Hz, Mode 7 at 68 Hz and Mode 28 at 132 Hz presented in Figure 11, causes an increase in the sound pressure level in the room, which leads to a decrease of the measured TL. This is not the case when the structure radiated in a semi-infinite domain. The sound reduction calculated in the two configurations becomes simil ar from 200 Hz, a value around the cutoff frequency of the acoustic cavity.

Figure 11 -Comparison of the TL between the two configurations in the narrow and 1/3 octave band The second aim of this section is to study the influence of the damping in the TL of the structure. For this, three analyses are made with the configuration (d). Results presented by the Figure 12 disclosed the fact that the damping smooths the curve of the TL and influences considerably the results above 100 Hz. At 500 Hz, the model without damping present a results 15 dB lower than the case with 1% damping. Below 100 Hz, the TL of all models present a similar behavior.

CONCLUSIONS AND PERSPECTIVES

In this work, the vibroacoustic behavior of the double-glazing glass is studied. First, the modal behavior of the structure is experimentally identified and validated with the finite element method. The established numerical model was used to carry out different acoustic analyses. The numerical examples show that the presence of the acoustic room has a significant influence on the TL in the very low frequencies due to its modal behavior below the Schneider frequency. In addition, the dynamic damping and the impedance don't have an influence below 100 Hz but improve the acoustic performances in the higher frequencies.

Within this framework, many future investigations and developments can be considered among which: the calculation of the TL with other configurations to study the influence of the acoustic cavities and their geometries; the influence of the structural and the acoustic boundary conditions; the recalibration of a model of one opening and then the hole window.

Figure 1 -

 1 Figure 1 -Double-glazing system

Figure 2 -

 2 Figure 2 -(a) Wooden window, (b) Components of one opening

Figure 3 -

 3 Figure 3 -Details of the double-glazing

Figure 5 -

 5 Figure 5 -Suspended double-glazing

Figure 6 -

 6 Figure 6 -Verification of the reciprocity For analyzing the experimental data, a home-made EMA toolbox based on the pLSCF (polyreference Least Squares Complex Frequency domain) (8) has been implemented. Figure7shows the comparison of the first ten natural frequencies of the structure obtained from the experimental and numerical analysis. For a such complex structure, results are considered satisfying since the maximum difference of natural frequency is 6.8 %, observed for Mode 8.

Figure 7 -

 7 Figure 7 -Comparison of the first ten experimental and numerical eigenfrequencies of the double-glazingTable1presents a comparison of modes shapes 1, 2 3 and 6 experimentally measured and found by finite element method. Once again, a good agreement between the two methods is observed. Table 1-Comparison of experimental and numerical mode shapes Mode 1 2 3 6

Figure 8 -

 8 Figure 8 -Experimental modal damping of the double-glazing

  the receiving room with a correction term considering the acoustic absorption 𝐴 (m)) and the area of the test specimen 𝑆 (m)): 𝑅 = 𝐿 𝑠 -𝐿 𝑟 + 10 log 10 (𝑆 𝐴 ⁄ )

Figure 10 -

 10 Figure 10 -Different configurations considered for problem resolution

Figure 12 -

 12 Figure 12 -Effect of the damping in the TL of the double-glazing (configuration d)

Table 1 -

 1 Comparison of experimental and numerical mode shapes

	Mode	1	2	3	6
	Experimental
		29.9 Hz	37.7 Hz	41.2 Hz	75.5 Hz
	Numerical
		29.8 Hz	36.6 Hz	38.6 Hz	73 Hz
	The modal damping is often the parameter extracted with the greatest uncertainty

Chaima.soussi@lecnam.net

ACKNOWLEDGEMENTS

The authors would like to express their thanks to CODIFAB (Comité professionnel de développement des industries françaises de l'ameublement et du bois) for its financial support.