Spectral nonlocal boundary conditions for the wave equation in moving media
Preprint, Inst. Appl. Math., the Russian Academy of Science


Ivan L.Sofronov, Olga V.Podgornova
( .., ..)

Russian Academy of Science, Keldysh Institute of Applied Mathematics
. ..

Moscow, 2004

This work was supported by RFBR grant 03-01-0567
( 03-01-0567)

Abstract

A spectral approach of generating low-reflecting boundary conditions for the wave equation in the moving media is proposed. Operator of boundary conditions is firstly derived in exact form for discrete equations, and then necessary approximation modifications are developed to obtain reasonable computational costs. The sum-of-exponentials representation of occurring temporal kernels is used as a key approach for such modifications.

- . , , . .


Introduction

 

The considered wave equation for the moving media occurs after the formal change of variables

 

(1.1)

 

in the wave equation

(1.2)

 

where is a given constant speed, is the speed of sound, . It reads

 

. (1.3)

 

Here is the axis looked to the right in the global (rest) coordinate system, and is the similar axis in the local system of coordinates uniformly moving to the left. The latter system is usually associated with the body (wing) immersed in the uniform flow.

 

Equation (1.3) describes propagation of the perturbations of the pressure or the velocity potential and can be immediately obtained from the Euler equations linearised about the uniform flow, cf. [Sofronov-JMAA] .

 

In order to numerically simulate acoustic waves governed by (1.3) high-order finite volume or finite difference methods are used. These methods require non-reflecting boundary conditions on open boundaries such that they could have really small reflections*. At least the error arising because of spurious reflections should not be greater than the approximation error in the interior. One of the best choices is exact (transparent) boundary conditions. The correspondent operators have been obtained and implemented for the wave equation in case of , spherical boundary, see [Sofronov-DAN] , [Sofronov-EJAM] , [Grote-Keller SIAM] , [Hagstrom-AN] , and in case of , channel, see [BBS-AIAA] . In both cases a spectral approach is used to derive analytically the desired operators of the boundary conditions. Namely the Fourier method on the open boundaries: spherical functions or imaginary exponentials for 3D or 2D spherical boundary, and cosines for inflow/outflow cross-sections of the channel, respectively. The larger number of the basis functions being taken into account the higher accuracy of a discrete counterpart (i.e. smaller amount of reflections). Each Fourier coefficient which is a function of time and the normal spatial variable is treated separately by recurrence formulae with respect to time.

 

Possibility of use of the Fourier method is a key feature in the construction of the abovementioned operators of boundary conditions. Evidently this spectral approach permits to tune the accuracy of required discretization to the approximation error in the interior. In both spherical and plain cases the Fourier method was used owing to the fact that the governing equations have uniform coefficients on the transversal coordinate surfaces: spherical or polar coordinates for Eq. (1.2) outside 3D/2D sphere, and Cartesian coordinates for Eq. (1.3) in the channel to the left from the inflow cross-section and to the right from the outflow cross-section, see Figure 1, shaded regions extended to infinity (white regions correspond to computational domains).

 

Figure 1 : Spherical (left) and plain-channel (right) geometries

 

 

 

 

Unfortunately an immediate treatment of Eq. (1.3) with in the spherical geometry, see Figure 1 left, by the spectral approach similar to the case of what is very desirable for the aeroacoustics in open domains is not possible because of variable coefficients with respect to the azimuth angle: the Fourier method does not work. In this paper, we propose a way to find an approximate solution to this challenge. Idea consists of using a discrete counterpart of the problem from the beginning with successive derivation and efficient approximation of the spectral boundary operator in terms of a discrete Fourier basis.

 

The outline of the paper is as follows. In Section 1 we formulate the problem, the two-dimensional case is considered for simplicity (polar coordinates). Section 2 describes main steps of the algorithm of generating boundary conditions correspondent to the homogeneous case of zero velocity , a bridge to the inhomogeneous case is made. The latter is considered in Section 3. Numerical examples demonstrating accuracy of the approach are given in Section 4. Section 5 contains several conclusions.

 

Note that idea to use discrete governing equations outside domain of interest was proposed by V. S. Ryabenkii and it has been explored, in particular in [RTT-JCP] to construct non-local ABCs for 3D wave equation. A principal discrimination of this and our approaches consists of ways of approximations of obtained operators that originally require too large computational resource.

The method [RTT-JCP] is based on the property of 3D wave equation to have lacuna, while our approach develops approximation of boundary operator by sum-of-exponentials. The latter is more generic from the view of applications; at least we can treat 2D wave equation where the method [RTT-JCP] does not work.

 

1.    Problem formulation and governing equations

 

 

We omit the prime in Eq. (1.3) hereafter for a convenience, and restrict ourselves by the two-dimensional case; the approach is generalized straightforwardly for the three-dimensional case as well.

Let us consider Cauchy problem for the equation

(2.1)

supposing that exciting data functions are concentrated inside a finite domain:

 

The original problem consists of constructing artificial boundary conditions ABCs on such that waves propagate through without reflection.

A disk with a radiusis taken here as the domain:

 

.

 

Remark. As usual in formulation of problem of generating ABCs one needs to point exactly the governing equations outside only. Consequently the aim is to replace these governing equations by proposed ABCs. No concretization of equations inside is required as a rule.

We introduce the polar coordinates

 

(2.2)

and rewrite (2.1) in the form:

 

(2.3)

 

 

Equation (2.3) or more precisely some its difference counterpart outside the disk will be the main equation in our analysis. The desired low-reflecting boundary conditions will be generated numerically.

 

We will need also a simple local boundary condition for (2.3) on . To generate it let us take the well-known condition

 

for the 2D wave equation and make the change of variables (1.1) . After some algebra, putting in the time-dependent coefficients, we obtain the desired local condition:

 

(2.4)

 

 

2.    Case a=0

 

We reproduce here main elements of the approach [Sofronov-EJAM] of generating analytical transparent boundary conditions for the Eq. (2.3) , . This will be a background to make a generalization for the case .

Let us consider first the following auxiliary extended IBVP for a function on

(3.1)

Here is the basis of imaginary exponentials on , is the Diracs delta function.

 

The problem has analytical solution expressed in the form

(3.2)

where is the modified Bessel function (see for example [A-S] ), denotes the inverse Laplace transform .

 

Evidently, the solution of the IBVP with arbitrary Dirichlect boundary data ,

(3.3)

is written down as

(3.4)

where denotes convolution with respect to the time variable, and are the Fourier coefficients defined from the decomposition of

on the boundary .

 

Notice that the convolution kernel is written in the factorization from

(3.5)

with

Coming back to the interior IBVP we propose to use (and we do use) the formula (3.4) to calculate function on the open boundary while developing a numerical algorithm for solving the reduced problem in . Let us clarify this on the example of an explicit difference scheme. Denote by last three -grid points of the polar mesh in . Suppose the solution is already known for the time-layers with . Then using a second-order finite-difference scheme one can update the solution on the time layer for all points except the boundary point . The solution at point is calculated by (3.4) taking Dirichlect data at as . Figure 2 schematically represents the algorithm. Thus we obtain the transition operator from the layer to the .

 

It is important to emphasize that the convolution kernel is handled by the sum-of-exponentials approximations:

 

This representation allows the recursive evaluation of the convolution operator in (3.4) and dramatically reduces; see details in [Sofronov-EJAM] .

 

Figure 2: Schematic representation of the update algorithm.

3.    Case a>0

 

Now consider the equation (2.3) for . Similarly to (3.1) we have the following auxiliary IBVP

(4.1)

where denotes the wave operator (2.3) in moving media.

Auxiliary elementary kernels

 

Evidently there is no simple analytical formula for the solution in this case. Therefore let us consider the discrete counterpart for (4.1) :

(4.2)

 

I.e. we introduce the polar grid in

and suppose that we are able to calculate solution of (4.2) - grid function with . The details of the finite-difference scheme will be discussed below.

Evidently we have discrete problems (4.2) since the discrete basis on consists of discrete functions .

 

First, similarly to (3.3) we consider the discrete problem

(4.3)

with arbitrary Dirichlect data . Its solution can be expressed in terms of the solution :

(4.4)

where are the Fourier-coefficients of in the basis , i.e.

and denotes the discrete convolution operator defined by the following rule

Next we introduce the elementary kernels which are the Fourier-components of in the basis numerated so that

(4.5)

 

here can have the values .

 

The following matrix notation clarifies the formula (4.5)

(4.6)

 

Remark. In case owing to the separation of variables the matrix in (4.6) is diagonal, i.e. if , cf. (3.5) .

Each elementary kernel depends now on temporal index only (at fixed radial index ).

 

Thus (4.4) can be rewritten in the form

(4.7)

 

Formula (4.7) will serve us to generate low-reflecting boundary conditions, cf. (3.5) .

Numerical aspects of the algorithm

 

At first we say some words about the finite-difference scheme for (2.3). All derivatives are approximated by central second order finite differences. The scheme is implicit in time because of the mixed derivatives and at each time step we have to solve the linear system , where is a solution on the current time step . Matrix has the form , where is the identical matrix, corresponds to the derivatives, corresponds to the derivatives. To inverse the matrix we use the simple iterations in form

with , is approximation of , is an iterative parameter. On each iteration step we have to inverse two three-diagonal matrices that are handled by the sweep method.

 

We define the basis by imaginary exponentials on the equidistant grid:

.

Discrete delta function is given simply by

According to (4.7) we must calculate the kernels for all time steps such that , where is a calculation time. However, similarly to the update algorithm shown in Figure 2, it is enough to keep functions only for single value of . Nevertheless these calculations of elementary kernels in (4.6) are very expensive. It requires also large memory resources to keep as well as large computational costs to calculate the convolution in (4.7) .

 

That is why we have developed set of modifications to (4.7) in order to sharply reduce the computational costs. First we subdivide the passing waves onto low and high frequencies (with respect to spatial grid size). Therefore we decrease the summation limits in (4.7) . Only low-frequency harmonics with are treated accurately with the non-local discrete boundary condition. For high-frequency harmonics, discretization of the local boundary condition (2.4) is used. The new limits correspond to the truncation of the matrix , i.e. instead of matrix we consider matrix.

 

Next we introduce restriction on the summation index : let belong to the interval simply throwing away any others .

 

We will see in the examples of numerical simulation that such approximations of the full matrix in (4.6) do have a sense: it is sufficient to take small enough to produce accurate results.

 

Thus we really need only a band submatrix in (4.6) :

 

Finally, and this is the most valuable modification to reduce computational costs, we use a technique developed in [AES-CMS] and approximate each discrete convolution kernel by sum of exponentials:

(4.8)

here is the power in the last term.

 

This representation allows for the recursive evaluation of the convolutions in (4.7) .

In practice we use for large enough computational time and therefore we need calculate in advance and keep only about complex numbers to represent each elementary kernel. So the cost of our approximation to (4.7) is not too large: more exactly the requirements on memory are estimated by of real values and the computational cost is estimated by operations per time step, .

 

Incorporation of the modified formula (4.7) into a difference scheme for interior problem in order to update solution at the external open boundary is made in the same manner as described in the previous section. The only discrimination is that we must treat a band matrix of elementary kernels (width) instead of simply diagonal one (the parameter)

 

According to the algorithm described in [AES-CMS] the approximation (4.8) can be obtained by knowledge of at . Thus the extended auxiliary problems are computed only for several first time steps.

4.    Numerical examples

 

In order to avoid singularities in the origin we consider the annular domain . We impose homogeneous Dirichlect boundary conditions at and our discrete non-local boundary conditions at . The velocity and .

Two equidistant meshes are used: coarse one with , and fine one with .

 

In the simulations we consider the equation (2.3). The initial data is taken to zero and the source is introduced as a right-hand side in equation (2.3) having the form

Here is so-called Ricker signal with the central frequency , see Figure 3 (left)

the source distribution is on Figure 3 (central) with

and the frequency dependence of the source is on Figure 3 (right)

 

Figure 3: time dependency, Ricker function (left); distribution on r variable (central); -distribution (right).

 

We compare calculated solutions with the reference solution obtained on the extended area and on the very fine mesh so that this discrete solution can be identified with the exact.

 

Below we represent the results in continuous norm measured over our annular domain . Note that the errors for -norm have the same orders and behavior.

 

In Figure 4 (top) we represent the relative errors of the solutions obtained on the extended areas, i.e. the errors that are due to the approximations of the difference scheme on our grids.

 

Then in Figure 4 (bottom) we represent the relative error of the solutions with low-reflecting boundary conditions in form (4.7) compared with the solution computed on the extended domain on the same mesh. We set for the coarse mesh and for the fine one. The results are pretty well: boundary errors are much less than the approximation errors (20 to 30 times) and dont affect the resulting error.

 

Note that if we use the local boundary condition (2.4) at then the errors have the values compared with the solution, i.e. the errors are about 100%.

 

The demonstrated results confirm that can be small enough compared to . In Figure 5 - norm of in logarithm scale is shown. We take here -norm is calculated with respect to , correspondent to 5 seconds. One can observe a sharp peak near .

 

 

Figure 4: relative errors of the solution calculated on extended domain, dashed line is for the coarse mesh (), solid line is for the fine mesh (). Top figure corresponds to the reference discrete solution on the extended region:

, bottom to the solution with our boundary condition at :

We summarize the influence of the parameters on relative errors in the tables below. Table 1 results correspond to the coarse mesh, Table 2 to the fine one.

The reader can compare these values with those on the Figure 4, right: coarse and fine grids, respectively.

 

 

5.5E-02

2.7E-02

2.2E-03

1.9E-03

5.5E-02

2.8E-02

1.2E-03

3.5E-04

Table 1: Relative errors for the coarse mesh for the different sizes of the matrix band.

 

 

5.6E-02

2.9E-02

2.4E-03

2.4E-03

5.6E-02

3.0E-02

7.6E-04

6.3E-04

Table 2: Relative errors for the fine mesh for the different sizes of the matrix band.

 

It is impotent to notice that we use the approximation representation (4.8) to reconstruct the kernels for where is large enough. According to the algorithm for finding coefficients in (4.8) we need function on a short time interval only, i.e. , where . Of course such construction is not correct for arbitrary medium. For example it is obvious that we cannot apply such procedure for the medium with some inhomogeneous in some distance from the external boundary. But our medium has no obstacles and we dont expect some impulse arrived from outside.

 

Another difficulty with usage of (4.8) occurs while considering large values of . If is small enough the kernel looks like one presented in Figure 6 (top.). Such kernels are approximated very well by the sum-of-exponentials (4.8) . But if increases the kernel becomes like one from Figure 6 (bottom). It is impossible here to construct the approximation (4.8) with decaying exponentials at short time. Fortunately for the case of our coarse and fine meshes we dont need to deal with such bad kernels. Pretty well accuracy is achieved without considering kernels of these types.

 

If we need finer meshes we must consider bad kernels as well. Let us discuss two possible ways how to avoid the difficulties with the approximation.

Evidently the nature of this oscillation behavior is owing to the delta-function Dirichlect boundary data while calculating the elementary kernels, see (4.2) . Therefore the first way is to work with submeshes. I.e. we can try to find the kernels on finer sub meshes with smooth delta functionoriginated from the main grid. Thus the kernels will be smoother and could permit desired approximations. The second way consists in using more sophisticated finite-difference scheme in (4.2) that gives smaller oscillations for discontinuous initial data.


Figure 5: - norm of , , versus distance . Velocity a=0.2 (top) and a=0.7 (bottom).

Figure 6: Amplitude of elementary kernels for ; (top) and (bottom).

 

5.    Conclusions

 

In this paper we have introduced the novel approach of constructing discrete transparent boundary condition for the wave equation in the moving media. Necessary approximation modifications of exact formula leading to low-reflecting boundary conditions are proposed. These modifications permit to rapidly calculate the boundary operator. Numerical examples show that the error due to reflections is much less than the error due to finite-difference scheme.

 

Also the described algorithm may be considered as a generic method to construct low-reflecting boundary conditions for the different kind of equations and boundary shapes. We already have some results concerning the wave equation in the layered media and we think about another applications.

 

As mentioned above there are some open questions while approximating the kernels. They required more detailed investigation and this will be a part of our future work.

 

References

 

[Sofronov-JMAA] Sofronov, I. L. Non-reflecting inflow and outflow in wind tunnel for transonic time-accurate simulation, J. Math. Anal. Appl., V. 221, (1998) 92115.

 

[Sofronov-DAN] Sofronov, I. L. Conditions for complete transparency on a sphere for the three-dimensional wave equation, Russ. Acad. Sci. Dokl. Math. Vol. 46, No.2 (1993) 397401.

 

[Sofronov-EJAM] Sofronov, I. L. Artificial boundary conditions of absolute transparency for two- and three-dimensional external time-dependent scattering problems, Euro. J. Appl. Math., V.9, No.6 (1998) 561588.

 

[Grote-Keller SIAM] M.J.Grote and J.B.Keller, Exact nonreflecting boundary conditions for the time dependent wave equation, SIAM J.Appl.Math. 55 (1995), 280-297.

 

[Hagstrom-AN] T.Hagstrom, Radiation boundary conditions for the numerical simulation of waves,

Acta Numerica 8 (1999), 47-106, Cambridge: Cambridge University Press, 47-106.

 

[BBS-AIAA] Ballmann J.; Britten G.; Sofronov I. Time-accurate inlet and outlet conditions for unsteady transonic channel flow, AIAA Journal, Vol. 40 (2002), No. 2., 17451754.

 

[RTT-JCP] V. S. RYABEN'KII, S. V. TSYNKOV, AND V. I. TURCHANINOV, Global Discrete Artificial Boundary Conditions for Time-Dependent Wave Propagation, J. Comput. Phys., 174 (2001) pp. 712-758.

 

[A-S] 1. Abramovitz M., Stegun I. A. Handbook of Mathematical Functions, National Bureau of Standards, Applied Math. Series #55. Dover Publications, 1965.

 

[AES-CMS] Arnold A; Ehrhardt M.; Sofronov I. Discrete transparent boundary conditions for the Schroedinger equation: Fast calculation, approximation, and stability, Comm. Math. Sci. 1 (2003), 501-556.

 

 

 



* The term non-reflecting boundary conditions is an ideal used often in the literature for majority of proposed boundary conditions that do have reflections, in fact. In this sense the term low-reflecting boundary conditions clarified in the next sentence seems to be more relevant.