Stretched grids as non-reflecting boundary conditions aim at gradual under-resolution and dissipation of sinusoidal waves and are computationally costly and ineffective if not implemented properly. This defines a need for insight to enable efficient implementation. Two methods of predicting the behavior of sinusoidal waves propagating through stretched grids are presented, an initial value and a boundary value approach. The initial value problem is an approximation and has the advantage that it can be applied with many higher order accurate finite difference schemes. The boundary value problem is more exact but is not applicable for higher order methods. Both methods are verified using a two-dimensional linearized Navier-Stokes solver, which is implemented in the framework of this study. Further, acoustic black hole layers are introduced as an alternative to stretched grids. Using both methods in combination reduces the number of grid points necessary in a boundary domain drastically. All methods are improved with respect to available parameters and their performance is compared using benchmark flows. Explicit recommendations for implementation are developed.
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 underresolution 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.
-
-
-
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X.