Computational Mathematics and Modeling, Vol. 18, No. 3, 2007

ARTIFICIAL BOUNDARY CONDITIONS FOR TWO-DIMENSIONAL EQUATIONS OF FLUID DYNAMICS. 1. CONVECTIVE WAVE EQUATION L. V. Dorodnitsyn

UDC 519.6:533.1,534.2

The article models external flow problems on artificially bounded regions. In the linear approximation we examine the reflection of acoustic waves in a moving medium, incident at various angles on a fixed boundary. We consider the construction of various boundary conditions and estimate their reflecting properties for plane waves and waves from point sources. The plane wave approximation is justified theoretically. A method is proposed for estimating the integral reflection coefficient for waves with a whole range of incidence angles.

Introduction Simulation of external gas flows, in particular, flow past bodies, generally requires mapping the problem from an infinite to a finite region. The finite region has an artificial boundary, on which boundary conditions (also called artificial) have to be defined. These artificial boundary conditions should ensure that the solutions of the problems in the infinite and the finite regions match or are at least close to one another. In physical terms, the artificial boundary should not reflect waves from inside the region to the outside. This gives rise to the concept of “nonreflecting boundary conditions”. In problems of subsonic fluid dynamics the formulation of artificial boundary conditions is not a simple undertaking. Difficulties are associated both with the construction of nonreflecting boundary conditions (this is discussed below) and with the large distortions of the solution when the boundary conditions are poorly chosen. Since the mid-1970s theoreticians have been working on the development of nonreflecting boundary conditions for various equations and systems, including the equations of fluid dynamics. The present author [1 – 3] has revised the known approaches and proposed a number of new techniques that cover certain effects observed in applications. However, these studies do not provide complete results for spatially multidimensional systems of equations. The difficulties with constructing a nonreflecting boundary condition for fluid-dynamic problems in multidimensional cases emerge already for linearized Euler equations. Here the main difficulty is created by the reflection of acoustic waves that satisfy the convective wave equation. We know that it is impossible to obtain local nonreflecting conditions for the multidimensional wave equations. Approximations to the nonreflecting condition, which are the subject of voluminous literature, converge nonuniformly and slowly. For this reason we have decided to devote the first article in the series to boundary conditions for the classical wave equation (in a stationary medium) and the convective wave equation (in a moving medium). Here we formulate the general principles and the procedure for the construction of artificial boundary conditions. We also provide typical examples of boundary conditions. Techniques are suggested for assessing the properties of the boundary conditions. The proposed approaches are supported by theoretical considerations, but the substantiation is far exhaustive (although, on the other hand, it extends beyond the limits of our main subject). The convective wave equation, unlike the classical wave equation, is not discussed widely in the literature. We emphasize the specific features of the problems for the wave equation in a moving medium as opposed to a stationary medium. In other cases, conversely, the presence of convection does not prevent the reduction of the problem to a simpler format. We anticipate that the results of our study will be applied to systems of nonlinear equations of fluid dynamics and moreover to finite-difference methods for the solution of these systems. Note that the oscillation generation Translated from Prikladnaya Matematika i Informatika, No. 24, pp. 76 – 110, 2006. 282

1046–283X/07/1803–0282

c 2007


Springer Science+Business Media, Inc.

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

283

mechanism in these models is not entirely clear, and the information about wave propagation is imprecise. The nonlinearity of the problems substantially obstructs the application of nonlocal boundary conditions, which are in any event quite complex and hard to use. Finally, different systems of main equations — ranging from the scalar wave equation to finite-difference schemes describing viscous gas dynamics — require different numbers of boundary equations. In subsequent articles, while increasing the complexity of the models, we will follow the “inheritance principle,” i.e., in each stage we will combine and adapt the methods and boundary conditions obtained in previous stages. The present article is a theoretical survey. The proposed boundary conditions will be tested numerically in later stages — after being applied to more complex models. This is the subject of future articles. 1. Procedure for Investigation of Fluid-Dynamic Problems In this article we consider linearized equations of fluid dynamics describing the behavior of small perturbations of the parameters against a constant background. Subsequently, some of the results will be successfully extended to essentially nonlinear problems. Let T  U = ρ
u
v
p
be the perturbation vector of the primitive variables. In linear models, perturbations are denoted by attaching a prime to the corresponding variable symbol, while fluid-dynamic background parameters are denoted by plain letters and are assumed constant. Here and in what follows ρ is the density, u, v are the velocity components, p = ρc2/γ is the pressure, c the velocity of sound, γ the adiabatic index. The linearized Euler equations have the form ∂U ∂U ∂U + Cx + Cy = 0, ∂t ∂x ∂y

(1.1)

where  u ρ          0 u  Cx =     0 0       0 ρc2

 0       0 1/ρ      ,   u 0        0 u 0

  v 0 ρ 0               0 v 0 0        Cy =  .       0 0 v 1/ρ             2 0 0 ρc v

Suppose that the original problem is posed in the region Ω, which is typically the entire plane excluding the immersed body. The problem is mapped onto a bounded simulation region Ω
⊂ Ω which encompasses the body and the most important flow elements. As the region Ω
we basically consider a halfplane with the artificial boundary Γ
= {x = const, −∞ < y < ∞}, generally {x = 0}. If we are interested in a right boundary condition, then Ω
= {x < 0}, in case of a left boundary condition, Ω
= {x > 0}. To allow for oscillation generation in nonlinear models, we adopt a simplified linear description. Assume that oscillation sources with certain known characteristics exist inside the region Ω
and only there. The linearized Euler equations thus differ from (1.1) by the presence of a nonhomogeneous right-hand side with a finite-support function F : ∂U ∂U ∂U + Cx + Cy = F (x, y, t), ∂t ∂x ∂y

supp(x,y) F ⊂ Ω
.

284

L. V. D ORODNITSYN

This system is formally obtained from (1.1) by introducing initial conditions and conditions on a rigid body through the apparatus of generalized functions [4], but we will not link the right-hand side F with a specific problem. Since there are no sources in the neighborhood of the boundary, the reflecting properties of the boundary conditions are determined by the behavior of the solutions of the homogeneous equation (1.1). A differential consequence of the system of linear Euler equations (1.1) is the acoustic wave equation p
tt + 2up
xt + 2vp
yt + u2 p
xx + 2uvp
xy + v 2 p
yy − c2 (p
xx + p
yy ) = 0,

(1.2)

which is called the convective wave equation. (Here and in what follows subscripts denote the derivatives of the scalar functions with respect to the corresponding variables.) This equation and its particular cases are the topic of our article. For the analysis of linear equations, we mainly consider their partial solutions in the form of plane harmonic waves p
= exp{iωt − ikx − iy},

ω ≥ 0,

k,  ∈ R.

(1.3)

All the harmonic parameters are assumed real, because Eq. (1.2) is without dissipation and the waves do not decay. Assume that the boundary condition is a homogeneous linear equation of general form Lp
x=0

= 0.

(1.4)

We primarily consider differential-algebraic operators L with constant coefficients. The substitution of the function (1.3) in the boundary conditions (1.4) produces the Fourier transform of the operator L, which we denote

, ω). L(k, 2. Reflection of Acoustic Waves in a Moving Medium Consider a medium described by linearized two-dimensional Euler equations and moving with a constant velocity. The plane waves are assumed to interact with the fixed natural boundary {x = const}. The specific form of the boundary condition (1.4) is immaterial: it is only important that the reflected waves are plane. Due to the purely computational formulation of our problem, it is not a classical problem of physics. The laws of reflection of plane acoustic waves from an impermeable boundary are presented in the literature [5 – 7], but they deserve further examination in our context. Consider the case of horizontal flow with velocity (u, 0). The acoustic wave equation (1.2) takes the form   p
tt + 2up
xt + u2 p
xx − c2 p
xx + p
yy = 0.

(2.1)

The flow is subsonic: 0 < u < c. Choose the right boundary, which is called the outlet boundary if the velocity u is positive. The Fourier transform of Eq. (2.1) is −ω 2 + 2uωk − u2 k 2 + c2 (k 2 + 2 ) = 0.

(2.2)

Hence for the y -component of the wave vector  we obtain the expression c c sin θ ≡s= , ω c + u cos θ

where (k, ) =

 k 2 + 2 (cos θ, sin θ).

(2.3)

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

285

θ1 θ∗r

B

r

β

α

A

r

O

✛

α


M

1/M

r

O

✲

sin θ1

1

s

θ2

cos θ1

Fig. 1. Angles of incident and reflected waves in a moving medium. Here θ ∈ (−π, π] is the angle between the normal to the wavefront and the x-axis. Deviating from the established tradition, we call it in what follows the angle of incidence or the angle of reflection. When a wave is reflected, neither its frequency ω nor the vertical wavenumber  change, and thus s remains unchanged. Denote the Mach number by M = u/c and consider the function (2.3): s(θ) =

sin θ ; 1 + M cos θ

|s(θ)| ≤ s(θ∗ ) = √

1 , 1−M2

θ∗ = arccos(−M).

(2.4)

When |s| < s(θ∗ ) Eq. (2.3) has two roots for θ : the angle of incidence θ1 and the angle of reflection θ2 . The geometrical interpretation of the relationships between the parameters of incident and reflected waves for s > 0 is shown in Fig. 1. On a unit circle centered at O we have marked the arcs spanning the angles θ1 , θ2 and θ∗ . Draw a tangent through the point θ∗ that intersects the abscissa axis at the point A. We will show that the points θ1 and θ2 lie on the secant through the point A. From the triangle AOθ∗ in Fig. 1 and formula (2.4) it follows that |OA| = 1/cos(π−θ∗ ) = 1/M. Choose the 1 . Erect the vertical OB. Then point θ1 and denote the angle β = OAθ tan β = M|OB| =

sin θ1 . + cos θ1

M−1

Hence, by (2.4), s(θ1 ) = |OB|. The same reasoning holds for the second intersection point of Aθ1 with the circle, namely θ2 . Thus, there exists a limiting angle θ∗ > π/2 such that the incidence angles fall to its right (−θ∗ < θ1 < θ∗ ) and the reflection angles to its left (θ∗ < |θ2 | ≤ π). The incidence angle is “not equal” to the reflection angle, i.e., θ1 = π − θ2 (for s > 0). If 1 < s < s(θ∗ ) both angles θ1 and θ2 are obtuse and the wavevectors of both waves point leftward (see [5, 8]). Using relationships (2.4), we introduce the “generalized cosine” D=

 |M + cos θ| 1 − (1−M2 )s2 = . 1 + M cos θ

(2.5)

Hence for the angles of the incident and the reflected wave we obtain (opening the absolute value in different ways) cos θ1 =

D−M , 1−MD

cos θ2 = −

D+M . 1+MD

(2.6)

286

L. V. D ORODNITSYN

From (2.6) and also directly from (2.4) we obtain a useful relationship linking the incidence and reflection angles: tan(θ1 /2) tan(θ2 /2) = tan2 (θ∗/2) =

1+M . 1−M

A similar formula is given in [7]. In a moving medium the direction of the normal to the wavefront differs from the direction of wave propagation. From (2.2) we obtain an expression for the group velocity (see [9]) (∂ω/∂k, ∂ω/∂) ≡ C ≡ |C|(cos α, sin α) = c (M + cos θ, sin θ) .

(2.7)

In Fig. 1 the group velocity vectors of the incident and the reflected waves normalized by the velocity of sound c are represented by the segments −−→ C1 /c = O
θ1 ,

−−→ C2 /c = O
θ2 .

An interesting property of the velocities of the two waves is that they form equal angles with the horizontal axis: α1 = π − α2 = α

(s > 0).

In this sense, the incidence angle is equal to the reflection angle. Indeed, applying successively formulas (2.7), (2.4), and (2.5), we obtain tan α =

sin θ 1 + M cos θ =s· = ±s/D. M + cos θ M + cos θ

The sign + obviously corresponds to the incident wave and − to the reflected wave: tan α1 = s/D,

tan α2 = −s/D.

(2.8)

A corollary of this result is that always, for any vertical wavenumber , one acoustic wave propagates to the right and the other to the left [6]. Equation (2.1) requires one boundary condition that does not reflect the right acoustic wave. The reflection of acoustic waves from the left (inlet) boundary obeys the same geometry, but the roles of the incidence and reflection angles θ1 and θ2 are interchanged. Now consider an arbitrary flow velocity (u, v), remaining inside the subsonic range: u2 + v 2 < c2 . The acoustic wave equation takes the form (1.2) instead of (2.1). The analysis reduces to the previous if we change to a coordinate system moving with the velocity (0, v) [5] or, equivalently, introduce the modified frequency ω − v in (2.2). Then formula (2.3) remains valid when s is replaced with sˆ ≡ c/(ω − v) = s (1 − vs/c)−1 =

sin θ . 1 + M cos θ

(2.9)

The correspondence between the parameters s and sˆ is one to one, because the entire flow is subsonic. A singularity arises when the flow is supersonic, while the boundary is subsonic, i.e., the velocity has a subsonic normal component, but this case is not considered here.

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

287

The presence of a vertical velocity reduces to the case v = 0 by substituting sˆ from (2.9) for s in all the relationships between the wave parameters. The restrictions (2.4) are thus defined for the function sˆ(θ), the formula (2.5) changes,  D = 1 − (1−M2 )ˆ s2 , (2.10) while the relationships between the angles and D in expressions (2.5) and (2.6) are conserved. The geometrical properties of the group velocity vector (2.7) change: it now includes a vertical component v. 3. Properties of the Boundary Conditions for the Wave Equation In the previous section we have considered the reflection of plane waves from a plane boundary: the wavevector of the incident wave uniquely determines the wavevector of the reflected wave. The amplitude ratio of the two modes follows from boundary condition (1.4) for a constant operator L. Let us discuss the techniques for the construction of boundary conditions and investigate their reflecting properties. For definiteness we focus on the right boundary of the region. We start with the wave equation in a stationary medium:   p
tt − c2 p
xx + p
yy = 0.

(3.1)

In this case the horizontal wavenumbers are given by k1 =

ω 1−s2 , c

k2 = −

ω 1−s2 , c

(3.2)

where s = c/ω. It is better to change in (3.2) to the dimensionless wavenumber ζ = ck/ω : ζ1 =



 ζ2 = − 1−s2

1−s2 ,

(3.3)

(ζj = ckj /ω.). The partial solution of Eq. (3.1) is taken as the following combination of incident and reflected waves: p
= exp{iωt − ik1 x − iy} + R exp{iωt − ik2 x − iy}.

(3.4)

The number R, called the (amplitude) reflection coefficient, is determined from the right boundary operator L as R=−

1 , , ω) L(k .

2 , , ω) L(k

(3.5)

Consider a narrower class of boundary conditions: equations in which all the derivatives are of the same order

, ω) is a homogeneous function of its arguments: m ≥ 0. In this case, the transform of the boundary operator L(k,

, ω) = ω m L(k/ω, /ω, 1) = ω m L
(ζ, s). L(k,

(3.6)

The expression for the reflection coefficient (3.5) is rewritten using (3.6) as R=−

L
(ζ1 , s) . L
(ζ2 , s)

(3.7)

288

L. V. D ORODNITSYN

The reflection coefficient is thus independent of frequency but is a function of s: R = R(s). The value of s in turn depends on the incidence angle θ (in this case, s = sin θ), and thus, R = R(θ). A key fact [10] is that it is impossible to pose an exact local nonreflecting condition for the wave equation (3.1), i.e., there is no differential or differential-algebraic operator for which the numerator of (3.7) or (3.5) is identically zero, while the denominator does√not vanish. This is so because the Fourier transform of a local operator may not contain an irrational expression 1−s2 . The boundary conditions [10, 11] are based on the approximation of the radical by a rational function,  pm (s) 1−s2 ∼ r(s) = , qn (s) where pm (s), qn (s) are polynomials of degrees m, n. This leads to the dispersion equation ζ = r(s), and the reflection coefficient (3.7) is given by the formula √ √ qn (s) 1−s2 − pm (s) 1−s2 − r(s) √ = . R(s) = √ 1−s2 + r(s) qn (s) 1−s2 + pm (s)

(3.8)

√ Since 1−s2 is actually a function of s2 , for r(s) we generally take an even function of s. However, some approaches [12, 13] use nonsymmetrical expressions for r(s) and assume a prevailing orientation of the incident waves. We will give the simplest approximations of the radical and the resulting boundary conditions. Taylor expansion in powers of s gives √ √

1−s2 = 1 + O(s2 );

(3.9)

1−s2 = 1 − 0.5s2 + O(s4 ).

(3.10)

Hence we obtain two versions of the right boundary condition, dropping in (3.9), (3.10) the residual terms, substituting the expressions in (3.2), and taking the inverse Fourier transform: p
t + cp
x = 0;

(3.11)

p
tt + cp
xt − 0.5c2 p
yy = 0.

(3.12)

The reflection coefficients (3.8) given (3.11), (3.12) are respectively √ 1 − 1−s2 √ = − tan2 (θ/2); R=− 2 1 + 1−s R=−

2 √ 1 − 1−s2 √ = − tan4 (θ/2). 1 + 1−s2

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

289

The boundary conditions (3.11) and (3.12) are thus exact nonreflecting boundary conditions for waves incident along the normal to Γ
, and their accuracy is respectively O(s2 ) and O(s4 ) for small s. For applied purposes we can use various equivalent forms of the boundary condition (3.12), substituting in it the main equation (3.1) that holds near the boundary. In particular, replacing the second time derivative p
tt gives the equation p
xt + cp
xx + 0.5cp
yy = 0,

(3.13)

p
tt + 2cp
xt + c2 p
xx = 0.

(3.14)

and replacing the term p
yy in (3.12) gives

The reflection coefficient for boundary conditions (3.12) and (3.14) is obviously the same. Note that the boundary operators from (3.11) and (3.14) generate a family of nth order differential operators


Lp ≡

∂ ∂ +c ∂t ∂x

n p
= 0.

(3.15)

The boundary condition (3.15) produces a reflection coefficient (3.7) equal to R=−

n √ 1 − 1−s2 √ = − tan2n (θ/2). 1 + 1−s2

An important issue is the well-posedness of the problem of constructing the boundary conditions in twodimensional problems [14]. For a problem to be well-posed it is necessary (but not sufficient) that the denominator of the reflection coefficient (3.5) or (3.7) does not vanish in the entire definition domain. Otherwise the boundary condition would allow an incoming wave with an arbitrary amplitude. High-order differential boundary conditions lead to unstable solutions in most cases. An exhaustive class of √ 2 well-posed approximations of the radical 1−s for the wave equation has been constructed in [11]. In other models, convection and additional factors create further difficulties. This theme will be discussed in more detail in a future article on linearized Euler equations. Another essential factor should now be mentioned. For none of the local boundary conditions is the reflection coefficient small for all incidence angles. In particular, we always have total reflection of tangential (“sliding”) waves: R(θ) → 1, θ → ±π/2. For the wave equation and also in other cases it is undesirable to have boundary conditions that allow the reflection coefficient to be greater than 1 for some wave parameters. This would imply wave reflection with an increase of the amplitude. Therefore the boundary condition is required to satisfy the constraint1 R(θ) ≤ 1

(3.16)

for all incidence angles θ. The boundary conditions written previously satisfy (3.16). Below we return to estimation of the reflecting properties of the boundary conditions in the entire range of incidence angles. 1

The relationship of this property with the well-posedness property is not at all trivial.

290

L. V. D ORODNITSYN

We do not attempt to achieve a substantial improvement of the approximation of the exact nonreflecting condition. Besides the previously mentioned factors, the main reason is that in real problems of fluid dynamics the number of such boundary conditions is essentially less than in linear theory. Let us consider the boundary conditions for the three-dimensional wave equation   p
tt − c2 p
xx + p
yy + p
zz = 0.

(3.17)

Take the halfspace Ω
= {x < 0, −∞ < y, z < ∞} with the artificial right boundary Γ
= {x = 0, −∞ < y, z < ∞}. The Fourier transform of (3.17) in terms of the functions p
= exp{iωt − ikx − iy y − iz z} is the dispersion equation   −ω 2 + c2 k 2 + 2y + 2z = 0. Solving this equation for the horizontal wavenumber k, we find ζ1 =



1−s2 ,

 ζ2 = − 1−s2 ,

where ζj = ckj /ω,

s2 =

 c2  2 y + 2z . 2 ω

We see that the problem of constructing a nonreflecting boundary condition reduces to the two-dimensional case (3.3). Approximating the radical by (3.9) and (3.10), we obtain the equations (3.11) and   p
tt + cp
xt − 0.5c2 p
yy + p
zz = 0

(3.18)

instead of (3.12). It is remarkable that replacing in (3.18) the expression p
yy + p
zz from (3.17) we return to Eq. (3.14). The expressions for the reflection coefficients determined by these two boundary conditions remain as previously. Let us consider the general case of the equation of two-dimensional waves in a moving medium (1.2) with 2 u + v 2 < c2 . Consider only the right boundary, which is the outlet boundary for u > 0 and the inlet boundary for u < 0. The dispersion relationships for the two acoustic modes and the corresponding expressions for the wavenumbers have the form (see (2.10)) ω = uk + v + ck k1 =

 1+2/k 2 ;

(ω − v)(−u + cD) c2 − u2

ω = uk + v − ck k2 = −

or

ζ1 =

s D−M ; sˆ 1−M2

(3.19)

 1+2/k 2 ;

(ω − v)(u + cD) c2 − u2

or

ζ2 = −

s D+M . sˆ 1−M2

(3.20)

As previously, we study wave reflection from the right boundary by examining the combination of harmonics (3.4) and the reflection coefficient (3.5).

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

291

The simplest form of the right boundary condition is the one-dimensional equation. Setting  = 0 in the dispersion relationship (3.19), we obtain k = ω/(c+u) or ζ = 1/(1 + M). Hence we obtain the condition on the right boundary p
t + (c+u)p
x = 0.

(3.21)

The reflection coefficient is calculated from (3.7): R=−

1 − D + (1−M)vˆ s/c . 1 + D + (1−M)vˆ s/c

(3.22)

R is of first order of smallness in sˆ (and also in s and in the angle θ). It is preferable to use approximations of the dispersion relationship (3.19) in which the radical D from (2.10) is expanded in powers of sˆ. The expansion D ∼ r(ˆ s) contains only even powers of sˆ. This property is also transmitted to the reflection coefficient (3.7), which by analogy with the stationary case (3.8) is expressed as R(ˆ s) =

D − r(ˆ s) . D + r(ˆ s)

(3.23)

In boundary conditions obtained by inverting the Fourier transform the derivatives with respect to t are necessarily accompanied by convective terms in the direction y : p
t + vp
y . The initial approximation D ∼ 1 substituted in (3.19) gives the dispersion relationship k=

ω − v c+u

and the boundary condition p
t + (c+u)p
x + vp
y = 0,

(3.24)

for which the reflection coefficient (3.23) is R=−

1−D 1−M =− tan2 (θ/2) = − cot2 (θ∗/2) tan2 (θ/2) 1+D 1+M

(see (2.4), (2.5)). Here R is of second order of smallness in θ, and it moreover is independent of the vertical flow component. Let us allow for the next term in the Taylor expansion of D(ˆ s) from (2.10). In (3.19) the expression for k1 leads to k=

ω − v 0.5c2 − . c+u ω − v

Hence we obtain the boundary condition (see (3.12)) p
tt + 2vp
yt + v 2 p
yy + (c+u)(p
xt + vp
xy − 0.5cp
yy ) = 0.

(3.25)

292

L. V. D ORODNITSYN

An alternative form of (3.25), analogous to (3.13), is obtained by replacing the term p
tt in (3.25) from the main equation (1.2): p
xt + (c+u)p
xx + vp
xy + 0.5cp
yy = 0.

(3.26)

The reflection coefficient (3.23) of the boundary conditions (3.25) and (3.26) is given by  R=−

1−D 1+D

2

 =−

1−M 1+M

2 tan4 (θ/2),

i.e., R = O(ˆ s4 ). It is easy to see that the boundary conditions (3.24) and (3.26) satisfy the constraint (3.16). The simplest boundary condition (3.21) is less reliable. Specifically, for an inlet boundary, if the absolute value of the background flow velocity is close to the velocity of sound and the vertical velocity component v is relatively large, then in some range of incidence angles θ the reflection coefficient (3.22) may increase to infinity or at least become greater than 1 in absolute value. The examination of this case is quite involved and is accordingly relegated to Appendix A. 4. Propagation and Reflection of Nonplane Waves In previous sections we have examined the properties of plane harmonic waves that satisfy the convective wave equation with certain boundary conditions. The theory assumed decomposition of every solution in such waves. However, in typical multidimensional external problems in practice we observe predominance of waves of other kinds from concentrated initial perturbations or from various sources (multipoles, such as monopoles, dipoles, quadrupoles, etc.). Decomposition in plane waves in principle exists, but the coefficients have a characteristic distribution over the modes, and this distribution is not easily amenable to theoretical analysis. For this reason we cannot ignore cylindrical and spherical (in the three-dimensional case) waves. The analysis of nonplane waves is much more complex than the analysis of plane waves. In this context we compare the findings for waves of different types and, to the extent possible, reduce the results to appropriately chosen plane waves. We consider cylindrical harmonics in a two-dimensional medium and, for completeness, also spherical three-dimensional waves. Note that in some settings the three-dimensional case is simpler than the two-dimensional one. We thus consider a nonhomogeneous wave equation in a stationary medium p
tt − c2 ∆p
= f (x, t),

x ∈ Rd ,

(4.1)

where x is a point in the d-dimensional space, ∆ is the Laplace operator. Assume that some point x0 is the source of perturbations. We will give some important examples of these problems for Eq. (4.1). The first case is the Cauchy problem with concentrated initial perturbation. Equation (4.1) is defined for t > 0, the right-hand side is f ≡ 0, the initial condition is p
= aδ(x − x0 ),

p
t ≡ 0,

t = 0,

(4.2)

where δ(x) is the multidimensional Dirac function. Another example is an isotropic source of harmonic oscillations. Equation (4.1) holds for t ∈ R, the right-hand side is f (x, t) = aδ(x − x0 )eiωt .

(4.3)

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

293

Let us proceed to solve these problems. We start with the three dimensional case, i.e., d = 3, x = (x, y, z). The point oscillation source problem (4.1), (4.3) takes the form   p
tt − c2 p
xx + p
yy + p
zz = aδ(x − x0 )eiωt ,

x ∈ R3 .

(4.4)

The solution of (4.4) is p
=

A exp{iωt − i¯ r}, r

(4.5)

where r is the distance between the points x and x0 , r¯ = ωr/c is the dimensionless radius. The solution of the point initial perturbation problem (3.17), (4.3) is expressed [4] in terms of generalized functions and is formally representable in the form p
= AV3 (r, t),

(4.6)

i.e., it is spatially isotropic. In the two-dimensional space (d = 2, x = (x, y)) the point source problem (4.1), (4.3) takes the form   p
tt − c2 p
xx + p
yy = aδ(x−x0 , y−y0 )eiωt ,

(x, y) ∈ R2 .

(4.7)

Its solution is 

a




∞ 

p
= a
eiωt H0 (¯ r) = √ exp{iωt − i¯ aj (i¯ r)−j  , r} 1 + r¯ j=1

(4.8)

r) is the zeroth Hankel function (of 2nd kind) representable in an asymptotic series. Below we will where H0 (¯ examine the short-wave and/or large-distance approximation r¯  1. The Cauchy problem (3.1), (4.3), as in (4.6), has an isotropic solution p
= AV2 (r, t), where the generalized function V2 is given, e.g., in [4]. For practical purposes, the following model appears to be most useful for solving applied problems and formulating the relevant boundary conditions. Suppose that we know the type of the point source, but have no information about its location x0 and intensity. We only assume that x0 is inside the numerical region. Consider a simple particular case of problems for which an exact nonreflecting condition exists: the propagation of isotropic waves in 3D space from the source or from an initial perturbation located at a given point. Then the solution depends on two variables (r, t). The wave operator should be rewritten in spherical coordinates, and in this case it decomposes into a product of two operators: p
tt






− c ∆p = 2

∂ ∂ c +c + ∂t ∂r r



∂ ∂ c −c − ∂t ∂r r



p
.

The first factor describes waves moving from the center to the periphery, while the second factor describes waves in the opposite direction. Only waves of the first type are physically meaningful, and therefore the solutions of the corresponding problems, for instance the function (4.5), satisfy the equation Lp
≡ p
t + c(p
r + p
/r) = 0.

(4.9)

294

L. V. D ORODNITSYN

Suppose that the problem is posed in a region with an artificial boundary, in particular in the halfspace Ω
= {x < 0, −∞ < y, z < ∞}. Here x0 < 0. Then Eq. (4.9) should be used as the boundary condition on Γ
= {x = 0}. Also recall that in problems without initial conditions, such as (4.1), (4.3), the unique solution is determined using the radiation condition [4, 15], which is obtained from (4.9) as r → ∞. In broader classes of problems for the 3D wave equation (with mixed-type sources), and also in all 2D cases, local exact nonreflecting conditions do not exist. A sequence of approximate nonreflecting conditions of nth differential order of accuracy O(¯ r−2n−1 ), n = 1, 2, . . . , has been constructed in [16] for spherical waves (see below on the accuracy of the boundary condition). Another difficulty is associated not with the type, but with the localization of the source. For instance, if a source of isotropic spherical waves is located at a point other than the assumed point x0 , then in boundary condition (4.9) we need to redefine the radius r, because in its initial form it is no longer exact. The series of boundary conditions from[16] can be applied here also, because in some external region the solution is decomposable in multipoles centered at x0 . Yet we lose the main advantage of [16], where the accuracy of each solution in the sequence increases with the increase of the oscillation frequency ω in (4.4). The source localization error is similar to the error in the determination of the incidence angle of plane waves. Based on this consideration we can apply “dispersionless” boundary conditions with operators of the form L
(ζ, s) to model waves from concentrated sources. Let us discuss the theoretical estimation of the accuracy of an artificial boundary condition. The answer is quite simple for plane harmonic waves interacting with a plane boundary, when the condition is constant on the entire boundary. The boundary condition is estimated by the amplitude reflection coefficient (3.5) or (3.7) (in the 3D case the formulas have a similar form). The reflection coefficient can be introduced from the following consideration. Assume that in some region with the boundary Γ we have the main equation and the boundary condition Lp
|Γ = 0.

(4.10)

Suppose that the solution of the main equation p
= p1 (x, t), conventionally called the incident wave, is such that (4.10) defines the solution of a differential problem of the form p
= p1 (x, t) + Rp2 (x, t)

(4.11)

with constant and unique R. The function p2 should be of the same order as p1 , e.g., satisfy a normalization condition of the form |p1 (xΓ , t)| = |p2 (xΓ , t)|, where xΓ is an arbitrary boundary point x ∈ Γ. Thus R = R[p1 ] = −

Lp1 (xΓ , t) Lp2 (xΓ , t)

(4.12)

is the sought reflection coefficient. Examples when these manipulations can be undertaken in practice are quite rare. For a given function p1 we have to construct explicitly the function p2 , preferably unique for the whole class of operators L from (4.10); alternatively we have to compare the boundary conditions with one another. This requirement is satisfied by plane waves and boundary operators with constant coefficients. The following example also admits the introduction of reflection coefficients.

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

295

Consider the isotropic wave equation in a region bounded by a sphere of radius r1 : p
tt

 −c

2

p
rr

2 + p
r r

 = 0,

0 ≤ r0 < r < r1 .

Following (4.11), we write the solution in the form (see (4.5)) p
= p1 + Rp2 ≡

 eiωt  exp{−i(¯ r − r¯1 )} + R exp{i(¯ r − r¯1 )} , r

(4.13)

where r¯1 = ωr1 /c. The value of R depends on the boundary condition for r = r1 . If at this point we specify Eq. (4.9), then obviously R = 0. We seek reflection coefficients for the case of boundary conditions from Sec. 3. In this case, the derivative with respect to x, i.e., with respect to the normal to the boundary, corresponds to the derivative with respect to r. For condition (3.11) Lp
≡ p
t + cp
r = 0 we obtain from (4.13) Lp1 = iωp1 + c

∂p1 cp1 =− , ∂r r

Lp2 = 2iωp2 −

cp2 . r

Hence, by (4.12), R=

i 1 + O(¯ r1−2 ). =− 2i¯ r1 − 1 2¯ r1

(4.14)

Condition (3.14) Lp
≡ p
tt + 2cp
rt + c2 p
rr = 0 yields Lp1 =

2c2 p1 , r2 R=

Lp2 = −

2¯ r12

 4  2 2 ω r + icωr − c2/2 p2 , 2 r

1 1 = 2 + O(¯ r1−3 ). + 2i¯ r1 − 1 2¯ r1

(4.15)

Formulas (4.14) – (4.15) demonstrate that the order of smallness of the reflection coefficient increases with the increase of the order of accuracy in s of the boundary conditions for a multidimensional wave equation. For most problems and most functions p1 (x, t) there is no explicit expression for the reflected wave p2 (x, t) from (4.11) or else this wave is represented as a decomposition (distribution) in known functions, where the decomposition coefficients depend on the choice of the boundary condition. In this case, we can calculate the numerator of the fraction in the right-hand side of (4.12), but the denominator becomes undefined. If the boundary condition is well-posed, then we can assert in a certain sense that the denominator in (4.12) does not vanish and conserves the prior order of magnitude. Thus, the boundary operator should primarily produce a small quantity when acting on a “trial” incident wave p1 The order of accuracy of the boundary condition is the ratio of the orders of magnitude of the numerator and the denominator of the fraction in (4.12).

296

L. V. D ORODNITSYN

We apply these considerations to investigate the reflecting properties of operators (3.6) in the problem of interaction of isotropic cylindrical waves (4.7) with the right boundary {x = 0}. Assume for definiteness that the source is at the point (x0 , 0), where x0 < 0. Consider the numerator of (4.12). In this case p1 = p
from (4.8). We assert that   Lp
(0, y, t) = ω m L
(ζ1 , s)p
+ O ω m r¯−1 p
, where ζ1 = −x0 /r,

(4.16)

 s = y/r,

r=

x20 + y 2

(it easy to verify that the first equality in (3.3) holds). In other words, the main contribution to the action of the operator L on a cylindrical wave is associated with the term describing the action of the same operator on the corresponding plane wave. In the neighborhood of the point (0, y1 ) the solution (4.8) is close to the function p
 A(y1 ) exp{iωt − ik1 x − iy},

where k1 = −

ωx0 , cr

=

ωy1 . cr

The bound (4.16) follows from the theorem in Appendix B. With small incidence angles (s  1) this result apparently can be strengthened. The order of accuracy of
L (ζ1 , s) in s is increased by increasing the order of smallness in r¯ of the residual term in (4.16). This is observed for the case of operators (3.15). To bound the reflection coefficient we need information about the denominator of the fraction in (4.12), i.e., the properties of the reflected wave p2 . Its scale (characteristic amplitude) obviously changes along the boundary, because (4.16) depends on s, and for typical boundary conditions the value of Lp
is small with normal incidence. We can accept the traditional physical conjecture that allows local replacement of both incident and reflected waves by effective plane waves. In this context, it is helpful to pass from the constant reflection coefficient R = const to the same reflected wave amplitude |p2 | for x = 0 and the variable coefficient R = R(y). Formula (4.16) suggests a stronger assumption. Every boundary operator from the class (3.6) acts on the function p2 approximately in the same way as on a plane harmonic wave with a wavevector corresponding to the laws of reflection of the incident wave at the given point. As a result, we may take the reflection coefficient from (4.12) close to R(s) defined in (3.7). Waves in a moving medium constitute a more complicated topic. Assume that spherical waves in a medium moving with the velocity u = (u, 0, 0) with |u| < c are generated by a stationary point source: p
tt + 2up
xt + u2 p
xx − c2 (p
xx + p
yy + p
zz ) = aδ(x − x0 )eiωt ,

x ∈ R3 .

(4.17)

Problem (4.17) has the solution p
= A

exp{iωt − i¯ r} , r(1 + M cos θ)

(4.18)

where (for simplicity we take y0 = z0 = 0) x − x0 r= = M + cos θ

 y2 + z2 . sin θ

(4.19)

The phase surfaces of (4.18) are the eccentric spheres shown in Fig. 2. The wave geometry is also illustrated by Fig. 1, where the radius r = 1 is measured from the center of the circle O, while the source is at the point O
. Note that the wave amplitude in (4.18) changes with the direction of radiation.

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

297

Fig. 2. Phase lines of waves from a point source in a moving medium.

Now consider waves radiated by a stationary cylindrical source with u = (u, 0):   p
tt + 2up
xt + u2 p
xx − c2 p
xx + p
yy = aδ(x−x0 , y−y0 )eiωt ,

(x, y) ∈ R2 .

(4.20)

Problem (4.20) can be solved by the descent method [15], i.e., by integrating the 3D solution (4.18) over the coordinate z. Making a change of variables, we reduce the problem to the stationary medium case (4.8): p
= a
H0



 r} 1 + M cos θ A exp{iωt − i¯ . r ¯ exp{iωt − iuωx/c2 }  √ √ 1 − M2 r 1 + M cos θ

(4.21)

Wave radiation is anisotropic. The Cauchy problems with a point initial perturbation in 2D and 3D spaces are transformed into similar problems for a stationary medium (4.1), (4.3) if we pass to a coordinate system moving with velocity u. The solution of the equation p
tt + 2up
xt + u2 p
xx − c2 ∆p
= 0 with the initial condition (4.3) has the form p
= AVd (|x − x0 − ut|, t),

d = 2, 3,

(4.22)

where Vd are the generalized functions defined for the stationary case. We thus observe isotropic waves from a moving center. We will not examine in detail the properties of boundary conditions for problems for waves from a source in a moving medium. Exact local nonreflecting conditions apparently do not exist in this case. Note that the approach of [16] is extended in [17] to the construction of boundary conditions in convective problems for spherical and cylindrical waves. Finally, local substitution of effective plane waves is also admissible for the investigation of wave interactions with a plane boundary {x = 0}. For boundary operators without low-order derivatives we can probably use formula (4.16). Its proof becomes more complex than for a stationary medium, or possibly requires generalization of the procedure. In what follows we estimate the quality of the boundary conditions and optimize them in the framework of the approximation involving plane waves, plane boundaries, and constant coefficients.

298

L. V. D ORODNITSYN

5. Integral Reflection Coefficients To estimate the reflecting properties of a boundary condition, we use the reflection coefficients defined by this condition. However, as we have seen before, a universal reflection coefficient does not exist in most problems. It depends on the choice of the incident wave, and is moreover local in many cases. For practical purposes we need to have a general measure of reflection of waves with various frequencies and incidence angles. It is thus necessary to introduce an integral reflection coefficient. Let Ω
= {x < 0} be the simulation region, Γ
= {x = 0} the boundary, and the space is two- or threedimensional. For plane waves and boundary conditions from the class (3.6), the reflection coefficient depends only on the incidence angle, R = R(θ). The integral reflection coefficient can be obtained by averaging R(θ) over θ from the entire range of incidence angles. The main question is how to choose the weighting function, which is determined by the wave type. We consider waves from concentrated perturbations and sources, replacing them with local plane waves in the neighborhood of the boundary. We allow for the dependence of the amplitude of these waves on the traversed distance and possibly the radiation angle. For averaging it makes sense to use the energy characteristics, and not the amplitudes of the modes; for instance, the energy or the energy flux [5], which is not the same. We will use energy density. For the incident wave the energy density is proportional to the squared absolute value of the pressure perturbation |p
|2 , and the energy density of the reflected wave in turn is proportional to |R|2 . Both parameters in the final analysis are expressed in terms of the incidence angle θ, and the result depends on the wave type. 5.1. Reflection of Waves in a Stationary Medium. We will apply the principles developed above to specific examples of differential equations. We start with an isotropic point source for a two-dimensional wave equation (4.7). The point (x0 , y0 ) is inside the simulation region, i.e., x0 < 0. Note that the specific location of the source does not change the result. The perturbation amplitude |p
| follows by (4.8) the asymptotic law    c A |p
| = √ 1 + O . ωr r On the right boundary {x = 0} the radius is r = |x0 |/cos θ and √ |p
|  |p
(θ)| = A
cos θ.

(5.1)

Define the integral reflection coefficient RS for problem (4.1), (4.3) using (5.1):

RS2

1 = π

π/2 −π/2

|p
(θ)|2 |R(θ)|2
dθ |p (0)|2

1 = π

π/2 |R(θ)|2 cos θdθ.

(5.2)

−π/2

The normalization of the integral is the ratio of intensities of oblique waves and waves incident at a right angle. In this way we distinguish between cylindrical and plane waves. If the reflection coefficient is R(θ) ≡ 1 (or if we estimate incident instead of reflected waves), then the integral reflection coefficient is less than 1: by (5.2), RS2 = 2/π,

RS = 0.7979.

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

299

We write out the integral reflection coefficients for known examples of boundary conditions with their R(θ): R(θ) = − tan2 (θ/2);

RS2 =

38 − 4, 3π

RS = 0.1787,

R(θ) = − tan4 (θ/2);

RS2 =

2642 − 8, 105π

RS = 0.09635.

Here we treat the integral reflection coefficient as a functional of the function R(θ) defined on an interval, i.e., RS [R]. The mean energy density of reflected waves divided by the mean energy density of incident waves equals RS2 [R]/RS2 [1]. When the reflection coefficient R is expressed in terms of the dimensionless wavenumber s = sin θ, the integral (5.2) transforms to RS2

1 = π

1 |R(s)|2 ds. −1

For the boundary condition based on a rational approximation of the radical, formula (3.8) gives 1 RS2 = π

1  −1

2 √ qn (s) 1−s2 − pm (s) √ ds. qn (s) 1−s2 + pm (s)

(5.3)

√ In the literature the radical 1−s2 is sometimes approximated by a rational function using the least squares method [11]. The problem is to find polynomials pm (s) and qn (s) of given degrees from the condition 1  min

pm ,qn : qn (0)=1 −1

2  qn (s) 1−s2 − pm (s) ds.

We see from (5.3) that the optimization problem should be posed somewhat differently, although minimization of (5.3) is a technically difficult undertaking. However, this is not included in our agenda. The integral reflection coefficient has been defined for a single point source (4.3). This measure can be defined with equal success for several or even distributed sources given that they are not coherent. Since the contribution of each source to the reflected wave energy is proportional to the integral (5.2), we obtain the same expression (5.2) regardless of the specific weights of the sources in the sum. Formula (5.2) is true not only for an isotropic source (4.3), but also for a randomly oriented multipole. Consider the cylindrical dipole f (x, y, t) = a

∂δ (x−x0 , y−y0 )eiωt , ∂xα

(5.4)

where ∂δ/∂xα is the derivative of the delta function with respect to the direction α, with the axis uniformly distributed in the angle α ∈ [−π/2, π/2]. The asymptotic behavior of |p
| is such that (see (5.1)) √ |p
|  |p
(θ, α)| = A
cos θ cos(θ−α).

(5.5)

300

L. V. D ORODNITSYN

Generalizing (5.2), we define the integral reflection coefficient for the source (5.4) as 1 π

RS2 =





π/2

dα

−π/2



π/2

−π/2 π/2


−π/2

|R(θ)|2 |p
(θ, α)|2 dθ .

|p (0, α)| dα 2

Substituting |p
| from (5.5) and changing the order of integration in the double integral, we obtain

RS2 =

1 π



π/2

−π/2

 |R(θ)| cos θdθ 2



π/2

−π/2

π/2

−π/2

cos2 (θ−α)dα

cos2 αdα

1 = π

π/2 |R(θ)|2 cos θdθ. −π/2

We again end up with formula (5.2) for the integral reflection coefficient. If instead of a constant-power source we consider the Cauchy problem (3.1), (4.3), the result remains the same. The perturbation amplitude p
at the instant of crossing the boundary has the asymptotic expression (5.1) (for initial functions weakly convergent to δ(x−x0 )), and the definition (5.2) of the integral reflection coefficient thus remains valid. In practice, the problem is posed not on a halfplane, but in a bounded region, for instance a rectangle whose right boundary is of finite length in y. This means that waves propagating at angles θ close to ±π/2 do not always reach the right boundary. Such incidence angles are not excluded, because the source (x0 , y0 ) may be arbitrarily close to the boundary, but their probability is smaller than the probability of small θ. It is hardly possible to allow for this effect precisely, but for most functions R(θ) the integral reflection coefficient (5.2) may be used as an upper bound on the true reflecting properties of the boundary condition. We now define the integral reflection coefficient for the general case. Following the principle used in the derivation of (5.2), we define the integral reflection coefficient by the formula RS2

1 = |Ω+ |

 |R(Ω)|2 Ω+

|p
(Ω)|2 dΩ. |p
(Ω0 )|2

(5.6)

Here Ω+ is the (solid) angle of incident waves, |Ω+ | is the magnitude of this angle, Ω is the current direction, Ω0 is the direction of normal incidence; p
is taken on the boundary of the region. Three-dimensional waves from a point isotropic source are described by Eq. (4.4). The right boundary Γ
is the plane {x = 0, −∞ < y, z < ∞}. Instead of the two-dimensional asymptotic expression, we have a different equality that follows from (4.5): |p
| = A/r = A
cos θ.

(5.7)

Let us find an expression for the integral reflection coefficient, applying formula (5.6) with Ω+ = {Ω : 0 ≤ θ ≤ π/2}, |Ω+ | = 2π. Typical boundary conditions give R = R(θ), and substituting (5.7) in (5.6) we obtain

RS2

π/2 = |R(θ)|2 cos2 θ sin θdθ. 0

(5.8)

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

301

We give the integral reflection coefficients for the boundary conditions discussed previously: R(θ) = 1;

RS2 = 1/3,

R(θ) = − tan2 (θ/2);

RS2 =

R(θ) = − tan4 (θ/2);

RS2 = 61 − 88 ln 2,

RS = 0.5774,

25 − 12 ln 2, 3

RS = 0.1248, RS = 0.05521.

The differences from the 2D case have obvious reasons. Comparing (5.2) and (5.8) we note that in the 3D space the role of both large (θ ≈ π/2) and small angles (θ ≈ 0) is reduced. 5.2. Wave Reflection in a Moving Medium. Consider the simplest case, when the medium moves with velocity u ∈ (−c, c) parallel to the x-axis. Let us elucidate how the waves from a stationary spherical source (4.17) behave near the right boundary. From (4.18) and (4.19) we have |p
| =

M + cos θ A = A
. r(1 + M cos θ) 1 + M cos θ

(5.9)

Let the reflection coefficient be R = R(θ), 0 ≤ θ ≤ θ∗ . Define the integral reflection coefficient by (5.6) with solid angle Ω+ = {Ω : 0 ≤ θ ≤ θ∗ } of magnitude |Ω+ | = 2π(1+M). Substituting (5.9) in (5.6) we obtain 1 RS2 = 1+M



θ∗ |R(θ)|

2

M + cos θ 1 + M cos θ

2 sin θdθ.

(5.10)

0

Table 1 lists the values of RS for typical examples of reflection coefficients R(θ) for various Mach numbers M. For M = 0 the results naturally coincide with the case of a stationary medium (5.8). The fact that for R(θ) ≡ 1 the value of RS varies with the Mach number M is explained by the nonequilibrium intensity of wave radiation over the direction (see (4.18)). The situation is different if we consider the Cauchy problem (4.17) with zero right-hand side and (4.3). By (4.22) the generalized amplitude, as in a stationary medium, depends on the radiation angle and has the asymptotic expression |p
|  A/r = A
(M + cos θ). In the end, the formula for the integral reflection coefficient differs from (5.10):

RS2

1 = (1+M)3

θ∗ |R(θ)|2 (M + cos θ)2 sin θdθ.

(5.11)

0

√ For R(θ) ≡ 1 we obtain the same integral reflection coefficient for all Mach numbers, RS = 1/ 3. The results for main R(θ) are given in Table 1. The dependence of RS on the Mach number is more pronounced than with a point source (5.10). Boundary conditions of the same type “perform” less well on subsonic inlet than on subsonic outlet.

302

L. V. D ORODNITSYN

Table 1. Integral Reflection Coefficients Rs[R] for Various Wave Types and Mach Numbers R(θ) Model

Source 2D (5.14) Cauchy 2D (5.13) Source 3D (5.14) Cauchy 3D (5.14)

M

θ∗ θ tan4 2 2

1

1−M θ tan2 1+M 2

0

0.7979

0.1787

0.09635

1/2

0.8372

0.1608

0.08268

−1/2

0.7734

0.1878

0.1043

0

0.7979

0.1787

0.09635

1/2

0.6090

0.1281

0.06250

−1/2

0.6540

0.2129

0.1217

0

0.5774

0.1248

0.05521

1/2

0.6744

0.1167

0.04661

−1/2

0.5185

0.1258

0.05907

0

0.5774

0.1248

0.05521

1/2

0.5774

0.07982

0.02836

−1/2

0.5774

0.1572

0.07800

cot4

Our analysis of 2D waves in a moving medium starts with the Cauchy problem (4.2) for the homogeneous equation (2.1). The perturbation amplitude has the asymptotic expression |p
|  A

√

√ r = A
M + cos θ.

(5.12)

The reflection coefficient R(θ) is defined for −θ∗ ≤ θ ≤ θ∗ . Substituting (5.12) in the formula (5.6) for the integral reflection coefficient, we obtain a generalization of (5.2):

2 RS2 = RSC

1 = ∗ 2θ

θ∗ −θ∗

|p
(θ)|2 1 |R(θ)|2
dθ = ∗ 2 |p (0)| 2θ (1+M)

θ∗ |R(θ)|2 (M + cos θ)dθ.

(5.13)

−θ∗

The values of RS [R] are summarized in Table 1. The integral reflection coefficient again essentially depends on the Mach number. We now consider waves radiated by a stationary cylindrical source (4.20). From (4.21) and the first equality in (4.19) it follows that on the right boundary A |p |  √ √ = A
r 1 + M cos θ




M + cos θ 1 + M cos θ

1/2 .

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

303

Instead of (5.13) we now have the following formula for the integral reflection coefficient:

2 RS2 = RS1

1 = ∗ 2θ

θ∗ |R(θ)|2 −θ∗

M + cos θ dθ. 1 + M cos θ

(5.14)

The values of RS for various coefficients R(θ) show a fairly weak dependence on the Mach number (Table 1). It follows from the preceding discussion that, for a variety of reasons, the integral reflection coefficient does not describe accurately the extent of wave reflection in a particular problem. At the same time RS [R] provides a universal measure of the reflecting properties of a given boundary condition. Assume for instance, that in some two-dimensional problem we have both an initial perturbation and an oscillation source, which moves more slowly than the flow and decays over time. Then we can assume with a high degree of confidence that the relative energy 2 2 from (5.14). density of the reflected waves falls between the values of the expressions RSC from (5.13) and RS1 Two variants of the integral reflection coefficients thus provide a lower and an upper bound for wave reflection in a more general case. Appendix A. Analysis of Reflection Coefficients for the Simplest Boundary Condition Let us determine for what background flow parameters the boundary condition (3.21) may or may not be a priori well-posed or may produce wave reflection with increasing amplitude. The expression for the reflection coefficient (3.22) is rewritten using formulas (2.9) and (2.5), and in new notation: R(θ) = −

(1−Mx )(1 − cos θ + My sin θ) , (1+Mx )(1 + cos θ) + (1−Mx )My sin θ

(A.1)

where Mx = u/c, My = v/c. We also introduce the scalar Mach number M≡



M2x + M2y < 1.

The coefficient (A.1) depends on the incidence angle and also on the two flow parameters M x and M y : R = R(θ; Mx , My ). Classifying these parameters, we identify three properties of the function R(θ): 1. For Mx ≥ 0, the coefficient R(θ) from (A.1) is defined for −θ∗ ≤ θ ≤ θ∗ and satisfies property (3.16). This implies that at the outlet boundary the boundary condition (3.21) is well-posed and guarantees that the wave amplitude decreases during reflection. To prove this property, we first show that the denominator in (A.1) does not vanish on [−θ∗ , θ∗ ]. We present the proof for the case My ≥ 0. Noting that cos θ ≥ −Mx , we obtain a lower bound on the denominator: (1+Mx )(1 + cos θ) + (1−Mx )My sin θ ≥ (1+Mx )(1−Mx ) − (1−Mx )My = (1−Mx )(1+Mx −My ). For Mx ≥ 0 the second factor is positive (as well as the first), and the overall expression is therefore also positive. The case Mx < 0 requires a different approach and is considered below.

304

L. V. D ORODNITSYN

If the denominator is positive on the entire interval, the real function R(θ), is continuous on this interval. We will show that R ≥ −1. To this end, subtracting the denominator of (A.1) from the numerator we obtain −2(Mx + cos θ) ≤ 0, and equality is attained only at the endpoints θ = ±θ∗ . The inequality proves the required property. It remains to check for what parameter values the coefficient R may take the value +1. The equation R = 1 is equivalent to the expression 1 + Mx cos θ + (1−Mx )My sin θ = 0.

(A.2)

The left-hand side of (A.2) is lower bounded by 1 − M2x − (1−Mx )My = (1−Mx )(1+Mx −My ), which is strictly positive for Mx ≥ 0. Hence it follows that R < 1. We have thus shown that R from (A.1) satisfies property (3.16). 2. For Mx < 0, a necessary and sufficient condition for the existence of an angle θ ∈ [−θ∗ , θ∗ ] such that in its neighborhood |R(θ)| → ∞ is given by the inequality M2y ≥

1 + Mx . 1 − Mx

(A.3)

The values of M x and M y satisfying (A.3) exist if and only if M ≥ M∗ ≈ 0.7625.

(A.4)

In other words, with Mach number M < M∗ the denominator in (A.1) never vanishes and the necessary wellposedness condition does not break down. The denominator in (A.1) vanishes if and only if tan(θ/2) ≡

sin θ 1 + Mx < 0. =− 1 + cos θ (1−Mx )My

(A.5)

The unique root of this equation inside the definition region exists when tan2 (θ/2) ≤ tan2 (θ∗/2) =

1 + Mx . 1 − Mx

Squaring the right-hand side of (A.5), we obtain (A.3). Let us find for what values of M x the condition (A.3) may hold. Combining (A.3) with the constraint on M, we obtain the double inequality 1 + Mx ≤ M2y < 1 − M2x . 1 − Mx

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

305

Since equality is attainable, this expression gives a necessary and sufficient inequality for M x . It is easily solved: −1 < Mx < 0, i.e., no new constraints arise. This completes the proof of (A.3). We will now show that inequality (A.3) is attainable only for sufficiently large M. Expression (A.3) is equivalent to the following: M2 ≥ M2x +

1 + Mx ≡ Φ(Mx ). 1 − Mx

Finding min Φ(x) ≈ 0.5814 = (M∗ )2

x∈[−1,0]

and noting that inequality (A.3) is unimprovable, we conclude that it holds when (A.4) is satisfied. 3. If M < M∗ ≈ 0.7603,

(A.6)

then property (3.16) holds a fortiori, i.e., for moderate M the reflected wave amplitude always decreases. It follows from the preceding that for M < M∗ the function R(θ) in (A.1) is continuous. We have shown that in this case R(θ) ≥ −1. We will now elucidate for what parameter values the function may take the value +1. We again obtain Eq. (A.2) and note that it has a solution only if the sum of the squares of the coefficients of cos θ and sin θ is not less than 1. Otherwise, it is guaranteed that R does not attain +1. Let  Mx = Mξ, My = M 1−ξ 2 . The sum of squares of the coefficients in the left-hand side of (A.2) is M2x + (1−Mx )2 M2y = M2 − 2M3 ξ(1−ξ 2 ) + M4 ξ 2 (1−ξ 2 ) ≡ F (ξ; M). Separately maximizing the coefficients of different degrees of M over all ξ ∈ [−1, 0], we obtain √ 4 3 3 M4 F (ξ; M) ≤ M + M + = Ψ(M). 9 4 2

The inequality Ψ(M) < 1 is in turn equivalent to (A.6). Thus, (A.6) is a coarse bound, but it does not make much sense to refine the maximum admissible Mach number because M ∗ is close to M ∗ from (A.4). Appendix B. General Property of Operators with Constant Coefficients We will prove the assertion from Sec. 4 about the action of dispersionless operators on a cylindrical harmonic at the point (x, y). For convenience we change to dimensionless variables.

306

L. V. D ORODNITSYN

Consider an nth order operator from the class (3.6):

Lu = L(n) u =

n  k=0

Lk

∂ n−k u , ∂tn−k

Lk =

k 

ajk

j=0

∂k . ∂xj ∂y k−j

(B.1)

The transform of the operator L(n) in application to the functions u(x, y, t) = exp{it − iζx − isy} is L
(n) (ζ, s) =

n 

(i)n−k L
k (ζ, s),

where L
k (ζ, s) = (−i)k

k 

ajk ζ j sk−j .

j=0

k=0

Here and below subscripts without parentheses identify a homogeneous polynomial of the corresponding degree, while subscripts in parentheses identify a polynomial of general form. The cylindrical wave is represented by the function u(x, y, t) = eit H0 (r),

where r =

 x2 + y 2 ,

(B.2)

H0 is the zeroth Hankel function. We rewrite (4.16) in the form L(n) u(x, y, t) = L
(n) (ζ1 , s)u + O(r−1 u), r → ∞,

(B.3)

where u is defined by (B.2), ζ1 =

x , r

s=

y . r

Instead of proving directly the bound (B.3), we prove by induction two equalities that hold for even and odd n and imply (B.3). An operator (B.1) of order n+1 is representable in the form L(n+1) u =

∂ L u + Ln+1 u, ∂t (n)

Ln+1 = L(1) n

∂ ∂ + L(2) . n ∂x ∂y

(B.4)

The first term (differentiation with respect to t) obviously does not affect the bound (B.3) for the functions (B.2). It is only necessary to verify the assertion for the homogeneous spatial operator Ln+1 consisting of two similar terms. Below we use the first Hankel function H1 (r), which satisfies the following identities: d H0 (r) = −H1 (r), dr ∂r x = , ∂x r

d 1 H1 (r) = H0 (r) − H1 (r), dr r ∂r y = . ∂y r

(B.5)

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

307

Note also the relationship between Hankel functions of 2nd kind:   H1 (r) = iH0 (r) + O r−1 H0 (r) .

(B.6)

Also note that L
k (ζ1 , s) = r−k L
k (x, y). We first show that formula (B.3) holds for first-order operators. For L = ∂/∂x, using the identities (B.5) and the asymptotic expression (B.6), we obtain     ∂ x ix H0 (r) = − H1 (r) = − H0 (r) + O xr−2 H0 (r) = −iζ1 H0 (r) + O r−1 H0 (r) . ∂x r r Since L
(ζ, s) = −iζ, we obtain (B.3). Theorem. The following relationships hold: Ln H0 (r) = r−n L
n (x, y)H0 (r) + r−2n+2 P(2n−4) (x, y)H0 (r) + r−2n+1 Q(2n−2) (x, y)H1 (r), n = 2m,

(B.7)

m = 1, 2, . . .

Ln H0 (r) = −ir−n L
n (x, y)H1 (r) + r−2n+2 P(2n−3) (x, y)H0 (r) + r−2n+1 Q(2n−3) (x, y)H1 (r), n = 2m + 1,

(B.8)

m = 1, 2, . . .

Proof. Using formulas (B.5), we easily verify that (B.7) holds for n = 2. Then the proof is by the following scheme. Given the truth of (B.7) for n = 2m, we obtain equality (B.8) for n+1 = 2m+1: Ln+1 H0 (r) = −ir−n−1 L
n+1 (x, y)H1 (r) + r−2n P(2n−1) (x, y)H0 (r) + r−2n−1 Q(2n−1) (x, y)H1 (r),

(B.9)

and then (B.7) for n+2 = 2m+2: Ln+2 H0 (r) = r−n−2 L
n+2 (x, y)H0 (r) + r−2n−2 P(2n) (x, y)H0 (r) + r−2n−3 Q(2n+2) (x, y)H1 (r).

(B.10)

Relying on (B.4), consider the operator Ln+1 u =

∂ Ln u. ∂x

Differentiating (B.7), we obtain dH0 ∂r dH0 ∂r ∂ Ln H0 (r) = r−n L
n (x, y) + r−2n+2 P(2n−4) (x, y) ∂x dr ∂x dr ∂x + r−2n+1 Q(2n−2) (x, y) + r−2n+2

∂ dH1 ∂r + r−n L
n (x, y)H0 (r) dr ∂x ∂x

∂ ∂ P(2n−4) (x, y)H0 (r) + r−2n+1 Q(2n−2) (x, y)H1 (r) ∂x ∂x

(B.11)

308

L. V. D ORODNITSYN

− nr−n−1

∂r
∂r Ln (x, y)H0 (r) − (2n−2)r−2n+1 P(2n−4) (x, y)H0 (r) ∂x ∂x

− (2n−1)r−2n

∂r (x, y)H1 (r). Q ∂x (2n−2)

Substituting the first three formulas from (B.5) and introducing new notation for the derivatives of polynomials with respect to x, we obtain ∂ Ln H0 (r) = −r−n−1 xL
n (x, y)H1 (r) − r−2n+1 xP(2n−4) (x, y)H1 (r) ∂x + r−2n xQ(2n−2) (x, y)H0 (r) − r−2n−1 xQ(2n−2) (x, y)H1 (r) + r−n M
n−1 (x, y)H0 (r) + r−2n+2 A(2n−5) (x, y)H0 (r) + r−2n+1 B(2n−3) (x, y)H1 (r) − nr−n−2 xL
n (x, y)H0 (r) − (2n−2)r−2n xP(2n−4) (x, y)H0 (r) − (2n−1)r−2n−1 xQ(2n−2) (x, y)H1 (r).

(B.12)

For n = 2 the term with A(2n−5) is missing. Let us discuss the first term in the right-hand side of (B.12). The transform of the operator (B.11) is L
n+1 (ζ, s) = −iζL
n (ζ, s). Therefore −r−n−1 xL
n (x, y)H1 (r) = −ir−n−1 L
n+1 (x, y)H1 (r).

(B.13)

Let us transform the terms with H0 (r): n/2

r−n M
n−1 (x, y) = r−2n (x2 +y 2 )

M
n−1 (x, y).


(x, y). We The number n/2 = m is an integer, and the function multiplying r−2n is therefore a polynomial N2n−1 do the same with the other terms. In the end,

r

−2n

 n/2 xQ(2n−2) (x, y) + (x2 +y 2 ) M
n−1 (x, y) + (x2 +y 2 )A(2n−5) (x, y) − n(x +y ) 2

2

n−2 2

xL
n (x, y)

= r−2n P(2n−1) (x, y)H0 (r).

 − (2n−2)xP(2n−4) (x, y) H0 (r) (B.14)

A RTIFICIAL B OUNDARY C ONDITIONS FOR T WO -D IMENSIONAL E QUATIONS OF F LUID DYNAMICS . 1

309

Finally, the terms in (B.12) containing H1 (r) are reduced to the form r

−2n−1

  2 2 2 2 (x +y )xP(2n−4) (x, y) − 2nxQ(2n−2) (x, y)H1 (r) + (x +y )B(2n−3) (x, y) H1 (r)

= r−2n−1 Q(2n−1) (x, y)H1 (r).

(B.15)

The right-hand sides of (B.13), (B.14), and (B.15) sum to the sought expression (B.9). The action of the operator Ln+1 u =

∂ Ln u ∂y

is reduced to the required form by the same method as (B.11), applying the fourth formula in (B.5). The representation (B.9) thus holds for every 1st order spatial operator applied to the operator Ln . This completes the proof for n + 1 = 2m + 1. Expression (B.10) is derived from (B.9) by the same procedure as in the previous stage. Q.E.D. Equality (B.7) directly implies the validity of the bound (B.3) for even n. Formula (B.8) for odd n also leads to (B.3) using relationship (B.6). The study has been supported by the Russian Foundation for Basic Research (grant 0.6-01-00293). REFERENCES 1. L. V. Dorodnitsyn, “Acoustic properties of continuous and discrete fluid-dynamic models,” Prikladnaya Matematika i Informatika, No. 6, 39 – 62 (2000). 2. L. V. Dorodnitsyn, “Nonreflecting boundary conditions for fluid-dynamic systems,” Zh. Vychisl. Matem. i Mat. Fiziki, 42, No. 4, 522 – 549 (2002). 3. L. V. Dorodnitsyn, “Artificial boundary conditions for simulation of subsonic gas flows,” Zh. Vychisl. Matem. i Mat. Fiziki, 45, No. 7, 1251 – 1278 (2005). 4. V. S. Vladimirov, Equations of Mathematical Physics [in Russian], Nauka, Moscow (1985). 5. A. T. Fedorchenko, “Reflection of a plane sound wave from a permeable surface in the presence of a normal gas flow,” Akust. Zh., 35, No. 5, 951 – 953 (1989). 6. M. B. Giles, “Nonreflecting boundary conditions for Euler equation calculations,” AIAA J., 28, No. 12, 2050 – 2058 (1990). 7. F. Q. Hu, “On absorbing boundary conditions for linearized Euler equations by a perfectly matched layer,” J. Comp. Phys., 129, No. 1, 201 – 219 (1996). 8. L. V. Dorodnitsyn, Nonreflecting Boundary Conditions: From a Conception to Algorithms [in Russian], Preprint, MAKS Press, Moscow (2002). 9. D. I. Blokhintsev, Acoustics of a Nonhomogeneous Moving Medium [in Russian], Nauka, Moscow (1981). 10. B. Engquist and A. Majda, “Absorbing boundary conditions for the numerical simulation of waves,” Math. Comp., 31, No. 139, 629 – 651 (1977). 11. L. N. Trefethen and L. Halpern, “Well-posedness of one-way wave equations and absorbing boundary conditions,” Math. Comp., 47, No. 176, L 421 – 435 (1986). 12. M. A. Il’gamov, “On nonreflecting conditions on boundaries of a simulation region,” in: Dynamics of Shells in a Flow [in Russian], Proceedings of a Seminar, No. 18, KFTI KF AN SSSR, Kazan’ (1985), pp. 4 – 76. 13. P. Luchini and R. Tognaccini, “Direction-adaptive nonreflecting boundary conditions,” J. Comp. Phys., 128, No. 1, 121 – 133 (1996). 14. R. L. Higdon, “Initial boundary-value problems for linear hyperbolic systems,” SIAM Rev., 28, No. 2, 177 – 217 (1986). 15. A. N. Tikhonov and A. A. Samarskii, Equations of Mathematical Physics [in Russian], Izd. MGU, Moscow (1999). 16. A. Bayliss and E. Turkel, “Radiation boundary conditions for wave-like equations,” Comm. Pure Appl. Math., 33, No. 6, 707 – 725 (1980). 17. T. Hagstrom, S. I. Hariharan, and D. Thompson, “High-order radiation boundary conditions for the convective wave equation in exterior domains,” SIAM J. Sci. Comp., 25, No. 3, 1088 – 1101 (2004).