Nicolas Vanspranghe

Francesco Ferrante

Christophe Prieur

Control of a Wave Equation with a Dynamic Boundary Condition

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The aim of this paper is to study a wave equation in a bounded one-dimensional medium supplied with a dynamic (or kinetic) boundary condition at one end. The system is actuated via a Neuman-type control at the other end. Dynamic boundary conditions involve second-order time derivative and are typically obtained in physical models for which the momentum of the boundary cannot be neglected. A prime example is an infinite-dimensional model for the propagation of mechanical vibrations along drilling rods. In that case, the control is the torque applied to one extremity, and the kinetic boundary condition is given by the behavior of the rock-tip contact, which is subject to nonlinear friction. In particular, stick-slip phenomena can occur at the rock-tip interface and generate unwanted vibrations that might jeopardize the plant. A mechanical analysis of the rock-tip dynamics is given in [START_REF] Thiago | Non-linear dynamics of a drill-string with uncertain model of the bitrock interaction[END_REF]. The problem of stabilization and regulation of the velocity at the rock-tip contact has sparked interest in the control commmunity, see e.g. [START_REF] Saldivar | The control of drilling vibrations: A coupled PDE-ODE modeling approach[END_REF]. In engineering applications, the goal is to minimize the stick-slip vibrations through a suitable control. Various boundary control strategies are proposed to address this problem in [START_REF] Roman | Boundary control of a wave equation with indomain damping[END_REF], including backstepping design -see also [START_REF] Roman | Backstepping control of a wave PDE with unstable source terms and dynamic boundary[END_REF]. In [START_REF] Mlayeh | Backstepping boundary observer based-control for hyperbolic PDE in rotary drilling system[END_REF], an observer-based boundary control design is proposed. In [START_REF] Terrand-Jeanne | Regulation of Inhomogeneous Drilling Model With a P-I Controller[END_REF], stabilization and regulation using a proportional integral boundary controller is investigated. In [START_REF] Smyshlyaev | Boundary control of an antistable wave equation with anti-damping on the uncontrolled boundary[END_REF] and [START_REF] Bresch | Output-feedback adaptive control of a wave PDE with boundary anti-damping[END_REF], a similar problem is considered, but with a boundary anti-damping only involving first-order derivatives. However, the aforementioned papers deal with linearized models. We should also mention [START_REF] Hastir | Well-posedness of infinite-dimensional linear systems with nonlinear feedback[END_REF] and [START_REF] Ramirez | Stabilization of infinite dimensional port-Hamiltonian systems by nonlinear dynamic boundary control[END_REF] where different classes of boundary nonlinearities for distributed parameter systems are considered.

Let us introduce the specific dynamical model under study in this paper. Let L > 0 and Ω (0, L) ⊂ R. We deal with Nicolas Vanspranghe, Francesco Ferrante and Christophe Prieur are with Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, 38000 Grenoble, France. Email: name.surname@gipsa-lab.fr. This work has been partially supported by MIAI @ Grenoble Alpes (ANR-19-P3IA-0003). the following control system:

     ∂ tt u -∂ xx u = 0 on Ω × R + , ∂ tt u(0, t) -∂ x u(0, t) = F (∂ t u(0, t)) for all t, ∂ x u(L, t) = U (t) for all t, (1)
where F is a scalar function that models the nonlinear boundary friction, and U (t) is the control input.

Our first goal is to prove the well-posedness of the control system (1) supplied with the following collocated feedback:

U (t) = -g(∂ t u(L, t)), (2)
where g is a continuous nondecreasing scalar function.

The output considered here is the velocity at the boundary where the actuator lies, meaning that an observer is not required to implement the controller. Note that this class of feedback laws includes nonlinearities of interest in control applications, such as saturations or deadzones -see e.g. [START_REF] Prieur | Wave Equation With Cone-Bounded Control Laws[END_REF].

In this paper, we prove that, under appropriate assumptions, closing the loop in (1) with (2) leads to a wellposed dynamical system. A precise variational formulation of the problem is given and the regularity of the solutions is rigorously investigated. The underlying control problem is the stabilization of a possibly non-dissipative boundary by an action on the other boundary. We prove that, under suitable conditions, exponential stability can be achieved using a proportional controler. To the best of our knowledge, the stabilization of such system in the presence of nonlinear boundary anti-damping has not been investigated so far.

This paper is organized as follows. Section II introduces the variational formulation of control system (1) and contains the first main result, namely the well-posedness of the closedloop dynamics. Section III introduces the natural energy associated with [START_REF] Bresch | Output-feedback adaptive control of a wave PDE with boundary anti-damping[END_REF] and states an exponential stability result under appropriate assumptions on the nonlinear map F . This is the second main contribution. Section IV contains some illustrative numerical results. Section V gives the proof of the well-posedness result. Section VI contains concluding remarks.

Notation: Given a Banach space E, we denote its norm by • E and we use the duality bracket φ, x E to write φ(x) for any x ∈ E and φ ∈ E , where E is the topological dual of E. If E is also a Hilbert space, its scalar product is written (•, •) E . The space of infinitely differentiable functions on Ω with compact support is denoted by C ∞ c (Ω). Also, for T > 0, we denote by W1,p (0, T ; E) the subspace of L p (0, T ; E) composed of (classes of) E-valued functions f such that, for some g in L p (0, T ; E), f (t) = f (0) + t 0 g(s) ds for a.e. t in (0, T ). Such class f is identified with its continuous representative and we say that f = g in the sense of Evalued distributions.

II. VARIATIONAL FORMULATION AND WELL-POSEDNESS

In this section, we establish the framework in which we analyze system (1), and state our well-posedness result. We start by introducing the energy spaces associated with (1) as well as some notation. First, define

H L 2 (Ω) × R. (3)
We endow H with the usual product Hilbertian structure. Define now the following subspace of H:

V (u, u(0)) : u ∈ H 1 (Ω) , (4)
which is equipped with the scalar product inherited from

H 1 (Ω) × R.
As stated in Section V, V is also a Hilbert space. We also introduce the state space

X V × H (5)
endowed with the product Hilbertian structure. For the sake of clarity, we use parenthesis to denote elements of V or H, and brackets to denote elements of X , as in

X = [u, v] ∈ X , u = (u, u(0)) ∈ V , etc. Now, let us define the bilinear continuous symmetric form a on V × V by a(u 1 , u 2 ) Ω ∂ x u 1 (x)∂ x u 2 (x) dx. (6)
Finally, we denote by δ L the linear form mapping u ∈ V into u(L), which belongs to V since H 1 (Ω) is continuously embedded into C( Ω) in the one-dimensional case.

Assumption 1. The scalar function F is globally Lipschitz.

We define the nonlinear operator B on H associated with the first-order boundary pertubation:

∀v = (v, γ) ∈ H, B(v) (0, F (γ)). ( 7
)
This operator is globally Lipschitz by Assumption 1. Let us precise the meaning of weak solutions to (1).

Definition 1. A weak solution to (1) on (0, T ) is any u in L ∞ (0, T ; V ) ∩ W 1,∞ (0, T ; H) verifying T 0 -(u (t), φ (t)) H + a(u(t), φ(t)) dt = T 0 (B(u (t)), φ(t)) H + U (t) δ L , φ(t) V dt (8)
for all test-functions φ in C c (0, T ; V ) ∩ C 1 c (0, T ; H). Definition 1 is motivated by simple formal calculations in which one multiplies the wave equation by some smooth test-function φ(x, t) and integrates it over Ω × (0, T ) using integration by parts and the boundary conditions. Note that Definition 1 makes sense if, say, U belongs to L2 (0, T ). Closing the loop, a hidden regularity property of the solutions is needed to ensure all terms are defined. 1 Assumption 2. The function g : R → R is continuous and nondecreasing. 2 We show that the closed-loop system generates a dynamical system on X by defining the operators S t for t ≥ 0 as follows:

∀X = [u 0 , u 1 ] ∈ X , S t (X) [u(t), u (t)] ∈ X , (9)
where u is the unique solution associated with initial data

[u 0 , u 1 ].
Theorem 1 (Hadamard well-posedness). Let [u 0 , u 1 ] ∈ X . Under Assumptions 1 and 2, there exists a unique (weak) solution u ∈ C(R + , V )∩C 1 (R + , H) to (1) with feedback (2) and initial data [u 0 , u 1 ]. The solution u enjoys the following hidden regularity property: for all T > 0,

u(L, •) ∈ H 1 (0, T ) and g(∂ t u(L, •)) ∈ L 2 (0, T ). (10)
Moreover, we can associate with (1) and control law (2) the semigroup S = {S t } t≥0 of nonlinear continuous operators on X as defined in (9). Theorem 1 is proved in Section V. In the following proposition, we give some additional regularity when the initial datum is smooth and verifies a compatibility condition. We write W

[H 2 (Ω) × R] ∩ V , equipped with the scalar product inherited from H 2 (Ω) × R. Proposition 1 (Strong solutions). Let T > 0. If [u 0 , u 1 ] belongs to W × V and verifies ∂ x u 0 (L) = -g(u 1 (L)), then, u ∈ L ∞ (0, T ; W ), u ∈ L ∞ (0, T ; V ). (11)
Moreover, any weak solution u is the limit of a sequence of strong solutions [u n , u n ] in C([0, T ], X ), and

∂ t u n (L, •) → ∂ t u(L, •) in L 2 (0, T ).
Proposition 1 is a byproduct of the proof of Theorem 1 and is used to justify computations performed in Section III.

III. STABILITY ANALYSIS OF THE CLOSED-LOOP SYSTEM

Next, we analyze the stability of the closed-loop system when the feedback is linear. We introduce the energy functional E defined on X by

E(u, v) 1 2 v 2 H + a(u, u) , (12)
which is, in the context of abstract wave equations, the natural "mechanical" energy. Here, for all u in V and v = (v, γ) in H,

E(u, v) = 1 2 Ω |v(x)| 2 + |∂ x u(x)| 2 dx + 1 2 |γ| 2 . ( 13
)
For any ρ and µ in H 1 (Ω), we define a weighted energy functional (or Lyapunov function candidate):

Γ ρ,µ (u, v) 1 2 Ω v ∂ x u µ ρ ρ µ v ∂ x u + 1 2 µ(0)|γ| 2 . ( 14
)
It follows from the boundedness of the weights ρ and µ that there exists a constant M ρ,µ > 0 such that

∀[u, v] ∈ X , Γ ρ,µ (u, v) ≤ M ρ,µ E(u, v). (15)
A sufficient condition for the converse inequality to hold is the existence of C 1 , C 2 > 0 such that

µ(x) ≥ C 1 and µ(x) 2 -ρ(x) 2 ≥ C 2 (16)
holds for all x in Ω. We also remark that E and Γ are continuous on X . We note that the position does not appear in the energy functionals and without Poincaré-type inequalities, one cannot directly infer an estimate of the norm of the position. Stabilization shall be investigated with respect to the following invariant set, which is the line composed of constant solutions:

A {[θ1, 0] ∈ X : θ ∈ R}, (17)
where 1 ((x → 1), 1). Observe that A is exactly the subset of X where the mechanical energy E vanishes. For any

[u, v] in X , dist([u, v], A) 2 ≤ CE(u, v), (18)
where C is some positive constant. This is a consequence of the Poincaré-Wirtinger inequality -see [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF] for instance.

From now on, we are interested in stabilizing system (1) using a proportional feedback, i.e.

U (t) = -α∂ t u(L, t), ( 19
)
where α is a positive gain to be tuned.

Assumption 3. The function F is q-Lipschitz for some q < 1/2. Also, F (0) = 0.

The condition F (0) = 0 is quite natural when dealing with friction. Assumption 3 is meant to be in force in the context of nonlinear anti-damping. 3 Theorem 2 (Exponential stability). Under Assumption 3 and with α = 1, A uniformly and exponentially attracts the bounded sets of X . More precisely, there exist two positive constants M and η such that for all solutions u to (1),

E(u(t), u (t)) ≤ M E(u 0 , u 1 ) exp(-ηt). (20)
Moreover, for each solution u, u(t) converges in V to some constant function u ∞ when t → +∞.

Proof. We pick a solution u and, for the sake of concision, we denote by Γ u ρ,µ the function t ∈ R + → Γ ρ,µ (u(t), u (t)), which is continuous. Take ρ, µ ∈ H 1 (Ω). Assume temporarily that u is a strong solution. For all τ ≥ 0, denoting by 3 Alternatively, if F has a stabilizing effect, e.g. F (0) = 0, F nonincreasing, it can be relaxed to F locally Lipschitz -with appropriate modifications in the proof of well-posedness.

Q τ the rectangle Ω × (0, τ ), we have the following identity:

Γ u ρ,µ (t) τ 0 = - 1 2 Qτ ∂ t u ∂ x u ρ µ µ ρ ∂ t u ∂ x u + 1 2 ρ(L)(1 + α 2 ) -2αµ(L) τ 0 |∂ t u(L, t)| 2 dt - ρ(0) 2 τ 0 |∂ t u(0, t)| 2 + |∂ x u(0, t)| 2 dt +µ(0) τ 0 F (∂ t u(0, t))∂ t u(0, t) dt. (21)
Equation ( 21) is obtained by multiplying the wave equation

∂ tt u(x, t) -∂ xx u(x, t) = 0 by φ(x, y) = µ(x)∂ t u(x, t) + ρ(x)∂ x u(x, t) ∈ H 1 (Q τ )
, integrating over Q τ , and performing a few integrations by parts. If ρ and µ are nonnegative, writing |F (s)| ≤ q|s| where q is the Lipschitz constant in Assumption 3, we have the following inequality:

Γ u ρ,µ (t) τ 0 ≤ - 1 2 Qτ ∂ t u ∂ x u ρ µ µ ρ ∂ t u ∂ x u + 1 2 ρ(L)(1 + α 2 ) -2αµ(L) τ 0 |∂ t u(L, t)| 2 dt + 1 2 [2qµ(0) -ρ(0)] τ 0 |∂ t u(0, t)| 2 dt. ( 22
)
Take µ(x) = 1. For ( 16) to hold, it suffices to have ρ(x) ≤ 1for some > 0. Now, we derive some sufficient conditions for the energy to decay exponentially. It suffices to have ρ(0) ≥ 2q + , and ρ (x) ≥ a.e. as well as

ρ(L) ≤ 2α 1 + α 2 - (23)
for some > 0. Since q < 1/2, there exists an increasing affine function ρ such that ρ(0) > 2q and ρ(L) < 1. Let α = 1 so that (23) holds. As a result, by Grönwall's lemma, we obtain the following: with this particular choice of ρ and µ, there exists a positive constant (solution independent) η such that

∀t ≥ 0, Γ u ρ,µ (t) ≤ Γ u ρ,µ (0) exp(-ηt). (24)
By a density-continuity argument, the uniform estimate (24) holds for weak solutions as well -see Proposition 1. Since (16) holds, then there exists M > 0 such that (20) is verified by any solution. As for the second statement of Theorem 2, let u be a solution. Take an increasing sequence of nonnegative real numbers t n such that t n → +∞ when n → +∞. Then {u(t n )} n≥0 is a Cauchy sequence in V . Indeed, for any m ≥ n, writing the variation of u(t) between t n and t m , we have

u(t m ) -u(t n ) H ≤ tm tn u (s) H ds ≤ M +∞ tn exp(-ηs/2) ds, (25)
which converges to 0 when n → +∞; also, we already know that a(u(t), u(t)) → 0 when t → +∞, thus u(t m )u(t n ) V can be arbitrarily small. Using a similar argument, one verifies that the limit does not depend on the sequence {t n }. Therefore, u(t) converges to some u ∞ ∈ V and u ∞ is constant by (18).

IV. NUMERICAL SIMULATIONS We provide some numerical computations for illustrative purposes. We discretize (1) using finite elements over space and finite differences over time, on the basis of the functional formulation of the problem. Figures 1 and2 are obtained with the following parameters: we take α = 1, F (x) = qx with q = 0.1 and L = 1. The domain Ω is discretized into 100 points and the time step is set to 0.001. Further computations suggest that the condition on q derived in Theorem 2 is nearly sharp as taking q = 0.5 leads to an exponentially unstable system. Computations also suggest that the proportional feedback is not robust to (numerical) errors: taking values of q slightly below 0.5 induces unclear situations where an exponential decay is not easily identifiable.

V. PROOF OF THE HADAMARD THEOREM

The proof of Theorem 1 relies on nonlinear semigroup techniques and appropriate energy estimates. We begin this section with some remarks on the functional settings of the problem. The following result states some useful properties of the spaces H and V . Lemma 1. V is a separable Hilbert space isomorphic to H 1 (Ω). Moreover, V is a dense subset of H.

Since u 2

V = u 2 H + a(u, u), we see that V is continuously embedded into H. We denote H the topological dual of H. By Riesz representation theorem, we make the identification H H . Because V is a dense subset of H, the latter can be identified as a dense subset of V , leading to the classical injection chain V → H H → V , where each space is dense and continuously embedded into the following one. We denote by A the continuous operator from V into V defined by Au 1 , u 2 V a(u 1 , u 2 ) for all u 1 , u 2 in V . Next, define an unbounded (nonlinear) operator A g on X by we obtain, after several integrations by parts, the following boundary estimate:

τ 0 |g(∂ t u m (L, t)) -g(∂ t u n (L, t))| 2 + |∂ t w mn (L, t)| 2 dt = 2 L Ω x∂ t w mn (x, t)∂ x w mn (x, t) dx τ 0 + 1 L Qτ |∂ t w mn | 2 + |∂ x w mn | 2 , (43)
From (43), using Cauchy-Schwarz and Young inequalities, and the Lipschitz property of B, we obtain the following:

τ 0 |g(∂ t u m (L, t)) -g(∂ t u n (L, t))| 2 + |∂ t w mn (L, t)| 2 dt ≤ M sup t∈[0,τ ] w mn 2 V + w mn 2 H , (44)
for all τ ∈ [0, T ] and some M > 0. As a consequence, {∂ t u n (L, •)} and {g(∂ t u n (L, •))} are Cauchy sequences in L 2 (0, T ). Since u n (L, •) converges (in particular) in L 2 (0, T ), it follows that {u n (L, •)} is a Cauchy sequence in H 1 (0, T ), meaning that u(L, •) belongs in fact to H 1 (0, T ), which is the desired hidden regularity property. By continuity of g and unicity of the limit, g(∂ t u n (L, •)) converges to g(∂ t u(L, •)) in L 2 (0, T ) -one has to consider a subsequence converging for a.e. t ∈ (0, T ).

Step 3: Existence of weak solutions. Each u n verifies the distributional identity ( 8) for all test-functions φ in C c (0, T ; V ) ∩ C 1 c (0, T ; H). We pick φ and we let n → +∞.

As [u n , u n ] converges to [u, u ] in C([0, T ], X ), ∂ t u n (L, •) converges to ∂ t u(L, •) in L 2 (
0, T ), and g(∂ t u n (L, •)) converges to g(∂ t u(L, •)) in L 2 (0, T ), we obtain that u verifies (8) as well and is therefore a weak solution to [START_REF] Bresch | Output-feedback adaptive control of a wave PDE with boundary anti-damping[END_REF].

Step 4: Unicity and Hadamard continuity. Let w u 1 -u 2 , where u 1 and u 2 are two arbitrary weak solutions. The following energy identity holds: for any 0 ≤ τ . It is obtained by multipling the wave equation by ∂ t w(x, t) and integrating over Q τ . From there, using the fact that g is nondecreasing and B is Lipschitz, since we also have w(t) 2 H | τ 0 ≤ 2 τ 0 w H w H , we can apply the standard Grönwall argument to obtain unicity and continuity with respect to initial conditions, proving that the operators S t introduced in (9) are well-defined and continuous. At this point, Proposition 1 is also proved.

VI. CONCLUSION AND PERSPECTIVES

In this paper, the well-posedness of a wave equation with a nonlinear dynamic boundary condition and nonlinear Neuman-type feedback is proved, in a variational framework. Using a proportional feedback law, exponential stabilization is achieved under suitable assumptions.

In the linear case, an infinite-dimensional frequencydomain approach would be interesting to tackle the stabilization problem and obtain sharp conditions. Besides, adding an integral action to the feedback law should be considered. Also, the experimental device of [START_REF] Roman | Boundary control of a wave equation with indomain damping[END_REF] simulating drilling dynamics may be used to test the model in real experiments.

Fig. 1 .

 1 Fig. 1. Evolution of the boundary position u(0, t) over time. It obeys a second-order differential equation coupled with the wave equation.

Fig. 2 .

 2 Fig.2. Evolution of the mechanical energy E(u(t), u (t)) over time. The uncontrolled system would be unstable due to the boundary anti-damping.

  u 1 (t)) -B(u 2 (t)), w (t)) H -[g(∂ t u 1 (L, t)) -g(∂ t u 2 (L, t))]∂ t w(L, t) dt (45)

In the language of systems theory, one may say that the Neuman input operator is unbounded with respect to the state space X .

It is a specificity of the one-dimensional case that no additional growth assumption on g is required to obtain weak solutions for such feedback.

One should verify that the second coordinate is indeed the weak derivative of the first coordinate.

.

(26) Note that the domain D(A g ) need not be a subspace. We start with the following first-order abstract Cauchy problem: Ẋ(t) + A g (X(t)) = F(X(t))

where F is the nonlinear perturbation operator on X defined by F(X) [0, B(v)] for all X = [u, v] ∈ X . We see that F is Lipschitz, since the H-valued mapping B is Lipschitz. We wish to prove that (27) is a Lipschitz perturbation of an evolution equation with maximal monotone generator, hence the following result.

Proposition 2. The unbounded operator A g + id is maximal monotone.

Proof. We start with the monotonicity, and then we shall prove the surjectivity of [A g + id] + id.

Step 1: Monotonicity.

by u (resp. v), and also X 1 -X 2 by X. We have

. Moreover, by integration by parts, we have

. Thus, from (28) we obtain

From (29) we then obtain

which is the desired result.

Step 2: Surjectivity.

Replacing u with (f 1 + v)/2 in (32) and using the condition

and also

We wish to apply Lemma 3 given in Appendix to Θ.

As a Hilbert space, V is reflexive; it is also separable -see Lemma 1. For all

H ≥ 0, which proves the monotonicity. Items ( 2) and ( 3) are easily verified using continuity arguments. Finally, let w V ; we have L, w V ,V ≤ L V w V and Θ(w), w V ,V ≥ w 2 V , so that Θ(w), w V ,V -L, w V ,V → +∞ when w V → +∞, which implies the desired property. We deduce from Lemma 3 that there exists v ∈ V such that (35) holds. Then, let u (f 1 + v)/2 ∈ V . For all w ∈ V , we have

Since f 2 belongs to L 2 (Ω), taking φ = (φ, 0) in (36), where φ is an arbitrary element of C ∞ c (Ω), allows us to prove that u ∈ H 2 (Ω). Now, integrating by parts in (36) gives

holding for all w ∈ V . Pick a triangular function ρ such that ρ(0) = 1 and ρ(L) = 0. Then (ρ, 1) ∈ V . Evaluating (37) with ρ(x) and then ρ(L-x) yields -∂ x u(0)+2v(0) = f 2 (0) and ∂ x u(L) = -g(v(L)), hence X = (u, v) ∈ D(A g ) and A g (X) + 2X = Y, which concludes the proof.

In order to consider any initial datum in the state space X , we complete Proposition 2 with the following lemma.

We can finally prove the Hadamard theorem. Considering the evolution equation 27 with a sequence of approximate initial data in D(A g ) provides strong solutions that converge in C([0, T ], X ). Then, one has to prove that the limit is a weak solution. For that purpose, using a multiplier method, we retrieve a standard hidden regularity property of the wave equation.

Proof of Theorem 1. We split the proof into four steps.

Step 1: Approximate solutions. Since D(A g ) is dense in X , we can pick a sequence of vectors X 0 n = [u 0 n , u 1 n ] ∈ D(A g ) that converges to [u 0 , u 1 ] in X . For each n ≥ 0, there exists a (unique) strong solution X n ∈ W 1,∞ (R + , X ) to the following reformulation of (27):

Indeed, Proposition 2 states that A g + id is maximal monotone, and F + id is still Lipschitz. Existence is given by [START_REF] Edwin | Monotone operators in Banach space and nonlinear partial differential equations[END_REF]Corollary 4.1]. For all n ≥ 0,

in the sense of strong differentiation, for a.e. t > 0, with X n taking values in D(A g ). We pick w in V and take the scalar product of (40) with w. With an integration by parts, and using the definition of D(A g ), we obtain

for a.e. t > 0, which means that, w being arbitrary,

∞ (0, T ; H) and u n ∈ W 1,2 (0, T ; V ), a standard timeregularization argument -omitted here -allows us to infer from 41 that the distributional identity [START_REF] Roman | Boundary control of a wave equation with indomain damping[END_REF] holds for all φ ∈ C([0, T ], V ) ∩ C 1 ([0, T ], H), meaning that u n is indeed a weak solution to [START_REF] Bresch | Output-feedback adaptive control of a wave PDE with boundary anti-damping[END_REF]

Thus, u n (0, •) belongs to H 2 (0, T ) with the desired weak derivatives. It can be inferred from the distributional identity [START_REF] Roman | Boundary control of a wave equation with indomain damping[END_REF] and the additional regularity given by (42) that u n (x, t) belongs to H 2 (Q T ) and verifies the wave equation on Q T .

Step 2: Uniform convergence. The convergence of the sequence {X n } in C([0, T ], X ) for any T > 0 follows from the standard monotonicity argument. We denote its limit by X = [u, u ]. Coming back to the sequence of approximate strong solutions, for all (m, n) in N 2 , we write w mn u mu n . We also write w 0 mn w mn (0) and w 1 mn w mn (0).

Multiplying the wave equation by x∂

APPENDIX

We use the following existence result, which is an application of the Galerkin method and the Brouwer fixedpoint theorem in finite dimension -see [11, Lemma 2.1 and Theorem 2.1] for the proof. Lemma 3. Let E be a separable reflexive Banach space and f ∈ E . Assume Θ : E → E verifies the following assumptions:

1) Θ is monotone i.e. Θ(x 1 ) -Θ(x 2 ), x 1 -x 2 E ≥ 0 for all x 1 , x 2 ∈ E; 2) Θ is bounded i.e. S ⊂ E bounded implies Θ(S)

bounded in E ; 3) For all x 1 , x 2 ∈ E, the scalar function t → Θ(x 1 + tx 2 ), x 2 E ,E is continuous. If for some ρ > 0, x E > ρ implies Θ(x), x E ,E > f, x E ,E , then there exists x ∈ E such that Θ(x) = f .


Jean-Loup Christen
email: jean-loup.christen@notk.org

Mohamed Ichchou

Abdelmalek Zine

Bernard Troclet

Wave Finite Element Formulation of the Acoustic Transmission Through ComplexInfinite Plates

Afinite element-based derivation of the transmission loss (TL) of anisotropic layered infinite plates is presented in this paper. The wave-finite element method (WFE) isused to represent the plate with afinite element model of asingle unit cell. The incident acoustic field is aknown plane wave,and the reflected and transmitted pressures are supposed to be plane waveswith unknown amplitudes and phases. The periodicity conditions on the unit cell allow tofind asimple matrix equation linking the amplitudes of the transmitted and reflected fields as afunction of the incident one. This approach is validated for several cases against classical analytical models for thin plates and sandwich constructions, where the results agree perfectly for areasonable mesh size. The method is then used to study the effect of stacking order in alaminated composite plate. The main interest of the method is the use of finite elements, which enables arelative easy modelling since most packages readily include different formulations, compared to analytical models, where different formulations have to be implemented for every kind of material.

Introduction

Sound transmission modelling is av ery important topic in the industry,y et it is often still difficult to calculate for complexg eometries or materials. Several methods have been proposed for the computation of transmission loss through infinite plate-likes tructures. Analytic models based on the constitutive equations for infinite thin plates [START_REF] Fahy | Sound and structural vibration[END_REF], sandwich constructions [START_REF] Kurtze | Neww all design for high transmission loss or high damping[END_REF][START_REF] Narayanan | Sound transmission through adamped sandwich panel[END_REF][START_REF] Kumar | Predicting the sound transmission loss of honeycomb panels using the wave propagation approach[END_REF], and even double plates with periodical stiffeners [START_REF] Lin | Sound transmission through periodically framed parallel plates[END_REF] have been proposed in the past. The Transfer Matrix Method (TMM) [START_REF] Brouard | A llard: Ag eneral method of modelling sound propagation in layered media[END_REF][START_REF] Llard | Modelling multilayered systems with porous materials using the transfer matrix method[END_REF] is ag eneral framework used to compute acoustic transmission features of infinite plane isotropic layerings comprising elastic solids, poroelastic materials and air gaps. It relies on an analytical formulation of the problem in each layer,thus making it necessary to derive newequations for every kind of material that can be encountered (isotropic or anisotropic, elastic or poroelastic).

Finite Element (FE) modelling has the practical advantage overanalytical models that most constitutive laws and couplings between materials are already implemented in commercial codes. Howeverm eshing al arge and simple structure such as ap late can lead to high computational times. The Wave-Finite Element (WFE)m ethod is ar e-sponse to this issue in the case of large periodic structures. It couples finite elements and the periodic structure theory (also called Bloch'stheory), that has been successfully used in the past for the analysis of free wave propagation in mono-and bi-dimensional waveguides [START_REF] Droz | Ahybrid wavemode formulation for the vibro-acoustic analysis of 2D periodic structures[END_REF]. The WFE wasalso used to compute vibroacoustic indicators such as group velocities and modal densities for use in Statistical Energy Analysis (SEA) [START_REF] Cotoni | Shorter: Astatistical energy analysis subsystem formulation using finite element and periodic structure theory[END_REF][START_REF] Chronopoulos | A unified approach for the broadband vibroacoustic response of composite shells[END_REF]. Extensions of this method towards the calculation of forced responses of waveguides have also been proposed [START_REF] Renno | Calculating the forced response of two-dimensional homogeneous media using the wave and finite element method[END_REF], without considering fluidstructure coupling. An application of 1D WFE formulation wasused by Serra et al. [START_REF] Serra | On the use of transfer approaches to predict the vibroacoustic response of poroelastic media[END_REF] for finite multilayers comprising poroelastic materials, butthe transfer approach necessitated meshing the whole surface of the plate. An assessment of wave methods can be found in [START_REF] Orrenius | Wave modelling in predictive vibro-acoustics: Applications to rail vehicles and aircraft[END_REF], comparing 3 methods, namely TMM, the WFE+SEA procedure of [START_REF] Cotoni | Shorter: Astatistical energy analysis subsystem formulation using finite element and periodic structure theory[END_REF] and aRayleigh-Ritz procedure.

Since the periodic structure theory deals primarily with infinite structures, the WFE can be used to compute the transmission loss of infinite plane structures excited by a superposition of plane wavesd irectly,w ithout having to resort to SEA. The derivation of this application is proposed in this paper,m aking use of the WFE framework for aplate with homogeneous properties in its plane. This enables to model an arbitrary layering with the finite element model of as ingle unit cell, which is just 1e lement wide in each of the plane'sdimensions, as in [START_REF] Manconi | Modelling wave propagation in two-dimensional structures using aw ave/finite element technique[END_REF] and [START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF]. This paper is structured as follows: the WFE model under forced plane wave is presented in Section 2, with the derivation of the formulation of the TL and absorption coefficient. The method is validated on several cases in Section 3: an isotropic thin plate, an orthotropic plate, asandwich construction with as oft viscoelastic core and as tiff sandwich plate with honeycomb core. Adiscussion on the element size closes this section. Finally,Section 4presents an application case in which the proposed method is used to study the effect of the order of the layers in acomposite plate.

Forced response to an incident plane wave

We consider amultilayered plate with one surface lying in the plane z = 0. Each layer is considered homogeneous, while this is of course not the case for the assembly.A s ah omogeneous structure is periodic for anyp eriod, the unit cell is ap arallelepiped. It is modelled with standard (8 nodes)b rick elements. The FE discretisation has only asingle element in each of the x and y directions, and as manyasneeded in the z direction.

Equilibrium formulation

Let D denote the dynamic stiffness matrix of the unit cell. This matrix is computed from the finite element stiffness and mass matrices K, M with the relationship

D = K -ω 2 M. (1)
Viscoelastic damping can be taken into account through a complexpart of the stiffness matrix K.

The equilibrium equation of the unit cell then writes

Du = f + e, (2)
where f represents the forces due to the interaction with the neighbouring cells, and e the forces imposed on the plate by the acoustic pressure fields on each side. As can be seen in Figure 1, the displacement vector u can be partitioned as

u =   u 1 u 2 u 3 u 4   . ( 3
)
We nowsearch for solutions as plane waveswith the same wave characteristics as the projection of the incident wave on the plate, hence where λ x = exp(ik x d x )and λ y = exp(ik y d y ). The reduced variable u 1 can be used with the transformation

u 2 = λ x u 1 , u 3 = λ y u 1 , u 4 = λ x λ y u 1 ,
u =   I λ x I λ y I λ x λ y I   u 1 = Tu 1 , ( 4
)
where I is the identity matrix.

The imposed pressure is condensed to forces exerted on the nodes, and can be partitioned in the same way, e = Te 1 .

(

The internal forces exerted on each vertical segment are derivedf rom the stress field applied to the normal to the twofaces related to that segment. The linearity of the stress tensor with respect to displacement field u and the previous periodicity relations lead to

f 1 + λ -1 x f 2 + λ -1 y f 3 + λ -1 x λ -1 y f 4 = 0, (6)
or Tf = 0, where T = I,λ -1 x I ,λ -1 y I ,λ -1 x λ -1 y I .Equation (2) can then be premultiplied by T and rewritten as

TDTu 1 = TTe 1 = 4e 1 . (7)
Let'sd enote A = 1 4 TDT.W efi nally obtain the reduced equation,

Au 1 = e 1 , ( 8
)
linking the forces and the displacements on one single segment instead of the initial four.

Acoustic coupling

We are concerned with an oblique incident plane wave impinging the plate on one side, with aknown amplitude p I . The interaction between this incident wave and the plate create ar eflected wave on the same side, and at ransmitted wave on the other one, whose respective amplitudes p R and p T we want to calculate. The sound field on the incident side is then the superposition of the incident and reflected waves,

P I = p(x, y, z)| z>0 = p I exp (-ik z,I z) + p R exp (+ik z,I z) , (9)
where the common factor exp(i(ωtk x xk y y)has been omitted for legibility.T he wavenumber components satisfy the relationship k 2 x + k 2 y + (k z,I ) 2 = k 2 = (ω/c I ) 2 and c I is the speed of sound in the fluid on the incident side. In the same way, the sound field on the transmission side contains only the transmitted wave propagating in the same direction as the incident wave, 2 ,where c T is the speed of sound in the fluid on the transmission side. Potential phase difference between the pressure fields are accounted for through the fact that the amplitudes p R and p T may be complex. The wavenumbers k x and k y are conserved across the plate, so only the k z component may vary with the nature of the fluid. We will consider in the following that the fluid is the same on both sides, thus ρ I = rho T = ρ 0 and c I = c T = c 0 .

P T = p(x, y, z)| z<0 = p T exp (-ik z,T z)) , ( 10
)
with the relationship k 2 x + k 2 y + (k z,T ) 2 = (ω/c T )
In this case, the incidence angle θ is the angle between the wave vector and the normal of the plate, while the azimuthal angle φ givesi ts orientation in the plane. We therefore have

k x = k sin θ cos φ, k y = k sin θ sin φ and k z = k cos θ.
The load imposed on the plate can be written from these twop ressure fields lumped on the nodes of the finite element model. As the pressure force is exerted along the normal to the structure, the only non-zero terms in the vector e 1 of forces imposed on segment 1will be those relative to the DOFs in z-direction on the incident side, e I ,and on the transmission side e T .For an elastic material, these quantities are scalar.Let u O and e O denote the displacement and force vectors on all other degrees of freedom in the segment. The force exerted on the nodes of segment 1isthen, as shown in Figure 1 e

1 =   e I e O e T   =   S(p I + p R ) 0 Sp T   , ( 11
)
where S is the free surface of the element, which is identical on both sides. All other vectors and matrices can be written following the same decomposition, allowing to rewrite [START_REF] Droz | Ahybrid wavemode formulation for the vibro-acoustic analysis of 2D periodic structures[END_REF] as

  a II A IO a IT A OI A OO A OT a TI A TO a TT     u I u O u T   =   S(p I + p R ) 0 Sp T   . ( 12
)
The second line of ( 12) leads to

u O = -A -1 OO (A OI u I + A OT u T ) , ( 13
)
which allows to condense this equation into

b II b IT b TI b TT u I u T = S(p I +p R ) Sp T . ( 14
)
The fluid-structure interaction is characterized by the continuity of normal particle velocity at the interface. This writes for the incident side,

ρ 0 ω 2 u I = ∂P I ∂z = -ik z (p I -p R ) , ( 15
)
and for the transmission side:

ρ 0 ω 2 u T = ∂P T ∂z = -ik z p T . ( 16
)
Introducing acoustic admittance Y 0 = cos(θ)/(iωρ 0 c 0 ), we obtain

u I = Y 0 (p I -p R ) , u T = Y 0 p T .
These expressions for u I and u T can be re-injected into [START_REF] Manconi | Modelling wave propagation in two-dimensional structures using aw ave/finite element technique[END_REF],leading to twoscalar equations linking the unknowns p R and p T and the incident pressure p I ,

  b II + S Y 0 -b IT b TI -b TT - S Y 0   p R p T = p I b II - S Y 0 b TI . (17)
Solving this equation givesthe acoustic transparency τ = |p T /p I | 2 and the absorption coefficient

α = 1 -|p R /p I | 2 .
The transmission loss for aplane wave is then finally TL = -10 log 10 τ. [START_REF] Christen | Global sensitivity analysis of analytical vibroacoustic transmission models[END_REF] This leads also to the diffuse field transmission loss by integrating overa ll possible incidences between minimum and maximum bounds θ ∈ [θ min ,θ max ]a nd directions φ ∈ [0, 2π]. It is recommended to avoid the full range [0, π 2 ]for θ because numerical errors occur close to grazing incidences θ = π/2. However, θ min = 0c auses no issue. The diffuse field transparencyisthen

τ d = 2π 0 θ max 0 τ(ωθ,φ)sin θ cos θdθdφ 2π 0 θ max 0 sin θ cos θdθdφ , ( 19
)
and the diffuse field transmission loss

TL d (ω) = -10 log 10 τ d (ω). ( 20
)

Validation results

In order to assess the validity of the method, acomparison with various analytical models is performed in this section. Adiscussion on the size of the elements to be used is proposed in Section 3.5.

Isotropic plate

The first case is an isotropic aluminium plate, under two different excitations: a45 o plane wave and adiffuse field. The material parameters are the Yo ung modulus E = 70 GPa, the Poisson ratio ν = 0.33, the density ρ = 2700 kg.m -3 .T he loss factor is η = 0.5%, which leads to ac omplexY oung'sm odulus Ẽ = E(1 + iη). The unit cell consists of 10 elements through the thickness, with an overall thickness of 2mmand lengths d x = d y = 0.1mm.

Twok inds of analytical models are used for comparison. The first one is based on the Kirchhoff-Love plate assumption, while the other is the Transfer Matrix Method (TMM)d escribed in [START_REF] Llard | Modelling multilayered systems with porous materials using the transfer matrix method[END_REF]. All three model agree very well with each other,w ith ad i ff erence of less than 1dBo ver most of the frequencyrange, except in atinyregion around the coincidence frequency, which is slightly underestimated in the thin plate model. Because the TMM and the current model agree in this region, this difference must be due to the fact that the Kirchhoff-Love assumptions neglect shear,while the other twouse a3Delasticity theory.

The diffuse field integration has been performed over the range [0, 0.999π/2] to avoid numerical issues with grazing incidence. The range wasd iscretised with 1000 points. The three models are found in this case to agree up to 1dBo vert he whole frequencyr ange, the small difference at coincidence disappearing because of the averaging over1000 incidence angles.

Orthotropic plate

The diffuse field transmission loss of an orthotropic plate is compared against an analytical model [START_REF] Guyader | Transmission of reverberant sound through orthotropic viscoelastic multilayered plates[END_REF]. The material parameters are giveni nT able I, where E stands for Yo ung'sm odulus, ν for Poisson'sr atios and G for shear modulus. The TL results are presented in Figure 4for frequencies ranging between 100 Hz and 20 kHz. The mesh of the unit cell has 10 elements in the z direction and lengths d x = d y = 0.1mm. Here again, the analytical and current models agree well with each other,w ith an error Table I. Parameters of the orthotropic plate.

E x E y = E z ν xy = ν xz ν yz G xy = G xz (GPa)( GPa) (-)( -) (GPa) 224
6.9 0.25 0.3 56.58 smaller than 1dBoverthe whole range. There are twocritical frequencies due to the very different Yo ung'sm oduli in the x and y directions. An expression for these frequencies can classically be derivedfrom the analytical model, and is givenbythe formula

G yz ρη h (GPa)( kg.m -3 )( -) (mm
f crit,x = c 2 0 2π 12(1 -ν xy ν yx )ρ E x h 2 , ( 21
)
where ν yx = (E y /E x )ν xy .T he other critical frequency f crit,y is obtained by replacing E x by E y in (21).I nthis case, the frequencies obtained by these formulas agree well with the dips observed in the TL curve.

Stiff sandwich plate

The 6 th order constitutive equation for the bending of a sandwich beam wasgiven by Mead and Markus [START_REF] Mead | The forced vibration of at hreelayer damped sandwich beam with arbitrary boundary conditions[END_REF] and later extended to aplate by Narayanan and Shanbhag [START_REF] Narayanan | Sound transmission through adamped sandwich panel[END_REF].

The main assumptions are that the predominant deformation of the core is due to shear,while that of the skins is due to pure bending. This implies that the Yo ung modulus of the core is high compared to the out-of-plane shear modulus. The model has fivem ain parameters, namely skin bending stiffness D t ,overall bending stiffness B,damping η,surface mass m and shear parameter g.The constitutive equation reads

D t (1 + iη)∇ 6 w -g(D t + B)(1 + iη)∇ 4 w +m ∂ 2 w ∂t 2 - mB N ∂ 2 ∂t 2 ∇ 2 w = (∇ 2 -g)(p L -p R ), ( 22
)
where w = v/iω is the normal displacement of the plate.

In the considered frame where af orced wave is imposed on the plate with aw avenumber k = ω c sin θ,t he spatial derivative operator ∇ can be replaced by -ik.T his leads to the following expression of the impedance:

Z(ω, θ) = D t k 6 + g(D t +B)k 4 -mω 2 k 2 + g(1-ν 2 s ) iω(k 2 + g) . ( 23
)
Fors andwiches made of isotropic materials and identical skins, the skin bending stiffness is 2 /2 and the shear parameter is g = 2G c /(E s h s h c ), where G c = G c,xz = G c,yz is the core transverse shear modulus, E s and ν s = ν s,xy are the in-plane Yo ung'sm odulus and Poisson ratio of the skins, and h s (resp. h c )isthe thickness of one skin (resp. the core). In-plane core shear modulus and all Yo ung'smoduli of the core are neglected in the analytical model.

D t = E s h 3 s /[6(1 -ν 2 s )] the overall bending stiffness is B = E s h 2 c h s (1 + h s /h c )
The twom odels are compared between 1kHz and 20 kHz on as andwich plate with isotropic skins and orthotropic core impinged by a45 o plane wave.The parameters of the sandwich construction are giveni nT able II, and the TL results shown in Figure 6. With amesh using 3e lements for each skin and 45 for the core, and again d x = d y = 0.1mm, the twom odels differ also by less than 1dBi nt he whole range. Since the analytical model takes shear into account, the agreement is also excellent at coincidence.

Sandwich plate with soft core

The same calculation wasm ade with at hick sandwich plate made of ar ather soft isotropic core and thick skins. The material properties are described in Table III. The transmission loss curvei ss hown on Figure 6. Again this correlates well with the analytical TMM, with amaximum error of 1dBoverthe considered range. Three dips in TL in the core, while the third one around 5620 Hz is due to coincidence with the bending wave.The dispersion curves for these twowaves are presented in Figure 7.

The displacement field can be calculated back from [START_REF] Orrenius | Wave modelling in predictive vibro-acoustics: Applications to rail vehicles and aircraft[END_REF].T heya re displayed at several frequencies on Figure 8. The first one is related to the symmetric motion of the skins (Figure 8a), while the second is related to a global bending motion of the sandwich, or antisymmetric motion of the skins (Figure 8b). In both cases, it can be observed that ashearing motion of the skins is involved, due to their rather high thickness. In the antisymmetric case, the cross-section remain straight and vertical, indicating that only shear is at stake, with practically no bending motion.

The field represented on Figure 8c is at afrequencybetween the twocoincidences, showing adecoupling of the twoskins. On the latter example, it may be noted that the amplitude of the z-displacement is higher on the transmission side than on the incident side, yet the amplitude of the transmitted wave is of course smaller than that of the incident wave.This is due to the fact that the pressure field on the incident side is result of interference between incident and reflected waves, while the transmitted field consists of asingle plane wave.

Discussion on element size

In the validation cases presented in this section, the mesh parameters are chosen so as to ensure convergence with the TMM model, which is based on 3D isotropic elastic theory and therefore exact in this case. Ac oarser mesh leads to errors in the coincidence frequencyregion. Of all the wavesconsidered in acoustic transmission, the acoustic wave in the fluid has the shortest wavelength, at the highest frequency( 20 kHz)c onsidered in all these cases. Since it is imposed on the structure at grazing incidence, it is relevant to compare the mesh size to it. The element sizes in x and y directions must therefore satisfy about 100 elements per wavelength, which is much stricter than the usual 10 elements per wavelength used for finite elements. Since there is only one element in the in-plane directions, this has no impact on the number of degrees of freedom in the mesh. However, the results are also found to be quite sensitive to the aspect ratio of the elements, which means that the element size in the z-direction must be about the same as that in the plane. Forthick layers, asmall d x may thus lead to ahigh number of elements through the thickness. On the other hand, it also sets amaximum size for d x and d y ,since there must be at least one element per layer. Above 10 kHz, this latter limit is usually less restrictive than the "100 element per wavelength" criterion for plates with thickness in the order of 1cm or less.

Influence of stacking on the transmission loss

The current approach can be used to study the influence of the stacking sequence of amultilayer plate on the TL. We consider a1.6 mm thick composite plate made of 24 layers of the same transverse isotropic material described in Table IV in order to study the effect of the stacking sequence on the TL. Three different stacking sequences have been studied, which are described in Table V. All of these have the same proportion of layers with fibres oriented towards each direction, at -45, 0, 45 or 90 degrees of the global coordinate system, orientation 0being aligned with the x axis. The layers are supposed to be perfectly bonded to each other,sothe FE model ensures continuity of displacements at the interfaces. The three configuration are aclassical symmetric quasi-isotropic stacking ("initial" configuration in Table V),arandom permutation of it ("random" case), and agrouping of the same orientation in ascending order ("sorted" case).

The TL wasc omputed for these three sequences with the method described above,u sing 2e lements per layer and amesh size of 0.1 mm in the x and y directions. The results are presented in Figure 9.

Concerning the lowfrequencybehaviour,itcan be seen that the stacking sequence has practically no influence for frequencies below2kHz, and little up to 4kHz. This is due to the fact that belowthis frequency, the plate behaves according to the mass law, and the actual stiffness has no effect on TL.

In high frequencies (above 8kHz), the random and initial stacking sequences have the same asymptotic behaviour.T his is due to the dominant effect of damping on the high-frequencyT L, above the highest critical frequency. These effects have been studied in [START_REF] Christen | Global sensitivity analysis of analytical vibroacoustic transmission models[END_REF]. The sorted case exhibits the same asymptotic behaviour,which is only visible for much higher frequencies, above 20 kHz and thus not represented on Figure 9 Between these, all configurations present acoincidence region, characterized by alower TL. Because of the presence of different layer orientations, the coincidence zone of all plates is strictly contained between the twoextreme critical frequencies of ap late where all layers are identically oriented. These bounds are givenb y( 21), which yields for this material f crit,x = 3770 Hz and f crit,y = 19 kHz. The random and initial configurations exhibit only one dip, indicating that theybehave closely to an isotropic plate in this respect. This is due to agood repartition of the layer orientations throughout the thickness. On the other hand, the sorted stacking exhibits aw ide coincidence re- gion with twocritical frequencies, abehaviour that is similar to that of an orthotropic plate.

The coincidence phenomenon is illustrated in Figure 10, which shows the free propagating wavenumbers as afunction of direction φ at afrequencyof7kHz, which is in the coincidence region for all 3cases. These results were obtained with aWFE calculation using the same mesh. The acoustic trace wavenumber k a = 2πf/c 0 is represented for ag razing incidence, and does not depend on φ.T he random stacking k-space has ashape close to acircle, and lies mostly below k a .The initial stacking has amore distorted k-space, partly above and partly belowt he acoustic wavenumber,b ut the deviation to the circle remains small, hence the narrowc oincidence region. Finally,t he sorted case has av ery stretched k-space, close to an ellipse with major axis oriented around -π/6radians, which also crosses the k a circle. This behaviour is related to the strongly inhomogeneous nature of this stacking case: for aw avew ith direction φ = 0, the stiffest layer (oriented along the fibers)i so nt he outer part of the plate, while a wave with φ = π/4" sees" the stiffest layer close to the central plane.

Conclusion

An ew method for computing the transmission loss of infinite plates made of arbitrarily layered materials has been presented in this paper.T he originality of the approach is that it couples the periodic structure theory to finite elements as in the WFE method, while taking into account the coupling with an acoustic fluid analytically. The method is validated against classical analytic Transfer Matrix Method computation for isotropic, orthotropic and sandwich constructions. The main interest of the method is that it natively takes into account complexbehaviour of the structure, such as symmetric motion of the skins with respect to the neutral plane. Besides, the use of finite elements makes it easy to retrieve the displacement field inside the structure. The computation time is much higher than that of analytical methods, butremains low, typically up to about 10ms per frequencyfor aplane wave.Apossible use of this method is the validation of analytical model for multilayers.

An application case is proposed to study the effect of the stacking sequence of alayered composite plate on the diffuse field TL. It is shown that the stacking sequence has an effect on the width of the coincidence range, and little influence outside it, the lowfrequencies being governed by mass lawand the higher by damping effects. The most extreme case of asorted stacking shows awide coincidence zone limited by twocritical frequencies, while the random permutation exhibits ad ip in an arrower zone around the lowest critical frequency( 7kHz in this case). In applicat ions where noise control is important at frequencies corresponding to the coincidence range, stacking sequence should therefore be considered for acoustic criteria too.

Figure 1 .

 1 Figure 1. FE model of au nit cell and notation for the applied forces.

Figure 2 .

 2 Figure 2. Transmission loss of an infinite steel plate for an oblique plane wave with θ = π/4o :a nalytic model (plate), +:TMM, -: current method, d x = 0.1mm.

Figure 3 .

 3 Figure 3. Diffuse field TL of an infinite steel plate +:analytical model (thin plate), -: current method, d x = 0.1mm.

Figure 4 .

 4 Figure 4. Diffuse field TL of an orthotropic plate + analytic,current method. The difference is less than 1dBo nt he whole range (100 Hz, 20 kHz).

Figure 5 .Figure 6 .

 56 Figure 5. TL of asandwich plate impinged by a45 o plane wave. +:Analytical model, -current method.

Figure 7 .

 7 Figure 7. Dispersion curves of the out-of-plane wavesofthe soft sandwich structure and trace wavenumber of the incident acoustic plane wave at 45 o .

Figure 8 .

 8 Figure 8. Displacement fields at several frequencies. The color indicates the absolute value of the z-displacement, in meters, for an incident amplitude of 1Pa. The incident wave comes from the top. (a) Symmetric coincidence at 5472 Hz. (b) Antisymmetric coincidence at 5623 Hz. (c) Displacement field at 5548 Hz.

FrequencyFigure 9 . 6 -5 π / 6 π / 3 -2 π / 3 π / 2 -π / 2 2 π / 3 -π / 3 5 π / 6 -Figure 10 .

 966332233610 Figure 9. Diffuse field transmission loss for 3di ff erent stacking sequences.

Table II .

 II Parameters of the stiff sandwich plate.

		E x (GPa)	E z (GPa)	ρ (kg.m -3 )	ν
	Skin	70	70	2700	0.33
	Core	0.08	88	0.1
		G xy (MPa)	G xz (MPa)	G yz (MPa)	h (mm)
	Skin	26310	26310	26310	2
	Core	0.04	44	30

Table III .

 III Parameters of the soft sandwich plate.

		E (MPa)	νρ (kg.m -3 )h (mm)
	Skin	8300	0.15	629.9	6.35
	Core		8.3	01 63 8.1
		1000 10	2000	4000	10000 20000

Table IV .

 IV Material parameters of asingle layer.F or atransverse isotropic material, we have E y = E z , ν xz = ν xy and G xz = G xy .

	E x (GPa)( GPa) E z	ν xy (-)( -) ν yz	G xy (GPa)
	175.0	6.90	0.2542	0.4689	4.18
	ρη (GPa)( kg.m -3 )( %) G yz	h (µm)
	2.35	1520	0.1	67
	Table V. The three stacking configurations.
	Initial Sorted Random [0; 90; -45; 90; 45; 0 3 ;-45; 90; -45;45;90; [45 3 ; -45 3 ;0 3 ;90 3 ] s [-45 6 ;0 6 ;45 6 ;90 6 ] -45; 0;90; -45; 45 2 ;90; -45; 45; 0; 45]

Acknowledgements

The authors gratefully acknowledge Airbus Defence and Space for their financial support. This work wasp artly conducted under the SMART-NEST project funded by the European Commission.