Y Boubendir

X Antoine

C Geuzaine

A Quasi-Optimal Non-Overlapping Domain Decomposition Algorithm for the Helmholtz Equation

This paper presents a new non-overlapping domain decomposition method for the Helmholtz equation, whose effective convergence is quasi-optimal. These improved properties result from a combination of an appropriate choice of transmission conditions and a suitable approximation of the Dirichlet to Neumann operator. A convergence theorem of the algorithm is established and numerical results validating the new approach are presented in both two and three dimensions.

Introduction

In this paper, we are interested in non-overlapping Domain Decomposition Methods (DDMs) for the Helmholtz equation. Such methods were introduced by Lions [START_REF] Lions | On the Schwarz Alternating Method III: A Variant for Non Overlapping Subdomains, Third International Symposium on Domain Decomposition Methods for Partial Differential Equations[END_REF] for the Laplace equation and extended to the Helmholtz equation by Després [START_REF] Després | Méthodes de Décomposition de Domaine pour les Problèmes de Propagation d'Ondes en Régime Harmonique[END_REF][START_REF] Després | Domain Decomposition Method and the Helmholtz Problem[END_REF][START_REF] Després | Domain Decomposition Method and the Helmholtz Problem (part II)[END_REF]. Essentially, the method consists in combining the continuity conditions (of the field and its normal derivative) on the artificial interfaces between subdomains, in order to obtain Robin boundary conditions and to solve the overall problem by iterating over the subdomains [START_REF] Quarteroni | Domain Decomposition Methods for Partial Differential Equations[END_REF][START_REF] Toselli | Domain Decomposition Methods Algorithms and Theory[END_REF][START_REF] Stupfel | A Hybrid Finite Element and Integral Equation Domain Decomposition Method for the Solution of the 3-D Scattering Problem[END_REF]. Robin conditions (also called absorbing or impedance boundary conditions) are chosen to couple the subdomains because using the natural conditions leads to divergent iterative algorithms [START_REF] Bendali | A FETI-like Domain Decomposition Method for Coupling FEM and BEM in Large-Size Problems of Acoustic Scattering[END_REF].

Improving the convergence properties of the iterative process constitutes the key in designing effective algorithms, in particular in mid and high frequency. The optimal convergence is obtained by defining the transmission conditions, on each interface, using the Dirichlet-to-Neumann (DtN) operator corresponding to the problem sharing the same interface [START_REF] Roux | Optimization of Interface Operator Based on Algebraic Approach[END_REF][START_REF] Roux | Approximation of Optimal Interface Boundary Conditions for Two-Lagrange Multiplier FETI Method[END_REF]. This however leads to a very expensive procedure in practice. The first contribution following the original method of Després, using particular non-local transmission conditions, was made in [START_REF] Collino | Domain Decomposition Method for Harmonic Wave Propagation: a General Presentation[END_REF] with a relative success regarding the effective convergence. A great variety of techniques based on local transmission conditions have also been proposed to improve the convergence: these include the class of FETI-H methods [START_REF] Farhat | A Two-Level Domain Decomposition Method for the Iterative Solution of High Frequency Exterior Helmholtz Problems[END_REF][START_REF] Farhat | Two-Level Domain Decomposition Methods with Lagrange Multipliers for the Fast Iterative Solution of Acoustic Scattering Problems[END_REF][START_REF] Farhat | A Domain Decomposition Method for Discontinuous Galerkin Discretizations of Helmholtz Problems with Plane Waves and Lagrange Multipliers[END_REF][START_REF] De La Bourdonnaye | A Non Overlapping Domain Decomposition Method for the Exterior Helmholtz Problem[END_REF], the optimized Schwartz approach [START_REF] Gander | Optimized Schwarz Methods without Overlap for the Helmholtz Equation[END_REF], and the evanescent modes damping algorithm [START_REF] Boubendir | An Analysis of the BEM-FEM Non-Overlapping Domain Decomposition Method for a Scattering Problem[END_REF][START_REF] Boubendir | Coupling of a Non-Overlapping Domain Decomposition Method for a Nodal Finite Element Method with a Boundary Element Method[END_REF][START_REF] Bendali | A FETI-like Domain Decomposition Method for Coupling FEM and BEM in Large-Size Problems of Acoustic Scattering[END_REF]. However, the related impedance operators do not accurately approximate the exact DtN operator on all the modes of the solution, which makes the resulting iterative methods suboptimal.

In this paper, we propose a new square-root based transmission condition, localized using Padé approximants, which accurately approximates the DtN operator and allows to design an algorithm with quasi-optimal convergence properties. Indeed, we will show that the rate of convergence is optimal on the evanescent modes and is significantly improved compared to current techniques for the remaining modes. In the effective convergence, this results in a DDM independent of the wavenumber as well as the mesh discretization. Moreover, we will see how the resulting approximate DtN operator is easy to implement in a basic finite element solver.

Γ Σ ij = Σ ji Γ ∞ n i Ω i n j Ω j
The paper is organized as follows. In Section 2, we introduce the scattering problem as well as the nonoverlapping DDM. We present in the third Section a non-local square-root operator which approximates the exact DtN transmission operator. Section 4 develops a convergence analysis for this approximate transmission condition on a model problem. Section 5 details the complex Padé approximation of the square-root operator to get a local representation. Section 6 presents the finite element implementation of the resulting DDM. Numerical results on both two-and three-dimensional problems are presented in Section 7.

Scattering Problem and Non-Overlapping Domain Decomposition Method

Let us consider the three-dimensional time-harmonic scattering problem of an incident acoustic wave by an obstacle K. We want to compute the scattered field u solution to the exterior Helmholtz equation with a Dirichlet boundary condition1 :

     ∆u + k 2 u = 0 in R 3 \K, u = f on Γ = ∂K, lim |x|→∞ |x|(∂ |x| u -ıku) = 0. (1)
The boundary data f is fixed by a plane wave: f = -e ıkα•x , with x = (x 1 , x 2 , x 3 ) ∈ R 3 and ı = √ -1. The incidence angle α is normalized on the unit sphere (|α| = 1) and k denotes the wavenumber, related to the wavelength λ of the incident wave through k = 2π/λ. The last equation of system (1) is the Sommerfeld radiation condition at infinity, which imposes that the scattered wave is outgoing.

To solve (1), we combine Absorbing Boundary Conditions (ABCs) with Lions-Després' non-overlapping domain decomposition method. The ABC method consists in truncating the infinite domain by introducing a fictitious boundary Γ ∞ to get a bounded computational region. Then, system (1) is approximated by

     ∆u + k 2 u = 0 in Ω, u = f on Γ, ∂ n u + Bu = 0 on Γ ∞ , ( 2
)
where Ω is the bounded domain enclosed by the fictitious boundary Γ ∞ and Γ (see Figure 1) and where the operator B represents an approximation of the DtN operator (for example B = -ık) on Γ ∞ . The vector n is the outwardly directed unit normal to Γ ∞ . Let us remark that all what follows can be adapted to the case of other truncation techniques, like e.g. perfectly matched layers (PMLs) [START_REF] Bermudez | An Exact Bounded Perfectly Matched Layer for Time-Harmonic Scattering Problems[END_REF][START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF][START_REF] Singer | A Perfectly Matched Layer for the Helmholtz Equation in a Semi-Infinite Strip[END_REF][START_REF] Turkel | Absorbing PML Boundary Layers for Wave-Like Equations[END_REF]. The first step of the Lions-Després domain decomposition method [START_REF] Lions | On the Schwarz Alternating Method III: A Variant for Non Overlapping Subdomains, Third International Symposium on Domain Decomposition Methods for Partial Differential Equations[END_REF][START_REF] Després | Méthodes de Décomposition de Domaine pour les Problèmes de Propagation d'Ondes en Régime Harmonique[END_REF] consists in splitting Ω into several subdomains Ω i , i = 1, . . . , N dom , such that (see Figure 1):

• Ω = N dom i=1 Ω i (i = 1, . . . , N dom ), • Ω i ∩ Ω j = ∅, if i = j, (i, j = 1, . . . , N dom ), • ∂Ω i ∩ ∂Ω j = Σ ij = Σ ji (i, j = 1, . . . , N dom )
is the artificial interface separating Ω i and Ω j as long as its interior Σ ij is not empty.

The second step consists in reducing the solution of the initial problem (2) by solving the local transmission problems

       ∆u (n+1) i + k 2 u (n+1) i = 0 in Ω i , u (n+1) i = f i on Γ i , ∂ ni u (n+1) i + Bu (n+1) i = 0 on Γ ∞ i , (3a)
∂ ni u (n+1) i + Su (n+1) i = g (n) ij on Σ ij , (3b)
and then in forming the quantities to be transmitted through the interfaces

g (n+1) ij = -∂ nj u (n+1) j + Su (n+1) j = -g (n) ij + 2Su (n+1) j on Σ ij , (4)
where

u i = u| Ωi , n i (resp. n j ) is the outward unit normal to the boundary of Ω i (resp. Ω j ), i = 1, . . . , N dom , j = 1, . . . , N dom , Γ i = ∂Ω i ∩ Γ, Γ ∞ i = ∂Ω i ∩ Γ ∞
and S an invertible operator. Let us note that the boundary condition on Γ i (resp. Γ ∞ i ) does not take place if the interior of ∂Ω i ∩ Γ (resp. ∂Ω i ∩ Γ ∞ ) is the empty set. We will assume in all that follows that the DDM is well-posed, in the sense that each subproblem (3a)-(3b) is itself well-posed, i.e., away from interior resonances.

Solving at each step all the local transmission problems through (3)-( 4) may be recast as one application of the iteration operator A :

× N dom i,j=1 L 2 (Σ ij ) → × N dom i,j=1 L 2 (Σ ij ) [11] defined by g (n+1) = Ag (n) + b, (5)
where g (n) is the set of boundary data (g

(n)
ij ) 1≤i,j≤N dom , and b is given by the Dirichlet boundary condition. Therefore, (3)-( 4) can be seen as an iteration of the Jacobi method (or fixed point iteration) applied to the linear system:

(I -A)g = b, ( 6
)
where I is the identity operator. Following this idea, any Krylov solver could also be applied to solve this equation. Indeed, we will see that the GMRES algorithm significantly improves the iterative process over the successive approximations procedure (5).

Square-Root Based Transmission Condition

It is well-known that the convergence of the domain decomposition method for scattering problems strongly depends on the choice of the transmission operator S. Indeed, to each choice of operator S corresponds an iteration operator A with particular spectral properties. Several techniques were developed to improve the convergence by modifying the original algorithm introduced by Després in [START_REF] Després | Méthodes de Décomposition de Domaine pour les Problèmes de Propagation d'Ondes en Régime Harmonique[END_REF]. In this last reference, the low-order approximation of the DtN operator

S 0 u = -ıku (7)
is used. For this choice of operator S, it can be shown [START_REF] Boubendir | An Analysis of the BEM-FEM Non-Overlapping Domain Decomposition Method for a Scattering Problem[END_REF] that the resulting iteration operator only acts on the part of the spectrum corresponding to the propagating modes while the eigenvalues related to the evanescent modes have unit modulus [START_REF] Boubendir | An Analysis of the BEM-FEM Non-Overlapping Domain Decomposition Method for a Scattering Problem[END_REF]. This directly impacts the convergence (divergence) properties of the resulting iterative scheme. Two families of techniques have been proposed to overcome this problem. First, algorithms based on the optimization of the rate of convergence were introduced by Gander et al. [START_REF] Gander | Optimized Schwarz Methods without Overlap for the Helmholtz Equation[END_REF], where improved local second-order approximations of the DtN map with optimized coefficients are used. In this case, the local transmission operators take the following general form for a generic transmitting boundary Σ:

S oo2 u = -αdiv Σ (A∇ Σ u) + au, (8)
where A and a are a surface tensor (which can be zero, leading to the usual Després approach) and a surface function, respectively, that can be optimized according to the rate of convergence. The operator div Σ in ( 8) is the surface divergence of a tangent vector field on Σ and ∇ Σ is the tangential gradient of a surface field. Second, the Evanescent Modes Damping Algorithm (EMDA) was introduced by Boubendir in [START_REF] Boubendir | An Analysis of the BEM-FEM Non-Overlapping Domain Decomposition Method for a Scattering Problem[END_REF][START_REF] Boubendir | Coupling of a Non-Overlapping Domain Decomposition Method for a Nodal Finite Element Method with a Boundary Element Method[END_REF], with the explicit aim to damp the evanescent modes:

S X u = -ıku + X u (9)
where X is a self-adjoint positive operator. We only consider here the usual case where X is a real-valued positive coefficient.

In this paper we propose a new "square-root" transmission operator [START_REF] Antoine | An Improved Surface Radiation Condition for High-Frequency Acoustics Scattering Problems[END_REF] that takes the following form:

S sq,ε u = -ık 1 + div Σ ( 1 k 2 ε ∇ Σ )u, (10)
where

k ε = k + ıε (11)
is a complexified wavenumber, and where the square-root √ A of an operator A is classically defined through the spectral decomposition of A [START_REF] Taylor | Pseudodifferential Operators[END_REF]. Furthermore, the notation √ z designates the principal determination of the square-root of a complex number z with branch-cut along the negative real axis.

We will show in what follows that

The choice of this square-root operator is motivated by developments of Absorbing Boundary Conditions (ABC) for scattering problems, where the goal is to minimize reflection at the fictitious surface Γ ∞ . In the eighties, the Engquist-Majda [START_REF] Engquist | Absorbing Boundary Conditions for the Numerical Simulation of Waves[END_REF] and Bayliss-Turkel [START_REF] Bayliss | Boundary Conditions for the Numerical Solutions for Elliptic Equations in Exterior Regions[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Condition on Surfaces of Arbitrary Shape[END_REF][START_REF] Bayliss | Radiation Boundary-Conditions for Wave-Like Equations[END_REF] ABCs were developed, which are written as

∂ n u + M u = 0 on Γ ∞ , ( 12
)
where M is a local approximation of the DtN operator that takes a form similar to [START_REF] Antoine | Phase Reduction Models for Improving the Accuracy of the Finite Element Solution of Time-Harmonic Scattering Problems I: General Approach and Low-Order Models[END_REF]. While these conditions are widely used for practical computations and most particularly for engineering purposes, it was recently shown in [START_REF] Antoine | An Improved Surface Radiation Condition for High-Frequency Acoustics Scattering Problems[END_REF] that they also do not model the evanescent and grazing modes, resulting in a loss of accuracy of the method. In order to solve this problem, a high-order local ABC was introduced in [START_REF] Antoine | An Improved Surface Radiation Condition for High-Frequency Acoustics Scattering Problems[END_REF][START_REF] Kerchroud | Numerical Accuracy of a Padé-Type Non-Reflecting Boundary Condition for the Finite Element Solution of Acoustic Scattering Problems at High-Frequency[END_REF], which uses M = S sq,ε in [START_REF] Bendali | A FETI-like Domain Decomposition Method for Coupling FEM and BEM in Large-Size Problems of Acoustic Scattering[END_REF] to model all the scattering modes: propagating, evanescent as well as (in an approximate way) grazing. This ABC can be localized with complex Padé approximants, and the coefficient ε in [START_REF] Bendali | Non-Overlapping Domain Decomposition Method for a Nodal Finite Element Method[END_REF] can then be chosen to minimize spurious reflections at the boundary. Without giving too much details about the derivation of the operator S sq,ε in [START_REF] Bayliss | Radiation Boundary-Conditions for Wave-Like Equations[END_REF], let us briefly explain where it comes from (see [START_REF] Antoine | An Improved Surface Radiation Condition for High-Frequency Acoustics Scattering Problems[END_REF] for more details). To obtain [START_REF] Bayliss | Radiation Boundary-Conditions for Wave-Like Equations[END_REF], we compute for a general surface Σ the principal symbol σ 1 (= -ık 1 -|ξ| 2 /k 2 ) [START_REF] Taylor | Pseudodifferential Operators[END_REF] of the exact DtN operator Λ [START_REF] Antoine | Bayliss-Turkel-like Radiation Condition on Surfaces of Arbitrary Shape[END_REF]. (The covariable ξ = (ξ 1 , ξ 2 ) is roughly the Fourier variable locally to a point of the surface after a local mapping [START_REF] Antoine | Bayliss-Turkel-like Radiation Condition on Surfaces of Arbitrary Shape[END_REF], with

|ξ| 2 = ξ 2 1 + ξ 2 2 .
) If we denote by σ the total symbol [START_REF] Taylor | Pseudodifferential Operators[END_REF] of Λ, it can be proved by using a factorization theorem [START_REF] Antoine | Bayliss-Turkel-like Radiation Condition on Surfaces of Arbitrary Shape[END_REF] at the microlocal level that we have an asymptotic expansion: σ ≈ σ 1 + σ R , where σ R is a remaining symbol of order less than one and therefore has an associated pseudodifferential operator of the same order. Then, we only retain the principal symbol σ 1 and do the approximation: Λ ≈ Op(σ 1 ), where Op(σ 1 ) designates the pseudodifferential operator with symbol σ 1 . In practice, for a general surface, it is difficult to write more precisely the corresponding operator. However, it can be proved [START_REF] Antoine | A Construction of Beam Propagation Methods for Optical Waveguides[END_REF] that, up to lower order symbols, the principal symbol of the square-root operator -ık 1 + div Σ ( 1

k 2 ε ∇ Σ ) is equal
to σ 1 , for ε = 0. Therefore, since we do the approximation modulo a lower order operator, we can approximate Λ by -ık 1 + div Σ ( 1 k 2 ε ∇ Σ ) and get the corresponding square-root transmission boundary condition. The damping parameter ε is then introduced to regularize the square-root operator with the aim to model more correctly the frequencies within the transition zone (called glancing zone in microlocal analysis [START_REF] Antoine | Bayliss-Turkel-like Radiation Condition on Surfaces of Arbitrary Shape[END_REF]). Let us remark that the derivation is similar to Beam Propagation Methods (BPMs) for bent optical waveguides [START_REF] Antoine | A Construction of Beam Propagation Methods for Optical Waveguides[END_REF][START_REF] Lu | Some Techniques for Computing Wave Propagation in Optical Waveguides[END_REF][START_REF] Ho | Improving Beam Propagation Method for TM polarization[END_REF], and that the theoretical factorization theorem used in the preceeding derivation allows to split the outgoing and incoming parts of the wave field through transmission boundary conditions [START_REF] Antoine | Bayliss-Turkel-like Radiation Condition on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A Construction of Beam Propagation Methods for Optical Waveguides[END_REF].

The Després and EMDA transmission operators [START_REF] Antoine | A Construction of Beam Propagation Methods for Optical Waveguides[END_REF] and ( 9) are particular cases of Impedance Boundary Conditions (IBC) [START_REF] Senior | Approximate Boundary Conditions in Electromagnetics[END_REF]. The second order [START_REF] Antoine | Phase Reduction Models for Improving the Accuracy of the Finite Element Solution of Time-Harmonic Scattering Problems I: General Approach and Low-Order Models[END_REF] and the new square-root operator (10) (as well as its Padé localization) fall into the realm of so-called Generalized Impedance Boundary Conditions (GIBC) [START_REF] Senior | Approximate Boundary Conditions in Electromagnetics[END_REF]. In what follows we will thus refer to the transmission conditions related to the Després and EMDA algorithms by IBC(0) and IBC(X ), respectively, and to the square-root transmission condition by GIBC(sq,ε).

Convergence Analysis for a Model Problem

In order to study the convergence of the proposed DDM with square-root transmission operator, we analyze the model problem depicted in Figure 2. This model problem couples two subdomains: a diskshaped bounded subdomain Ω 1 of radius R 0 and an unbounded domain

Ω 0 = R 2 \ Ω 1 : Ω 0 := {x ∈ R 2 , |x| > R 0 }, Ω 1 := {x ∈ R 2 , |x| < R 0 }, (13)
with ∂Ω 0 = ∂Ω 1 = Σ. We study the spectral properties of the iteration operator obtained from the domain decomposition algorithm coupling these two subdomains. Studying the coupling of bounded and unbounded subdomains will allow us to understand the main properties that one could not analyze by considering two bounded (e.g. a square domain divided in two) or two unbounded (e.g. two half-planes) subdomains. The considered model problem essentially contains the main difficulties encountered when solving exterior scattering problems. Note that a closely related problem where the unbounded domain is replaced by an annulus-shaped bounded domain with an ABC on the exterior boundary could also have been selected. However, this would only lead to more complex expressions without providing additional insight.

n 1 Ω 0 R 0 Σ Ω 1 n 0

Convergence Theorem

For the considered model problem, the domain decomposition method with square-root transmission condition consists in solving separately at each iteration the following problems

       ∆u 0 + k 2 u 0 = 0 in Ω 0 , ∂ n0 u 0 + S sq,ε u 0 = g 0 on Σ, lim |x|→∞ |x| 1/2 (∂ |x| u 0 -ıku 0 ) = 0, (14)
and

∆u 1 + k 2 u 1 = 0 in Ω 1 , ∂ n1 u 1 + S sq,ε u 1 = g 1 on Σ, ( 15
)
where g 0 and g 1 are defined as in equation ( 4). The convergence analysis can be developed by studying the spectral properties of the iteration operator A defined by [START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF], where

A := ΠT , (16)
and where the two operators T and Π are defined as [START_REF] Boubendir | An Analysis of the BEM-FEM Non-Overlapping Domain Decomposition Method for a Scattering Problem[END_REF]:

T g (n) := T 1 0 0 T 0   g (n) 1 g (n) 0   and Π := 0 I I 0 , (17)
with

T 1 g (n) 1 := -g (n) 1 + 2S sq,ε u (n+1) 1 , T 0 g (n) 0 := -g (n) 0 + 2S sq,ε u (n+1) 0 . ( 18
)
Since we have a circular geometry, we can analyze the error mode-by-mode by using a Fourier-Hankel series expansion in the polar coordinates system (r, θ). Let us set

u ℓ (r, θ) = +∞ m=-∞ u ℓ,m (r)e ımθ and g ℓ (r, θ) = +∞ m=-∞ g ℓ,m (r)e ımθ , ℓ = 0, 1. (19)
Since {e ımθ } m∈Z defines an orthonormal basis of L 2 (Σ), problems ( 14) and ( 15) lead to solve the decoupled problems

         1 r ∂ r (r∂ r u 0,m ) - m 2 r 2 u 0,m + k 2 u 0,m = 0 for r > R 0 , -∂ r u 0,m + S sq,ε m u 0,m = g 0,m for r = R 0 , lim r→+∞ r 1/2 (∂ r u 0,m -ıku 0,m ) = 0, ( 20
) and    1 r ∂ r (r∂ r u 1,m ) - m 2 r 2 u 1,m + k 2 u 1,m = 0 for r < R 0 , ∂ r u 1,m + S sq,ε m u 1,m = g 1,m for r = R 0 . (21)
In the polar coordinate system the surface divergence and the tangential gradient appearing in (10) reduce to the curvilinear derivative R -1 0 ∂ θ . Hence, for a mode m, we get the analytical expression of the square-root transmission operator:

S sq,ε m = -ık 1 - m 2 k 2 ε R 2 0 := -ık(R m + ıX m ). ( 22
)
The solution of the first and last equation of the exterior problem [START_REF] Després | Domain Decomposition Method and the Helmholtz Problem[END_REF] is given by u 0,m (r) = α m H Writing the modal decomposition of the operators T 0 and T 1 as T 0 = +∞ m=-∞ T 0,m e ımθ and T 1 = +∞ m=-∞ T 1,m e ımθ , one can show using the second equations of ( 20) and ( 21) that:

T ℓ,m = -kZ ℓ,m + S sq,ε m kZ ℓ,m + S sq,ε m = -Z ℓ,m + ı(R m + ıX m ) Z ℓ,m + ı(R m + ıX m ) , ℓ = 0, 1, (23)
with

Z 0,m = H (1) ′ m (kR 0 ) H (1) m (kR 0 ) and Z 1,m = - J ′ m (kR 0 ) J m (kR 0 ) . ( 24
)
Since the convergence properties of the DDM are strongly related to the spectrum of the iteration operator A, we study the spectrum of the modal matrices

A m = 0 T 0,m T 1,m 0 (25)
for each mode, where A = +∞ m=-∞ A m e ımθ . Let us begin by considering that the iteration equation ( 6) is solved using the Jacobi algorithm. Then, according to Theorem 1 below (see [START_REF] Boubendir | An Analysis of the BEM-FEM Non-Overlapping Domain Decomposition Method for a Scattering Problem[END_REF]), convergence occurs if and only if the spectral radius ρ(A) = max m∈Z ρ(A m ) is smaller than one:

Theorem 1 Let s ∈ R be such that (g 0 , g 1 ) ∈ V := H -s (Σ) × H -s (Σ).
The domain decomposition algorithm converges in V if and only if for all m, ρ(A m ) < 1, with ρ(A m ) the largest modulus of the two eigenvalues of the 2 × 2 matrix A m .

In our case, the two eigenvalues of A m are: λ ± m = ± T 0,m T 1,m . Hence, it is sufficient to prove that |T 0,m T 1,m | < 1 for each mode m to guarantee the convergence of the iterative method. Before proving this inequality, let us first establish the following lemma.

Lemma 2 Consider the function F defined by

F(R, X , Z) = -Z + ı(R + ıX ) Z + ı(R + ıX ) , ( 26
)
where the real coefficients R and X are positive (R > 0, X > 0) and where Z is a complex number, with ℜ(Z) := -x and ℑ(Z) := y. Then, under the conditions x > 0 and y > 0, we have

|F(R, X , Z)| < 1. ( 27
)
Proof. By writing Z = -x + ıy, we obtain

|F(R, X , Z)| 2 = (x -X ) 2 + (R -y) 2 (x + X ) 2 + (R + y) 2 < 1, (28) since
(R -y) 2 < (R + y) 2 and (x -X ) 2 ≤ (x + X ) 2 .
To prove the convergence of the algorithm, we first study separately the two coefficients T 0,m and T 1,m , related to the exterior and interior problems respectively. The behavior of the coefficients depends on the choice of the damping parameter ε appearing in the definition of the square-root operator [START_REF] Dular | A General Environment for the Treatment of Discrete Problems and its Application to the Finite Element Method[END_REF] through [START_REF] Bendali | Non-Overlapping Domain Decomposition Method for a Nodal Finite Element Method[END_REF].

Proposition 3 For ε > 0 and any m ∈ Z, we have:

|T 0,m | < 1. ( 29
)
If ε = 0 and m = ±kR 0 , then the inequality (29) also holds.

Proof. For any mode m, the coefficient Z 0,m can be written as

Z 0,m = H (1) ′ m (kR 0 )H (1)
m (kR 0 ) |H (1) m (kR 0 )| 2 , ( 30
)
with ℑ(H

(1) ′ m (kR 0 )H (1)
m (kR 0 )) > 0 and ℜ(H

(1) ′ m (kR 0 )H (1)
m (kR 0 )) < 0 (see [START_REF] Boubendir | An Analysis of the BEM-FEM Non-Overlapping Domain Decomposition Method for a Scattering Problem[END_REF] for the proof). Let us first assume that ε = 0. If one has a propagating mode, i.e., |m| < kR 0 , then S sq,0 m = -ıkR m . For the evanescent modes such that kR 0 < |m|, S sq,0 m = kX m . In these two cases, we deduce that |T 0,m | < 1 by using lemma 2 if there is no mode m such that m = ±kR 0 .

Let us now consider the damped version ε > 0 of the transmission operator S sq,ε m . We first recall that if a and b are two real numbers such that b = 0, then √ a + ıb = p + ıq, where p and q are two real numbers given by

p = 1 √ 2 a 2 + b 2 + a, q = sign(b) √ 2 a 2 + b 2 -a. (31)
Expanding S sq,ε m = -ık(R m + ıX m ) as

S sq,ε m = -ık R 2 0 ((k 2 -ε 2 ) 2 + 4k 2 ε 2 ) -m 2 (k 2 -ε 2 ) + 2ım 2 kε R 2 0 ((k 2 -ε 2 ) 2 + 4k 2 ε 2 ) , (32)
and knowing that k > 0 and ε > 0, we get that R m > 0 and X m > 0. Using [START_REF] Ianculescu | Parallel Iterative Solution for the Helmholtz Equation with Exact Non-Reflecting Boundary Conditions[END_REF] and Lemma 2, this allows us to conclude that |T 0,m | < 1, for any m ∈ Z and ε > 0.

Concerning the interior problem, we have the following result.

Proposition 4 If ε ≥ 0 and |m| > kR 0 , we have

|T 1,m | < 1. ( 33
)
If ε = 0 and |m| ≤ kR 0 , we have

|T 1,m | = 1. ( 34
)
Proof. We observe that Z 1,m is real. Let us first consider the case when ε = 0. If |m| ≤ kR 0 , then X m = 0 and we trivially have |T 1,m | = 1. When |m| > kR 0 , we study the sign of Z 1,m . Since

Z 1,m = - J ′ m (kR 0 ) J m (kR 0 ) = - u ′ 1,m (R 0 ) ku 1,m (R 0 )
, by using the variational formulation of problem ( 21) we obtain

ℜ R 0 u ′ 1,m (R 0 )u 1,m (R 0 ) = R0 0 r|u ′ 1,m | 2 + ( m 2 r 2 -k 2 )|u 1,m | 2 r dr.
Hence, for all |m| > kR 0 , we deduce that

ℜ(R 0 u ′ 1,m (R 0 )u 1,m (R 0 )) > 0. Furthermore, R m = 0 but X m = 0, meaning that |T 1,m | < 1.
Let us now consider the case ε > 0. If |m| > kR 0 the same procedure can be used, in which case

R m > 0 and |T 1,m | < 1.
¿From the preceding propositions we can prove the following result.

Theorem 5 For ε = 0 and if m is not a cut-off mode, i.e. |m| = kR 0 , we have the inequality

|T 0,m T 1,m | < 1. ( 35
)
If ε > 0, then there exists ε max such that for 0 < ε < ε max the inequality (35) is also satisfied.

Proof. Let ε = 0 and let us assume that there is no cut-off mode. Then, the result is a consequence of Propositions 3 and 4. If ε > 0, the inequality ( 35) is true if m corresponds to an evanescent mode. However, we cannot conclude directly for the propagating ones because of the operator T 1 . These propagating modes are in finite number and for ε = 0 we have [START_REF] Kerchroud | Performance Study of Plane Wave Finite Element Methods with a Padé-Type Artificial Boundary Condition in Acoustic Scattering[END_REF]. Therefore, we can then deduce by continuity on ε that there exists ε max such that for 0 < ε < ε max , the inequality ( 35) is satisfied.

The preceding theorem ensures convergence of the Jacobi domain decomposition algorithm with the square-root operator when 0 < ε < ε max . The analysis is not valid for the cut-off mode |m| = kR 0 when ε = 0: indeed, the operator S sq,ε is then not invertible, which makes the domain decomposition algorithm ill-posed. In the actual implementation of the DDM it is thus important to consider ε > 0 in order to guarantee the well-posedness of the DDM. In addition, as we will show, ε > 0 contributes in improving the convergence in the transition zone between the propagating and evanescent modes.

Quasi-Optimality

Since

H (1) m (t) ≈ - 2 m (m -1)! iπt m and J m (t) ≈ - t m 2 m m! for |m| ≫ |t|, (36)
we get, for |m| ≫ kR 0 , Z 0,m ≈ -m kR 0 and Z 1,m ≈ m kR 0 [START_REF] Lu | Some Techniques for Computing Wave Propagation in Optical Waveguides[END_REF] and thus lim This shows that there is a clustering of the eigenvalues for |m| ≫ kR 0 . Obviously, the limit is zero if ε = 0. Therefore, the rate of convergence is optimal for the evanescent modes. This property remains true for ε > 0 since k ≫ ε.

m-→∞ λ ± m = ± ıε -2k -ıε for ε > 0. (38)
To demonstrate the efficiency of the algorithm on the rest of the modes, let us choose k = 6π, R 0 = 1 and a maximal number of modes m max = [10kR 0 ] (where [10kR 0 ] denotes the integer part of 10kR 0 ), and let us compare the new algorithm with the Després and EMDA algorithms. We report in Figure 3(a) the modal spectral radius ρ(A m ) with respect to m for the transmitting boundary conditions IBC(0) (Després), IBC(k/2) (EMDA) and GIBC(sq, 0). As is well-known, IBC(0) has a spectral radius equal to 1 for the evanescent modes, which is improved by IBC(k/2). For the new transmission condition with ε = 0, we clearly observe an optimal convergence rate for the evanescent part, which is coherent with our previous comment. We also observe a significant improvement over the Després and EMDA algorithms on the propagating modes.

The damping parameter ε can be optimized to further improve the spectrum of the iteration operator corresponding to the modes in the transition zone. If we denote by ρ sq,ε the spectral radius of A m when the GIBC(sq,ε) transmission condition is used, we look for ε opt such that

ρ sq,εopt = min ε>0 ρ sq,ε . (39)
As mentioned above, the square-root operator is a suitable approximation of the DtN operator for the propagating and the evanescent modes. By construction, it appears that the maximal error between the exact DtN operator and its square-root approximation S sq,ε is obtained for the modes m satisfying m = ±[kR 0 ] [6]. These specific modes are better known as grazing modes in the context of scattering problems. For such problems it was shown in [START_REF] Antoine | An Improved Surface Radiation Condition for High-Frequency Acoustics Scattering Problems[END_REF] that the optimal value of the damping parameter is ε opt = 0.4k 1/3 H 2/3 , where H is the mean curvature on Σ. The resulting modal spectral radius for GIBC(sq,ε opt ) is shown in Figure 3(a), where the improvement around m = kR 0 is clearly visible. We can see in Figure 3(b) that the spectral radius of the iteration operator is indeed minimum for ε = ε opt . From Figure 3(b) one can also observe that X opt = k/2 (the value chosen for EMDA in [START_REF] Boubendir | An Analysis of the BEM-FEM Non-Overlapping Domain Decomposition Method for a Scattering Problem[END_REF][START_REF] Boubendir | Coupling of a Non-Overlapping Domain Decomposition Method for a Nodal Finite Element Method with a Boundary Element Method[END_REF]) is a fairly good choice to minimize the spectral radius for IBC(X ). In this example, X opt = 3k/4 could also be a good choice.

Krylov Subspace Solvers

Now, let us analyze the solution of (6) by using a Krylov subspace method like GMRES. It is well-known that fast convergence of the GMRES solver is related to the existence of clustering of its eigenvalues. Let us recall that the equation that we solve is given by the operator: (I -A). Therefore, its eigenvalues are: µ ± m = 1 ∓ λ ± m for each mode m ∈ Z. We report in Figure 4 the spectrum of the iteration operator for IBC(0), IBC(k/2) and GIBC(sq,ε opt ). We consider again kR 0 = 6π for a maximal number of modes m max = [10kR 0 ]. For all transmission operators the spectrum lies in the right half-plane, which makes the GMRES converging. However, for IBC(0) many eigenvalues spread out in the complex plane. A slightly better clustering occurs for IBC(k/2), while there is an excellent clustering of the eigenvalues for GIBC(sq, ε opt ). In particular, only a few eigenvalues related to the propagating modes do not cluster but are very close to (1, 0). Furthermore, the eigenvalues related to the evanescent modes seem to cluster at (1, 0). According to [START_REF] Milinazzo | Rational Square-Root Approximations for Parabolic Equation Algorithms[END_REF], we show that the clustering is at point

C * = k + ıε opt k + ıε opt /2 (40)
(in our case C * = 1.0008+0.0282ı). Since the eigenvalues cluster for the evanescent modes, we can expect that in the numerical implementation the GMRES convergence rate will be independent of the density of discretization points per wavelength n λ [START_REF] Antoine | Alternative Integral Equations for the Iterative Solution of Acoustic Scattering Problems[END_REF][START_REF] Antoine | Generalized Combined Field Integral Equations for the Iterative Solution of the Three-Dimensional Helmholtz Equation[END_REF][START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF]. This will be indeed the case as seen in Section 7.

Localization of the Square-Root Operator by Complex Padé Approximants

The square-root operator ( 10) is a nonlocal operator (a pseudo-differential operator of order +1). Therefore, it is impractical in a finite element setting since it would lead to consider full matrices for the transmission boundaries. Fortunately, a localization process of this operator can be efficiently done by using partial differential (local) operators to have a sparse matrix representation. In [START_REF] Kerchroud | Performance Study of Plane Wave Finite Element Methods with a Padé-Type Artificial Boundary Condition in Acoustic Scattering[END_REF][START_REF] Kerchroud | Numerical Accuracy of a Padé-Type Non-Reflecting Boundary Condition for the Finite Element Solution of Acoustic Scattering Problems at High-Frequency[END_REF][START_REF] Antoine | An Improved Surface Radiation Condition for High-Frequency Acoustics Scattering Problems[END_REF], this is realized by a rotating branch-cut approximation of the square-root and then applying complex Padé approximants of order N p [START_REF] Milinazzo | Rational Square-Root Approximations for Parabolic Equation Algorithms[END_REF]:

1 + div Γ ∞ ( 1 k 2 ε ∇ Γ ∞ )u ≈ R α Np (div Γ ∞ ( 1 k 2 ε ∇ Γ ∞ ))u = C 0 u + Np ℓ=1 A ℓ div Γ ∞ (k -2 ε ∇ Γ ∞ )(1 + B ℓ div Γ ∞ (k -2 ε ∇ Γ ∞ )) -1 u, (41)
which corresponds to the complex Padé approximation

√ 1 + z ≈ R α Np (z) = C 0 + Np ℓ=1 A ℓ z 1 + B ℓ z . ( 42
)
Let us mention some similar ideas in the approximation of the square-root with the approach described in [START_REF] Guddati | Continued Fraction Absorbing Boundary Conditions for Convex Polygonal Domains[END_REF]. The complex coefficients C 0 , A ℓ and B ℓ are given by

C 0 = e ı α 2 R Np (e -ıα -1), A ℓ = e -ıα 2 a ℓ (1 + b ℓ (e -ıα -1)) 2 , B ℓ = e -ıα b ℓ 1 + b ℓ (e -ıα -1) . ( 43
)
In [START_REF] Senior | Approximate Boundary Conditions in Electromagnetics[END_REF] α is the angle of rotation, (a ℓ , b ℓ ), ℓ = 1, ..., N p , are the standard real Padé coefficients

a ℓ = 2 2N p + 1 sin 2 ( ℓπ 2N p + 1 ) , b ℓ = cos 2 ( ℓπ 2N p + 1
), [START_REF] Singer | A Perfectly Matched Layer for the Helmholtz Equation in a Semi-Infinite Strip[END_REF] and R Np is the real Padé approximant of order N :

√ 1 + z ≈ R Np (z) = 1 + Np ℓ=1 a ℓ z 1 + b ℓ z . ( 45
)
For a variational representation, the approximation of the Padé-localized square-root transmission operator ( 10)-( 45) is realized by using auxiliary coupled functions [START_REF] Kerchroud | Performance Study of Plane Wave Finite Element Methods with a Padé-Type Artificial Boundary Condition in Acoustic Scattering[END_REF][START_REF] Kerchroud | Numerical Accuracy of a Padé-Type Non-Reflecting Boundary Condition for the Finite Element Solution of Acoustic Scattering Problems at High-Frequency[END_REF][START_REF] Antoine | An Improved Surface Radiation Condition for High-Frequency Acoustics Scattering Problems[END_REF]]

S Np,α,ε u = -ık(C 0 u + Np ℓ=1 A ℓ div Γ ∞ ( 1 k 2 ε ∇ Γ ∞ ϕ ℓ )) on Γ ∞ , ( 46
)
where the functions ϕ ℓ , ℓ = 1, .., N p , are defined on Γ ∞ as the solutions of the following surface PDEs:

(1 + B ℓ div Γ ∞ ( 1 k 2 ε ∇ Γ ∞ ))ϕ ℓ = u. ( 47
)
The resulting transmitting boundary condition is again a Generalized Impedance Boundary Condition, and is denoted by GIBC(N p , α, ε) for the Padé approximation with N p auxiliary functions, for an angle of rotation α and a damping parameter ε. The lowest-order approximation would be S Np,α,ε = -ıkI which also corresponds to the original approximation S 0 of S, i.e., IBC(0).

Let us analyze on the model problem of Section 4 the effect of the Padé approximation on the rate of convergence. For a mode m, the transmission operator acts as

S Np,α,ε m = -ıkR α Np (- m 2 k 2 ε R 2 0 ). ( 48
)
In the following the angle of rotation α is always taken equal to π/4, which was found to be the optimal choice through numerical experiments. We report in Figure 5 the spectrum of GIBC(sq, ε opt ), GIBC(4, π/4, ε opt ) and GIBC(8, π/4, ε opt ). As previously noticed, there is an excellent clustering of the eigenvalues for GIBC(sq, ε opt ). As expected, the larger N p , the better the approximation of the spectrum of the square-root. Moreover, N p allows to adjust the spectrum accuracy for large modes m (evanescent modes). Hence, it can be conjectured that N p must be taken with respect to the density of discretization points per wavelength n λ . We will see in the numerical simulations (see Section 7) that relatively small values of N p (N p = 2, 4, 8) give optimal convergence results.

Finite Element Implementation

Let us now describe the implementation of the domain decomposition algorithm with Padé-type transmission condition in a finite element context. Recall first that the iterative method consists in solving problems (3) and computing the transmitted quantities (4). We restrict this description only to a problem posed in a domain, noted here Ω, with no exterior interface nor interface connected to the scatterer. This simply implies that Γ i = ∅, Γ ∞ i = ∅ and ∂Ω is an artificial interface. The extention to the other cases is direct. For the sake of clarity, we abbreviate the notations u

(n+1) i , g (n) ij and g (n+1) ij
, used in problem (3) and equation ( 4), by u, g in and g out respectively. An iteration then leads to solving first the interior boundary value problem in terms of u with coupling auxiliary functions {ϕ ℓ } ℓ=1,...,Np according to the use of Padé approximation of order N p :

             ∆u + k 2 u = 0 in Ω, ∂ n u -ıkC 0 u -ık Np ℓ=1 A ℓ div ∂Ω ( 1 k 2 ε ∇ ∂Ω ϕ ℓ ) = g in on ∂Ω, -u + B ℓ div ∂Ω ( 1 k 2 ε ∇ ∂Ω ϕ ℓ ) + ϕ ℓ = 0, ℓ = 1, ..., N p , on ∂Ω, (49)
where C 0 , A ℓ and B ℓ are given by [START_REF] Senior | Approximate Boundary Conditions in Electromagnetics[END_REF]. The auxiliary functions {ϕ ℓ } ℓ=1,...,Np are introduced to deal with the use of the Padé approximation of order N p . Considering now test-functions v ∈ H 1 (Ω) and v ℓ ∈ H 1 (∂Ω), ℓ = 1, ..., N p , we get the coupled variational formulation

                         Ωi (∇u • ∇v -k 2 uv) dΩ -ıkC 0 ∂Ω uvd∂Ω +ık Np ℓ=1 A ℓ ∂Ω 1 k 2 ε ∇ ∂Ω ϕ ℓ • ∇ ∂Ω v d∂Ω = - ∂Ω g in v d∂Ω, - ∂Ω uv ℓ d∂Ω -B ℓ ∂Ω 1 k 2 ε ∇ ∂Ω ϕ ℓ • ∇ ∂Ω v ℓ d∂Ω + ∂Ω ϕ ℓ v ℓ d∂Ω = 0, ℓ = 1, ..., N p . (50)
Following the solution of (50), we update the boundary data through the relation

g out = -g in -2ıkC 0 u -2ık Np ℓ=1 A ℓ div ∂Ω ( 1 k 2 ε ∇ ∂Ω ϕ ℓ ), on ∂Ω. (51)
Let us now consider a covering Ω h of Ω using n t tetrahedral (or triangular in two-dimensions) finite elements with n v vertices. Parameter h is the usual finite element notation for the maximal side length of the tetrahedrals. All the notations are extended with h as subscript for the discrete version of the domains as well as unknowns. We choose here linear finite element approximations but the extension to higher-order elements is straightforward. For scattering problems, it is common to introduce the density of discretization points per wavelength: n λ = λ/h. The local density of discretization points is therefore: n λ = λ/h. Let us denote by S Ω h and M Ω h respectively the stiffness and mass matrices for linear elements associated with the domains Ω h . The matrices have size n v × n v . Furthermore, we introduce S ∂Ω h and M ∂Ω h as the respective stiffness and mass matrices related to the transmitting surface ∂Ω h . If these correspond to a generalized stiffness matrix for a surface function β, then it is quoted S ∂Ω h β . All these matrices have a size n ∂ t × n ∂ t . Let us denote by u ∈ C nv the local unknown vector and ϕ ℓ ∈ C n ∂ t the surface unknown auxiliary vectors obtained with linear finite elements. The discrete test-vectors and right hand side are also bold typed. Then, the discretization of the variational problem (49) leads to the solution of the coupled linear system

       (S Ω h -k 2 M Ω h -ıkC 0 M ∂Ω h )u + ık Np ℓ=1 A ℓ S ∂Ω h k -2 ε ϕ ℓ = -M ∂Ω h g in -M ∂Ω h u -(B ℓ S ∂Ω h k -2 ε -M ∂Ω h )ϕ ℓ = 0, ℓ = 1, ..., N p . (52)
The size of this linear system is (

n v + (1 + N p )n ∂ t ) × (n v + (1 + N p )n ∂ t )
. Each augmented system is solved directly by a LU factorization method since each matrix has a moderate size. The boundary data update can be done by vector addition using relation (4).

Numerical Results

This section is devoted to numerical simulations validating the new domain decomposition algorithm using the Padé-localized square-root transmission condition. We use the finite element method with linear (P1) basis functions implemented in GetDP/Gmsh [START_REF] Dular | A General Environment for the Treatment of Discrete Problems and its Application to the Finite Element Method[END_REF][START_REF] Geuzaine | Gmsh: a Three-Dimensional Finite Element Mesh Generator with Built-in Pre-and Post-Processing Facilities[END_REF] to approximate the local problems (3), as explained in Section 6. The convergence criterion for all the presented examples is identical: the iterations are stopped when the initial residual has decreased by a factor of 10 -6 .

Two-Dimensional Examples

The two-dimensional examples concern the scattering of a plane wave by a unit sound-soft circular cylinder. The second-order Bayliss-Turkel [START_REF] Bayliss | Boundary Conditions for the Numerical Solutions for Elliptic Equations in Exterior Regions[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Condition on Surfaces of Arbitrary Shape[END_REF][START_REF] Bayliss | Radiation Boundary-Conditions for Wave-Like Equations[END_REF] artificial boundary condition is set on a fictitious N wavenumber for N dom = 8 (right) when using the Jacobi or GMRES algorithm, for the "circle-concentric" decomposition.

circular boundary with radius 4 (see Figure 6).

Jacobi vs. GMRES

We begin by analyzing the behavior of the iterative methods (successive approximations, i.e. Jacobi, or GMRES [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]) when the domain is partitioned into N dom concentric subdomains. We call this partitioning the "circle-concentric" decomposition (see Figure 6(a)). The "size" L i of the subdomain Ω i , defined as the difference between its exterior and interior radius, is equal to 3/N dom . Le us start by considering a wavenumber k = π and a density of discretization points per wavelength n λ = 12. We report in Table 1 (left) the number of iterations required using Jacobi and GMRES. Concerning the complex Padé approximation, we fix N p = 8 auxiliary equations and θ = π/4. The damping parameter ε is optimized according to ε opt = 0.6k 1/3 R -2/3 j (where R j is the radius of the j-th interface). Let us precise that the "size" of the overall domain is L = 3 = 1.5λ. This means that in the case of N dom = 10 (resp. N dom = 15) subdomains, the size of one subdomain is about 0.15λ (resp. about 0.1λ) which is very small. For these cases, the EMDA algorithm solved by the Jacobi method diverges while our approach converges. As it can be seen, the algorithm scales according to N dom , and as expected, the GMRES outperforms the Jacobi method.

Consider now a partition of the initial domain into N dom = 8. The discretization is again fixed by n λ = 12. The same parameters (N p = 8, θ = π/4 and ε opt = 0.6k 1/3 R -2/3 j ) are kept to approximate the transmission operator given by [START_REF] Roux | Optimization of Interface Operator Based on Algebraic Approach[END_REF]. We report in Table 1 (right) the number of iterations for converging according to the wavenumber k. As we can see, the convergence is quasi independent of the wavenumber k. This is known to be a very difficult goal to attain when designing iterative schemes for scattering problems. In particular, this means that our Padé-type transmission condition with GMRES DDM is well-suited for high frequency scattering problems-we will analyze this behaviour further in Section 7.1.2. ¿From these first numerical results we clearly see that the GMRES algorithm outperforms the successive approximations (Jacobi), especially when the number of subdomains is large. In what follows we will thus only report results for the GMRES algorithm. Let us consider N dom = 5 subdomains. Figure 7(a) shows the number of iterations to reach convergence with respect to the wavenumber k, for two densities of discretization points per wavelength n λ . For EMDA with IBC(k/2), the number of iterations increases with respect to k and n λ . Concerning GIBC(N p , π/4, ε opt ), we can see that the convergence rate is independent of both the wavenumber and density of discretization points per wavelength. Furthermore, in this case, N p = 2 provides already an optimal convergence compared with N p = 8. Let us recall that one of the important points here is that we use complex approximants. Using real Padé approximants would lead to much higher values of N p , hence penalizing the computational cost of the overall method. Let us also mention that increased values of n λ are not only necessary for accuracy purposes but are also imposed to avoid pollution effects in the calculation of the solution at high wavenumbers [START_REF] Ianculescu | Parallel Iterative Solution for the Helmholtz Equation with Exact Non-Reflecting Boundary Conditions[END_REF][START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF][START_REF] Ihlenburg | Dispersion Analysis and Error Estimation of Galerkin Finite Element Methods for the Helmholtz Equation[END_REF][START_REF] Antoine | Phase Reduction Models for Improving the Accuracy of the Finite Element Solution of Time-Harmonic Scattering Problems I: General Approach and Low-Order Models[END_REF]. Figure 7(b) gives the number of iterations with respect to the density of discretization points per wavelength n λ , for two wavenumbers k. A convergence independent of the mesh size can be achieved provided that N p is sufficiently large, although a very small value (N p = 2) already provides a quasi-optimal result. In Figure 7(c) we report the number of iterations of the GMRES DDM with respect to the number of subdomains N dom . We directly see that all the methods scale, for IBC as well as GIBC transmission conditions. Again, for GIBC(N p ,π/4,ε opt ), the scaling does not depend on k (and N p ), while this is not the case for IBC(k/2).

IBC(k/2) k = 2π IBC(k/2) k = 6π GIBC(2,π/4,ǫopt) k = 2π GIBC(2,π/4,ǫopt) k = 6π GIBC(8,π/4,ǫopt) k = 2π GIBC(8,π/4,ǫopt) k = 6π
As a partial conclusion, we can state that using the GMRES algorithm with GIBC(N p ,π/4,ε opt ) in the DDM leads to a quasi-optimal convergence: the convergence rate is independent of both the wavenumber k and the density n λ . Furthermore, the method scales with the number N dom of subdomains. These conclusions are confirmed by analyzing the behavior of our algorithms for the "circle-pie" decomposition (see Figure 6(b)). Since we have straight interfaces we cannot consider our optimized value of ε opt anymore. However, numerical simulations show that ε opt = k/4 is a suitable choice. In figures 8(a)-8(c), we represent the same three sets of curves as for the "circle-concentric" decomposition. We get a convergence independent of the wave number k and density of discretization n λ (by taking N p large enough but relatively small). Again, the algorithm scales with respect to N dom . More generally, the convergence results are even better for the "circle-pie" than for the "circle-concentric" decomposition.

Let us remark here that in the context of the Krylov iterative solution of integral equations for acoustic scattering, quasi optimal convergence independent of the wavenumber and density of discretization points per wavelength was obtained in [START_REF] Antoine | Alternative Integral Equations for the Iterative Solution of Acoustic Scattering Problems[END_REF][START_REF] Antoine | Generalized Combined Field Integral Equations for the Iterative Solution of the Three-Dimensional Helmholtz Equation[END_REF][START_REF] Antoine | Integral Equations and Iterative Schemes for Acoustic Scattering Problems[END_REF] by considering Generalized Combined Field Integral Equations including the square-root operator.

Computational Cost

The overall computation time of the DDM depends on the number of iterations and the cost of solving (52) for each subdomain. As explained in Section 6, we use a direct sparse linear solver, and the main cost in solving (52) is thus the initial computation of the LU factorization of the augmented matrices of size N = (n v + (1 + N p )n ∂ t ) × (n v + (1 + N p )n ∂ t ). We report in Table 2 the number of unknowns N and the factorization time for one subdomain in the case of the "circle-pie" decomposition, with N dom = 5, k = 8π and n λ = 10. The factorization is performed using MUMPS [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF] on a 2.66 GHz MacBook Pro laptop computer. For the typical values of N p required to obtain quasi-optimal convergence (N p = 2 or N p = 8), we see that the computational overhead is small (less than 10%).

Three-Dimensional Example

We now test our algorithm on a more realistic three-dimensional example: the scattering by a submarine. The geometry is represented in Figure 9. For comparison purposes, we choose a unit length submarine. The incident wave is a plane wave with α = (cos π/8, 0, sin π/8) T . We consider two high wave numbers, k = 14π and k = 28π, a density of discretization points per wavelength n λ = 10 and we fix N dom = 5 subdomains. The scattered field computed with the GMRES DDM and Padé-type transmission boundary condition is represented on Figure 9. Figure 10 reports the residual decay for different transmission conditions and parameter values. Since we have a plane interface here, we cannot consider our optimized value of ε opt . As for the "circle-pie" decomposition analyzed in Section 7.1, numerical simulations show that ε opt = k/4 is a suitable choice. We see that GIBC(N p , π/4, ε opt ) leads to a very good convergence rate of the GMRES, and is clearly better than EMDA (with IBC(k/2)). For this three-dimensional case, an optimal value of N p to get the quasi-optimality is N p = 8.

Conclusion

We have proposed a quasi-optimal non overlapping domain decomposition algorithm for the Helmholtz equation. It is based on a suitable approximation of the DtN operator. A convergence analysis for a model problem has been developed, showing the main features of our algorithm. Several numerical tests in both two-and three dimensions validates the numerical quasi-optimality of the proposed algorithm. Extension to Maxwell's equations is currently being developed.

Figure 1 :

 1 Figure 1: Example of a two-dimensional non-overlapping domain decomposition method.

Figure 2 :

 2 Figure 2: Model problem with two subdomains and a circular interface.

( 1 )

 1 m (kr), where H (1) m denotes the Hankel function of the first kind and order m. For the interior problem (21) the solution of the first equation reads u (1) m (r) = β m J m (kr), where J m is the Bessel function of order m.

  Spectral radius of the modal iteration operator Am vs. the mode m.

  Spectral radius of the iteration operator A vs. the damping parameter ε.

Figure 3 :

 3 Figure 3: Spectral radius of the iteration operators and optimal choice of the damping coefficient ε.

Figure 4 :

 4 Figure 4: Eigenvalues distribution in the complex plane for (I -A) and different transmission operators.

Figure 5 :

 5 Figure 5: Eigenvalue distribution in the complex plane for the exact and Padé-localized square-root transmission operator of order 4 (left) and 8 (right).

Figure 6 :

 6 Figure 6: Two-dimensional test cases: reconstruction of the scattered field on the global domain after a DDM computation with k = 4π and N dom = 5.

7. 1 . 2

 12 Influence of the Wavenumber, the Discretization Density and the Number of Subdomains

  , Ndom = 5 IBC(k/2) nλ = 10 IBC(k/2) nλ = 20 GIBC(2,π/4,ǫopt) nλ = 10 GIBC(2,π/4,ǫopt) nλ = 20 GIBC(8,π/4,ǫopt) nλ = 10 GIBC(8,π/4,ǫopt) nλ = Number of iterations vs. wavenumber nλ GMRES iterations circle-concentric, Ndom = 5 IBC(k/2) k = 2π IBC(k/2) k = 6π GIBC(2,π/4,ǫopt) k = 2π GIBC(2,π/4,ǫopt) k = 6π GIBC(8,π/4,ǫopt) k = 2π GIBC(8,π/4,ǫopt) k = Number of iterations vs. mesh density Ndom GMRES iterations circle-concentric, nλ = 10

  Number of iterations vs. number of subdomains

Figure 7 :

 7 Figure 7: Convergence for the "circle-concentric" decomposition.

Figure 8 :

 8 Figure 8: Convergence for the "circle-pie" decomposition.

Figure 9 :

 9 Figure 9: Submarine problem with 5 subdomains: iso-surfaces of the real part of the scattered field for k = 14π.

Table 1 :

 1 Number of iterations vs. number of subdomains for k = π (left) and Number of iterations vs.

	dom Jacobi GMRES	k	Jacobi GMRES
	2	12	9	π	102	41
	5	45	24	2 π	62	36
	10	152	52	3 π	52	38
	15	+500	74	4 π	64	38

Table 2 :

 2 Number of degrees of freedom and computational time for the initial factorization for one subdomain of the "circle-pie" decomposition, with N dom = 5, k = 8π and n λ = 10.

		Number of Unknowns Factorization Time (s)
	IBC(k/2)	20812	0.37
	GIBC(1, π/4, k/4)	21053	0.37
	GIBC(2, π/4, k/4)	21294	0.38
	GIBC(4, π/4, k/4)	21776	0.39
	GIBC(8, π/4, k/4)	22740	0.40
	GIBC(16, π/4, k/4)	24668	0.45
	GIBC(32, π/4, k/4)	28524	0.53

† Institut Elie Cartan Nancy (IECN), Nancy University, INRIA Corida Team, B.P. 239, F-54506 Vandoeuvre-lès-Nancy Cedex, France. ‡ University of Liège, Department of Electrical Engineering and Computer Science, Montefiore Institute B28, B-4000 Liège, Belgium.

The Dirichlet boundary condition models a sound-soft obstacle; Neumann (sound-hard) or Fourier-Robin (impedance) boundary conditions may also be set, which would not fundamentally change the rest of the paper.

the nonlocal operator S sq,ε can be accurately localized using complex Padé approximants, and suitably combined with finite element methods;

the convergence of the resulting DDM is quasi-optimal: the rate of convergence corresponding to the evanescent modes is zero, and significantly improved for the remaining modes. This results in an effective solution where the iterative procedure is quasi independent of both the wavenumber and the level of accuracy of the spatial discretization (mesh refinement).

Acknowledgments Y. Boubendir gratefully acknowledges support from NSF through grant No. DMS-1016405. X. Antoine gratefully acknowledges support from the Agence Nationale pour la Recherche (Ref : ANR-09-BLAN-0057-01) (http://microwave.math.cnrs.fr/) and the Fondation de Recherche pour l'Aéronautique et l'Espace (IPPON Project). C. Geuzaine gratefully acknowledges support from the Belgian Science Policy (IAP P6/21), Belgian French Community (ARC 09/14-02) and Walloon Region (WIST3 No 1017086 "ONELAB").

Residual

submarine, Ndom = 5, nλ = 10


Elisabeth Lahalle

Gilles Fleury

Rawad Zgheib

New non-uniform transmission and ADPCM coding system for improving both signal to noise ratio and bit rate

Keywords: Adaptive transmission, adaptive reconstruction, audio-signal compression

Here we address the problem of adaptive digitaltransmission systems. New systems based on a nonuniform transmission (NUT) principle are proposed, utilizing a recently proposed algorithm for adaptive identification and reconstruction of AR processes subject to missing data. We propose a new adaptive sampling (nonuniform transmission) method combined with the adaptive reconstruction algorithm. A new NUT-ADPCM coding-decoding system is designed. The proposed system is demonstrated for audio-signal compression and compared to the ADPCM G.726 standard. The new system yields improvements in both signal-to-noise ratio and average bit rate.

I. INTRODUCTION

I N the design of digital transmission systems, it is important to find a good tradeoff between a low bit rate and a high signal-to quantization-noise ratio [START_REF] Rabiner | Digital processing of speech signals[END_REF], [START_REF] Wade | Codage et traitement du signal[END_REF]. In a previous study [START_REF] Mirsaidi | Reducing quantization error using prediction/non uniform transmission[END_REF], a new concept for the design of digital coding systems was introduced, based on nonuniform transmission of signal samples. The idea is to avoid transmitting a sample if it can be efficiently predicted. The use of the nonuniform transmission principle in an ADPCM coding system was first suggested in [START_REF] Mirsaidi | An ADPCM-like system based on non uniform signal transmission[END_REF]. An LMS-like algorithm [START_REF] Mirsaidi | LMS Like AR modeling in the case of missing observations[END_REF] was suggested for the prediction of samples that were not transmitted. However, this algorithm converges toward biased estimations of the model parameters and does not use an optimal predictor in the leastmean-squares sense [START_REF] Zgheib | New fast recursive algorithms for simultaneous reconstruction and identification of AR processes with missing observations[END_REF].

Recently, we proposed two new adaptive algorithms for the identification and reconstruction of nonstationary AR processes subject to missing data using a Kalman filter for the prediction. The first algorithm is based on a pseudolinear RLS algorithm for the identification [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF], and the second [START_REF] Zgheib | Lattice algorithm for adaptive stable identification and robust reconstruction of non stationary AR processes with missing observations[END_REF] is based on a lattice filter structure, thus guaranteeing stability and robustness [START_REF] Makhoul | Linear prediction: A tutorial review[END_REF], [START_REF] Makhoul | Adaptive lattice analysis of speech[END_REF]. Both algorithms are fast and offer an optimal reconstruction error in the least-meansquares sense. They showed good performance in terms of quadratic reconstruction error when applied to speech-signal reconstruction with a high probability of missing samples. Indeed, nonstationary AR processes can model a large number of signals in practical situations such as speech signals [START_REF] Rabiner | Digital processing of speech signals[END_REF], [START_REF] Makhoul | Adaptive lattice analysis of speech[END_REF]. We here propose to use these algorithms in a nonuniform transmission (NUT) system, previously introduced in [START_REF] Mirsaidi | Reducing quantization error using prediction/non uniform transmission[END_REF], to E. Lahalle, G. Fleury and R. Zgheib are with E3S -SUPELEC Systems Sciences / Dept. of Signal Processing and Electronic Systems, Gif-sur-Yvette, France e-mail: (firstname.lastname@supelec.fr), tel: +33(0)1 69 85 14 27, fax:+33(0)1 69 85 14 29. improve its performance. We then combine them with adaptive reconstruction algorithms to design a new NUT-ADPCM coding-decoding system.

In the following, we begin by recalling the nonuniform transmission idea introduced in [START_REF] Mirsaidi | Reducing quantization error using prediction/non uniform transmission[END_REF], then the adaptive algorithm for the reconstruction of AR processes with missing observations [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF]. The performance of NUT transmission systems using [START_REF] Mirsaidi | LMS Like AR modeling in the case of missing observations[END_REF] and [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF] identification algorithms are compared. Finally, a new NUT-ADPCM system is designed. It is compared in the last section to the ADPCM [1] coding system through simulations of speech signals.

A. NUT Coding 1) Adaptive sampling:

In a simplified digital-transmission system, assuming that the channel is perfect, the reconstruction error is equal to the quantization error. Mirsaidi et al. [START_REF] Mirsaidi | Reducing quantization error using prediction/non uniform transmission[END_REF] proposed a novel nonuniform transmission system that reduced quantization error by using a parametric signal-modeling approach for the prediction. In the receiver, a sample at time n is predicted using the estimated model parameters at time n-1 and the available quantized samples. The key idea of the system they proposed is that at time n an efficient prediction may provide a smaller error than a quantization. In this case, it is better to replace the transmitted quantized sample with its prediction in the receiver. Thus, using an efficient prediction method for signals subject to missing samples, the number of transmitted samples is considerably reduced and the reconstruction performance is also improved.

2) Kalman Pseudo-linear RLS Algorithm: Let {x n } be an AR process of order L with parameters {a k }, and { n } the corresponding innovation process of variance σ 2 . The loss process is modeled by an i.i.d. binary random variable {c n }, where c n = 1 if x n is available; otherwise, c n = 0. Let {z n } be the reconstruction of the process {x n }. If x n is available z n = x n ; otherwise, z n = xn , the prediction of x n . To identify the AR process subject to missing data in real time, we use the pseudolinear RLS algorithm [START_REF] Sanchis | Recursive identification under scarce measurements-convergence analysis[END_REF]. In the case of missing observations, the regression vector Ψ n = [x n-1 . . . x n-L ] cannot be constructed with only available data. Missing data are replaced by their predictions, i.e., Ψn = [z n-1 . . . z n-L ] . Thus, Ψn depends on the available parameters. We show [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF] that the least-squares estimate of the parameters is unbiased if the predictor used is optimal in the least-mean-squares sense. A Kalman filter is used for recursive optimal prediction. Let x n = [x n . . . x n-L+1 ] be the state vector. The predicted and filtered estimates are denoted by xn+1|n and xn+1|n+1 , respectively. P n+1|n and P n+1|n+1 are, respectively, the prior and posterior prediction-error covariance matrices. K n+1 is the Kalman filter gain. The Kalman pseudolinear RLS algorithm is summarized, at time n+1, as follows :

A n+1 =      ân,1 . . . . . . ân,L 1 0 0 . . . . . . 0 1 0      (1)
P n+1|n = AP n|n A + R , ( 2
)
xn+1|n = Ax n|n , ( 3
)
Ψn = xn|n , (4)
If x n+1 is available, i.e. c n+1 = 1, the predicted state vector xn+1|n is filtered by the Kalman filter (( 5),( 6),( 7)) and the AR parameters are updated using the RLS-like algorithm (( 8),( 9),( 10)):

K n+1 = P n+1|n c n+1 (c n+1 P n+1|n c n+1 ) -1 , ( 5
)
P n+1|n+1 = (I d -K n+1 c n+1 )P n+1|n , (6)
xn+1|n+1 = xn+1|n + K n+1 (y n+1 -c n+1 xn+1|n ), (7)
γ n+1 = c n G n Ψn+1 λ + Ψ n+1 G n Ψn+1 , ( 8
)
ân+1 = ân + γ n (y n+1 -c n+1 xn+1|n ), (9)
G n+1 = 1 λ (I d -γ n+1 Ψ n+1 )G n (10)
If x n+1 is missing, i.e. c n+1 = 0, the predicted state vector is not filtered by the Kalman filter and the AR parameters are not updated using the RLS-like algorithm :

K n+1 = 0, (11)
P n+1|n+1 = P n+1|n , ( 12
) xn+1|n+1 = xn+1|n , ( 13
)
γ n+1 = 0, (14)
ân+1 = ân , ( 15
)
G n+1 = 1 λ G n (16)
3) Comparison: In this section, the periodic transmission system, the nonuniform transmission system using the LMSlike algorithm proposed in [START_REF] Mirsaidi | LMS Like AR modeling in the case of missing observations[END_REF], and the Kalman pseudolinear RLS Algorithm are compared; the methods are referred to hereafter as Methods 1, 2 and 3, respectively. The test signal used was"Mary had a little lamb, its fleece was white as snow" sampled at F s = 8kHz. The forgetting factor used in the pseudolinear RLS algorithm [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF] is λ = 0.999 and the stepsize used in the LMS-like algorithm [START_REF] Mirsaidi | LMS Like AR modeling in the case of missing observations[END_REF] is µ = 7.10 -5 . In both parametric methods, the speech signal is modeled by an AR process of order 10. The three methods are evaluated in terms of mean bit rate and the signal-to-reconstruction error ratio given by:

(SN R) = 10log 10 N i=1 x 2 i N i=1 (x i -x i ) 2 .
In the case of a transmitted sample, the reconstruction error is the same as the quantization error; in the case of a missing sample, it is the same as the prediction error. For the periodic transmission system, the bit rate is given by r = BF s , where B is the number of bits per sample. For other two methods, a onebit flag is transmitted for each nontransmitted sample. We consider the average bit rate, given by r = (T B +1)F s , where T stands for the proportion of samples transmitted. For comparison purposes, Table I shows the results of the three methods. They are compared with respect to the signal-to-noise ratio (SNR) and the bit-rate performance of the G.726 ADPCM codec at 32 kbits. For almost the same SNR, Method 2 shows a higher average bit rate than for periodic transmission with B = 5, which is of course due to the onebit flag. In contrast, Method 3, using an optimal predictor [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF], outperforms all other methods in both SNR and average bit rate. Consequently, a new nonuniform transmission ADPCM coding system using an optimal predictor [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF] [12] is proposed below.

B. NUT-ADPCM Coding

The ADPCM codec, defined by the ITU-T Recommendations G.726 [1], is based on the differential coding principle (DPCM). Because the difference between the signal and its prediction has a lower variance than the original samples, it can be accurately quantized with fewer bits than the original samples. The performance of the codec is aided by using adaptive prediction and quantization, so that the predictor and difference quantizer adapt to the changing characteristics of the signal being coded. The adaptive predictor used is based on an ARMA(2,6) model of the signal.

1) system design: When the prediction method used is efficient, the prediction error of some samples might be of a magnitude less than the minimal quantization step provided by the adaptive quantizer. We propose avoiding transmitting the quantized prediction errors in that case and so expect to reduce the number of transmitted bits without sacrificing reconstruction quality. Therefore, the adaptive predictor used in the transmission and reception must solve the problem of the reconstruction of signals with missing samples. In addition, the transmission of a one-bit flag to indicate the transmission mode is required. The proposed system differs from the ADPCM codec mainly in the adaptive predictor used. Additionally, it contains a transmission-decision element for the quantized prediction error and for the one-bit flag.

Adaptive predictor: We then use the Kalman Recursive Least-Squares Lattice algorithm (KRLSL) [START_REF] Zgheib | Lattice algorithm for adaptive stable identification and robust reconstruction of non stationary AR processes with missing observations[END_REF] algorithm as an adaptive predictor in the NUT-ADPCM coding system instead of the Kalman RLS [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF] to ensure the stability of the system. The signal is modeled by a nonstationary AR(L) process of order 10 instead of an ARMA [START_REF] Makhoul | Linear prediction: A tutorial review[END_REF][START_REF] Mirsaidi | An ADPCM-like system based on non uniform signal transmission[END_REF]. Recall here that the KRLSL algorithm is based on an adaptive Burg algorithm for the identification of the AR parameters, thus guaranteeing at each time point the stability of the identified filter [START_REF] Zgheib | Lattice algorithm for adaptive stable identification and robust reconstruction of non stationary AR processes with missing observations[END_REF].

Additionally, the KRLSL algorithm uses a Kalman filter for the prediction of the samples. When a sample is available, the state, consisting of the last L samples, is updated by the Kalman filter proportionally to the prediction error of the sample [START_REF] Zgheib | Lattice algorithm for adaptive stable identification and robust reconstruction of non stationary AR processes with missing observations[END_REF]. However, in the case of NUT-coding systems, the receiver has only quantized versions of the prediction errors. Therefore, even if the true prediction error of a sample is known at the transmitter, the quantized prediction error is used with the Kalman filter to update the state. This reproduces the same sample prediction as is done in the receiver. Hence, when a quantized prediction error is received, the sample reconstruction is computed as in equation xn|n = xn|n-1 +e Q n,P . The sample thus reconstructed is used to initialize the latticeidentification algorithm ((17)).

In the KRLSL algorithm, the recursive equations of the RLSL algorithm ((17),( 18),( 19),( 20),( 21),( 22),( 23)) are used instead of the equations of the RLS algorithm ((8),( 9),( 10)) at each available sample, x n+1 . They are used to estimate the reflection coefficients k(l) t at each time t since the last available sample. The AR parameters, at time n + 1, are then deduced from the reflection coefficients k(l) n+1 using the Durbin-Levinson recursions.

f (0) t = b(0) t = xt|n+1 , , (17)
k(0) t = 1, (18)
Ĉ(l) t = λ Ĉ(l) t-1 -2 f (l-1) t b(l-1) t-1 , (19)
D(l) t = λ D(l) t-1 - f (l-1)2 t b(l-1)2 t-1 , (20)
k(l) t = - Ĉ(l) t D(l) t , (21)
f (l) t = f (l-1) t - k(l) t b(l-1) t-1 , (22) b(l) t = b(l-1) t - k(l) t f (l-1) t-1 (23)
Transmission decision: Once the sample is predicted via the Kalman filter, the prediction error is quantized using the same adaptive quantizer as described in the recommendations of the G.726 specifications [1]. If the sample is perfectly predicted (i.e., the prediction error is quantized to zero), the code corresponding to the quantized prediction error is not transmitted. A one-bit flag must also be transmitted to the receiver. In the transmitter, the adaptive predictor must be informed about the transmission decision of the sample code to execute the appropriate steps, which are also executed in the receiver. In this scheme, to reduce the number of bits transmitted, the number of codes transmitted should compensate for the extra bit flag. Thus, let N be the size of a signal coded using an ADPCM codec at B bits. The total number of bits transmitted is then equal to B ADP CM = N B bits. The nonuniform transmission ADPCM system also uses B bits (the same number as the ADPCM codec to which it is compared) to code a quantized prediction error. The total number of bits transmitted by the NUT-ADPCM system is then equal to B N U T = pN B + N f , where p is the ratio of transmitted samples and N f is the number of the transmitted flags. A flag D n = 1 is transmitted at each new transmission following a sequence of non-transmitted samples. A flag D n = 0 is transmitted when the transmission stops. The condition for a profitable use of the NUT-ADPCM system is thus pN B + N f < N B. A simplified diagram of the NUT-ADPCM codec described above is presented in Figure 1. 2) Simulations: The nonuniform transmission ADPCM system described above is compared here to the ADPCM through simulations on the same speech signal as above. Different values of the number of bits are used to code the prediction error B. The different methods are coded using MATLAB and executed on a Pentium 4 PC with a 3GHz processor. The test results are presented in Table II Referring to Table II, the NUT-ADPCM yields in all cases an improvement in both average transmission bit rate and SNR compared to ADPCM. Moreover, listening tests show a better quality using the NUT-ADPCM method than using the ADPCM for the same number of bits. The performance improvement is nevertheless obtained at the expense of a slight increase in computation time.

II. CONCLUSION

New nonuniform transmission coding systems were proposed based on nonuniform transmission methods using adaptive parametric prediction methods [START_REF] Zgheib | New fast algorithm for simultaneous identification and optimal reconstruction of non stationary AR processes with missing observations[END_REF], [START_REF] Zgheib | Lattice algorithm for adaptive stable identification and robust reconstruction of non stationary AR processes with missing observations[END_REF]. Nonuniform transmission was studied in the differential coding case. A nonuniform transmission method for quantized error prediction was proposed. As the prediction error is lower than the quantization error, this method can be considered to be a nearlossless compression method (or even a lossless compression method disregarding the quantization effect. A NUT-ADPCM coding system based on this nonuniform transmission method was described using the same adaptive quantization method as the ADPCM described in the G.726 [1] recommendations. In contrast to the ADPCM, the adaptive-prediction method used in the NUT-ADPCM is the one proposed in [START_REF] Zgheib | Lattice algorithm for adaptive stable identification and robust reconstruction of non stationary AR processes with missing observations[END_REF], thus the signal was modeled by an AR process. In addition, transmission decisions for the quantized prediction errors and the flags were introduced. The NUT-ADPCM coding system was compared to the ADPCM as described in the G.726 recommendations through simulations with speech signals. Improvements in both SNR and average transmission bit rate were observed. Thus, when the prediction errors are coded with four bits (the ADPCM working at 32kbits/s), the NUT-ADPCM coding offers an improvement of 10% in the average bit rate and an improvement of 6dB in the SNR compared to the ADPCM coding system.

Fig. 1 .

 1 Fig. 1. Simplified diagram of a non-uniform transmission ADPCM coder decoder.

  . The term CPU in the table indicates the computation time in seconds required by each of the methods.


Michel Bellet

Andres Rodriguez-Villa

Jean-François Agassant

Finite element and automatic remeshing methods for the simulation of complex blow molded polymer components

This paper presents a three dimensional finite element model of the extrusion blow molding process. The code Tform3 has the following characteristics: membrane formulation, linear triangle elements, updated Lagrangian implicit formulation, viscoelastic differential constitutive equations. The paper presents a brief recall of the formulation and then addresses three key issues of the simulation: automatic identification of constitutive equation parameters, automatic remeshing, coupling between gas pressure and inflation. An example of application to the extrusion blow molding of a bottle is presented.

INTRODUCTION

The extrusion blow-molding process consists in blowing an extruded cylindrical thin tube, or parison, of molten polymeric material (polyethylene, polyester, polyvinylchloride) into a mold. The process is applied to the production of hollow components such as bottles and tanks.

One of the key problems in product development is the fulfilment of a prescribed thickness distribution of the final part. Obviously numerical simulation can help in process optimisation and the finite element software Tform3 has been developed in that sense. After a brief recall of its basic concepts, the paper will focus on three key issues: automatic identification of rheology parameters, automatic remeshing, and coupling between internal gas pressure and inflation kinetics.

TFORM3. FEM MEMBRANE MODEL

Due to the reduced thickness of inflated parisons with respect to their curvature radii, it is possible to consider them in a first approach as geometrical surfaces, or membranes. This is equivalent to neglecting the mechanical power dissipated by flexion and shear through thickness with respect to the power of tangential expansion. In other words, the stress state is considered planar in the tangent plane to each point of the membrane.

Membrane assumption is done by most of the authors for polymer blowing finite element analysis (thermoforming or extrusion blow-molding process). Let us quote [START_REF] Warby | Finite element model of viscoelastic membrane deformation[END_REF][START_REF] Charrier | Free and constrained inflation of elastic membranes in relation to thermoforming[END_REF][START_REF] Lorenzi | Finite element simulation of thermoforming and blow molding[END_REF][START_REF] Kouba | Computer simulation of thermoforming in complex shapes[END_REF][START_REF] Shrivastava | Large deformation analysis of nonlinear viscoelastic membranes with reference to thermoforming[END_REF][START_REF] Vantal | Numerical simulation of the thermoforming of multi-layer polymer sheets[END_REF][START_REF] Rodriguez-Villa | Finite element simulation of the extrusion blow-molding process[END_REF][START_REF] Santhanam | A design tool for blow-molding and thermoforming process[END_REF][START_REF] Zhang | Process simulation of polymer blow molding by 3-D rigid viscoplastic FEM[END_REF][START_REF] Laroche | Process modelling of complex blow-molded parts[END_REF]. It should be noted that some authors, in this context, use an explicit formulation based on a centered finite differences scheme (equivalent to a resolution for the acceleration variable) [START_REF] Rachik | Simulation numérique du soufflage et du thermoformage des plastiques[END_REF][START_REF] D'oria | Progress in numerical modelling of the thermoforming process[END_REF][START_REF] Verron | Contribution expérimentale et numérique aux procédés de moulage par soufflage et de thermoformage[END_REF], whereas others prefer an implicit formulation, using the velocity or displacement as primitive variable. As the inertia terms are generally weak or moderate in extrusion blow-molding processes, the software Tform3 has been based on an implicit scheme.

The alternatives to the membrane assumption are on one hand the shell formulation, not frequently used in this field, probably because of its low increase of precision as compared to the membrane one, and, on another hand, the volumic formulation. Let us quote [START_REF] Song | Finite element analysis of inflation of an axisymmetric sheet of finite thickness[END_REF] for thermoforming, [START_REF] Debbaut | Numerical simulation of the blow molding process[END_REF][START_REF] Hartwig | Finite element simulation of the stretch blow molding process[END_REF][START_REF] Schmidt | Viscoelastic simulation of PET stretch/blow molding process[END_REF] for blow-molding. These authors stress the point that this formulation is interesting in presence of significative shear through the thickness, due to the tooling (stretching rod, punch). Another reason is the heat transfer analysis, although this can be accounted for in membrane formulation [START_REF] Vantal | Numerical simulation of the thermoforming of multi-layer polymer sheets[END_REF].

In our formulation, we define a system of two curvilinear coordinates onto the inflated membrane. The mechanical equilibrium is expressed by the principle of virtual power according to the following equation [START_REF] Bellet | Numerical simulation of thin sheet forming processes by the finite element method[END_REF]):

0 dS * . p dS e * v * p 3 d j i ij = - σ ∀ ∫ ∫ Ω Ω v g v (1)
in which e is the local thickness, p d is the differential pressure applied to the region Ω p of the deformed parison Ω, σ σ σ σ is the Cauchy stress tensor of contravariant components σ ij , g 3 is the normal vector to the membrane, v* is any virtual velocity field of covariant derivative components v* i | j . According to the membrane assumption, i,j=1,2.

In this equation, gravity and friction effects have been neglected. Hence a sticking contact is considered, which is representative of actual contact between hot polymer and cold metallic tools.

At each time increment [t,t+∆t], given an initial balanced configuration Ω t , the problem consists in solving (1) at time t+∆t. An implicit time integration scheme for the velocity is used and therefore, using the constitutive and incompressibility equations, (1) can be solved for the single unknown velocity field v t+∆t .

Discretisation with linear triangles yields a set of non-linear equations which is solved by the Newton-Raphson method. Then the configuration is updated.

A contact algorithm computes the position of each node with respect to the tooling surfaces, which are discretised with triangles. When a node penetration has been detected, the algorithm computes the normal velocity to be prescribed during next time increment in order to bring the node back onto the tool surface. According to the sticking contact assumption, the relative velocity of contacting nodes with respect to the tool is prescribed to zero.

CONSTITUTIVE EQUATION AND IDENTIFICATION

The pioneer works on blow-molding use hyperelastic constitutive equations, invoking short forming times. The most frequent models are Mooney-Rivlin and Ogden ones: [START_REF] Charrier | Free and constrained inflation of elastic membranes in relation to thermoforming[END_REF][START_REF] Lorenzi | Finite element simulation of thermoforming and blow molding[END_REF]. This approach is still popular: [START_REF] Rachik | Simulation numérique du soufflage et du thermoformage des plastiques[END_REF][START_REF] Hartwig | Finite element simulation of the stretch blow molding process[END_REF][START_REF] D'oria | Progress in numerical modelling of the thermoforming process[END_REF][START_REF] Santhanam | A design tool for blow-molding and thermoforming process[END_REF][START_REF] Laroche | Process modelling of complex blow-molded parts[END_REF] However, molten polymer rheology is viscoelastic, and most authors naturally use an extension of hyperelastic laws, i.e. integral viscoelastic models. The pioneer work of [START_REF] Warby | Finite element model of viscoelastic membrane deformation[END_REF] in this field has been followed by [START_REF] Kouba | Computer simulation of thermoforming in complex shapes[END_REF]

Maxwell viscoelastic model

The upper-convected Maxwell model is expressed by:

σ σ σ σ = -+ p v I T
(2)

T T v v t + = λ δ δ µ 2 & ε ε ε ε (3)
where p is an arbitrary pressure, T v the extra-stress tensor, λ the relaxation time, µ the dynamic viscosity and δT v / δt the upper-convected Oldroyd's derivative:

[ ] [ ] δ δ T T Gra d v T T Gra d v v v v v T t d dt = - - (4)
It is to be noticed that this upper-convected derivative reduces in curvilinear coordinates to the total derivative of covariant components [START_REF] Bellet | Numerical simulation of thin sheet forming processes by the finite element method[END_REF][START_REF] Rodriguez-Villa | Finite element simulation of the extrusion blow-molding process[END_REF].

Identification by inverse method

This method consists in the association of an instrumented blow molding test with a procedure of automatic identification using the finite element simulation of the test (Rodriguez-Villa 1997).

The mold shape is shown at figure 1. The inflation is recorded through a lateral window, using a video recorder at 25 pictures.s -1 (the typical inflation time being in the range 0.6 to 3 s.). A recording example is given at figure 2. It provides the measurement of the radius at midheight of the parison as a function of time R(t). In addition, the internal air pressure p a (t) is recorded (fig. 3). The reproducibility of such measurements has been checked for different inflation rates and temperatures. The identification step has been carried out as follows. Each blow molding test is simulated using the finite element software Tform3, applying to the parison the pressure load which has been effectively measured during the test. Then, by successive iterations, the set of parameters P (=(λ,µ) for Maxwell law) is modified until the error between the measured and the computed curves R(t) is reduced to a minimum. This error is defined as:

( ) ( ) E ( ) ( ) P P = - = = ∑ ∑ R R R i ca l i exp i 1 n exp i exp 2 i 1 n exp 2 (5)
in which nexp is the number of sampling points. The minimisation of E is achieved by a Gauss-Newton method, needing the first derivatives of E, and therefore the derivatives of R i ca l with respect to each parameter of P. Because of the very low computation time of the direct FEM simulation of the test, those derivatives are approached by finite differences. Figure 4 shows the comparison, for a single test, between experimental and computed (after optimization) curves R(t). It is found, as expected, that neither the linear elastic nor the Newtonian laws are able to account for the actual radius evolution under the measured pressure load. This clearly demonstrates the viscoelastic character of molten polymer during blow molding. For the three viscoelastic laws which have been tested (Maxwell, Oldroyd-B, Phan-Tien & Tanner), it has always been possible to exhibit a set of parameters resulting in a low error E. Table 1 gives the parameter values identified for a Maxwell constitutive equation on four blow molding tests for different inflation rates and temperatures. The dependence of viscosity on strain rate is shown in figure 5, with reference to measurements issued from shear and tensile tests. The viscosity is found higher in blowing conditions (biaxial stretching) than in shear conditions and comparable to the elongational viscosity found in uniaxial tensile conditions at similar strain levels. The dependence of λ and µ on strain rate is similar to the one found in shear conditions. It can be concluded that the Maxwell law is representative -at least qualitatively -of actual phenomena. However, a much better quantitative agreement should be obtained on the whole set of tests using a generalized Maxwell model.

T (°C) ε eq (s -1 ) µ (Pa.s) λ (s) E (

AUTOMATIC AND ADAPTIVE REMESHING

During blow molding, the material encounters large deformations. Hence the use of a simple updated Lagrangian formulation, consisting in convecting the mesh with the material, would result in a loss of precision due to mesh coarsening or mesh distortion. This can be avoided by regenerating the mesh from time to time.

The problem consists in remeshing a warped triangulated surface. The algorithm is initially issued from the remeshing procedure of Forge3 software [START_REF] Coupez | Parallel adaptive remeshing in 3D moving mesh finite element numerical grid generation[END_REF]. The basic principle is to extract around each node or each mesh edge a sub-topology composed of the elements including the considered node or edge, as explained on figure 6. Alternative topologies are then built and evaluated using a quality criterion (see hereunder). The best topology is then fixed and reinserted in place of the initial one and another node or edge is examined. The procedure is stopped when no more topological changes are found necessary. The quality Q τ of a sub-topology τ is taken equal to the minimum quality of its elements e:

Q Min Q e e τ τ = ∈ (6)
where Q e denotes element quality, which takes into account both shape ratio and size.

The shape quality is defined by the ratio:

( )

Q f S P sh e e e = 2 (7)
where S e is the element surface, P e its perimeter and f a coefficient such that the shape quality of equilateral triangles is equal to 1. The size quality is defined by considering the deformation between the equilateral triangle of aimed size h obj and the actual element (fig. 7). in which it can be seen that the quality tends to 1 as C:C tends to 1 (which results in fact in producing elements smaller than the expected size h obj ). In addition, e sz Q is all the more discriminant as e sh Q is high. Finally, the element quality is determined by the product:

( ) ( ) Q Q Q f S P Q e sh
Let us define now the aimed size h obj . It results from an estimator of the geometric error associated with the current discretisation of the membrane. Such a problem has been addressed by [START_REF] Bonet | Error estimators and enrichment procedures for the finite element analysis of thin sheet large deformation processes[END_REF] and the present approach is similar. It is based first on a nodal smoothing of the tangent basis:

where the sum is extended to the elements e surrounding the considered node n. Then the deformation between the current and the smoothed surface can be measured by the Green-Lagrange type tensor:

( )

L g g ij ij ij = - 1 2 ~ (11)
This permits to define the local error as:

It can be shown that this error varies quadratically with element size h in the case of a regular cylindrical surface (Rodriguez-Villa 1997). Generalizing the same dependency, the aimed size h obj can be defined by:

E E h h obj obj = (13)
where E obj and E are the objective error and the current local error respectively.

Finally, modifications of local topology are controlled by examining their impact on mean curvature variations. It can be noticed that the measured pressure is considerably lower (by a ratio ranging from 30 to 100) than the nominal one in the upstream air feeding system. Moreover, the pressure is far from being constant during the inflation stage. The main pressure peak corresponds to the beginning of parison inflation: the radius significantly increases and at the same time the pressure decreases. This clearly demonstrates the coupling between inflation kinetics and pressure evolution.

In order to account for this phenomenon, a thermodynamic model has been developed [START_REF] Rodriguez-Villa | Etude théorique et expérimentale de l'extrusion-soufflage de corps creux en polymère[END_REF]. The specific air flow rate q e (mass per time unit) is assumed constant, the air being in a given thermodynamic state (pressure p e , specific mass ρ e and temperature T e ). The energy variation of the air located inside the parison results from the energy input due to incoming air, from the energy dissipated by external forces and from the heat exchanged with the surroundings. This energy change is split into variations of internal energy, kinetic energy and potential energy: Assuming that the blowing power e e q e remains constant, it is then possible to determine it by a simple test of air blowing in a reservoir of constant volume V c . In this case, dV = 0 and we have: The blowing power can be estimated by measuring the initial pressure rate dp a /dt, denoted & P 0 . This parameter is a characteristic of the blowing machine. The pressure evolution during an extrusion blow molding on the same machine will be given by:

γ γ p dV dt p dV dt V dp dt V P a a t m a c -- + = ( ) & 1 0 (20)
A validation of this model has been carried out by coupling it to the computation of inflation kinetics in Tform3. Figure 8 gives an example -for one of the previous identification tests -of the good agreement between the measured pressure and the predicted one by the coupled model. From this result, it can be concluded that this coupling is definitely necessary to simulate extrusion blow molding in a realistic manner.

APPLICATION

Figure 9 shows a simulation of the extrusion blow molding of a detergent bottle made of high density polyethylene. The part encounters severe thickness gradients around its hollow handle. The initial geometry of the parison is: radius 60 mm, height 500 mm, thickness 1.5 mm. Due to the vertical symmetry, only half of the parison is computed. It is assumed to be fixed at the upper end (output of extrusion die) and a zero vertical displacement is prescribed at the lower end.

On this figure, the progressive onset of contact with the mold can be seen, during the closure of the mold, before the blowing itself. The efficiency of adaptive remeshing procedures is demonstrated.

Figure 10 presents the final computed thickness distribution. It can be seen that the parison has been laterally displaced with respect to the extrusion axis. This has been found necessary to prevent a folding initiating during the handle blowing. This has been confirmed experimentally and a good compromise has to be found between this displacement and the associated low thickness values obtained at the bottom of the bottle, on opposite side. The computing time is 1130s on a single processor of a Silicon Graphics R10000 origin 200 machine (412 time increments, maximum node number: 5000).

  with K-BKZ model, Shrivastava & Tang 1993, Verron 1997 with Christensen model. Others prefer differential viscoelastic models: Debbaut et al. 1993 with Giesekus law, Rachik et al. 1994 with White-Metzner model, Schmidt et al. 1996 and Rodriguez-Villa et al. 1995, 1997 with Maxwell, Oldroyd-B and Phan-Tien & Tanner models. Those last three constitutive equations are implemented in Tform3.

Figure 1 .

 1 Figure 1. Instrumented mold.

Figure 2 .

 2 Figure 2. Video acquisition during inflation.

  Figure

Figure 4 .

 4 Figure 4. Automatic identification of the parameters of several constitutive equations on a blow molding test.

Figure 5 .

 5 Figure 5. Comparison of viscosities measured in shear test, uniaxial tensile test and in the present blowing test.

Figure 6 .

 6 Figure 6. Extraction of sub-topology around an edge or a node. Choice among alternative sub-topologies and reinsertion.

FFigure 7 .

 7 Figure 7. Transformation between aimed and actual element. Denoting C the Cauchy dilatation tensor of transformation F, i.e. C F F = T , the size quality is defined by the formula: ( ) Q Q sz e sh e = +

  the parison has been recorded as a function of time during identification tests (fig.3).

  potential energy variations, and assuming transformation adiabaticity, (14) can be written, expressing the work of external forces: of incoming air is denoted e e , its internal energy u e and its velocity v e . c v is the specific heat at constant volume, we have: of specific heats. The equation yielding the air pressure evolution p a is then easily obtained:

Figure 8 .

 8 Figure 8. Comparison of measured and predicted pressure.

Figure 9 .

 9 Figure 9. Simulation of the blow molding of a detergent bottle.

Figure 10 .

 10 Figure 10. Final configuration of the simulated bottle. Deformed mesh and thickness profile (mm).

Table 1 .

 1 Maxwell viscosity and relaxation time identified by inverse method on four different blowing tests.

	%)