DTIC ADA453168: Artificial Boundary Conditions for the Numerical Simulation of Unsteady Acoustic Waves

34 Pages·0.47 MB·English

by  Defense Technical Information Center

Artificial Boundary Conditions for the
Numerical Simulation of Unsteady Acoustic
Waves
S. V. Tsynkov1
Department of Mathematics and Center for Research in Scientific Computation
(CRSC), North Carolina State University, Box 8205, Raleigh, NC 27695, USA.
Abstract
Weconstructnon-localartificialboundaryconditions(ABCs)forthenumericalsim-
ulation of genuinely time-dependent acoustic waves that propagate from a compact
source in an unbounded unobstructed space. The key property used for obtaining
the ABCs is the presence of lacunae, i.e., sharp aft fronts of the waves, in wave-type
solutions in odd-dimension spaces. This property can be considered a manifestation
of the Huygens’ principle. The ABCs are obtained directly for the discrete formula-
tion of the problem. They truncate the original unbounded domain and guarantee
the complete transparency of the new outer boundary for all the outgoing waves. A
centralfeatureoftheproposedABCsisthattheextentoftheirtemporalnonlocality
is fixed and limited, and it does not come at the expense of simplifying the original
model. It is rather a natural consequence of the existence of lacunae, which is a fun-
damental property of the corresponding solutions. The proposed ABCs can be built
for any consistent and stable finite-difference scheme. Their accuracy can always
be made as high as that of the interior approximation, and it will not deteriorate
even when integrating over long time intervals. Besides, the ABCs are most flexible
from the standpoint of geometry and can handle irregular boundaries on regular
grids with no fitting/adaptation needed and no accuracy loss induced. Finally, they
allow for a wide range of model settings. In particular, not only one can analyze the
simplest advective acoustics case with the uniform background flow, but also the
case when the waves’ source (or scatterer) is engaged in an accelerated motion.
Key words: Time-dependent sound waves, unbounded domains, lacunae,
non-deteriorating method, long-term computation, limited temporal nonlocality.
Email address: [email protected] (S. V. Tsynkov).
URL: http://www.math.ncsu.edu/∼stsynkov (S. V. Tsynkov).
1 TheauthorgratefullyacknowledgessupportbyAFOSR,GrantF49620-01-1-0187,
and that by NSF, Grant DMS-0107146.
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it

does not display a currently valid OMB control number.
1. REPORT DATE 3. DATES COVERED
2003 2. REPORT TYPE 00-00-2003 to 00-00-2003
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Artificial Boundary Conditions for the Numerical Simulation of Unsteady
5b. GRANT NUMBER
Acoustic Waves
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
North Carolina State University,Center for Research in Scientific REPORT NUMBER
Computation,Raleigh,NC,27695-8205
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
see report
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 33
unclassified unclassified unclassified
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
1 Introduction
Artificial Boundary Conditions (ABCs) is a common name for a group of
methods employed for solving infinite-domain problems on a computer. ABCs
facilitate truncation of the original unbounded domain, and provide the re-
quired closure for the resulting truncated formulation. The literature on the
subject of ABCs is broad, and we refer the reader to the review papers [1–3].
Two opposite trends can be identified in building the ABCs. High accuracy
requirements typically imply that the ABCs should be nonlocal. In partic-
ular, the exact ABCs are always nonlocal in multidimensional settings. For
unsteady problems, this also means nonlocality in time. Clearly, nonlocality
of the ABCs may result in high computational costs and elaborate implemen-
tation strategies. As such, local alternatives obtained either independently or
as a approximations of nonlocal ABCs become viable. These ABCs are usually
inexpensive and easy to implement, but may lack computational accuracy.
In the current paper we propose the ABCs for the numerical simulation of
non-dispersive unsteady acoustic waves. We assume that there is a (complex)
phenomenon/process confined to a bounded region that manifests itself by the
radiation of acoustic waves in the far field. The waves subsequently propagate
across the unbounded and unobstructed space, which is assumed isotropic
and homogeneous. The proposed ABCs are nonlocal and guarantee that the
accuracy of the boundary treatment can always be made at least as high as
that of the interior discretization. However, the key feature of the proposed
ABCs is that the extent of their temporal nonlocality is limited, and does not
increase as the time elapses. The bound on temporal nonlocality comes as a
consequence of the presence of lacunae, i.e., sharp aft fronts of the waves, in
the solutions of the linearized Euler equations. In general, existence of the
lacunae is a fundamental property of wave-type solutions in odd-dimension
spaces. It is often referred to in connection with the Huygens’ principle.
The ABCs are constructed by decomposing the original problem into the inte-
riorandauxiliarysub-problems.Thelatterislinearandhomogeneousthrough-
out the entire space, and is integrated with the help of lacunae. In so doing,
the interior solution is used to generate sources for the auxiliary problem, and
the auxiliary solution is used to provide the missing boundary data, i.e., the
require closure, for the interior problem.
Several other nonlocal ABCs’ methodologies for unsteady waves have been
recently advocated in the literature, most notably [4–10]; we also mention the
survey paper [11] and the bibliography there. Compared to these techniques,
ourapproachhasanumberofdistinctivecharacteristics.Itshighaccuracyand
restrictedtemporalnonlocalityhavealreadybeenmentioned.Besides,thepro-
2
posed lacunae-based ABCs are obtained directly for the discrete formulation
of the problem, and can supplement any consistent and stable finite-difference
scheme. In other words, they bypass the common stages of first deriving and
then approximating the continuous ABCs, see [2]. The lacunae-based ABCs
are especially designed to withstand long-term numerical integration with no
deterioration of accuracy. Moreover, neither they are restricted to any par-
ticular shape of artificial boundary nor require grid fitting. They also enable
analysis of a variety of problem formulations, e.g., the case of sound sources
moving with acceleration, which is equivalent to advective acoustics with un-
steady background flow. We note that the source motion does not have to be
considered relative to an ambient medium only. Via the appropriate Galileo
transform, it can be considered relative to a given mean flow as well.
The ABCs that we derive in the current paper for acoustics generalize and
further extend our previous work [12,13], in which a similar methodology was
introduced for the scalar wave equation. The current approach also fits into
the general theoretical framework of [14].
2 Lacunae of the Wave Equation
Consider a three-dimensional wave equation, x = (x ,x ,x ):
1 2 3
1 ∂2ϕ ∂2ϕ ∂2ϕ ∂2ϕ
−(cid:32) + + (cid:33) = f(x,t), t ≥ 0, (1a)
c2 ∂t2 ∂x2 ∂x2 ∂x2
1 2 3
∂ϕ
ϕ(cid:12)(cid:12)(cid:12)(cid:12)t=0 = ∂t(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)t=0 = 0, (1b)
with homogeneous initial conditions. For every (x,t), the solution ϕ = ϕ(x,t)
of problem (1a), (1b) is given by the Kirchhoff integral (see, e.g., [15]):
1 f (ξ,t−(cid:37)/c)
ϕ(x,t) = dξ, (2)
(cid:90)(cid:90)(cid:90)
4π (cid:37)
(cid:37)≤ct
where ξ = (ξ ,ξ ,ξ ), (cid:37) = |x − ξ| = (x −ξ )2 +(x −ξ )2 +(x −ξ )2,
1 2 3 (cid:113) 1 1 2 2 3 3
and dξ = dξ dξ dξ . If we assume that the right-hand side (RHS) f(x,t)
1 2 3
of equation (1a) is compactly supported in space-time on the domain Q ⊂
R3 ×[0,+∞), then formula (2) immediately implies that
ϕ(x,t) ≡ 0 for (x,t) ∈ (x,t) |x −ξ| < c(t−τ), t > τ . (3)
(cid:92) (cid:110) (cid:12) (cid:111)
(cid:12)
(ξ,τ)∈Q (cid:12)
The region of space-time defined by formula (3) is called lacuna of the solution
ϕ = ϕ(x,t). This region is obtained as the intersection of all characteristic
3
cones of equation (1a) once the vertex of the cone sweeps the support of the
RHS: suppf ⊆ Q. From the standpoint of physics, lacuna corresponds to
that part of space-time, on which the waves generated by the sources f(x,t),
suppf ⊆ Q, have already passed, and the solution has become zero again.
The phenomenon of lacunae is inherently three-dimensional. The surface of
the lacuna represents the trajectory of aft (trailing) fronts of the waves. The
existence of aft fronts in odd-dimension spaces is known as the Huygens’ prin-
ciple, as opposed to the so-called wave diffusion which takes place in even-
dimension spaces, see, e.g., [15]. Let us also note that the aft fronts and the
lacunae would still be present in the solution of the wave equation (1a) if the
homogeneous initial conditions (1b) were replaced by some inhomogeneous
initial conditions with compact support.
The notion of lacunae (or lacunas) was first introduced and studied by
Petrowsky in [16] (see also an account in [17, Chapter VI]) for a variety of
hyperbolic equations and systems; and general characterization of their coef-
ficients was provided that would guarantee existence of the lacunae. However,
since work [16] no constructive examples of either equations or systems have
been obtained for which lacunae would be present in the solutions, besides the
actual wave equation (1a), as well as those equations that either reduce to or
are derived from, the wave equation.
3 The Acoustics System of Equations
The acoustics (linearized Euler’s) system in its simplest form governs the
propagation of sound in an ambient compressible fluid [18, Chapter VIII]:
1 ∂p
+ρ ∇·u = ρ q ,
c2 ∂t 0 0 vol
(4)
∂u
ρ +∇p = b .
0 vol
∂t
In (4), c is the speed of sound, and the variables u = u(x,t) and p = p(x,t)
represent small perturbations of the velocity and pressure, respectively. Adi-
abatic law p = c2ρ has been applied that relates p to the density perturba-
tion ρ. The quantity ρ is the constant background density. The source term
0
q = q (x,t) that alters the balance of mass in the system is known as vol-
vol vol
ume velocity per unit volume, and the source term b = b (x,t) that alters
vol vol
the balance of momentum is known as force per unit volume (see, e.g., [19]).
Proposition 1 Assume that the initial velocity field is conservative:
u(x,0) = ∇ϕ (x). Then, the velocity potential exists in the solutions of
0
system (4) for all t ≥ 0 if and only if the forcing b is also conservative:
vol
b (x,t) = ∇ψ(x,t). In this case, the potential satisfies the wave equation.
vol
4
PROOF. Let us assume that the force field is conservative, and integrate the
second equation of system (4) with respect to time starting from t = 0:
t t
ρ u(x,t) = ρ u(x,0)− ∇p(x,τ)dτ + b (x,τ)dτ
0 0 (cid:90) (cid:90) vol
0 0
t t
= ρ ∇ϕ (x)− ∇p(x,τ)dτ + ∇ψ(x,τ)dτ
0 0 (cid:90) (cid:90)
0 0 (5)
t t
= ∇(cid:20)ρ0ϕ0(x)−(cid:90) p(x,τ)dτ +(cid:90) ψ(x,τ)dτ(cid:21)
0 0
d=ef ρ ∇ϕ(x,t).
0
Therefore, we obtain:
∂ϕ(x,t)
u(x,t) = ∇ϕ(x,t), ρ = −p(x,t)+ψ(x,t), (6)
0
∂t
and by differentiating the second relation in (6) w.r.t. time and subsequently
substituting it into the first equation of system (4), we conclude that the
potential ϕ = ϕ(x,t) defined by (5) will satisfy the wave equation (1a) with
the RHS given by f(x,t) = −q (x,t)+ 1 ∂ψ(x,t).
vol ρ0c2 ∂t
Conversely, let us assume that the velocity field is conservative: u(x,t) =
∇ϕ(x,t).Then,themomentumequationinsystem(4)canberecastasfollows:
∂ϕ
∇(cid:20)ρ0 +p(cid:21) = bvol,
∂t
which means that the forcing must be conservative: b (x,t) = ∇ψ(x,t). (cid:50)
vol
Proposition1alsoimpliesthatifb (x,t) ≡ 0thenthepotentialalwaysexists.
vol
In this case, the second equality of (6) reduces to the conventional relation
between the potential and pressure: ρ ∂ϕ = −p, see [18, Chapter VIII].
0∂t
Applying the gradient operator to the wave equation for the potential, we
conclude that the velocity vector also satisfies the wave equation:
1 ∂2u 1 ∂b
−∆u = −∇q + vol. (7)
c2 ∂t2 vol ρ c2 ∂t
0
Next, by differentiating the first equation of (4) with respect to time, taking
the divergence of the second equation of (4), and substituting the result into
the first one, we arrive at the wave equation for p:
1 ∂2p ∂q
−∆p = −∇·b +ρ vol. (8)
c2 ∂t2 vol 0 ∂t
5
The key consideration of interest is that if the RHSs of system (4) are com-
pactly supported in space and time on some domain Q ⊂ R3 ×[0,+∞), then
the RHSs of both equation (7) and equation (8) will also be compactly sup-
ported on the same domain Q. Consequently, solutions of equations (7) and
(8) will have lacunae of the same shape as prescribed by formula (3). Thus,
we have arrived at the following
Proposition 2 Let the RHSs of system (4) be compactly supported in space
and time: suppq (x,t) ⊆ Q and suppb (x,t) ⊆ Q, where Q ⊂ R3 ×
vol vol
[0,+∞). Let also u(x,t) = ∇ϕ(x,t). Then, solutions of the acoustics sys-
tem (4) with homogeneous initial conditions will have lacunae of the same
geometry as provided by formula (3).
In other words, we have shown that existence of the velocity potential is suf-
ficient for the acoustics system (4) to have lacunae in its solutions; and that
the geometry of the lacunae is determined by that of the support of the RHSs,
as in the case of the wave equation. We do not know whether this condition
is also necessary for having the lacunae. However, in the view of the comment
provided in the end of Section 2, we will hereafter be using the foregoing suf-
ficient condition as the only reliable indication of the presence of lacunae in
the solutions of system (4).
The group of numerical algorithms for acoustics that we describe hereafter
will be essentially based on the presence of lacunae. However, the Kirchhoff
formula (2) will never be used as an actual part of the algorithm construction,
it will only be needed at the theoretical stage, for determining the shape of
the lacunae that will later be incorporated into a purely finite-difference con-
text. Previously, the idea of using the Huygens’ principle for constructing the
ABCs was promoted by Ting & Miksis [20] and Givoli & Cohen [21]. Both
papers, however, have suggested to use numerical quadratures to approximate
the integral (2), and then couple it with the interior solution. Moreover, the
approach of [20] has never been actually implemented in a practical compu-
tational setting, whereas the approach of [21] required artificial dissipation to
be added to the scheme to fix the arising instabilities.
4 Lacunae-Based Integration of the Acoustics System
We will be looking for an irrotational solution of system (4) on a bounded
domain S = S(t) ⊂ R3. We assume that this is a domain of fixed shape
that has a finite diameter d for all t ≥ 0. We also assume that ∀t ≥ 0 :
{suppq (x,t)∩R3} ⊆ S(t) & {suppb (x,t)∩R3} ⊆ S(t). In other words,
vol vol
we want to compute the acoustic field on a given domain that also contains
all the field sources. While always remaining finite in size, this domain S(t) is
6
allowed to move across the space R3 according to a general law:
t
u = u (t), x = x (t) = u (τ)dτ, (9)
0 0 0 0 (cid:90) 0
0
where u and x are the velocity vector and coordinates of a given point in-
0 0
side S(t), respectively. The only limitation that we impose is that motion (9)
be subsonic: max |u (t)| = k < c. Let now q˜ and b˜ be some station-
t 0 vol vol
ary acoustic sources (RHSs to system (4)), i.e., compactly supported func-
tions of x on the fixed domain S(0) for all t ≥ 0: {suppq˜ (x,t) ∩ R3} ⊆
vol
S(0) & {suppb˜ (x,t)∩R3} ⊆ S(0). Then, we can incorporate the transla-
vol
tional motion (9) into the acoustics system (4) as follows, see [12]:
1 ∂p(x,t)
+ρ ∇·u(x,t) = ρ q˜ (x −x (t),t) ≡ ρ q (x,t),
c2 ∂t 0 0 vol 0 0 vol
(10)
∂u(x,t)
ρ +∇p(x,t) = b˜ (x −x (t),t) ≡ b (x,t).
0 vol 0 vol
∂t
With no loss of generality, the initial conditions for integrating system (10)
will always be assumed homogeneous. It is easy to see from (10) that the time-
dependent nature of the acoustic field is caused by the motion of the sound
sources, on which the genuine unsteadiness of the sound generating mecha-
nisms may be “superimposed.” The domain of integration S(t) has a fixed
size/shape and traces the motion of the sources. As shown in [12, Appendix]
using the Galileo transform, system (10) is equivalent to advective acoustics
with stationary sources and unsteady background flow −u = −u (t).
0 0
Let us emphasize that the setup we have introduced is quite general and also
includes the case of the waves’ sources that move relative to a given mean flow
(rather than relative to the ambient fluid only). Indeed, a similar argument
based on the Galileo transform would imply that having the mean flow v and
0
the source motion u on top of it is equivalent to either stationary sources in
0
the mean flow v − u or alternatively, sources that move with the velocity
0 0
−v +u through the ambient medium. The latter setup is obviously the same
0 0
as (9), as long as the condition max |−v (t)+u (t)| = k < c is met.
t 0 0
The case of primary interest for our analysis will be that of continuously
operating sources in (10), t ∈ [0,+∞). Let us, however, assume for the time
being that not only the RHSs in system (10) are compactly supported in
space, but also that their “lifespan” in time is finite: suppq (x,t) ∈ Q and
vol
suppb (x,t) ∈ Q, where R3 ×[0,+∞) (cid:51) Q = {(x,t) x ∈ S(t),t ≤ t ≤ t }.
vol (cid:12) 0 1
(cid:12)
Then, Proposition 2 implies that no later than (cid:12)
d+(t −t )(c+k)
t = t ≡ t + 1 0 ≡ t +T (11)
2 0 c−k 0 int
all of the domain S(t) will fall into the lacuna defined by formula (3) (see
7
[12,13] for a more detailed argument), and will remain inside the lacuna con-
tinuously thereafter, i.e., for all t ≥ t . In other words, for the solution of
2
system (10), we will have p(x,t) = 0 and u(x,t) = 0 when x ∈ S(t) and
t ≥ t , where t is defined by formula (11).
2 2
Next, we realize that during the time interval T no wave can travel in
int
space further away than the distance cT from the boundary of the domain
int
S(t ). This means that we will also have p(x,t) = 0 and u(x,t) = 0 for
0
dist(x,S(t )) > cT and t ≤ t ≤ t . As such, instead of the free unob-
0 int 0 2
structed space outside S(t) we may consider outer boundaries with arbitrary
(reflecting) properties. As long as none of these boundaries is located closer
than cT to S(t ), the solution of (10) inside S(t) is not going to feel their
int 0
presence for t ≤ t ≤ t .
0 2
In fact, the foregoing limitation for the location of outer boundaries can even
be relaxed if instead of requiring that no wave may reach an outer boundary
before t = t we introduce a weaker requirement that no reflected wave may
2
reach S(t) before t = t . The latter consideration easily translates into the
2
following estimate for the minimal distance between the domain S(t ) and the
0
allowed location of any reflecting boundary (see [12,13] for more detail):
c+k
Z = T . (12)
min int
2
We note that both estimates (11) and (12) are conservative in the sense that
they do not take into account the direction of the motion. In case the motion
of the sources is characterized by a specific or predominant direction, then the
quantity Z can be further reduced in the orthogonal directions.
min
Altogether we conclude that the solution of system (10) driven by the sources
compactly supported on the domain Q = {(x,t) x ∈ S(t),t ≤ t ≤ t } and
(cid:12) 0 1
(cid:12)
subject to the homogeneous initial conditions, can(cid:12) be obtained on S(t) for all
t ≥ 0 as follows. First, this solution is obviously equal to zero for 0 ≤ t < t .
0
Next, on the time interval t ≤ t ≤ t , see formula (11), system (10) should
0 2
be integrated on the auxiliary domain of size
Z = d+2Z = d+(c+k)T (13)
min int
centered around S(t ), which in any event is going to yield the correct solution
0
inside S(t). Finally, for all t ≥ t the solution on S(t) will be equal to zero
2
again because all the waves will have left the domain by t = t (lacuna).
2
Let us now address the case of continuously operating sources. For that pur-
8
pose, we introduce a partition of unity on the semi-infinite interval t ≥ 0:
∞
∀t ≥ 0 : Θ(t−σTj) = 1, (14)
(cid:88)
j=0
where T > 0 and 1 ≤ σ < 1 are two parameters, and Θ = Θ(t) is a smooth,
2
even, “hat”-type function with suppΘ(t) ⊆ [−T, T]:
2 2
0, t > T,
 2
Θ(t) = 11,−Θ(σT −t), 0(σ≤−t1≤)T(σ<−t ≤12)TT,, (15)
2 2
Θ(−t), t < 0.
Then, we define compactly supported sources for j = 0,1,2,...:
q(j)(x,t) = q (x,t)Θ(t−σTj), suppq(j)(x,t) ⊆ Q ,
vol vol vol j
(16)
b(j)(x,t) = b (x,t)Θ(t−σTj), suppb(j)(x,t) ⊆ Q ,
vol vol vol j
where according to (15):
1 1
Q = (x,t) x ∈ S(t), σj − T ≤ t ≤ σj + T , (17)
j (cid:110) (cid:12) (cid:16) 2(cid:17) (cid:16) 2(cid:17) (cid:111)
(cid:12)
(cid:12)
and consider a collection of sub-problems driven by the RHSs (16):
1 ∂p(j) ∂u(j)
+ρ ∇·u(j) = ρ q(j), ρ +∇p(j) = b(j),
c2 ∂t 0 0 vol 0 ∂t vol
(18)
p(j)(x,t) = 0, u(j)(x,t) = 0.
(cid:12)(cid:12)t=(σj−1)T (cid:12)(cid:12)t=(σj−1)T
(cid:12) 2 (cid:12) 2
As formula (14) is a partition of unity, we have:
∞ ∞
q (x,t) = q(j)(x,t) and b (x,t) = b(j)(x,t),
vol (cid:88) vol vol (cid:88) vol
j=0 j=0
and because of the linear superposition, the solution p(x,t), u(x,t) of system
(10) subject to the homogeneous initial conditions at t = 0 can be expanded
in terms of the solutions to systems (18):
∞ ∞
p(x,t) = p(j)(x,t) and u(x,t) = u(j)(x,t). (19)
(cid:88) (cid:88)
j=0 j=0
The series (19) are formally infinite. However, for any t > 0 and x ∈ S(t) each
will, in fact, contain only a finite fixed number of nonzero terms. Indeed, due
to the causality for a given t > 0 and all (σj−1)T > t, i.e., all j > (t +1)/σ,
2 T 2
wewillhavep(j)(x,t) = 0andu(j)(x,t) = 0 ontheentirespaceR3.Moreover,
9