Atmospheric sound propagation over large-scale irregular terrain

Atmospheric sound propagation over large-scale irregular terrain Martin Almquist

1 2

∗1

, Ilkka Karasalo

†2

and Ken Mattsson

‡1

Department of Information Technology, Uppsala University

Department of Hydroacoustics, Swedish Defence Research Agency, SE 17290 Stockholm, Sweden

September 16, 2013

Abstract

A benchmark problem on atmospheric sound propagation over irregular terrain has been solved using a stable fourth-order accurate nite dierence approximation of a high-delity acoustic model. A comparison with the parabolic equation method and ray tracing methods is made. The results show that ray tracing methods can potentially be unreliable in the presence of irregular terrain. Key words: wave propagation; irregular terrain; high-order nite dierence methods

1 Introduction High-delity simulations of sound propagation require an accurate treatment of the medium of propagation. In the atmosphere the speed of sound typically varies in space, which causes refraction of sound rays. As sound waves hit the surface of the earth they are partly reected; the degree and direction of reection depend on the topography and the type of ground. Measurements have shown that in order to accurately predict sound pressure levels (SPL), it is important to take the properties of both the atmosphere and the underlying terrain into account (see for example [40, 24]). In realistic applications, the atmospheric variations can be complex and the terrain is often irregular, which means both that the topography is non-trivial and that the type of ground can vary.

[email protected]; Corresponding author [email protected] ‡ [email protected]

∗ †

1

One high-delity model for sound propagation is based on the linearized Euler equations [41], which is most often solved with the staggered grid nite dierence time domain method originally presented in [48].

In [17] a

pseudo-spectral time domain method was used to solve the linearized Euler equations in urban courtyards.

This is an ecient approach when the

underlying geometry ts into a block-Cartesian geometry. Another high-delity model is the time-dependent acoustic wave equation. However, the wave equation has been considered too computationally demanding in realistic large-scale 3-D settings. One way to reduce the computational cost when the sound source consists of a single frequency is to Fourier transform the wave equation, which yields the Helmholtz equation. From the Helmholtz equation, it is possible to further reduce the computational cost via an assumption that reections at boundaries are negligible. The resulting model is called the parabolic equation (PE) model [25], often used in ocean acoustics.

Because of the assumption about reections,

the PE method is only valid for very moderate topographies. A signicant drawback with the Helmholtz equation (and the PE method) is that timedependent phenomena, such as a varying atmosphere or sources (for instance wind turbines) can not easily be approximated. However, despite the possible simplications, the above models are all generally regarded as too computationally expensive in realistic 3-D applications.

Hence, most numerical methods for sound propagation are based

on much simpler models.

A commonly used sound propagation model is

NORD2000 [1], based on ray tracing methods, which combine geometrical ray theory with the theory of diraction (for examples of applications see [46, 47, 23, 8]). Ray tracing methods are computationally cheap, but they incorporate neither a complex atmosphere nor irregular terrain properly. Unfortunately, due to a shortage of reliable measurements over largescale domains, it has been dicult to validate the sound propagation models currently used in practical applications. An alternative to measurements is validation through manufactured solutions (where one compares the numerical results with analytical solutions), but this technique has not yet been widely accepted outside the eld of applied mathematics. In this paper we shall solve a realistic sound propagation problem with the time-dependent acoustic wave equation as the underlying model.

As

mentioned above, this model is desirable from an accuracy perspective, but it is computationally costly.

It is therefore imperative to use an ecient

numerical method, to minimize the number of unknowns. It is well known that high-order (higher than second order) spatially accurate nite dierence schemes combined with high-order accurate time marching schemes are very well suited for wave propagation problems on large domains (see the pioneering paper by Kreiss and Oliger [22]). To guarantee an accurate approximation, it is necessary that the underlying numerical scheme can be proven stable, which is a non-trivial task using high-order nite dierence

2

methods. A robust and well-proven high-order nite dierence methodology that ensures stability of time-dependent partial dierential equations (PDEs) is the summation-by-partssimultaneous approximation term (SBP-SAT) method. The SBP-SAT method combines semi-discrete operators that satisfy a summation-by-parts (SBP) formula [21], with physical boundary conditions implemented using the simultaneous approximation term (SAT) method [4]. Examples of the SBP-SAT approach can be found in[36, 37, 38, 32, 34, 35, 39, 30, 44, 26, 7, 33, 14, 13, 16]. An added benet of the SBP-SAT method is that it naturally extends to multi-block geometries while retaining the essential single-block properties: stability, accuracy, and conservation [5]. Thus, problems involving complex domains or non-smooth geometries are easily amenable to the approach. References [30, 34, 15, 29] report applications of the SBP-SAT method to problems involving nontrivial geometries. Most of the published results for the SBP-SAT method are for rst order hyperbolic systems (and the NavierStokes equations). The extension of the SBP-SAT technique to the second order wave equation is found in [30, 31].

However, those studies were re-

stricted to 1-D concerning high-order nite dierence methods (3-D results for the nite volume technique were presented). The present study is a direct extension of the SBP-SAT method to the second order wave equation on multidimensional curvilinear domains, including non-trivial boundary conditions. In the rst part of this paper, the SBP-SAT method is applied to the benchmark problem introduced in [40]. The SBP-SAT method is here extended to the second order wave equation on a curvilinear 2-D domain with non-trivial boundary conditions. A fourth-order accurate SBP-SAT approximation is implemented and veried by a grid-convergence study against manufactured solutions.

In the second part of this paper, the benchmark

problem in [40] is used to compare the PE and ray tracing methods with the newly implemented SBP-SAT method. Conclusions as to the validity of the PE and ray tracing methods in the presence of a complex atmosphere and irregular terrain are drawn. The rest of the paper is organized as follows: In Section 2 we introduce some denitions and illustrate the SBP-SAT method by applying it to a 1-D problem.

In Section 3 we analyze a model problem in 2-D. We then

introduce the benchmark problem in Section 4, and describe how the SBPSAT technique has been adapted to this problem. The PE method and the ray interpolation methods are introduced in Sections 5 and 6, respectively. The implementation of the fourth-order accurate SBP-SAT method is veried in a series of convergence studies in Section 7.

In Section 8, we present

the results from the benchmark problem, comparing the dierent numerical methods. Conclusions are presented in Section 9. The nite dierence (SBP) operators are listed in Appendix.

3

2 The 1-D problem In this section we dene the SBP-SAT method.

To illustrate the power

and simplicity of the method we shall consider the following second-order hyperbolic problem:

0 ≤ x ≤ 1, x = 0, x = 1, 0 ≤ x ≤ 1,

autt = (bux )x , L0 u = g0 (t), L1 u = g1 (t), u = f1 , ut = f2 , where

a(x) > 0

and

b(x) > 0,

t ≥ 0, t ≥ 0, t ≥ 0, t = 0,

(1)

and the boundary conditions are given by

L0 u = α0 utt + α1 ut + α2 u + α3 ux = g0 (t), x = 0 L1 u = β0 utt + β1 ut + β2 u + β3 ux = g1 (t), x = 1. In the present study we restrict ourselves to the case

(2)

α3 6= 0, β3 6= 0, referred

to as mixed boundary conditions. For the treament of Dirichlet conditions we refer to [31]. Before we start employing the SBP-SAT method, some denitions are needed.

Let the inner product for real-valued functions

u, v ∈ L2 [0, 1]

R1

be

(u, v) = 0 u v c(x) dx, c(x) > 0, and let the corresponding norm kuk2c = (u, u). The domain (0 ≤ x ≤ 1) is discretized using the following N + 1 equidistant grid points:

dened by be

xi = i h, i = 0, 1, . . . , N, h = The approximate solution at grid point solution vector is

v = [v0 , v1 , . . . , vN ]T .

xi

1 N

is denoted

. vi ,

and the discrete

Similar to the continuous inner

product, we dene an inner product for discrete real-valued vector functions

u, v ∈ RN+1 by (u, v)Hc = uT HC v , where H is diagonal and positive denite and C is the projection of c(x) onto the diagonal. The corresponding 2 T norm is kvkH = v HC v . c

Remark

The matrix product

HC

denes a norm if and only if

metric and positive denite. This can only be guaranteed if

H

HC

is sym-

is a diagonal

matrix (see [43] for a detailed study on this). The following vectors will be frequently used:

e0 = [1, 0, . . . , 0]T , eN = [0, . . . , 0, 1]T .

(3)

2.1 The SBP-SAT method SBP operators are essentially central nite dierence stencils, closed at the boundaries with carefully chosen one-sided dierence stencils which mimic the underlying integration-by-parts formula in a discrete norm. In the present

4

paper we address the SBP operators by the accuracy of the central scheme and the type of norm which they are based on. A SBP operator is closed with

pth

2pth

order diagonal norm

order accurate one sided stencils (see [32]).

For rst order hyperbolic problems, this implies that the convergence rate (i.e., global convergence) drops to

(p + 1)th

order when using a

2pth

order

diagonal norm SBP operator. For strongly parabolic problems and second order hyperbolic problems the convergence rate instead drops to

(p + 2)th

order (see [10, 45] for more information on the accuracy of nite dierence approximations). To dene the SBP-SAT method, we present denitions 2.1-2.3 (rst stated in [35] and [30]). We here say that a scheme is explicit if no linear system of equations needs to be solved to compute the dierence approximation.

Denition 2.1 An explicit 2pth-order accurate nite dierence stencil with minimal stencil width of a Cauchy problem is called a 2pth-order accurate narrow-stencil. Denition 2.2 A dierence operator D1 = H −1 Q = H −1 (Q¯ − 12 e0 eT0 + 1 T 2 eN eN ) approximating ∂/∂ x, using a 2pth-order accurate narrow-stencil, is said to be a 2pth-order accurate narrow-diagonal rst-derivative SBP oper¯+Q ¯ T = 0. ator if H is diagonal and positive denite and Q Denition 2.3 Let D2(b) = H −1 (−M (b) − e0 b0 S0 + eN bN SN ) approximate ∂/∂ x ( b ∂/∂ x), where b(x) > 0, using a 2pth-order accurate narrow-stencil. (b) D2 is said to be a 2pth-order accurate narrow-diagonal second-derivative (b) is symmetric SBP operator, if H is diagonal and positive denite, M and positive semi-denite and S0 and SN approximate the rst-derivative operator at the boundaries. The superscript

(b)

emphasizes that

M (b)

and

(b)

D2

depend on

b(x).

The

explicit dependence can be found in [28]. For completeness we have included the fourth order SBP operators (used in the present study) in Appendix. The following denition introduced in [28] is also central in this paper:

Denition 2.4 Let D1 and D2(b) be 2pth-order accurate narrow-diagonal rst- and second-derivative SBP operators. If M (b) = D1T HBD1 + R(b) , and the remainder R(b) is positive semi-denite, D1 and D2(b) are called compatible. 2.1.1 Continuous analysis Multiplying the rst equation in (1) by

ut

and integrating by parts (referred

to as the energy method) leads to


d kut k2a + kux k2b = 2(bux ut )x=1 − 2(bux ut )x=0 , dt 5

(4)

where

(bux ut )x=1

means

(bux ut )

evaluated at

x = 1.

We also identify

E = kut k2a + kux k2b ,

(5)

as the total energy (kinetic and potential). Multiplying the rst equation in (1) by

ut ,

integrating by parts and

imposing the boundary conditions (2) leads to

d ¯ E = BTx=0 + BTx=1 , dt

(6)

where

¯ = kut k2a + kux k2 + b(0) (α0 u2t + α2 u2 )x=0 − b(1) (β0 u2t + β2 u2 )x=0 , E b α3 β3 and

 2 g02 b(0) g0 1 BTx=0 = + 2b(0)α u − t α3 2 α1 x=0 − 2 α1 α3  2 g 2 b(1) g1 1 u − + 21β1 β3 . BTx=1 = − 2b(1)β t β3 2 β1

(7)

(8)

x=1

Here we assume that

α3 6= 0

and

β3 6= 0.

The following Lemma is central in

the present study,

Lemma 2.5 Eq. (1) with boundary conditions (2) has a bounded energy in terms of initial and boundary data if α3 6= 0, β3 6= 0, α0 α3 ≥ 0, α2 α3 ≥ 0, β0 β3 ≤ 0, β2 β3 ≤ 0, α1 α3 < 0 and β1 β3 > 0 hold. Proof

¯ is non-negative and well E α2 α3 ≥ 0, β0 β3 ≤ 0, β2 β3 ≤ 0 hold.

dened if

α3 6= 0, β3 6= 0, α0 α3 ≥ 0,

By integrating (6) in time, we obtain

 2  2 ! 2b(1)β g (τ ) 2b(0)α g (τ ) 1 1 1 0 ¯ + − dτ = E(t) ut − ut − β3 2 β1 x=1 α3 2 α1 x=0 0  Z t 2 g1 (τ )b(1) g02 (τ )b(0) ¯ E(0) + − dτ . 2 β1 β3 2 α1 α3 0 Z

If

t

α1 α3 < 0

and

β1 β3 > 0

hold we have a strong estimate of

¯ E(t)

in terms

of initial and boundary data.

Remark

The implication of Lemma 2.5 is that we can have at most linear

time-growth of

kuk.

Linear growth (in t) does not violate well-posedness (see

[11]). However, linear time-growth of

kuk

β1 = 0,

f2 6= 0 α0 = α1 = β0 =

can only occur if we have

combined with pure Neumann boundary conditions, i.e,

so that we have a zero eigenvalue in the spectrum [27].

6

2.1.2 Semi-discrete analysis The discrete approximation of (1) using the SBP-SAT method is

(b)

Avtt = D2 v + τ0 H −1 e0 (L0 v − g0 ) + τ1 H −1 eN (L1 v − g1 ) ,

(9)

e0 and eN are dened in (3). (We assume the same initial conditions v = f1 , vt = f2 as in the continuous case). The matrix A has the values of a(x) injected on the diagonal. The semi-discrete boundary operators that

where

mimic (2) are given by

L0 v = α0 (vtt )0 + α1 (vt )0 + α2 v0 + α3 S0 v L1 v = β0 (vtt )N + β1 (vt )N + β2 vN + β3 SN v.

(10)

vtT H

and adding the

Applying the energy method by multiplying (9) by transpose leads to

d EH = − 2b0 (vt )0 S0 v + 2bN (vt )N SN v dt
+ 2τ0 α0 (vt )0 (vtt )0 + α1 (vt2 )0 + 2τ0 (α2 (vt )0 v0 + α3 (vt )0 S0 v − (vt )0 g0 )
+ 2τ1 β0 (vt )N (vtt )N + β1 (vt2 )N

(11)

+ 2τ1 (β2 (vt )N vN + β3 (vt )N SN v − (vt )N g1 ) , where

EH = kvt k2Ha + v T M (b) v .

Lemma 2.6 Eq. (9) with boundary operators (10) exactly mimics the continuous energy estimate (6) if τ0 = αb03 , τ1 = − bβN3 , and is thus stable if the conditions in Lemma 2.5 hold. Proof

insert

τ0 =

b0 α3 ,

τ1 = − bβN3

in (11) to obtain

d ¯ EH = BT0 + BTN , dt

(12)

where

¯H = kvt k2H + v T M (b) v + b0 (α0 (vt2 )0 + α2 (v 2 )0 ) − bN (β0 (vt2 )N + β2 (v 2 )N ) , E a α3 β3 and

 2 g2 b BT0 = + 2bα03α1 (vt )0 − 2gα01 − 2 α01 α0 3  2 g2 b BTN = − 2bβN3β1 (vt )N − 2gβ11 + 2 β1 1N β3 .

Equation (12) is a semi-discrete analogue to (6), and stability follows if the conditions in Lemma 2.5 hold.

7

Figure 1: The mapping between cartesian (left) and curvilinear (right) coordinates

3 Analysis in 2D In this section we analyze the scalar 2-D wave equation with mixed boundary conditions. To allow for complex domains, we transform the equation given on a curvilinear domain to an equation on a rectangular domain. We then derive an energy estimate for the continuous case.

After discretizing the

model in space with the SBP-SAT method, we prove stability by exactly mimicking the continuous energy estimate.

3.1 Denitions To make the notation more compact we introduce the Kronecker product:



···

c0,0 D . . .

 C ⊗D =

. . .

cp−1,0 D · · · C is a p × q matrix and D is (N + 1) × (N + 1) identity matrix.

where

an

 ,

Ω0 .

We will refer to

m×n

the rectangular domain.

(Nξ + 1)(Nη + 1)

matrix. We also let

Ω

IN

be the

we transform it to

as the physical domain and

Ω0

as

The rectangular domain is discretized using the

grid points:

 (ξi , ηj ) = The boundaries of

S

Ω

(13)

cp−1,q−1 D

If the problem is given on a curvilinear domain the unit square,



c0,q−1 D

i j , Nξ Nη Ω0

 ,

i = 0, 1, ..., Nξ ,

are denoted by

W

(west),

j = 0, 1, ..., Nη . N

(north),

E

(east) and

(south), respectively, as shown in Figure 1. The approximate solution at

(ξi , ηj ) is denoted by vij , and the discrete solution vector is v = [v00 , ..., v0Nη , v10 , ..., vNξ Nη ]T . The matrix RW is dened so that RW v is a vector with the same length as v and the same elements on the posi-

a grid point

tions corresponding to the west boundary, but zeros everywhere else. The

8

matrices

RN , RE

and

RS

are dened similarly for the north, east and south

boundaries, respectively. By

D1ξ

we denote the 2-D version of the narrow-stencil rst-derivative

  (b) ∂ ∂ ∂ D1 , approximating ∂ξ . Similarly, D2ξ approximates ∂ξ b ∂ξ . In the same manner, we let Hξ denote the 2-D version of the diagonal matrix (b) H , applied in the ξ -direction. D1η , D2η and Hη are dened similarly for the η -direction. To simplify the notation (without any restriction) we here assume Nξ = Nη = N . The 2-D operators can be neatly expressed in terms of the 1-D SBP operator

operators using the Kronecker product:

D1ξ

= D1 ⊗ IN , D1η = IN ⊗ D1

Hξ

= H ⊗ IN ,

Hη

= IN ⊗ H

EW

= e0 ⊗ IN ,

ES

= IN ⊗ e 0

RW

=

T , EW EW

RS

= ES EST

EE

= eN ⊗ IN , EN

RE e0

where the vectors

b

=

and

eN

T, EE EE

(14)

= IN ⊗ e N , T, = EN EN

RN

are dened in (3). Assuming that the coecient

is constant, we can also write

(b)

D2ξ

(b)

(b)

= D2 ⊗ IN , D2η

In the case of a variable coecient

b,

(b)

= IN ⊗ D2 .

(15)

however, (15) does not hold. To cover

also that case, we introduce the notation

(ξ)

bi (η) = b(ξi , η),

(η)

bi (ξ) = b(ξ, ηi ),

i = 0, 1, . . . , N.

(16)

We also dene



(b)

dij

       =      



D2

(η)

b0



!

       .      

i,j  ! (η) b1

D2

i,j ..

.

 ! (η) bN

D2

(17)

i,j We then have

(b) D2ξ



d11

(b)

d12

(b)

···

d1(N +1)

   =  

d21

(b)

d22

(b)

···

d2(N +1)

. . .

. . .

..

. . .

(b)

(b)

.

d(N +1)1 d(N +1)2 · · · 9

(b)

(b)



(b)

     

d(N +1)(N +1)

(18)

and

 (b)

D2η



 D2    =   

(ξ)

b0



   (ξ) b1

D2

..

.

  (ξ) bN

    .   

(19)

D2

3.2 The continuous problem We consider the following problem:

(x, y) ∈ Ω,

utt = b∆u

t≥0 (20)

γ1 u + γ2 ∇u · n + γ3 ut = 0 (x, y) ∈ ∂Ω, t ≥ 0 u = f1 , ut = f2 , (x, y) ∈ Ω, t = 0, where

b(x, y) > 0.

We have chosen homogeneous boundary conditions to

avoid unnecessary notation in the analysis. Similarly to the 1-D analysis in Section 2.1, the analysis holds for inhomogeneous conditions as well.

We

γ2 6= 0, which includes the important (γ1 = 0, γ2 = 1, γ3 = 0). We can add dissipation to (20) by adding a term b∇ · (σ∇ut ), σ(x, y) ≥ 0

also limit our present study to the case case of Neumann conditions

to the right hand side of the PDE. The added dissipation term will be used to create absorbing layers at articial boundaries in Section 4.1. Including the dissipation term, the problem reads

utt = b∆u + b∇ · (σ∇ut )

(x, y) ∈ Ω,

t≥0 (21)

γ1 u + γ2 ∇u · n + γ3 ut = 0 (x, y) ∈ ∂Ω, t ≥ 0 u = f1 , ut = f2 , (x, y) ∈ Ω, t = 0.

We now transform the problem to a rectangular domain. Assume that there is a smooth one-to-one mapping



from

Ω0

to

Ω.

The Jacobian

J

x = x(ξ, η) y = y(ξ, η),

of the transformation is

J = xξ yη − xη yξ . The

scale factors η1

η2 of the transformation are dened q q 2 2 η1 = xξ + yξ , η2 = x2η + yη2 .

and

as (22)

Since the mapping is one-to-one, the Jacobian is everywhere non-zero. By the chain rule, we have



uξ = ux xξ + uy yξ uη = ux xη + uy yη , 10

which is equivalent to

(

ux = uy =

u

Replacing

with

uxx =

1 J

1 J

uyy =

1 J

1 J

By adding

uxx

ux

1 J 1 J

(uξ yη − uη yξ ) = (uη xξ − uξ xη ) =

and

uy




(uξ xη − uη xξ ) xη uyy

(23)

in (23) yields

(uξ yη − uη yξ ) yη

and

1 J ((uyη )ξ − (uyξ )η ) 1 J ((uxξ )η − (uxη )ξ ) .

ξ


ξ

−

1 J

1 J

−

1 J

1 J

(uξ yη − uη yξ ) yξ




(uξ xη − uη xξ ) xξ

η


η

.

(24)

and rearranging terms, the rst equation in (21) can

be written as

˜ tt = ∆u ˜ +∆ ˜ σ ut , Ju

(ξ, η) ∈ Ω0 ,

(25)

where we have dened

˜ = (α1 uξ )ξ + (βuξ )η + (βuη )ξ + (α2 uη )η , ∆u ˜ σ u = (σα1 uξ )ξ + (σβuξ )η + (σβuη )ξ + (σα2 uη )η , ∆ α1 =



1 1 2 1 2 xη + yη2 , β = − (xη xξ + yη yξ ) , α2 = xξ + yξ2 , J J J

and

J J˜ = . b Using equation (23) to transform ∇u · n in the second equation in (21) yields the transformed boundary conditions:

 γ1 η2 u − γ2 (α1 uξ + βuη ) + γ3 η2 ut    γ1 η2 u + γ2 (α1 uξ + βuη ) + γ3 η2 ut  γ1 η1 u − γ2 (α2 uη + βuξ ) + γ3 η1 ut   γ1 η1 u + γ2 (α2 uη + βuξ ) + γ3 η1 ut

= 0, = 0, = 0, = 0,

(ξ, η) ∈ W (ξ, η) ∈ E (ξ, η) ∈ S (ξ, η) ∈ N.

(26)

The complete transformed problem is given by (25), (26) and the initial conditions stated in (21). Applying the energy method (here assuming that also the time derivative of the boundary condition (26) holds) leads to

d E =− dt

Z

γ3 + σγ1 η2 u2t dη − γ2

W

Z −

γ3 + σγ1 η1 u2t dξ − γ2

N

Z  −

Z E Z

γ3 + σγ1 η2 u2t dη γ2 γ3 + σγ1 η1 u2t dξ γ2

S

utξ utη

T    σα1 σβ utξ dΩ0 , σβ σα2 utη

Ω0

11

(27)

where

    Z  T  Z 1 ˜ 2 0 uξ α1 β uξ Jut dΩ + E= dΩ0 + BT  , uη β α2 uη 2

(28)

Ω0

Ω0 and

Z Z Z γ1 γ1 γ1 γ1 2 2 2 η2 u dη + η2 u dη + η1 u dξ + η1 u2 dξ BT = γ2 γ2 γ2 γ2 S Z W Z E ZN Z γ3 γ γ γ3 3 3 2 2 2 + σ η2 ut dη + σ η2 ut dη + σ η1 ut dξ + σ η1 u2t dξ. γ2 γ2 γ2 γ2 Z

W

E



N

(29)

S



α1 β 2 is positive denite since α1 > 0 and α1 α2 − β = β α2 (xξ yη − xη yξ )2 = J 2 > 0. Thus, the problem has an energy estimate if the The matrix

relations

γ1 ≥ 0, γ2

γ3 ≥0 γ2

(30)

hold. The last term in (27) implies damping of the energy for

σ > 0.

3.3 The semi-discrete problem In the semi-discrete setting we use the following notation for the matrices corresponding to the continuous variable coecients, for readability purposes: If

λ

denotes a variable coecient in the continuous setting, we here

denote the matrix with the values of the diagonal by

λ.

λ(ξ, η)

at the grid points injected on

There is no risk of confusion since it will always be clear

from context whether we are in a continuous or semi-discrete setting. The semi-discrete version of (26) is given by

 W L v = RW {γ1 η2 v − γ2 (α1 Sξ v + βD1η v) + γ3 η2 vt } = 0     LE v = R {γ η v + γ (α S v + βD v) + γ η v } = 0 1 2 2 1 ξ 1η 3 2 t E S  L v = RS {γ1 η1 v − γ2 (α2 Sη v + βD1ξ v) + γ3 η1 vt } = 0    N L v = RN {γ1 η1 v + γ2 (α2 Sη v + βD1ξ v) + γ3 η1 vt } = 0.

(31)

In the numerical scheme we also impose the time-derivative of the boundary conditions when

Lv = f ,

σ > 0.

we impose both

For instance, if we have the boundary condition

Lv = f

and

Lvt = ft

using the SAT technique.

The semi-discrete approximation of (25) and (26) using the SBP-SAT method is

˜ tt = D(α1 ) v + D1ξ βD1η v + D1η βD1ξ v + D(α2 ) v Jv 2η 2ξ (σα1 )

+ D2ξ

(σα2 )

vt + D1ξ σβD1η vt + D1η σβD1ξ vt + D2η

vt

+ τ1 Hξ−1 LW v + τ1 Hξ−1 LE v + τ2 Hη−1 LS v + τ2 Hη−1 LN v + τ3 σHξ−1 LW vt + τ3 σHξ−1 LE vt + τ4 σHη−1 LS vt + τ4 σHη−1 LN vt . 12

(32)

The rst main result of the present study is stated in the following theorem:

Theorem 3.1 The scheme (32) is stable if τ1 = τ2 = τ3 = τ4 − γ12 and (30) holds. Proof

vtT Hξ Hη

Applying the energy method by multiplying (32) by

and

adding the transpose leads to

d dt EH

= vtT (1 + τ1 γ2 )Hη α1 (−RW + RE )Sξ v

+

vtT (1 + τ2 γ2 )Hξ α2 (−RS + RN )Sη v

+

vtT (1

+ τ3 γ2 )Hη σα1 (−RW + RE )Sξ vt

+

vtT (1 + τ4 γ2 )Hξ σα2 (−RS + RN )Sη vt

+

vtT (1 + τ1 γ2 ) Hη (−RW + RE ) βD1η v

+

vtT (1 + τ2 γ2 ) Hξ (−RS + RN ) βD1ξ v

+

vtT (1 + τ3 γ2 ) Hη (−RW + RE ) σβD1η vt vtT (1 + τ4 γ2 ) Hξ (−RS + RN ) σβD1ξ vt

+

vtT (τ1 γ3 + τ3 σγ1 )Hη η2 (RW + RE ) vt

+

vtT (τ2 γ3 + τ4 σγ1 )Hξ η1 (RS + RN ) vt

+

(σα1 )

−vtT Hη Mξ

(σα2 )

vt − vtT Hξ Mη

vt

+

+

−2 (D1ξ vt )T σβHξ Hη (D1η vt ) , where

EH =

1 T ˜ 2 vt Hξ Hη Jvt 1 2 1 2 1 2

By choosing



+

 v + 2 (D1ξ v)T βHξ Hη (D1η v) +
−v T τ1 γ1 Hη η2 (RW + RE )v − v T τ2 γ1 Hξ η1 (RS + RN )v +
−vtT τ3 σγ3 Hη η2 (RW + RE )vt − vtT τ4 σγ3 Hξ η1 (RS + RN )vt . (α1 )

v T H η Mξ

(α2 )

v + v T Hξ Mη

τ1 = τ2 = τ3 = τ4 = − γ12

we obtain an energy estimate com-

pletely analogous to (27). If (30) holds, we have a non-growing energy.



4 Model problem 4.1 The continuous problem We consider sound propagation over the terrain shown in Figure 2. A source emitting spherical waves with a frequency of 50 Hz is placed at range height

z = 10

m.

13

r = 0,

100 90 80

Height z (m)

70 60 50 40 30 20 10 0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Range r (m)

Figure 2: The topography

The propagation of sound waves is governed by the acoustic wave equation

utt = b∆u, where

u

is the acoustic pressure and

b

(33)

is the square of the wave velocity. As

in Section 3.2, we introduce dissipation,

utt = b∆u + b∇ · (σ∇ut ).

(34)

(r, φ, z) and assuming φ-direction) results in the axi-

Expressing equation (34) in cylindrical coordinates symmetry in the azimuthal direction (the

symmetric two-dimensional restriction of (34),

1 1 1 utt = (rur )r + uzz + (σrur )r + (σuz )z . b r r

(35)

We construct the physical domain by introducing articial boundaries in the manner shown in Figure 3. In the simulations, the west boundary was

r = 1 m and the north boundary at heights ranging from z = 300 z = 750 m. The dot at the z -axis marks the location of the source, just

placed at m to

outside the west boundary of the domain. The boundary condition at the south boundary is a locally-reacting impedance condition given by (see [40])

 pω

0

c where

c

−

χ q u + ∇u · n + ut = 0, 2 c

(36)

n is the unit outward normal, ω0 is the angular χ is the curvature and p and q are real numbers that

is the wave speed,

frequency of the source, satisfy the relation

p + qi =

14

i , Zˆ

(37)

100 90

Artificial boundaries

80

Height z (m)

70 60

Computational domain

50 40 30 20 10 0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Range r (m)

Figure 3: Qualitative description of the physical domain.

The dot at the

z -axis

represents the point source

where

Zˆ is the normalized sound impedance, Zˆ = 18.3+17.5i.

χ

The curvature

is dened as

χ(r) = where

H(r)

−Hrr (r) (1 + Hr2 (r))3/2

,

denotes the height of the ground at horizontal position

(38)

r.

The north and east boundaries are articial boundaries that are introduced to truncate the unbounded domain. We need to avoid reections at these boundaries. Absorption of waves at articial boundaries is an important numerical issue. One possible approach is to apply an absorbing boundary condition (ABC), for example the rst order Engquist Majda ABC [9],

   √   ut + b∇u · n |z=zmax = 0,  √  √   ut + b∇u · n + b u |r=rmax = 0. 2r

(39)

However, a rst order Engquist Majda ABC is perfectly absorbing only at normal incidence. At

45◦

incidence on a straight boundary, 17% of the in-

coming wave is reected, and close to glancing the reection coecient tends to unity. Another, more ecient approach is to introduce an absorbing layer (AL) close to the articial boundary. creasing the dissipation coecient

σ

This can be done by smoothly infrom zero to some xed value.

The

dissipation must increase quickly enough that the waves are damped eciently inside the layer, but it must also increase smoothly enough that we avoid reections at the interface between the AL and the interior domain. The wider one can aord to make the layer, the easier it is to nd a function

σ

that fullls both criteria.

This AL is a special case of the method

presented by Appelö and Colonius in [3]. Appelö and Colonius also slowed

15

down waves inside the layer by stretching the grid at the boundaries and included higher-order dissipation operators for better performance.

Since

the focus in the present study is not on optimal absorbing layers, we here settle for a simple version of their AL. We will verify that the truncation of the domain does not aect the solution by placing the articial boundaries at dierent locations in the simulations. At the west boundary, the boundary condition is determined by the source.

Consider a point source with amplitude

a distance

r˜ from

(s)

u

f.

and frequency

   A r˜ (˜ r, t) = sin 2πf t − . r˜ c (r, z) = (0, z0 ).

Now let the source be located at source to a point

A

At

the source, the acoustic pressure is given by

(r, z)

(40)

The distance

r˜

from the

is given by

r˜ =

p r2 + (z − z0 )2 .

(41)

Combining (40) and (41) yields

p

A

u(s) (r, z, t) = p sin 2πf r2 + (z − z0 )2

t−

r2 + (z − z0 )2 c

We will impose this boundary data on the west boundary.

!! .

(42)

This can be

done using Dirichlet boundary condition by extending the SAT technique developed in [31].

However, in the present study we instead impose the

data using a mixed boundary condition (which allows for a stronger energy estimate in terms of boundary data),

√ ut +

√

√ √ b b (s) (s) (s) b∇u · n + u = ut + b∇u · n + u , 2r 2r

r = rmin .

(43)

(ξ, η) and perform a transformation (0 ≤ ξ ≤ 1, 0 ≤ η ≤ 1) onto the physical domain

We now introduce logical coordinates from the unit square

(r(ξ, η), z(ξ, η)) ∈ Ω.

Equation (35) transforms into

˜ tt = ∆u ˜ +∆ ˜ σ ut , Ju where we have dened

(ξ, η) ∈ Ω0 ,

rJ J˜ = , b

˜ = (α1 uξ )ξ + (βuξ )η + (βuη )ξ + (α2 uη )η , ∆u ˜ σ u = (σα1 uξ )ξ + (σβuξ )η + (σβuη )ξ + (σα2 uη )η , ∆ and

α1 =



r r 2 r 2 rη + zη2 , β = − (rξ rη + zξ zη ) , α2 = rξ + zξ2 . J J J 16

(44)

Table 1: Number of grid points corresponding to dierent resolutions. North boundary at

z = 500

m

Resolution (points per wavelength)

Nr

Nz

Nr · Nz

6 9 12

1801 2701 3601

451 676 901

0.81 · 106 1.83 · 106 3.24 · 106

The model that we solve is the equation (44) with the boundary conditions (36) and (43) at the south and west boundaries, and either the ABC (39) or the AL at the north and east boundaries. In the ABC approach, is identically zero. In the AL approach,

σ

σ

is non-zero close to the north and

east boundaries. The equation (44) has the same form as (25), and all the boundary conditions (36), (43) and (39) are of the mixed type analyzed in Section 3. Thus, the analysis performed in Section 3, proving well-posedness for the continuous problem and stability for the discrete scheme, holds for this model too.

4.2 Implementation details We have implemented a fourth-order SBP-SAT method of the model problem. The spatial discretization is thus fourth-order accurate in the interior scheme and second-order accurate in the boundary closures. For completeness we have included the operators (rst presented in [28]) in Appendix. The classical fourth-order accurate Runge-Kutta method was used for discretization in time. In order to apply the solver to the model problem, a computational grid must be constructed in the physical domain. Generating a good grid on a complex domain is not a trivial task. If the grid is not smooth enough, the convergence rate will decrease. In this case, the terrain prole has only two continuous derivatives, and hence we suspect that the grid will not support high-order accuracy. creating grids.

We have used Pointwise, a commercial software for

Figure 4 shows an example of a coarse grid generated in

Pointwise. Table 1 lists the number of gridpoints required when using 6, 9 and 12 points per acoustic wavelength, with the north boundary at a height of 500 m. In the simulations, we time-advanced until the solution became periodic in time and then computed the amplitude of the sound waves by measuring

|v|max ,

the maximum absolute value of the solution during one period,

one meter above ground.

The propagation loss

P

(measured in dB) was

computed as

 P = −20 · log

17

|v|max A

 ,

(45)

100 90 80

Height z (m)

70 60 50 40 30 20 10 0

0

200

400

600

800

1000 Range r (m)

1200

1400

1600

1800

2000

Figure 4: Coarse example grid

where

A

is the amplitude of the point source.

5 The parabolic equation method The PE method used in the model problem of Section 8 is described briey below. More details can be found in [20]. As an initial step, a smooth approximation

h(r)

of the ground height

as function of range is computed from the data, using an interpolating or a variance-reducing B-spline expansion [6, Ch. XI], the choice depending on the smoothness of the data. The geometry is then mapped from the physical

(r, z)

domain to a rectangle in the

(ξ, η)

plane by an orthogonal curvilinear

transform

r = r(ξ, η), such that

η

z = z(ξ, η),

rξ rη + zξ zη = 0,

(46)

is constant along the boundaries of the computational domain.

Assuming cylindrical symmetry, the Helmholtz equation for the complex pressure

u(ξ, η)

is

f (f −1 uξ )ξ + f (g −1 uη )η + k 2 a2 u = 0, where, using the unit dB/wavelength for the attenuation

f = ρa/rb, a = (rξ2 + zξ2 )1/2 ,

g = ρb/ra, b = (rη2 + zη2 )1/2 ,

k=

(47)

α, ω c (1

+ iα log(10) 40π ).

The PE approximations are derived by writing (47) in the form

T 2 u = (1 − L)u − k0−2 Ru,

18

(48)

where

R

k0

is a reference wavenumber,

T

and

L

are dierential operators, and

a function:

T u = −ik0−1 f 1/2 (f −1/2 u)ξ , Lu = R = Discarding the term

(49)

−k0−2 [f (g −1 (uη ))η + k 2 a2 u] (3fξ2 − 2f fξξ )/(4f 2 ).

Ru,

which is small since

+ u,

R ∼ ξ −2

(50) (51)

as

ξ

grows, equation

(48) is simplied to

T 2 u = (1 − L)u.

(52)

The PE schemes compute one-way solutions to (52) by solving

√

Tu = with the pseudo dierential operator of

1−Lu

√

(53)

1 − L replaced by a rational function

L, Tu =

analoguously with [2].

Pm

and

Qn

Pm (L) u, Qn (L)

(54)

are polynomials of degrees

m

and

n

in

the Padé approximation

√

1 − x = Pm (x)/Qn (x) + O(xm+n+1 ),

x → 0.

(55)

Thus the JEPE PE-approximations are

u = Qn (L)v, T Qn (L)v = Pm (L)v, n

Increasing the Padé order

(56)

m = max(n, 1),

n = 0, 1, 2, . . .

(57)

reduces the phase error as function of elevation

angle, but also increases the computational work. In practice, Padé orders

n = 0, 1, 2

are the most frequently used and correspond to the narrow angle

◦

◦

◦

(15 ), the wide angle (35 ) and the very wide angle (55 ) approximations, cf. [19, Sec. 6.2.4]. Equation (57) with initial conditions at ξ = ξ0 and boundary conditions η = 0 (the upper boundary) and η = −H (the ground) is solved using the method of lines. Thus, u, T , L and the boundary conditions are discretized at

vertically using a centered second-order nite dierence scheme, [42, Sec. 9].

The vertically discretized form of equation (57) is a system of ODEs

(omitting the indices

m

and

n)

d DQ(L)w = ik0 D[P (L) − Q(L)]w dξ

w(ξ0 ) = w0

(58)

for the scaled and wavenumber-shifted complex pressure

w(ξ) = e−ik0 ξ 1/2 v(ξ). 19

(59)

D

and

L

are diagonal and tri-diagonal matrix-valued functions of

tively, with

D

real and the imaginary part of

The initial prole

w0

L

ξ,

respec-

diagonal and non-positive.

is computed from the given source data (height and

vertical directivity), by low-pass ltering w.r.t. vertical wave number to the validity interval of the PE scheme. Equation 58 is then solved by a two-step fourth order A-stable second derivative method by Jeltsch (method J4 in [18]).

6 The ray interpolation method In the ray interpolation method the ground height is described by a smooth function of range identical to that in the PE method described in Section 5. A ray trajectory

(r(s), z(s))

where

s

is arc length, is a solution to the ODE

system [19, Sec. 3.2.1]

dr/ds = cos(φ) dz/ds = sin(φ)

(60)

dφ/ds = {sin(φ) ∂c/∂r − cos(φ) ∂c/∂z}/c. c = c(r, z) is the sound speed and φ = φ(s) the elevation angle of the ray.

In

the high-frequency limit, the ray trajectories are streamlines of the acoustic intensity eld i.e. propagation paths of acoustic energy. The waveeld at a

(r, z) is then a sum of contributions from all rays passing through (r, z) eigenrays from the source to (r, z)). Each eigenray contributes to the

point - all

sum with the eld inside an innitesimal tube surrounding the ray. With a

P0 at 1 m range, (r(s), z(s)) is

mono-frequency monopole source with amplitude of the ray-tube eld along an eigenray

n (s)

b P (f, s) = P0 α(s) ei2πf τ (s) eiπnc (s)/2 Πj=1 γj

the value

(61)

where

f

frequency cos φ0 α(s) = ray tube area factor rA(s) φ0 = φ(0) launch angle

(62)

A(s) = − sin(φ)∂r/∂φ0 + cos(φ)∂z/∂φ0

(65)

τ (s)

(63) (64)

travel time along ray

(66)

nc (s)

number of caustic points along ray

(67)

nb (s)

number of ground reflections along ray

(68)

γj The waveeld a source at

(0, zs )

0

reflection coefficient for j th ground reflection P (f, r)

at points

(r, h(r)) 0 < r < R

is computed in two steps.

20

(69)

on the ground from

First, a pre-dened number

K

of ray paths

(rj (s), zj (s))

are computed

by solving the ray ODEs with a fourth order Runge-Kutta method [12, p. 178] with variable stepsize and local error control.

The rays start from

(rj (0), zj (0)) = (0, zs ) with uniformly distributed vertical launch angles Φj , j = 1, ..., K , and the ODE system (60) is augmented with one equation each for the travel time τ (s) along the ray and the partial derivatives ∂r(s)/∂φ0 , ∂z(s)/∂φ0 , ∂φ(s)/∂φ0 with respect to launch angle φ0 . Rays are reected at the ground and terminated at the maximal range r = R or at the upper boundary z = Z of the computational region. The number and the source

the locations of caustic points and ground reections along each ray are determined. Then, for each point

(rj , h(rj ))

on a receiver grid, the eld

P (f, r)

is ob-

tained as a sum of contributions of the form (61) approximating the eigenrays by cubic interpolation to appropriate ray subsets.

7 Convergence results In this section we verify the implementation of the fourth-order SBP-SAT method and investigate the quality if the grid in a series of convergence studies. We will calculate the convergence rate

kvref − v (N2 ) kh kvref − v (N1 ) kh

q = log10 where

d

is the dimension (d

=2

is the discrete

l2

as

!

 / log10

vref with N

here),

the corresponding numerical solution

q

N1 N2

1/d ,

(70)

v (N ) is − v (N ) kh

is a reference solution, grid points and

kvref

norm of the error.

7.1 Convergence study without absorbing layer To separate the eects of the grid from the numerical method, we here present a convergence study on a smooth curvilinear grid.

We use the analytical

solution

(a)

u

A

sin 2πf (r, z, t) = p r2 + (z − z0 )2

!! p r2 + (z − z0 )2 t− c

(71)

f and A = 1, f = 12.5

which is the pressure eld created by a point source with frequency

A, located z0 = 10 m.

(r, z) = (0, z0 ).

amplitude

at

We here set

Hz and

In the computations we use the same setup as for

the benchmark problem, except that we here set the dissipation coecient

σ

to zero everywhere (in order to have an analytical solution) and impose

the analytical solution

u(a)

as initial and boundary data.

The setup with

the smooth domain and the initial condition is shown in Figure 5.

21

The

Figure 5: The analytical solution at time Table 2:

log(l2 − errors)

t=0

on a smooth grid

and convergence rates on a smooth grid, without

absorbing layer

Nr × Nz

log el2

q

81 × 21 161 × 41 321 × 81 641 × 161 1281 × 321

-0.26

0.00

-1.30

3.46

-2.62

4.45

-3.73

3.64

-4.88

3.80

convergence results are presented in Table 2. We note that we obtain the expected fourth order convergence rate. Next, we investigate the quality of the grid generated (using the commercial grid-generator Pointwise) for the benchmark problem by running a convergence study with exactly the same setup on that grid.

The results

are presented in Table 3. We note that we obtain approximately third order convergence on this grid and draw the conclusion that the grid, as expected, is not smooth enough to support high-order accuracy (higher than third order). When solving the benchmark problem, we will thus have to make do with third order convergence.

7.2 Convergence study with absorbing layer To verify also the implementation of the AL, we here perform a convergence

σ on the benchmark grid. The 3601×881 grid points was used as a reference solution.

study with a non-zero dissipation coecient solution obtained with

The results are presented in Table 4. Similar to Table 3, we obtain slightly less than third order convergence on this grid.

22

Table 3:

log(l2 − errors) and convergence rates on the grid generated for the

benchmark problem, without absorbing layer

Nr × Nz

log el2

226 × 56 451 × 111 901 × 221 1801 × 441

0.33

0.00

-0.53

2.87

-1.37

2.77

-2.19

2.72

Table 4:

q

log(l2 − errors) and convergence rates on the grid generated for the

benchmark problem, with absorbing layer

Nr × Nz

log el2

226 × 56 451 × 111 901 × 221 1801 × 441

0.46

0.00

-0.32

2.59

q

-1.09

2.56

-1.84

2.48

8 Computations We have solved the benchmark problem described in Section 4 for two different sound speed proles:

•

prole 1: Constant prole,

•

prole 2: Linear prole,

c = 340

m/s

c = c0 + kz

with

c0 = 340

m/s and

k = 0.1

s−1 . To guarantee a correct solution it is important to verify: 1) that the SPL at ground level is grid-converged, and 2) that reections at articial boundaries are negligible. We begin this section by investigating the eects of the articial boundary treatment, and then perform a grid-convergence study. The grid-converged results obtained with the SBP-SAT method are then compared with the results obtained with the PE and ray tracing methods.

8.1 Domain truncation To investigate the eect of the articial boundaries, we computed the propagation loss for dierent locations of the north boundary, for case 1. In Figure 6 we compare the rst order Engquist Majda ABC with the AL approach to truncate the domain at the north boundary. The eects of the reections using the rst order Engquist Majda ABC decrease as we move the north boundary higher, but even with the north boundary at a height of 2000 m the reections interfere with the interior waves and cause rapid oscillations

23

in the SPL at ground level. The spurious reections when using the AL are much smaller. After 500 m, the results do not change visibly. In the remaining computations we place the north boundary at

z = 500

m and employ

the AL approach to truncate the domain.

Remark

We also extended the domain and moved the east boundary further

to the right, but the location of the east boundary turned out to have no impact on the SPL at ground level.

8.2 Grid-convergence To verify that the discretization errors are negligible, we study how the computed propagation loss varies with grid renement. The results are shown in Figure 7. The curve obtained using 4.5 grid points per wavelength deviates signicantly from the others, while the curves corresponding to 6 and 12 grid points per wavelength are almost identical, i.e., indicating grid convergence. In the remaining simulations, grids with 12 points per wavelength were used.

8.3 Comparison of models Figure 8 shows the propagation loss, obtained with the SBP-SAT method, in the entire domain up to a height of 100 m. The eects of the refraction that occurs with prole 2 is most apparent far away from the source. In Figure 9 we compare the result obtained with the SBP-SAT method with the result published in [40] and the result obtained with the PE method described in Section 5, for prole 1. We observe that the results are in fairly good agreement. The maximum dierence between the SBP-SAT and the PE methods in Section 5 and in [40] are 8 dB and 4 dB, respectively. In Figure 10 we compare the ray tracing methods with the SBP-SAT method, for prole 1.

The dierence between the computed SPL using

the SBP-SAT method and the most accurate ray tracing method is always greater than 15 dB beyond 1100 meters. We also note that the results obtained with the dierent ray tracing methods dier signicantly from one another, and they all under-predict the SPL. In Figure 11 we compare the SBP-SAT method with the PE and ray tracing methods, with prole 2.

The PE method again shows reasonable

agreement with the SBP-SAT method. The ray tracing methods are here in better agreement with the SBP-SAT method than they were in Figure 10 (prole 1), but again they under-predict the SPL, except Ray interpolation that now over-predicts (except at the very end of the domain).

9 Conclusion The theory surrounding the SBP-SAT technique has been extended with a result that proves the stability of the SBP-SAT method for the second order

24

0 200 m 400 m 1000 m 2000 m

Propagation loss (dB)

20

40

60

80

100

120

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Range r (m)

(a) Engquist Majda ABC 0 300 m 500 m 750 m

10

Propagation loss (dB)

20 30 40 50 60 70 80 90

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Range r (m)

(b) AL Figure 6: Propagation loss measured 1 m above ground for dierent locations of the north boundary, using (a) rst order Engquist Majda ABC and (b) an AL

25

10 4.5 points/wavelength 6 points/wavelength 12 points/wavelength

Propagation loss (dB)

20 30 40 50 60 70 80 90

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Range r (m)

Figure 7: Convergence study.

The graphs show the propagation loss 1 m

above the ground for dierent levels of grid renement.

wave equation on a curvilinear 2-D domain with mixed boundary conditions. A fourth-order accurate SBP-SAT method has been applied to the benchmark problem on atmospheric sound propagation introduced in [40]. Since the SBP-SAT method is here applied to the full wave equation model, it can be used as a reference against which simpler (and computationally cheaper) methods can be validated. The present study has shown that, when applying the SBP-SAT method to sound propagation problems, the following should be considered:

•

The introduction of articial boundaries must not aect the solution. One way to achieve this is with carefully constructed absorbing layers.

•

A grid generated for a realistic topography might not support highorder accuracy.

•

The SPL must be grid-converged.

The results presented in Figures 10 and 11 show that ray tracing methods are not reliable for prediction of SPL in the case of irregular terrain. The PE methods show reasonable agreement with the SBP-SAT method, both with constant speed of sound and with a linear sound speed prole, which is expected since the topography in this problem is rather gentle. One would expect the PE methods, and the ray tracing methods in particular, to be more unreliable in the case of more pronounced topography. This is something we hope to address in a coming study.

26

(a) Prole 1

(b) Prole 2 Figure 8: Propagation loss (dB) for dierent sound speed proles, (a) prole 1 and (b) prole 2

27

0

Propagation loss (dB)

SBP−SAT Parakkal et al. I. Karasalo Topography 50

100

150

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Range r (m)

Figure 9: Propagation loss 1 m above ground for SBP-SAT and PE, with prole 1

Propagation loss (dB)

0 SBP−SAT Bell−shaped beams Gaussian beams Ray interpolation Topography 50

100

150

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Range r (m)

Figure 10: Propagation loss 1 m above ground for the SBP-SAT method and the ray tracing methods, with prole 1

28

Propagation loss (dB)

0 SBP−SAT PE Bell−shaped beams Gaussian beams Ray interpolation Topography

50

100

150

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Range r (m)

Figure 11: Propagation loss 1 m above ground for the SBP-SAT method, the PE and ray tracing methods, with prole 2

A Finite dierence operators For completeness we present the fourth order SBP operators. Here

h denotes

(b) ) are the standard the grid-spacing. The interior stencils (in D1 and M central

4th

order accurate nite dierence stencils.

At the boundaries we

use one-sided stencils that are formally second order accurate. The discrete norm

H

is dened:

 17 48

    H = h   

 59 48

43 48

       

49 48

1 ..

.

The rst derivative SBP operator is given by,

 24 − 17  −1  42 1  43 D1 =  3 h  98  0 

59 34

0 59 − 86 0 0

4 − 17

3 − 34 0

0 59 − 98

59 86

0 0 4 − 43

0 − 23

0

1 2

1 12

..

32 49

..

.

.

0 0 0 4 − 49 2 3

..

0 0 0 0 1 − 12 ..

.

The third-order accurate boundary derivative operator

S0 =

1  11 h 6

−3

3 2

29

− 13

0 0 ...



.

S0



..

        .

is given by,

The interior stencil of

−h M (b)

i

at row

is given by (i

mi, i−2 =

1 6

bi−1 −

1 8

bi−2 −

1 8

bi

mi, i−1 =

1 6

bi−2 +

1 6

bi+1 +

1 2

bi−1 +

1 2

1 mi, i = − 24 bi−2 −

5 6

bi−1 −

5 6

bi+1 −

1 24

bi+2 +

1 2

mi, i+1 =

1 6

mi, i+2 =

1 6

bi−1 +

1 6

bi+1 −

1 8

1 8

bi −

bi +

bi bi+2 −

3 4

bi

.

bi+1

bi+2

−h M (b)

The left boundary closure of

1 2

= 7 . . . N − 6):

6×6

(given by a

matrix) is given

by

m1, 1 =

12 17

m1, 2 =

− 59 68

m1, 3 =

2 17

m1, 4 =

3 68

m1, 5 =

49579087 10149031312

m1, 6 =

1 − 784

m2, 2 =

3481 3264

b1 +

59 192

27010400129 345067064608

b2 +

6025413881 21126554976

b1 −

b1 −

59 192

b1 −

1244724001 21126554976

b4 +

59 m2, 3 = − 408 b1 −

29294615794607 29725717938208

1328188692663 m2, 5 = − 37594290333616 b3 +

m2, 6 = m3, 3 =

1 51

m3, 4 =

1 136

m3, 5 =

− 10532412077335 42840005263888

m3, 6 =

960119 − 1280713392

m4, 4 =

3 1088

m4, 5 =

4959271814984644613 − 20965546238960637264

m4, 6 =

368395 − 2230716

m5, 5 =

8386761355510099813 128413970713633903242

m5, 6 =

35039615 − 213452232

m6, 6 =

3290636 80044587

b1 +

b3 +

8673 2904112

b1 −

b1 +

b3 −

b5 −

b4 −

b4 −

b4 +

b4 +

564461 13384296

1613976761032884305 7963657098519931984

3391 6692148

b5 +

b3 −

b3 +

b7 −

5580181 6692148

b4

b4

752806667 539854092016

1 6

b4

b4

4836340090442187227 5525802884687299744

b4 +

b4

60834186813841 1278205871342944

b3 −

507284006600757858213 475219048083107777984

b5 +

2944673881023 29725717938208

13777050223300597 26218083221499456

b2 +

125059 743572

236024329996203 1278205871342944

1328188692663 37594290333616

8673 − 2904112

59 192

b4

b4

b3 +

260297319232891 2556411742685888

59 b1 + m2, 4 = − 1088

b3

b3

9258282831623875 7669235228057664

b1 +

b4 +

b4

b4

2083938599 8024815456

752806667 21126554976

49579087 10149031312

b3 − 1 784

b3 −

b3 +

69462376031 2070402387648

537416663 7042184992

213318005 16049630912

b2 +

b3 +

1 6

b3 +

b6 −

b3 +

1063649 8712336

b7 +

1 24

b3 −

b8 +

b3

b4

b5

b3 1 24

b5 +

b4 +

2224717261773437 2763180339520776

5 6

564461 4461432

15998714909649 37594290333616

13091810925 13226425254392

b5 +

b3 +

1869103 2230716

378288882302546512209 270764341349677687456

17220493277981 89177153814624

b3 −

33235054191 26452850508784

b5 +

1 8

b6 +

b4 −

1950062198436997 3834617614028832

375177 743572

b4

b5

b6

b4 +

5 6

b6 +

1118749 2230716

1 24

b5 −

660204843 13226425254392

1 2

b7 +

280535 371786

b5

b6

b3 +

3 4

b6

The corresponding right boundary closure is obtained by replacing

bi →

bN +1−i for i = 1, . . . , 8 followed by a permutation of both rows and columns. (b) . The matrix M (b) is Let mi,j be the entry at row i and column j in M symmetric, which means that it is completely dened by the entries on and above the main diagonal, i.e.,

mj,i = mi,j , i = 1, . . . , N,

30

j = i, . . . , N .

References [1] Nord2000. comprehensive outdoor sound propagation model. part 2: Propagation in an atmosphere with refraction.

Technical Report

1851/00, Nordic Noise Group & Nordic Road Administration, 2006. [2] A. Bamberger, B. Engquist, L. Halpern, and P. Joly.

Higher order

paraxial wave equation approximations in heterogeneous media.

Journal of Applied Mathematics, 48:129154, 1988. [3] D. Appelö and T. Colonius.

SIAM

A high-order super-grid-scale absorbing

layer and its application to linear hyberbolic systems.

J. Comput. Phys.,

228:42004217, 2009. [4] M. H. Carpenter, D. Gottlieb, and S. Abarbanel. Time-stable boundary conditions for nite-dierence schemes solving hyperbolic systems: Methodology and application to high-order compact schemes.

put. Phys., 111(2), 1994.

J. Com-

[5] M. H. Carpenter, J. Nordström, and D. Gottlieb. A Stable and Conservative Interface Treatment of Arbitrary Spatial Accuracy.

Phys., 148, 1999. [6] Carl de Boor.

J. Comput.

A Practical Guide to Splines. Springer-Verlag, New York,

1978. [7] Peter Diener, Ernst Nils Dorband, Erik Schnetter, and Manuel Tiglio. Optimized high-order derivative and dissipation operators satisfying summation by parts, and applications in three-dimensional multi-block evolutions.

J. Sci. Comput., 32(1):109145, 2007.

[8] J. J. Embrechts. Broad spectrum diusion model for room acoustics raytracing algorithms.

The Journal of the Acoustical Society of America,

107(4):20682081, 2000. [9] B. Engquist and A. Majda.

Absorbing boundary conditions for the

numerical simulation of waves.

Mathematics of Computation,

31:629

651, 1977. [10] B. Gustafsson. The convergence rate for dierence approximations to general mixed initial boundary value problems.

SIAM J. Numer. Anal.,

18(2):179190, Apr. 1981. [11] B. Gustafsson, H.-O. Kreiss, and J. Oliger.

and dierence methods.

John Wiley & Sons, Inc., 1995.

[12] E. Hairer, S.P. Nørsett, and G.Wanner.

Equations I.

Time dependent problems

Springer Verlag, 1993.

31

Solving Ordinary Dierential

[13] J. Hicken and D. Zingg.

Superconvergent functional estimates from

summation-by-parts nite-dierence discretizations.

Scientic Computing, 33(2):893922, 2011.

SIAM Journal on

[14] J. E. Hicken and D. W. Zingg. Parallel newton-krylov solver for the euler equations discretized using simultaneous approximation terms.

Journal, 46(11):27732786, 2013/02/09 2008.

AIAA

[15] Jason E. Hicken and David W. Zingg. Aerodynamic optimization algorithm with integrated geometry parameterization and mesh movement.

AIAA Journal, 48(2):400413, 2013/02/09 2010. [16] J.E. Hicken.

Output error estimation for summation-by-parts nite-

dierence schemes.

Journal of Computational Physics,

231(9):3828 

3848, 2012. [17] Maarten Hornikx and Jens Forsse'n.

Modelling of sound propagation

to three-dimensional urban courtyards using the extended fourier pstd method.

Applied Acustics, 72:665676, 2011.

[18] R. Jeltsch. Multistep Methods Using Higher Derivatives and Damping at Innity.

Math. Comp., 31:124138, 1977.

[19] F.B. Jensen, W.A. Kuperman, M.B. Porter, and H. Schmidt.

tational Ocean Acoustics.

Compu-

AIP Press, New York, 1994.

[20] I. Karasalo and A. Sundström. JEPE - a high-order PE-model for rangedependent uid media. In

Proc. 3rd European Conference on Underwa-

ter Acoustics, pages 189194, Heraklion, Crete, Greece, 1996.

[21] H.-O. Kreiss and G. Scherer. Finite element and nite dierence meth-

Mathematical Aspects of Finite Elements in Partial Dierential Equations., Academic Press, Inc., 1974. ods for hyperbolic partial dierential equations.

[22] Heinz-Otto Kreiss and Joseph Oliger. Comparison of accurate methods for the integration of hyperbolic equations.

Tellus XXIV, 3, 1972.

[23] Samuli Laine, Samuel Siltanen, Tapio Lokki, and Lauri Savioja. celerated beam tracing algorithm.

Applied Acoustics,

Ac-

70(1):172  181,

2009. [24] Conny Larsson. Weather eects on outdoor sound propagation.

national Journal of Acoustics and Vibration, 5, 2000.

[25] E. Larsson and L. Abrahamsson.

Helmholtz and parabolic equation

solutions to a benchmark problem in ocean acoustics.

Am., 113:24462454, 2003.

32

Inter-

J. Acoust. Soc.

[26] L. Lehner, O. Reula, and M. Tiglio. Multi-block simulations in general relativity: high-order discretizations, numerical stability and applications.

Classical Quantum Gravity, 22:52835321, 2005.

[27] K. Mattsson. Boundary procedures for summation-by-parts operators.

Journal of Scientic Computing, 18:133153, 2003.

[28] K. Mattsson.

Summation by parts operators for nite dierence ap-

proximations of second-derivatives with variable coecients.

Scientic Computing, 51:650682, 2012.

Journal of

[29] K. Mattsson and M. Almquist. A solution to the stability issues with block norm summation by parts operators.

J. Comput. Phys., 253:418

442, 2013. [30] K. Mattsson, F. Ham, and G. Iaccarino. propagation in discontinuous media.

Stable and accurate wave

J. Comput. Phys., 227:87538767,

2008. [31] K. Mattsson, F. Ham, and G. Iaccarino. Stable boundary treatment for the wave equation on second-order form.

ing, 41:366383, 2009.

Journal of Scientic Comput-

[32] K. Mattsson and J. Nordström. Summation by parts operators for nite dierence approximations of second derivatives.

J. Comput. Phys.,

199(2):503540, 2004. [33] K. Mattsson and F. Parisi. Stable and accurate second-order formulation of the shifted wave equation.

Commun. Comput. Phys., 7:103137,

2010. [34] K. Mattsson, M. Svärd, M.H. Carpenter, and J. Nordström. High-order accurate computations for unsteady aerodynamics.

Computers & Fluids,

36:636649, 2006. [35] K.

Mattsson,

M.

Svärd,

and

M.

Shoeybi.

Stable

schemes for the compressible navier-stokes equations.

and

accurate

J. Comput. Phys.,

227(4):22932316, 2008. [36] J. Nordström and M. H. Carpenter. Boundary and interface conditions for high-order nite-dierence methods applied to the Euler and NavierStokes equations.

J. Comput. Phys., 148:341365, 1999.

[37] J. Nordström and M. H. Carpenter. High-order nite dierence methods, multidimensional linear problems, and curvilinear coordinates.

Comput. Phys., 173:149174, 2001. [38] J. Nordström and J. Gong. problems.

J.

A stable hybrid method for hyperbolic

J. Comput. Phys., 212:436453, 2006. 33

[39] J. Nordström, K. Mattsson, and R.C. Swanson. Boundary conditions for a divergence free velocity-pressure formulation of the incompressible navier-stokes equations.

J. Comput. Phys., 225:874890, 2007.

[40] S. Parakkal, KE. Gilbert, D. Xiao, and HE. Bass.

A generalized po-

lar coordinate method for sound propagation over large-scale irregular terrain.

J Acoust Soc Am., 128(5):25732580, 2010.

[41] T. Van Renterghem and D. Botteldooren. Prediction-step staggered-intime fdtd: An ecient numerical scheme to solve the linearised equations of uid dynamics in outdoor sound propagation.

Applied Acustics,

68:201216, 2007. [42] A. Sundström. Energy-conserving parabolic wave equations. FOA Report C20895-2.7, National Defence Research Establishment, Stockholm, 1992. [43] M. Svärd. On coordinate transformation for summation-by-parts operators.

Journal of Scientic Computing, 20(1), 2004.

[44] M. Svärd, M. H. Carpenter, and J. Nordström.

A stable high-order

nite dierence scheme for the compressible NavierStokes equations, no-slip wall boundary conditions.

J. Comput. Physics,

227:48054824,

May 2008. [45] M. Svärd and J. Nordström.

On the order of accuracy for dierence

approximations of initial-boundary value problems.

J. Comput. Physics,

218:333352, October 2006. [46] Henry Weinberg and Robert Burridge. Horizontal ray theory for ocean acoustics.

The Journal of the Acoustical Society of America, 55(1):63

79, 1974. [47] Chang-Fa Yang, Boau-Cheng Wu, and Chuen-Jyi Ko.

A ray-tracing

Antennas and Propagation, IEEE Transactions on, 46(6):907919, 1998. method for modeling indoor wave propagation and penetration.

[48] K.S. Yee.

Numerical solution of initial boundary value problems in-

volving Maxwell's equations in isotropic media.

Propag., 14:302307, 1966.

34

IEEE Trans. Antennas