Spectral nonlocal boundary conditions for the wave equation in moving media


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 twodimensional 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. Ryaben’kii and it has been explored, in particular in [RTTJCP] to construct nonlocal 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 [RTTJCP] is based on the property of 3D wave equation
to have lacuna, while our approach develops approximation of boundary operator
by sumofexponentials. The latter is more generic from the view of applications;
at least we can treat 2D wave equation where the method [RTTJCP] does not work.
We
omit the prime in Eq. (1.3)
hereafter for a convenience, and restrict ourselves
by the twodimensional case; the approach is generalized straightforwardly for
the threedimensional 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 radius_{}is 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
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
“lowreflecting” 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 wellknown condition
_{}
for
the 2D wave equation and make the change of variables (1.1)
. After some algebra, putting _{} in the timedependent
coefficients, we obtain the desired local condition:
We
reproduce here main elements of the approach [SofronovEJAM] 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 Dirac’s delta
function.
The problem has analytical
solution expressed in the form
_{}
(3.2)
where _{} is the modified
Bessel function (see for example [AS]
),_{} 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
timelayers with _{}. Then using a secondorder finitedifference 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
sumofexponentials approximations:
_{}
This
representation allows the recursive evaluation of the convolution
operator in (3.4)
and dramatically reduces; see details in [SofronovEJAM] .
Figure 2:
Schematic representation of the update algorithm.
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.
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 finitedifference 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
Fouriercoefficients 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
Fouriercomponents 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 lowreflecting boundary
conditions, cf. (3.5)
.
At
first we say some words about the finitedifference
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 threediagonal
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 lowfrequency harmonics with _{} are treated
accurately with the nonlocal discrete boundary condition. For highfrequency
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 [AESCMS] 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 [AESCMS] the approximation (4.8)
can be obtained by knowledge of _{} at _{}. Thus the extended auxiliary problems are computed only for
several first time steps.
In
order to avoid singularities in the origin we consider the annular domain _{}. We impose homogeneous Dirichlect boundary conditions at _{} and our discrete
nonlocal 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 righthand side in equation (2.3)
having the form
_{}
Here _{} is socalled 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 lowreflecting
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 don’t 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.5E02 
2.7E02 
2.2E03 
1.9E03 
_{} 
5.5E02 
2.8E02 
1.2E03 
3.5E04 
Table 1: Relative errors for the coarse
mesh for the different sizes of the matrix _{} band.

_{} 
_{} 
_{} 
_{} 
_{} 
5.6E02 
2.9E02 
2.4E03 
2.4E03 
_{} 
5.6E02 
3.0E02 
7.6E04 
6.3E04 
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 don’t 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 sumofexponentials (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 don’t 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 deltafunction
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”
function_{}originated from the main grid. Thus the kernels will be
smoother and could permit desired approximations. The second way consists in
using more sophisticated finitedifference 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).
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
lowreflecting 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 finitedifference scheme.
Also
the described algorithm may be considered as a generic method to construct
lowreflecting 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.
[SofronovJMAA] Sofronov, I. L. Nonreflecting
inflow and outflow in wind tunnel for transonic timeaccurate simulation,
J. Math. Anal. Appl., V. 221, (1998) 92—115.
[SofronovDAN] Sofronov,
I. L. Conditions for complete transparency on a sphere for the
threedimensional wave equation, Russ. Acad. Sci. Dokl. Math. Vol. 46, No.2
(1993) 397—401.
[SofronovEJAM] Sofronov, I. L. Artificial boundary conditions of absolute
transparency for two and threedimensional external timedependent scattering
problems, Euro. J. Appl. Math., V.9, No.6 (1998) 561—588.
[GroteKeller
SIAM] M.J.Grote and J.B.Keller, Exact nonreflecting
boundary conditions for the time dependent wave equation, SIAM J.Appl.Math.
55 (1995), 280297.
[HagstromAN]
T.Hagstrom, Radiation boundary
conditions for the numerical simulation of waves,
Acta Numerica 8 (1999), 47106, Cambridge: Cambridge University Press, 47106.
[BBSAIAA] Ballmann J.; Britten G.; Sofronov I. Timeaccurate
inlet and outlet conditions for unsteady transonic channel flow, AIAA
Journal, Vol. 40 (2002), No. 2., 1745—1754.
[RTTJCP] V. S. RYABEN'KII, S. V. TSYNKOV, AND V. I. TURCHANINOV,
Global Discrete Artificial Boundary Conditions for TimeDependent Wave
Propagation, J. Comput. Phys., 174 (2001) pp. 712758.
[AS] 1. Abramovitz M., Stegun I. A. Handbook of
Mathematical Functions, National Bureau of Standards, Applied Math. Series
#55. Dover Publications, 1965.
[AESCMS] Arnold A; Ehrhardt
M.; Sofronov I. Discrete transparent boundary conditions for the
Schroedinger equation: Fast calculation, approximation, and stability, Comm. Math. Sci. 1 (2003), 501556.
* The term “nonreflecting 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 “lowreflecting boundary conditions” clarified in the next sentence seems to be more relevant.