EFFICIENT CALCULATION OF THE THREE-DIMENSIONAL SOUND PRESSURE FIELD AROUND A NOISE BARRIER

Journal of Sound and Vibration (1996) 197(5), 547–571

EFFICIENT CALCULATION OF THE THREE-DIMENSIONAL SOUND PRESSURE FIELD AROUND A NOISE BARRIER D. D Ecole Nationale des Ponts et Chausse´es CERAM, 1 avenue Montaigne, 93167 Noisy Le Grand Cedex, FR (Received 4 November 1993, and in final form 25 March 1996) A numerical method is presented for calculating the sound pressure around a noise barrier of constant but arbitrary cross-section. To obtain the exact solution of this problem, the Helmholtz equation must be solved in the three-dimensional domain outside the barrier. It is shown in this article how to calculate this 3-D sound pressure from solutions of simpler problems defined on the two-dimensional domain outside a cross-section of the barrier. The numerical solution of a large three-dimensional problem is thus avoided and the efficiency of the calculation is considerably improved especially when a whole frequency spectrum is needed. By using the boundary element method to calculate the numerical solutions in the two-dimensional domains, examples are given of determinations of sound pressure fields created by a point source and by an incoherent line source. In this way the efficiency of barriers of different cross-sections can be compared by using the real sound pressure around them. From the frequency spectrum at a point, one can then calculate by Fourier transform the temporal variations of the sound pressure created by a noise source moving in a direction parallel to the barrier. 7 1996 Academic Press Limited

1. INTRODUCTION

A noise barrier is the most common protection used against noise propagation in outdoor acoustics. To be able to estimate the efficiency of a noise barrier of given geometry, one must calculate its excess attenuation, defined as the quotient of the sound pressure with the barrier over the sound pressure in free field. The calculation of excess attenuations of barriers is thus the main modelling problem and knowledge of acoustic fields around such obstacles is required. For a straight wall on a rigid ground, the acoustic field is well approximated by semi-analytical formulas derived from the Sommerfeld-MacDonald solution for a semi-infinite screen. Maekawa [1] used these results to propose a method which has long been used for noise barrier calculations. Pierce [2], Jonasson [3] and Tolstoy [4] improved the method by solving the diffraction problem for a two-dimensional angle or for a polygonal line. For more complex shapes, Keller’s geometrical theory of diffraction allows the determination of asymptotic behaviors for high frequencies in the shadow zone; see, for instance, the papers of Kurze and Anderson [5] and Kurze [6]. However more precise results for arbitrary shapes and boundary conditions can be derived from numerical solutions of the Helmholtz equation in the domain outside the noise barrier. They are 547 0022–460X/96/450547 + 25 $25.00/0

7 1996 Academic Press Limited

548

. 

usually obtained from the boundary element method for which one needs only to mesh the boundary of the fluid domain. Using such a method Seznec [7] and Hothersall et al. [8, 9] presented numerical results for two-dimensional diffraction problems. In their calculus the sound pressure is supposed to be created by a coherent line source, the sound pressure of which is given by p(r) = (i/4)H0 (kr),

(1)

where r is the vector between the sound source and the observation point, r is its norm, k is the wavenumber equal to 2pf/c, with f the frequency and c the sound speed. H0 is the Hankel function of the first kind. The Helmholtz equation is solved in a two-dimensional domain which represents the domain outside a cross-section of the wall. Unfortunately the experimental data and the real sound pressures usually consist of noise radiated from a three-dimensional point source, the sound pressure of which is given by p(r) = eikr/4pr

(2)

The two-dimensional models can give only a crude estimate of the real sound field created by such a point source. However, Daumas [10] and subsequent authors compared measured insertion losses for sound pressures created by point sources with results of two-dimensional calculations. They found a close agreement between the two values. So 2-D calculations give very useful results to predict the performances of 3-D barriers. One must notice, however, that the comparison is made on the amplitude of the sound pressure and no result is given for the phase. So the possibility of calculating interference effects is not clear. The second and more important remark is that these comparisons are made for situations where the source and the receiver are in the same plane perpendicular to the barrier. For reception points outside this plane, the 2-D calculations provide no way to estimate the pressure. So the solution of the true three-dimensional problem is desirable in order to estimate correctly the real attenuation everywhere outside the barrier and to compare it with experimental results usually obtained by using a rather point-like sound source like a loudspeaker. Furthermore, knowledge of the pressure in a straight line along the barrier allows one to calculate the attenuation for an incoherent line source, which is known to be a good description of road traffic noise. The three-dimensional sound pressure can be obtained by numerical methods like the boundary element method or the finite element method. Unfortunately these methods usually require meshing a two-dimensional surface for the boundary element method or a large part of the three-dimensional domain outside the barrier for the finite element method. So the problem to be solved becomes rapidly very large and only simulations for small walls of limited length and for low frequencies are usually done. Using the boundary element method Filippi [11] gave results for a noise source, assuming a finite wall with zero thickness and a wavelength of the same order as the dimensions of the wall. Kawai [12] and Antes [13] made calculations for finite length walls with pressure fields created by point sources. They discretized a wall of rectangular shape and made calculations for frequencies about 100 Hz for walls of some meters in height and width. For instance, for a real barrier with a length of 10 m and a height of 2 m, the discretization of the rectangular surface by a classical 3-D boundary element method with five nodes per wavelength needs approximately 500/l 2 unknowns where l is the wavelength. At 1000 Hz, l = 0·34 m and one must solve a full linear system with about 4 × 103 unknowns which leads to possible

3

   

549

but heavy computations. In comparison, the 2-D solution on a cross-section only requires 10/l 1 30 unknowns for the same frequency. To avoid such difficulties in 3-D calculations, it is proposed here to solve these problems by using a series of solutions of simpler two-dimensional problems in the frequency domain, for both real and complex frequencies. To solve the two-dimensional problems one can use any numerical or analytical method. For the examples described in this article results are presented which are obtained from a numerical solution based on the boundary element method on a cross-section of the barrier. The meshes are thus one-dimensional, the calculations are two-dimensional but the mathematical transformation presented will give the three-dimensional sound pressure. Furthermore for a broad band frequency analysis it is shown that the whole spectrum for a three-dimensional source can be deduced from the two-dimensional spectrum for real and imaginary frequencies. Then it is shown how this can be used to estimate the attenuation provided by noise barriers, especially in the case of an incoherent line source. Finally, by Fourier transformation, the time variations of the pressure are calculated for a source moving in a direction parallel to the barrier.

2. PROBLEM SPECIFICATION

It is supposed that the barrier has a constant cross-section, as for the example presented in Figure 1, and that it is infinite in length. In other words the geometry is unchanged by a translation in the z-direction. The sound pressure is supposed to be created by a single frequency point source in the fluid domain, with time behavior e−ivt. To calculate the pressure around the barrier one must solve the three-dimensional problem Dpdif + k 2pdif = 0

in V3 ,

1pdif /1n + 1pinc /1n = 0

pinc (r) = eik=r − rs =/4p=r − rs =,

on 1V3 ,

1pdif /1n − ikpdif = o(1/r),

(3)

where V3 is the exterior 3-D domain outside the barrier with boundary 1V3 . The sound pressure pinc and pdif are respectively the incident and the diffracted fields and the sum ptot = pinc + pdif is the total pressure. It is supposed that the boundary of the fluid domain 1V3 , consisting of the ground and the barrier, is rigid. The vector rs is the position of the sound source. The last relation is the radiation condition for outgoing waves at infinity.

Figure 1. Structure of constant cross-section.

. 

550

According to the limiting absorption principle (see reference [14]), one knows that the solution is also the limit when n:0+ of solutions of Dpdif + k 2(n)pdif = 0

in V3 ,

1pdif /1n + 1pinc /1n = 0

on 1V3

pinc (r) = eik(n)=r − rs =/4p=r − rs =,

(4)

with pdif $ L 2(V3 ), k(n) = zk 2 + in, n q 0 and Im(k(n)) q 0. To solve this problem by the boundary element method one must descretize the boundary or at least a large part of the wall. The discretization of the ground can be avoided if one uses a Green function with an image source in the ground to take account of the rigid boundary condition. Nevertheless the dimension of the discrete problem to be solved limits the practical solutions to walls of small lengths and to the low frequency domain. To avoid this complex computational scheme, one can start from the following relation in which the field of a three-dimensional point source in free space is calculated from the fields of two-dimensional sources by Fourier transform (see reference [15], formula 6·616): 2 + z2

eik(n)zr

=

4pzr 2 + z 2

1 2p

g

+a

−a

i H (rzk 2(n) − a 2) eiaz da. 4 0

(5)

Here 0 E arg (zk 2(n) − a 2) Q p and r = zx2 + y 2 q 0 is the radial distance in the (x, y) plane between the observation point and the source. For imaginary values of the variable one has, when u q 0, (i/4)H0 (iu) = (1/2p)K0 (u) where K0 is the modified Bessel function of zero order. The function (i/4)H0 (r zk 2(n) − a 2) is in fact the sound pressure of a line source with a possible extension to the imaginary axis. One can then extend this relation for point sources to solutions of diffraction problems and calculate the solutions of 3-D problems as Fourier transforms of solutions of 2-D problems. So, one can define a family of two-dimensional problems depending on a complex parameter m. One must calculate the diffracted qdif , incident qinc or total qtot sound pressures solutions of the following two-dimensional exterior problems with rigid boundary conditions on the ground, on the barrier and outgoing wave conditions at infinity: (D + m2)qdif = 0

in V2 ,

1qdif /1n + 1qinc /1n = 0

qinc (x, y, m) = (i/4)H0 (rm).

on 1V2 , (6)

V2 is the 2-D fluid domain outside a cross-section of the barrier and 1V2 is its boundary. The outgoing wave condition can be expressed as the existence of a function f with compact support such that qdif (x, m) =

g

V2

i H (m=x − z=)f(z) dz. 4 0

(7)

In the case of a real parameter m one can use the more familiar Sommerfeld radiation condition 1qdif /1n − imqdif = o(1/zr).

(8)

This problem has a unique solution for complex m such that Im(m) e 0 (see reference [14]).

3

   

551

Now the solution of the three-dimensional problem (4) can be recovered from the solutions of the previous two-dimensional problems by the formula p(x, y, z, k(n)) =

1 2p

g

+a

eiazq(x, y, zk 2(n) − a 2) da,

(9)

−a

where p and q denote the incident, diffracted or total sound pressure. The solution of the three-dimensional problem is thus the Fourier transform of two-dimensional solutions which are much easier to obtain. The solution of problem (3) can be found as p(x, y, z, k) = lim+ n:0

1 2p

g

+a

eiazq(x, y, zk 2(n) − a 2) da,

(10)

−a

Remark. When k Q =a= the solution q(x, y, zk2 − a 2) is exponentially decreasing because K0 (z) 0 zp/2z e−z

(11)

for large values of z. So the integral (9) is in fact exponentially decreasing in a when =a= q k and, generally, only a rather small range in the imaginary complex domain needs to be calculated. To estimate the gain in computation provided by formula (10), one notices that the boundary element solution of the two-dimensional problem requires only meshing a curve instead of a surface for the three-dimensional problem. Supposing for instance that the 2-D mesh is made of n nodes. For the same precision, the surface mesh of the 3-D problem needs a number of nodes proportional to n 2. If the computational time is proportional to the power three of the number of nodes, the times required are about n3 in 2-D and n 6 in 3-D. One must, however, make several 2-D calculations. To estimate the number of calculations needed, one notices that the 2-D solution oscillates approximately like eikr where r is the distance between the source and the receiver. Here rmax is the maximum distance of interest. To get five points per period one must sample the wavenumber with the spacing Dk = 2p/5rmax

(12)

Df = c/5rmax ,

(13)

and the frequency with

So the maximum number of calculations to get the 3-D pressure at a frequency f is, in the real domain, n = f/Df = 5frmax /c = 5rmax /l.

(14)

With approximately an equal number of points for the imaginary frequencies, one can estimate the number of frequency points by nf = 10rmax /l.

(15)

Taking again the example described in the introduction, one can estimate the number of calculations to be (500/l 2)3 in the 3-D case and (10/l)3 × 10rmax /l = 104rmax /l 4 in the 2-D case. In Table 1 these two values are compared for different frequencies, for rmax = 50 m.

. 

552

T 1 Comparison of the number of operations Frequency (Hz)

100

500

1000

5000

l (m) 2D case 3D case

3·40 4 × 103 8 × 104

0·68 2 × 106 109

0·34 4 × 107 8 × 1010

0·07 2 × 1010 1015

Although one needs to calculate for several frequencies in the 2-D case, the proposed method saves many calculations. The reduction is greater when the frequency is higher. In this example, at 1000 Hz, the 3-D calculation is about 2000 times longer than the evaluation of the 2-D series. This approach is still more interesting when a complete frequency spectrum over a band [0, fmax ] is needed. In the 3-D case one must calculate for each required frequency. With the proposed method only two-dimensional solutions over the band [0, fmax ] and in a limited part of the complex imaginary frequency axis are needed. The integral (9) has a numerical cost which is negligible compared to the solutions of problems (6). So the 3-D spectrum is calculated with a computational cost similar to the 2-D one.

3. NUMERICAL SOLUTION OF 2-D PROBLEMS

To solve the two-dimensional problem one can use any method. To illustrate the previous discussion the following examples will be calculated with the boundary element method. Let G k(y, x) be the Green function for a half-space with a rigid boundary condition: that is, Gk(y, x) = (i/4)H0 (kr) + (i/4)H0 (kr'),

(16)

where H0 is the Hankel function of first kind of order zero. The distances are defined by r = =y − x= and r' = =y − x'= where x is the source point, x' the image source with respect to the ground and y the observation point (see Figure 2). The total pressure qtot solution of equations (6) satisfies the Kirchoff relation on the boundary

c(y)qtot (y) −

g

1V

qtot (x)

1G k dx = qinc (y), 1nx

Figure 2. Definition of domains and normals.

(17)

3

   

553

with c(y) = u(y)/2p if y$1V ( if y is a smooth point). u(y) is the exterior angle at y and qinc (y) is the incident pressure created by sources in the exterior domain. For the special case in which y is one of the points on the boundary where the barrier intersects the ground surface one has in fact c(y) = u(y)/p. It is supposed that the normal at the boundary is directed into Vext . To simplify the expression the rigid boundary condition is used. Relation (17) can be transformed into the more regular expression 1 2

qtot (y) −

g

(qtot (x) − qtot (y))

1V

1G 0 dx − 1nx

g

0

qtot (x)

1V

1

1G k 1G 0 − dx = qinc (y); 1nx 1nx

(18)

where G0 is the static Green function given by G 0(y, x) = −(1/2p) log r − (1/2p) log r'.

(19)

To avoid possible problems associated with the singular frequencies for which equation (18) has more than one solution and which are the resonance frequencies of the interior domain, the Burton and Miller [16] formulation of the integral relation is used. This consists in adding to relation (18) its derivative multiplied by a complex parameter b with Im(b) $ 0. They showed that this new equation has a unique solution at every frequency. Lin [17] proved the validity of this approach. Transforming the derivative equation to improve its regularity one finally obtains qtot (y) −

g

(qtot (x) − qtot (y))

1V

−b

1G 0 dx − 1nx

g

0

qtot (x)

1V

g

(qtot (x) − qtot (y) − 9t qtot (y) · (x − y))

g

9t qtot (y) · nx

1V

−b

1V

= qinc (y) + b

1G 0 dx − b 1ny

g

1V

qtot (x)

0

1

1G k 1G 0 − dx 1nx 1nx 1 2G 0 dx 1nx 1ny

1

1 2G k 1 2G 0 − dx 1nx 1ny 1nx 1ny

1qinc (y), 1ny

(20)

where 9t means the tangential derivative. This relation involves weakly singular terms only. The boundary condition 1q/1n = 0 has been used to remove the corresponding terms in the equation. More details on the derivation of this relation can be found in reference [18]. In the following, the parameter b is chosen as b = i=k=/max(=k=2, 50);

(21)

thus it gives the classical value i/k at high frequency and a small value at low frequency. The discretization is made by quadratic three nodes elements with the interpolation r(j) = N1 (j)r1 + N2 (j)r2 + N3 (j)r3

(22)

for the geometry where r1 , r2 and r3 are the positions of the nodes and Ni (j) are the interpolation functions defined on [−1, 1]. One has a similar expression for the interpolation of the pressure. To solve equation (20) numerically the collocation method is used. To avoid possible problems with the non-uniqueness of the normal at the extremity nodes of the elements

. 

554

Figure 3. Nodes (Ni ) and collocation points (Ci ) inside an element.

the collocation points are chosen to be inside the element as indicated in Figure 3. In an element one defines three possible collocation points which have the co-ordinates −0·95, 0·05 and 0·95 in the j interval. To get a final square system one must associate a unique collocation point to each node. For the interior node only one choice is possible (C2 ). For an extremity node the collocation point is chosen among the two nearest collocation points. Two choices are possible because there is a collocation point in the two elements having this node respectively as first or third node. The final choice depends on the global numbering of the nodes. Which one is effectively chosen has very little influence on the final result. After the numerical evaluation of the integrals over each element one obtains the global system AQ = F

(23)

where

&'

q1 . Q= . . qn

and

F=

&

'

qinc (y1 ) + b 1qinc /1ny (y1 ) . . , . qinc (yn ) + b 1qinc /1ny (yn )

(24)

n is the total number of nodes, y1 , . . . , yn are the collocation points and q1 , . . . , qn the nodal values of the pressure. A results from the assembly of the elementary matrices of each element.

4. COMPARISON WITH KNOWN 3-D SOLUTIONS

To estimate the precision of the previous methods and their efficiency, the numerical results are compared in two cases for which the 3-D solution can be easily calculated. The main purpose is to study the influence of the discretization on the numerical errors. These discretization errors have two origins. A first part of the error comes from the numerical solution of the 2-D problems by the boundary element method. This study is classical and nothing new is provided here. Secondly the calculus is based on the evaluation of the 2-D solution at discrete frequencies, on finite intervals. To calculate the 3-D pressure for a frequency f, one must calculate the 2-D solution in the interval [0, f] and for the imaginary frequencies. So the precision of the result will depend on the interval chosen in the imaginary axis. With [0, i fmax ] denoting this interval, fmax must be chosen large enough to allow one to neglect the contribution of frequencies larger than fmax in the numerical evaluation of formula (10). Formula (11) can provide an estimate of this upper limit. To get a solution with an error of order e for points at distance larger than rmin one must have e−armin E e:

(25)

3

   

555

Figure 4. Point source in a half-space.

that is, (c/2prmin ) log (1/e) E fmax .

(26)

The solution is also strongly dependent on the number of frequency points chosen in each interval. Formula (14) can provide an upper bound on the number of points required in the real interval [0, f]. In the imaginary interval, the pressure is more regular and the number of points can be smaller. The first case is a point source over a rigid plane for which the analytical solution is derived from relation (5) and is given by eikzr2 + z2 4pzr 2 + z 2

+

eikzr'2 + z2 i = 4pzr'2 + z 2 8p

g

+a

(H0 (rzk 2 − a 2) + H0 (r'zk 2 − a 2)) eiaz da.

(27)

−a

r is the distance between the projections of the source and the reception points in the (x, y) plane. r' is the corresponding distance for the image source. The sound pressure is to be calculated at the three points defined in Figure 4. Their co-ordinates are given in Table 2. Point R1 is a rather difficult case, being very close to the sound source. In this example the two-dimensional sound pressure is given analytically by q(f) = (i/4)H0(2pfr/c) + (i/4)H0 (2pfr'/c),

(28)

with c = 343·5 m/s. This expression is evaluated for discrete frequencies values. To estimate the influence of the frequency discretization and the value of fmax , four series of calculations were made with the following values for the frequencies (all in Hz): (D1) f0 = 0·0001, f1 = 0·01, f2 = 0·1, f3 = 1, f4 = 5 and between f5 = 10 and fmax = f204 = 2000 in steps of 10; (D2) f0 = 0·0001, f1 = 0·01, f2 = 0·1, f3 = 1, f4 = 5 and between f5 = 10 and fmax = f403 = 2000 in steps of 5; (D3) f0 = 0·0001, f1 = 0·01, f2 = 0·1, f3 = 1, f4 = 5 and between f5 = 10 and fmax = f404 = 4000 in steps of 10; (D4) f0 = 0·001, f1 = 0·01, f2 = 0·1, T 2 Source and reception points Source R1 R2 R3

x (m)

y (m)

z (m)

0 0·1 10 10

0·8 0·8 1 1

0 0 0 10

. 

556

T 3 Frequency steps and errors associated with the three points Df (Hz) R1 R2 R3

e (fmax = 2000 Hz)

687 6·9 6·9

2·6 × 10 0 0

e (fmax = 4000 Hz)

−2

6·6 × 10−4 0 0

f3 = 1, and between f4 = 2·5 and fmax = f1603 = 4000 in steps of 2.5. For each discretization the pressures q(fn ) and q(i fn ) were calculated. Between the frequencies fn and fn + 1 the pressure was calculated by a quadratic interpolation on the values of q(fn − 1 ), q(fn ) and q(fn + 1 ), thus one could rebuild a continuous approximation of the exact function from these discrete frequency points. The value of Df calculated by formula (13) and the error e, obtained by taking equality in formula (26), are given in Table 3 for the three calculation points. The symmetry allows one to calculate i 4p

g

amax

(H0 (rzk 2 − a 2) + H0 (r'zk 2 − a 2)) cos (az) da,

(29)

0

2 where zamax − k 2 = 2pfmax /c. This integral is divided into four parts for a respectively in [0, 0·5k], [0·5k, k], [k, 1·5k] and [1·5k, amax ]. Each part is divided into sub-intervals and in each interval the numerical integration is made by a Gauss formula with seven points. Changes of variable can be used to improve the numerical integration near a = k. Table 4 gives the pressure at points R1, R2 and R3 for different frequencies and for the four frequency discretization schemes (D1), (D2), (D3) and (D4). The influence of the maximum frequency of calculation fmax is especially important for points close to the source point as one can see by examining the results for the point R1. For this point the decrease of the frequency step has no influence while the increase in fmax clearly improves the result.

T 4 Pressure in a half-space Frequency (Hz) R1 R1 R1 R1 R1

analytical numerical numerical numerical numerical

R2 R2 R2 R2 R2 R3 R3 R3 R3 R3

100

500

1000

2000

0·734 + 0·155i 0·725 + 0·155i 0·726 + 0·155i 0·733 + 0·155i 0·734 + 0·155i

0·461 + 0·673i 0·452 + 0·673i 0·452 + 0·673i 0·460 + 0·673i 0·461 + 0·673i

−0·228 + 0·726i −0·239 + 0·726i −0·239 + 0·726i −0·227 + 0·726i −0·227 + 0·726i

−0·717 − 0·350i −0·842 − 0·349i −0·842 − 0·349i −0·716 − 0·349i −0·716 − 0·349i

analytical (10−2) numerical (D1) (10−2) numerical (D2) (10−2) numerical (D3) (10−2) numerical (D4) (10−2)

1·431 − 0·626i 0·986 − 0·689i 1·330 − 0·654i 0·986 − 0·689i 1·420 − 0·630i

−0·558 − 1·041i −0·590 − 0·784i −0·589 − 1·010i −0·590 − 0·787i −0·563 − 1·040i

−0·096 + 0·162i −0·068 + 0·147i −0·096 + 0·163i −0·066 + 0·146i −0·096 + 0·163i

0·553 + 1·431i 0·584 + 1·100i 0·576 + 1·400i 0·587 + 1·100i 0·557 + 1·430i

analytical (10−2) numerical (D1) (10−2) numerical (D2) (10−2) numerical (D3) (10−2) numerical (D4) (10−2)

0·746 + 0·829i 0·741 + 0·534i 0·711 + 0·815i 0·741 + 0·534i 0·740 + 0·827i

−0·488 − 0·845i −0·530 − 0·566i −0·489 − 0·826i −0·532 − 0·565i −0·488 − 0·841i

−0·287 + 0·500i −0·154 + 0·483i −0·275 + 0·505i −0·154 + 0·484i −0·285 + 0·501i

0·256 + 0·461i 0·362 + 0·377i 0·268 + 0·442i 0·362 + 0·377i 0·257 + 0·460i

(D1) (D2) (D3) (D4)

3

   

557

Figure 5. Diffraction by a rigid cylinder.

The results for the points R2 and R3 show that density of the discrete frequencies is essential for points far from the source. On the contrary the value of fmax is not important, as expected. The results could be improved by using the fact that, for large distances, the pressure behaves like eikr and it would be possible to interpolate e−ikrp(r, k) instead of p(r, k). This could save many calculations in the 2-D solution especially for points at large distance for which, according to formula (13), the frequency step must be small. In the second case, the pressure of a point source diffracted by a rigid cylinder has been studied. A solution of this problem is given in the Appendix. It consists mainly of Fourier transforms of Bessel functions. This solution can be efficiently calculated with arbitrary precision by standard numerical routines. In the present case this formula has been calculated with high accuracy and it provides a reference solution to which the numerical results will be compared. So the global efficiency of the numerical method (boundary element solution followed by a numerical Fourier transform of the interpolated results) can be estimated by comparison with this reference solution. The cylinder has a radius of 1 m and is infinite in length. The calculations are for three different points which are defined in Figure 5. Their co-ordinates are given in Table 5. The source is on the x-axis so the solution is symmetrical with respect to the (x, z)-plane and is equal to the pressure diffracted by a rigid half-cylinder on a rigid plane. To estimate the influence of the frequency discretization two calculations were made with the following values for the frequencies (Hz): (D1) f0 = 0·001, f1 = 0·01, f2 = 0·1, f3 = 1, f4 = 5 and between f5 = 10 and fmax = f204 = 2000 in steps of 10; (D2) f0 = 0·0001, f1 = 0·01, f2 = 0·1, f3 = 1, f4 = 5 and between f5 = 10 and fmax = f403 = 2000 in steps of 5. The circular section of the half-cylinder was divided into N + 1 nodes of co-ordinates Nn = (cos (np/N), sin (np/N))

(30)

where 0 E n E N. N = 101 was taken to get approximately five nodes per wavelength at the maximum frequency. To calculate the numerical solutions of the two-dimensional problems a quadratic mesh was used with three node elements as described in Figure 3. T 5 Source and reception points around the cylinder Source R1 R2 R3

x (m)

y (m)

z (m)

2 4 10 10

0 0 0 1

0 0 0 10

. 

558

T 6 Pressure around a cylinder Frequency (Hz) −2

100

500

1000

2000

R1 appendix (10 ) R1 numerical (D1) (10−2) R2 numerical (D2) (10−2)

−2·621 − 1·187i −2·886 − 1·220i −2·618 − 1·203i

3·816 − 3·270i 3·740 − 3·279i 3·814 − 3·275i

0·912 − 4·533i 0·868 − 4·519i 0·908 − 4·533i

−2·736 − 1·890i −2·762 − 1·886i −2·767 − 1·887i

R2 appendix (10−2) R2 numerical (D1) (10−2) R2 numerical (D2) (10−2)

−0·197 + 0·522i −0·480 + 0·539i −0·228 + 0·521i

−1·071 − 0·902i −1·066 − 0·721i −1·081 − 0·876i

0·150 + 1·252i 0·172 + 1·097i 0·164 + 1·236i

−0·783 + 0·018i −0·722 + 0·264i −0·777 + 0·025i

R3 appendix (10−2) R3 numerical (D1) (10−2) R3 numerical (D2) (10−2)

0·209 − 0·526i −0·141 − 0·540i 0·173 − 0·532i

−0·448 − 0·809i −0·515 + 0·724i −0·463 − 0·660i −0·415 + 0·682i −0·451 − 0·797i −0·503 + 0·730i

−0·198 − 0·706i −0·270 − 0·639i −0·207 − 0·695i

Figure 6. Real part of the pressure at R1. —, Reference; ·,, calculated.

Figure 7. As Figure 6 but imaginary part of the pressure.

3

   

559

When the 2-D solutions have been found for these discrete frequencies, the 3-D solution is calculated numerically as in the half-plane case. In Table 6 the analytical and numerical results for four frequencies are compared. Both

Figure 8. As Figure 6 but at R2.

Figure 9. As Figure 7 but at R2.

Figure 10. As Figure 6 but at R3.

. 

560

Figure 11. As Figure 7 but at R3.

discretizations give an approximation of the solution. The discretization (D2) is clearly better than (D1), as expected. In Figures 6 and 7 the calculated solution is compared with the reference one for the real and imaginary components of the pressure at point R1 on the frequency band [0, 2000 Hz] with the frequency discretization (D2). Figures 8–11 show similar results for points R2 and R3. The comparison is quite good except perhaps at very low frequencies where some differences can be seen. The two preceding comparisons have thus proved the accuracy of the proposed method. The frequency discretization does not need to be very fine except perhaps for very low frequencies or for observation points far from the source point. For most practical cases one can calculate the pressure of the three-dimensional problem for frequencies as high as 2000 Hz without difficulty, which is not the case if a full three-dimensional problem were to be solved.

5. APPLICATIONS TO NOISE BARRIERS

5.1.   The excess attenuation is defined as the sound pressure level with a wall divided by the sound pressure in free field: At = 20 log10 (=pwall =/=pfree =).

(31)

This formula can be used directly for point sources or coherent line sources for which =pfree = = (1/4pr) or =pfree = = 14=H0 (kr)= respectively. Calculations for incoherent line sources parallel to the barrier will also be made. An incoherent line source is modelled as uncorrelated three-dimensional point sources located on a straight line parallel to the z-axis. The amplitude of such a source as function of the abscissa z can be represented as a random process with the cross-correlation function E(m(z)m*(z + u)) = n(z) d(u),

(32)

where n(z) e 0. In free field the pressure of the whole line at a point x = (x, y, 0) is given by

p(x) =

g

+a

−a

eikr(z) m(z) dz 4pr(z)

(33)

3

   

561

where r(z) = zx + y + z . The density of acoustic potential energy has the expectation 2

2

2

0 1

e=E

=

=p(x)=2 1 = 4rc 2 4rc 2

1 4rc2

g

+a

−a

g g +a

−a

+a

eikr(z) e−ikr(z') E(m(z)m*(z') dz dz' 4pr(z) 4pr(z')

−a

n(z) dz. [4pr(z)]2

(34)

(35)

If n is constant one can calculate the energy density for a line in free field by efree =

1 pn n = , 4rc 2 (4p)2r 64prc 2r

(36)

with r = zx 2 + y 2. In the presence of a wall, the pressure created at point x = (x, y, 0) by a unit source located in the plane z = z0 is, with respect to formula (10), P(x, z0 ) =

1 2p

g

+a

e−iaz0q(x, y, zk 2 − a 2) da.

(37)

−a

For a line source with the amplitude m(z) one has

g

p(x) =

+a

P(x, z)m(z) dz.

(38)

−a

The expectation of the density of acoustic potential energy is thus ewall =

1 4rc 2

g

+a

=P(x, z)=2n(z) dz.

(39)

−a

If the amplitude n is constant ewall =

n 4rc 2

g

=P(x, z)=2 dz

(40)

g

=q(x, y, zk 2 − a 2)=2 da.

(41)

+a

−a

Using Parseval’s relation yields ewall =

n 8prc 2

+a

−a

This integral is calculated by Gaussian quadratures as for the numerical evaluation of formula (29). Note that this formula needs only the knowledge of the 2-D solutions and can be efficiently calculated. The attenuation provided by the barrier is then given by At = 10 log10 (ewall /efree ).

(42)

5.2.     The previous results can now be used to compare the efficiency of noise barriers of different shapes. The attenuations for the four cross-sections of Figures 12–15 were

562

. 

Figure 12. Straight wall. Case A.

Figure 13. Left-bent top. Case B.

Figure 14. Right-bent top. Case C.

Figure 15. T-shape wall. Case D.

calculated with a sound source at 50 cm above the ground and 1·9 m to the right of the wall at point (2, 0·5). This source can be a point source, a coherent or an incoherent line source. The observation point is 2 m behind the wall and 50 cm above the ground at point (−2, 0·5). The two points are on the same cross-section (same z). The origin O is always at the bottom left of the wall. In the four cases, the thickness of the wall is e = 10 cm and the total height H is equal to 2 m. For case D the width of the top is 50 cm. For cases B and D the wall is bent from (3/4) H to H at an angle of 45°. The walls and the ground are assumed rigid.

3

   

Figure 16. Excess attenuation. Case A. —, Coherent line; ·····, point source; – – – , incoherent line.

Figure 17. As Figure 16 but Case B.

Figure 18. As Figure 16 but Case C.

563

564

. 

Figure 19. As Figure 16 but Case D.

Figures 16–19 show that a point source and a coherent line source give almost the same result. The attenuation is a complex curve especially for the cases B, C and D because of interference between the different waves diffracted by the barrier and reflected from the ground. The comparison between these curves to find the best shape is not easy. An incoherent line source gives, on the contrary, very different results with an attenuation increasing with frequency and a much lower average value. This result needs to be modelled three-dimensionally and cannot easily be deduced from the two-dimensional attenuation for a coherent line source. There is no interference phenomenon in this case. So the attenuation seems to be less dependent on a particular choice for the frequency or the position. This could provide a better estimate of the efficiency of the wall. In Figure 20 the attenuations for an incoherent line source for different forms of the top of the wall are compared. The differences are small but the T-form seems a little better by some decibels. In Figure 21 the results presented here are compared with the numerical results of reference [9]. The calculation is made for a source at point (−15, 0, 0) and a receiver at point (50, 0, 0). The wall has a T-shape with a width of 20 cm, the top of the T- is 1·5 m

Figure 20. Differences in attenuation for different forms of the top of the wall. —, A; ·····, B; – – – –, C; — — —, D.

3

   

565

Figure 21. Comparison with reference [9] results. —, Reference [9] results; ·····, coherent line results; – – – –, point source; — — —, incoherent line source.

wide and the global height is 2·72 m. The wall is symmetrical with respect to the y-axis. The comparison of the results shows a good agreement between the results of reference [9], calculated for a coherent line source at the octave center frequencies, and the calculations here for the 2-D coherent line source and the 3-D point source insertion losses. The insertion loss is defined as the quotient between the pressure without the wall over the pressure with the wall. The excess attenuation for an incoherent line source is also given, showing a lower attenuation at high frequencies.

6. MOVING SOURCES

The previous results were obtained in the frequency domain, supposing the source was harmonic in time. One can, however, use this to get results for non-stationary cases. More precisely, the pressure at a point created by a source moving in a direction parallel to the barrier can be calculated. 6.1.  The pressure is the solution of DP − (1/c 2)1 2P/1t 2 = −s(t) d(x − X(t))

(43)

if the source sends a signal with the amplitude s(t) at the position X(t) and time t. The source is supposed to move along the z-direction, so X(t) = (x0 , y0 , z(t)). With p denoting the Fourier transform of P in z and t, one obtains

P(x, y, z, t) =

1 (2p)2

−a

g g +a

p(x, y, a, v) =

g g +a

−a

+a

p(x, y, a, v) eiaz e−ivz da dv,

−a

+a

−a

P(x, y, z, t) e−iaz eivt dz dt.

(44)

. 

566

Taking the Fourier transform of equation (43) in z and t yields

g g +a

Dx,y p − a 2p + k 2p = −

−a

0g

+a

s(t) d(x − X(t)) e−iaz eivt dz dt

−a

1

+a

s(t) e−iaz(t) eivt dt d(x − x0 ) d(y − y0 )

=−

−a

= −S(a, v) d(x − x0 ) d(y − y0 ),

(45)

with S(a, v) =

g

+a

s(t) e−iaz(t) eivt dt.

(46)

−a

So p(x, y, a, v) = S(a, v)q(x, y, zk 2 − a 2)

(47)

where q is the previously defined solution of the two-dimensional problem. The pressure field is thus P(x, y, z, t) =

1 (2p)2

g g +a

−a

+a

S(a, v)q(x, y, zk 2 − a 2) eiaz e−ivt da dv.

(48)

−a

6.2.      For a uniform motion z(t) = z0 + Ut and with sˆ denoting the Fourier transform of s S(a, v) =

g

+a

s(t) e−ia(z0 + Ut) eivt dt = e−iaz0sˆ (v − aU).

(49)

−a

The speed of the source generates a shift in frequency, a Doppler effect. The pressure is then P(x, y, z, t) =

1 (2p)2

g g +a

−a

0 X

+a

e−iaz0sˆ (v − aU)q x, y,

−a

1

v2 − a 2 eiaz e−ivt da dv. c2

(50)

For a harmonic source of unit amplitude pulsating at the angular frequency v0 , one has s(t) = e−iv0t. Its Fourier transform is sˆ (v) = 2p d(v − v0 ) and P(x, y, z, t) =

1 −iv0t e 2p

g

+a

q(x, y, z(v0 + aU)2/c 2 − a 2) eia(z − z0 − Ut) da.

(51)

−a

The integral is no longer symmetric with respect to 0. One also has (v0 + aU)2/c 2 − a 2 = 0

(52)

for a+ = k/(1 − M),

a− = −k/(1 + M),

(53)

3

   

567

where M = U/c is the Mach number and k = v0 /c. The integral (51) is evaluated numerically for a q 0 over the intervals [0, 0·95a+ ], [0·95a+ , a+ ], [a+ , 1·05a+ ] and [1·05a+ , amax ] where amax is defined by 2 amax − (v0 + amax U)2/c 2 = (2pfmax /c)2.

(54)

Similarly for a Q 0 the integral is calculated on [0·95a− , 0], [a− , 0·95a− ], [1·05a− , a− ] and [amin , 1·05a− ] where amin is the negative solution of equation (54). Each interval is divided into 500 sub-intervals over which a seven points Gauss formula is used. 6.3.  The results are first compared with those of an analytical solution. In free field the retarded potential formula allows one to calculate the pressure created by a point source with a uniform motion and the amplitude s(t) = e−iv0t as (see reference [19]) P(x, t) =

e−iv0(t − =x − y/c) , 4p(=x − y= − (x − y) · U/c)

(55)

where U is the velocity of the source and y is its position at time t − =x − y=/c, which is the solution of y = X(t − =x − y=/c),

(56)

where X(t) is the position of the source at time t. To take into account the rigid ground one must add the pressure created by the image source. So P(x, t) =

e−iv0(t − = x − y=/c) e−iv0(t − =x−y'=/c) + 4p(=x − y= − (x − y) · U/c) 4p(=x − y'= − (x − y') · U/c)

(57)

where y' = X'(t − =x − y'=/c),

(58)

X' being the position of the image source. For the numerical example it is supposed that the source moves along the z-direction at (2, 0·5) in the (x, y)-plane with the speed 50 m/s and radiates a harmonic signal of frequency 200 Hz. The pressure is calculated at point (−2, 0·5, 0). In all the examples the source is supposed to be at z = 0 at time 0. So at this time the source is in front of the observation point. The source sends the signal s(t) = Re(e−iv0t). The signal at the observation point is Re(P(x, y, z, t)) where P is given by formula (57) or (51). P can also be written as an amplitude multiplying the source signal as P(t) = A(t) e−iv0t.

(59)

To avoid an over complex figure Re(eiv0tP(x, y, z, t)) = Re(A(t)) is presented instead of Re(P(x, y, z, t)). In Figure 22 the analytical and numerical results over the time interval [−0·5 s, 0·5 s] are compared. It can be seen that the numerical results give a very good estimate of the pressure. To take more complex examples the straight wall of section 5.2 is again used and the pressure at the observation point behind it is calculated. The source still moves along the z-direction at (2, 0·5) in the (x, y)-plane. Re(A(t)) with P given by formula (51) is again presented.

568

. 

Figure 22. Source moving in a half-space. —, Analytical; ,, numerical.

Figure 23. Pressure for the frequency 200 Hz and the source speed 50 m/s. —, With wall; ·····, without wall.

Figure 24. As Figure 23 but for pressure modulus.

3

   

569

Figure 25. As Figure 23 but for source speed 100 m/s.

For the results in Figure 23 it is supposed that the source has a speed of 50 m/s and sends a signal of frequency 200 Hz. The time history is presented in the time interval [−0·5 s, 0·5 s] with and without the wall. One notices an important decrease in pressure when the wall is added. The behavior is different for the negative and positive times because of the Doppler effect. Figure 24 presents the modulus of P(x, y, z, t) over the time interval [−5 s, 5 s] for the complex signal s(t) = e−iv0t and the same speed of the source. This has the advantage of removing the oscillations and it allows a better estimate of the duration of the noise and of the attenuation provided by the barrier. For Figure 25 the source is assumed to have a speed of 100 m/s. The duration of the noise at the reception point is smaller than in Figure 23 and the frequency appears consequently higher. A stronger Doppler effect is noticeable in this case. On the contrary the maximum of the pressure does not seem to be strongly affected by the speed of the source.

Figure 26. As Figure 23 but for the frequency 800 Hz.

570

. 

In Figure 26 the source has a speed of 50 m/s but sends a signal of frequency 800 Hz. The attenuation provided by the barrier is stronger than for a source at 200 Hz.

7. CONCLUSION

The proposed method allows an efficient determination to be made of excess attenuation of noise barriers of constant cross-sections. One has to solve only two-dimensional problems for real and imaginary frequencies with a very simple one-dimensional mesh. One can obtain the pressure field for a point source and for an incoherent line source. Comparisons can be made between different geometries in this interesting practical case. Furthermore by Fourier transformation one can obtain the time signal at a given point in space for a noise source moving in a direction parallel to the barrier and can show the influence of the frequency and of the velocity of the source. It would be interesting to improve the model to take into account an impedance boundary condition on the barrier and a finite impedance ground.

REFERENCES 1. Z. M 1965 Memoirs of the Faculty of Engineering, Kobe University 11, 29–53. Noise reduction by screens. 2. A. D. P 1974 Journal of the Acoustical Society of America 55, 941–955. Diffraction of sound around corners and over wide barriers. 3. H. G. J 1972 Journal of Sound and Vibration 25, 577–585. Diffraction by wedges of finite acoustic impedance with applications to depressed roads. 4. I. T 1989 Journal of the Acoustical Society of America 85, 661–669. Exact, explicit solutions for diffraction by hard sound barriers and seamonts. 5. U. J. K and G. S. A 1971 Applied Acoustics 4, 35–53. Sound attenuation by barriers. 6. U. J. K 1974 Journal of the Acoustical Society of America 55, 504–518. Noise reduction by barriers. 7. R. S 1980 Journal of Sound and Vibration 73, 195–205. Diffraction of sound around barriers: use of the boundary elements technique. 8. D. C. H, S. N. C-W and M. N. H 1991 Journal of Sound and Vibration 146, 303–322. Efficiency of single noise barriers. 9. D. C. H, D. H. C and S. N. C-W 1991 Applied Acoustics 32, 269–287. The performance of T-profile and associated noise barriers. 10. A. D 1978 Acustica 40, 213–222. Etude de la diffraction par un e´cran mince dispose´ sur le sol. 11. P. J. T. F 1977 Journal of Sound and Vibration 54, 473–500. Layer potentials and acoustic diffraction. 12. Y. K and T. T 1990 Applied Acoustics 31, 101–117. The application of integral equation methods to the calculation of sound attenuation by barriers. 13. H. A 1991 Boundary element methods in acoustics, 225–260 (chapter 11). Applications in environmental noise (Computational mechanics publications). Amsterdam: Elsevier Applied Sciences. 14. C. H. W 1975 Scattering theory for the d’Alembert equation in exterior domains (volume 442 of Lecture Notes in Mathematics). Berlin: Springer-Verlag. 15. I. S. G and I. M. R 1980 Table of integrals, series and products. New York: Academic Press. 16. A. J. B and G. F. M 1971 Proceedings of the Royal Society of London A323, 201–210. The application of integral equation methods to the numerical solution of some exterior boundary-value problems. 17. T. C. L 1984 Journal of Mathematical Analysis and Applications 103, 565–574. A proof for the Burton and Miller integral equation approach for the Helmholtz equation.

3

   

571

18. D. D 1994 Ph.D. Thesis, Ecole Nationale des Ponts et Chausse´es. L’acoustique des proble`mes couple´s fluide-structure: application au controˆle actif du son. 19. D. S. J 1989 Acoustic and electromagnetic waves. Oxford University Press (Science Publications).

APPENDIX: POINT SOURCE SCATTERING BY A RIGID CYLINDER

A rigid cylinder of radius b scatters the sound pressure created by a point source located at (x0 , y0 , z0 ). The pressure is the solution of Dp + k 2p = −d(x − x0 ) d(y − y0 ) d(z − z0 ),

(60)

with 1p/1n = 0 on the surface of the cylinder and the outgoing wave condition at infinity. pˆ denotes the partial Fourier transform of p in the z-direction: pˆ (x, y, kz ) =

g

+a

p(x, y, z) e−ikz z dz

(61)

−a

The Fourier transform is the solution of Dx,ypˆ + (k 2 − kz2 )pˆ = −e−ikz z0 d(x − x0 )d(y − y0 )

(62)

The solution of this problem in the domain outside the cylinder is (see reference [19]) +a

pˆ = s

−a

i iHm (mr0 ) (J (mr)H'm (mb) − Hm (mr)J'm (mb)) e−im(F − F0) e−ikz z0, 4 H'm (mb) m

(63)

with r0 = zx02 + y02, r = zx 2 + y 2 and m = zk 2 − kz2 . F and F0 are respectively the polar angles of points (x, y) and (x0 , y0 ) with respect to the x-axis. Jm and Hm are respectively the Bessel and Hankel functions of order m. This formula is valid for r Q r0 . In the other case one must exchange r and r0 . The sound pressure solution of (60) is thus p=

i +a −im(F − F0) se 8p −a

g

+a

−a

Hm (mr0 ) (J (mr)H'm (mb) − Hm (mr)J'm (mb)) eikz (z − z0) dkz . H'm (mb) m

(64)