Analyzing stretched grids and introducing acoustic black hole layers as non-reflecting boundary conditions

Contents

1 Introduction
1.1 Earlierinvestigations
1.2 Objectives and structure of the thesis

2 Governing equations and finite difference techniques
2.1 LinearizedEulerEquations
2.1.1 Continuousdifferential equations
2.1.2 Boundary conditions forthe governing equations
2.2 Finitedifferencemethods
2.2.1 Finitedifference scheme
2.2.2 1D scheme
2.2.3 Deriving coefficients for variable grid spacing
2.2.4 Time integration
2.2.5 Artificial damping and viscosity
2.2.5.1 Viscosity
2.2.5.2 Selective artificial damping
2.2.5.3 Deriving a new exact method to obtain SAD coefficients
2.2.5.4 Filtering
2.3 Wave propagation through finite difference schemes
2.3.1 Dispersion relation (modified wavenumber)
2.3.2 Growing or decaying solutions
2.3.3 Single-sided schemes for boundaries
2.3.4 Growing or decaying solutions through time integration schemes
2.3.5 Theshortwavecomponent

3 Wave propagation through stretched grids
3.1 Dispersion relation for unequal grid spacing
3.2 Phasevelocity
3.3 Groupvelocity
3.4 Reflection
3.5 Cut-offfrequency

4 Non-reflecting boundary conditions
4.1 Categorizing non-reflective boundary conditions
4.2 Thompson Boundary conditions
4.2.1 Deriving general equations
4.2.2 Non-reflecting subsonic outflow BC
4.2.3 Non-reflective inflow boundary conditions
4.2.4 Discussion and limitations of the Thompson approach
4.3 Stretchedgrids
4.3.1 Initialvalue problem
4.3.2 The boundary value problem
4.4 Acoustic black hole layers
4.4.1 Acousticblackholes
4.4.2 Deriving acoustic black holes for finite difference schemes
4.4.3 Numericalscheme
4.4.4 Stretched grid analogy
4.4.5 Hybrid and optimized BHL

5 Benchmark and optimization techniques
5.1 Benchmarktechniques
5.1.1 Testflows
5.1.1.1 Planewave
5.1.1.2 Gausspulse1D
5.1.1.3 Gausspulse2D
5.1.1.4 Obliquewave
5.1.1.5 Vortex
5.1.2 Measuring errors
5.2 Optimizationtechniques
5.2.1 Optimizationin1D and2D
5.2.2 Matlab optimization functions

6 2D solver
6.1 Generalsetupandcapability
6.2 Boundary conditions and initial values
6.3 Numericalscheme
6.4 Implementation
6.5 Discussion and limitations
6.6 1D solver

7 Verifying analytical expressions numerically
7.1 Stretchedgrids
7.1.1 Amplitude decay predicted through initial value problem (1D)
7.1.2 Amplitude decay predicted through initial value problem (2D)
7.1.3 The boundary value problem (1D)
7.1.4 The boundaryvalue problem (2D)
7.2 Black hole layers
7.2.1 Numericalproofoftheory
7.2.2 Stretched grid analogy
7.3 Conclusion:Analysis

8 Optimizing boundary conditions
8.1 Parameterstudy
8.1.1 Influence ofgrid stretching on reflections
8.1.2 Influence ofdamping and viscosity
8.1.3 Reflections from the black hole layer
8.2 Optimizing black hole layers
8.2.1 Optimization in one space dimension
8.2.2 Optimization in two space dimensions
8.3 3-zone acoustic black hole layer
8.4 Conclusion: Improved implementation

9 Benchmarkcalculations
9.1 Gauss pulse, no mean stream
9.2 Vortex with mean stream
9.3 Conclusion:Benchmark

10 Conclusion
10.1 Methodologicalsummary
10.2 Recommendations for implementation
10.3 Aspectsforfutureresearch

Appendix A Additionalmaterial
A.1 Finite difference coefficients
A.2 Amplitude decay through initial value problem

Appendix B Matlab code
B.1 solver_1D.m
B.2 optimize_main.m
B.3 run_solver_2D.m
B.4 solver_2D.m

Appendix C Literature

List of figures

Figure 2.1 Grid nodes of an N+M+1 stencil

Figure 2.2 Distribution of dl with respect to grid nodes

Figure 2.3 Fourier transform for 3-, 5- and 7-stencils

Figure2.4 Fourier transform for the new damping approach comparing to Tam’s optimized coefficients

Figure 2.5 Half logarithmic plot of the Fourier transform for the new damping approach comparing to Tam’s optimized coefficients

Figure 2.6 Dispersion relation of different numeric schemes

Figure 2.7 Dispersion relation for one-sided stencils

Figure 2.8 Temporal dispersion relation for different integration schemes

Figure 2.9 Stability of different integration schemes

Figure 3.1 Dispersion relation in a stretched grid with varying difference schemes

Figure 3.2 Dispersion relation in a stretched grid with varying stretching ratio

Figure 3.3 Dispersion relation in a stretched grid including viscosity

Figure 3.4 Group velocity in a stretched grid with varying difference schemes. The respective cut-offfrequency is indicated

Figure 3.5 Group velocity in a stretched grid with varying stretching ratio

Figure 3.6 Distribution of ElH with respect to grid nodes

Figure 4.1 Plate with an acoustic black hole attached [13]

Figure 5.1 Sinusoidal wave

Figure 5.2 Gauss pulse (1D) with a half width of one

Figure 5.3 Gauss pulse (2D)

Figure 5.4 Oblique wave with inflow, outflow and periodic boundary conditions

Figure 5.5 Distribution ofthe circumferential velocity

Figure 5.6 Vortex with a counter clockwise rotation leaving the domain

Figure 5.7 Setup for measuring the error or optimizing parameters within the boundary domain

Figure 6.1 Setup of the 2D solver

Figure 6.2 Finite difference operator

Figure 7.1 Amplitude decay predicted through the initial value problem method with varying damping

Figure 7.2 Amplitude decay predicted through the initial value problem method with varying grid stretching

Figure 7.3 Amplitude decay predicted through the initial value problem method with varying order of accuracy

Figure 7.4 Amplitude decay predicted through the initial value problem method using viscosity

Figure 7.5 Oblique waves with varying angle of incidence (0°, 20°, 40° and 60°)

Figure 7.6 Amplitude decay for oblique waves with varying incidence and constant mean stream

Figure 7.7 Analytical solution predicted with the boundary value problem method with varying viscosity and grid stretching s=1.03

Figure 7.8 Analytical solution predicted with the boundary value problem method with varying viscosity and grid stretching s=1.05

Figure 7.9 Analytical solution predicted with the boundary value problem method with varying viscosity and grid stretching s=1.1

Figure 7.10 Analytical solution predicted with the boundary value problem method with varying viscosity and grid stretching s=1.2

Figure 7.11 Analytical solution predicted through the boundary value problem method with mean stream and varying oblique incident angle

Figure 7.12 Proofofthe lineardependence ofthe time step and the wavelength

Figure 7.13 BHL with a stretched grid analogy

Figure 8.1 Cost function over grid stretching for various finite difference schemes

Figure 8.2 Cost function over damping for the 2nd and 4th order scheme and several grid-stretching factors

Figure 8.3 Cost function over viscosity for the 2nd and 4th order scheme and several

grid-stretching factors

Figure 8.4 Comparing different damping stencils with s = 1.05 and 4th order difference scheme

Figure 8.5 Cost function over grid stretching comparing stretched grids and black hole layers

Figure 8.6 Results of an optimization in 1D ofa black hole layer with 8 points

Figure 8.7 Results of an optimization in 2D of a black hole layer with 8 points

Figure 8.8 Results of an optimization in 2D of a black hole layer with 12 points

Figure 8.9 Zonal approach for an 8-point BHL

Figure 8.10 Parameter study for the 8-point BHL

Figure 8.11 Parameter study for the 12-point BHL

Figure8.12 Parameterstudyforthe 16-pointBHL

Figure 9.1 Comparing stretched grids with varying grid stretching for the 2D Gauss pulse benchmark

Figure 9.2 Comparing black hole layers with the stretched grid analogy and hybrid layers for the 2D Gauss pulse

Figure 9.3 Thompson approach for the 2D Gauss pulse

Figure 9.4 Comparing different black hole layers with the 3-zone approach

Figure 9.5 Comparing a stretched grid, hybrid layer and black hole layer with s = 1 as well as a 16-point 3-zone BHL for the 2D Gauss pulse benchmark

Figure 9.6 Comparing a stretched grid, hybrid layer and black hole layer with s = 1 as well as a 8-point 3-zone BHL for the 2D Gauss pulse benchmark

Figure 9.7 Spatial error distribution for different boundary conditions at four different times

Figure 9.8 Comparing stretched grids with varying grid stretching for the vortex benchmark

Figure 9.9 Comparing black hole layers with the grid stretching analogy and hybrid layers for the vortex benchmark

Figure 9.10 Thompson approach for the vortex benchmark

Figure 9.11 Comparing a stretched grid, hybrid layer and black hole layer with s = 1 for the vortex benchmark

Figure9.12 Spatial error distribution for different boundary conditions at four times The vortex propagates from the middle to the right

List of tables

Table 2.1 Coefficientsfortime integration

Table 2.2 Coefficientsfordifferentdamping schemes

Table 2.3 Damping coefficients obtained with the new approach

Table 5.1 Data for oblique waves with different incident angles

Table 6.1 Available boundary conditions for the 2D solver

Table 6.2 Available conditions for initial values for the 2D solver

List of symbols

illustration not visible in this excerpt

1 Introduction

Especially in urban and suburban areas, noise pollution has become one of the key objections for infrastructural projects and few external circumstances have such a negative influence on human well-being and the quality of living. Regulations are becoming more and more strict with regard to noise pollution and large measures are taken to protect human hearing organs at working places.

Efforts of reducing noise require accurate methods of predicting sources and propagation and convection of sound. One of the challenges in computational aero-acoustics are very small amplitudes compared to large flow gradients in turbulent flows and the short wave lengths of acoustic waves compared to typical propagation distances. Both result in the necessity of using highly accurate numerical methods. In the majority of acoustic solvers, finite difference methods are used, since they are easily extended to higher order accuracy. [1]

Besides the requirement of high accuracy, one of the most critical aspects of CAA is the implementation of robust and stable boundary conditions. Especially free field as well as in- and outflow boundary conditions are challenging. Inflow boundary conditions should allow for defining incoming waves at the boundary and outgoing waves to leave the domain without causing reflections or inducing flows. Free field boundary conditions should ideally not reflect any outgoing waves. Outflow boundary conditions should be capable of not reflecting any waves while flows such as strong flow gradients or vortices leave the domain. [1]

In common CFD solvers (e.g. Star-CCM+, Ansys), so-called Dirichlet boundary conditions are often used, which specify the correct number of variables at the boundaries such that the equations are well-posed. At the inlet, the velocity can be specified while at the outlet often the pressure is implemented as boundary condition. These conditions are however reflecting acoustic waves.

Many different inflow, outflow and free field boundary conditions have been developed, for example linearized boundary conditions such as Giles boundary conditions, non-linear Thompson boundary conditions, different kinds of buffer zones such as PMLs as well as the use of grid stretching. The latter method aims at stretching the grid close to the boundary in such a way that outgoing waves are gradually under-resolved and eventually dissipated, making non-reflecting boundary conditions in principle unnecessary. Therefore, stretched grids can also be implemented to obtain non-reflecting conditions in solvers where only Dirichlet boundary conditions are available. Stretched grids are on the one hand popular among acoustic engineers, since they are easy to use and obviate the need for implementing modifications of the principle finite difference code, but on the other hand they are inefficient if not used properly. This thesis aims at analyzing, developing and improving stretched grids as non-reflecting boundary conditions.

The conflict when implementing stretched grids as boundary conditions is that grid stretching factors should be small to minimize reflections caused by the stretching. Small grid stretching ratios require however large and computationally costly boundary zones. This defines a need for insight and optimization of available parameters to be able to implement the computationally most efficient method with a given accuracy. The results of this thesis are prepared fordirect implementation in a numerical code.

Another objective of this study is to introduce an alternative implementation of free field boundary conditions, called black hole layer, since it is based on the principle of acoustic black holes. This technique has been applied on physical plates to reduce the bending waves reflected from the edges by reducing the wave propagation speed locally inside the plates. This method is transferred to the finite difference method by reducing the time step locally inside a boundary domain and aims similarly to stretched grids at gradual under­resolution of sinusoidal waves.

1.1 Earlierinvestigations

Kreiss and Efraimsson [2] presented two approaches for analyzing the behavior of sinusoidal waves entering stretched grids. The first one is an initial value problem, which is equal to the one that is commonly used to calculate the dispersion relation of a given finite difference scheme. The other approach is a boundary value problem equaling a Fourier transform in time. Both techniques have successfully been used to predict the amplitude decay of waves propagating through stretched grids.

The analysis in [2] is an extension of the works done by Vichnevetsky [3],[4] who has developed a method for analyzing the first order semi-discretized advection equation. He presented also a way of calculating reflections resulting from non-equally spaced grids, which has been used in this study.

1.2 Objectives and structure of the thesis

This thesis aims at analyzing stretched grids as boundary conditions using the boundary value approach presented by Kreiss and Efraimsson and an extended initial value approach. Both techniques are verified using a numerical code.

Further, the concept of acoustic black hole layers is derived for the finite difference method and theiropportunities and limitations as non-reflecting boundary conditions are discussed.

The code used for the numerical study is developed in the framework of this thesis. The used linearized Navier-Stokes equations and necessary finite difference techniques are described in chapter 2. This chapter introduces also the main characteristics necessary to describe waves propagating through finite difference schemes such as the dispersion relation and group velocity.

The latter analysis is extended to variable grid spacing in chapter 3, since these quantities are necessary to describe waves travelling through stretched grids. A technique is presented that can in contrast to preexisting methods be applied to many higher order accurate finite difference schemes.

Chapter 4 applies these quantities in a method to predict the amplitude decay of waves travelling through stretched grids. The boundary value method is also presented. Further, an overview of existing non-reflecting boundary conditions is given and equations for Thompson boundary conditions are derived for the current solver. Equations describing the black hole layer are derived.

As mentioned previously, the analytic correlations are to be verified in a numerical code. Appropriate benchmark flows are described in 5.

The capabilities and basic characteristics of the 2D solver developed are outlined in chapter 6.

The following chapters apply the numeric code in different applications. Chapter 7 verifies the correlations describing the amplitude decay and black hole layer theory. Chapter 8 presents a numerical study analyzing the amount of reflections obtained for different finite difference techniques. In addition, black hole layers are optimized numerically with respect to grid stretching and time scaling.

The number of presented non-reflecting boundary conditions promotes the question of which of them is performing best under challenging conditions found in practice. Chapter 9 presents benchmark calculations for a two dimensional vortex and Gauss pulse evaluating all approaches discussed.

Finally, the conclusions of this study are presented in chapter 10.

2 Governing equations and finite difference techniques

As part of this thesis, a two-dimensional acoustic solver was implemented. The governing equations used in the solver are presented (section 2.1.1) in chapter 2.1. General requirements for boundary conditions are outlined in paragraph 2.1.2. Non-reflecting boundary conditions are discussed in chapter 4.

In addition to the governing equations, the numerical scheme (sections 2.2.1 and 2.2.2) and further finite difference methods are necessary for implementing a numerical code. Finite difference parameters are derived for variable grid spacing in section 2.2.3. Time integration is discussed in section 2.2.4. Further finite difference techniques are damping, viscosity and filtering (2.2.5.2-2.2.5.4), which play an important role in reducing the amplitude of waves within a boundary domain and are necessary to dissipate the short wave component of finite difference schemes (2.3.5). In addition to existing techniques, a new, improved implementation of selective artificial damping is derived. All presented techniques are implemented in the one- and two-dimensional acoustic solver developed in the framework of this thesis.

Wave propagation through finite difference schemes is one of the central aspects of this thesis. In order to examine the characteristics of propagating waves for constant grid spacing, the dispersion relation of discrete schemes is discussed in 2.3.1. The dispersion relation is also a tool for investigating stability (2.3.2) and the behavior of one-sided schemes, which are required for boundary closures (section 2.3.3).

2.1 Linearized Euler Equations

The governing equations as well as suitable boundary conditions are outlined in the following.

2.1.1 Continuous differential equations

The governing differential equations used in the solver are shown in equations (2.1) to (2.3) as described in Gustafsson, Kreiss and Oliger [5]. They are: Continuity equation:

illustration not visible in this excerpt

Momentum equations:

illustration not visible in this excerpt

The Euler equations are shown here and viscosity terms will be included later in the context of different damping implementations. In addition to the governing equations, another condition is required for closure. We assume a calorically perfect gas, yielding that the isentropic state condition can be used as relation between density and pressure:

illustration not visible in this excerpt

The equations may be linearized around a mean flow with the velocity u0 in x-direction due to the small magnitudes of perturbations. The ansatz is made, where the exact variable is indicated with an upper case letter. The variables u, v and p are defined as small perturbations of the mean flow u0 and the mean pressure p0. Further, it is assumed that the mean flow in y-direction is zero (v0 = 0). These simplifications lead to the following set of equations and are basis for the numerical finite difference scheme

illustration not visible in this excerpt

2.1.2 Boundary conditions for the governing equations

Equation (2.5) is an initial-boundary-value problem and requires being well-posed in order to have a unique solution [6]. Requirements for different flow conditions such as inflow, outflow and wall boundary conditions are presented in [5]. According to this analysis, the number of boundary conditions must be equal to the number of negative eigenvalues of the coefficient matrices to achieve well-posed semi bounded conditions. We have for Supersonic inflow, c < u0: All variables have to be specified at the boundary, for example

illustration not visible in this excerpt

Subsonic inflow, 0 < u0 < c: Two variables have to be specified, for example using a forcing function f, which is dependent on y and t:

illustration not visible in this excerpt

Neitherof these conditions represents a non-reflective boundary condition. An appropriate choice to avoid reflections is presented in 4.2. Rigid wall, u0 = 0: One condition has to be specified, which naturally is the normal velocity component:

illustration not visible in this excerpt

Another possibility to model walls is the use of ghost points [6]. Both approaches have been used in the current study.

Subsonic outflow, —c < u0 < 0: One boundary condition is necessary. The choice of this condition determines if the outflow is reflective. A suitable choice for non-reflecting boundary conditions is presented in chapter4.2.

Supersonic Outflow, u0 < —c: No boundary conditions have to be specified.

2.2 Finite difference methods

Implementing a code requires a numerical scheme as well as explicit techniques and coefficients, e.g. for random grid spacing, artificial damping and viscosity. The finite difference methods used in the framework of this thesis are presented in the following paragraph.

2.2.1 Finitedifferencescheme

Starting from equation (2.5), an explicit finite difference scheme is developed. A similar scheme has for example been used by Tam and Webb [7] to analyze their radiation boundary conditions. The scheme consists of a multi-step time integration and can be combined with many types of spatial discretization methods, for example with the fourth order 7-stencil DRP scheme developed by the same authors. Another method for deriving coefficients for spatial discretization is shown in chapter 2.2.3.

From equation (2.5), the following abbreviations are introduced:

illustration not visible in this excerpt

Considering the matrix multiplication, E and F are both column vectors. With these definitions, the governing equations can be split into the two equations

illustration not visible in this excerpt

The index l represents the grid points in x-direction, the index m the grid points in y-direction and n the respective time step. The variables N and M define the size of the stencil of the finite difference scheme and are visualized in Figure 2.1. The variable R indicates the number of steps involved in the multi-step time integration method. The initial conditions for this scheme are (see [6]):

illustration not visible in this excerpt

For internal stability considerations, it is referred to Tam and Webb [7], who analyzed a similar scheme.

Close to the boundaries, a one-sided scheme may be used to implement the conditions presented in chapter2.1.2.

2.2.2 IDscheme

Besides the two-dimensional scheme, a simpler and more efficient one-dimensional code is desired. The above differential equations and numeric scheme may easily be simplified to a one-dimensional case. The corresponding differential equations are

As for the two-dimensional case, the matrix expression is abbreviated with

illustration not visible in this excerpt

The initial conditions forthis scheme are equivalent to the ones ofthe two-dimensional code.

2.2.3 Deriving coefficients for variable grid spacing

Changing the grid spacing inside the buffer zone is one of the central points of this thesis. General equations for calculating finite difference coefficients are therefore derived in the following. A similar derivation for constant grid spacing is presented in [6]. In more general numeric codes, variable grid sizes are usually implemented using curvilinear coordinates [1]. The direct formulation of coefficients for variable grid spacing simplifies analysis of the numeric scheme and is therefore preferred in the context ofthis study.

A common ansatz for calculating the derivative of a continuously differentiable function f with respect to x is

illustration not visible in this excerpt

where Ax is the distance between two grid points and aj a set of unknown coefficients. The finite difference stencil ranges from the grid points l - N to l + M and the corresponding grid spacings are defined such as in Figure 2.2. Now, the grid spacing is described with dimensionless coefficients hi, that are defined as the distance between two grid points divided by Ax:

illustration not visible in this excerpt

Figure 2.2 Distribution of d, with respect to grid nodes.

In order to calculate the coefficients aj, equation (2.12) may be expanded in a Taylorseries:

illustration not visible in this excerpt

This linear system of equations is of the dimension N + M + 1 and has to be solved separately for every grid point in forehand. This technique has two main advantages for random grid spacing:

- The accuracy of the finite difference scheme is preserved for non-equally distributed grids.

- No significant negative influence on the computation speed in comparison to an equally spaced grid is expected since the coefficients are calculated once and stored for each grid point.

For a finite difference code requiring several dimensions, this technique is applied to all dimensions independently.

A number of commonly used finite difference stencils is included in appendix A.1.

Tam and Webb [7] optimized a symmetric 7-stencil scheme with 4th order accuracy resulting the following coefficients for constant grid spacing.

illustration not visible in this excerpt

2.2.4 Time integration

The multi-step time integration scheme presented in [7] is considered. The scheme is based on the ansatz [6]

The coefficients bj can be calculated using the following system of equations:

Solutions of this system of equations are included in Table 2.1. Tam and Webb [7] also presented an optimized version of the four-step integration scheme, which is included in the table.

illustration not visible in this excerpt

Table 2.1 Coefficients for time integration.

2.2.5 Artificial damping and viscosity

Damping and dissipation can be implemented using many different techniques. The full Navier-Stokes equations contain viscosity terms that are presented in section 2.2.5.I. Another common approach is selective artificial damping (2.2.5.2 and 2.2.5.3) which is computationally less costly and is often used in finite difference schemes. The aim of adding damping in acoustic solvers is usually to dissipate the short wave component present in the solutions (see 2.3.5). Therefore, only high frequencies should be dissipated. Since common approaches are not satisfactory, a new way of obtaining coefficients for selective artificial damping is presented aiming at reducing the dissipation of well-resolved waves. Finally, the concept of filtering is discussed in the context of stretched grids as boundary conditions (2.2.5.4).

2.2.5.1 Viscosity

In industry, full non-linear Navier-Stokes equations are used frequently for aero-acoustic simulations. These solvers use viscosity terms inherent in the full Navier-Stokes equations to damp out the high frequency content. Although linearized, viscosity terms can be added to the set of equations used in this study, which results in the following changes:

The numerical scheme is accordingly modified with an additional term, which is for the u- component of the equations:

illustration not visible in this excerpt

The velocity component in y-direction changes similarly and the last equation does not feature any viscosity, since it is derived from the continuity equation. The viscosity e is defined such that Ax2 cancels out, which has advantages for the analytical analysis of the behavior of the scheme. The coefficients cj correspond to the second derivative of the respective velocity component and can be calculated with the same coefficient matrix as the first derivatives (2.15). The system of equations is (2.19)

2.2.5.2 Selective artificial damping

Another way of adding damping to the equations is the concept of selective artificial damping (SAD) [6]. The numerical scheme is changed in the following way:

illustration not visible in this excerpt

D. and Dm are weighting factors dependent on the node location. Only the values that are differentiated in the respective direction are taken into account, i.e. u and p in x-direction and v and p in y-direction.

The coefficients dj are commonly obtained from a Fourier transform (for details see [6] and [7]) which describes the frequency range of the damping

illustration not visible in this excerpt

where the damping coefficients are assumed to be symmetric (dj _ d-j). For a three-stencil scheme, the number of unknown damping coefficients is two, for a five-point stencil there are three unknown coefficients and for a seven-point scheme, there are four unknown coefficients. In the following, two conditions are defined using the above Fouriertransform.

Since long waves should not be dissipated, a condition can be defined

illustration not visible in this excerpt

In order to damp out spurious waves most effectively, the Fourier transform should be maximum at the wave number of spurious waves:

illustration not visible in this excerpt

From these two conditions, the damping coefficients for a three point stencil can be calculated directly as d0 _ 0.5 and di _ -0.25.

For a five-point stencil, one condition is missing and the coefficients become d0 _ 1 ~^2 and di _ -0.25.

For certain values of the possible parameters, the Fourier transform becomes negative with a minimum between 0 and n. Since negative damping causes growing solutions, it has to be determined which values d0 and d2 can be taken such that the function does not become negative.

for the position of the minimum. Since no minimum is desired, this expression has to be smaller than -1 or larger than 1:

For the first condition, do should be larger than 0.375 and for the second condition, do should be smaller than 0.625. That is, there is a maximum between 0 and n for values of do larger than 0.625. Several of the resulting curves are plotted in Figure 2.3, of which the upper boundary curve is brown and the lower boundary curve is red.

For the 7-stencil scheme, two conditions for determining the coefficients are missing and the common approach is to fit the Fourier transform to a Gauss function under the constraints given above [6],[7]. A suitable cost function is for example 2

illustration not visible in this excerpt

with the constant 0 < b < n. The latter can be used to improve the quality of the resulting Fourier transform. Results of the optimization are shown in Figure 2.3 where b = 2 was used for all optimizations. It can be seen that the curve for a = 0.2n is partly negative resulting in an amplification instead of damping, which is not desirable. Therefore, these coefficients should not be used.

illustration not visible in this excerpt

Table 2.2 Coefficients for different damping schemes.

In the framework of this thesis, these damping coefficients are used even for stretched grids. Numerical studies showed that damping coefficients calculated using grid spacing information resulted in less effective damping if the grid is stretched. Therefore, the coefficients for constant grid spacing are used also for stretched grids.

2.2.5.3 Deriving a new exact method to obtain SAD coefficients

There are several drawbacks when the coefficients of the Fourier transform are determined by curve fitting to a Gauss function. That is, that the curves fit quite well to the Gauss function for larger values of aAx, but often show a oscillatory behavior around D = 0 for small values of aAx. The second problem is, that the Gauss curves are not decaying quickly enough for aAx going towards zero, resulting in excessive dissipation of the mean flow.

Therefore, better damping coefficients are required for the Fourier transform, which only dissipate the short wave component. A new approach is proposed, which is based on the idea of moving all roots of D in the range of 0 < aAx < n to aAx = 0. Further, all higher derivatives of D should only have roots at the position aAx = 0, which guarantees the quickest and smoothest possible decay of D for aAx tending towards zero. Considering again the Fourier transform

illustration not visible in this excerpt

all uneven higher order derivatives are set to zero, resulting in the following set of equations:

illustration not visible in this excerpt

The even higher order derivatives are zero in aAx = 0 by default due to the symmetric character of the equation. More specific, the derivatives should have a double root at aAx = 0. By default, one of the roots of the derivatives is zero in aAx = 0. Requiring a double root moves the second possible root between 0 < aAx < n to aAx = 0.

For example for a 9-point daming stencil, there are 5 unknowns, d0, d1( d2, d3 and d4. Thus, the first, third and fifth order derivative have to be evaluated. In the following, the coefficients will be derived in detail for this stencil.

For a 9-point stencil, the first derivative is

illustration not visible in this excerpt

Using the addition theorems we have

illustration not visible in this excerpt

Now, the equation is factorized and the wavenumber is set to zero (aAx = 0):

illustration not visible in this excerpt

It can be seen, that the equation has two roots, the sinus term and the sum in the parenthesis. Thus, the slope of the Fourier transform is zero at aAx = 0, independently of the damping coefficients. Since no maxima or minima are desired between 0 and n, the second root is also moved to aAx = 0 and thus

illustration not visible in this excerpt

Since the second derivative is due to the symmetric damping stencil equal to zero, the third derivative can also be calculated and set to zero. This equation reads:

illustration not visible in this excerpt

After applying the above addition theorems, setting the equation to zero and simplification, the final equation is

illustration not visible in this excerpt

The fifth derivative is calculated in the same way and simplified, resulting in

illustration not visible in this excerpt

For smaller systems, the lowest row and last columns of the Matrix may be removed. The resulting coefficients are included in Table 2.3.

illustration not visible in this excerpt

Table 2.3 Damping coefficients obtained with the new approach.

illustration not visible in this excerpt

Figure 2.4 Fourier transform for the new damping approach comparing to Tam’s optimized coefficients.

The coefficients in Table 2.3 are evaluated against the optimized 7-stencil scheme with a = 0.3n in Figure 2.4 and Figure 2.5. It can be seen, that the 7-stencil curve of the new approach is equivalent to the optimized one in the range of short wave components while the curve featuring the new coefficients converges much quicker for small wave numbers. This results in less dissipation for large wavelengths.

A logarithmic plot in Figure 2.5 shows the detailed behavior of the new damping curves in the frequency range of physical waves and helps choosing the right stencil based on maximum damping limits desired.

illustration not visible in this excerpt

Figure 2.5 Half logarithmic plot of the Fourier transform for thenew damping approach comparing to Tam’s optimized coefficients.

2.2.5.4 Filtering

In ordinary flows, the aim of damping is to reduce the short wave length content, i.e. dispersive and parasite waves. As can be seen in Figure 2.3, the presented approach of selective artificial damping dissipates more than necessary. Another method of filtering high frequency content is the use of a Taylor filter, which filters only high frequency content [6].

In the current study, it is desirable that damping dissipates acoustic waves in the buffer layer. As a conclusion, selective artificial damping was used within the buffer zones instead of Taylorfiltering such as described in [6].

In the physical domain, no damping or filtering was used, since the parasite waves reflected from the boundary should be measured and not removed by a filter. In a solver aiming at physical problems, a filtershould be used within the physical domain, however.

2.3 Wave propagation through finite difference schemes

Finite difference schemes do not solve the underlying partial differential equations exactly but are an approximation. The differences of the approximate to the exact solution are discussed in this section for constant grid spacing.

2.3.1 Dispersion relation (modified wavenumber)

A dispersive medium is characterized in general by the dependence of the phase velocity of the frequency [3]. Central finite difference schemes yield dispersive solutions, just as dispersive media. The dispersive terms are especially important in the context of stretched grids, since the wave travelling characteristics are dependent on the resolution of the wave. As a consequence, the relation between frequency and the wavenumber is dependent on the grid spacing. The relation between the frequency and wave number is in this thesis called the dispersion relation and will here be analyzed in some detail due to its importance for the whole topice.