Denys Dutykh
email: denys.dutykh@ucd.ie

Olivier Goubet
email: olivier.goubet@u-picardie.fr

DERIVATION OF DISSIPATIVE BOUSSINESQ EQUATIONS USING THE DIRICHLET-TO-NEUMANN OPERATOR APPROACH

Keywords: Boussinesq equations, viscosity, dissipation, dispersive waves, boundary layer 2010 Mathematics Subject Classification. 76B25 (primary), 76B07, 65M70 (secondary)

The water wave theory traditionally assumes the fluid to be perfect, thus neglecting all effects of the viscosity. However, the explanation of several experimental data sets requires the explicit inclusion of dissipative effects. In order to meet these practical problems, the theory of visco-potential flows has been developed (see P.-F. Liu & A. Orfila [42] and D. Dutykh & F. Dias [25]

). Then, usually this formulation is further simplified by developing the potential in an entire series in the vertical coordinate and by introducing thus, the long wave approximation. In the present study we propose a derivation of dissipative Boussinesq equations which is based on two various asymptotic expansions of the Dirichlet-to-Neumann (D2N) operator. Both employed methods yield the same system at the leading order by different ways.

Introduction

The water wave problem has a long history which counts more than two hundred years [START_REF] Craik | The origins of water wave theory[END_REF]. Traditionally the fluid under consideration is supposed to be perfect which allowed to reveal several useful variational structures [START_REF] Petrov | Variational statement of the problem of liquid motion in a container of finite dimensions[END_REF]44,67,8]. However this assumption is not valid for all time scales and for all physical situations. The explanation of several laboratory experiments and some modeling problems require the explicit inclusion of weak dissipative effects. In the prominent work by J. [START_REF] Bona | An Evaluation of a Model Equation for Water Waves[END_REF] [START_REF] Bona | An Evaluation of a Model Equation for Water Waves[END_REF] one can find the following passage emphasizing the importance of dissipative effects: ". . . it was found that the inclusion of a dissipative term was much more important than the inclusion of the nonlinear term, although the inclusion of the nonlinear term was undoubtedly beneficial in describing the observations."

Another more recent experimental study emphasizes more specifically the boundary layer effects [START_REF] Hammack | Solitary wave collisions[END_REF]: "A third source of experimental differences is residual boundary layer motions that are left behind as a solitary wave propagates in the channel. In both the co-and counterpropagating binary collision experiments solitary waves encounter the boundary-layer wakes of the other wave. This wake does have a small effect on wave speeds that can be significant in our data analysis. Detailed measurements of these boundary layer motions were not made; hence, it is not known how reproducible they are."

The complete model of viscous free-surface flows should consider two phases (air and water) separated by an interface. Each phase is described by incompressible Navier-Stokes equations with respective viscosities. Such numerical models currently exist [START_REF] Raval | Numerical simulation of viscous, nonlinear and progressive water waves[END_REF][START_REF] Coifman | Nonlinear harmonic analysis and analytic dependence[END_REF], however they allow to simulate only a few wavelengthes in reasonable times. Consequently, several approximate models have been developed to take into account weak viscous effects. One of the first simplifications consists in replacing the upper light fluid by a free surface. This research direction was successfully pursued in several studies mentioned below.

Historically various dissipative terms have been incorporated in long wave equations such as Burgers [START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF], KdV [START_REF] Ott | Damping of Solitary Waves[END_REF]34,[START_REF] Grimshaw | Damping of large-amplitude solitary waves[END_REF] or some Boussinesq type systems [36,37,[START_REF] Chen | Long-time asymptotic behaviour of dissipative Boussinesq systems[END_REF][START_REF] Dutykh | Dissipative Boussinesq equations[END_REF][START_REF] Dutykh | Mathematical modelling of tsunami waves[END_REF]35,[START_REF] Dutykh | Mathematical modeling in the environment[END_REF][START_REF] Chen | Decay of solutions to a water wave model with a nonlocal viscous dispersive term[END_REF]. Some attempt has been also made to incorporate ad-hoc dissipative terms into the complete water wave problem as well [45,54,55,[START_REF] Joseph | Potential Flow of Viscous Fluids: Historical Notes[END_REF][START_REF] Wu | A note on stabilizing the Benjamin-Feir instability[END_REF]38].

Recently, the theory of visco-potential flows has been developed for weakly dissipative regimes [START_REF] Liu | Viscous effects on transient long-wave propagation[END_REF][START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF][START_REF] Liu | Boundary layer flow and bed shear stress under a solitary wave[END_REF][START_REF] Dutykh | Visco-potential free-surface flows and long wave modelling[END_REF][START_REF] Dutykh | Group and phase velocities in the free-surface visco-potential flow: new kind of boundary layer induced instability[END_REF]. The main idea consists in writing the Helmholtz decomposition of the velocity field for linearized Navier-Stokes equations with free surface. These equations can be identically satisfied by a potential velocity field but not the boundary conditions. Then, the rotational part is expressed asymptotically (and hence, approximatively) in terms of the velocity potential and free surface elevation. In this way we introduce some corrections to boundary conditions on the free surface and/or on the bottom. In deep water case this question was studied by F. [START_REF] Dias | Theory of weakly damped free-surface flows: a new formulation based on potential flow solutions[END_REF] [START_REF] Dias | Theory of weakly damped free-surface flows: a new formulation based on potential flow solutions[END_REF]. In the finite depth case the situation is more complex since a boundary layer at the bottom should be taken into account [START_REF] Liu | Viscous effects on transient long-wave propagation[END_REF][START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF]. This physical effect modifies the bottom boundary condition by adding a nonlocal term in time which represents mathematically a half-order integral. Consequently, a memory effect is introduced into the model in the finite depth. Some attempt to the efficient discretization of this fractional integral has been performed during last years [START_REF] Torsvik | An efficient method for the numerical calculation of viscous effects on transient long-waves[END_REF][START_REF] Dutykh | Visco-potential free-surface flows and long wave modelling[END_REF][START_REF] Chen | Decay of solutions to a water wave model with a nonlocal viscous dispersive term[END_REF][START_REF] Chehab | Numerical study of a family of dissipative KdV equations[END_REF].

The main scope of the present study is to combine the visco-potential approach with Dirichlet-to-Neumann (D2N) operator formulation of the water wave problem [START_REF] Petrov | Variational statement of the problem of liquid motion in a container of finite dimensions[END_REF]67,[START_REF] Craig | Nonlinear modulation of gravity waves: a rigorous approach[END_REF]12] which can be successfully used for water wave modeling purposes [10,11,[START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF]. Using this novel visco-potential formulation and an asymptotic expansion of the D2N operator, we derive a long wave model of Boussinesq type. We underline that our derivation remains at the formal level. More theoretical study will be needed to justify our current developments.

The present study is organized as follows. In Section 2 we present the visco-potential formulation with D2N operator, while dissipative Boussinesq equations are derived in Section 3 following two different approaches presented in Sections 3.1 and 3.2. Finally, this study is ended by outlining main conclusions in Section 5. In Section 4 we discuss some results related to dissipative Boussinesq system and the dispersion relation analysis.

Mathematical modeling

Consider an ideal incompressible fluid of constant density ρ. The vertical projection of the fluid domain Ω is a subset of R 2 . The horizontal independent variables are denoted by x = (x, y) and the vertical one by z. The origin of the cartesian coordinate system is chosen such that the surface z = 0 corresponds to the still water level. The fluid is bounded below by the bottom z = -h(x, t) and above by the free surface z = η(x, t). We assume that the total depth H(x, t) := h(x, t) + η(x, t) remains positive H(x, t) ≥ h 0 > 0 at all times t ∈ [0, T ]. The sketch of the physical domain is shown in Figure 1.

Remark 1. We make the classical assumption that the free surface is a graph z = η(x, t) of a single-valued function. It means in practice that we exclude some interesting phenomena, (e.g. wave breaking phenomena) which are out of the scope of this modeling paradigm.

The governing equations of the classical water wave problem are the following [39, 56, 48, 64]:

∇ 2 φ + ∂ 2 zz φ = 0, (x, z) ∈ Ω × [-h, η], (2.1)
∂ t η + ∇φ • ∇η -∂ z φ = 0, z = η(x, t), (2.2)
∂ t φ + 1 2 |∇φ| 2 + 1 2 (∂ z φ) 2 + gη = 0, z = η(x, t), (2.3)
∂ t h + ∇φ • ∇h + ∂ z φ = 0, z = -h(x, t), (2.4)
with φ the velocity potential (by definition the velocity field (u, v) = ∇ x,z φ), g the acceleration due to gravity force and ∇ = (∂ x , ∂ y ) denotes the gradient operator in horizontal Cartesian coordinates. Incompressibility leads to Laplace's equation for φ. The main difficulty of the water wave problem lies on the nonlinear free boundary conditions. Equations (2.2) and (2.4) express the free-surface kinematic condition and bottom impermeability respectively, while the dynamic condition (2.3) expresses the free surface isobarity. Remark 2. Surface tension effects can also be included in the water wave problem. In this case, the Bernoulli condition (2.3) has to be modified (see [18] for a general review of this topic):

∂ t φ + 1 2 |∇φ| 2 + 1 2 (∂ z φ) 2 + gη = σ ρ ∇ • ∇η 1 + |∇η| 2 , z = η(x, t),
where σ is the surface tension coefficient. However, this effect is not of primary importance in our study and, consequently, will be ignored.

Below we will need the unitary exterior normals to the fluid domain. It is straightforward to obtain the following expressions for the normals at the free surface and bottom respectively:

nf = 1 1 + |∇η| 2 -∇η 1 , nb = 1 1 + |∇h| 2 -∇h -1 .
Russian physicists A. Petrov (1964) [START_REF] Petrov | Variational statement of the problem of liquid motion in a container of finite dimensions[END_REF] and V. [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF] [67] proposed a more compact formulation of the water wave problem (2.1) -(2.4) based on the Hamiltonian structure. Let us introduce the trace of the velocity potential at the free surface:

ϕ(x, t) := φ(x, η(x, t), t).
This variable plays the role of generalized momentum in the Hamiltonian description of water waves [67,15]. The second canonical variable is the free surface elevation η. Another important ingredient is the normal velocity at the free surface v n which is defined as:

v n (x, t) := 1 + |∇η| 2 ∂φ ∂n f z=η = (∂ z φ -∇φ • ∇η)| z=η . (2.5)
Using all these elements, the boundary conditions (2.2) and (2.3) on the free surface can be rewritten in terms of ϕ, v n and η [START_REF] Craig | Nonlinear modulation of gravity waves: a rigorous approach[END_REF]12,[START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF] as following:

∂ t η -D η (ϕ) = 0, ∂ t ϕ + 1 2 |∇ϕ| 2 + gη - 1 2(1+|∇η| 2 ) D η (ϕ) + ∇ϕ • ∇η 2 = 0, (2.6)
where we introduced the so-called classical Dirichlet-to-Neumann (D2N) D η operator [9, 12], which maps the velocity potential at the free surface ϕ to the normal velocity v n :

D η : ϕ → v n = 1 + |∇η| 2 ∂φ ∂ nf z=η ∇ 2 φ + ∂ 2 zz φ = 0, (x, z) ∈ Ω × [-h, η], φ = ϕ, z = η, 1 + |∇h| 2 ∂φ ∂n b = ∂ t h, z = -h.
This operator takes the free surface shape η(x, t) as a parameter and the velocity potential trace ϕ(x, t) as an argument (Dirichlet boundary condition). As a result it returns the normal velocity

v n = 1 + |∇η| 2 ∂φ ∂ nf z=η
at the free surface (Neumann data) which explains the common name of this operator.

Remark 3. The function h(x, t) represents the ocean bathymetry (depth below the still water level, see Figure 1) and is assumed to be known. The dependence on time is included in order to take into account the bottom motion during, for example, an underwater earthquake [16, [START_REF] Dutykh | Water waves generated by a moving bottom[END_REF][START_REF] Dutykh | Linear theory of wave generation by a moving bottom[END_REF]. Since in the present study we are mainly interested in quantifying a solitary wave amplitude decay under the viscous damping, we will take hereafter the even bottom case h(x, t) ≡ h = const.

Now we are going to introduce weak dissipative effects directly into the Petrov-Zakharov formulation (2.6) in the spirit of the visco-potential flows theory [START_REF] Liu | Viscous effects on transient long-wave propagation[END_REF][START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF][START_REF] Dutykh | Group and phase velocities in the free-surface visco-potential flow: new kind of boundary layer induced instability[END_REF]:

∂ t η -Dη (ϕ) = 2ν∇ 2 η, ∂ t ϕ + 1 2 |∇ϕ| 2 + gη - 1 2(1+|∇η| 2 ) Dη (ϕ) + ∇ϕ • ∇η 2 = 2ν∇ 2 ϕ, (2.7)
where ν is the kinematic (molecular or eddy) viscosity and the modified D2N operator D is defined as: Dη :

ϕ → v n = 1 + |∇η| 2 ∂φ ∂ nf z=η ∇ 2 φ + ∂ 2 zz φ = 0, (x, z) ∈ Ω × [-h, η], φ = ϕ, z = η, ∂ z φ = -ν π t 0 φzz √ t-τ dτ, z = -h.
The local dissipative terms in (2.7) are due to viscous effects in the bulk of the fluid while the integral term in the bottom kinematic condition takes into account boundary layer effects at the bottom. For more details on the modeling and derivation of visco-potential flows formulation we refer to [START_REF] Liu | Viscous effects on transient long-wave propagation[END_REF][START_REF] Dutykh | Mathematical modelling of tsunami waves[END_REF][START_REF] Dutykh | Visco-potential free-surface flows and long wave modelling[END_REF]. In the next Section we will simplify dynamic equations (2.7) by applying a long wave approximation.

Remark 4. The bottom kinematic boundary condition can be equivalently rewritten using the Laplace equation which is still valid on boundaries under the assumption of their regularity:

∂ z φ = ν π t 0 ∇ 2 φ √ t -s ds, z = -h.
The notation can be greatly simplified if we apply the Laplace transform to both sides of the last identity [START_REF] Widder | The Laplace transform[END_REF][START_REF] Gradshteyn | Tables of Integrals, Series, and Products[END_REF]:

∂ z φ = ν τ ∇ 2 φ,
where τ is the Laplace transform parameter (to shorten the notation the images of the Laplace transform are denoted by the same symbol). We stress that the use of the Laplace transform is not fundamental to the subsequent derivation. Alternatively we could denote the half-integral by another symbol. We find it natural to use the apparatus of the operational calculus.

Dimensionless variables

In order to apply an asymptotic expansion in the following Section, we have to introduce dimensionless variables to reveal small parameters [2, [START_REF] Dutykh | Dissipative Boussinesq equations[END_REF]19,[START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF].

Unscaled independent variables are defined as follows:

x ′ = x ℓ , z ′ = z h , t ′ = t ℓ/ √ gh ,
where ℓ is a typical wavelength, h is the average water depth and √ gh gives the phase speed of infinitely long shallow water waves. The scaling of dependent variables needs one more characteristic length -the typical wave amplitude a:

η ′ = η a , φ ′ = φ ℓag/ √ gh .
In dimensionless variables the governing equations (2.7) take the following form (we dropped out the primes for the sake of convenience):

∂ t η - 1 µ 2 Dεη (ϕ) = 2ν∇ 2 η, (2.8) ∂ t ϕ + ε 2 |∇ϕ| 2 + η - ε µ 2 Dεη (ϕ) + εµ 2 ∇ϕ • ∇η 2 2(1 + ε 2 µ 2 |∇η| 2 ) = 2ν∇ 2 ϕ, (2.9)
and the D2N operator Dεη in unscaled variables becomes:

Dεη : ϕ → v n = 1 + ε 2 µ 2 |∇η| 2 ∂φ ∂ nf z=εη = ∂ z φ -εµ 2 ∇φ • ∇η z=εη ∇ 2 φ + 1 µ 2 ∂ 2 zz φ = 0, (x, z) ∈ Ω × [-1, εη], φ = ϕ, z = εη, ∂ z φ = - ν πµ 2 t 0 φzz √ t-s ds, z = -1.
Remark 5. In the light of Remark 4, the bottom boundary condition can be more conveniently rewritten in dimensionless variables using the Laplace transform:

∂ z φ = - 1 µ ν τ φ zz ≡ µ ν τ ∇ 2 φ.
We introduced above three classical scaling parameters: the nonlinearity ε := a h , the dispersion µ 2 := h ℓ 2 and the inverse Reynolds number also denoted by ν for the sake of convenience ν → ν √ gh ℓ . First two parameters are classical in long wave theory [51, 1, 3, 63, 49, 2, 24] while the last one is less classical in this context and it is present due to viscous effects. In the so-called Boussinesq regime it is additionally assumed that nonlinearity and dispersion are of the same order, i.e. the Stokes-Ursell number [START_REF] Ursell | The long-wave paradox in the theory of gravity waves[END_REF] is of the order of unity:

St := ε µ 2 = O(1).
(2.10)

The last assumption St ∼ 1 will be systematically used in the following Section. It will allow us to neglect higher order terms since:

ε 2 ∼ εµ 2 = O(µ 4 ).

D2N operator asymptotic expansions

In the current Section we simplify the water wave problem (2.7) by expanding asymptotically the D2N map. We apply two different methods which lead to the same result. The first method is called the shape derivative method and consists in expanding the D2N map in the nonlinearity parameter ε. The second method presented in Section 3.2 expands the velocity potential in powers of the dispersion parameter µ 2 . It is worth to point out that if the Stokes-Ursell number is of the order of one (the so-called Boussinesq regime [19]), i.e. ε ∼ µ 2 , both methods lead to the same asymptotic model.

To shorten the notation and make the derivation as clear as possible, we restrict our attention to the 2D case. The horizontal variable is denoted simply by x hereafter. Consequently, we will have to compute an asymptotic approximation to the following expression:

Dεη (ϕ) = φ z -εµ 2 φ x η x z=εη .
(3.1)

Shape derivative method

In this section we follow in great lines the study [40]. However, the original derivation procedure was greatly simplified to fit our particular purposes.

We know that the velocity potential φ(x, z) solves the following elliptic boundary-value problem:

φ xx + 1 µ 2 φ zz = 0, (x, z) ∈ R × [-1, εη] (3.2) φ = ϕ, z = εη, (3.3)
φ z = µ ν τ φ xx , z = -1. (3.4)
We will find a solution φ(x, z; t) to this problem in the form:

φ(x, z; t) = φ (0) + εφ (1) + o(ε).

(3.5)

In the following the dependence on time will be omitted to shorten the notation. Now we have to construct zero and first order approximations φ (0) and φ (1) .

Solution in the fixed domain

It can be easily seen that the function φ (0) satisfies the following set of equations:

φ (0) xx + 1 µ 2 φ (0) zz = 0, (x, z) ∈ Ω 0 = R × [-1, 0] φ (0) = ϕ, z = 0, φ (0) z = µ ν τ φ (0) xx , z = -1.
Since the domain is a strip, we can employ the Fourier transform in the horizontal variable x in order to solve this problem. Throughout this section by f (k) we will denote the Fourier transform of a function f (x) ∈ L 2 (R). In Fourier variable, we then solve the ODE φ(0) zz -µ 2 k 2 φ(0) = 0, supplemented with boundary conditions φ(0) = φ at z = 0 and φ(0

) z = -µk 2 ν τ φ(0) at z = -1.
After some simple computations we find

φ(0) (k, z) = Rν (k, z)[ φ(k)],
where the pseudo-differential operator R ν (x, z)[•] is defined through its symbol which reads:

Rν (k, z)[•] := cosh µ|k|(z + 1) -|k| ν τ sinh µ|k|(z + 1) cosh(µ|k|) -|k| ν τ sinh(µ|k|) •
For our purposes it is sufficient to use an approximate operator ∂ z Rν (k, z)[•] at z = 0 to the order O(µ 2 ):

Rν (k, z)[•] = k 2 1 - k 2 µ 2 3 + µk 2 - 1 µ ν τ + 4µ 2 k 4 3 -k 2 ν τ -µk 4 ν τ 3 2 + o(µ 2 ) ,
which leads to the following term appearing in the D2N map Dεη (ϕ):

1 µ 2 φ (0) z (x, 0) = -ϕ xx - µ 2 3 ϕ xxxx + 1 µ ν τ (ϕ xx + µ 2 ϕ xxxx ) - ν τ ϕ xxxx + µ 2 3 ϕ xxxxxx + µ ν τ 3 2 ϕ xxxxxx + o(µ 2 ).
In the last expression we can neglect higher order viscous terms to obtain a simplified representation:

1 µ 2 φ (0) z (x, 0) = -ϕ xx - µ 2 3 ϕ xxxx + 1 µ ν τ (ϕ xx + µ 2 ϕ xxxx ) - ν τ ϕ xxxx + o(ν + µ 2 ).

The shape derivative

The next step consists in computing the first order approximation φ (1) (x, z). The main problem is that in a difference with the previous case, the function φ (1) is defined on a perturbed domain. In order to proceed, we will make a change of variables to work on a fixed domain:

(x, z) ∈ Ω ε -→ (x, y = r(x, z) = z -εη 1 + εη ) ∈ Ω 0 . (3.6)
It is easy to see that the Jacobian J of the last transformation is equal to J = r ′ z = 1 1+εη . Throughout this section all functions defined on the straightened domain Ω 0 will be denoted by tildes, e.g. φ(x, y(x, z)) = φ(x, z).

Let Φ ∈ C ∞ 0 (Ω ε ) be a test function. Then, equation (3.2) can be rewritten in a weak form:

1 µ 2 Ωε φ z Φ z dx dz + Ωε φ x Φ x dx dz = 0,
or after performing the change of variables (3.6):

1 µ 2 Ω0 φy ψy dx dy 1 + εη + Ω0 (1 + εη) φx -(y + 1)εη x φy ψx - (y + 1)εη x 1 + εη dx dy = 0.
We can substitute into the last integral identity the asymptotic expansion (3.5) recasted in new variables: φ = φ(0) + ε φ(1) + o(ε).

At the lowest order it will give us the weak form of the Laplace equation (3.2) for φ(0) which is satisfied identically. At the first order in ε we obtain the following integral identity:

1 µ 2 Ω0 φ(1) y ψy dx dy + Ω0 φ(1)
x ψx dx dy = = Ω0 η µ 2 φ(0) y ψy + (y + 1)η x ( φ(0)

x ψy + φ(0) y ψx ) -η φ(0)

x ψx dx dy.

Switching back to the differential form and performing some simplifications yield this inhomogeneous elliptic PDE for the first order approximation φ(1) function:

1 µ 2 φ(1) yy + φ(1) xx = 2η µ 2 φ(0) yy + 2(y + 1)η x φ(0) xy + η xx (y + 1) φ(0) y , (x, y) ∈ R × [-1, 0]. (3.7)
This equation should be completed by appropriate boundary conditions which are obtained by expanding corresponding conditions (3.3), (3.4) of the complete problem:

φ(1) = 0, y = 0, (3.8) φ(1) y -µ ν τ φ(1) xx = η φ(0) y , y = -1. (3.9)
In order to solve this problem we will consider the following ansatz for the solution: φ(1) (x, y) ≈ Φ(x, y) = (y + 1)η(x) φ(0) y (x, y).

(3.10)

where the dependence on time can be postponed until considering the dynamic free surface boundary conditions. Substituting expansion (3.13) into equations (3.2) -(3.4) leads the following infinite sequence of elliptic problems:

φ (n) zz + φ (n-2) xx = 0, (x, z) ∈ R × [-1, εη] φ (n) = δ 0n ϕ, z = εη, φ (n) z = ν τ φ (n-1) xx , z = -1,
where φ (-2) = φ (-1) ≡ 0 by definition and δ 0n is the usual Kronecker symbol (equal to 1 if n = 0 and equal to 0 otherwise). Since we are going to neglect all the terms of the order O(µ 4 ) and higher, we have to compute only first three terms (φ (0) , φ (1) and φ (2) ) in the asymptotic expansion (3.13). After some algebraic computations that we omit here, we obtain the following expressions of the first three solutions in the infinite hierarchy:

φ (0) (x, z) = ϕ(x), (3.14)
φ (1) (x, z) = ν τ (z -εη)ϕ xx , (3.15) φ (2) (x, z) = - 1 2 (z -εη)(z + εη + 2)ϕ xx - ν τ ν(z -εη) (1 + εη)ϕ xx xx . (3.16)
After substituting relations (3.14) -( 3.16) into the definition of the D2N map (3.1) and neglecting higher order terms, we derive exactly the same expression (3.12) for Dεη (ϕ) as we obtained in the previous Section 3.1.2. Consequently, both derivation methods yield the same result at least for asymptotic orders considered in this study. We have to note that the second method seems to lead to the result in a more direct way. However, the theoretical and practical potential of the shape derivative method may not be fully explored yet. Each method allows to derive models that are relevant in the range of corresponding small parameters: the first one is suitable to derive a family of models with respect to the nonlinearity parameter ǫ, while the second one is relevant for the dispersion parameter µ at the price to solve a hierarchy of D2N problems.

The choice of the velocity variable is far from being unique [START_REF] Bona | A model for the two-way propagation of water waves in a channel[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF]2], but we prefer not to enter into the debate.

The asymptotic representation (3.14) -(3.16) can be averaged over the water depth to provide us the corresponding expression for the depth-averaged velocity potential:

φ = ϕ -µ ν τ ϕ xx + µ 2 3 ϕ xx + O(µ 3 + νµ 2 ).
In order to invert this relation and, thus, express the velocity potential at the free surface ϕ in terms of φ, it is better to work with Fourier transforms [2]:

φ = 1 + µ ν τ k 2 - µ 2 3 k 2 φ + O(µ 3 ).
The last relation can be easily inverted within the same asymptotic accuracy:

φ = 1 + µ ν τ k 2 - µ 2 3 k 2 -1 φ + O(k 4 ) = 1 -µ ν τ k 2 + µ 2 3 k 2 φ + O(k 4 ).
In physical space we obtain the requested relation between ϕ and φ:

ϕ = φ + µ ν τ φxx - µ 2 3 φxx + O(µ 3 ).
The direct substitution of the last relation into governing equations (4.2), (4.1) and inverting the Laplace equation leads the following set of Boussinesq type equations which represent a viscid generalization of the celebrated Peregrine system over the flat bottom [51]:

η t + (1 + εη)ū x = 2νη xx + 1 µ ν π t 0 (ū x -µ 2 3 ūxxx )(x, s) √ t -s ds, (4.3) ūt
+ εūū x + η x - µ 2 3 ūxxt = 2ν ūxx - ∂ ∂t 1 µ ν π t 0 µ 2 ūxx (x, s) √ t -s ds , (4.4)
where ū := φx is the depth-averaged velocity variable. To our knowledge this is the first dissipative Boussinesq system which contains a half-order derivative and a half-order integral in the same time.

Dispersion relation analysis

Let us consider the linearized version of the viscous Peregrine system (4.3), (4.4) with dimensionless parameters ε, µ 2 set for convenience to the unity (it is always possible by choosing the appropriate scaling):

η t + ūx = 2νη xx + ν π t 0 (ū x -1 3 ūxxx )(x, s) √ t -s ds, (4.5) ūt
+ η x - 1 3 ūxxt = 2ν ūxx - ∂ ∂t ν π t 0 ūxx (x, s) √ t -s ds . (4.6)
Equations (4.5), (4.6) can be analyzed if we apply the Fourier transform in space x → k and the Laplace transform [START_REF] Widder | The Laplace transform[END_REF] in time t → τ reducing these PDEs to a system of two linear equations. This analysis is relatively classical, consequently, we do not enter into algebraic details. The solvability condition for this system gives the required dispersion relation between the frequency τ and the wavenumber k:

(τ + 2νk 2 ) τ 1 + k 2 3 + 2νk 2 - √ τ νk 2 + k 2 1 - ν τ 1 + k 2 3 = 0. (4.7)
The last relation is a polynomial equation of the degree five if we make a change of variables τ = z 2 . There is no general closed-form solution to this equation. However, it can be efficiently solved numerically for any given combination of parameters.

For comparison, the dispersion relation of the free surface Euler equations (the fluid is supposed to be perfect) takes the following classical form [START_REF] Mei | The applied dynamics of ocean surface waves[END_REF]64]:

τ 2 + k tanh k = 0. ( 4.8)
Below we make a comparison between the last relation with the dispersion relation (4.7) to viscous Peregrine equations (4.3), (4.4).

Numerical results and discussion

In this section we assess the performance of the derived above viscous Peregrine system in approximating the dispersive properties of the full Euler equations for a perfect fluid with a free surface. The dispersion relation analysis allows also to study the dissipative properties of the proposed model. For this purpose we consider another quantity of interest which is called the phase velocity defined as:

c(k) := τ (k) k = c r (k) + ic i (k),
where the imaginary part c i (k) is responsible for the wave propagation phenomena and the real part c r (k) describes the dissipative properties of the system (small perturbations are damped if c r < 0 and amplified otherwise).

The wave frequency τ is found for each value of the wavenumber k by solving numerically (Newton-type iterations initialized with relation (4.8)) equation (4.7). In this way we obtained numerical results presented below. For the full Euler equations we use the exact analytical relation (4.8). Finally, the phase speed is obtained after a simple division operation by the wavenumber k. Remark 6. The choice of the eddy viscosity for practical simulations of water waves is not obvious. However, in various experimental and theoretical studies researchers independently come to the conclusion that a value of ν ≈ 10 -3 m 2 /s fits very well available data [60, [START_REF] Bona | An Evaluation of a Model Equation for Water Waves[END_REF][START_REF] Tuck | Sea wave pattern evaluation[END_REF]. The most recent experimental study of plunging breakers confirms this value again [START_REF] Tian | Energy dissipation in two-dimensional unsteady plunging breakers and an eddy viscosity model[END_REF]. Consequently, we retain this value for our numerical illustrations as well. For comparison, the molecular viscosity of the water is of the order ν ≈ 10 -6 m 2 /s which is too small to model the energy dissipation phenomena in a laboratory wave tank, for example. Hence, this molecular viscosity should be replaced by its effective value.

On Figure 2 we show the imaginary part of the phase velocity for Euler and viscous Peregrine equations. For our numerical computations we fixed the viscosity parameter ν to be equal to 10 -3 (see Remark 6). A zoom on long waves region is shown on Figure 2(a) where we can see that the wave propagation is slightly slown down by viscous effects. Mathematically this effect is ascribed to nonlocal terms which are more important in magnitude for small wavenumbers kh. We can conclude that for kh < 2 the viscous Peregrine system provides a reasonnable approximation to full Euler equations. The behaviour of the phase velocity for shorter waves is shown on Figure 2(b). The real part of the wave frequency τ is depicted on Figure 3. On the left image in Figure 3(a) one can see a zoom on the long wave region. The apparent transient behaviour at kh ≈ 1 corresponds to the change between the region where the nonlocal terms are dominant (kh < 1 for the present choice of ν) to the short wave region dominated by local diffusion (τ ≈ -2νk 2 for kh > 1).

Conclusions

In the present study we investigate further the derivation methods of dissipative Boussinesq equations [START_REF] Dutykh | Dissipative Boussinesq equations[END_REF][START_REF] Dutykh | Visco-potential free-surface flows and long wave modelling[END_REF]. First, we recast the visco-potential formulation [START_REF] Liu | Viscous effects on transient long-wave propagation[END_REF][START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF] using the normal velocity at the free surface which is given by the so-called Dirichlet-to-Neumann map in the spirit of the Petrov-Zakharov formulation [START_REF] Petrov | Variational statement of the problem of liquid motion in a container of finite dimensions[END_REF]67]. Then, two completely different approaches have been employed to simplify this complete visco-potential formulation. The first method is called the shape derivative method [40] which expands the D2N map in powers of the nonlinearity parameter ε. The second method is more classical [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF][START_REF] Madsen | Higher-Order Boussinesq-Type Equations for Surface Gravity Waves: Derivation and Analysis[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF] and consists in expanding the velocity potential φ and, consequently, the D2N map in powers of the dispersion (or shallowness) parameter µ 2 . Both methods lead to the same asymptotic expression for the D2N map and thus, to the same dissipative Boussinesq type system. The shape derivative method seems to involve more computational work. However, we prefer to keep both methods in our exposition. The main reason is that in the absence of the Stokes number assumption (ε ∼ µ 2 ), two methods are fundamentally different since they involve expansions in different small parameters (nonlinearity and dispersion). Consequently, the shape derivative method is fully dispersive, while the long wave expansion is fully nonlinear if no additional assumptions are made on the magnitude of the second parameter (dispersion and nonlinearity, respectively). Moreover, we came up with two new Boussinesq type systems, one of them being a viscid generalization of classical Peregrine system [51]. Finally, the dispersion relation properties of the derived system were discussed as well.

In future investigations the decay rate of solutions to new systems at infinite times t → +∞ will be studied. These asymptotic estimations can be supported by several numerical simulations of a solitary wave decay rate. The effect of new viscous terms is to be revealed in future studies. visity to Amiens. The authors would like to thank Hervé Le Meur and Professors Min Chen and Frédéric Dias for helpful discussions on viscous Boussinesq equations.

Figure 1 .

 1 Figure 1. Sketch of the physical domain.

Figure 2 .Figure 3 .

 23 Figure 2. Imaginary part of the phase velocity with zoom on long waves region on the left image. The blue solid line corresponds to the dispersion relation of the full Euler equations while the black dashed line represents viscous Peregrine system.

Acknowledgement D. Dutykh acknowledges the support from French Agence Nationale de la Recherche, project MathOcean (Grant ANR-08-BLAN-0301-01) along with the support from ERC under the research project ERC-2011-AdG 290562-MULTIWAVE. Also D. Dutykh would like to acknowledge the support from the University of Picardie Jules Verne during his

By direct substitution we can check that equation (3.7) and bottom boundary condition (3.9) are satisfied identically. The free surface condition (3.8) poses however some problems, since Φ(x, y = 0) = η(x) φ(0) y (x, 0) = 0. Fortunately, this shortcoming can be easily corrected since the problem (3.7) -(3.9) is linear in velocity potential. Consequently, the ansatz (3.10) can be corrected to satisfy the free surface condition (3.8):

where the function Ψ(x, y) is a solution to the following problem:

The solution to these equations has been constructed in Section 3.1.1 (with viscous parameter ν = 0) and an application of operator R 0 concludes the construction of the shape derivative:

]. Now we are able to obtain the asymptotic expansion of the D2N map (3.1). It can be performed directly in straightened variables (x, y), since one has:

After substituting expansion (3.5) for the velocity potential φ into the rescaled version of D2N map definition (3.11) and truncating the expansion at O(ε + µ 2 ) terms, we obtain the following asymptotic expression which concludes this section:

where we retained only the leading order viscous effects O(ν

2 ) and O(ν).

Long wave expansion method

The developments presented in this section are mainly inspired by the study [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF]. We recall that we expand here the velocity potential in powers of the dispersion parameter µ 2 :

Dissipative Boussinesq equations

In this section we derive a Boussinesq-type evolution system for viscous long water waves using the asymptotic expansion (3.12) of the D2N map derived above. Substituting expression (3.12) into evolution equations (2.8), (2.9) and neglecting higher order terms we obtain the following system of Boussinesq-type equations [START_REF] Dutykh | Viscous potential free-surface flows in a fluid layer of finite depth[END_REF][START_REF] Dutykh | Visco-potential free-surface flows and long wave modelling[END_REF]:

2) where u := ϕ x is the velocity defined at the free surface.

We can further improve dispersive properties of the derived Boussinesq system (4.1), (4.2) if we change the velocity variable to the depth averaged velocity [51], for example.


Lâm Laurent Lim
email: lam.lim@gi-etu.grenoble-inp.fr

Simme Douwe Flapper

Jean-Philippe Gayon
email: jean-philippe.gayon@grenoble-inp.fr

Optimal control of a hybrid production/remanufacturing system using one shared resource

Keywords: remanufacturing, shared server, optimal policy, heuristics, Markov processes, quasi-birth-death processes

Introduction

Reusing product returns

There are several reasons why companies consider reusing products: the growth of environmental consciousness, legal obligation, costs reduction, open potential new markets etc. Returns complicate the inventory control and integrating the return flow of used products into the materials planning is a key issue in reverse logistics.

Hybrid production and remanufacturing systems

In our problem, a company can use both a remanufacturing facility to sell returned products and a classical production facility to satisfy the demand that cannot be fulfilled by returns. Such systems are called hybrid production/remanufacturing systems and have received a growing attention in recent years. To perform production and remanufacturing, a company can use dedicated or shared resources. Sharing resources increases the flexibility of the system but complicates the production/remanufacturing control. The study of [START_REF] Tang | Economic Lot Scheduling Problem with Returns[END_REF] is one example of a concrete application of a hybrid system with shared resource in the automotive industry.

Model formulation

We consider a single-product production-inventory system with independent product returns. Returned items require remanufacturing operations (always successful) before being sellable again. Production and remanufacturing operations are performed by a single shared resource. Switching times and costs are neglected. Customer demand can be fulfilled without distinction by either newly produced or remanufactured products. Two inventories are considered: the remanufacturable inventory composed of product returns, and the serviceable inventory composed of new and remanufactured products. We use a continuous-time model with Markovian assumptions: demands and returns follow independent Poisson processes, production and remanufacturing times are exponentially distributed. We further assume linear holding and backlogging costs. We want to find the optimal policy that minimizes costs (under discounted and average cost criteria).

Structure of the optimal policy

Optimal policy in the general situation

We show that the optimal policy structure seems to be complex in the general case. However, we prove that if the remanufacturing rate is larger than the production one and the remanufacturing holding cost is larger than the serviceable one, then it is optimal to repair returned products as soon as possible: in this case, the optimal policy belongs to the class of policies with priority on returns.

A characterization of the optimal policy with priority on returns

Using Markov decision processes [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF], we characterize the optimal policy for the class of policies where returns have priority over the production. For the discounted cost criteria, the optimal policy is the solution of a fixed-point equation [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]. Based on structural results for queueing systems [START_REF] Koole | Structural results for the control of queueing systems using event-based dynamic[END_REF], we show that the optimal policy is a base-stock policy: the server produces new products if and only if the serviceable inventory is below a certain threshold. Then, this result is extended to the average cost criteria, as explained in [START_REF] Weber | Optimal Control of Service Rates in Networks of Queues[END_REF].

Algorithm to compute the optimal policy parameters and heuristics

Stochastic dynamic programming can be used to compute the optimal policy parameters for both discounted and average cost criteria. However, for policies with priority on returns, we provide a faster method. When considering a base-stock policy, we show that the system evolves according to a quasi-birth-death process [START_REF] Neuts | Matrix Geometric Solutions in Stochastic Models[END_REF]. Because of the special structure of the generator matrix Q associated to this process, we use a matrix-geometric approach [START_REF] Neuts | Matrix Geometric Solutions in Stochastic Models[END_REF] to compute the steadystate probabilities. The advantage of using this method is to reduce the sizes of matrices that are required for computations. Because computations may be time consuming, we provide heuristics for approximating the optimal policy parameter. To derive these heuristics, we modify the exact formula of the optimal base stock level for a make-to-stock queue with direct returns [START_REF] Gayon | A make-to-stock queue with product return[END_REF].


General Integral Formulation for the 3D Thin Shell Modeling

Tung Le-Duc , Gérard Meunier , Olivier Chadebec , Jean-Michel Guichon , and Joao Pedro A. Bastos Grenoble Electrical Engineering Laboratory, University of Grenoble, Grenoble INP/UJF/CNRS UMR 5269, Grenoble, France GRUCAD/EEL/CTC/UFSC, Florianópolis, SC 88040-900, Brazil

In order to compute electromagnetic problems in presence of thin conductive and magnetic shells, a new integral formulation is proposed. Based on the principle of shell element, the formulation is general and enables the modeling of various problems for any skin-depth, avoiding the air region meshing. The formulation is validated thanks to an axisymmetric FEM and compared with another shell formulation implemented in a 3D FEM. Advantages and drawbacks of this new formulation are discussed.

Index Terms-Eddy currents, integral equation method, quasi-static, skin effect, thin shell.

I. INTRODUCTION

T HIN conductive magnetic shells are frequent in electro- magnetic devices (such as planar conductors, electromagnetic shielding…) and the analysis of their eddy-current distributions is a complex problem in engineering.

The geometry of a thin region is characterized by a high ratio between length and thickness. Thus, using a volume mesh for the thin region leads to a large number of elements. Furthermore, when the frequency is high the skin depth becomes much thinner than the thickness . It increases the size of the needed mesh if this effect is properly modeled. Shell element formulations have been developed to overcome such difficulties, e.g., with the boundary element method (BEM) [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], with the finite-element method [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF]- [START_REF] Thomas | A three-dimensional electromagnetic shell finite element for coupled vector-scalar potential formulations[END_REF] and with an integral method recently proposed by the authors of this paper [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF].

In [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF], a general shell element formulation for thin conductive regions modeling has been proposed ( or or

). As in [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF]- [START_REF] Thomas | A three-dimensional electromagnetic shell finite element for coupled vector-scalar potential formulations[END_REF], this formulation takes into account the field variation through the depth due to the skin effect. The shell is modeled by an integral formulation avoiding the air region meshing. Based on a simple surface discretization of the shell averaged surface, the number of unknowns is considerably reduced, since its depth is not meshed. However, this already proposed shell element formulation was only suitable for conductive materials and could not be applied to shells with both magnetic and conductive properties.

This paper presents a new integral formulation which allows the modeling of magnetic and conductive thin shells in the general case ( or or ). This method is a natural extension of the method proposed in [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF]. Moreover, we will see that an additional advantage is its implementation simplicity avoiding the complex computation of volume singular integrals as in [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF]. Let use notice that the formulation is currently limited to linear cases but this assumption is often reliable for shielding effects computations.

In the second part of the paper, the integral equation governing the thin shell regions is presented. In the third part, this equation will be coupled with an integral formulation. Finally, a numerical example is proposed in the last part. Results obtained with our formulation are compared with results obtained with the Finite Element Method (FEM).

II. THIN SHELL EQUATION

Let be a thin magnetic conductive shell with thickness (conductivity and linear permeability ) and , boundaries of the shell with the air region (see Fig. 1.). We consider that a given skin depth is associated to the shell. The field variation of the tangential component across the thickness of shell can be approximated by the analytical solution of the problem for an infinite plane [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF]:

(1) where , and are the tangential values of the magnetic field on both sides of the shell. Applying Galerkin's method on the Maxell-Faraday equation for the side "1" of the shell, we have [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF], [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF]:

(2)
where , ; w is a set of nodal surface weighting functions; is the normal vector corresponding to the side "1" of the shell and is the line region delimiting the surface. The symmetrical equation corresponding to the other side of the shell is obtained by simply exchanging the indices "1" and "2":

(3) 0018-9464/$31.00 © 2013 IEEE where is the normal vector corresponding the side "2" of the shell and is the line region delimiting the surface. Equations ( 2) and (3) represent the electromagnetic behavior of the conductive shell itself. In order to take into account the surrounding air region, several authors choose the Finite Element Method [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF]- [START_REF] Thomas | A three-dimensional electromagnetic shell finite element for coupled vector-scalar potential formulations[END_REF]. Thus, the air region has to be also discretized.

In [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF], we developed an integral approach based on the coupling of ( 2) and (3) with a volume integral equation. The obtained formulation presents an interesting compactness treating different air volumes as a single physical region.

However, this formulation still suffered from drawbacks. The first one is its application only on non magnetic materials. Another difficulty is the computation of volume integrals which can become very complicated to accurately consider the exponential decreasing of the current across the thickness, especially when the influence of an element on itself is focused. The following section proposes an approach for solving these issues.

III. INTEGRAL FORMULATION

Equations ( 2) and ( 3) are written on the averaged surface of the thin region. Subtracting (3) and (2) leads to:

(4) where is the averaged induction and is the normal vector corresponding to the surface (see. Fig. 1

.).

A reduced magnetic scalar potential is now introduced on each side of the shell. Because of the jump of the magnetic field tangential component through the thickness, a double layer surface element is necessary [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF]. The expressions related to the tangential magnetic fields on both sides and outside of the shell are: [START_REF] Krähenbühl | Homogenization of lamination stacks in linear magnetodynamics[END_REF] where and are the tangential source fields generated by inductors at each side of the shell and are the both associated scalar reduced magnetic potentials.

Using (4) and ( 5) and assuming small variations of through the shell, we get: [START_REF] Dang | Subproblem approach for thin shell dual finite element formulations[END_REF] where is the scalar magnetic potential discontinuity, i.e. the potential jump of the double layer.

The magnetization law of the material is linear and associated to a constant relative permeability . The (6) can thus be rewritten as: [START_REF] Guérin | 3-D magnetic scalar potential finite element formulation for conducting shells coupled with an external circuit[END_REF] where is the averaged magnetization. On the side "1" of the shell, the total magnetic field is the sum of the inductor field , the field created by magnetization and the field created by eddy currents flowing in the shell. Hence an integral equation linking the local field to the magnetization and the eddy current can be written [START_REF] Chadebec | A review of magnetostatic moment method[END_REF]:

(8) where is the vector linking the integration point on to the point where the field is expressed and is the current density.

The integration of tangential magnetization through the depth can be expressed by [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF]: [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF] where . The equivalent shell current distribution can be also expressed by [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF]: [START_REF] Chadebec | A review of magnetostatic moment method[END_REF] Using ( 9) and ( 10), [START_REF] Thomas | A three-dimensional electromagnetic shell finite element for coupled vector-scalar potential formulations[END_REF] becomes: (11) The equation on the side "2" of shell is obtained by permuting the subscripts "1" and "2" as [START_REF] Abakar | 3D modeling of shielding structures made by conductors and thin plates[END_REF] If we consider that has no variation through the thickness of the shell, we obtain: [START_REF] Nguyen | Ships magnetic anomaly computation with integral equation and fast multipole method[END_REF] Equations ( 7) and [START_REF] Nguyen | Ships magnetic anomaly computation with integral equation and fast multipole method[END_REF] have to be discretized. The easiest way is to mesh the averaged surface , located halfway between both boundaries of the shells (see Fig. 1.), into triangular elements. They are associated to a uniform tangential component of the eddy current and uniform magnetization (meaning 1-order shape functions for the potential jump and 0-order shape functions for the magnetization). The algebraic linear system obtained has complex unknowns (a single complex magnetization vector per element and a single complex magnetic scalar potential per node; the mesh being composed of nodes). Let us notice that shape functions with higher orders could be also used.

Equation ( 7) can be written with the following matrix form: [START_REF] Buchau | Parallelization of a fast multipole boundary element method with cluster OpenMP[END_REF] where: is a complex vector of dimension , is a complex vector of dimension .

is a sparse matrix, expressed as follows:

(15) is a sparse matrix and can be written as:

(16)

A matrix system representing ( 13) can also be written thanks to a collocation approach or a most expensive Galerkin projection. Here a simple collocation at elements centroids can be used to get [START_REF] Chadebec | A review of magnetostatic moment method[END_REF]:

(17)
where represents the identity matrix; is the source fields projected at element centroids, is a dense matrix and is a also dense matrix. The coefficients of these matrices are expressed as follows:

(18) where denotes the classical scalar product operator and is the vector basis of cell element i and is the vector linking the integration point on to the centroid of element i.

(19) Matrix terms ( 18) and ( 19) present singularities classically obtained with integral formulations so particular treatments are needed especially for the computation of the element influence on itself. Coefficients (18) are computed thanks an equivalent magnetic charge distribution located on edge of elements like proposed in [START_REF] Chadebec | A review of magnetostatic moment method[END_REF]. The integral is computed numerically. This well-known approach avoids the superposition of Gauss points and points where the field has to be computed and thus suppress the singularity. For (19), we consider that the averaged internal field created by a conductor carrying a uniform current on itself is globally very low. Thus, the coefficient C(i,i) is imposed to zeros while other coefficients are computed with standard numerical integration. Let us observe that in the proposed formulation, integrals are computed only thanks to a surface numerical integration. This is an important difference with [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF] where accurate and uneasy volume integration across the thickness was necessary. The method is then considerably simpler and reliable as shown in the following part.

IV. VALIDATION

The formulation has been implemented for a 3D geometry and has been validated by two test problems. The first case is a thin magnetic conductive disk placed in a magnetic field created by a current loop (Fig. 2.). This example is modeled by three different methods on a standard computer (Intel CPU 5160@3 GHz; 16 GB). The first one is the axisymmetric FEM using FLUX software [11] (100,000 elements, time of resolution 70 s). The second one is a shell element formulation implemented in 3D FEM code [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF], [START_REF] Guérin | 3-D magnetic scalar potential finite element formulation for conducting shells coupled with an external circuit[END_REF], [11] with 300,000 volume elements. The last one is the proposed integral method, the disk being meshed in 800 surface elements. We focus on the magnetic field computed in the air region for different frequencies at the point A (0.25, 0, 0.1) (Fig. 3).

If we consider the axisymmetric FEM method as our reference, the integral method leads to an error of 1.6% for the magnetic field at and 0.2% at (Fig. 4). The integral formulation leads to more accurate results than the same implementation of the shell element but considering the air region treated with the 3D FEM. Let us notice that computation times are similar (around 300 s for 3D FEM and 240 s for the integral formulation).

The second test problem is the modeling of a real device proposed by EDF (Electricité de France) [START_REF] Abakar | 3D modeling of shielding structures made by conductors and thin plates[END_REF]. In order to limit short-circuit currents, the power station "Folies" is equipped with three reactances in series with a power transformer secondary. The current in each reactance is 1000 A phase-shifted with 120 degree and the frequency is 50 Hz. However, in nominal conditions, these three reactances create a leakage induction in the neighborhood. In order to minimize this electromagnetic disturbance, an electromagnetic shielding is added between the sources and the area to be protected. The skin depth is smaller than the shielding thickness. The Fig. 6 shows the current distribution computed by the proposed integral method and a more classical FEM 3D. Results provided by our coupling are very encouraging, the convergence being reached with few elements (less than 1000). This is the strong point of integral formulations which lead to accurate field computation in air regions with very low number of unknowns. However, some matrices obtained are fully populated and compression algorithms must be used if there are a large numbers of elements. In [START_REF] Nguyen | Ships magnetic anomaly computation with integral equation and fast multipole method[END_REF], [START_REF] Buchau | Parallelization of a fast multipole boundary element method with cluster OpenMP[END_REF], authors chose to use the well-known "Fast Multipole Method" (FMM) to compress matrices and .

V. CONCLUSION

In this paper, we have presented a new integral formulation using shell elements in order to model thin conductive and magnetic regions. The formulation is general and various skin depths ( or or ) can be accurately modeled. It main advantage is its simplicity of implementation and its reliability. It drawbacks in comparison with our previous formulation proposed in [START_REF] Le-Duc | A new integral formulation for eddy current computation in thin conductive shells[END_REF] is that results are slightly less accurate for the field near the shell than in and we cannot compute the losses. However, the formulation remains quite efficient if problems where the magnetic field in the air and global shell current distribution must be evaluated. According to us, the coupling of this formulation with compression algorithms does not suffer from any difficulty.

Fig. 1 .

 1 Fig. 1. Thin conductive and magnetic region and associated notations.

Fig. 2 .

 2 Fig. 2. Thin magnetic conductive disk ( , , , ,).

Fig. 3 .

 3 Fig. 3. Magnetic field in the point A (0.25, 0, 0.1), calculated by 3 methods versus .

Fig. 4 .

 4 Fig. 4. Difference in % with the axisymmetric FEM solution for magnetic field computation versus .

Fig. 5 .

 5 Fig. 5. Three reactances and electromagnetic shielding.

Fig. 6 .

 6 Fig. 6. Surface distribution of current (A/m) in the plaque for test case "Folies". a) FEM 3D with shell elements; b) Integral Method.

ACKNOWLEDGMENT

This work was supported in part by the French National Research Agency (ANR/project E-CEM) and the Associated International Laboratory "James Clerk Maxwell" between France and Brazil (LIA 917-CNRS/CNPQ).