UPTEC F 16046
Examensarbete 30 hp Augusti 2016
Modeling Sound Propagation from Wind Turbines using Linearized 3D Euler Equations
Ylva Rydin
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
Abstract
Modeling Sound Propagation from Wind Turbines using Linearized 3D Euler Equations
Ylva Rydin
In this report sound propagation from wind turbines is modeled using the 3D linearized Euler equations which allows varying atmospheric fields containing wind.
The Euler Equations are linearized and a stable high order finite difference
approximation of the equations in 1D and 3D is obtained. Stability and convergence of the problem are proven and numerically verified.
Ämnesgranskare: Ken Mattsson Handledare: Jonatan Werpers
Sammanfattning
Globalt byggs det fler och större vindkraftverk idag än någonsin tidigare. När vindkraftsparker byggs är det viktigt att de placeras så att bullernivåerna inte blir ett störningsmoment i närliggande bostadsområden. Att förutspå bullerut- bredning från vindkraftparker noggrant är avancerat eftersom ljudutbredningen beror av många olika faktorer. Till exempel vind, atmosfäriska variationer av tryck och densitet samt atmosfärisk dämpning. I några av de tidigare pro- jekt där finita differenser har använts för att beräkna ljudutberdning från vind- kraftverk har effekten från vind inte tagits med i beräkningarna [4],[5]. I detta arbete har denna effekt inkluderats genom att beräkna ljudutbredning från vin- dkraftverk med Eulers ekvationer. Dessa ekvationer är ett system av första ordningens ickelinjära partiella differentialekvationer som beskriver rörelsen av ickeviskösa fluider. Begreppet fluid innefattar även gaser inom fluidmekanik och därför kan luftflöden beräknas med dessa ekvationer.
Ekvationerna har först lineariserats och sedan diskretiserats i rummet med finita differenser av hög ordning då dessa är väl anpassade för att lösa vågut- bredingsproblem. Vi har visat att problemet är välställt och metoden stabil vilket betyder att data i simuleringen garanterat inte kommer växa på ett icke fysikaliskt sätt samt att den numeriska lösningen konvergerar.
Fokus har legat på ljudnivåer på avstånd i storleksordningen kilometer ifrån vin- dkraftverken. Då dessa avstånd är stora i förhållande till vindkraftverkens stor- lek har vindkraftverken approximerats med punktkällor som avger ljud. Ljudet har i sin tur modellerats som tryckvariationer i punktkällan. Målet med detta ar- bete har varit att linearisera Eulers ekvationer samt ta fram och impelementera en stabil diskretisering i 3D. Vi har visat att den numeriska implementationen är stabil och konvergerar som förväntat. Metoden bör i ett framtida projekt utvidgas för att kunna hantera icke-reguljära geometrier så att mer realistiska vindkraftpaker kan modelleras. Om koden dessutom paralleliseras kommer stora och nogranna beräkningar på verkliga vindkraftspaker kunna genomföras.
Contents
1 Introduction 1
2 Euler Equations 2
2.1 Linearization . . . 3
2.2 Symmetrization . . . 5
3 Modeling Atmospheric Properties 6 3.1 Background Field . . . 6
3.2 Atmospheric Absorption . . . 8
3.3 Point Source . . . 8
3.4 Boundary Conditions . . . 9
4 The 1D Problem 10 4.1 Continuous Analysis . . . 10
4.2 Boundary Analysis . . . 12
4.3 Semi-discrete Stability Analysis . . . 14
4.3.1 Definitions . . . 14
4.3.2 Discretization . . . 15
4.3.3 Stability Analysis . . . 15
5 The 3D Problem 17 5.1 Continuous Analysis . . . 17
5.2 Semi-discrete Stability Analysis . . . 18
6 Computations 1D 21 6.1 Point Source . . . 21
6.2 Convergence . . . 22
6.2.1 Analytic Solution . . . 22
6.2.2 Reference Solution . . . 23
6.3 Stiffness of Semi-discrete Problem . . . 24
6.4 Comparison between Background Fields and Absorption . . . 25
7 Computations 3D 27 7.1 Convergence . . . 27
7.1.1 Analytic Solution . . . 27
7.2 Manufactured Solution . . . 28
7.3 Comparison of Background Fields and Boundary Conditions . . . 29
8 Conclusions and Further Improvements 32
A Analytic Solution to 1-Dimensional Linearized Euler with Forc-
ing 35
1 Introduction
Wind power is an important part of the renewable energy development. There- fore, the global market has been growing rapidly in recent years and this devel- opment is predicted to continue [1]. An issue when building wind farms is that there are restrictions on noise levels in residential areas. Hence, it is important to accurately predict the sound propagation from the wind farms on beforehand so that the turbines can be placed at a safe distance from these areas [2]. To model sound propagation from wind turbines is a difficult task since it depends upon factors of the surroundings. For example, it depends on the topology of the area, atmospheric properties such as wind, variations in air density, air hu- midity and atmospheric damping.
Sound propagation from wind turbines has in earlier projects been modeled by solving the acoustic wave equation with finite difference methods [4], [5]. A problem with the acoustic wave equation is that it does not take effects of wind into account. To obtain an extension of these earlier models including wind the Euler equations on linearized form are used to model sound propagation in this project. To get a stable and accurate approximation a high order finite difference discretization with Summation-By-Parts (SBP) properties is used as spatial discretization method in this project. This method also provides useful tools for stable boundary treatment. Since the focus in this work is far field simulations the wind turbines are modeled as point sources which is accurate far away from the source. The goals of this thesis has been:
• Linearizing the Euler equations in a way that allow varying background fields.
• Writing down a well-posed mathematical model of sound propagation from wind turbines based on the linearized Euler equations with physically cor- rect boundary conditions.
• Discretize the model in a stable way using finite differences.
• Numerically verify convergence and stability of the finite difference scheme.
The Euler equations are first introduced and linearized in Section 2. In Section 3 the atmospheric properties, such as wind and varying speed of sound, are discussed and a method to find atmospheric background fields consistent with the Euler equations is explained. Further, a point source that models wind tur- bines is introduced and the required boundary conditions to model the ground and the domain boundaries towards the atmosphere are discussed. In Section 4, well-posedness is shown for the 1D linearized Euler equations, a finite dif- ference scheme for a numerical solution is introduced and stability is shown for the scheme. In Section 5, a well-posed 3D scheme is introduced. In Section 6 and 7, numerical experiments verifying the stability and convergence of the numerical method are presented. Further, experiments investigating the impact of different background properties are performed in 1D and 3D.
2 Euler Equations
In this section the Euler Equations are introduced on conservative form. To simplify the analisys and improve performance of the numerical implementation the Euler equations are linearized and symmetrized. To do this, the equations are rewritten on primitive form and the atmosphere is assumed to behave as an ideal gas, using the ideal gas law as the equation of state. Finally, the equations are linearized considering perturbations around a solution of the full non-linear equations.
The Euler Equations are the governing equations of non-viscous fluid flow. The equations can be obtained from the well known Navier-Stokes equations by ne- glecting viscous effects. They can be written on conservative form as:
∂U
∂t +∂F
∂x +∂G
∂y +∂H
∂z = ˜f, (1)
where U is given by (ρ, ρu, ρv, ρw, w)T, ˜f is a forcing function dependent on time and space but not the solution,
F =
ρu ρu2+ p ρuv + p ρuw + p (e + p)u + pv + pw
, G =
ρv ρvu + p ρv2+ p ρvw + p (e + p)v + pu + pw
,
and
H =
ρw ρuw + p ρvw + p ρw2+ p (e + p)w + pu + pv
.
Here ρ is the density, u, v, w, velocity components in the x, y and z directions, pthe pressure and e the internal energy. These equations represents the conser- vation of mass, momentum and internal energy respectively. The ideal gas law p = (γ − 1)(e −12ρ(u2+ v2+ w2))relates the pressure velocity, internal energy and density. Here γ = 1.4 is the heat capacity ratio for air. In [3] it is shown that by using the ideal gas law as the equation of state the Euler equations (1) can be rewritten on non-conservative form as
∂v
∂t + ˆA(v)∂v
∂x+ ˆB(v)∂v
∂y + ˆC(v)∂v
∂z = ˆf, (2)
where the primitive variables v = (ρ, u, v, w, p)T, and the coefficient matrices are given by,
A(v) =ˆ
u ρ 0 0 0
0 u 0 0 1ρ
0 0 u 0 0
0 0 0 u 0
0 γp 0 0 u
, B(v) =ˆ
v 0 ρ 0 0
0 v 0 0 0
0 0 v 0 1ρ
0 0 0 v 0
0 0 γp 0 v
,
and
C(v) =ˆ
w 0 0 ρ 0
0 w 0 0 0
0 0 w 0 0
0 0 0 w 1ρ
0 0 0 γp w
.
Some of the analysis and numerical experiments in this work are performed in 1D since the 1D equations are easier to work with and the 1D results gen- eralize well to 3D. The 1D equations can be obtained by assuming no flow or variation of the solution in the y and z directions. This corresponds to setting the wind speeds v, w = 0 and assuming ∂v∂y = ∂v∂z = 0. The 1D case the Euler equations are reduced to following system of equations,
∂v
∂t + ˆA(v)∂v
∂x = ˆf, (3)
where v = (ρ, u, p)T.The 1D coefficient matrix is
A =ˆ
u ρ 0
0 u 1ρ 0 γp u
.
2.1 Linearization
To linearize the Euler Equations (2) we consider a small perturbation v1around a known solution to the full non-linear equations, v0. This is done by splitting the field into two parts as v = v0+ εv1, where v0 will be referred to as the background flow and εv1 is a small perturbation around the background flow on which the computations will be performed. The forcing function is also split as ˆf = ˆf0+ εˆf1.The linearization is done the same way for all three coefficient matrices. Therefore, the linearization will be demonstrated for the 1D problem,
∂v
∂t + ˆA(v)∂v
∂t = ˆf . (4)
In this section the Einstein summation notation [10] will be used which allows rewriting (4) as
∂vi
∂t + ˆAij(v)∂vj
∂t = ˆfi. (5)
where,
A =ˆ
Aˆ11 Aˆ12 · · · Aˆ1m
Aˆ21 Aˆ22 · · · Aˆ2m ... ... ... ...
Aˆm1 Aˆm2 · · · Aˆmm
, v =
v1
v2
...
vm
and i, j = 1, 2, . . . , m.
This notation is used to suppress the summation sign by letting indexes appear- ing twice in a term indicate summation over that index. For example matrix- vector multiplication can be written as
Av =ˆ X
j
Aˆijvj= ˆAijvj,
using this notation. To obtain the linearization, the coefficient matrix, ˆAij(v) is Taylor expanded around v0in ε as:
Aˆij(v0+ εv1) = ˆAij(v0) + ε∂ ˆAij(v0)
∂vk
v1k+ O(ε2). (6) Substituting v = v0+ εv1 and the Taylor expanded coefficient matrix in (6) into equation (5) gives,
∂vi0
∂t + ε∂vi1
∂t + ˆAij(v0) + ε∂ ˆAij(v0)
∂vk
v1kh∂v0j
∂x + ε∂vj1
∂x i
+ O(ε2) = ˆfi. which can be rearranged as:
ε ∂vi1
∂t + ˆAij(v0)∂vj1
∂x +∂ ˆAij(v0)
∂vk
∂vj0
∂xvk1+
∂v0i
∂t + ˆAij(v0)∂vj0
∂x − ˆfi0
| {z }
=0
+ε2∂ ˆAij(v0)
∂vk
∂vj1
∂xvk1= ε ˆfi1. (7)
Since the background field v0is assumed to solve the non linear Euler equations (4), equation (7) can be written as
∂v1i
∂t + ˆAij(v0)∂vj1
∂x +∂ ˆAij(v0)
∂vk
∂v0j
∂xv1k+ ε∂ ˆAij(v0)
∂vk
∂v1j
∂xvk1= ˆfi1. (8) By neglecting the comparably small ε terms, the following linear approximation of the Euler Equations is obtained
∂v1
∂t + ˆA(v0)∂v1
∂x + ˆE(v0)v1= ˆf1, (9) where
E(vˆ 0) =∂ ˆAij(v0)
∂vk
∂v0j
∂x =
u0x p0x 0
−(ρp00x)2 u0x 0 0 p0x γu0x
.
The 3D linearized Euler equations are easily obtained by linearizing ˆB and ˆC the same way as ˆA, the resulting equations are
∂v1
∂t + ˆA(v0)∂v1
∂x + ˆB(v0)∂v1
∂z + ˆC(v0)∂v1
∂y + ˆE(v0)v1= ˆf1 (10) where
E(vˆ 0) = ∂ ˆAij(v0)
∂vk
∂vj0
∂x +∂ ˆBij(v0)
∂vk
∂vj0
∂y +∂ ˆCij(v0)
∂vk
∂vj0
∂z
=
u0x+ v0y+ w0z p0x p0y p0z 0
−(ρp00x)2 u0x u0y u0z 0
− p
0 y
(ρ0)2 v0x vy0 v0z 0
−(ρp00z)2 w0x w0y w0z 0 0 p0x p0y p0z γ(u0x+ v0y+ w0z)
.
2.2 Symmetrization
The coefficient matrices ˆA, ˆB and ˆC can be symmetrized simultaneously by a change of variables derived in [3]. The change is described by v1= S−1q, where
S =
√γρ0
c0 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
ρ0c0
√γ 0 0 0 qγ−1
γ ρ0c0
.
The new symmetric equations are written as
qt+ Aqx+ Bqy+ Cqz+ Eq = f (11) where q = S−1v = (ρ0c√0
γρ1, u1, v1, w1, − c0
ρ0√
γ(γ−1)ρ1+
√γ c0ρ0√
γ−1p1)T and f = S−1ˆf. The new coefficient matrices are
A(v0) = S−1A(vˆ 0)S =
u0 √c0γ 0 0 0
c0
√γ u0 0 0 qγ−1
γ c0
0 0 u0 0 0
0 0 0 u0 0
0 q
γ−1
γ c0 0 0 u0
,
B(v0) = S−1B(vˆ 0)S =
v0 0 √c0γ 0 0
0 v0 0 0 0
c0
√γ 0 v0 0 qγ−1
γ c0
0 0 0 v0 0
0 0 qγ−1
γ c0 0 v0
,
C(v0) = S−1C(vˆ 0)S =
w0 0 0 √c0γ 0
0 w0 0 0 0
0 0 w0 0 0
c0
√γ 0 0 w0 qγ−1
γ c0
0 0 0 qγ−1
γ c0 w0
and,
E(v0) =
L √ρ0xγρc00
ρ0yc0
√γρ0
ρ0zc0
√γρ0 0
−
√γp0x
ρ0c0 u0x u0y u0z 0
−
√γp0y
ρ0c0 v0x v0y vz0 0
−
√γp0z
ρ0c0 wx0 wy0 w0z 0
L√
γ − 1 −p0(ρρ0x0−p)2c0x0ρ0a −p
0ρ0y−p0yρ0
(ρ0)2c0 a −p0(ρρ0z0−p)2cz00ρ0a γL
,
where L = u0x+ vy0+ w0z, the constant a =q γ
γ−1 and c0= qγp0
ρ0 is the speed of sound.
The 1D linearized Euler equations,
qt+ Aqx+ Eq = f (12)
are obtained in the same way by multiplication from the left by the inverse of the matrix
S =
√γρ0
c0 0 0
0 1 0
ρ0c0
√γ 0 qγ−1
γ ρ0c0
. Here q = (ρ0c√0
γρ1, u1, − c0
ρ0√
γ(γ−1)ρ1+
√γ c0ρ0√
γ−1p1)T, f = S−1ˆf,
A =
u0 √c0γ 0
c0
√γ u0 qγ−1
γ c0 0 qγ−1
γ c0 u0
,
and
E =
u0x √ρ0xγρc00 0
−
√γp0x
ρ0c0 u0x 0
u0x√
γ − 1 −p0(ρρ0x0−p)2cx00ρ0
q γ
γ−1 γu0x
.
3 Modeling Atmospheric Properties
When modeling sound propagation from wind turbines atmospheric properties such as speed of sound, wind and absorption has great impact. For example, you might have experienced that it is easier to hear someone who screams to you in the wind direction than someone screaming towards it. In this section a set of background fields that solves the non-linear Euler equations are de- rived. Further, an absorbing damping term is constructed and added to the Euler equations. Moreover, the point source is introduced and the boundary conditions used within the project are discussed.
3.1 Background Field
The main reasons for modeling sound propagation with linearized Euler equa- tions instead of the simpler wave equation is that it allows more complex back- ground fields with for example, wind and varying pressure and density. In this section a set of non-constant background fields are derived to allow simulating different types of atmospheric conditions. This is done first in 1D and later in 3D.
In the linearizion of the Euler equations (10), it is assumed that the background field is a solution to the full Euler Equations. Hence the 1D background field v0 must be a solution to
dv0
dt + ˆA(v0)∂v0
∂x = ˆf0, (13)
where the forcing function ˆf0 is restricted to be zero in this study. To sim- plify the computations, only time independent background fields are considered which means that dvdt0 = 0. Hence, solutions to following system of differential equations are studied
(ρu)x= 0 (14a)
uux+px
ρ = 0 (14b)
γpux+ upx= 0. (14c)
If all parameters are non-zero the solution to equation (14a) is ρu = k1 where k1 is a constant. This allows rewriting (14b) as (k1u + p)x = 0, which gives p = k2− k1u. Inserting this into (14c) gives u = k3. Hence the background field must be constant if all parameters are non-zero.
If instead no wind (u = 0) is assumed the remaining part of equation (14b) is pρx = 0. Now a possible background field is
v0=
ρ(x)
0 k
,
where k > 0 is a constant and ρ(x) > 0 is any strictly positive function.
In the 3D case the background field has to solve the 3D Euler equations (2). If only time independent solutions are considered it has to solve
(ρu)x+ (ρv)y+ (ρw)z= 0, uux+ vuy+ wuz+pρx = 0, uvx+ vvy+ wvz+pρy = 0, uwx+ vwy+ wwz+pρz = 0, γp(ux+ vy+ wz) + upx+ vpy+ wpz= 0.
In the model it is possible to impose a variety of background fields. In the present study background fields are assumed vary only in the z-direction, which corresponds to background fields varying with altitude. Hence the background fields are solutions to
(ρw)z= 0, wuz= 0, wvz+pρy = 0, wwz+pρz = 0, γpwz+ wpz= 0.
As in the 1D case this system requires w = 0 for other solutions than a constant background field. A possible solution is
v0=
ρ(z) u(z) v(z) 0 k
,
where k > 0 is constant and ρ(z) > 0, u(z) and v(z) are real valued functions.
3.2 Atmospheric Absorption
The linearized Euler equations does not take atmospheric absorption into ac- count. In earlier projects, [5] the sound propagation has been modeled with the acoustic wave equation,
ptt= c2pxx− βpt. (15)
In this case, damping has been added to the model by introducing the term
−βpt. The sound attenuation coefficient β governs the magnitude of the damp- ing, it depends on the frequency of the sound, the air humidity and temperature [8]. To add a similar damping to the linearized and symetrized Euler equations the matrix E can be extended with the damping term
ED= E +
α1 0 0 0 α2 0 0 0 α3
, (16)
where α1, α2 and α3, are damping parameters. This term damps variations in both velocity u and pressure p. This parameter has to be adjusted after the modeled conditions.
3.3 Point Source
In the Euler equations sound generators such as wind turbines are commonly modeled with a point source term. In 1D, following time dependent source term that triggers the pressure is us used:
f = p0g(t)δ(x − x0)
1
√ 0 γ − 1
.
This source is located at x0, p0 is the source amplitude, δ(x) the Dirac distri- bution and g(t) a source function which describes the pressure variations at the source. The pressure component of the solution to (12) with a Gaussian source function g(t) = √12e−(t−1.8)22σ2 with σ = 3.51 is displayed in Figure (1).
0 500 1000 1500 2000 2500 3000 3500 4000 4500 x
1 2 3 4
t
Figure 1: The time evolution of a Gaussian pulse generated from a point source at x0= 2500
3.4 Boundary Conditions
Two types of boundary conditions are needed in the model. Firstly a non- reflecting boundary condition that lets the sound leave the domain is required for the boundaries towards the atmosphere. Secondly a boundary condition to model sound reflection at the ground is needed. In this section both types of boundary conditions are introduced in 1D. These boundary conditions will be used in the later analysis and numerical experiments in both 1D and 3D. The non-reflecting boundary conditions used in this work are the so called charac- teristic boundary conditions,
A+q = 0, x = a
A−q = 0, x = b. (17)
Here the matrices A+ and A−, derived in Section 4.2 are constructed so that the ingoing characteristic variables are set to zero at the boundaries. These boundary conditions are non-reflecting since they do not let any information enter the domain without preventing information from leaving. To model sound reflection on the right boundary (x = b) the following boundary conditions,
A+q = 0, x = a
0 1 0 q = 0, x = b. (18)
are imposed. The boundary condition at x = b is a wall-boundary condition that reflects the sound. Here the velocity u1 = 0 is imposed on the boundary which gives reflection since the sound can not leave the domain. To display the difference between the boundary conditions Figure 2 shows a simulation with two Gaussian pulses going in opposite directions in a domain with the boundary conditions in equation (18). Here, it can be observed that the left pulse is passing through the boundary, while the right pulse is reflected.
0 500 1000 1500 2000 2500 3000 3500 4000 4500 x
0 0.5 1 1.5 2 2.5
t
Figure 2: The time evolution of two Gaussian pulses moving in opposite direc- tions. The domain has characteristic boundary conditions on the left boundary and a wall condition to the right.
4 The 1D Problem
In this section, the linearized 1D Euler equations:
qt= −A(x)qx− E(x)q + f, a ≤ x ≤ b, t > 0,
Laq = ga, x = a, t ≥ 0,
Lbq = gb, x = b, t ≥ 0,
q = F, a ≤ x ≤ b, t = 0.
(19)
are analyzed. Firstly, well-posedness for the continuous problem is shown using the energy method. Further, a semi discrete finite difference scheme will be presented and a stability analysis of the scheme that mimics the continuous analysis is performed.
4.1 Continuous Analysis
In this section well-posedness of (19) is shown. A problem is well-posed if it has an existing and unique solution that depends continuously on the problem’s data (initial condition, boundary condition and forcing) [9]. In this report it is assumed that the problem has an existing and unique solution. To show that the problem depends continuously on data the following definition of a well-posed initial-boundary value problem is used,
Definition 1 An initial-boundary value problem is well-posed if for every f ∈ C∞ there exists a smooth and unique solution such that,
||u(·, t)|| ≤ Keα(t−t0)(||F(·)|| + max
t0≤τ ≤tkf (·, τ )k),
where α and K are constants independent of the forcing function f , the initial condition F and t0.
The energy method will be employed to show well-posedness. Before continuing with the analysis the following inner product is introduced.
Definition 2 Let u = [u(1), u(2), ..., u(m)]T and v = [v(1), v(2), ..., v(m)]T be real valued vector functions and let [u,v] ∈ L2[a, b]. The inner product is defined as (u,v) =Rb
a u · v dx
Lemma 1 The 1D linear Euler equations,
qt= −A(x)qx− E(x)q + f, a ≤ x ≤ b, t > 0,
Laq = ga, x = a, t ≥ 0,
Lbq = gb, x = b, t ≥ 0,
q = F , a ≤ x ≤ b, t = 0.
are well-posed according to definition 1 if the boundary conditions are chosen so that
−qTAq|x=bx=a≤ 0.
Proof. To obtain the energy estimate equation (19) is multiplied by qT from the left,
qTqt= −qTAqx− qTEq + qTf.
By using the chain rule, (Aq)x= Axq + Aqxand integration by parts,
(q, qt) =1
2((Aq)x, q) +1
2(qx, Aq) +1
2(q, Axq) − qTAq|x=bx=a
− (q, Eq) + (q, f)
(20)
is obtained, which has the transpose
(qt, q) = −1
2(qx, Aq) +1
2((Aq)x, q) + 1
2(qT, Axq)
− (q, ETq) + (f, q).
(21)
Adding (20) and (21) gives, d
dt||q||2= (q, (Ax− E − ET)q) + (f, q) + (q, f) − qTAq|x=bx=a. If the boundary conditions are chosen so the term,
− qTAq|x=bx=a≤ 0, (22)
the energy estimate d
dt||q||2≤ (q, (Ax− E − ET)q) + (f, q) + (q, f) (23) is obtained. Since Ax and E are known and bounded the energy estimate in (23) can be rewritten as
d
dtkqk2≤ 2αkqk2+ 2kqkkf(·,t)k. (24)
where,
2α = max
x |Ax− E − ET|.
Employing the chain rule and dividing (24) by 2kqk the estimate,
d
dtkqk ≤ αkqk + kf(·, t)k is given. By introducing the variable change u = e−αtkqk,
eαtdu dt + αu
≤ αeαtu + kf(·, t)k is obtained, which gives
du
dt ≤ e−αtkf(·, t)k.
By integrating over time,
ku(t)| ≤ ku(0)k + Z t
0
e−ατkf(·, τ )kdτ, is obtained. Changing the variables back to q gives
||q|| ≤ eαt kF k + Z t
0
e−ατkf(·, τ )kdτ
≤ Keαt kF k + max
0≤τ ≤tkf(·, τ )k, (25)
where K is a constant. Hence the problem is well-posed according to Definition
1.
4.2 Boundary Analysis
In earlier analysis the boundary conditions have been assumed to give a well- posed energy estimate, this means that the condition,
− qTAq|x=bx=a≤ 0, (26)
must hold. In this section it is shown that the two boundary conditions used in this project gives well-posed energy estimates.
To state the non-reflecting characteristic boundary conditions the coefficient matrix A will be split into two parts A = A++ A−. Since A is symmetric it can be diagonalized as A = VADAVA−1where DA,is a diagonal matrix with the eigenvalues of A on the diagonal. By using this the matrix A is split as
A± = VA DA± |DA| 2
VA−1, where
VA=
−√
γ − 1 √γ−11 √γ−11
0 √1
2 −√1
2
√1 γ
qγ−1 2γ
qγ−1 2γ
, and DA=
u 0 0
0 u + c 0
0 0 u − c
.
The characteristic boundary conditions can now be stated as A+q = 0, x = a
A−q = 0, x = b. (27)
With these boundary conditions the boundary term can be written as,
−qTAq|x=bx=a= − qTAq|x=b+ qTAq|x=a
= −(qTA+q + qTA−q)|x=b+ (qTA+q + qTA−q)|x=a
= −qTA+q|x=b+ qTA−q|x=a ≤ 0.
This is negative semi definite since A+is positive semi definite and A−negative semi definite.
The wall boundary condition can only be set on boundaries with exactly one outgoing characteristic variable. This is true at both boundaries when u0= 0. In the case when u06= 0this holds at the boundary where the wind is going out of the domain. So when 0 < u0< c0the reflecting boundary condition can only be set on the right boundary while when −c0< u0< 0it can only be imposed on the left boundary. The boundary analysis is the same in both of these cases.
Hence, the analysis will be performed for the case where 0 ≤ u0 < c0 with a wall boundary condition at the right boundary and a characteristic condition at the left boundary as
A+q = 0, x = a
0 1 0 q = 0, x = b. (28)
With these boundary conditions the boundary term can be written as
−qTAq|x=bx=a = − qTAq|x=b+ qTAq|x=a
= −qTAq|x=b+ (qTA+q + qTA−q)|x=a
= −qTAq|b+ qTA−q|x=a
≤ qTAq|x=b.
(29)
By using a variable change to the characteristic variables w = [wT+ wT−]by the transform q = VAw equation (29) can be rearranged to
− qTAq|x=bx=a≤ w+T−|u0| 0
0 |c0− u0| − |c0+ u0|
w+|x=b. (30) This expression is negative definite for all 0 ≤ u0 < c0. This is the case for all situations modelled in this project. Hence both boundary conditions fulfils the inequality in (26) and therefore a well posed problem is given with these boundary conditions.
4.3 Semi-discrete Stability Analysis
In this section a stable semi-discrete finite difference scheme for the 1D linear Euler equations (3) will be introduced. To ensure stability finite difference operators with the Summation By Parts (SBP) property will be used. Further the boundary conditions will be imposed with a technique called Simultaneous- Approximation-Term (SAT ) which imposes the boundary conditions weakly.
This means that the boundary points are not forced to an exact value at the boundaries, instead a penalty term is added that penalizes any deviation from the imposed boundary conditions. By utilizing the properties of the SBP-SAT finite difference scheme stability will be shown with an estimate that mimics the continuous energy estimate.
4.3.1 Definitions
Before the semi discrete scheme is presented some notations and definitions must be introduced.
The computational domain x ∈ [a, b] is discretized by N + 1 grid points as
x = hi, i = 0, 1, 2, ..., N, h = 1 N. The discrete solution vector at each time is the m(N + 1) vector
q = [q0(1), q1(1), ..., qN(1), q0(2), ..., qN(2), q0(3), ..., qN(2)]T.
To mimic the continuous case a discrete inner product between two vectors u, v ∈ Rm(N +1) is defined by (u, v)H = uTHv with the implied norm ||u||H = uTHu. The matrix H must be Hermitian and positive definite.
The following m(N + 1) vectors will be used frequently e0= [1, 0, ..., 0]T and eN = [0, ..., 0, 1]T. Further the Kronecker product defined as,
C ⊗ D =
c0,0D . . . c0,q−1D ... ... ...
cp−1,0D . . . cp−1,q−1D
, where C is a p × q matrix and D a m × n matrix, will be used.
Let In denote the n × n identity matrix. To allow a more compact notation, let any operator D be extended as
D = Im⊗ D.
Further, let the coefficient matrices be extended as
A = A ⊗ IN +1 and, E = E ⊗ IN +1. (31) Finally, the notation
v0= eT0v, vN = eTNv, will be used for any (N + 1) vector v.
Definition 3 A pth-order difference operator D1= H−1(Q −12e0eT0 +12eNeTN) approximating δ/δx is a pth-order accurate first derivative SBP operator if H = HT > 0 and Q + QT = 0.
4.3.2 Discretization
Using the SBP operator the semi discrete version of (19) is obtained, qt= −AD1q + Eq + f, t > 0,
Laq0= ga, t ≥ 0, LbqN = gb, t ≥ 0,
q = F, t = 0.
(32)
Using SAT technique to impose the boundary conditions gives the finite differ- ence scheme,
qt= −AD1q + Eq + f + H−1e0τ0(Laq0− ga)
+ H−1eNτN(LbqN − gb), (33) where τ0and τN are the penalty parameters which can chosen to obtain a stable estimate.
4.3.3 Stability Analysis
The matrix E, boundary data, initial data and the forcing function will be put to zero in the stability analysis since they will not impact the stability of the problem [9]. The following definition of stability will be used:
Definition 4 A semi-discrete approximation is stable if there are constants K and α such that for all t0, f and F
kq(t)k ≤ eα(t−t0)(kF k + max
t0≤τ ≤tkf k). (34) In the analysis below a variable change to the characteristic variables has been performed on the boundary terms as
w+0= Raw−0 t ≥ 0,
w−N = Rbw+N t ≥ 0, (35)
where
Ra= −(LaVA+)−1(LaVA−), Rb= −(LbVA−)−1(LbVA+).
Here the matrix VA has been split into the two blocks VA= [VA+ VA−]corre- sponding to the characteristics.
Lemma 2 The semidiscrete finite difference scheme in (33) is stable if the penalty parameters
τ0= VA−|DA+| 0
VA−1 and τN = VA
0
−|DA−|
VA−1,
are imposed and the boundary conditions are chosen as in (27) or (28).
Proof. Multiplying (33) by qTH and adding the transpose leads to the semi discrete energy estimate
d
dt||q||2H=q0TAq0− qTNAqN + q0Tτ0Laq0+ qNTτNLbqN. It is possible to choose the penalty parameters in different ways. Here,
τ0= VA
−|DA+| 0
VA−1 and τN = VA
0
−|DA−|
VA−1,
are chosen which, in terms of characteristic variables gives the energy estimate d
dt||q||2H = − w0+T |DA+|w0+− w0−T |DA−|w0−+ 2w0+T |DA+|Raw0−
− wN −T |DA−|wN −− wN +T |DA+|wN ++ 2wTN −|DA−|RbwN +. After some manipulations,
d
dt||q||2H= − (w0+− Raw0−)T|DA+|(w0+− Raw0−) + wT0−(RTa|DA+|Ra− |DA−|)w0−
− (wN −− Rbw0+)T|DA−|(wN −− RbwN +) + wTN +(RTb|DA−|Rb− |DA+|)wN +,
(36)
is obtained. For the characteristic boundary conditions in (27) the energy esti- mate
d
dt||q||2H = − (w0+− Raw0−)T|DA+|(w0+− Raw0−)
− (wN −− Rbw0+)T|DA−|(wN −− RbwN +) ≤ 0
(37)
is given. This is stable according to Definition 4. For the boundary conditions in (28) the semi discrete energy estimate is
d
dt||q||2H = − (w0+− Raw0−)T|DA+|(w0+− Raw0−)
− (wN −− Rbw0+)T|DA−|(wN −− RbwN +) + wN +T (−|u0| 0
0 |c0− u0| − |c0+ u0|
⊗ IN +1) wN +≤ 0,
(38)
which negative and for all 0 ≤ u0 < c0. Hence the numerical scheme is stable
for the boundary conditions in (27) and (28).
Figure 3: Computational domain Ω in 3D
5 The 3D Problem
In this section the 3D linearized Euler equations (10) are analyzed. The 3D analysis is similar to 1D since each direction is treated separately in the same way as in the 1D case. In this case there are 6 boundaries, 2 per dimension.
They are denoted W (west), E (east), S (south), B (bottom) and T (top) as shown in Figure 3. The 3D domain ΩE,N,TW,S,B is W < x < E, S < y < N and B < z < T. To simplify the 3D analysis the background field will be considered constant and data put to zero. The 3D problem is stated as
qt+ Aqx+ Bqy+ Cqz+ Eq = f q ∈ Ω t > 0,
q = F t = 0 (39)
with the general boundary conditions
LWq = gW q ∈ W, LEq = gE q ∈ E, LSq = gS q ∈ S, LNq = gN q ∈ N, LBq = gB q ∈ B, LTq = gT q ∈ T.
(40)
5.1 Continuous Analysis
In this section well-possessedness is shown for the linearized 3D Euler equations (39) with general the boundary conditions in equation (40). In the analysis the following inner product is used.
Definition 5 Let the vector functions u,v ∈ L2[Ω1,1,10,0,0], where u = [u(1), u(2), ..., u(m)]T and v = [v(1), v(2), ..., v(m)]T. Further, let the inner product be defied as (u,v) = R1
0
R1 0
R1
0 uTv dxdydz, with corresponding norm ||u||2= (u,v).
Multiplying (39) by qT and putting data and E to zero gives,
qTqt= −qTAqx− qTBqy− qTCqz, (41)
which has the hermitian transpose
qTtq = −qTxAq − qTyBq − qTzCq. (42) By adding (41) and (42), integrating over the domain and assuming homoge- neous boundary data the energy estimate
∂
∂t||q||2= BTW + BTE+ BTS+ BTN+ BTB+ BTT, (43) is obtained where
BTW = Z
W
qTAqdS, BTE= Z
E
qTAqdS, BTS =
Z
S
qTBqdS, BTN = Z
N
qTBqdS, BTB =
Z
B
qTCqdS, BTT = Z
T
qTCqdS.
By treating these boundary terms as in the 1D case, it can be shown that this estimate is well posed.
5.2 Semi-discrete Stability Analysis
In this section the 3D SBP operators are defined and stability for the 3D-problem is shown. The continuous domain is discretized by (Nx+ 1)(N y + 1)(N z + 1) points where
(xi, yj, zk) = ( i Nx, j
Ny, k
Nz), i = 1, 2..., Nx, j = 1, ..., Ny, k = 1, ..., Nz. The approximate numerical solution at each time step is represented by the (Nx+ 1)(Ny+ 1)(Nz+ 1)mvector
q = [q000(1), q(1)001..., q00N(1)
z, q(1)010, q..., q0N(1)
yNz, ..., q(1)N
xNyNz, ..., q(m)N
xNyNz]T. The matrices EW, EE, ES, EN, EB and ET are designed so the vectors E∗q contains the elements from corresponding boundary. The difference operator D1xis an extension of the 1D operator D1and approximates the derivative ∂/∂x.
Further Hx is the extension in x-direction of the discrete norm H. Analogous D1y, Hy, D1zand Hzare the derivatives and norms in corresponding directions.
Dx= Im⊗ D1⊗ INy⊗ INz, Hx= Im⊗ H ⊗ INy ⊗ INz, EW = Im⊗ e0x⊗ INy⊗ INz, EE = Im⊗ eN x⊗ INy⊗ INz, Dy= Im⊗ INx⊗ D1⊗ INz, Hy= Im⊗ INx⊗ H ⊗ INz, ES = Im⊗ INx⊗ e0y⊗ INz, EN = Im⊗ INx⊗ eNy⊗ INz, Dz= Im⊗ INx⊗ INy ⊗ D1, Hz= Im⊗ INx⊗ INy ⊗ H, EB = Im⊗ INx⊗ INy ⊗ e0z, ET = Im⊗ INx⊗ INy. ⊗ eNz, Hxy= HxHy Hxz = HxHz
Hyz= HyHz Hxyz= HxHyHz
