J Albella
email: jorge.albella@usc.es

H Ben Dhia

S Imperiale

J Rodríguez

Mathematical and numerical study of transient wave scattering by obstacles with the Arlequin Method

Keywords: Wave propagation, Domain Decomposition, Stability Analysis

In this work we extend the Arlequin method to overlapping domain decomposition technique for transient wave equation scattering by obstacles. The main contribution of this work is to construct and analyze from the continuous level up to the fully discrete level some variants of the Arlequin method. The constructed discretizations allow to solve wave propagation problems while using non-conforming and overlapping meshes for the background propagating medium and the surrounding of the obstacle respectively. Hence we obtain a flexible and stable method in terms of the space discretization -an inf-sup condition is proven -while the stability of the time discretization is ensured by energy identities.

Introduction

In this work we are interested in the mathematical analysis of the Arlequin method for wave propagation simulation in the context of transient wave scattering by bounded obstacles.

Considering a discretization process by continuous finite elements, we assume that meshes of the neighborhoods of each obstacle or inclusion and meshes of the background medium are provided independently of one another. We aim at proposing a method that enables to compute the solution of the wave scattering without involving a well-adapted global mesh of the configuration but instead involving a coupling of the solutions computed on each provided meshes in a non-conformal way. The advantage of such approach is twofold:

• No global mesh of the configuration is required. This is especially useful if no embedded mesh generator is available. Moreover, only the quality of the meshes independently of one another, have to be checked and optimized.

• Local space refinement can be intrinsically done. In the non destructive testing configuration described above, it is optimal to use a coarse regular mesh for the background medium. This mesh has to be adapted to the smallest wavelength of the emitted waves, whereas the mesh representing the defect could be finer and welladapted to the geometry only.

Among the class of methods that achieve the two mentioned advantages we focus on methods that

• guaranty discrete energy preservation: this is specifically useful for the numerical stability of the time scheme.

• allow the use of efficient explicit time-stepping method: the idea is to use an explicit time discretization method for the computation of the solution in the background medium with a time step not modified by the coupling procedure (one can think to locally implicit time integration as in [START_REF] Doleana | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF][START_REF] Descombes | Locally implicit time integration strategies in a discontinuous Galerkin method for Maxwell's equations[END_REF] or to local time stepping method for the computation in the fine mesh as in [START_REF] Becache | Space-time mesh refinement for elastodynamics[END_REF][START_REF] Derveaux | Effective computational methods for wave propagation[END_REF]).

• is compatible with high order discretization in space: high order methods for the computation of transient waves propagation have proven to be really efficient (see [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF][START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF][START_REF] Maday | Spectral element methods for the incompressible Navier-Stokes equations, State-of-the-art surveys on computational mechanics[END_REF][START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF]).

The fictitious domain method presented in [START_REF] Bécache | A Fictitious Domain Method with Mixed Finite Elements for Elastodynamics[END_REF] is designed to take into account scattering by impenetrable obstacle. It has some of the properties mentioned above: it preserves a discrete energy and allows the use of a time step adapted to the discretization properties of the background medium. However in [START_REF] Bécache | Convergence results of the fictitious domain method for a mixed formulation of the wave equation with a Neumann boundary condition[END_REF] it is shown that it can not be compatible with high order space discretization. The space-time refinement method presented in [START_REF] Becache | Space-time mesh refinement for elastodynamics[END_REF][START_REF] Derveaux | Effective computational methods for wave propagation[END_REF] is based upon local time stepping and boundary coupling by mortar elements. It satisfies all the mentioned properties but is a non-overlapping domain decomposition method: it requires conformity between the geometry of the sub-domain boundaries, this can not be guaranteed in a generic optimization procedure where the position of a defect is not known exactly and changes between each simulations.

Following the work of [START_REF] Ghanem | Arlequin framework for multi-model, multi-time scale and heterogeneous time integrators for structural transient dynamics[END_REF], in this article we extend the Arlequin method -originally developed in [START_REF] Dhia | Problèmes mécaniques multi-échelles: la méthode Arlequin[END_REF] (see also [START_REF] Dhia | Further insights by theoretical investigations of the multiscale arlequin method[END_REF]) for static problem, in [START_REF] Dhia | Level-sets and arlequin framework for dynamic contact problems[END_REF] and [START_REF] Dhia | Level-sets fields, placement and velocity based formulations of contact-impact problems[END_REF] for dynamic problems and recently integrated in an industrial computational platform in [START_REF] Qiao | Implementation of the Arlequin method into ABAQUS: Basic formulations and applications[END_REF] -to the transient wave scattering problem while ensuring discrete energy preservation, the use of an explicit time stepping method in the coarse region and a convergence behaviour compatible with continuous quadratic finite elements. Moreover we present some variants of the Arlequin method adapted to transient wave scattering that allow less constraints in the meshes generation. Let us illustrate our approach by considering the scattering of waves by a crack (Figure 1a). A global mesh of the cracked domain is represented Figure 1b and should be generated for each crack position, size and orientation modifications. The standard Arlequin approach would consist in, first constructing a coarse homogeneous mesh of the uncracked domain, as in Figure 2a, then constructing a local fine mesh of the neighborhood of the crack, called the patch, as shown in Figure 2b, remove unnecessary elements of the coarse mesh and finally applying a matching on the region, denoted ω, as represented in Figure 2b and Figure 2c. The matching is applied by imposing the equality of the fields in a weak variational sense in ω, partitioning the different involved variational formulation in the overlap and by correcting these formulations by means of a Lagrange multiplier which implies the projection of one mesh on the other (see [START_REF] Gander | Algorithm 932: PANG: software for nonmatching grid projections in 2D and 3D with linear complexity[END_REF] for a projection algorithm of linear complexity). The obtained method is a mixed dynamic weak formulation. Therefore, as for mortar element methods for instance, after space discretization, the well-posedness of the underlying semi-discrete problem is related to the provability of a so called inf-sup condition (see [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Dhia | Mathematical analysis of the mixed Arlequin method[END_REF]): let V 1,h and V 2,h be the Hilbert discretization spaces for the solution computed in the coarse and fine mesh respectively and M h the Hilbert discretization space for the Lagrange multiplier, then the inf-sup condition is satisfied if there exists a positive constant C such that

inf m h ∈M h sup (v 1,h ,v 2,h )∈V 1,h ×V 2,h |b(v 1,h -v 2,h , m h )| m h M h v 1,h 2 V 1,h + v 2,h 2 V 1,h 1/2 ≥ C,
where the bilinear form b refer to the coupling operator. Remark that in the case illustrated in Figure 2, the boundary of the overlapping region is compatible with internal and/or external elements boundaries of the coarse and fine meshes. This allows to have a systematic approach to choose an appropriate space in which to look for the Lagrange multiplier so that the inf-sup condition is automatically satisfied: M h can be chosen as the restriction of the functions in V 1,h or V 2,h to ω (as suggested and proved in [START_REF] Dhia | Mathematical analysis of the mixed Arlequin method[END_REF]).

(a) Construction of a coarse homogeneous mesh of the domain.

(b) Fine mesh adapted to the crack and conform with ω.

(c) Reduction of the coarse mesh presented. Also conform with ω. One of the main interest of the Arlequin method is that it allows, either to move the patch or to create a new one with enhanced flexibility: the global coarse mesh could not be changed during these operations. To do so, the procedure described above must be adapted. More precisely, in the configuration depicted in Figure 3a, the crack is translated and rotated, if the same mesh of the patch is used as before then the coupling region boundary ∂ω is no more compatible with the elements' boundaries of the coarse mesh. A systematic simple choice for the space M h that guarantees that the inf-sup condition holds is no longer possible. Hence the stability and accuracy of the discretization process is not easily guaranteed. There are two possible strategies to circumvent this shortcoming:

• The coarse mesh elements can be cut or modified so that conformity between ∂ω and elements' boundaries is recovered. This corresponds to a remeshing procedure that may lead to a poor mesh in terms of quality that will require eventually a global remeshing. In the end, it is likely that small or distorted elements have to be introduced in the coarse mesh, thus penalizing the overall discretization process.

• The strategy classically used in the Arlequin method (see [START_REF] Dhia | Problèmes mécaniques multi-échelles: la méthode Arlequin[END_REF][START_REF] Dhia | Mathematical analysis of the mixed Arlequin method[END_REF][START_REF] Dhia | Further insights by theoretical investigations of the multiscale arlequin method[END_REF][START_REF] Dhia | Application of the Arlequin method to some structures with defects[END_REF] ) is to first consider that the coupling region ω is a part of the whole fine mesh, as illustrated in Figure 3b, then use a smooth degenerate partition of the energy in order to "disregard" finite elements of the coarse mesh in the zone of interest (see Figure 3c). The discretization space for the Lagrange multiplier is then chosen as a subspace of the discretization space generated from the fine mesh (i.e M h = V 2,h | ω in the case depicted). Although this approach introduces some consistency error than can be controlled, we show in section 5.2 that for our transient problem, if explicit time stepping is used, the degenerescence of the partition of the energy may imply degenerescence of the time step restriction, thus penalizing the overall scheme.

The strategy we propose will allow to circumvent the problems mentioned above. The idea is to couple the fine and coarse meshes only close to their boundaries and to remove unnecessary elements of the coarse mesh (those elements interacting with the zone of interest of the patch). The coupling is then done either using boundary coupling, i.e. imposing by mortar elements (see [START_REF] Wohlmuth | Discretization Methods and Iterative Solvers Based on Domain Decomposition[END_REF]), that the trace of the solutions are equal at the boundaries of the patch or using volume coupling, i.e. imposing that the solution are equal on the volume (as shown in Figure 4). To these constraints we associate two Lagrange multipliers that will be chosen as the trace (in case of coupling on the boundary) or the restriction (in case of coupling in a volume) of the corresponding coarse or fine finite elements spaces. With this choice it turns out that the semi-discrete inf-sup condition is automatically guaranteed for a sufficiently fine mesh refinement. The 1D convergence analysis presented in section 6 suggests that, for stability reasons, one should always prefer volume coupling on the regions where the physical coefficients vary, however this may not be possible, since the two coupling regions must not intersect. Therefore, a common scenario is to use a boundary coupling on the coarse mesh and a volume coupling on the fine mesh as illustrated by Figure 4, this allows to choose the volume coupling in a larger region.

The article is organized as follows • In section 2, we give a mathematical analysis of the Arlequin formulation at the continuous level for wave equation (proofs of this section are given in Appendix B).

• Section 3 is dedicated to the introduction of some variants of the Arlequin formulation that leads to the strategy depicted in Figure 4.

• In section 4 we briefly described the discretization in space of the methods and shows why the inf-sup condition can be verified for sufficiently small discretization parameters.

• In section 5, the construction and analysis of a conservative locally implicit coupling for the time discretization is presented and an estimate on the influence of the Arlequin coupling on the time step restriction is provided.

• Finally, in Appendix A we study analytically the influence of the Arlequin method to the time step restriction in a 1D configuration.

Formulation of the Arlequin methodology at the continuous level

The scattering problem

We are interested in the transient simulation of the propagation and scattering of a wave due to a local defect, a hole in our context. The non-defected domain Θ ⊂ R d is assumed to be homogeneous and Lipschitz regular, having defined a bounded open domain O compactly embedded in Θ, the defected domain is defined as

Ω = Θ \ O.
We assume that O may be empty (this is useful in a unidimensional setting). Being given some functions in L ∞ (Ω) denoted (ρ, µ) with ρ(x) ≥ ρ -> 0 and µ(x) ≥ µ -> 0, a.e. x ∈ Ω, we look for the solution u(x, t), for times t ∈ [0, T ] of the wave equation with regular enough source term, f ∈ W 1,1 (0, T ; L 2 (Ω)), and, for the sake of simplicity homogeneous Neumann boundary condition:

ρ ∂ 2 t u -∇ • µ ∇u = f, in Ω, t ∈ [0, T ], ∇u • n = 0, on ∂Ω, t ∈ [0, T ], (1)
with vanishing initial data and, without loss of generality, vanishing source term at the initial time

u(•, 0) = 0, ∂ t u(•, 0) = 0, f (•, 0) = 0 in Ω. (2)
The variational formulation associated to equations ( 1) and ( 2) is obtained by multiplying (1) with test function v ∈ H 1 (Ω) and integrating over Ω.

Find u(•, t) ∈ H 1 (Ω), for t ∈ (0, T ) such that

(ρ ∂ 2 t u, v) L 2 (Ω) + (µ ∇u, ∇v) L 2 (Ω) = (f, v) L 2 (Ω) ∀ v ∈ H 1 (Ω), (3)
where

(•, •) L 2 (Ω)
is the standard inner product in L 2 . It is well known that a solution of (3) exists is unique and satisfies

u ∈ C 2 ([0, T ]; L 2 (Ω)) ∩ C 1 ([0, T ]; H 1 (Ω)).
We assume now given two overlapping open sub-domains Ω 1 and Ω 2 such that the domain Ω is the union of these subdomains (see figure 5),

Ω j ⊂ Ω j ∈ {1, 2}, Ω = Ω 1 ∪ Ω 2 .
The overlapping region between Ω 1 and Ω 2 is denoted ω, it corresponds to a matching region and is assumed non-empty

ω = Ω 1 ∩ Ω 2 , ω = ∅.
Moreover, we assume that the boundary Γ of the hole O is a subset of ∂Ω 2 (see Figure 5) and that the overlapping region do not intersect Γ:

Γ = ∂O ⊂ ∂Ω 2 , ∂ω ∩ Γ = ∅.
Later in the discretization process, we will define finite element spaces on Ω 1 and Ω 2 independently and ω will be a matching region. To make things easier we assume in the following that f is compactly supported in Ω 1 \ Ω 2 . Following the Arlequin methodology developed in [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF] we look for a suitable continuous formulation that will allow a non uniform discretization process of the two sub-domains Ω 1 and Ω 2 . This is the objective of the following section.

Formulation of the wave propagation problem in overlapping domains

First we define the following Hilbert space

V = {v = (v 1 , v 2 ) | v 1 ∈ H 1 (Ω 1 ), v 2 ∈ H 1 (Ω 2 ), v 1 = v 2 in ω},
We introduce two couples of bounded functions,

α j ∈ L ∞ (Ω j ), β j ∈ L ∞ (Ω j ), j = 1 or 2, ⌦ 1 ! ⌦ ⌦ 2 ! O O Figure 5
: Typical configuration of a domain including a hole (a defect) and some heterogeneities (due to damaging) close to it. The heterogeneities and the defect are captured by Ω 2 whereas the background medium is captured by Ω 1 . The coupling domain ω is the overlapping region between Ω 1 and Ω 2 that represent partitions of the unity. Then, we define the Arlequin formulation of the problem (3):

       Find u(•, t) ∈ V, for t ∈ (0, T ), such that ∀ v ∈ V, 2 j=1 (α j ρ ∂ 2 t u j , v j ) L 2 (Ω j ) + 2 j=1 (β j µ ∇u j , ∇v j ) L 2 (Ω j ) = (f, v 1 ) L 2 (Ω 1 \Ω 2 ) , (4)
completed with zero initial conditions u j (•, 0) = 0, ∂ t u j (•, 0) = 0, for i = 1, 2 and f (•, 0) = 0.

We assume that the couples of functions (α i , β i ) satisfy the properties given below Assumption 1.

In the non-overlapping region we have,

α j = 1, β j = 1, in Ω j \ ω, j ∈ {1, 2}. (5)
In the matching domain we assume that

2 j=1 α j = 1, 2 j=1 β j = 1, in ω. ( 6
)
Moreover we assume that there exists

(α 0 , β 0 ) ∈ R 2 such that inf x∈Ω j α j (x) ≥ α 0 > 0, inf x∈Ω j β j (x) ≥ β 0 > 0, j ∈ {1, 2}. (7)
We are now in position to present a theorem for existence/uniqueness of the problem [START_REF] Bécache | A Fictitious Domain Method with Mixed Finite Elements for Elastodynamics[END_REF][START_REF] Dhia | Problèmes mécaniques multi-échelles: la méthode Arlequin[END_REF] and guarantee its equivalence with respect to problem (2, 3):

Theorem 2.1. [proved in Appendix B.1]
If Assumption 1 is satisfied then problem (2, 4) has a unique solution

(u 1 , u 2 ) ∈ C 2 ([0, T ]; L 2 (Ω 1 ) × L 2 (Ω 2 )) ∩ C 1 ([0, T ]; V ), (8)
moreover if u is the solution of (2, 3) then u = ũ where ũ ∈ H 1 (Ω) is defined by

ũ =        u 1 in Ω 1 \ ω, u 2 in Ω 2 \ ω, u j in ω, j ∈ {1, 2}.
The discretization of problem ( 4) is rather difficult since the basis function in the space V must be constructed and therefore the equality constraint of u 1 and u 2 over ω imposed strongly. Such constrained can be imposed weakly using a Lagrangian multiplier, this is the object of the next section.

Remark 2.2. It appears that condition [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF] is too strong, since by assumption we have

α 1 + α 2 = 1, β 1 + β 2 = 1, in ω.
As shown in Appendix Appendix B the energy E(t) is always positive, meaning that vanishing or negative weighting coefficient could be considered. However it will be clear with the discrete analysis presented later that degenerate partition of unity may lead to unstable explicit time discretization. Now we introduce from [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF] a reformulation of the wave propagation problem, where the coupling is imposed by means of a Lagrangian multiplier:

                     Find (u 1 (•, t), u 2 (•, t), (•, t)) ∈ H 1 (Ω 1 ) × H 1 (Ω 2 ) × H 1 (ω), for t ∈ [0, T ], such that for all (v 1 , v 2 , m) ∈ H 1 (Ω 1 ) × H 1 (Ω 2 ) × H 1 (ω) : 2 j=1 (α j ρ ∂ 2 t u j , v j ) L 2 (Ω j ) + 2 j=1 (β j µ ∇u j , ∇v j ) L 2 (Ω j ) + (v 1 -v 2 , ) H 1 (ω) = (f, v 1 ) L 2 (Ω 1 \Ω 2 ) , (9a)
(u 1 -u 2 , m) H 1 (ω) = 0, (9b)
with zero initial condition.

Theorem 2.3. [proved in Appendix B.2]
If Assumption 1 is satisfied, there exist a unique solution (u 1 , u 2 , ) of (2, 9) and it satisfies

(u 1 , u 2 , ) ∈ 2 j=1 C 2 ([0, T ]; L 2 (Ω j )) ∩ C 1 ([0, T ]; H 1 (Ω j )) × C 0 ([0, T ]; H 1 (ω)) , moreover (u 1 , u 2 ) ∈ V is also solution of (4).
Finally, to conclude this section we give a PDE interpretation of the Lagrangian multiplier . In what follows, such interpretation (equation ( 10) below) is used to construct different approximations at the discrete level of system (9).

Proposition 2.4. [proved in Appendix B.3]

If Assumption 1 is satisfied, the unique solution (u 1 , u 2 , ) of (2, 9) must verify, in a distributional sense, the following equality in the overlapping region ω:

2 -∆ = α 2 ρ ∂ 2 t u 2 -∇ • β 2 µ ∇u 2 -α 1 ρ ∂ 2 t u 1 + ∇ • β 1 µ ∇u 1 in D (ω). (10)
Note that no smoothness assumption on the β j holds therefore we can not ensure that the term ∇ • β j µ ∇u j is regular (L 2 for instance). Finally we deduce from proposition 2.4 an important result for the construction of the alternative formulation presented in section 3. To state the result we need to introduce some notations: we first define γ i and γ e the interior and exterior boundary of ω as well as two disjoint open subdomains ω i and ω e (by disjoint we mean that they have no common boundary). These set are either empty set or satisfy (see Figure 6)

γ i ⊂ ∂ω i and γ e ⊂ ∂ω e .
We denote by ω c the complementary set defined by Then we consider the following assumption:

ω c = ω \ ω i ∪ ω e . i e ! c ⌘ ! (a) Boundary-Boundary ! i e ! c (b) Volumic-Boundary ! e ! i ! c (c) Volumic-Volumic
Assumption 2. When ω c = ∅ we assume that (α 1 , β 1 ) and (α 2 , β 2 ) satisfy

α 1 = β 1 = C and α 2 = β 2 = 1 -C in ω c ,
where C is a strictly positive scalar.

Under this assumption, we obtain the following corollary from Proposition 2.4:

Corollary 2.5. [proved in Appendix B.4] If Assumption 1 and Assumption 2 are satisfied, then -∆ = 0 in D (ω c ). (11)

Alternative formulation of the wave propagation problem

In this section we suggest alternative formulations of problem [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF]. Although all the representations we give are equivalent at the continuous level they will give different discrete schemes. These alternative formulations are obtained using the property satisfied by namely equation [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF].

Variant (a): Boundary-Boundary coupling (ω i = ∅ and ω e = ∅). Using equation [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF] we have for all w ∈ H 1 (ω)

( , w) H 1 (ω) = w, ∇ • n γ i + w, ∇ • n γe
with n the outward unitary normal to ω. Then one can choose to introduce λ γ,i ≡ ∇ • n and λ γ,e ≡ ∇ • n in H -1/2 (γ i ) and H -1/2 (γ e ) respectively and substitute in (9a)

(v 1 -v 2 , ) H 1 (ω) ≡ v 1 -v 2 , λ γ,i γ i + v 1 -v 2 , λ γ,e γe . ( 12
)
For symmetry reasons we want to modify (9b) accordingly, more precisely we replace (9b) by the equation

u 1 -u 2 , µ γ,i γ i + u 1 -u 2 , µ γ,e γe = 0, (µ γ,i , µ γ,e ) ∈ H -1/2 (γ i ) × H -1/2 (γ e ). (13)
It obviously implies that u 1 = u 2 on ∂ω however it is not clear that this implies u 1 = u 2 in ω and this should be verified for consistency reasons. The consistency is preserved because of Assumption 2: choosing test functions in (9a) as

v 1 = w/C ∈ H 1 0 (ω) and v 2 = -w/(1 -C) ∈ H 1 0 (ω) , we obtain ρ ∂ 2 t (u 1 -u 2 ), w L 2 (ω) + µ ∇(u 1 -u 2 ), ∇w L 2 (ω) = 0,
hence u 1u 2 is identically 0 in ω since it satisfies an homogeneous wave equation with homogeneous Dirichlet data and zero initial data.

Variant (b): Volumic-Boundary coupling (ω i = ∅ and ω e = ∅). Again, from equation [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF], we have for all w ∈ H 1 (ω)

(w, ) H 1 (ωc) = w, ∇ • n γe + w, ∇ | ωc • n ∂ωc\γe . ( 14
)
where n is the unitary outward normal to ω c . We introduce the function i ∈ H 1 (ω i ) as the unique solution of

(w, i ) H 1 (ω i ) = (w, ) H 1 (ω i ) + w, ∇ | ωc • n ∂ω i \γ i , w ∈ H 1 (ω i ). ( 15
)
where n is again the unitary outward normal to ω c . Remark that i is well defined since the normal trace of ∇ on ∂ω i \ γ i ≡ ∂ω c \ γ e is well defined (because of ( 11)). We can substitute the last term of ( 14) using ( 15), we obtain

(w, ) H 1 (ω) = w, ∇ • n γe + (w, i ) H 1 (ω i ) .
Such equality can be used, to introduce a set of unknowns i ∈ H 1 (ω i ) and again λ γ,e ≡ ∇ • n in H -1/2 (γ e ) to substitute in (9a)

(v 1 -v 2 , ) H 1 (ω) ≡ (v 1 -v 2 , i ) H 1 (ω i ) + v 1 -v 2 , λ γ,e γe . ( 16
)
Choosing (v 1 , v 2 , m) such that v 1|ω = w ∈ H 1 0 (ω \ ω i ), v 2 = 0 and m = 0 in (9)
we have that u 1 satisfy a wave equation with no source term in ω c . The same is true for u 2 and therefore for the difference u 1u 2 . We can then replace (9b) by

(u 1 -u 2 , m i ) H 1 (ω i ) + u 1 -u 2 , µ γ,e γe = 0, (m i , µ γ,e ) ∈ H 1 (ω i ) × H -1/2 (γ e ), (17)
which implies u 1 = u 2 in ω.

Remark 3.1. In the previous configuration we have arbitrarily choose to introduce the new unknown i ∈ H 1 (ω i ), however the same treatment could have be done using a function e ∈ H 1 (ω e ). The equations (9a, 9b) can then be modified accordingly

(v 1 -v 2 , ) H 1 (ω) ≡ (v 1 -v 2 , e ) H 1 (ωe) + v 1 -v 2 , λ γ,i γ i . (18)
and

(u 1 -u 2 , m e ) H 1 (ωe) + u 1 -u 2 , µ γ,i γ i = 0, (m e , µ γ,i ) ∈ H 1 (ω e ) × H -1/2 (γ i ). (19)
Variant (c): Volumic-Volumic coupling (ω i = ∅ and ω e = ∅). Finally the case [START_REF] Qiao | Implementation of the Arlequin method into ABAQUS: Basic formulations and applications[END_REF][START_REF] Rylander | Stability of Explicit-Implicit Hybrid Time-Stepping Schemes for Maxwell's Equations[END_REF] and [START_REF] Wathen | An analysis of some element-by-element techniques[END_REF][START_REF] Doleana | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF]) can be combined. The equations (9a, 9b) can then be modified accordingly using the following substitutions

(v 1 -v 2 , ) H 1 (ω) ≡ (v 1 -v 2 , i ) H 1 (ω i ) + (v 1 -v 2 , e ) H 1 (ωe) (20)
and

(u 1 -u 2 , m i ) H 1 (ω i ) + (u 1 -u 2 , m e ) H 1 (ωe) = 0, (m i , m e ) ∈ H 1 (ω i ) × H 1 (ω e ). (21)

Space discretization

We present a numerical method to handle the problem (9) in the specific case of volumicvolumic coupling (equation ( 20) and ( 21)). The space discretization relies on a standard Galerkin approach. We use in practice the first and second order Lagrange elements. However the proposed approach is fully compatible with high order spectral elements (see [START_REF] Maday | Spectral element methods for the incompressible Navier-Stokes equations, State-of-the-art surveys on computational mechanics[END_REF], [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF], or [START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF] for more details on spectral elements) that provides mass lumping using some quadrature rule. However to simplify the analysis we do not consider the use of quadrature and assume that all the bilinear forms introduced can be exactly evaluated. Moreover, in order to show that the discrete inf-sup condition holds, we consider as an illustrative example, that for every h one can construct some standard, quasi-uniform, finite elements triangulations T 1 of Ω 1 and T 2 of Ω 2 and restrict our study to those cases when ω i and ω e are chosen to be conform with those triangulations in the sense that there exist T i ⊂ T 1 triangulation of ω i and T e ⊂ T 2 triangulation of ω e (see Figure 7). We introduce some conforming continuous finite element spaces based on the aformentionned triangulations

V 1,h ⊂ H 1 (Ω 1 ), V 2,h ⊂ H 1 (Ω 2 ), L i,h ⊂ H 1 (ω i ), L e,h ⊂ H 1 (ω e ).
The semi-discrete problem reads:

                         Find (u 1,h (•, t), u 2,h (•, t), i,h (•, t), e,h (•, t)) ∈ V 1,h × V 2,h × L i,h × L e,h such that for all (v 1,h , v 2,h , m i,h , m e,h ) ∈ V 1,h × V 2,h × L i,h × L e,h : 2 j=1 (α j ρ ∂ 2 t u j,h , v j,h ) L 2 (Ω j ) + 2 j=1 (β j µ ∇u j,h , ∇v j, ) L 2 (Ω j ) + (v 1,h -v 2,h , i,h ) H 1 (ω i ) + (v 1,h -v 2,h , e,h ) H 1 (ωe) = (f, v 1,h ) L 2 (Ω 1 \Ω 2 ) , (22a)
(u 1,h -u 2,h , m i,h ) H 1 (ω i ) + (v 1,h -v 2,h , m e,h ) H 1 (ωe) = 0. (22b)
We deduce the algebraic version of the semi discrete variational problem by defining the vectors U 1,h , U 2,h , and Λ h representing respectively the decomposition of u 1,h , u 2,h and ( i,h , e,h ) in their Lagrangian basis. We obtain

             M 1,h d 2 dt 2 U 1,h + S 1,h U 1,h + B 1,h Λ h = M 1,h F 1,h , M 2,h d 2 dt 2 U 2,h + S 2,h U 2,h -B 2,h Λ h = 0, B T 1,h U 1,h -B T 2,h U 2,h = 0, (23)
where {M j,h } are definite symmetric mass matrices and {S j,h } are semi-definite symmetric stiffness matrices. The complete error analysis of the semi-discrete and discrete problem is out of the scope of this article, we refer to [START_REF] Dhia | Mathematical analysis of the mixed Arlequin method[END_REF] for a general convergence result in the static case. Let us only mention that, following the arguments presented in [START_REF] Derveaux | Effective computational methods for wave propagation[END_REF], existence of solution of ( 23) is related to the invertibility of the matrix

B T 1,h M -1 1,h B 1,h + B T 2,h M -1 2,h B 2,h ,
which is guaranteed if and only if

ker(B 1,h ) ∩ ker(B 2,h ) = {0}.
As pointed out in [START_REF] Derveaux | Effective computational methods for wave propagation[END_REF] this last condition can be expressed as a discrete inf-sup condition. Theorem 4.1. We have for h sufficiently small inf

(m i,h ,m e,h ) sup (v 1,h ,v 2,h ) |(v 1,h -v 2,h , m i,h ) H 1 (ω i ) + (v 1,h -v 2,h , m e,h ) H 1 (ωe) | m i,h 2 H 1 (ω i ) + m e,h 2
H 1 (ωe) 1/2 v 1,h 2 H 1 (Ω 1 ) + v 2,h 2
H 1 (Ω 2 ) 1/2 ≥ C,
where C does not depend on h.

Proof. Let us consider any (m i,h , m e,h ) ∈ L i,h × L e,h and remember that L i,h and L e,h are build from the triangulations T i ⊂ T 1 and T e ⊂ T 2 . Then to build adequate v 1,h ∈ V 1,h and v 2,h ∈ V 2,h we need to assume that (which is clear for h sufficiently small) there exist

{κ 1,c k } = T 1,c ⊂ T 1 and {κ 2,c k } = T 2,c ⊂ T 2 such that γ ci := ∂ω c ∩ ∂ω i ⊂ K 1 c := k κ 1,c k ⊂ ω c , γ ce := ∂ω c ∩ ∂ω e ⊂ K 2 c := k κ 2,c k ⊂ ω c ,
as shown in Figure 7. Now we introduce the finite functional spaces of Lagrange elements:

W 1,h ⊂ H 1 (K 1 c ), W 2,h ⊂ H 1 (K 2 c ),
and notice that the traces of W 1,h|γ ci and W 2,h|γce coincide with the traces of L i,h|γ ci and L e,h|γce . Moreover, as γ ci and γ ce are closed, we can introduce m 0 i,h and m 0 e,h as the continuous extension by zero of m i,h|γ ci and m e,h|γce to the whole boundaries ∂K 1 c and ∂K 2 c . Then considering continuous lifting operators

L 1 h (•), L 2 h (•) (as introduced in [26]) we define v * 1,h = L 1 h (m 0 i,h ), v * 2,h = L 2 h (m 0 e,h ) to finally take (v 1,h , v 2,h ) ∈ V 1,h × V 2,h as: v 1,h :=    m i,h in ω i v * 1,h in K 1 c 0 elsewhere (even ω e ) and v 2,h :=    -m e,h in ω e -v * 2,h in K 2 c 0 elsewhere (even ω i )
which by definition and continuity of the lifting and the trace operator satisfy, for C 1 , C 2 ∈ R + independent of h:

(v 1,h -v 2,h , m i,h ) H 1 (ω i ) + (v 1,h -v 2,h , m e,h ) H 1 (ωe) = m i,h 2
H 1 (ω i ) + m e,h 2
H 1 (ωe) , v 1,h 2 H 1 (Ω 1 ) = v * 1,h 2 H 1 (K 1 c ) + m i,h 2 H 1 (ω i ) ≤ C 1 m i,h 2 H 1 (ω i ) , v 2,h 2 H 1 (Ω 2 ) = v * 2,h 2 H 1 (K 2 c ) + m e,h 2
H 1 (ωe) ≤ C 2 m e,h 2
H 1 (ωe) .
Then we can conclude the proof by

|(v 1,h -v 2,h , m i,h ) H 1 (ω i ) + (v 1,h -v 2,h , m e,h ) H 1 (ωe) | m i,h 2 H 1 (ω i ) + m e,h 2
H 1 (ωe) 1/2 v 1,h 2 H 1 (Ω 1 ) + v 2,h 2
H 1 (Ω 2 ) 1/2 ≥ 1 max(C 1 , C 2 ) .

Time discretization

We choose a uniform finite difference discretization of the time derivative, therefore we define a fixed time step ∆t, the unknown

(U 1,h , U 2,h , Λ h ) is approximated at each t n = n∆t by (U n 1,h , U n 2,h , Λ n h ).
The time discretization we suggest relies on the conservative Neumark scheme family. We use centered second order approximation of the second order time derivative as well as centered approximation for the term S 1,h U 1,h and S 2,h U 2,h . These approximations are parametrized by θ 1 and θ 2 respectively and corresponds to so-called theta-schemes. We obtain for all n > 0 the algebraic equations:

                                     M 1,h U n+1 1,h -2U n 1,h + U n-1 1,h ∆t 2 + S 1,h θ 1 U n+1 1,h + (1 -2θ 1 )U n 1,h + θ 1 U n-1 1,h + B 1,h Λ n+1 h + Λ n-1 h 2 = M 1,h F n 1,h , (24a)
M 2,h U n+1 2,h -2U n 2,h + U n-1 2,h ∆t 2 + S 2,h θ 2 U n+1 2,h + (1 -2θ 2 )U n 2,h + θ 2 U n-1 2,h -B 2,h Λ n+1 h + Λ n-1 h 2 = 0, (24b)
B T 1,h U n+1 1,h -B T 2,h U n+1 2,h = 0. (24c)
In these equations the parameters θ 1 and θ 2 are chosen respectively to be 0 (explicit case) or strictly positive (implicit case). In the explicit case we expect that the chosen time step will be enforced by a CFL type condition for the corresponding domain spatial discretization, for the implicit case such time step-restriction condition is relaxed but for an higher computational cost.

In the situation where mesh near the defect O is fine compared to the supposedly coarse mesh of Ω 1 it is preferable to use the time step imposes by the coarse mesh and use an implicit discretization in the fine one and a natural choice would be instance θ 1 = 0 and

θ 2 ≥ 1/4.
In what follows we present a detailed stability analysis of System [START_REF] Dhia | Level-sets and arlequin framework for dynamic contact problems[END_REF]. It is non trivial because the parameters α j and β j jump at the boundary of the coupling domain and as shown in what follows such jumps may reduce the maximum time step allowed in the case of explicit discretization.

Energy estimate

Multiplying (24a) by

U n+1 1,h -U n-1 1,h /2∆t and (24b) by U n+1 2,h -U n-1
2,h /2∆t, summing and using (24c) to get rid of the term involving the lagrangian multipliers we obtain

2 j=1 M j,h U n+1 j,h -2U n j,h + U n-1 j,h ∆t 2 + S j,h θ j U n+1 j,h + (1 -2θ j )U n j,h + θ j U n-1 j,h • U n+1 j,h -U n-1 j,h 2∆t = M 1,h F n 1,h • U n+1 1,h -U n-1 1,h 2∆t
It is then clear that the derivation of the discrete energy relation is independent of the coupling procedure, therefore with some standard arguments (see [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF] for instance) we obtain

2 j=1 E n+1/2 j -E n-1/2 j ∆t = M 1,h F n 1,h • U n+1 1,h -U n-1 1,h 2∆t , (25)
where E n+1/2 j is defined by

E n+1/2 j = 1 2 M j,h - (1 -4 θ j ) 4 ∆t 2 S j,h U n+1 j,h -U n-1 j,h 2∆t • U n+1 j,h -U n-1 j,h 2∆t + S j,h 2
U n+1 j,h + U n-1 j,h 2 • U n+1 j,h + U n-1 j,h 2 .
It is possible to show that stability of the schemes is guaranteed if E n+1/2 j is positive, this is true if the classical CFL conditions are fulfilled :

M j,h - (1 -4 θ j ) 4 ∆t 2 S j,h ≥ 0 if θ j < 1/4 . ( 26
)
Remark 5.1. Note also that equation [START_REF] Dhia | Level-sets fields, placement and velocity based formulations of contact-impact problems[END_REF] do not directly implies the stability of the problem in a well-defined norm, in [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF] it is shown that the preservation of E n+1/2 j indeed guaranteed the stability of the solution in a classical sense. The extension of such results is not obvious in the cases presented here and may be wrong if the parameters α j and β j degenerate.

Remark 5.2. The cases 0 < θ j < 1/4 corresponds to conditionally stable implicit scheme. In [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF] variants of such schemes are studied in a simpler configuration, it is shown that the choice θ j = 1/12 provide a fourth order scheme and that the consistency error constant growths with values of θ j larger and larger than 1/4. Such observations motivate the choice θ j = 1/4: it yields the unconditional stable scheme with the minimum consistency error (among the family of second order two time-steps schemes).

Estimate of the CFL Condition of the explicit/implicit scheme

In what follows we consider the specific case of θ 1 = 0 and θ 2 ≥ 1/4, which corresponds to an explicit/implicit coupling. In that case we give an analysis of the stability condition starting from [START_REF] Bernardi | A local regularization operator for triangular and quadrilateral finite elements[END_REF], which corresponds to check the positivity of E n+1/2 1 only. It seems nonoptimal since, theoretically the positivity of

E n+1/2 1 + E n+1/2 2
has to be ensure. However since the solutions are impose to be equal in a coupling domain in a weak way only some pathological may be constructed.

Local estimate of the CFL Condition. From the previous discussion it could seems that the time step restrictions are computed independently from one another. However they are coupled in a more subtle way. As mentioned before, the parameters α j and β j jump inside a single mesh element and this may decrease the CFL condition compared to the maximum time step allowed in the homogeneous case. The condition θ 2 ≥ 1/4 ensures that E n+1/2 2 is positive a sufficient condition of stability is the positivity of

E n+1/2 1 which is equivalent to M 1,h - ∆t 2 4 S 1,h ≥ 0. ( 27
)
This corresponds to check the positivity of

(α 1 ρ v 1,h , v 1,h ) L 2 (Ω 1 ) - ∆t 2 4 (β 1 µ ∇v 1,h , ∇v 1,h ) L 2 (Ω 1 ) , v 1,h ∈ V 1,h ⊂ H 1 (Ω 1 ), (28)
this is the CFL condition which acts a time step restriction.

We now present a local analysis of the CFL Condition based upon local estimate first introduced in [START_REF] Wathen | An analysis of some element-by-element techniques[END_REF]. We introduce a standard finite element triangulation {κ k } = T 1 of Ω 1 (the union of the κ k recovers Ω 1 and they do not intersect). We then defined element wise scalar ∆τ k as

∆τ k := sup δτ ∈ R + | (ρ v 1,h , v 1,h ) L 2 (κ k ) - δτ 2 4 (µ ∇v 1,h , ∇v 1,h ) L 2 (κ k ) ≥ 0, ∀ v 1,h ∈ V 1,h .
The estimation of the CFL condition that one should obtain without any coupling by the Arlequin method then reads ∆t ≤ min

k {∆τ k }. ( 29
)
We now defined the element wise scalar ∆t k as

∆t k := max δt ∈ R + | (α 1 ρ v 1,h , v 1,h ) L 2 (κ k ) - δt 2 4 (β 1 µ ∇v 1,h , ∇v 1,h ) L 2 (κ k ) ≥ 0, ∀ v 1,h ∈ V 1,h ,
and relates ∆t k to ∆τ k . This gives a local estimate of the influence of the coupling on the CFL. To do so, remark that, for all k and all positive δt,

(α 1 ρ v 1,h , v 1,h ) L 2 (κ k ) - δt 2 4 (β 1 µ ∇v 1,h , ∇v 1,h ) L 2 (κ k ) ≥ inf x∈κ k α 1 ρ 1 2 v 1,h 2
L 2 (κ k ) -sup x∈κ k β 1 δt 2 4 µ 1 2 ∇v 1,h 2 L 2 (κ k ) ,
from which we deduce that

∆t k ≥ ∆τ k inf x∈κ k α 1 sup x∈κ k β 1 . ( 30
)
An estimate of the CFL condition for the coupled problem with the Arlequin method is therefore

∆t ≤ min k    ∆τ k inf x∈κ k α 1 sup x∈κ k β 1   
Remark 5.3. Although our estimate (30) is not optimal as shown by the one-dimensional analysis done in Appendix A it shows that the time step degenerates towards 0 if α 1 vanishes, this justifies [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF] in Assumption 1. Moreover the estimate [START_REF] Gander | Algorithm 932: PANG: software for nonmatching grid projections in 2D and 3D with linear complexity[END_REF] performs well in a one-dimensional setting however numerical evidence shows that it can be sub-optimal in higher dimensions for meshes with a high level of connectivity (i.e. one element have a high number of neighboring elements).

Remark 5.4. It is possible to obtain estimate of the stability condition independent of α 1 and β 1 by using an implicit scheme only for the penalizing elements. This method is presented in [START_REF] Rylander | Stability of Explicit-Implicit Hybrid Time-Stepping Schemes for Maxwell's Equations[END_REF] for the case of the Maxwell equation, it can be seen as θ-scheme where the parameter θ depends on the space variable and allowed to implicit locally the scheme.

Remark 5.5. Note also that quadrature formulae enable the use of mass lumping techniques which leads to the construction of a diagonal mass matrix. This is no longer possible close to the intersection with the coupling domain ω. This latter drawbacks is not of significant importance since the computation of the Lagrange multipliers already involved the inversion of elliptic operator in ω.

1D convergence results

Description of the numerical configuration

In this section we present 1D numerical space-time convergence analysis of the schemes [START_REF] Dhia | Level-sets and arlequin framework for dynamic contact problems[END_REF] with the different variants presented section 3.

Continuous equations. The domain we consider is Ω = (0, [START_REF] Bécache | Convergence results of the fictitious domain method for a mixed formulation of the wave equation with a Neumann boundary condition[END_REF] with

Ω 1 = 0, 9 4
,

Ω 2 = 3 4 , 3 . ω = 3 4 , 9 4
, γ e = 3 4 , γ i = 9 4 .

The coupling domain ω i and ω e are defined up to an element of the different meshes considered:

ω e = 3 4 , 21 16 + O(h 2 ) , ω i = 27 16 + O(h 1 ), 9 4
,
where h 1 and h 2 stands for the space-step of the space discretization in Ω 1 and Ω 2 . See Figure 8 for a representation of the various domains. In these domains we solve the onedimensional scalar wave equations (1) given by

ρ ∂ 2 t u -∂ x • µ ∂ x u = 0, x ∈ Ω, t ∈ [0, T ],
with inhomogeneous Dirichlet boundary condition at x = 0 and outgoing boundary conditions at x = 3. More precisely we set

u(0, t) = 1 [0,7.5) 200 • 12 2 (3 - 2
( t+ 9 2 )( t-15 2 ) +200•4
and

∂ x u(3, t) + ∂ t u(3, t) = 0.
initial conditions are set to 0, the final time of simulation is T = 15. This configuration represents a smooth pulse traveling from left to right. Space discretization. The domain Ω 1 is assumed to represent a coarse region and therefore we use a uniform mesh of N 1 elements with space step h 1 = |Ω 1 | N 1 (N 1 will take values in {12, 24, 48, 96, 192, 384} to study the convergence of our algorithms). The domain Ω 2 represents a region where heterogeneities are supported in. To take into account these properties we use non uniform meshes built as follow: given a refinement factor R (in what follows R = 5) we start from a uniform mesh of N 2 = N 1 × R elements, the space-step being defined as

h 2 = |Ω 1 |
N 2 , then every vertices x h 2 of this mesh are slightly shifted in order to avoid effects related on mesh-uniformity. More precisely we set

x h 2 ← x h 2 + h 2 3 sin N 2 + 2 2 π x h 2 -min x h 2 max x h 2 -min x h 2
so that the end points of Ω 2 are not modified. Moreover, because we considered a 1D setting, all the variants of the Arlequin methods defined section 3 can be used therefore we are able to compare the properties of the suggested variants. For the convergence analysis we construct on the mesh of Ω 1 and Ω 2 the same P 2 standard Lagrangian finite element (with exact integration). Time discretization. We use the algorithm presented section 5 with

θ 1 = 0, θ 2 = 1 4 .
As explained section 5 the time step restriction depends only the properties of the matrices M 1,h and S 1,h relative to the space discretization in the coarse region. Therefore the time step ∆t used is computed with the maximum eigenvalue of M -1 1,h S 1,h thanks to relation [START_REF] Bernardi | A local regularization operator for triangular and quadrilateral finite elements[END_REF]. Note that we always used a safety margin by choosing a time step which is 0.95 the theoretical (computed) maximum time step allowed.

Space-time convergence analysis. We study the convergence behavior of the volume unknown u 1 and u 2 in the following norm

e 1 := u 1 -u L 2 (0,T ;H 1 (Ω 1 )) , e 2 := u 2 -u L 2 (0,T ;H 1 (Ω 2 )) ,
where the exact solution u is either computed analytically in the homogeneous case or numerically on a fine grid for the inhomogeneous case.

Homogeneous case. In a preliminiary example we consider the homogeneous case where ρ = µ = 1 as well as the simplest choice of partition of unity

α 1 = α 2 = β 1 = β 2 = 1/2.
Inhomogeneous case. Our objective is now to take into account a case where there is some heterogeneities in the background medium that are slowly oscillating with a small amplitude around a medium value. Close to a given position (here x = 9/4), we assume that the speed of waves sharply increases up to approximately three times its base value.

Such configuration is supposed to model a damage area closed to a defect. Therefore we choose ρ = 1 -sin( 8πx 3 ) 10

, µ = 1 + 10 e 500•6 2

( x+ 3 4 )( x-21 4 ) +4•500 , (31)
and a partition of unity given by

α 1 = 3 4 ρ in ω i ∪ ω e 1 2 in ω c , β 1 = 1 4 µ in ω i ∪ ω e 1 2 in ω c , α 2 = 1 -α 1 , β 2 = 1 -β 1 , (32)
as illustrated Figure 9. This choice of coefficient is driven by the objective to release any time step restriction for the explicit scheme.

Figure 9: Parameters of the 1D experiment in the heteogeneous case.

Convergence results

Numerical results presented in Figure 10 and Figure 11 show that all variants keep the appropriate convergence rate (see section 3 for the definition of the coupling variants). A comparative between the convergence curves shows that a slightly increase of the error for new variants may happen. However, this lost of accuracy is justified by the computational time saved due to the flexibility on the mesh generation and the reduction of degrees of freedom of the Lagrange multiplier used to impose the matching on the overlapping region. Notice that for the case of inhomogeneous physical coefficients the convergence curve for the "boundary-boundary" coupling (Var. (a)) is not presented. This is due to the CFL condition (28) which, on those cases when the strong variation of the physical coefficients happens inside of ω c , then leads to a much smaller time step (as shown in Section 6.3) and therefore to a convergence curve which is not comparable to the ones obtained for the other variants. Observe that the standard Arlequin strategy shows a better behavior in the homogenous case but this does not appear to be a general behavior as shows the Figure 11.

Numerical illustration of the CFL condition estimation in the inhomogeneous case

Figure 12 shows computations of the ∆τ k and ∆t k (used in local estimate of the CFL without and with coupling respectively) as detailed in section 5.2 for inhomogeneous case and second order elements. One can see that in both cases the local estimate of the CFL will be driven by the elements at the boundary of the coupling region. Global numerical computations of the maximal time step allowed show that the coupling is improving the situations: the maximum time step allowed is circa 0.010 if no coupling technique is used (for the case N 1 = 24) but is around to 0.023 when the volume/volume coupling strategy is used (the coefficient of the configuration are given by (31) and (32)).

Finally note, although the "boundary-boundary" coupling (Var. (a)) offers the greatest flexibility in term of mesh generation (as explained in the introduction of this article) and is also the less computational intensive, it may not be used because of heterogeneities that would restrict to much the maximum time step allowed. Therefore the most robust strategy is the "volumic-boundary" coupling (Var. (b)) which is a reasonable alternative in term of stability and computational cost. Clearly we should tend to this situation when γ → 0 (α 1 = α 2 and h is irrelevant), when h → 0 (we have α = 1 in almost all the element) and when h → 1 whenever γ = 1 (in this situation we have α(x) = a > 0 in almost all the element). A singular situation is encountered when γ = 1 (that is, a = 0) since for h → 1 the mass matrix M vanishes. More precisely we have 0.8 0.8 0.8 0.8 1 1 0.9 0.9

lim h →0 λ(γ, h) = 12, lim γ →0 λ(γ, h) = 12, lim h →1 λ(γ, h) =    12, if γ ∈ [0, 1) +∞, if γ = 1 lim γ →1 λ(γ, h) = 12 (h -1) 2 .
(a) The graph of λ(γ, h) for (γ, h) ∈ (0, 0.9) × (0, 1). (b) The graph of λ 2 (1/2, h) for h ∈ (0, 1).

. In the particular case in which γ = 1/2 (that is of our interest; see figure A.13b) we observe that λ(1/2, h), h ∈ [0, 1] lays on the interval [10.5, 14]. The CFL condition would not be very penalized (Worst case scenario would be to have a time step equal to 12/14 0.926 times the maximal time step allowed without coupling.

moreover if u is the solution of (2, 3) then u = ũ where ũ ∈ H 1 (Ω) is defined by

ũ =        u 1 in Ω 1 \ ω, u 2 in Ω 2 \ ω, u j in ω j ∈ {1, 2}.
Proof. It is rather standard to obtain existence and uniqueness result for the the problem (2, 4) and the proof will not be reproduced here. Note however that part of the analysis is based upon the obtention of the following energy identity: by setting v i = ∂ t u i in (4) we find, 1 2

d dt E(t) = (f, ∂ t u 1 ) L 2 (Ω 1 \Ω 2 ) ⇒ E(t) ≤ ρ -1 - t 0 f (s) L 2 (Ω 1 ) ds (B.2)
where the energy E is defined by

E(t) = 2 j=1 (α j ρ ∂ t u j , ∂ t u j ) L 2 (Ω j ) + 2 j=1 (β j µ ∇u j , ∇u j ) L 2 (Ω j ) , (B.3)
which is also, taking into account relations [START_REF] Dhia | Mathematical analysis of the mixed Arlequin method[END_REF][START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF] of Assumption 1:

E(t) = 2 j=1 (ρ ∂ t u j , ∂ t u j ) L 2 (Ω j \ω) + 2 j=1
(µ∇u j , ∇u j ) L 2 (Ω j \ω)

+ ρ ∂ t u 1 , ∂ t u 1 L 2 (ω) + µ∇u 1 , ∇u 1 L 2 (ω) .
Note that precise energy estimate of the solution with respect to the source term can easily be deduced from (B.2). The extra regularity given in (B.1) is obtained by differentiating once with respect to time formulation (4): thanks to (2) (and since f (•, 0) = 0), the unknown ∂ t u is solution of the same formulation with zero initial condition and source term ∂ t f . The existence and uniqueness proof above can be repeated, it leads to the required extra regularity property. Now remark that ũ ∈ H 1 (Ω) by construction. We now prove that ũ satisfies the variational formulation (3) and conclude by uniqueness of the solution of the standard wave equation [START_REF] Bécache | Convergence results of the fictitious domain method for a mixed formulation of the wave equation with a Neumann boundary condition[END_REF]. Using the decomposition of the solution into subdomains and the value of ũ on each of these subdomains, we have, for all v ∈ H 1 (Ω), and any k ∈ {1, 2}

(ρ ∂ 2 t ũ, v) L 2 (Ω) + (µ ∇ũ, ∇v) L 2 (Ω) = (ρ ∂ 2 t u k , v) L 2 (ω) + (µ ∇u k , ∇v) L 2 (ω) + 2 j=1 (ρ ∂ 2 t u j , v) L 2 (Ω j \ω) + 2 j=1 (µ ∇u j , ∇v) L 2 (Ω i \ω) . (B.4) j=1 (α j ρ ∂ 2 t u j , v) L 2 (Ω j \ω) + 2 j=1
(β j µ ∇u j , ∇v) L 2 (Ω j \ω)

+ (α 1 + α 2 ) ρ ∂ 2 t u k , v L 2 (ω) + (β 1 + β 2 ) µ ∇u k , ∇v L 2 (ω) = (f, v) L 2 (Ω 1 \Ω 2 ) . (B.5)
We can use the assumptions ( 5) and ( 6) on the α j 's and β j 's to simplify (B.5), then the right hand side of (B.4) can be substituted: we find that, for all v ∈ H 1 (Ω),

(ρ ∂ 2 t ũ, v) L 2 (Ω) + (µ ∇ũ, ∇v) L 2 (Ω) = (f, v) L 2 (Ω 1 \Ω 2 ) ,
and so, we recover (3), which implies u = ũ by uniqueness of the solution of the wave equation with the same source term and zero initial conditions (obviously ũ(•, 0) = 0 and ∂ t ũ(•, 0) = 0 in Ω).

Theorem Appendix B.2. [proof of Theorem 2.3]

If Assumption 1 is satisfied, there exist a unique solution (u 1 , u 2 , ) of (2, 9) and it satisfies

(u 1 , u 2 , ) ∈ 2 j=1 C 2 ([0, T ]; L 2 (Ω j )) ∩ C 1 ([0, T ]; H 1 (Ω j )) × C 0 ([0, T ]; H 1 (ω)) , moreover (u 1 , u 2 ) ∈ V is also solution of (4).
Proof. The last part of the statement is easily proved. From the second equation of [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF] we have that (u 1 , u 2 ) ∈ V and choosing test function (v 1 , v 2 ) ∈ V the first equation of ( 9) gives (4). The rest of the proof is dedicated to proving existence/uniqueness/regularity of solutions of problem [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF].

The notation (•, •) L 2 (Ω) denotes the complex hermitian inner product in L 2 (Ω).

Extension of the source term. We denote by f e ∈ W 1,1 (R + , L 2 (Ω 1 )) the extension of the source term f (•, t) for t > T such that f e (•, t) = 0 for t > 2T . Note that with either f e chosen or f the solutions should coincide for t ≤ T . Therefore existence/uniqueness and estimates will be obtained with source term f e . To obtain a preliminary result we assume that f e ∈ C ∞ 0 (R + , L 2 (Ω 1 )).

Existence/uniqueness/estimations in Laplace domain. We introduce the Laplace transform L of any time dependent function h vanishing at t = 0 as

L(h(t)) = h(s) = +∞ 0 h(t) e -st dt, s = j ξ + η, η > 0.
Note that, the source f e , the Lagrange multiplier as well as u 1 and u 2 and their first two derivatives vanish at t = 0. The variational formulation of (2, 9) written for t ∈ R + becomes after transformation in Laplace domain: find

( u 1 , u 2 , ) ∈ H 1 (Ω 1 ) × H 1 (Ω 2 ) × H 1 (ω) such that for all ( v 1 , v 2 , m) ∈ H 1 (Ω 1 ) × H 1 (Ω 2 ) × H 1 (ω):              s s 2 2 j=1 (α j ρ u j , v j ) L 2 (Ω j ) + s 2 j=1 (β j µ ∇ u j , ∇ v j ) L 2 (Ω j ) +s ( v 1 -v 2 , ) H 1 (ω) = s ( f e , v 1 ) L 2 (Ω 1 \Ω 2 ) , s ( u 1 -u 2 , m) H 1 (ω) = 0, (B.6)
To proceed with the analysis we introduce the bilinear forms a(s) :

H 1 (Ω 1 ) × H 1 (Ω 2 ) × H 1 (Ω 1 ) × H 1 (Ω 2 ) → C, b(s) : H 1 (Ω 1 ) × H 1 (Ω 2 ) × H 1 (ω) → C,
where a(s) is defined by

a s ; ( u 1 , u 2 ), ( v 1 , v 2 ) = s s 2 2 j=1 (α j ρ u j , v j ) L 2 (Ω j ) + s 2 j=1 (β j µ ∇ u j , ∇ v j ) L 2 (Ω j ) . (B.7) and b(s ; ( v 1 , v 2 ), µ) = s ( v 1 -v 2 , m) H 1 (ω) .
We can then recast our problem into a mixed form:

   a s ; ( u 1 , u 2 ), ( v 1 , v 2 ) + b(s ; ( v 1 , v 2 ), ) = s ( f e , v 1 ) L 2 (Ω 1 \Ω 2 ) , b(s ; ( u 1 , u 2 ), m) = 0. (B.8)
To guarantee existence/uniqueness of the solution we need to check the inf-sup condition for an appropriate norm. As in [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF] we define a s-dependent H 1 -like norm:

v j 2 Ω j = |s| 2 v j 2 L 2 (Ω j ) + ∇ v j 2 L 2 (Ω j ) .
We have that a(s) is continuous for the composed norm

|a s ; ( u 1 , u 2 ), ( v 1 , v 2 ) | ≤ C |s| ( u 1 2 Ω 1 + u 2 2 Ω 2 ) 1 2 ( v 1 2 Ω 1 + v 2 2 Ω 2 ) 1 2 , (B.9)
with C a positive scalar depending only on the L ∞ norm of ρ, µ, α j and β j . The coercivity of a(s) is also guaranteed

|a s ; ( u 1 , u 2 ), ( u 1 , u 2 ) | ≥ c η ( u 1 2 Ω 1 + u 2 2 Ω 2 ), (B.10)
where c is a positive scalar that depends only on α 0 , β 0 , ρ and µ. Moreover b(s

; ( v 1 , v 2 ), m) is continuous |b(s ; ( v 1 , v 2 ), m)| ≤ |s| ( v 1 2 Ω 1 + v 2 2 Ω 2 ) 1 2 m H 1 (ω) . (B.11)
In [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF], Theorem 2.25, a statement of existence/uniqueness results is given for mixed problems of the form (B.8). It relies on: the continuity properties (B.9) and (B.11), the coercivity property (B.10), and an inf-sup condition in complex domain that still need to be proved. More precisely the inf-sup condition is: there exists k(s) a strictly positive constant such that

sup ( v 1 , v 2 )∈H 1 (Ω 1 )×H 1 (Ω 2 ) |b(s ; ( v 1 , v 2 ), m)| ( v 1 2 Ω 1 + v 2 2 Ω 2 ) 1 2 ≥ k(s) m H 1 (ω) . (B.12)
Adequate choice of v 1 and v 2 in (B.12) will enable us to prove this inequality. Following the study of [START_REF] Dhia | Mathematical analysis of the mixed Arlequin method[END_REF], we set v 2 = 0, and

v 1 = L 1 ( m) ∈ H 1 (Ω 1 )
, where L 1 is a continuous lifting operator (see [START_REF] Mclean | Strongly Elliptic systems and Boundary Integral equations[END_REF] theorem 4.10):

v 1 = m in ω, -∆ v 1 = 0 in Ω 1 \ ω, ∇ v 1 • n = 0 on ∂Ω 1 \ ∂ω.
For this specific choice of ( v 1 , v 2 ), inequality (B.12) becomes

|s| m 2 H 1 (ω) v 1 Ω 1 ≤ sup ( v 1 , v 2 )∈H 1 (Ω 1 )×H 1 (Ω 2 ) |b(s ; ( v 1 , v 2 ), m)| ( v 1 2 Ω 1 + v 2 2 Ω 2 ) 1 2
.

Since the Lifting operator is continuous we have

v 1 H 1 (Ω 1 ) ≤ L m H 1 (ω)
, where L is a positive scalar depending only on Ω 1 . As a consequence Existence and uniqueness for problem (B.8) is therefore guaranteed. Moreover, using standard results on mixed problems (see [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]), it can be shown that the solution ( u 1 , u 2 , ) of (B.6) satisfies the estimates Using Plancherel theorem we find that the unique solution (u 1 , u 2 , ) of (2, 9) satisfies e -η t ∂ n+1 t u j (t) L 2 (Ω j ) ∈ L 2 (R + ), e -η t ∂ n t ∇u j (t) L 2 (Ω j ) ∈ L 2 (R + ), e -η t H 1 (ω) ∈ L 2 (R + ), for n = 0 and n = 1. This implies in particular ∂ t u j ∈ L 2 0, T ; L 2 (Ω j ) , u j ∈ L 2 0, T ; H 1 (Ω j ) and ∈ L 2 0, T ; H 1 (ω) .

v 1 Ω 1 ≤ L 1 + |s| 2 m H 1 (ω) , (
( u 1 2 Ω 1 + u 2 2 Ω 2 ) 1 2 ≤ 1 c η f e L 2 (
as well as ∂ 2 t u j ∈ L 2 0, T ; L 2 (Ω j ) and ∂ t u j ∈ L 2 0, T ; H 1 (Ω j ) . By injection one can deduce (see [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF] Chap. XVIII) that

(u 1 , u 2 ) ∈ 2 j=1 C 1 ([0, T ]; L 2 (Ω j )) ∩ C 0 ([0, T ]; H 1 (Ω j )) .
By repeating the same arguments and observing that e -η t ∂ 2 t f e (t) L 2 (Ω 1 ) ∈ L 2 (R + ) after differentiating in time the problem [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF], we can also show that

(u 1 , u 2 ) ∈ 2 j=1
C 2 ([0, T ]; L 2 (Ω j )) ∩ C 1 ([0, T ]; H 1 (Ω j )) and ∈ H 1 0, T ; H 1 (ω) .

Note that the last relation implies ∈ C 0 [0, T ]; H 1 (ω).

Source terms with minimal regularity. We now show that there exists a unique solution of problem [START_REF] Bécache | A Fictitious Domain Method with Mixed Finite Elements for Elastodynamics[END_REF][START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF] with the adequate regularity without assuming that f e is in C ∞ 0 (R + ; L 2 (Ω 1 )). By density, there exists a sequence of compactly supported functions f m e ∈ C ∞ 0 (R + ; L 2 (Ω 1 )) such that f m e → f e in W 1,1 (R + ; L 2 (Ω 1 )) with f e (•, 0) = 0. The associated solutions are denoted (u m 1 , u m 2 , m ), we have, as show previously

(u m 1 , u m 2 ) ∈ 2 j=1 C 2 ([0, T ]; L 2 (Ω j )) ∩ C 1 ([0, T ]; H 1 (Ω j ))
as well as the energy estimate (obtained using energy techniques on smooth solutions and Gronwall's lemma in time domain) ∂ k+1 t u j (t) L 2 (Ω j ) + ∇∂ k t u j (t) L 2 (Ω j ) + (t) H 1 (ω) .

We have We have that (u m 1 , u m 2 , m ) → (u 1 , u 2 , ) in W . Then one can show that (u 1 , u 2 , ) is indeed solution of (9) by writing [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF] for (u m 1 , u m 2 , m ) then going to the limit m → +∞.

(u m 1 -u n 1 , u m 2 -u n 2 , m -n ) W ≤ C sup

  (a) Cracked domain.(b) Adapted mesh.

Figure 1 :

 1 Figure 1: Global adapted mesh for a cracked domain.

Figure 2 :

 2 Figure 2: Standard Arlequin decomposition for a scattering problem.

  (a) New crack position, but same mesh presented in Figure2a.(b) Mesh presented in Figure2badapted to new crack position.(c) There is no valid coupling region ω conform with one of the meshes (In red the region where the coefficient degenerates).

Figure 3 :

 3 Figure 3: The standard Arlequin strategy applied to wave scattering.

  (a) Fine adapted mesh. It can be used for any new crack position. (b) Reduction of the coarse mesh presented in Figure 2a.

Figure 4 :

 4 Figure 4: New decomposition strategy that can be adapted to any position of the crack.

Figure 6 :

 6 Figure 6: Representation of the typical domain decompositions considered for the overlapping region ω.

(a) Triangulation T 1

 1 of Ω 1 and its sub-triangulations T i of ω i and T 1,c in the proof of Theorem 4.1 (b) Triangulation T 2 of Ω 2 and its sub-triangulations T e of ω e and T 2,c in the proof of Theorem 4.1 (c) Triangulations in ω. Notice that T 1,c and T 2,c could intersect. However T i and T e can not.

Figure 7 :

 7 Figure 7: Standard triangulations so that the inf-sup condition is satisfied.

Figure 8 :

 8 Figure 8: Representation of the configuration and geometry of the problem considered.

Figure 10 :

 10 Figure 10: Convergence curves for the case of homogeneous coefficients. The optimal time step is the same for all the variants. Second order finite elements are used.

Figure 11 :

 11 Figure 11: Convergence curves for the case of inhomogeneous coefficients. Var. (a) results are not comparable due to a much more restrictive CFL condition. Second order finite elements are used.

Figure 12 :

 12 Figure 12: Histogram representing the local estimation of the CFL for each element of the mesh with second order elements in a inhomogeneous case (µ varies strongly in Ω 2 ).

Figure A. 13 :

 13 Figure A.13: Graph of the second eigenvalue λ(γ, h) when a = b.

2 H 1 L 1 + |s| 2 ≤ |s| m 2 H 1 L 1 + |s| 2 .

 21122112 B.13) and from inequality (B.13) we can deduce the following results|s| m (ω) (ω) v 1 Ω 1 ⇒ (B.12) holds with k(s) = |s|Remark that inequality (B.12) also holds with k(s) ≡ k(η) = L -1 min(η/2, 1/2) since |s| 1 + |s| 2 ≥ min(|s|/2, 1/2) ≥ min(η/2, 1/2).

0 ∂ 2 j=1C 2 ( 1 k=0 2 j=1

 02212 (t) L 2 (Ω j ) + ∇∂ k t u m j (t) L 2 (Ω j ) ≤ C e (s) L 2 (Ω 1 ) ds,where, in what follows, C is a positive scalar depending only on T , the coefficients (ρ, µ) and the lifting operator L 1 . This implies the following estimate on the Lagrange multiplier (simply choose v 1 = L 1 m and v 2 = 0 in (9) as done above in the Laplace domain)m (t) H 1 (ω) ≤ C f m e (t) L 2 (Ω 1 ) + C t t f m e (s) L 2 (Ω 1 ) ds.We now construct a sequence(u m 1u n 1 , u m 2u n 2 , mn ) in the Banach space W = [0, T ]; L 2 (Ω j )) ∩ C 1 ([0, T ]; H 1 (Ω j )) × C 0 ([0, T ]; H 1 (ω))equipped with the norm(u 1 , u 2 , ) W = sup t∈[0,T ]

L 2 (Ω 1 ) + t 0 ∂

 20 t (f m ef n e )(s) L 2 (Ω 1 ) ds , (B.16)which implies that our sequence is a Cauchy sequence since the right hand side of the previous equation is also a Cauchy sequence (f m e → f e in W 1,1 (R + ; L 2 (Ω 1 ))) and supt∈[0,T ] (f m ef n e )(t) L 2 (Ω 1 ) ≤ C f m ef ne W 1,1 (R;L 2 (Ω 1 ))

  Ω 1 ) , H 1 (ω) ≤Existence/uniqueness in time domain. Estimate (B.[START_REF] Maday | Spectral element methods for the incompressible Navier-Stokes equations, State-of-the-art surveys on computational mechanics[END_REF] is the key estimate to obtain existence and uniqueness of the solution in time domain. Following the standard arguments of[START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF] Chap. XVI, observe that for any causal time dependent function h(t) which n first derivatives vanish at the origin we have s n h = F(e -η t ∂ n t h(t)) where F is the Fourier transform from the time variable t to the frequency variable ξ and η is assume fixed and strictly positive. Now, since by assumption f e is smooth we havee -η t f e (t) L 2 (Ω 1 ) ∈ L 2 (R + ), e -η t ∂ t f e (t) L 2 (Ω 1 ) ∈ L 2 (R + )and since f e and ∂ t f e vanish at the origin we have by application of Plancherel theorem Therefore from estimate (B.14) and (B.15) we see that u j (s) Ω j and H 1 (ω) are square integrable functions of ξ:

	1 k(η) f e L R 1 + C|s| c η f e 2 L 2 (Ω 1 ) + s f e 2 L 2 (Ω 1 ) dξ < +∞.	(B.15)

2

(Ω 1 ) (B.

[START_REF] Maday | Spectral element methods for the incompressible Navier-Stokes equations, State-of-the-art surveys on computational mechanics[END_REF]

where C and c are respectively the same constant as in (B.9) and (B.10).

R (|s| 2 + |s| 4 ) u j (s) 2 L 2 (Ω j ) dξ < +∞, R (1 + |s| 2 ) ∇ u j (s) 2 L 2 (Ω j ) dξ < +∞ and R (s) 2 H 1 (ω) dξ < +∞.

Appendix A. Local analysis of the CFL in a 1D setting

In this section we exhibit the impact of the Arlequin methodology on the stability condition of the method. To do so we consider the one dimensional case. By the local arguments used in section 5.2 we only need to consider a reference element: the interval [0, 1]. For h ∈ [0, 1] we consider the functions

This configuration correspond to a situation where the coupling region boundary lies inside an element. The coefficients α j and β j jump and therefore modify the properties of the mass and stiffness bilinear forms hence the CFL condition when an explicit scheme is used. If using P 1 continuous elements, the basis functions are given by

To estimate the global CFL condtion we want to find the largest scalar ∆τ such that

Such problem can be recast as a generalized eigenvalue problem. We have ∆τ = 2/ √ max λ where the λ are defined as the eigenvalues of the pair

where M and K are respectively 2 × 2 local mass and stiffness matrix. We have

The first eigen-pair of (A.2) is given by ( 0, [1, 1] T ) while the second eigenvalue is given by

A particular case that is of our interest occurs when a = b. In that case we set γ = 1a, we get

In figure A.13a we represent the graph of λ(γ, h). It is worth to note that the larger the eigenvalue is, the smaller the time step needs to be chosen. Moreover, setting α(x) = 1 we would have found λ(0, h) = 12.

Proposition Appendix B.3. [proof of Proposition 2.4]

If Assumption 1 is satisfied, the unique solution (u 1 , u 2 , ) of (2, 9) must verify, in a distributional sense, the following equality in the overlapping region ω:

Proof. By setting v 1 = w| Ω 1 and v 2 = -w| Ω 2 in ( 9) with w ∈ D(ω) (extended by 0 in Ω), we find

From equation (B.18) we can deduce the partial differential equation satisfied by (t) in ω.

We see that

where ω c is included in ω (the functions being extended by 0), v 2 = 0 and m = 0 in (9) we have

and setting v 2 = w ∈ H 1 0 (ω c ), v 1 = 0 and µ = 0 we find

Adding the last two equations (B.21) and (B.22) we get

We have shown that: If u 1 = u 2 in any domain ω c ⊂ ω and Assumption 1 is satisfied then