Two-Scale Simulation of Acoustic Equations

(1)

c

TU Delft, The Netherlands, 2006

TWO-SCALE SIMULATION OF ACOUSTIC EQUATIONS

S´ebastien Jund∗ and Eric Sonnendr¨ucker∗

∗ _{Institut de Recherche Mathématique Avancée} Université Louis Pasteur

7, rue Ren´e Descartes 67084 Strasbourg Cedex, France.

e-mail : jund@math.u-strasbg.fr, sonnen@math.u-strasbg.fr

Key words: acoustic equations, domain decomposition, multi-scale resolution, Raviart-Thomas finite-elements.

Abstract. We develop a numerical method for solving the linearized acoustic equations on a grid involving zones with cells of very different sizes, in order for example to compute sources coming from the flow on a fine grid. The method is based on domain decomposition techniques which lead us to introduce two auxiliary problems and show theoretically how they allow us to calculate the solution of the initial problem. These two auxiliary problems are discretized using Raviart-Thomas finite elements [5, 3] on two different scales which introduce some errors that we correct by setting to zero an operator we know has to be zero in the theoretical study.

1 INTRODUCTION

An important issue when dealing with near to mid field acoustics is to obtain a good approximation of the sources from the numerical computation of the flow. The space scale involved in the flow solver are a lot smaller than those needed for the acoustics solver. Therefore it is natural to perform the acoustics computation on a mesh with two zones involving largely different cells sizes, small cells around the flow where the sources are computed and larger cells away from the flow where only the larger wavelength acoustic waves need to be computed. A naive coupling of the discretizations on the two zones can not only lead to inaccurate results, but often generate numerical instabilities due to short wavelength modes which are not resolved on the coarse scale bouncing back and forth on the fine grid.

(2)

The paper is organized as follows: First we introduce a continuous problem based on two auxiliary problems with an artificial interface separating the two regions of interest. The continuous problem enables us to find an operator on the interface on which to iterate to obtain the continuity of the normal component. Then we introduce the Raviart-Thomas finite element that we shall use and the discretized problem. In the following section we present the algorithm that is used to compute an approximation of the solution in the steady state case and after that we extend the method to the time-dependent case before presenting numerical results.

2 CONTINUOUS PROBLEM

Let Ω ⊂ IR2 be an open, bounded, regular domain of IR2 with boundary Γ = ∂Ω and let ω ⊂ Ω be another open regular domain with boundary γ. We denote by −→n the unit normal vector outward to the domain ω on γ.

n Ȧ ī Ȗ Figure 1: Domain.

We consider the following problem:

U − c∇(∇ · U) = f in Ω, (1)

U · −→n = 0 on Γ. (2)

where f ∈ (L2_(Ω))2 _{is the sum of two functions f}

1, f2 ∈ (L2(Ω))2, f = f1 + f2 with

supp(f2) ⊂ ω, and c is a positive constant.

(3)

         V − c∇(∇ · V ) = f1 in Ω, V · −→n = 0 on Γ, [V · −→n]= 0 on γ, [(∇ · V )−→n]= −λ on γ (3) and ₍ W− c∇(∇ · W ) = f1+ f2 in ω, W· −→n = V · −→n on γ (4)

where, [ψ] = ψ+ _{− ψ}− _{denotes the jump of ψ on γ , ψ}+ _{and ψ}− _{are respectively the}

restriction of ψ to Ω\ω and ω.

Since supp(f2) ⊂ ω, and with regard to (3) and (4), it is sufficient to remark that

Proposition 2.1 If λ is chosen such that

λ = (∇ · V−)−→n − (∇ · W )−→n , (5) then the solution U of the initial problem (1) is given by

U = V+.χ(Ω\ω) + W.χ(ω) (6)

Proof. The proof is based on the following : Take U verifying (6). Then ∀φ ∈ H0(−→div, Ω),

Z ΩU.φ dX+ c Z Ω(∇.U)(∇.φ) dX = Z ΩU.φ dX + c Z Ω\ω(∇.V +_{)(∇.φ) dX + c}Z ω(∇.W )(∇.φ) dX (7) where Z Ω\ω(∇.V +_{)(∇.φ) dX = −}Z Ω\ω(∇(∇.V +_{)).φ dX +}Z Γ∪γ((∇.V +₎₋→_m_{).φ dσ} (8) and Z ω(∇.W )(∇.φ) dX = − Z ω(∇(∇.W )).φ dX + Z γ((∇.W )− →_n_{).φ dσ} ₍₉₎

where −→m is the outbound normal on ∂(Ω\ω), within particular −→m = −−→n on γ. Hence combining (8) and (9), equation (7) becomes:

(4)

Using (3) and (5) we have λ:= (∇.V−)−→n − (∇.V+)−→n = (∇.V−)−→n − (∇.W )−→n , hence Z γ{((∇.W )− →_n_{) − ((∇.V}+₎₋→_n_{)}.φ dσ = 0.} Finally, we obtain Z ΩU.φ dX + c Z Ω(∇.U)(∇.φ) dX = Z Ω(f1+ f2).φ dX.

Therefore, defining the operator

T λ= ((∇.W )−→n) − ((∇.V+)−→n), it is sufficient to find λ such that T λ = 0. We can also prove that Lemma 2.2 The operator T is given by

T λ= λ + ((∇.W )−→n) − ((∇.V+)−→n) (11) where V and W satisfy respectively

( V − c∇(∇.V ) = f1 in Ω, V .−→n = 0 on Γ, (12) and ₍ W − c∇(∇.W ) = f1+ f2 in ω, W .−→n = V .−→n on γ. (13)

Proof. We give a sketch of the proof. We introduce Wf = W − V− _and _Wf _{= W − V}−

and notice that they both verify

(

E− c∇(∇.E) = f2 in ω,

E.−→n = 0 on γ. (14)

The uniqueness of the solution gives (∇.Wf)−→n = (∇.Wf)−→n which means that

(∇.W )−→n − (∇.V−)−→n = (∇.W )−→n − (∇.V−)−→n .

Then the result follows by adding and subtracting (∇.V+₎₋→_n _{in the left hand side and}

(5)

(∇.V+)−→n − (∇.V−)−→n = −λ, and

((∇.W )−→n) − ((∇.V+₎₋→_n_{) = T λ.}

Thus, it is trivial that T λ = 0 is equivalent to λ = ((∇.V+)−→n) − ((∇.W )−→n).

But if this is true when we consider the continuous problems (12) and (13), it is sure that solving the corresponding discretized problems and defining

λH = ((∇.V +

H)−→n) − ((∇.Wh)−→n) (15)

introduces errors that we will have to correct. But still, the approach is quite clear at this moment : by solving the discretized problems of (12) and (13) we define λH =

((∇.V+H)−→n) − ((∇.Wh)−→n) and then we can solve the discretized problems of (3) and (4)

and finally define the solution of the initial problem (1) as in (6).

3 NUMERICAL APPROXIMATION

We propose to use Raviart-Thomas finite-elements of order one, two and three. We recall that the ith _{order finite-element is given by the following triplet (K, P, Σ) where}

• K denotes a quadrilateral, • P = Qi,i−1

Qi_−1,i

!

denotes a polynomial space where Qm,n = hxiyj,0 ≤ i ≤ m, 0 ≤ j ≤ ni,

• Σ = {σ_Z kl,1 ≤ l ≤ 2(i + 1), 0 ≤ k ≤ i − 1} is the set of linear forms σkl(f ) = Γl

Lkl(f.τl)dσ, where {Γl}1≤l≤2(i+1) denotes a regular subdivision of the

quadrilat-eral K (two times i + 1 lines, including the two times two opposite edges of K), τl

is a unit tangential vector associated to Γl and Lkl is the kth normalized Legendre

polynomial calculated on Γl.

We give for example the four basis functions (i.e. the four polynomial functions Pi

in P which satisfies σi(Pj) = δij ∀(i, j) ∈ {1..4}2) associated to a reference quadrilateral

[0, h] × [0, h] for the first order elements: P1 = 1 √ h(1 − y h) 0 ! , P2 = y h√h 0 ! , P3 = 0 1 √ h(1 − x h) ! , P4 = 0 x h√h ! , and the twelve basis functions for the second order elements (the six first for L0,l and the

(6)

P1 = (h2 −3yh+2y2₎ √ hh2 0 ! , P2 = −4y(−h+y)_√ hh2 0 ! , P3 = y_(−h+2y) √ hh2 0 ! , P4 = 0 h2 −3xh+2x2 √ hh2 ! , P5 = 0 −4x(−h+x)_√ hh2 ! , P6 = 0 x(−h+2x)_√ hh2 ! , P7 =   √ 3(−h3 +2xh2 +3yh2 −6xyh−2y2 h+4xy2 √ hh3 0  , P8 =   −4 √ 3y(h2 −2xh−yh+2xy) √ hh3 0  , P9 =   √ 3y(h2 −2xh−2yh+4xy)_√ hh3 0  , P10 =   √ 0 3(−h3_+3xh2_+2yh2 −6xyh−2x2_h+4x2_y) √ hh3  , P11 =   ₍₋₄√ 0 3x(h2 −2yh−xh+2xy) √ hh3  _, _P₁₂ ₌   √ 0 3x(h2 −2yh−2xh+4xy)_√ hh3  _, Γ Γ Γ Γ 1 2 3 4 n n n n 1 2 3 4 K_ref Γ Γ Γ Γ 1 2 3 4 n n n n1 6 3 4 K ref Γ Γ5 6 n n 2 5

Figure 2: The reference quadrilateral, its four and twelve (six times two per Γl) degrees of freedom for the first and second order elements.

Once the two meshes are fixed, we define two spaces of finite-element approximation PH and Ph respectively associated to the coarse mesh over the whole domain Ω and to

the fine mesh over ω.

Assuming that the boundary γ coincides with the union of the edges of quadrilaterals defining the coarse mesh, we also define a finite-element approximation space ∆H which

corresponds to the normal trace on γ of any function in PH. More precisely, using the ith

order finite elements implies that ∆H is a 1D IPi−1 finite-element space on γ discretized

(7)

moment integrals over each edge. It is in this last space that λ has to be searched. The idea is to define an operator TH : ∆H → ∆H approximating appropriately T λ and

then to search for λH such that THλH = 0 in order to correct the λH defined in (15).

Thus, an approximation of the solution U of the problem (1) is given by

UHh = VH+.χ(Ω\ω) + Wh.χ(ω) (16)

where VH+ is the restriction on Ω\ω of the solution of the following equation

Z ΩVH.φ dX+ c Z Ω(∇.VH)(∇.φ) dX = Z Ωf1.φ dX+ c Z γ λH.φ dσ ∀φ ∈ PH, (17)

and Wh is the solution of

Z ω Wh.φ dX + c Z ω(∇.Wh)(∇.φ) dX = Z ω(f1+ f2).φ dX ∀φ ∈ Ph (18)

with boundary conditions Wh.−→n = VH.−→n on γ.

One can see that the normal restriction ((∇.V_H+)−→n) on γ is by definition in ∆H,

but it is not the case of ((∇.Wh)−→n). This means that we can not directly define

TH = ((∇.Wh)−→n) − ((∇.VH+)−→n).

Then, the most natural way to construct a consistant approximation of T λ is to find δH ∈ ∆H such that Z γ δH.Φ dσ = Z γ((∇.Wh)− →_n_{) − ((∇.V}+ H)−→n).Φ dσ ∀Φ ∈ ∆H, (19)

and then define THλH = δH. Unfortunately this simple definition does not work well in

practice.

A more sophisticated way to define TH is the following : rather than integrating directly

Z

γ((∇.Wh)−

→_n_{).Φ dσ and} Z

γ((∇.V +

H)−→n).Φ dσ in equation (19) which will be a very coarse

approximation because we only use information on the boundary γ, we use the fact that if V − c∇(∇.V ) = f1 in Ω, (20) then Z Ω\ωV + H.φ dX + c Z Ω\ω(∇.V + H)(∇.φ) dX = Z Ω\ωf1.φ dX− c Z γ((∇.V + H)−→n).φ dσ ∀φ ∈ PH, (21) and combining this with equation (17) we obtain

(8)

for any Φ in ∆H, where ˜Φ is an extension of Φ in PH.

In the same way we obtain another expression of c Z γ((∇.Wh)− →_n_{).Φ dσ =}Z ωWh. ˜Φ dX + c Z ω(∇.Wh)(∇. ˜Φ) dX − Z ω(f1+ f2). ˜Φ dX. (23) We finally propose to calculate THλH = δH where δH verifies : ∀Φ ∈ ∆H,

Z γδH.Φ dσ = Z γλH.Φ dσ + 1 c Z ωWh. ˜Φ dX + Z ω(∇.Wh)(∇. ˜Φ) dX − 1 c Z ω(f1+ f2). ˜Φ dX − 1 c Z ω V_H−. ˜Φ dX + Z ω(∇.V − H)(∇. ˜Φ) dX − 1 c Z ω f1. ˜Φ dX . (24) 4 SOLUTION ALGORITHM

As said at the end of section 2 the use of numerical resolution introduces some errors that have to be corrected. This means that defining λH by (15) does not ensure that

THλH = 0. This is the reason why we decide to apply a fixed point method on (I − TH)

as follows :

1. Define λH = 0 identically as an initialisation.

2. Solve _         VH − c∇(∇.VH) = f1 in Ω, VH.−→n = 0 on Γ, [VH.−→n]= 0 on γ, [(∇.VH)−→n]= −λH on γ (25)

This resolution gives VH.−→n on γ.

3. Solve ₍

Wh − c∇(∇.Wh) = f1+ f2 in ω,

Wh.−→n = VH.−→n on γ.

(26)

4. Define THλH = δH using equation (24).

5. Finally define λH = λH − THλH.

(9)

One can see that at the first iteration, with λH = 0, solving (25) and (26), defining

THλH = δH and redefining λH = λH − THλH corresponds exactly to what we said at

the end of section 2, i.e solving the discretized problems of (12) and (13) and defining λH as in (15), so that defining UHh as in (16) with VH and Wh calculated at the second

iteration should be a good approximation of the solution of the initial problem. The further iterations correspond only to the correction of λH which lead to THλH = 0.

5 TIME DEPENDENT PROBLEM

We now consider the following time dependent problem:

(

∂_t2U − ∇(∇ · U) = f in Ω × [0, T ],

U · −→n = 0 on Γ × [0, T ]. (27)

where T is a given final time, and f is also supposed to be the sum of two functions f1, f2 ∈ (L2(Ω) × [0, T ])2 with supp(f2(., t)) ⊂ ω at each time t ∈ [0, T ]. We introduce the

two auxiliary problems:

         ∂_t2V − ∇(∇ · V ) = f1 in Ω × [0, T ], V · −→n = 0 on Γ × [0, T ], [V · −→n]= 0 on γ × [0, T ], [(∇ · V )−→n]= −λ on γ × [0, T ] (28) and ₍ ∂_t2W − ∇(∇ · W ) = f1+ f2 in ω × [0, T ], W · −→n = V · −→n on γ × [0, T ] (29)

After fixing a time step ∆t, introducing the time discretization tn = n∆t, and discretizing

the second derivative in time by the standard centered second order scheme, these two problems become:          Vn+1− ∆2 t∇(∇ · Vn+1) = ∆2tf1n+1+ 2Vn− Vn−1 in Ω, Vn+1· −→n = 0 on Γ, [Vn+1_{· −}→_n_{]= 0 on γ,} [(∇ · Vn+1₎₋→_n_{]= −λ}n+1 _{on γ} (30) and ₍ Wn+1_{− ∆}2 t∇(∇ · Wn+1) = ∆2t(f1n+1+ f2n+1) + 2Wn− Wn−1 in ω, Wn+1· −→n = Vn+1_{· −}→_n _{on γ} (31)

which have to be solved at each time step n∆t. Defining

(10)

we automatically have ˜

f2 = ∆2tf n+1

2 + 2(Wn− Vn) − (Wn−1− Vn−1).

Then, and most importantly, we have to give an interpretation of the source terms ˜

f1 and ˜f2. Indeed, ˜f1 and ˜f2 make appear terms like Vn and Vn−1 which have to be

thought as the solution of the initial problem (27) on the coarse grid respectively at time n∆t and (n − 1)∆t which are not the solutions of problem (30) : we recall that problem

(30) only gives the solution of the initial problem (27) in Ω\ω and that this solution is given by Wn _{(solution of problem (31)) in ω so that the terms V}n _{and V}n−1 _{in ˜}_f

1 and ˜f2

must be defined as Vn _{and V}n₋₁ _{(solution of problem (30) at time n∆}

t and (n − 1)∆t) in

Ω\ω and as the projection on the coarse space discretization of Wn _{and W}n−1 _(solution

of problem (31) at time n∆t and (n − 1)∆t) in ω.

Taking this into account, we can solve the problems (30) and (31) at each time step using the two scale method previously described in this paper.

6 NUMERICAL RESULTS

In this section we give numerical results which highlight the efficiency of the two-scale method. In the following, the solution of the initial problem using our two-scale method is called numerical solution, and the solution of the initial problem calculated directly on the extension of the fine grid over the whole domain is refered to be the reference solution. As we are only interested in the loss of precision when we use our method rather than a direct resolution on a uniformly refined mesh, we will only consider the error between these two solutions; which will be given in L2_{(Ω) and H(}−→_{div, Ω) norms. We will split}

these errors into interior errors and exterior errors which respectively correspond to the errors in ω and Ω\ω. The first two test cases highlight the efficiency of the method on the stationary problem and the third and fourth test cases deal with the time dependent problem.

6.1 Test case 1.

In this test case, we consider Ω = [−π 2, π 2] × [− π 2, π 2], ω = [− π 6, π 6] × [− π 6, π 6], f1 = 2 cos(y) 2 cos(x) !

as background source and f2 = 10 cos(9x) cos(9y)χ(ω) 1

1

!

.

We give in Figure 3 the first field of the reference solution, the numerical solution af-ter one iaf-teration and the numerical solution afaf-ter two iaf-terations. One can see that the algorithm produces exactly what we expected : what we see after the first iteration is not supposed to be an approximation of the solution but something that allows us to cal-culate λH which is an approximation of λ. Then, at the second iteration, the algorithm

(11)

0 X Y Z 0 X Y Z 0 X Y Z

Figure 3: First field of the reference solution, the numerical solution after one iteration and the numerical solution after two iterations for the third test case.

6.2 Test case 2.

For our second test case we consider Ω = [−π 2, π 2]×[− π 2, π 2], ω = [− 13π 54, 13π 54]×[− 13π 54 , 13π 54], f1 = 2 cos(y)_{2 cos(x)} !

as background source and f2 = 10 cos(9x) cos(9y)χ([−π₆,π₆]×[−π₆,π₆]) 1₁

!

. One can see that the only difference between this test case and the third test case is that we take ω slightly greater than supp(f2). We give in Figure 4 the first field of the

reference solution and the numerical solution after two iterations. We give in Table 3 the errors calculated on a 27 × 27 coarse grid and a 3 time finer fine grid.

6.3 Test case 3.

For this first time dependent test case we take Ω = [−12, 12] × [−12, 12], ω = [−4, 4] × [−4, 4], f1 = f2 = 0

0

!

. We also take U(x, y, 0) =

(12)

First order finite elements

Iteration 1 2 5 10

L2 _error _{0.218627E+01 0.156232E+01 0.155576E+01 0.155576E+01}

L2 exterior error 0.150718E+01 0.115932E+01 0.117306E+01 0.117308E+01 L2 interior error 0.158372E+01 0.104729E+01 0.102192E+01 0.102190E+01 H(div) error 0.221645E+01 0.159050E+01 0.158382E+01 0.158382E+01 H(div) exterior error 0.154440E+01 0.119459E+01 0.120785E+01 0.120787E+01 H(div) interior error 0.158980E+01 0.105007E+01 0.102450E+01 0.102448E+01

Second order finite elements

Iteration 1 2 5 10

L2 _error _{0.238895E+01 0.596396E+00 0.593638E+00 0.593639E+00}

Third order finite elements

Iteration 1 2 5 10

L2 error 0.239262E+01 0.182636E+00 0.182854E+00 0.182855E+00 L2 exterior error 0.165784E+01 0.149684E+00 0.150459E+00 0.150459E+00 L2 _{interior error} _{0.172516E+01 0.104646E+00 0.103913E+00 0.103913E+00}

H(div) error 0.240816E+01 0.183206E+00 0.183419E+00 0.183420E+00 H(div) exterior error 0.167243E+01 0.150357E+00 0.151127E+00 0.151127E+00 H(div) interior error 0.173269E+01 0.104677E+00 0.103939E+00 0.103939E+00

Table 1: Numerical errors on a 9 × 9 grid for the first test case.

0 0

!

, where η = 10 and ε = 0.5.

We give in Figure 5 the first field of the reference solution and the numerical solution after forty-five time iterations (i.e. at time t=6.25). We give in Table 4 the errors calcu-lated on a 27 × 27 coarse grid and a 2 time finer fine grid at time t=4, t=8 and t=12.

6.4 Test case 4.

(13)

Iteration 1 2 5 10

L2 _error _{0.236381E+01 0.490444E+00 0.485887E+00 0.485887E+00}

Iteration 1 2 5 10

L2 _error _0.238689E+01 _0.708898E-01 _0.707343E-01 _0.707344E-01

L2 exterior error 0.165395E+01 0.649181E-01 0.650094E-01 0.650095E-01 L2 interior error 0.172096E+01 0.284781E-01 0.278769E-01 0.278767E-01 H(div) error 0.240243E+01 0.712253E-01 0.710700E-01 0.710700E-01 H(div) exterior error 0.166854E+01 0.652830E-01 0.653735E-01 0.653736E-01 H(div) interior error 0.172848E+01 0.284811E-01 0.278791E-01 0.278789E-01

Iteration 1 2 5 10

L2 error 0.238695E+01 0.819654E-02 0.821060E-02 0.821062E-02 L2 exterior error 0.165399E+01 0.760862E-02 0.762701E-02 0.762702E-02 L2 _{interior error} _0.172100E+01 _0.304830E-02 _0.304017E-02 _0.304020E-02

H(div) error 0.240249E+01 0.822539E-02 0.823938E-02 0.823940E-02 H(div) exterior error 0.166858E+01 0.763965E-02 0.765796E-02 0.765797E-02 H(div) interior error 0.172852E+01 0.304841E-02 0.304024E-02 0.304027E-02

Table 2: Numerical errors on a 27 × 27 grid for the first test case.

∂tU(x, y, 0) =

0 0

!

, where η = 10 and ε = 0.5.

We give in Figure 6 the first field of the reference solution and the numerical solution at time t=12, and in Table 5 the errors calculated on a 27 × 27 coarse grid and a 2 time finer fine grid at time t=4, t=8 and t=12.

6.5 Interpretation

(14)

X Y Z X Y Z

Figure 4: First field of the reference solution and the numerical solution after two iterations for the second test case.

few iterations; the interior error, which is considerably smaller, really decreases iteration after iteration so that we can conclude that the only error we add by using our two-scale decomposition comes from the fact that we use a coarser finite-element space in Ω\ω which is unable to fit the solution as well as the finer space defined on the extension of the fine grid over the whole domain Ω. Indeed, it is not only necessary to ensure that the finite-element space is fine enough to capture the source term (what we did by defining Ph

in ω) but also that the finite-element space is fine enough to represent the solution over the whole domain Ω, which is not trivial as we use low order finite-elements, knowing that even if a peaking source term is located in ω, it also perturbates, in general, the solution in Ω\ω. We highlight this consideration by simply taking ω slightly larger than supp(f2),

as in the second test case, and see how much better the results are, compared to the first test case.

(15)

Iteration 1 2 5 10

L2 _error _0.469249E+00 _0.685708E-01 _0.680016E-01 _0.680016E-01

L2 exterior error 0.312402E+00 0.626100E-01 0.629453E-01 0.629453E-01 L2 interior error 0.350143E+00 0.279633E-01 0.257316E-01 0.257315E-01 H(div) error 0.503046E+00 0.105062E+00 0.104682E+00 0.104682E+00 H(div) exterior error 0.350356E+00 0.101262E+00 0.101467E+00 0.101467E+00 H(div) interior error 0.360979E+00 0.279987E-01 0.257457E-01 0.257456E-01

Iteration 1 2 5 10

L2 _error _0.477248E+00 _0.821211E-02 _0.820462E-02 _0.820462E-02

L2 exterior error 0.318481E+00 0.765013E-02 0.765598E-02 0.765600E-02 L2 interior error 0.355437E+00 0.298568E-02 0.294986E-02 0.294984E-02 H(div) error 0.504537E+00 0.839975E-02 0.839238E-02 0.839239E-02 H(div) exterior error 0.346950E+00 0.785110E-02 0.785678E-02 0.785679E-02 H(div) interior error 0.366311E+00 0.298599E-02 0.295009E-02 0.295007E-02

Iteration 1 2 5 10

L2 error 0.477250E+00 0.740001E-03 0.742053E-03 0.742056E-03 L2 exterior error 0.318483E+00 0.695763E-03 0.697112E-03 0.697114E-03 L2 _{interior error} _0.355439E+00 _0.252023E-03 _0.254318E-03 _0.254322E-03

H(div) error 0.504537E+00 0.744640E-03 0.746679E-03 0.746682E-03 H(div) exterior error 0.346948E+00 0.700692E-03 0.702032E-03 0.702033E-03 H(div) interior error 0.366312E+00 0.252032E-03 0.254324E-03 0.254328E-03

Table 3: Numerical errors on a 27 × 27 grid for the second test case.

CONCLUSION AND PERSPECTIVES

We have developed a numerical method based on a domain decompostion idea which can handle well two regions of the computation domain which are discretized with different meshes using Raviart-Thomas finite elements. Note that care must be taken in order to define the numerical interface operator TH. A simple definition using only interface

information like (19) does not work well in practise. It is necessary to use the more sophisticated definition (24).

(16)

X Y Z

X Y

Z

Figure 5: First field of the reference solution and the numerical solution after forty-five time iterations for the first time dependent test case.

X Y

Z

X Y

Z

(17)

Time t 4 8 12

L2 error 0.581116E+00 0.719711E+00 0.635122E+00 L2 _{exterior error} _{0.531007E+00 0.715346E+00 0.634705E+00}

L2 _{interior error} _0.236067E+00 _0.791401E-01 _0.229959E-01

H(div) error 0.189652E+01 0.181267E+01 0.135066E+01 H(div) exterior error 0.183644E+01 0.180503E+01 0.135005E+01 H(div) interior error 0.473557E+00 0.166193E+00 0.403889E-01

Time t 4 8 12

L2 error 0.924122E-01 0.860510E-01 0.432997E-01 L2 _{exterior error} _0.886897E-01 _0.855219E-01 _0.432661E-01

L2 interior error 0.259643E-01 0.952856E-02 0.170669E-02 H(div) error 0.420503E+00 0.318012E+00 0.156138E+00 H(div) exterior error 0.418192E+00 0.316847E+00 0.156076E+00 H(div) interior error 0.440241E-01 0.272051E-01 0.439168E-02

Time t 4 8 12

L2 _error _0.139498E-01 _0.103088E-01 _0.449157E-02

L2 exterior error 0.137516E-01 0.102902E-01 0.449112E-02 L2 interior error 0.234338E-02 0.618670E-03 0.640346E-04 H(div) error 0.141093E+00 0.606411E-01 0.250898E-01 H(div) exterior error 0.141044E+00 0.606128E-01 0.250890E-01 H(div) interior error 0.369926E-02 0.185147E-02 0.199354E-03

(18)

Time t 4 8 12

L2 error 0.223674E+00 0.779898E+00 0.105697E+01 L2 _{exterior error} _{0.199908E+00 0.753085E+00 0.104881E+01}

L2 _{interior error} _{0.100334E+00 0.202743E+00 0.131048E+00}

H(div) error 0.481584E+00 0.149183E+01 0.218860E+01 H(div) exterior error 0.469498E+00 0.147570E+01 0.218318E+01 H(div) interior error 0.107218E+00 0.218813E+00 0.153950E+00

Time t 4 8 12

L2 error 0.267213E-01 0.686831E-01 0.874091E-01 L2 _{exterior error} _0.255068E-01 _0.680402E-01 _0.867683E-01

L2 interior error 0.796434E-02 0.937600E-02 0.105641E-01 H(div) error 0.575619E-01 0.195273E+00 0.228699E+00 H(div) exterior error 0.569467E-01 0.195035E+00 0.228447E+00 H(div) interior error 0.839323E-02 0.964150E-02 0.107185E-01

Time t 4 8 12

L2 _error _0.275997E-02 _0.615208E-02 _0.775471E-02

L2 exterior error 0.269564E-02 0.611029E-02 0.771919E-02 L2 interior error 0.592397E-03 0.715845E-03 0.741333E-03 H(div) error 0.160975E-01 0.203000E-01 0.235309E-01 H(div) exterior error 0.160859E-01 0.202872E-01 0.235192E-01 H(div) interior error 0.610490E-03 0.720756E-03 0.742916E-03

(19)

REFERENCES

[1] Y. Achou, Y. Maday The mortar element method with overlapping subdomains., SIAM J. Numer. Anal. 40(2), 601-628 (2002)

[2] R. Glowinski, J. He, J. Rappaz & J. Wagner A multi-domain method for solving numerically multi-scale elliptic problems , R. Glowinski and al., C. R. Acad. Sci. Paris, Ser. I 338 (2004).

[3] J.C. N´ed´elec A new family of mixed finite elements in IR3 , Numer. Math., 50:57-81, 1986.

[4] A. Quateroni, A. Valli Domain Decomposition Methods for Partial Differen-tial Equations, Numerical Mathematics and Scientific Computation. Oxford Science publications, Oxford, 1999.