Category COMPUTATIONAL AEROACOUSTICS

Nowadays, acoustic treatment is invariably used on the inside surface of commercial aircraft jet engines for fan noise reduction. Acoustic treatment panels or liners, when properly toned, are extremely effective noise suppressors. Because of structural integrity requirements, the acoustic liners are usually of the Helmholtz resonator type. Figure 10.1 shows a single-layer acoustic liner. The liner consists of a face sheet with holes at a regular pattern. Underneath the face sheet are cavities. The cavities control the frequency range within which damping is most effective. The dominant damping mechanism of this type of liner is attributed to the dissipation arising from vortex shedding at the mouths of the holes (see Melling, 1973; Tam and Kurbatskii, 2000; Tam et al., 2001). The vortex-shedding process converts acoustic energy into the rotational kinetic energy of the shed vortices. These vortices are subsequently dissipated by molecular viscosity. In the absence of vortex shedding, which is the case at a low sound pressure level, the principal damping mechanism is associated with viscous dissipation of the oscillatory flow at the mouths of the Helmholtz resonators. The flow and acoustic fields around the Helmholtz resonators of acoustic liners are exceedingly complicated, especially when there is a mean flow over the liner.

For engineering purposes, a gross macroscopic description of the effects of the acoustic liners on incident acoustic waves is definitely preferred over a more demanding microscopic description of the actual phenomenon. In the aeroacoustic community, it is an accepted practice to characterize the macroscopic properties of an acoustically treated surface by a single quantity Z called the impedance. Impedance is defined as the ratio of the acoustic pressure p to the acoustic velocity component normal to the treated surface vn (positive when pointing into the surface). That is,

P = ZVn. (10.1)

Impedance is a complex quantity, Z = R – iX (e-lmt time dependence is assumed). The use of a complex quantity is needed to account for the damping and phase shift imparted on the sound waves by the acoustically treated surface. The acoustic resis­tance R and the acoustic reactance X are generally frequency dependent. They also vary with the intensity of the incident sound waves and the adjacent mean flow veloc­ity. These quantities are usually measured empirically, although some semiempirical formulas are available for their estimates. It has been found experimentally that R is positive and does not vary much with frequency. On the other hand, X may be

Figure 10.1. A single-layer Helmholtz resonator type acoustic liner.

both positive or negative depending on the frequency. For many acoustic liners, the dependence of X on angular frequency ы can be represented adequately by a simple analytical expression of the following form (Tam and Auriault, 1996):

X _ X-1 + X

—- ———— + Xi ы.

p a0 ы

where X-1 and X1 are two parameters, p is the gas density, and a0 is the speed of sound. Figure 10.2 shows an accurate representation of the data of Motsinger and Kraft (1991) by the two-parameter formula. The values of X-1 and X1 obtained by mean-least-square fit are -13.48 and 0.0739, respectively, where ы is measured in kiloradians per second.

The impedance boundary condition (10.1) is basically a boundary condition established for frequency-domain analysis. As it is, it cannot be used in time – domain computation. For time marching computation, a suitable equivalent of the impedance boundary condition in the time domain is needed. One significant advan­tage of time-domain methods over frequency-domain methods is that broadband noise problems can be handled relatively easily, almost without extra effort. For

broadband noise problems, frequency-domain methods are computationally inten­sive and laborious. In addition, if the problem is nonlinear, frequency-domain meth­ods would be exceedingly difficult to carry out.

Time-domain problems can be solved only if they are well posed and stable. In the presence of a mean flow, standard formulation of the impedance boundary condition is known to support spurious unstable solution of the Kelvin-Helmholtz type. It turns out that, for acoustic liner problems, such instability is rather weak. It can be suppressed by the imposition of artificial selective damping at the liner surface.

In many internal aeroacoustics problems, the mean flow is driven by a pressure gradient. For example, maintaining the nozzle flow in Figure 9.17 would require a pressure difference imposed at x1 and x2. Now, if it is required to impose at x = x1 the boundary condition p = p1 and at x = x2 the boundary condition p = p2(p is the mean flow pressure), then it would be difficult to impose additional incoming wave or outgoing wave or absorbing boundary conditions at x = x1 and x = x2. To circumvent this pressure gradient boundary condition problem, it is recom­mended to perform a pressure gradient transformation to absorb the pressure gra­dient into the governing equations. In this way, it will remove the need for enforc­ing a static pressure boundary condition and thus allow the imposition of inflow/ outflow or absorbing boundary conditions at the two ends of the computational domain.

Figure 9.17. Nozzle flow problem re­quiring the imposition of a static pres­sure difference at x1 and x2.

Consider the nozzle flow problem of Figure 9.17. Mathematically, the pressure gradient field may be represented by

The pressure gradient is

d = в.

ax

Now, the following transformation of variables will move the pressure gradient boundary condition into the differential equations. Let

v = V, p = p’, p = P + —.

Substitution of transformation (9.72) into the Navier-Stokes equations, the gov­erning equations for (V, p’, —) in Cartesian tensor notation are found to be

dp + d(P'[12] j) = 0 91 d x ;

9 v – ‘ 9 v­

—- + v■ —- dt 1 9 x j

+v і d—+vi5 i1e+y (p+—)^=0

The velocity boundary conditions are unchanged by the transformation. The mean pressure boundary conditions for — are

x = x1; — = 0; x = x2, — = 0.

EXERCISES

Figure 9.18. Lines of constant phase associated with plane waves in the computational plane.

Consider plane wave solutions of the form,

u u V = Re V gKax+fty-mt) p p_

(2)

where Re{} is the real part of. Substitution of (2) into (1) gives

(m – MTOa)[(« – Mma}2 – a2 – в2] = 0.

(4)

The mode shape is given by the eigenfunction of Eq. (3). Upon replacing a and в by Eq. (7), the vorticity wave solution may be written as

The total wave number of vorticity wave in Eq. (8) is equal to <w/MOsin в. For convenience, let us set the total wave number to 2n or ю = 2n MOsin в. This is the same as assigning the length scale L to be equal to the wavelength.

Use a 200 x 200 mesh as the computation domain. Choose a suitable size Ax and Ay. Propagate a plane vorticity wave across the computation domain with в = 60° and sv = 10—5. You may add artificial selective damping if needed. Compare your numerical result with the exact solution. Find the maximum error over the entire computation domain over a period of oscillation.

9.2. This is a continuation of Problem 9.1. Now suppose one is interested in acoustic wave propagation. The dispersion relation is

(ю — MOa)2 — a2 — в2 = 0

or

ю = MO a + (a2 + в2)2 (1)

(the negative sign of the square root gives a second acoustic mode).

Let V and ф be the velocity and direction of propagation. Then, from the group velocity, it is easy to find

da = M° + (a2+в2 )1/2 = V cos ф (2)

If the direction of propagation ф is specified (see Figure 9.19), then the solution of Eqs. (1) to (3) gives

V = Mm cos ф + (1 — MJ sin2 ф)2 (4)

oj(V cos ф — MTO)

1 + MTO (V cos ф — MTO)

юУ sin ф

1 + MTO (V cos ф — MTO)

The eigenfunction for the acoustic wave is given by the solution of Eq. (3) of Problem (9.1). It may be written in the following form:

(7)

where V, a, and в are given by Eqs. (4), (5), and (6), and sa is the amplitude of the acoustic wave. It is easy to show that the total wave number is equal to ш /[1 + MTO(Vcosф – MTO)]. Thus, by taking ш = 2n[1 + MTO(Vcos<p – MTO)], the acoustic wavelength is chosen as the length scale.

Use a 200 x 200 mesh as the computation domain. Choose a suitable mesh size, propagate a plane acoustic wave train across the computation domain with <p = 30° and ca = 10-5. Compare your numerical results with the exact solution. Find the error distribution across the computation domain.

9.3. Many flow and acoustic problems are quasi-two-dimensional. This happens when the flow spreads out very slowly. A good example is the case of a jet. Locally, it makes sense to assume that the jet flow is parallel, that is, the velocity is parallel and independent of the axial coordinate x.

Consider the scattering of an incident plane acoustic wave by a jet. Let the mean axial velocity and the density of the jet be u(y, z) and p(y, z), locally parallel flow assumed. p is a constant. The linearized Euler equations are

d и _9 и d и d и 1 9 p

+ и + v + w + = 0

dt dX dу dz p dx

9 v _9v 19p

+ и + = 0

91 dx p у

d w _9 w 1 9 p

+ и + = 0

dt dx p z

dp dp

+ U + kp 91 dx

Let the incident plane acoustic wave be

i – (cos ©x+sin © sin фу+sin © cos фz—aJt)

p = ea0

where (©, ф) are the angular coordinates of the direction of propagation with respect to a spherical polar coordinate system centered at the jet exit, and a0 is the ambient sound speed.

Since the coefficients of Eq. (1) are independent of x, the x dependence of the scattered field must be the same as the incident wave. By factoring out exp[i(Q/a0)(cos©)x], the governing equations for the scattered sound field has spa­tial dependence on y and z only. Outside the jet (see Figure 9.20), where й = 0, the equations reduced to

dй iQ cos ©

+ p = 0

dt a0 P0

w + 1 3p = 0

dt P0 dy

д w 1 d p „

+ = 0 d t p0 d z

д v d w

й + — + = 0.

dy dz

It is easy to show that (3) forms a dispersive wave system. The asymptotic solution of a dispersive wave system is not suitable for constructing radiation boundary condition.

Develop a set of PML absorbing boundary condition for Eq. (3) to be used in the boundary region of a computation domain as shown in the figure. The PML is to absorb the scattered waves. Show that your PML equations are stable.

9.4. This problem is intended to test the discrete formulation and computation of sound transmitted through a surface of discontinuity. Here, the surface of discon­tinuity is formed by the interface of two fluids with different densities and sound speeds as shown in Figure 9.21. An incident acoustic wave at an angle of incidence в impinges on the interface. Part of the wave is transmitted and part of it is reflected. For computation purposes, we will use the following length, velocity, time, pressure,

Figure 9.21. Schematic diagram showing the incident, the reflected, and the trans­mitted waves.

and density scales. Subscripts 1 and 2 indicate fluids 1 and 2. length scale = L

velocity scale = a1 (sound speed in region 1) L

time scale = —

a

density scale = p 1 (density of ‘uid in region 1) pressure scale = p1a1

P1 Y1

The governing equations for small-amplitude disturbances in fluid 1 are

дР1 + ( дЩк + = 0

дt дх ду ‘

The governing equations for small-amplitude disturbances in fluid 2 are

Pi = P2’ vi = v2-

Now, consider a plane wave at an incident angle 0 and frequency ы given below:

where Re is the real part of. Determine the intensity and direction of the transmitted and reflected waves for the two cases with 0 = 20° and 65°. The frequency and other parameters are ы = 0.7, a = 0.694, and X = 1. Plot contours of p at intervals of 0.25 at the beginning of a cycle.

The exact solution is given in Section 1.4. At 0 = 65°, there is a total internal reflection. The transmitted wave amplitude is complex indicating that there is no transmitted wave radiated to the far field.

(Hint: In this problem, although both p and v are continuous at the vortex sheet their y derivatives are, however, discontinuous. Because of the discontinuous derivatives, one may use the ghost point method to enforce the continuity of p and v.)

9.5. For the purpose of computing sound radiation from a conical surface with a half cone angle 5 as shown in Figure 9.22, a set of natural coordinates to use is (f, n, Ф), where f and n are related to the cylindrical coordinates (r, х, ф) by

f = x/ cos 5, n = r – x tan 5.

(f, n) are, in fact, the body-fitted oblique Cartesian coordinates in the ф = constant plane. By assuming the solution to have a еітф dependence, show that, when еітф is

Figure 9.22. Oblique Cartesian coordi­nates (f, n) and a PML on a conical surface.

Figure 9.23. Profile of damping function a(z).

A subscript m is used to denote the solution associated with. (um, vm, wm) are the velocity components in the cylindrical coordinates. For computing sound radiation, it is proposed to use PML as boundary condition. The PML encloses the computational domain as shown in Figure 9.22.

Show that an appropriate set of PML equation is

where a (z) is a damping function with a profile as shown in Figure 9.23. The distri­bution of this function in the computational plane is shown in Figure 9.24.

Implement the PML equations to test their effectiveness. You may use an oscil­lating point monopole located at the center line of the conical surface as the source of sound.

In CAA, one often has to compute acoustic reflections or propagation over surfaces that are made up of different materials. At the junction of two materials, such as shown in Figure 9.14, the incident sound waves encounter a line discontinuity. Computationally, the grid size is, inevitably, finite. Because of the finite size mesh, experience indicates that spurious reflected waves are usually produced. To minimize spurious reflected waves, it is necessary to use a very fine computational mesh. This is sometimes a very costly remedy. It would be ideal to have a method that minimizes such spurious reflected waves without having to refine the computational mesh.

In jet engines, acoustic liners are, invariably, installed on the inside surface of the inlet and exhausted ducts, (see Motsinger and Kraft, 1991; Eversman, 1991). Acoustic liners are the most effective device for the suppression of engine fan noise at the present time. The junction between two acoustical liners with differ­ent impedance are examples of surface discontinuity. To facilitate the installation and maintenance of acoustic liners, hard wall splices are introduced. These splices are installed between two pieces of acoustic liners. Their presence creates additional surface discontinuities. To minimize the scattering of acoustic waves by hard wall splices, engine manufacturers resort to the use of very narrow splices. This poses a real dilemma for duct mode computation. In the absence of narrow splices, the use of a computational mesh with 7 to 8 mesh points per wave length in conjunction with a large-stencil CAA algorithm will have enough resolution for the computation. However, with narrow splices, the splice width is sometimes less than the size of a mesh spacing. This is as shown in Figure 9.15. In this case, the hard wall splice is invisible to the computation. Traditional wisdom will call for the use of much finer mesh. A factor of 10 or more mesh size reduction is often deemed necessary. This is to allow 8 to 10 mesh points to lie on the hard wall splice. This is, undoubtedly, very expensive. Not only the number of grid points has to increase substantially, but also the size of the time step used in the computation has to be significantly reduced. It
would, therefore, be extremely useful to have a method that can perform the nar­row splice-scattering computation using the original coarse grid without mesh size reduction.

In this section, a method capable of computing acoustic wave propagation and reflection in the presence of surface discontinuities as shown in Figures (9.14) and (9.15) using a large-size computational mesh is discussed. The method is known as the “wave number truncation method” (see Tam and Ju, 2009).

Now consider the wall boundary condition shown in Figure 9.14. On the hard wall side of the junction, the rigid wall boundary condition is

v = 0, (9.61)

where v is the velocity component normal to the wall. Let к be the elastic constant of the elastic wall and Z be the vertical displacement of the surface (positive in the direction of outward-pointing normal), then the elastic wall boundary condition is

(9.62)

Figure 9.15. A narrow splice with width smaller than a grid spacing.

where p is the pressure and pref is the reference pressure for zero displacement. By differentiating Eq. (9.62) with respect to t, it is found that
dZ = v = – Kdp d t V dt
(9.63)
Boundary conditions (9.61) and (9.63) may be combined into a single boundary condition by means of the unit step function, H(x). Let the wall junction be at x = 0, then the unified boundary condition at the wall is v = – кдР H (x). (9.64) dt It is useful to consider boundary condition (9.64) in wave number space. The Fourier transform of boundary condition (9.64) is
f H(a — k) — (k) dk. 2nj d t — X
V(a)	(9.65)
The right side of Eq. (9.65) is a convolution integral. The Fourier transform of the unit step function is
H (a) = ——^—. (9.66) 2na Thus, H(a) involves all wave numbers, including wave numbers in the short and ultrashort wave number range of the computational scheme (see Chapter 2, Section 2.2, and Figure 2.1). The convolution integral of Eq. (9.65) would, therefore, generate wave numbers outside the long wave range. In other words, because H(a) contains spurious waves with respect to the computational scheme, the surface discontinuity would scatter off spurious waves in an acoustic wave reflection computation. A way to minimize the generation of spurious waves is to remove all the wave numbers in the short and ultrashort wave number range from H (a); i. e., wave number higher than the cutoff. Thus, a modified boundary condition is to replace H(a) by H(a) [H(a + ac) — H(a — ac)] in Eq. (9.65), where ac is the cutoff wave number (see Figure 2.1). In physical space, this is tantamount to replacing H(x) in Eq. (9.64) by H(x), where H(x) is the inverse Fourier transform of H(a) [H(a + ac) — H(a — ac)]; i. e.,
11 = – +~ Si(acx). 2 П
(9.67)
In Eq. (9.67), the Sl(z) function is defined as (see Abramowitz and Stegun, 1964 [Chapter 5, Section 5.2]),
z

Figure 9.16 shows H(x) as a function of acx. Now, it is proposed to use the following modified boundary condition instead of the exact boundary condition (9.64):

v =-к pH(x). (9.69)

at

That is, H (x) is replaced by H (x). The effectiveness and accuracy of the wave number truncation method and its extension (narrow scatterer) will be demonstrated by several examples in the last section of Chapter 10.

In three dimensions, the linearized Euler equations may be written in a matrix form

as

d u 9 u „ 9 u

+ A + В

dt dx dy

where

A straightforward application of the method of Section 9.5.1 yields the following PML equations:

u + (ay + az )ql + ayazq2 u + (az + ax )ql + azaxq2

U + (ax + ay )ql + axayQ2

(ax + ay + az)U + (ayaz + azax + axay )4l + axayaz02-

The auxiliary variables ql and q2 are given by

Notice that the absorption coefficients ax, a and az are arbitrary functions of x, y, and z, respectively. Again, as in the two-dimensional case, the auxiliary variables, q1 and q2 need to be computed only in the PML regions.

The PML absorbing boundary condition is especially useful in a ducted envi­ronment. It is very effective for simulating a long-duct termination. Figure 9.13 shows a computational domain inside a circular duct with rigid walls. In cylindrical coordinates, (x, r, ф), the linearized Euler equations may be written as

where

The corresponding PML equation is

where

In applying PML to a ducted domain, the wall boundary condition for the PML is the same as that for the Euler domain. For instance, for the rigid wall circular duct problem of Figure 9.13, the same rigid wall boundary condition, namely, v = 0, at r = D/2 is to be used in the PML regions.

In implementing PML, the damping coefficients, say, ctx(x), cty(y), are often taken as smooth functions. A good practice is to set these functions equal to zero at the interface with the Euler domain. It is also a good practice to let ctx(x) increase smoothly to a constant level in the main part of the PML. At the termination of the PML, the computed variables are usually exponentially small so that reflection is not a concern. However, it is a good practice to impose a standard radiation boundary condition that is of minimal cost to the computation.

Perfectly matching between the solutions of the PML equations and the linearized Euler equations is one of the most important properties of a PML. It is, however, not obvious from the derivation that the governing PML equation actually possesses this unusual property. Here, it is intended to confirm perfect matching. For this purpose, consider the interface between the Euler computational domain and the vertical PML on the right side of Figure 9.9. If will now be assumed that the computational domain is large compared with the size of the outgoing disturbances, so that the interface may be regarded as infinite in extent in the y direction as shown in Fig­ure 9.10. For convenience, a Cartesian coordinate system centered at the interface will be used. In the Euler domain, the linearized Euler equation (9.35), when written out in full, becomes

dp dp du dv

— + M— I——— 1–

dt dx dx dy

9 u 9 u 9 p

—– + M—– + —

dt dx dx

9 v 9 v 9 p

—– + M—– + —

dt dx dy

9 P d P du dv

— + M— I——— 1–

dt dx dx dy

To allow the most general incident disturbances to impinge on the interface, one may choose to apply Fourier transform in y and Laplace transform in t to Eq. (9.44). This reduces Eq. (9.44) to a system of ordinary differential equations in x in the form (~ over a variable is used to indicate a Fourier-Laplace-transformed variable) as follows:

. _ . dp du _______

— iojp + M——I—-—+ ifiv =
dx dx

. _ , dU dp)

—io)u + M—— + —— = dx dx

dv)

—iov + M— + їв p = dx

dp) du)

—io p + M— + — + ipv =
dx dx

In Eq. (9.45), в is the Fourier transform variable and o is the Laplace trans­form variable. Eq. (9.45) is a fourth-order differential system in x. Four linearly

interface

independent solutions can easily be found. They are given as follows:

(a) Entropy waves

й = v = p = 0, /5 = A(p, oi)ei{a/M)x. (9.46)

The amplitude function A (в, o) is arbitrary.

(b) Vorticity waves

p = p = 0, й = Б(в, ш)еі{ш/М)х, v = – щ B(P, oj)ei(a/M)x

The amplitude function Б (в, o) is arbitrary.

(c) Acoustic waves

(9.48)

The amplitude function С(в, o) is arbitrary and

-oM + [o2 – в2 (1 – M2)]2
1 – M2

The branch cuts of the square root function in the o plane are shown in Figure 9.11. Eq. (9.48) gives the acoustic waves that propagate or spread to the right. The left propagating acoustic wave solution is obtained by replacing the square root function by its negative value. Note that, for acoustic waves that propagate to the far field, X+ must be real. This requires o2 > в2 (1 – M2). The components for which o2 < в2 (1 – M2) are nonpropagating acoustic near field. They decay exponentially with increasing x.

The most general disturbances that exit the Euler domain on the right side to enter the PML are represented by a combination of entropy, vorticity, and acoustic

Figure 9.11. Branch cuts for the square root function [o2 — в2 (1 — M2 )]1/2 in the complex o plane.

waves. It turns out that the PML equation supports a similar set of wave solutions. The PML equations, Eq. (9.42), when written out in full, are

dp dp d u 9

~dt + M dx + ax + ay(v + U-q-) + U-p + 1 — м2

9 u „9 u 9 p orM

—– + M——- + + uxu + t

dt dx 9 X 1 — M2

9 v 9 v 9 p 9 q4 „x

+ M + + oX —+ axv + – 2 v — 0

dt 9- 9y x dy 1 — M2

9 p d p du dv 9 q3 °xM

+ M + + + U – – + <Jxp + – 2 (u + ^^) = 0

dt dx dx 9 y 9 y 1 — M2

The Fourier-Laplace transform of Eq. (9.49) in y and t is

dp du v иx

шр + M + + ївv v — иx + – 2 (p + Mtt) — 0

dx dx xio 1 — M2

du dp ux

—lou + M + + – 2 (u + Mp) — 0

dx dx 1 — M2

. ,dv „ в ux

—lov + +Ів p – U- op + T—Mv = 0

, dp du в ux _ ___

—lop + M + + ifiv — ux v + – 2 (p + Mu) = 0.

d – d – x o 1 — M2

It should be pointed out that at the interface – = 0, the y and t dependence of the transmitted waves must be the same as those of the incident waves. This requires

that в and ш in Eq. (9.50) be the same as those in Eq. (9.45). Eq. (9.50) is a fourth – order differential system in x. The four linearly independent solutions, which can be verified easily, are as follows (a) Entropy waves

	x M
й = i? = p = 0,	(9.51)

where А(в, ш) is arbitrary. By comparing Eq. (9.46) and Eq. (9.51), it is clear by setting А (в, ш) = А (в, ш) that the entropy wave solution in the PML will match perfectly the entropy wave solution of the Euler domain at the interface x = 0. There is no reflected wave. In addition, Eq. (9.51) indicates that any entropy waves transmitted into the PML will be damped spatially with a damping rate of ax/[M(1 – M2)] in the x direction.

(b) Vorticity waves

By comparing Eq. (9.47) and Eq. (9.52), it is obvious that by setting B(fi, ш) = В(в, ш), the vorticity waves in the Euler domain are completely trans­mitted into the PML without reflection at the interface x = 0. In the PML, the spatial damping rate for the vorticity waves is the same as that of the entropy waves.

(c) Acoustic waves

(9.53)

where С(в, ш) is arbitrary and X = X+ + m(l—xM2)[ш2 — в2(1 — M2)]2.

By comparing Eq. (9.48) and Eq. (9.53), it is possible to conclude that there is perfect matching between the acoustic waves incident on the interface at x = 0 and the transmitted acoustic waves in the PML. Perfect matching is achieved by setting С(в, ш) = С(в, ш). This is true not only for propagating acoustic waves for which ш2 > в2(1 – M2), but it is also true for nonpropagating near-field components with ш2 < в2(1 – M2). The spatial damping rate in the PML is equal to ш(1—xM)2 [ш2 — в2 (1 — M2 )]2 in the x direction.

It was mentioned previously that the first version of PML equations were unsta­ble. Now, it is easy to show that PML equation Eq. (9.49) does not support unstable solutions. To prove this assertion, one is required to show that there is no unstable acoustic wave solution in the PML. The simplest way to complete the proof is to take the Fourier transform of Eq. (9.50) in x. Eq. (9.50) is a homogeneous system. It is straightforward to find the condition for nontrivial solution is given by setting the determinant of the coefficient matrix to zero. This gives, after some algebraic simplifications,

в2

1 — M2’

Figure 9.12. Mapping of the upper half of the ю plane on to the F plane. where

ю2 / Мю 2 (ю + iax )2 1 — M2

(9.55)

F (ю)

In Eq. (9.55), a is the Fourier transform variable for x. Eq. (9.54) is the dispersion relation for the acoustic waves supported by the PML equations. To prove that there is no unstable acoustic wave solution, it is sufficient to show that dispersion relation (9.54) has no roots in the upper half ю plane for arbitrary a and в. This can be done by mapping the upper half ю plane into the complex F plane. Figure 9.12 gives the image in the F plane. Notice that there is no value in the upper half ю plane that makes F real and positive. Therefore, Eq. (9.54) cannot be satisfied because the right side of the equation is real and positive. Since both the entropy and vorticity waves are damped waves, PML Eq. (9.49) would only have damped wave solutions.

The PML equation for the linearized Euler equation (linearized over a uniform mean flow in the x direction) may be derived following three basic steps. To make the discussion of each of these steps as simple as possible, two-dimensional problems in Cartesian coordinates will first be considered. Let L be the length scale, a0 (sound speed of the uniform mean flow) be the velocity scale, L/a0 be the time scale, p0 (gas density of mean flow) be the density scale, and p0a0 be the pressure scale. The dimensionless linearized Euler equation may be written in the following form:

where

M is the Mach number of the mean flow.

The first step is to make a change of variables from (x, y, t) to (X, y, t). The variables are related by

This transformation is necessary to make the PML equation stable. The lin­earized Euler equation now becomes

M d U d U 2 1 d U

I +- TA I -=■ + a + (1 – M2)1В— = 0­1 – m2 J dt dx dy

In Eq. (9.37), I is the identity matrix. Assuming, for the time being, a time dependence of e-i ®?, Eq. (9.37) may be rewritten as [11]

so that

dX + aB Yy

Eqs. (9.42) and (9.43) are the PML equations for rectangular computational domain as shown in Figure 9.9.

It is worthwhile to point out that the auxiliary variable q is needed only in the PML domains. This is because the spatial derivative dq/dx disappears for the PML equation when ay = 0 and, similarly, 9q/dy disappears when ax = 0. Thus, 9q/dx is required only in a horizontal PML and 9q/9y is required only a vertical PML. In short, there is no need to compute or store q in the Euler computation domain.

Perfectly matched layer (PML) was invented by Berenger (1994,1996) as an absorb­ing boundary condition for computational electromagnetics. The idea is to enclose a computational domain by a PML. In the PML, a modified set of governing equations are used. The modified equations, as designed, have the unusual characteristics that, when an outgoing disturbance impinges on the interface between the computational domain and the absorbing layer, no wave is reflected back. In other words, all out­going disturbances are transmitted into the absorbing layer where they are damped out.

Hu (1996) was the first to apply PML successfully to aeroacoustic problems governed by the linearized Euler equations. In the beginning, the linearization was done over a uniform mean flow. Recently, Hu (2005) had extended his work to nonuniform mean flow. Tam et al. (1998) analyzed the original version of the PML equations. They pointed out that those equations supported spatially growing unsta­ble solutions. This is because, in the presence of a mean flow, the PML equations give rise to dispersive waves. In the original version of PML, a small band of these waves has phase velocity in the opposite direction to the group velocity. Spatial damping in the PML is associated with the direction of the phase velocity. As a result, these waves grow in amplitude as they propagate across the PML instead of being damped. Hu (2001,2004) has since resolved this instability problem. In this chapter, the more recent PML equations developed by Hu are presented and analyzed.

Let the (x, r, ф) components of the velocity field be denoted by (u, v, w). The r derivative of u at points near the cylindrical axis may be treated exactly the same as if u is a scalar. However, care must be exercised in forming the finite difference approximation for dr and fw.

In a cylindrical coordinate system, v is positive in the positive direction of r. Now, to form a finite difference quotient to approximate the r derivative of v for a point at r = Ar in the direction of ф1, a 7-point symmetric stencil will extend to the other side of the origin at r = 0 as shown in Figure 9.7. For the points beyond r = 0, positive v is in the opposite direction. So, for the purpose of forming a directional derivative, we must use – v for these points. Hence,

Similar reasoning leads to the formula for dr as follows:

9.4.2.5 The Values of v and w at r = 0

The cylindrical axis is a singular line of the cylindrical coordinates. All the other points have a v and w velocity component. However, at r = 0, there is no radial or any azimuthal velocity component. In the r – ф plane, there is only a velocity vector vr=0 at r = 0.

At r = 0, the cylindrical coordinate system is singular. It is not recommended to compute the velocity vr=0 by means of the governing equations of motion. A good approximation is to determine vr=0 by interpolation from values of v and w of the first few rings of points around the origin r = 0 (see Figure 9.8). For mesh

Figure 9.8. Determination of v at r = 0 by interpolation using values of (v, w) at mesh points in the first or first two rings.

points at angle ф = ф1 and r = nAr, the velocity components in the directions of the cylindrical coordiates and in the y and z direction are related by

vy = v cos ф1 – w sin ф1; vz = v sin ф1 + w cos ф1.

Now, the value of vy and vz at r = 0 may be taken, as a first approximation, to be equal to the average of the values of vy and vz at all the mesh points of the first ring. A more accurate value of these two quantities may be obtained by applying the multidimensional interpolation procedure of Chapter 13 to the values at the mesh points on the first two rings of the grid.

In using cylindrical coordinates for two – or three-dimensional computation, there are two basic problems in computing the solution at and near r = 0. The first problem is how to approximate d/dr by a finite difference quotient at mesh points close to r = 0. The second problem is how to calculate the solution at r = 0 where there is an apparent singularity. This section addresses these two issues in the general situation.

9.4.2.1 Directional Derivative

Let es be a unit vector in the і direction as shown in Figure 9.5. To approximate З/Зі, one may use a 7-point finite difference quotient on a line in the es direction as shown; i. e., 1

Figure 9.6. Approximating the r deriva­tive at ф = ф1, r = Ar as directional derivative by a 7-point stencil finite dif­ference quotient.

Here, ЭФ/Эз = es – VФ is the directional derivative of Ф in the es direction. j = 0 is the point at which the derivative is to be computed.

9.4.2.2 r Derivative of a Scalar Variable Near r = 0

The r derivative of a variable Ф at ф = фх and r = Ar (see Figure 9.6) may be regarded as a directional derivative. That is,

A finite difference approximation to the directional derivative of Eq. (9.31) may be formed in the same way as Eq. (9.30). Thus,

where the origin r = 0 is at j = -1. In this way, except for r = 0, there is no problem in computing the solution to the discretized Navier-Stokes equations on a cylindrical coordinate mesh.

9.4.2.3 Scalar Variables at r = 0

Because the Navier-Stokes equations in cylindrical coordinates have apparent sin­gularities at r = 0, the values of scalar flow variables cannot be advanced in time by means of the discretized form of the equations. To find the values of the scalar variables at r = 0, a simple way is to compute first the solution at all other mesh points in the computational domain. After this is done, the values of the flow variables at r = 0 may be found by using high-order multidimensional optimized interpolation based on the values of all the mesh points on the first three rings of the cylindri­cal mesh. Just an average of all the values of the variable on the first ring is often a good approximation. For this purpose, the multidimensional optimized interpolation method of Chapter 13 would be very useful.

There are many aeroacoustics problems for which the use of a cylindrical coordinate system is the most natural choice for computing the solution. A simple example is the case of computing the noise from a circular jet. Another example is the numerical calculation of the propagation of acoustic waves inside a circular duct. The governing equations are the Navier-Stokes equations. When written in cylindrical coordinates, some terms of these equations with 1 /r coefficient apparently could blow up at the axis. In other words, there is an apparent singularity (although not a true one) at the axis (r = 0) of the coordinate system. In this section, a way to treat such apparent singularity is discussed.

9.4.1 Linear Problem Involving a Single Azimuthal Fourier Component

Acoustic problems are often linear. If the problem is linear and involves only a single azimuthal Fourier component, say the nth mode, then usually it is advantageous to separate out the azimuthal dependence by assuming that all variables have the dependence of єтф. On factoring out єтф, the reduced set of equations effectively becomes two-dimensional. Let the mean flow be axisymmetric. Near the x-axis, the mean flow may be assumed to be locally uniform; i. e., u = U, v = w = 0, p = p, p = p, and T = T. The linearized equations of motion are

+u +pv ■ v =0 (912)

¥ + = -1 vP + v v2v (9’13)

+4X) + pv- v=*tv2 T (914)

p = pRT + p RT, (9.15)

where vt and kt are the turbulent kinematic viscosity and thermal conductivity, if the flow is turbulent; otherwise, the molecular values should be used.

It can easily be shown that when Eqs. (9.12) to (9.15) are written in cylindrical coordinates (r, ф, x), it is an eighth-order differential system in r. Four of the solutions are bounded at r = 0. The other four are unbounded. Here, interest is confined to the bounded solutions alone.

In general, the velocity vector may be represented by a scalar potential Ф and vector potential A, that is,

v = vФ + v x A. (9.16)

For Eqs. (9.12) to (9.15), the solutions for the scalar and vector potentials can be obtained separately. It turns out that the viscous solutions are related only to the

vector potential. For convenience, the vector potential solutions will be referred to as viscous solutions.

Substituting Eq. (9.16) into Eqs. (9.12) to (9.15), it is easy to find that the governing equations for the scalar and vector potentials are

and

The viscous solutions are given by p = p = T = 0, v = Vx A.

9.4.1.1 Scalar Potential Solutions

General solutions of the scalar potential, Ф, may be found by applying Fourier – Laplace transforms to x and t and expanding ф dependence in a Fourier series, that is,

This reduces Eqs. (9.17) to (9.20) to a fourth-order differential system in r. The two solutions bounded at r = 0 may be expressed in terms of Bessel functions of order n. The complete solutions when written out in full are as follows:

un = ikJn (l■±r)- Vn = drrJn (X±r)

pn = p D’ (to — uk) — vt (4 + k2)] Jn (^±r)

where

ikt (to — uk)/R — p — cvpT + icVpvt [(to — uk)/R]
kt [vt/R + iT/(to — uk)]

—cvp(to — uk)2
Rkt [vt/R + iT/(to — uk)]

9.4.1.2 Viscous Solutions

There are two sets of viscous solutions. The first set is found by letting

A = ф ex, p = p = T = 0, (9.24)

where ex is the unit vectors in the x direction. Again, by the use of Fourier-Laplace transforms and Fourier expansion, the solution that is bounded at r = 0 is

p n pn Tn un 0

On following these steps, the bounded solution, after some algebra, is found to be

pn = pn = Tn = 0

9.4.1.3 Analytic Continuation into the r < 0 Region

The general solution is a linear combination of the scalar and vector potential solu­tions. The Bessel functions of these solutions can be continued analytically into the nonphysical region r < 0. For positive r, the analytic continuation formula for integer-order Bessel function is

Jn (—£r) = (—1)nJn (£r). (9.98)

By means of Eq. (9.98) and the preceding general solutions, it is straightforward to establish that

Pn (—r, x, t) = ( —1)n Pn (r, x, t)

p. n (—r, x, t) = ( — 1)npn (r, x, t)

Tn (—r, x, t) = (—1)n Tn (r, x, t) un (—r, x, t) = (—1)n un (r, x, t)

Vn (—r, x, t) = ( — 1)n+1 Vn (r, x, t) Wn (—r, x, t) = ( — 1)n+1 Wn (r, x, t),

Figure 9.5. Approximating directional derivative in the es direction at і by a 7-point stencil finite difference quotient.

where pn, pn, etc. are the amplitude functions of the Fourier series expansions in ф.

Now, Eq. (9.29) may be used to extend the solution into the nonphysical negative r region to facilitate the computation of high-order large-stencil finite difference in r for points adjacent to the jet axis. That is, when a finite difference stencil extends into the negative r region, the values of the variables are found by Eq. (9.29). It is noted, however, that for n = 1,

lim v1 (r, x, t) = 0, lim w1 (r, x, t) = 0,

r^0 r^0

in general. For this reason, terms such as v1/r and w1/r cannot be computed at the jet axis r = 0. It is recommended that the values of v and w and all the other variables not be computed directly by the finite difference marching scheme at r = 0. Instead, the numerical solution at a new time level for all the other mesh points is first computed. The values of all the variables on mesh points in the region r < 0 is then calculated by analytical continuation. This leaves the values at r = 0 still to be determined. Here, it is suggested that they are to be found by high-order symmetric interpolation (see Chapter 11). Once the values of all physical variables at the axis r = 0 are found, the computation at the new time level is completed.