Category COMPUTATIONAL AEROACOUSTICS

There is a fundamental difference between space and time. Space extends from negative infinity to positive infinity, whereas time only increases in one direction. This inherent difference requires finite difference approximation for time derivatives to be different from spatial derivatives.

In general, there are two ways to form a discretized approximation of a time derivative. They are the single time step method and the multilevel time discretiza­tion method. In this chapter, the single time step method will be discussed briefly. The multilevel discretization method will be discussed in greater detail. One major difference between the single time step method and multilevel methods is that, for wave propagation problems, the latter, when properly implemented, would lead to dispersion-relation-preserving (DRP) schemes (see Tam and Webb, 1993). DRP property is very desirable for computing wave propagation problems. DRP schemes will be discussed in detail in the next two chapters.

In this section, the coefficients of several large optimized stencils are provided. The range of optimization in wave number space is aAx = n.

9- Point Stencil

П = 1.28, a0 = 0.0

al = —a-1 = 0.8301178834769906875382633360472085 a2 = — a—2 = —0.23175338776901819008451262109655756 a3 = —a—3 = 0.05287205020483696423592156502901203 a4 = — a—4 = —6.306814638366300019250697235282424E-3

11-Point Stencil

П = 1.45, a0 = 0.0

ax = —a—1 = 0.8691451973307874495001324546317289 a2 = — a—2 = —0.28182159562075193452800172953580692 a3 = —a—3 = 0.08707110821545964532454798738037454 a4 = —a—4 = —0.019510858728038347594594712914321635 a5 = —a—5 = 2.265620835298174792121178791209576E-3

13-Point Stencil

П = 1.63, a0 = 0.0

al = —a—1 = 0.897853870484230495572085208841816 a2 = —a—2 = —0.32269821467978701682118324993041657 a3 = —a—3 = 0.12096287073505875036103692011036235 a4 = —a—4 = —0.037989102193448210970017384251701696 a5 = —a—5 = 8.526107608989087816184470859634284E-3 a6 = —a—6 = —1.0033637668308847022803811005724351E-3

The a Ax versus aAx curves of these stencils are shown in Figure 2.11. Figure 2.12 shows the variation of the slope ^Ax) versus aAx. Figure 2.13 shows an enlarged portion of Figure 2.11 near the maxima of the curves.

0.995

EXERCISES

2.1. Dimensionless variables with respect to length scale = Ax, velocity scale = c (speed of sound), and timescale = Ax/a are to be used.

Find the exact solution of the simple convective wave equation

d U d U

+ = 0 d t dx

Note that the entire solution propagates to the right at a dimensionless speed of unity.

Discretize the x derivative by the 7-point optimized finite difference stencil. Solve the initial value problem computationally. You may use any time marching scheme with a small time step At. Demonstrate that the numerical solution separates into two pulses. One pulse propagates to the right at a speed nearly equal to unity. The other pulse, composed of high wave numbers, propagates to the left. Explain why the solution separates into two pulses. Find the speed of propagation of each pulse from the j(aAx) versus a Ax curve and compare with your computed results.

2.2. Solve the following initial value problem computationally on a grid with Ax = 1

d U d U

+ = 0. dt dx

u (x, 0) = e ax cos (px)

Use the standard sixth-order central difference scheme.

Take a = 0.035 and consider three cases:

1. д = 0.8

2. д = 1.6

3. д = 2.1

For each case, plot u (x, t) at t = 300 and compare with the exact solution. Explain your results.

2.3. Solve the initial value problem:

— + — = 0 at t = 0, u = H(x + 50) – H(x — 50), dt dx

computationally on a grid with Ax = 1 , where H (x) is the unit step function,

Plot u(x, t) at t = 300. Determine from your numerical solution

(a) the most negative wave velocity (waves propagating to the left)

(b) the wave number of the wave with zero group velocity Compare the wave speed and wave number of (a) and (b) with the values from the ^ curve.

Do the above using

(1) the standard sixth-order central difference scheme and

(2) the optimized scheme (7-point stencil optimized over —1.1 < a Ax < 1.1).

Near the boundary of a computation domain, a symmetric stencil cannot be used; backward difference stencils are needed. Figure 2.10 shows the various 7-point back­ward difference stencils. Unlike symmetric stencil, the wave number of a backward difference stencil is complex. As will be shown later, a stencil with a complex wave

Figure 2.10. 7-point backward differ­ence stencils for mesh points shown as black circles.

number would lead to either damping or amplification depending on the direction of wave propagation. For numerical stability, artificial damping must be added. The subject of artificial damping is discussed in Chapter 7. The following is a set of stencil coefficients that has been obtained through optimization taking into consideration both good resolution and minimal damping/amplification requirements. These back­ward difference stencils are also useful near a wall or surface of discontinuity.

Note: For the three sets of stencil coefficients aj*, a51, and aj2 subscript j = 0 is the stencil point where the derivative is to be approximated by a finite difference quotient. The first superscript is the number of mesh points ahead of the stencil point in the positive direction and the second superscript is the number of mesh points behind.

By now, it is clear that, if one wishes to resolve short waves using a fixed size mesh, it is necessary to use schemes with large stencils. However, there is limitation to this. Beyond the optimized 7-point stencil (n = 1.1) with resolved wave number range up to a Ax = 0.95 (seven mesh spacings per wavelength), it will take a very large increase in the stencil size to improve the resolution to five or four mesh spacings per wavelength. So the cost goes up quite fast.

Consider a very large stencil, say N = 7 or a 15-point stencil. The stencil coefficients can be determined in the same way as before by minimizing the inte­grated error E.

n

E = J a Ax — aAx|2 d (a Ax).

-n

For n = 1.8 (a good overall compromise), the coefficients are

a0 = 0.00000000000000000000000000000 + 0 a3 = —a—1 = 9.1942501110343045059277722885,0 — 1 a2 = —a—2 = —3.55829599268352687556676424010 — 1 a3 = —a—3 = 1.52515016084064924691049286790 — 1 a4 = — a—4 = —5.94630408297157726668285968990 — 2 a5 = —a—5 = 1.90107527095082986598491679880 — 2 a6 = —a—6 = —4.38086492973364818511370009070 — 3 a7 = — a—7 = 5.38961218686233846596929558780 — 4

For comparison purposes, the coefficients of the standard fourteenth-order central difference scheme are

a0 = 0.00000000000000000000000000000 + 0 ax = —a—1 = 8.750000000000000000000000000000 — 1 a2 = —a—2 = —2.91666666666666666666666666670 — 1 a3 = —a—3 = 9.72222222222222222222222222220 — 2 a4 = —a—4 = —2.65151515151515151515151515010 — 2 a5 = —a—5 = 5.3030303030303030303030303030300 — 3 a6 = —a—6 = —6.79875679875679875679875679870 — 4 a7 = —a—7 = 4.16250416250416250416250416250 — 5

The a Ax versus a Ax and the dal da versus a Ax curves are shown in Figures 2.7 and 2.8, respectively. The choice of n is based on the enlarged da/da curve. For very large stencils, these curves exhibit small-amplitude oscillations. The n = 1.8 case gives a good balance between the resolved band width and constancy of group velocity. This stencil can resolve waves as short as 3.5Ax.

To demonstrate the effectiveness of the use of a large computation stencil to resolve short waves, again consider the initial value problem of Eqs. (2.15) and (2.16). The initial disturbance is again taken to be a Gaussian

f (x) = he— ln2( b)2, (2.23)

Figure 2.7. aAx versus aAx curve for the optimized 15-point stencil (n = 1.8) and the standard fourteenth-order cen­tral difference scheme………………………………………….. , opti­mized scheme; , standard

fourteenth-order scheme.

where h is the initial maximum amplitude of the pulse and b is the half-width. The Fourier transform of (2.23) is

hb _ а2ь2

f = – г Є 4ln2.

2(n ln2) 2

Notice that the smaller the value of b, the sharper is the pulse waveform in physical space. Also, the corresponding wave number spectrum of the pulse is wider. This makes it more difficult to produce an accurate numerical solution. The numerical results correspond to h = 0.5 and b = 3,2, and 1 using the standard fourteenth-order

Figure 2.8. d(aAX) versus aAx for the

—– 15-point optimized scheme (n =

1.8) and the……… standard fourteenth-

order central difference scheme.

Figure 2.9. Comparisons between the computed and the exact solutions of the convective wave equation. t = 400, h = 0.5. Dispersion-relation-preserving scheme is the optimized scheme, exact solution, computed solution.

scheme and the optimized scheme (n = 1.8) are given in Figure 2.9. With b > 2, the optimized scheme gives good waveform. However, for an extremely narrow pulse, b = 1, none of these schemes can provide a solution free of spurious waves.

There is no question that a good finite difference scheme should have a Ax equal to a Ax over as wide a bandwidth as possible. However, from the standpoint of wave propagation, it is also important to ensure that the group velocity of the finite difference scheme is a good approximation of that of the original equation. It will be shown later that the group velocity of a finite difference scheme is determined by da /da. da /da, the slope of the a Ax — a Ax curve, should be equal to 1 if the scheme is to reproduce the same group velocity or speed of propagation of the original partial differential equation.

Figure 2.2 shows the da/da curve of the optimized 7-point stencil scheme as a function of a Ax. It is seen that for 0.4 < aAx < 1.15, da/da is larger than 1. This will lead to numerical dispersion, because some wave components would propagate faster than others and faster than the wave speed of the original partial differential equation. Around aAx = 0.8 to 1.0, the slope is greater than 1.02. Thus, the wave component in this wave number range will propagate at a speed 2 percent faster (assuming time is computed exactly). That is to say, after propagating over 100 mesh spacings, the group of fast waves has reached the 102 mesh point. This is too much dispersion for long-range propagation.

One way to reduce numerical dispersion is to reduce the range of optimization. Instead of the range of aAx = 0 to n/2, one may use the range of 0 to n. That is, the integrated error is

n

E = j |aAx — aAx|2 d (aAx).

0

By setting n to a specified set of values, the coefficients ax, a2, and a3 of the 7-point central difference stencil may be found as before. Upon examining the corresponding da/da curve, it is recommended that the case n = 1.1 has the best overall wave

Figure 2.3. aAx versus a Ax.————- ,

optimized fourth-order scheme;——- ,

standard sixth-order scheme;.. , stan­dard fourth-order scheme; , stan­

dard second-order scheme.

propagation characteristics. The stencil coefficients are

a0 = 0 a1 = – a-1 = 0.77088238051822552

a2 = – a-2 = -0.166705904414580469 a3 = – a-3 = 0.02084314277031176.

Figure 2.3 shows the a Ax versus a Ax relations for the standard sixth-order, fourth-order, and second-order central difference scheme and the optimized scheme (n = 1.1). The values of aAx at approximately the point of deviation from the a Ax – aAx straight-line relation for these schemes, are

where X is the wavelength. Thus, a high-order finite difference scheme can better resolve small-scale features.

Figure 2.4 is a comparison of the curves. Figure 2.5 is the same plot but with an enlarged vertical scale. The increase in group velocity near aAx = 0.7 is around 0.3 percent. This means that the scheme has very small numerical dispersion. The resolved bandwidth for | dal da -1.01 < 0.003 is about 15 percent wider than that of the standard sixth-order scheme. This improvement is achieved at no cost because both schemes have a 7-point stencil.

example. Consider the initial value problem governed by the one-dimensional convective wave equation. Dimensionless variables with Ax as the length scale,

Ax/c (c = sound speed) as the timescale will be used. The mathematical problem is

The exact solution is

u = f (x —1).

Figure 2.5. Enlarged d(fAX) versus aAx

curves.—— , optimized scheme and the

…… , standard sixth-order central differ­ence scheme.

In other words, the solution consists of the initial disturbance propagating to the right at a dimensionless speed of 1 without any distortion of the waveform. In particular, if the initial disturbance is in the form of a Gaussian function,

with a half-width of three mesh spacings, then the exact solution is

(2.18)

On using a (2N + 1)-point stencil finite difference approximation to d/dx in Eq. (2.15), the original initial value problem is converted to the following system of ordinary differential equations:

t = 0 ul = f (l ).

l is the spatial index.

For the Gaussian pulse, the initial condition is

The system of ordinary differential equations (2.19) and initial condition (2.21) can be integrated in time numerically by using the Runge-Kutta or a multistep method. The solutions using the standard central second-order (N = 1), fourth-order (N = 2), sixth-order (N = 3), and the optimized 7-point (N = 3, n = 1.1) finite difference schemes at t = 400 are shown in Figure 2.6. The exact solution in the form of a Gaussian pulse is shown as the dotted curves. The second-order solution is in the form of an oscillatory wave train. There is no resemblance to the exact solution. The numerical solution is totally dispersed. The fourth-order solution is better. The sixth – order scheme is even better, but there are still some dispersive waves trailing behind the main pulse. The wavelength of the trailing wave is about 7, which corresponds to a wave number of 0.9, which is beyond the resolved range of the sixth-order scheme. The optimized scheme gives very acceptable numerical results.

To understand why the different finite difference scheme performs in this way, it is to be noted that the Fourier transform of the initial disturbance is

It is easy to verify that a significant fraction of the initial wave spectrum lies in the short-wave range (a Ax > ac Ax) of the a Ax versus a Ax curve of the standard second – and fourth-order finite difference schemes. (Note: Ax = 1 in this problem.) These wave components propagate at a slower speed (smaller than the group velocity). This causes the pulse solution to disperse in space and exhibits a wavy tail.

Another way to understand the dispersion effect is to superimpose the f (a) curve on the dal da versus a curve. The standard sixth-order scheme maintains a wave speed nearly equal to 1 (da/da = 1.0) up to a Ax = 0.85. Beyond this value the wave speed decreases. The optimized scheme maintains an accurate wave speed up to a Ax = 1.0. Now, for the given initial condition, there is a small quantity of wave energy around a Ax = 0.9. This component, having a slower wave speed, appears as trailing waves in the solution of the standard sixth-order scheme. The wavelength of the trailing wave is X = (2n|a) ~ 7.0. This is consistent with the numerical result shown in Figure 2.6. For the optimized scheme, there are trailing waves with wave number ~ 1.0 or X ~ 6.3, but the amplitude is much smaller and is, therefore, not so easily detectable. Basically, none of these schemes can resolve waves with wavelengths less than six mesh spacings really well.

Transform of f (x + X) = e>aX f (a).

Now, by applying the Fourier transform to Eq. (2.6) and making use of the Derivative and Shifting theorems, it is found that

where

a on the left of (2.8) is the wave number of the partial derivative. By comparing the left and right sides of this equation, it becomes clear that a on the right side is effectively the wave number of the finite difference scheme. a Ax is a periodic function of a Ax with a period of 2n. To ensure that the Fourier transform of the finite difference scheme is a good approximation of that of the partial derivative over the wave number range of, say, waves with wavelengths longer than four Ax or a Ax < (n/2), it is required that aj be chosen to minimize the integrated error, E, over this wave number range, where

2

E = j |aAx — aAx[3]d (a Ax). (2.10)

— 2

In order that a is real, the coefficients aj must be antisymmetric; i. e.,

a0 = 0, a—j = —aj. (2.11)

On substituting a from Eq. (2.9) into Eq. (2.10) and taking Eq. (2.11) into account, E may be rewritten as,

Eq. (2.13) provides N equations for the N coefficients aj, j = 1,2,…, N. Now consider the case of N = 3 or a 7-point stencil.

There are three coefficients a1, a2, and a3 to be determined. It is possible to combine the truncated Taylor series method and the Fourier transform optimization method.

Figure 2.1. a Ax versus a Ax for the 7- point fourth-order optimized central dif­ference scheme.

Let Eq. (2.14) be required to be accurate to O(Ax4). This can be done by expanding the right side of Eq. (2.14) as a Taylor series and then equating terms of order Ax and (Ax)3. This gives

9 4a1

20 – T 2 a1

-15 + У’

On substitution into Eq. (2.12), the integrated error E is a function of a1 alone. The value a1 may be found by solving

— — 0. d ax

This gives the following values

a0 — 0 a1 — —a-1 — 0.79926643

a2 — — a—2 — —0.18941314 a3 — — a—3 — 0.02651995.

The relationship between a Ax and a Ax over the interval 0 to n for the optimized 7-point stencil is shown graphically in Figure 2.1. For a Ax up to about 1.2, the curve is nearly the same as the straight line a Ax — a Ax. Thus, the finite difference scheme can provide reasonably good approximation for a Ax < 1.2 or X > 5.2Ax, where X is the wavelength.

For a Ax greater than 1.6, the a (a) Ax curve deviates increasingly from the straight-line relationship. Thus, the wave propagation characteristics of the short waves of the finite difference scheme are very different from those of the partial differential equation. There is an in-depth discussion of short waves later.

Figure 2.2. d(ffX) versus aAx for N = 3, fourth-order scheme optimized over the range of —n /2 to n/2.

In this chapter, how to form finite difference approximations to partial derivatives of the spatial coordinates is considered. The standard approach assumes that the mesh size goes to zero, i. e., Ax ^ 0, in formulating finite difference approximations. However, in real applications, Ax is never equal to zero. It is, therefore, useful to realize that it is possible to develop an accurate approximation to d/dx by a finite difference quotient with finite Ax. Many finite difference approximations are based on truncated Taylor series expansions. Here, a very different approach is introduced. It will be shown that finite difference approximation may be formulated in wavenumber space. There are advantages in using a wave number approach. They will be elaborated throughout this book.

2.1 Truncated Taylor Series Method

A standard way to form finite difference approximations to d/dx is to use Taylor series truncation. Let the x-axis be divided into a regular grid with spacing Ax. Index I (integer) will be used to denote the 1th grid point. On applying Taylor series expansion as Ax ^ 0, it is easy to find (xt = I Ax) as follows:

9 u 1 d2u 2 1

“’+■ = u(x’ + Ax) = u‘ + jїї),Ax + 2 j a?) Ax + 3!

j du) 1 j d 2u) 2 1

"’-■ = u(xt – Ax) = ,,t – – tAx + 2 ^ – j!

Thus,

„ j 9 u ) 1 j 9 3u ) 3

“<+1-u’-1 = 2Ax Vx, + 3 ,Ax +

On neglecting higher-order terms, a second-order central difference approximation is obtained,

(2.2)

This is a 3-point second-order stencil. By keeping a larger number of terms in (2.1) before truncation, it is easy to derive the following 5-point fourth-order stencil and the 7-point sixth-order stencil as follows:

From the order of the truncation error, one expects that the 7-point stencil is probably a more accurate approximation than the 5-point stencil, which is, in turn, more accurate than the 3-point stencil. However, how much better, and in what sense it is better, are not clear. These are the glaring weaknesses of the Taylor series truncation method.

2.2 Optimized Finite Difference Approximation in Wave Number Space

Suppose a central (2N + 1)-point stencil is used to approximate 9/9x; i. e.,

At one’s disposal are the stencil coefficients a; j = – N to N. These coefficients will now be chosen such that the Fourier transform of the right side of Eq. (2.5) is a good approximation of that of the left side, for arbitrary function f. Fourier transform is defined only for functions of a continuous variable. For this purpose, Eq. (2.5), which relates the values of a set of points spaced Ax apart, is assumed to hold for any similar set of points Ax apart along the x-axis. The generalized form of Eq. (2.5), applicable to the continuous variable x, is

Eq. (2.5) is a special case of Eq. (2.6). By setting x = lAx in Eq. (2.6), Eq. (2.5) is recovered.

The Fourier transform of a function F(x) and its inverse are related by

(Ю

F (a) = E [ f (x)e-iaxdx 2n J

— TO TO

F (x) =/ F (27)

— TO

The following useful theorems concerning Fourier transform are proven in Appendix A.

2.2.1 Derivative Theorem

Transform of 8

d x

It turns out that the discretized interface problem consisting of finite difference equations (1.61) to (1.64) and boundary conditions (1.65) and (1.66) supported nontrivial homogeneous solutions. These are eigensolutions. They will be referred to as boundary modes. These solutions are bounded spatially, but they may grow or decay in time. A numerical solution is stable only if all boundary modes decay in time. If there is one growing boundary mode, then, in time, this mode will dominate the entire solution leading to numerical instability.

To find the boundary modes, let them be of the form,

pfl = AeiaiAx+ie1mAy, m = -1, 0,1, 2,… (1.80)

p® = BeiaiAx-ie2mAy, m = 1, 0, – 1, – 2,… (1.81)

where a, the wave number, is real and arbitrary. вi and в2 may be complex. For spatially bounded boundary modes, it is required that

Im(e1) > 0; Im(e2) > 0. (1.82)

IfIm(e1 )= 0,thenRe(e1 )> 0 and, similarly, if Im(e2 )= 0,thenRe(e2 )> 0 to satisfy the outgoing wave condition.

Substitution of Eq. (1.80) into Eq. (1.61) yields

Similarly, substitution of Eq. (1.81) into Eq. (1.63) leads to

The corresponding v{pm and vf^ may easily be found by solving Eqs. (1.62) and (1.64). They are found to be

,2(1) = Л 8ш(в1Ду> еіаІДх+і0,тДу vl, m = A -(1Ь Є

шр(1> Ду

д(2> ЕрШ(в2Ду> ^іаІДх-ів2тДу

Vl, m = B -(2Ь Є ‘

шр(2> Ду

Upon imposing the dynamic and kinematic boundary condition (1.65) and (1.66) on solutions (1.80), (1.81), (1.85), and (1.86), it is straightforward to obtain

A – B = 0

For nontrivial solutions, the determinant, D, of the coefficient matrix of the above 2 x 2 linear system for A and B must be equal to zero; i. e.,

The eigenvalues, «Дх/й(1>, of the boundary modes are the roots of Eq. (1.87) for a given аДх.

To determine the eigenvalues numerically, the following grid search method may be used. To start the search, the complex «Дх/а(1> plane is divided into a grid. The intent is to calculate the complex value of D at each grid point. Once the values of D on the grid are known, a contour plot is used to locate the curves Re(D) = 0 and Im(D) = 0 (where Re() and Im() denote the real and imaginary part). The intersections of these curves give D = 0 or the locations of the boundary mode frequencies. To improve accuracy, one may use a refined grid or use a Newton or similar iteration method.

To calculate the value of D at a grid point in the complex «Дх/а(1> plane for a prescribed density and mesh size ratio and аДх, the values of в 1Ду and в2Ду are first computed by solving Eqs. (1.83) and (1.84). All these values are then substituted into Eq. (1.87) to compute D. Figure 1.7 shows the contours Re(D) = 0 and Im(D) = 0 for the case p(1>/p(2> = 2.0 , Дх = Ду, and аДх = n/10. These two curves intercept each other on the real «Дх/а(1> axis at «Дх/а(1> = (0.536, 0.0). Since «Дх/а(1> is real, this root corresponds to a neutral mode. With the addition of artificial selective damping to the finite difference scheme, which is a standard procedure and will be

Figure 1.7. Location of the roots of D in the complex «Ax/5[1] [2] plane. p(1)/p(2) =

2.0, Ax = Ay, and aA = n/10.———— ,

Re(D) = 0.0;…….. , Im(D) = 0.0.

discussed in a later chapter, the mode will be damped. Thus, the finite difference scheme is not subjected to numerical instability arising from boundary modes.

EXERCISES

where l, m are the spatial indices in the x and y directions (I = 0,1,2,…, N; m = 0, 1,2,…, M). The boundary conditions for the finite difference equations are

at l = 0, l = N, m = 0, m = M, Ul = 0.

l, m

Now, suppose the simple wave equation is discretized by using a second-order central difference approximation. This gives

d2ue c2

~dF = (Xx? (u’+1 _ 1 +111-1

where l is the spatial index.

Consider an incident wave of the form

ut (t) = Ає-і(кІАх+юі

where к Ax = cos-1 (1 – ^2-). Note: the value of к, the wave number, is determined by the requirement that the incident wave satisfies the governing finite difference equation.

Find the reflected wave of the discretized system and show that the reflected wave has the same amplitude as the incident wave. Notice that the wave number к and, hence, the phase velocity, «/к, and group velocity, йш/йк, of the waves supported by the finite difference approximation are not the same as those of the original simple wave equation. Is the approximation good at high or low frequency?

1.3. Acoustic waves in a medium at rest are governed by the linearized Euler and energy equations. Let (u, v) be the velocity components in the x and y directions and p0 and p0 be the static density and pressure. The governing equations are

du 1 d p

+ = 0 d t p0 dx

A +1 dJP = 0

d t P0 dy

where y is the ratio of specific heats. Upon eliminating u and v, it is easy to find that the governing equation for p is

д 2 p = c2 = m

dt2 c дx2 + дy2 ’ c p0

Consider plane acoustic waves propagating at an angle в with respect to the y-axis incident on a wall lying on the x-axis. The incident wave is given by

p = Ae – — (X sin в +y cos в +ct)

p incident Ae

The boundary condition at the wall is

У = о, v = 0.

By the у-momentum equation, the wall boundary condition for pressure p is

At у = 0, — = 0.

У ду

In order to satisfy the wall boundary condition, a reflected wave is generated at the wall. This reflected wave cancels the pressure gradient of the incident wave in the у direction at the wall, namely,

It is easy to find that the solution of Eq. (A) that satisfies Eq. (C) is,

„ _ aJc (-x sin в +у cos в-ct)

preflected

By formulas (B) and (D), it is straightforward to show that the following well-known properties of acoustic wave reflection are true.

(i) The angle of reflection is equal to the angle of incidence.

(ii) The reflected wave has the same intensity as the incident wave.

(iii)

There is a pressure doubling at the wall, that is,

where the overbar denotes a time average.

Suppose the problem is solved by finite difference approximation with a mesh size of Ax and Ду. On using a second-order central difference approximation, Eqs. (A) and (C) become

where l, m are the spatial indices in the x and у directions and the wall is located at m = 0.

Find the incident wave with frequency ю and angle of incidence в. Find the ref­lected wave so that boundary condition (F) is satisfied. Determine whether the acoustic wave reflection properties (i), (ii), and (iii) are satisfied. Would you be able to find the incident and reflected waves if a fourth-order central difference scheme is used?

Consider a plane acoustic wave of angular frequency ы incident on the interface at an angle of incidence в as shown in Figure 1.4. The appropriate solution of Eqs. (1.38) to (1.40) may be written in the following form:

(1.47)

where A is the amplitude and Re{} is the real part of. The reflected wave in region (1) has a form similar to Eq. (1.47), which may be written as

(1.48)

where R is the amplitude of the reflected wave. The transmitted wave in region (2) must have the same dependence on x and t as the incidence wave. Let

(1.49)

By substituting Eq. (1.49) into Eqs. (1.41) to (1.43), it is easy to find after some simple elimination,

On solving Eq. (1.50), the transmitted wave with an amplitude T may be written as

(2) (2) p(2)

9v(2) _ 1 9p(2)

dt p(2) 9 y

A rectangular mesh with mesh size Ax and Ay, as shown in Figure 1.5, will be used for the finite difference model. Let I and m be the mesh indices in the x and y directions. The fluid interface is at m = 0.

The time dependence of all variables is in the form e-lmt. This may be factored out from the problem. For simplicity, let

ui1m (t)=Re{ pp е-ш (1.60)

and similarly for all the other variables. Suppose the standard second-order central difference is used to approximate the spatial derivatives in Eqs. (1.56) to (1.59), it is straightforward to obtain

+ 1 (nd) – 2 f,(1) + nd) ^

+ (Ay)2 lm+1 2pl, m + pl, m-1)

Figure 1.5. Computation mesh.
(1.62)
+ — d(2) — 2 d(2) + d(2) + (Ay)2 Pl, m+1 2 pl, m + pl, m—1
= 0	(1.63)
(1.64)

Now, because a 3-point central difference stencil is used, Eqs. (1.61) and (1.62) should be applicable to values of m > 1. However, if it is insisted that they are to hold true at m = 0, just above the interface discontinuity, then a problem arises. These equations will include the nonphysical variable pf— 1. Similarly, Eqs. (1.63) and (1.64) should be applicable to m < -1. But if they are to be satisfied at m = 0, just below the interface discontinuity, then there is the problem of nonphysical variable pfl appearing in the equations. pf— 1 and p® will be referred to as ghost values and the corresponding row of points at m = — 1 and m = +1 are the ghost points. In a time – domain computation, the values pf ‘0 , v^O are given by the time marching scheme of Eqs. (1.61) and (1.62). Similarly, pfO and v(20 are given by those of Eqs. (1.63) and (1.64). Boundary conditions (1.44) and (1.45), however, require the continuity of pressure and vertical velocity component at the interface. This would not be possible in general. Here, it is important to recognize that, although the pressure is continuous at the interface, the pressure gradient is not. In fact, the pressure gradient should be discontinuous. The jump at the discontinuity is determined by boundary conditions (1.44) and (1.45). In the finite difference model, one may accept pf—1 and p® , the ghost values, as the variables controlling the pressure gradients. These ghost values

of pressure are to be chosen such that the discretized version of boundary conditions (1.44) and (1.45), namely,

f(1) – f ® p 1,0 — p,0

yd) — y(2)

v 1,0 — vl,0

are satisfied. Notice that there are two boundary conditions. They can be enforced by choosing the two ghost values properly.

The finite difference model of sound transmission as formulated, consisting of Eqs. (1.61) to (1.66), can be solved exactly analytically. The exact solution may be found as follows.

The incoming wave at an angle of incidence в and having a total wave number к may be written in the form

p(incoming) ^^/k(sinвДхІ+cosвДym)

cos(k sin в Ax) + cos(k cos вДу)
&
— 0. (1.68)
+
(Ax)2	(Ay)2

The cosine functions are symmetric with respect to к, so that if к is a positive real root, then – к is also a root. The positive real root is the total wave number of the incident wave. If & is complex, then к is also complex. It is to be determined by analytic continuation in the complex plane. The reflected wave has the same dependence on I but propagates in the positive m direction. It is straightforward to find that the total wave field in region (1), being the sum of the incident and the reflected waves, may be written as

p(l) Ae_ ‘к (sin eAxl+cos eAym) + Rg_’^(sin eAxl—cos eAym)

(m — 0,1, 2,…).

R is the reflected wave amplitude. By substituting Eq. (1.69) into Eq. (1.62), the vertical velocity component is found, as follows:

5Іп(кcoseAy) г Ag—^(sineAxl+coseAym) + r ^Ік(— sineAxl+coseAym)] m 12 3

&p(1)Ay L ’ ’ ’ * ‘

(1.70)

The transmitted wave in region (2) must have the same dependence on I as the incident and reflected waves. Thus, let

f.(2) 7^ „_ік sin в Axl—ів Aym

pl, m —

(m — 0, _ 1, _ 2, _ 3,…).

The transmitted wave must satisfy (1.63). On substituting (1.71) into (1.63), it is found that the equation is satisfied if в is the root of

For the outgoing wave, в is the positive real root. If ш is complex, then в is obtained by analytic continuation.

The vertical velocity component of the transmitted wave is given by the substi­tution of (1.71) into (1.64). This yields

It should be obvious that the ghost value of pressure must have the same dependence on I as the incident, reflected, and transmitted waves. Hence, we may write

^(1) ^ ,,—iksineAxl

pl, —1 = p— 1Є

^(2) ^ „—iksineAxl

pl,1 = p1e.

There are four boundary conditions that need to be satisfied at m = 0. The first two are from Eqs. (1.61) and (1.63). By setting m = 0, these equations give

(1.76)

The other two boundary conditions are the dynamic and kinematic boundary con­ditions (1.65) and (1.66). The incident, reflected, and transmitted wave solutions as found consist of Eqs. (1.69), (1.70), (1.71), and (1.73). There are four constants in the solution; they are R, T , p—1 , and p 1 . By substituting the solution into the four boundary conditions, it is straightforward to find that the transmitted and reflected wave amplitudes are given by

T _ 2

A 1 + і р^ 1 sin (в Ay)

+ ур(2) / sin(k cos eAy)

R = T — A.

It is easy to show that, in the limit Ax, Ay ^ 0, k ^ ш/а(1), в ^ ш/а(2) [1 — (р(1)/р(2)) sin2 в]1/2, so that Eq. (1.78) is identical to the exact solution (1.54) of the continuous model.

Figure 1.6 shows the transmitted wave amplitude computed by Eq. (1.78) for the case р(1)/р(2) = 1.2 with &>Ax/a(1) = n/6 and n/4, Ax = Ay as a function of the angle of incidence. There is good agreement with the exact solution (1.54). For the given density ratio, total internal refraction occurs at в = 65.91°. For an

Figure 1.6. Transmitted wave amplitude as a function of the angle of inci­dence. , exact solution; ,

&>Ax/a(1) = n/6;….. , &>Ax/a(1) = n/4.

incident angle greater than this value, the incident wave is totally reflected at the interface and there is no transmitted wave. Tbecomes complex, and the transmitted wave becomes damped exponentially in the negative y direction. The cutoff angle as determined by Eq. (1.78) is found to be very close to the exact value.

How best to transform a boundary value problem governed by partial differential equations into a computation problem governed by difference equations is not always obvious. The task is even more difficult if the original problem involves a surface of discontinuity. There is a lack of discussion in the literature about how to model a discontinuity in the context of finite difference. The purpose of this section is to show how one such model maybe set up. At the same time, this model will demonstrate that the finite difference formulation of boundary value problems may support spurious boundary modes. These modes might have no counterpart in the original problem.

Figure 1.4. Schematic diagram showing (a) the incident, reflected, and transmitted sound waves at a fluid interface, (b) the slightly deformed fluid interface.

For static equilibrium, the pressure balance condition is

p(1) = p(2) or P(1)(fl(1)}2 = P(2)(fl(2)}2.