- Email: [email protected]
Sound propagation over ground under the influence of a sound speed profile in the atmosphere
Journal of Sound and Vibration (1990) 139(l),
71-81
SOUND PROPAGATION OVER GROUND UNDER THE INFLUENCE OF A SOUND SPEED PROFILE IN THE ATMOSPHERE K. B. RASMUSSEN TheAcoustics Laboratory, Building 352, Technical University of Denmark, DK-2800 Lyngby, Denmark (Received 7 February 1989, and in revised form 11 August 1989)
Sound propagation from a monopole point source located over a semi-infinite porous ground is investigated. The atmosphere is described as a generalized horizontally stratified medium in order to take vertical variations in the sound speed into account. In each layer a linearly varying sound speed is permitted, hence making it possible to simulate a piecewise linear vertical velocity variation. The mathematical problem is formulated in terms of the Hankel transform, and the sound field is determined by numerical integration of the inverse transform. However, the numerical approach is substantially different from the one used in previous work in this field [l-3]. Results of calculations for different atmospheric conditions are shown and discussed, and compared to measured data.
1. INTRODUCTION Determination of the influence of wind and temperature gradients on the propagation of sound over plane ground is an important part of environmental prediction schemes. This is due to the substantial influence from such gradients on the sound levels measured at distances exceeding 100 m. From a practical point of view, something has to be done in order to take these influences into account, and since the theoretical side of the problem is still largely unsolved, one possibility is to incorporate results from extensive measurements supported by ray-tracing calculations, for which the degree of accuracy is uncertain. Attempts to address the problem more directly by using numerical solutions of the wave equation have been made, however [2-41, with encouraging results, except for serious numerical stability problems. Yet another numerically oriented solution approach is presented in what follows. In principle this approach is very similar to previous approaches, but it does not suffer from any stability problems. The method is to implement a direct numerical solution to the wave equation in a medium (air) having piecewise linear velocity variations. The ground is described as being a porous material, characterized by the formulas provided by Delany and Bazley in their now famous paper. This choice of ground model has gained widespread acceptance. The difference from the works by previous authors rests in the choice of formulation of the mathematical problem. In previous papers dealing with the subject [2,3], methods similar to the so-called Thomson-Haskell solution for a horizontally layered structure have been used. However, it is believed that more efficient and stable results may be obtained by employing the more direct implementation originating from Schmidt and Jensen [5]. This implementation is known in seismic applications to be approximately ten times faster than the comparable Thomson-Haskell algorithm. Hence, the problem addressed in the present paper is fundamentally one of a multi-layer structure having horizontal stratification. The lowest layer is always a porous medium, 71 0022460x/90/1coO71+
11 $03.00/O
@ 1990 Academic FVess Limited
12
K. R.
RASMUSSEN
which is described as a fluid behaving according to Delany and Bazley formulas. The usual impedance description has not been applied: i.e., local reaction has not been incorporated explicitly in the equations. The calculations presented in this paper are, however, indistinguishable from local reaction results, since the losses and the propagation velocity in the ground are such that differences disappear. This has been verified for the case of zero gradient. (Local reaction results are obtained when the sound speed in the ground is approaching zero, and/or the losses approach infinity.) The upper layers representing the atmosphere may each have a linear sound speed variation as a function of height. In principle, the number of atmospheric layers is unlimited, but in practice the number is limited due to considerations of calculation time and memory size. For calculations involving many layers, numerical problems could occur since a linear equations system involving 2N - 2 equations for N layers is to be solved. The rest of the paper deals with the mathematical basis for the implementation, and some appropriate comments on how various numerical problems have been solved. In the last part of the paper some numerical results are shown, the principle purpose of which is to demonstrate the influence of sound speed gradients on atmospheric near-to-theground propagation. The velocity gradients are interpreted as equivalent wind speed gradients, and the numerical results are compared with measured data. 2. THEORY The mathematical description given in this section is a combination of the multi-layerproblem solution given by Schmidt and Jensen [5] and as used by Rasmussen and Vistisen [6] with Airy function solutions to problems with linear velocity variations [l, 71. VELOCITY LAYERS 2.1. CONSTANT For a horizontally layered structure composed of solely fluid layers as shown in Figure 1, the acoustic field is expressed directly in terms of pressure p, defined as a solution to the wave equation
V2p + k2p = 0, where k is the wavenumber.
The vertical displacement w =
(-
dp)-’
(1)
component
w is given by
ap/az,
(2)
where o is angular frequency and p is density. The eiW’time convention is used. Introducing the Hankel transform, which is suitable due to the symmetry around the vertical z-axis, one obtains [5] 3: K,(s, z) Jo(sr)s ds, P(C z) = (3) I0 where the transformed
pressure is given by the kernel K,(s, z) = A-(s) e-‘*‘“‘+
Figure 1. Horizontal
stratification.
Layer no. 1 is the ground.
A+(s) e+zu’s),
The source is assumed
(4)
to be located
on the z-axis.
OUTDOOR
where A+ and A- are determined layers. a is given by
SOUND
according
73
PROPAGATION
to the boundary
conditions
between the
@y(S)= (sZ- I?)“*,
(5)
where k is the wavenumber for the layer in question. The wavenumber is complex due to the losses in the medium, and is thus defined as k = w/c - ia, where a represents the attenuation. Since the wavenumber is complex, it should be specified that the square root in equation (5) is understood to have positive real part. From the Hankel transform it is seen that s is the horizontal wavenumber component and that in the transformed domain the coefficients A’ and A- do not depend on z. For the vertical displacement one obtains, from equations (2) and (3), X K,. J,)(sr)s ds, w(r,z)=(-w?p)-’ (6) I0 where the kernel is given by K,, = aAEquation (1) is for the source-free of the relation for a point source,
emZU(‘)+aAf(S) et”““.
(7)
field. Including the source field is done by means X
p*(r, z) =
I0
i,, JO(sr)s ds,
(8)
where the kernel is given by & = e -I;-=,l*/(y,
(9)
where z, is the source height. The asterisk denotes source contribution. The field given by equation (8) is equal to eeikR/ R, which is the appropriate expression for a point source when the eiw’time convention is used. For the vertical displacement one obtains i?(r, z) = (-
w”p)-’
cc J k,.
JO(rs)s ds,
(10)
0
where the kernel is given by
k,=-sign(z-~,)e~(~-‘.(~.
(11)
In each layer the kernels given in equation (4) are defined separately, and (Yas well as A+ and A- are layer-dependent parameters, which are more accurately written as (Y,, AZ and A,, where n is the layer number. Finding AZ and Ai as a function of s is the central part of the field calculation. The integration from zero to infinity is done numerically and the necessary integration interval is limited, since the kernel approaches zero for large and small values of the horizontal wavenumber. For each s-value used in the numerical integration, the value of the A-parameters are found by means of solving 2N -2 linear equations, where N is the number of layers; see Figure 1. The linear equations are determined from continuity of w and p at the boundaries. Equations are formulated for each layer boundary on the basis of the kernels (4) and (7). For the boundaries of the layer containing the source, contributions from the kernels (9) and (11) are added. Hence, each layer n gives rise to equations of the type K W, n,” -K W,fl+l,l = k,.,,
- L.,+ ,.,r
K,.,, - &,n+,,,= &,“.U - &.“+l.,,
(12, 13)
where the subscripts n and n + 1 denote the layer number, whereas the subscripts I and u denote the lower and upper interface respectively for the layer. The source contributions on the right side will be zero except in the source layer.
K. B. RASMUSSEN
74
It is an advantage to scale equations (12) and (13) in order to obtain coefficients of the same order of magnitude. Hence, displacement equations are multiplied by l/k, where k is the wavenumber in the source layer. A local z co-ordinate is introduced in each layer, in such a way that z is zero on the lower interface. This choice saves computation time, since exponential functions having zero argument may be left out. Furthermore, the choice of the local z co-ordinate is a part of the precautions ensuring numerical stability of the equations system. If the equations are arranged in a matrix structure, as shown schematically in Figure 2, then a band matrix emerges, which may readily be solved by means of Gaussian elimination. The A- and A’ values for the layer containing the observation point are used to determine the contribution to the pressure integral, equations (3) and (4).
2 2 Figure 2. Schematic last layer A+ is zero.
representation
of matrix for a three-layer
problem.
In the first layer A-
is zero; in the
2.2. LINEAR VELOCITY VARIATION In the previous section the stratification consisted solely of homogenous layers. It is possible to introduce layers having a linear sound speed variation as a function of the vertical z co-ordinate. The solution is based upon the Airy function. The wavenumber in a layer having a velocity gradient y is given by
VI),
k(z) = wI(co[l+
(14)
leading to (valid for (yz(
k,= o/co.
(15)
Inserting the approximation in equation (15) in the wave equation leads to Airy’s differential equation, and hence to Airy function solutions. The solution may be stated as [I,%71 CO KP J,,(sr)s ds, (16) P(C z) = I0 where K,=A-(s)V(~+y)+A+(s)w(~+y), 7=(s2-
k;)l’,
y=z/4
I=
sign (r)((rl2ki)-““,
(17)
and V and W are related to the Airy function as follows [7,8]: V(z) = 7~“’ Ai (z),
The displacement,
W(z) = 2q7’12eirr/’ Ai (z e’*“/‘).
(18919)
w, is determined through equation (2) as before. The kernel becomes K,=A-(s)V’(~+y)/l+A+(s)W’(~+y)/l,
(20)
OUTDOOR
SOUND
PROPAGATION
75
where V’and W’ denote differentiation with respect to the height z. Hence, determining the displacement involves calculating the derivatives of the V and W functions. For large values of the real part of 7 the behaviour of W and V and their derivatives are dominated by exponential functions, g(x) = e-Y, f(x) = e:““‘, (21322) where x is the real part of 7. f represents the asymptotic behaviour of W and W’, and g the behaviour of V and V’. Hence, for numerical reasons it is advisable to multiply W and W’by l/f(x) whenever the real part of 7 is positive, and to multiply V and V’by l/g(x) for positive real part of 7. In fact, this modification ensures numerical stability of the equations system (see Figure 2), because very large and very small numerical values are removed by the modification. The legality of the process is ensured by the fact that 7 is independent of the variable z. Source contributions are not included, since-it has been assumed in this work that the source is not located in a gradient layer. This is not a problem as far as the application of the computer model is concerned, since the source may always be located in a thin homogenous layer. Calculating the field for a problem involving one or more gradient layers means setting up the linear equations as described in section 2.1, but in gradient layers the kernels are expressed by equations (17) and (20). Due to the assumption that the source is always located in a homogeneous layer, the source contribution is always determined from the expressions in section 2.1. When gradient layers are present, for numerical reasons the problem should always be formulated so that the source height is greater than the receiver height. When this condition is not met, the principle of reciprocity is employed: the source and receiver are interchanged and the signs of the gradients are also changed. For the situation dealt with in this paper, this is in fact equivalent to an interchange of source and receiver heights. 2.3. THE GROUND LAYER The ground layer is described as a fluid layer of infinite depth and having heavy losses. Wavenumber and effective density should be determined according to models describing the acoustic properties of porous materials and, in the present work, the formulas from Delany and Bazley [9] are used: pc=p,,c,[l+0~0511(f/a)-0’7’-i0~0768(f/a)-0’73], k= Ic,Jl+0~0858(f/a)-0’70-iO~1749(~/a)-0’59],
(23) (24)
where u denotes the specific resistance of the ground in Ns/m4. po, co and k. denote density, speed of sound and wavenumber for air, respectively. Note that the density becomes complex due to the friction. For a homogeneous atmosphere it has been established that the results from what is an extended reaction calculation with use of equations (23) and (24) are indistinguishable from the results obtained from a local reaction calculation. 2.4. REMARKS CONCERNING THE IMPLEMENTATION Since the above formulas for the sound pressure are not suited for analytical evaluation, but are evaluated by numerical approaches, and since the numerical evaluation involves a number of important decisions, some remarks on the subject are appropriate. The integration path could be deformed away from the real axis, in order to improve convergence when performing numerical integration. Generally speaking such a step would, however, necessitate tracing the location of the poles, since allowance must be
76
K. 9.
RASMUSSEN
made for the passage of poles during deformation. In the present case, the integration is performed along the real axis, and the presence of losses (mainly in the ground layer) ensures that no poles are located on the axis. If the near field is without interest, the integrals to be evaluated may be rewritten by employing the asymptotic representation for the Bessel function. Such an approach brings the formulation into a form which permits numerical integration by means of the efficient Fast Fourier Transform algorithm [4-61. Furthermore, approaches of this kind are suitable for implementation on vector processors and, hence, very efficient processing is the result. On the other hand, if an adaptive integration algorithm is employed, varying the step size according to the more or less rapid variations in the kernel, then the near field may be retained while saving integration points in regions where the kernel is slowly varying as a function of the horizontal wavenumber. Retaining the unapproximated Bessel function also means avoiding the introduction of an extra pole in the kernel due to the asymptotic Bessel function expression. This approach has been employed in the present work. The integration is performed with a density of points determined by two criteria: the first criterion is that the approximate Bessel function periodicity is taken into account by at least six points per period; the second criterion is that the contribution to the integral from one segment must be smaller than 0.03 JO(rs). For each point in the integration, Bessel and Airy functions are evaluated. The Bessel function of zero order is evaluated by means of rational approximations found in reference [lo], and the Airy functions and their derivatives are evaluated as described in the Appendix. The kernel in the integrand is found as a solution to linear equations having band structure and as defined in the previous sections. The linear equations are solved by means of Gaussian elimination with pivoting, and the band structure is taken into account in order to reduce the calculation time. The computation time is determined by the number of kernel points to be determined, since each kernel value is found from the solution of the linear equations system. However, the computation time may be reduced in the following way. For large s-values the corresponding kernel contribution represents near-horizontal propagation and evanescent waves. It may be shown from analytic solutions to three-layer problems that when s approaches infinity (the evanescent region) the kernel depends only on s and source and receiver height. The layer-dependent parameters disappear and the calculation is essentially for an infinite medium. This means that the kernel may be approximated by the kernel corresponding to a simplified problem for large s values. In a similar but more intuitive manner, the problem may be simplified for s approaching zero. Small values of s represent near-vertical propagation contributions. For vertical contributions a horizontal sound speed gradient is without effect except for a phase shift as the contribution propagates through areas of different propagation velocity. The introduction of layers in the atmosphere does lead to reflection of sound when the layer boundaries are passed. These reflected contributions are, however, due to the discontinuity in the derivative of the sound speed gradient, and hence a second order effect which may safely be ignored. In fact they represent a weakness in the problem formulation rather than a true effect. The result is that the kernel may safely be approximated by the kernel for a ground homogenous air interface for small s values. The Thomson-Haskell approach is the usual one to apply to this class of problems. It is a matrizant method involving propagation of the field through all layers in the problem. This,leads to numerical problems for large values of the horizontal wavenumber. In the present approach (originating from that of Schmidt and Jensen [S] for non-gradient layers) the continuity at the layer boundaries is established by the linear equations system,
OUTDOOR
SOUND
77
PROPAGATION
which is arranged in such a manner that the most important contributions are close to the diagonal in the matrix equation. This means that the Gaussian elimination program removes the numerical stability problems.
3. NUMERICAL
RESULTS
The numerical tool presented in the previous sections may be applied to the simulation of various outdoor sound propagation conditions over plane ground. Possibilities of influence from extended reaction phenomena in connection with realistic ground models may be investigated, but that is not dealt with here. The main reason for developing the program was to investigate the theoretical importance of the presence of a wind speed gradient in the lower atmosphere. A model with a linear sound speed variation had been developed previously [I], but the present work makes it possible to establish the influence from the fact that the actual gradient is not adequately described by a simple linear velocity variation. Strictly speaking, the mathematical model deals with sound speed variations and not wind speed variations, but the effective sound speed may be viewed as the superposition of the physical sound speed and the wind speed component in the direction of sound propagation. In the following, results for a source height of 1.45 m and a receiver height of 0.5 m are shown for various sound speed conditions. The ground is specified by the Delany and Bazley relations, equations (23) and (24), for a flow resistance, a, given by 200x 10’ Ns/m4. The results are compared with measured one-third octave data from grasscovered ground previously published in references [l] and [ 111. The measured data were obtained for wind speeds of 2-2.5 m/s measured at a height of 10 m. Atmospheric absorption has been taken into account in the calculations by introducing a loss of O-002 dB per wavelength. Measured and calculated data are valid for a distance of 80 m except for that used in Figure 3 which is valid for 20 m. Figure 3 is included in order to demonstrate the deviation between theoretical and experimental values for a short distance where the wind gradient has only small influence. Calculations have been performed for the positive gradients shown in Figure 4 and for negative gradients having the same shape and size. One gradient is a simple linear variation
-20
-30
I
I
I
I 250
I 500
I 1000
-
125
Frequency
2000
(Hz)
Figure 3. Sound pressure level relative to free field over a typical grassy field. Horizontal distance 20m; source height 1.45 m; receiver height 0.5 m. --, Measured data, downwind [l]; - - -, measured data, upwind [I]; --, predicted for zero gradient.
78
K. R. RASMUSSEN
up to 20 m, and the others are bilinear up to 20 m. The bilinear cases simulate the meteorological fact that for a wind speed gradient and gradient is steeper close to the ground. All the gradients in Figure 4 are, however, terminated at a height of 20 m. This termination is without influence on the results but ensures that only a well defined range of sound velocities enter the calculations. Since the source must be in a layer with zero gradient, such a thin (O-1 m) layer is included in the calculations. This layer is not visible in Figure 4. An indication is given in Figure 5 that for a positive gradient some improvement of the similarity between measured and calculated data may be obtained by the inclusion of a bilinear velocity variation in comparison with the simple linear variation. It is indicated in Figure 6 that the bilinear cases overestimate the shadow-effect occurring at high frequencies due to the negative gradient. The size of this deviation is dependent
20-
343
344
3 IS
C (m/a)
Figure 4. Sound speed profiles used in calculations. 1, Simple gradient; 2, bilinear variation crossing simple gradient at 10 m; 3, bilinear variation crossing simple gradient at 20 m.
I
I
I
Frequency
(Hz1
Figure 5. Sound pressure level relative to free field for a positive velocity gradient over a typical grassy field. Horizontal distance 80 m; source height 1.45 m; receiver height 0.5 m. -, Measured data from reference [ 11; . . . . predicted for simple gradient: - - -, predicted for bilinear case no. 2; - --, predicted for bilinear case no. 31
OUTDOOR
%
SOUND
79
PROPAGATION
-10 -
-2o-
-30 125
250
500
Frequency (Hz) Figure 6. Sound pressure level relative to free field for a negative velocity gradient over a typical grassy field. Horizontal distance 80 m; source height 1.45 m; receiver height 0.5 m. -, Measured data from reference [ 11; . . I . . predicted for simple gradient; - - -, predicted for bilinear variation, case no. 2; --, predicted for bilinear Lariation, case no. 3.
upon the choice of parameters for the bilinear curve; the lower part of the sound speed variation is especially critical for high-frequency propagation. However, the deviation could also be severely influenced by the fact that turbulence effects are unaccounted for in the simulation. In fact, the shape of the high-frequency part of the measured curve has the approximate 20 dB per decade slope which is to be expected from turbulence considerations [ 111. It should be noted that the measured difference in level between the positive and negative gradient cases is more than 10 dB at 2 kHz. From a practical point of view, this large difference is well accounted for by the simple linear sound speed profile simulation. From a theoretical point of view the use of a linear profile is, of course, not
2-
OF
I
01
I
7,877
I
I
I
la.052
la.228
la.403
Horzontol
I
la.579
I
wovenumber
Figure 7. Pressure kernel for positive bilinear variation, case no. 3, for 1000 Hz. -, imaginary part.
Real part; --,
80
K. B. RASMUSSEN
-II 17.701
I I7.877
I I8.052
I
18~228
Hor~zontol
Figure 8. Pressure imaginary part.
kernel
for negative
bilinear
I 18.403
I 18.579
I 19-754
wavenumber
variation,
case no. 3, for 1000
Hz. -,
Realpart; --,
satisfactory. Ideally, a very accurate sound speed profile (consisting of a considerable number of line segments) should be used in connection with turbulence considerations. The numerical integration process is illustrated in Figures 7 and 8, showing the central part of the kernel for increasing and decreasing sound speed versus height at 1 kHz. As might be expected the kernel, which represents the wavenumber spectrum, has a more pronounced mode structure for increasing sound speed versus height than for decreasing sound speed versus height. This is a consequence of the concentration of sound energy along the surface for a positive gradient. On the basis of the material in this paper it seems that a simple linear sound speed profile could be a reasonable representative of the combined effect of turbulence and sound speed gradients in the atmosphere. Further measured data representing point source propagation under well defined meteorological conditions are needed in order to test this preliminary conclusion. Furthermore, indications that an improved prediction may be obtained for a more realistic sound speed profile are present, but such a more realistic profile makes it more urgent to introduce turbulence in the calculations.
REFERENCES 1. K. 9. RASMUSSEN 1986 Journal ofSound and Vibration 104,321-335. Outdoor sound propagation under the influence of wind and temperature gradients. KUESTER,D.C.CHANG,W.F.RICHARDS,R.GILBERT~~~ N. 2. R.RAsPET,S.W.LEE,E. BONG 1985 Journal of the Acoustical Society of America 77, 345-352. A fast field program for sound propagation in a layered atmosphere above an impedance ground. 3. S. W. LEE, N. BONG, W. F. RICHARDSand R. RASPET1986 Journal of the Acoustical Society of America 79, 628-634. Impedance formulation of the fast field program for acoustic wave propagation in the atmosphere. 4. T. L. RICHARDSand K. AI-I-ENBOROUGH1986 Journal ofSound and Vibration 109,157-167. Accurate FFT-based Hankel transforms for predictions of outdoor sound propagation. 5. H. SCHMIDT and F. 9. JENSEN 1985 Journal ofrhe Acoustical Society ofAmerica 77,813-825. A full wave solution for propagation in multilayered viscoelastic media with application to Gaussian beam reflection at fluid-solid interfaces. and 9. VIS-~ISEN1986 Repot-f86.149 Qdegaard & Danneskiold-Samstie 6. K. 9. RASMUSSEN Kroghsgade 1 DK-2100, Copenhagen, Denmark. Calculation of sound propagation under water by means of FIT-model. (In Danish.) I. A. D. PIERCE 1981 Acoustics: An Introduction to its Physical Principles and Applications. New York: McGraw-Hill. 8. V. A. FOCK 1965 Electromagnetic Diflraction and Propagation Problems. Oxford: Pergamon Press.
OUTDOOR
SOUND
81
PROPAGATION
and E. N. BAZLEY 1970 Applied Acoustics 3, 105-116.Acoustical properties of fibrous absorbent materials. 10. M. ABRAMOWITZ and I. A. STEGUN 1970 Handbook cf Mathematical Functions. New York: Dover. 11. K. B. RASMUSSEN 1985 Danish Acoustical Institute Report. The effect of wind and temperature gradients on sound propagation outdoors.
9. M. E. DELANY
.APPENDIX:
ON
THE
CALCULATION OF AIRY DERIVATIVES
FUNCTIONS
AND
THEIR
The W and V functions defined in terms of the Airy function in equations (18) and (19) are complex functions of complex argument. They may be evaluated on the basis of formulas from the book by Abramowitz and Stegun [lo] as follows:
Ai
= c,.f(z) - c,g(z), cZ = - Ai’(0) = 0.258819,
c, = Ai(0) = 0.355028, f(z)
=
(Al)
Z0 1+“‘ + z6 2.3 2.3.5.6+2.3.5.6.8.9 Z
+...,
V and
W are then determined
(A3)
IO
+. . .
g(z)=z+L+z:
3.4
(A2)
(A4)
3.4.6.7+3.4.6.7.9.10
from W(z) = 2.,.r’/’ e’njh Ai (z e’2n/‘).
V(z) = T”’ Ai (z),
(A5, A6)
The derivatives V’ and W’ are calculated from the similar expressions obtained when equations (A3) and (A4) are differentiated. Equation (Al) is slow in computation due to the series (A3) and (A4). For large positive and negative arguments (real part) asymptotic expressions may be used [8, lo]. For Re (z) >>1 one obtains (including three terms) U(z) = z -“4 e”( 1 + a,/x
+ u,/x” + a,/~~),
U’(z) = z’14 e”( 1 - b,/x - b,/x’-
(A7)
b,/x3),
(A8)
V( 2) = o*5z-“4
e-“( 1 - a,/x
+ a,/.~‘-
a,/~“),
(A9)
V’(z) = -0.5~“~
e-‘( 1+ b/x
- b,/x’+
b,/x3),
(AlO)
where x = (2/3)z”* and a, = O-069444, a, = O-037133, a3 = 0.037993, b, = O-097222, b2 = 0.043885 and b, = 0.042463. W and W’ are determined from equations (A7)-(AlO) with W(z)=
For negative U(-z) U’( -z)
values
one similarly
U(z)+iV(z).
obtains
= zz”4 cos (x + 7r/4)(1-
(using
(All) the same definitions
UJX’) + zz”J sin (x+ 7r/4)(a,/x
= z”4 sin (x + 7r/4)( 1 + b2/x2) + z”’ cos (x + r/4)(
V(-z)
= zP1’4 sin (x + r/4)(
V’( -z)
= - z’/4 cos (~+~/4)(1+b~/x’)+z”~sin
of CJ and x) - a,/~‘),
(A12)
b,/x - b,/x3),
(A13)
1 - a2/x2) - 2-“4 cos (x + ~/4)(a,/x
- a,/~~),
(x+n/4)(b,/x-b,/x’).
(A14) (A15)
71-81
SOUND PROPAGATION OVER GROUND UNDER THE INFLUENCE OF A SOUND SPEED PROFILE IN THE ATMOSPHERE K. B. RASMUSSEN TheAcoustics Laboratory, Building 352, Technical University of Denmark, DK-2800 Lyngby, Denmark (Received 7 February 1989, and in revised form 11 August 1989)
Sound propagation from a monopole point source located over a semi-infinite porous ground is investigated. The atmosphere is described as a generalized horizontally stratified medium in order to take vertical variations in the sound speed into account. In each layer a linearly varying sound speed is permitted, hence making it possible to simulate a piecewise linear vertical velocity variation. The mathematical problem is formulated in terms of the Hankel transform, and the sound field is determined by numerical integration of the inverse transform. However, the numerical approach is substantially different from the one used in previous work in this field [l-3]. Results of calculations for different atmospheric conditions are shown and discussed, and compared to measured data.
1. INTRODUCTION Determination of the influence of wind and temperature gradients on the propagation of sound over plane ground is an important part of environmental prediction schemes. This is due to the substantial influence from such gradients on the sound levels measured at distances exceeding 100 m. From a practical point of view, something has to be done in order to take these influences into account, and since the theoretical side of the problem is still largely unsolved, one possibility is to incorporate results from extensive measurements supported by ray-tracing calculations, for which the degree of accuracy is uncertain. Attempts to address the problem more directly by using numerical solutions of the wave equation have been made, however [2-41, with encouraging results, except for serious numerical stability problems. Yet another numerically oriented solution approach is presented in what follows. In principle this approach is very similar to previous approaches, but it does not suffer from any stability problems. The method is to implement a direct numerical solution to the wave equation in a medium (air) having piecewise linear velocity variations. The ground is described as being a porous material, characterized by the formulas provided by Delany and Bazley in their now famous paper. This choice of ground model has gained widespread acceptance. The difference from the works by previous authors rests in the choice of formulation of the mathematical problem. In previous papers dealing with the subject [2,3], methods similar to the so-called Thomson-Haskell solution for a horizontally layered structure have been used. However, it is believed that more efficient and stable results may be obtained by employing the more direct implementation originating from Schmidt and Jensen [5]. This implementation is known in seismic applications to be approximately ten times faster than the comparable Thomson-Haskell algorithm. Hence, the problem addressed in the present paper is fundamentally one of a multi-layer structure having horizontal stratification. The lowest layer is always a porous medium, 71 0022460x/90/1coO71+
11 $03.00/O
@ 1990 Academic FVess Limited
12
K. R.
RASMUSSEN
which is described as a fluid behaving according to Delany and Bazley formulas. The usual impedance description has not been applied: i.e., local reaction has not been incorporated explicitly in the equations. The calculations presented in this paper are, however, indistinguishable from local reaction results, since the losses and the propagation velocity in the ground are such that differences disappear. This has been verified for the case of zero gradient. (Local reaction results are obtained when the sound speed in the ground is approaching zero, and/or the losses approach infinity.) The upper layers representing the atmosphere may each have a linear sound speed variation as a function of height. In principle, the number of atmospheric layers is unlimited, but in practice the number is limited due to considerations of calculation time and memory size. For calculations involving many layers, numerical problems could occur since a linear equations system involving 2N - 2 equations for N layers is to be solved. The rest of the paper deals with the mathematical basis for the implementation, and some appropriate comments on how various numerical problems have been solved. In the last part of the paper some numerical results are shown, the principle purpose of which is to demonstrate the influence of sound speed gradients on atmospheric near-to-theground propagation. The velocity gradients are interpreted as equivalent wind speed gradients, and the numerical results are compared with measured data. 2. THEORY The mathematical description given in this section is a combination of the multi-layerproblem solution given by Schmidt and Jensen [5] and as used by Rasmussen and Vistisen [6] with Airy function solutions to problems with linear velocity variations [l, 71. VELOCITY LAYERS 2.1. CONSTANT For a horizontally layered structure composed of solely fluid layers as shown in Figure 1, the acoustic field is expressed directly in terms of pressure p, defined as a solution to the wave equation
V2p + k2p = 0, where k is the wavenumber.
The vertical displacement w =
(-
dp)-’
(1)
component
w is given by
ap/az,
(2)
where o is angular frequency and p is density. The eiW’time convention is used. Introducing the Hankel transform, which is suitable due to the symmetry around the vertical z-axis, one obtains [5] 3: K,(s, z) Jo(sr)s ds, P(C z) = (3) I0 where the transformed
pressure is given by the kernel K,(s, z) = A-(s) e-‘*‘“‘+
Figure 1. Horizontal
stratification.
Layer no. 1 is the ground.
A+(s) e+zu’s),
The source is assumed
(4)
to be located
on the z-axis.
OUTDOOR
where A+ and A- are determined layers. a is given by
SOUND
according
73
PROPAGATION
to the boundary
conditions
between the
@y(S)= (sZ- I?)“*,
(5)
where k is the wavenumber for the layer in question. The wavenumber is complex due to the losses in the medium, and is thus defined as k = w/c - ia, where a represents the attenuation. Since the wavenumber is complex, it should be specified that the square root in equation (5) is understood to have positive real part. From the Hankel transform it is seen that s is the horizontal wavenumber component and that in the transformed domain the coefficients A’ and A- do not depend on z. For the vertical displacement one obtains, from equations (2) and (3), X K,. J,)(sr)s ds, w(r,z)=(-w?p)-’ (6) I0 where the kernel is given by K,, = aAEquation (1) is for the source-free of the relation for a point source,
emZU(‘)+aAf(S) et”““.
(7)
field. Including the source field is done by means X
p*(r, z) =
I0
i,, JO(sr)s ds,
(8)
where the kernel is given by & = e -I;-=,l*/(y,
(9)
where z, is the source height. The asterisk denotes source contribution. The field given by equation (8) is equal to eeikR/ R, which is the appropriate expression for a point source when the eiw’time convention is used. For the vertical displacement one obtains i?(r, z) = (-
w”p)-’
cc J k,.
JO(rs)s ds,
(10)
0
where the kernel is given by
k,=-sign(z-~,)e~(~-‘.(~.
(11)
In each layer the kernels given in equation (4) are defined separately, and (Yas well as A+ and A- are layer-dependent parameters, which are more accurately written as (Y,, AZ and A,, where n is the layer number. Finding AZ and Ai as a function of s is the central part of the field calculation. The integration from zero to infinity is done numerically and the necessary integration interval is limited, since the kernel approaches zero for large and small values of the horizontal wavenumber. For each s-value used in the numerical integration, the value of the A-parameters are found by means of solving 2N -2 linear equations, where N is the number of layers; see Figure 1. The linear equations are determined from continuity of w and p at the boundaries. Equations are formulated for each layer boundary on the basis of the kernels (4) and (7). For the boundaries of the layer containing the source, contributions from the kernels (9) and (11) are added. Hence, each layer n gives rise to equations of the type K W, n,” -K W,fl+l,l = k,.,,
- L.,+ ,.,r
K,.,, - &,n+,,,= &,“.U - &.“+l.,,
(12, 13)
where the subscripts n and n + 1 denote the layer number, whereas the subscripts I and u denote the lower and upper interface respectively for the layer. The source contributions on the right side will be zero except in the source layer.
K. B. RASMUSSEN
74
It is an advantage to scale equations (12) and (13) in order to obtain coefficients of the same order of magnitude. Hence, displacement equations are multiplied by l/k, where k is the wavenumber in the source layer. A local z co-ordinate is introduced in each layer, in such a way that z is zero on the lower interface. This choice saves computation time, since exponential functions having zero argument may be left out. Furthermore, the choice of the local z co-ordinate is a part of the precautions ensuring numerical stability of the equations system. If the equations are arranged in a matrix structure, as shown schematically in Figure 2, then a band matrix emerges, which may readily be solved by means of Gaussian elimination. The A- and A’ values for the layer containing the observation point are used to determine the contribution to the pressure integral, equations (3) and (4).
2 2 Figure 2. Schematic last layer A+ is zero.
representation
of matrix for a three-layer
problem.
In the first layer A-
is zero; in the
2.2. LINEAR VELOCITY VARIATION In the previous section the stratification consisted solely of homogenous layers. It is possible to introduce layers having a linear sound speed variation as a function of the vertical z co-ordinate. The solution is based upon the Airy function. The wavenumber in a layer having a velocity gradient y is given by
VI),
k(z) = wI(co[l+
(14)
leading to (valid for (yz(
k,= o/co.
(15)
Inserting the approximation in equation (15) in the wave equation leads to Airy’s differential equation, and hence to Airy function solutions. The solution may be stated as [I,%71 CO KP J,,(sr)s ds, (16) P(C z) = I0 where K,=A-(s)V(~+y)+A+(s)w(~+y), 7=(s2-
k;)l’,
y=z/4
I=
sign (r)((rl2ki)-““,
(17)
and V and W are related to the Airy function as follows [7,8]: V(z) = 7~“’ Ai (z),
The displacement,
W(z) = 2q7’12eirr/’ Ai (z e’*“/‘).
(18919)
w, is determined through equation (2) as before. The kernel becomes K,=A-(s)V’(~+y)/l+A+(s)W’(~+y)/l,
(20)
OUTDOOR
SOUND
PROPAGATION
75
where V’and W’ denote differentiation with respect to the height z. Hence, determining the displacement involves calculating the derivatives of the V and W functions. For large values of the real part of 7 the behaviour of W and V and their derivatives are dominated by exponential functions, g(x) = e-Y, f(x) = e:““‘, (21322) where x is the real part of 7. f represents the asymptotic behaviour of W and W’, and g the behaviour of V and V’. Hence, for numerical reasons it is advisable to multiply W and W’by l/f(x) whenever the real part of 7 is positive, and to multiply V and V’by l/g(x) for positive real part of 7. In fact, this modification ensures numerical stability of the equations system (see Figure 2), because very large and very small numerical values are removed by the modification. The legality of the process is ensured by the fact that 7 is independent of the variable z. Source contributions are not included, since-it has been assumed in this work that the source is not located in a gradient layer. This is not a problem as far as the application of the computer model is concerned, since the source may always be located in a thin homogenous layer. Calculating the field for a problem involving one or more gradient layers means setting up the linear equations as described in section 2.1, but in gradient layers the kernels are expressed by equations (17) and (20). Due to the assumption that the source is always located in a homogeneous layer, the source contribution is always determined from the expressions in section 2.1. When gradient layers are present, for numerical reasons the problem should always be formulated so that the source height is greater than the receiver height. When this condition is not met, the principle of reciprocity is employed: the source and receiver are interchanged and the signs of the gradients are also changed. For the situation dealt with in this paper, this is in fact equivalent to an interchange of source and receiver heights. 2.3. THE GROUND LAYER The ground layer is described as a fluid layer of infinite depth and having heavy losses. Wavenumber and effective density should be determined according to models describing the acoustic properties of porous materials and, in the present work, the formulas from Delany and Bazley [9] are used: pc=p,,c,[l+0~0511(f/a)-0’7’-i0~0768(f/a)-0’73], k= Ic,Jl+0~0858(f/a)-0’70-iO~1749(~/a)-0’59],
(23) (24)
where u denotes the specific resistance of the ground in Ns/m4. po, co and k. denote density, speed of sound and wavenumber for air, respectively. Note that the density becomes complex due to the friction. For a homogeneous atmosphere it has been established that the results from what is an extended reaction calculation with use of equations (23) and (24) are indistinguishable from the results obtained from a local reaction calculation. 2.4. REMARKS CONCERNING THE IMPLEMENTATION Since the above formulas for the sound pressure are not suited for analytical evaluation, but are evaluated by numerical approaches, and since the numerical evaluation involves a number of important decisions, some remarks on the subject are appropriate. The integration path could be deformed away from the real axis, in order to improve convergence when performing numerical integration. Generally speaking such a step would, however, necessitate tracing the location of the poles, since allowance must be
76
K. 9.
RASMUSSEN
made for the passage of poles during deformation. In the present case, the integration is performed along the real axis, and the presence of losses (mainly in the ground layer) ensures that no poles are located on the axis. If the near field is without interest, the integrals to be evaluated may be rewritten by employing the asymptotic representation for the Bessel function. Such an approach brings the formulation into a form which permits numerical integration by means of the efficient Fast Fourier Transform algorithm [4-61. Furthermore, approaches of this kind are suitable for implementation on vector processors and, hence, very efficient processing is the result. On the other hand, if an adaptive integration algorithm is employed, varying the step size according to the more or less rapid variations in the kernel, then the near field may be retained while saving integration points in regions where the kernel is slowly varying as a function of the horizontal wavenumber. Retaining the unapproximated Bessel function also means avoiding the introduction of an extra pole in the kernel due to the asymptotic Bessel function expression. This approach has been employed in the present work. The integration is performed with a density of points determined by two criteria: the first criterion is that the approximate Bessel function periodicity is taken into account by at least six points per period; the second criterion is that the contribution to the integral from one segment must be smaller than 0.03 JO(rs). For each point in the integration, Bessel and Airy functions are evaluated. The Bessel function of zero order is evaluated by means of rational approximations found in reference [lo], and the Airy functions and their derivatives are evaluated as described in the Appendix. The kernel in the integrand is found as a solution to linear equations having band structure and as defined in the previous sections. The linear equations are solved by means of Gaussian elimination with pivoting, and the band structure is taken into account in order to reduce the calculation time. The computation time is determined by the number of kernel points to be determined, since each kernel value is found from the solution of the linear equations system. However, the computation time may be reduced in the following way. For large s-values the corresponding kernel contribution represents near-horizontal propagation and evanescent waves. It may be shown from analytic solutions to three-layer problems that when s approaches infinity (the evanescent region) the kernel depends only on s and source and receiver height. The layer-dependent parameters disappear and the calculation is essentially for an infinite medium. This means that the kernel may be approximated by the kernel corresponding to a simplified problem for large s values. In a similar but more intuitive manner, the problem may be simplified for s approaching zero. Small values of s represent near-vertical propagation contributions. For vertical contributions a horizontal sound speed gradient is without effect except for a phase shift as the contribution propagates through areas of different propagation velocity. The introduction of layers in the atmosphere does lead to reflection of sound when the layer boundaries are passed. These reflected contributions are, however, due to the discontinuity in the derivative of the sound speed gradient, and hence a second order effect which may safely be ignored. In fact they represent a weakness in the problem formulation rather than a true effect. The result is that the kernel may safely be approximated by the kernel for a ground homogenous air interface for small s values. The Thomson-Haskell approach is the usual one to apply to this class of problems. It is a matrizant method involving propagation of the field through all layers in the problem. This,leads to numerical problems for large values of the horizontal wavenumber. In the present approach (originating from that of Schmidt and Jensen [S] for non-gradient layers) the continuity at the layer boundaries is established by the linear equations system,
OUTDOOR
SOUND
77
PROPAGATION
which is arranged in such a manner that the most important contributions are close to the diagonal in the matrix equation. This means that the Gaussian elimination program removes the numerical stability problems.
3. NUMERICAL
RESULTS
The numerical tool presented in the previous sections may be applied to the simulation of various outdoor sound propagation conditions over plane ground. Possibilities of influence from extended reaction phenomena in connection with realistic ground models may be investigated, but that is not dealt with here. The main reason for developing the program was to investigate the theoretical importance of the presence of a wind speed gradient in the lower atmosphere. A model with a linear sound speed variation had been developed previously [I], but the present work makes it possible to establish the influence from the fact that the actual gradient is not adequately described by a simple linear velocity variation. Strictly speaking, the mathematical model deals with sound speed variations and not wind speed variations, but the effective sound speed may be viewed as the superposition of the physical sound speed and the wind speed component in the direction of sound propagation. In the following, results for a source height of 1.45 m and a receiver height of 0.5 m are shown for various sound speed conditions. The ground is specified by the Delany and Bazley relations, equations (23) and (24), for a flow resistance, a, given by 200x 10’ Ns/m4. The results are compared with measured one-third octave data from grasscovered ground previously published in references [l] and [ 111. The measured data were obtained for wind speeds of 2-2.5 m/s measured at a height of 10 m. Atmospheric absorption has been taken into account in the calculations by introducing a loss of O-002 dB per wavelength. Measured and calculated data are valid for a distance of 80 m except for that used in Figure 3 which is valid for 20 m. Figure 3 is included in order to demonstrate the deviation between theoretical and experimental values for a short distance where the wind gradient has only small influence. Calculations have been performed for the positive gradients shown in Figure 4 and for negative gradients having the same shape and size. One gradient is a simple linear variation
-20
-30
I
I
I
I 250
I 500
I 1000
-
125
Frequency
2000
(Hz)
Figure 3. Sound pressure level relative to free field over a typical grassy field. Horizontal distance 20m; source height 1.45 m; receiver height 0.5 m. --, Measured data, downwind [l]; - - -, measured data, upwind [I]; --, predicted for zero gradient.
78
K. R. RASMUSSEN
up to 20 m, and the others are bilinear up to 20 m. The bilinear cases simulate the meteorological fact that for a wind speed gradient and gradient is steeper close to the ground. All the gradients in Figure 4 are, however, terminated at a height of 20 m. This termination is without influence on the results but ensures that only a well defined range of sound velocities enter the calculations. Since the source must be in a layer with zero gradient, such a thin (O-1 m) layer is included in the calculations. This layer is not visible in Figure 4. An indication is given in Figure 5 that for a positive gradient some improvement of the similarity between measured and calculated data may be obtained by the inclusion of a bilinear velocity variation in comparison with the simple linear variation. It is indicated in Figure 6 that the bilinear cases overestimate the shadow-effect occurring at high frequencies due to the negative gradient. The size of this deviation is dependent
20-
343
344
3 IS
C (m/a)
Figure 4. Sound speed profiles used in calculations. 1, Simple gradient; 2, bilinear variation crossing simple gradient at 10 m; 3, bilinear variation crossing simple gradient at 20 m.
I
I
I
Frequency
(Hz1
Figure 5. Sound pressure level relative to free field for a positive velocity gradient over a typical grassy field. Horizontal distance 80 m; source height 1.45 m; receiver height 0.5 m. -, Measured data from reference [ 11; . . . . predicted for simple gradient: - - -, predicted for bilinear case no. 2; - --, predicted for bilinear case no. 31
OUTDOOR
%
SOUND
79
PROPAGATION
-10 -
-2o-
-30 125
250
500
Frequency (Hz) Figure 6. Sound pressure level relative to free field for a negative velocity gradient over a typical grassy field. Horizontal distance 80 m; source height 1.45 m; receiver height 0.5 m. -, Measured data from reference [ 11; . . I . . predicted for simple gradient; - - -, predicted for bilinear variation, case no. 2; --, predicted for bilinear Lariation, case no. 3.
upon the choice of parameters for the bilinear curve; the lower part of the sound speed variation is especially critical for high-frequency propagation. However, the deviation could also be severely influenced by the fact that turbulence effects are unaccounted for in the simulation. In fact, the shape of the high-frequency part of the measured curve has the approximate 20 dB per decade slope which is to be expected from turbulence considerations [ 111. It should be noted that the measured difference in level between the positive and negative gradient cases is more than 10 dB at 2 kHz. From a practical point of view, this large difference is well accounted for by the simple linear sound speed profile simulation. From a theoretical point of view the use of a linear profile is, of course, not
2-
OF
I
01
I
7,877
I
I
I
la.052
la.228
la.403
Horzontol
I
la.579
I
wovenumber
Figure 7. Pressure kernel for positive bilinear variation, case no. 3, for 1000 Hz. -, imaginary part.
Real part; --,
80
K. B. RASMUSSEN
-II 17.701
I I7.877
I I8.052
I
18~228
Hor~zontol
Figure 8. Pressure imaginary part.
kernel
for negative
bilinear
I 18.403
I 18.579
I 19-754
wavenumber
variation,
case no. 3, for 1000
Hz. -,
Realpart; --,
satisfactory. Ideally, a very accurate sound speed profile (consisting of a considerable number of line segments) should be used in connection with turbulence considerations. The numerical integration process is illustrated in Figures 7 and 8, showing the central part of the kernel for increasing and decreasing sound speed versus height at 1 kHz. As might be expected the kernel, which represents the wavenumber spectrum, has a more pronounced mode structure for increasing sound speed versus height than for decreasing sound speed versus height. This is a consequence of the concentration of sound energy along the surface for a positive gradient. On the basis of the material in this paper it seems that a simple linear sound speed profile could be a reasonable representative of the combined effect of turbulence and sound speed gradients in the atmosphere. Further measured data representing point source propagation under well defined meteorological conditions are needed in order to test this preliminary conclusion. Furthermore, indications that an improved prediction may be obtained for a more realistic sound speed profile are present, but such a more realistic profile makes it more urgent to introduce turbulence in the calculations.
REFERENCES 1. K. 9. RASMUSSEN 1986 Journal ofSound and Vibration 104,321-335. Outdoor sound propagation under the influence of wind and temperature gradients. KUESTER,D.C.CHANG,W.F.RICHARDS,R.GILBERT~~~ N. 2. R.RAsPET,S.W.LEE,E. BONG 1985 Journal of the Acoustical Society of America 77, 345-352. A fast field program for sound propagation in a layered atmosphere above an impedance ground. 3. S. W. LEE, N. BONG, W. F. RICHARDSand R. RASPET1986 Journal of the Acoustical Society of America 79, 628-634. Impedance formulation of the fast field program for acoustic wave propagation in the atmosphere. 4. T. L. RICHARDSand K. AI-I-ENBOROUGH1986 Journal ofSound and Vibration 109,157-167. Accurate FFT-based Hankel transforms for predictions of outdoor sound propagation. 5. H. SCHMIDT and F. 9. JENSEN 1985 Journal ofrhe Acoustical Society ofAmerica 77,813-825. A full wave solution for propagation in multilayered viscoelastic media with application to Gaussian beam reflection at fluid-solid interfaces. and 9. VIS-~ISEN1986 Repot-f86.149 Qdegaard & Danneskiold-Samstie 6. K. 9. RASMUSSEN Kroghsgade 1 DK-2100, Copenhagen, Denmark. Calculation of sound propagation under water by means of FIT-model. (In Danish.) I. A. D. PIERCE 1981 Acoustics: An Introduction to its Physical Principles and Applications. New York: McGraw-Hill. 8. V. A. FOCK 1965 Electromagnetic Diflraction and Propagation Problems. Oxford: Pergamon Press.
OUTDOOR
SOUND
81
PROPAGATION
and E. N. BAZLEY 1970 Applied Acoustics 3, 105-116.Acoustical properties of fibrous absorbent materials. 10. M. ABRAMOWITZ and I. A. STEGUN 1970 Handbook cf Mathematical Functions. New York: Dover. 11. K. B. RASMUSSEN 1985 Danish Acoustical Institute Report. The effect of wind and temperature gradients on sound propagation outdoors.
9. M. E. DELANY
.APPENDIX:
ON
THE
CALCULATION OF AIRY DERIVATIVES
FUNCTIONS
AND
THEIR
The W and V functions defined in terms of the Airy function in equations (18) and (19) are complex functions of complex argument. They may be evaluated on the basis of formulas from the book by Abramowitz and Stegun [lo] as follows:
Ai
= c,.f(z) - c,g(z), cZ = - Ai’(0) = 0.258819,
c, = Ai(0) = 0.355028, f(z)
=
(Al)
Z0 1+“‘ + z6 2.3 2.3.5.6+2.3.5.6.8.9 Z
+...,
V and
W are then determined
(A3)
IO
+. . .
g(z)=z+L+z:
3.4
(A2)
(A4)
3.4.6.7+3.4.6.7.9.10
from W(z) = 2.,.r’/’ e’njh Ai (z e’2n/‘).
V(z) = T”’ Ai (z),
(A5, A6)
The derivatives V’ and W’ are calculated from the similar expressions obtained when equations (A3) and (A4) are differentiated. Equation (Al) is slow in computation due to the series (A3) and (A4). For large positive and negative arguments (real part) asymptotic expressions may be used [8, lo]. For Re (z) >>1 one obtains (including three terms) U(z) = z -“4 e”( 1 + a,/x
+ u,/x” + a,/~~),
U’(z) = z’14 e”( 1 - b,/x - b,/x’-
(A7)
b,/x3),
(A8)
V( 2) = o*5z-“4
e-“( 1 - a,/x
+ a,/.~‘-
a,/~“),
(A9)
V’(z) = -0.5~“~
e-‘( 1+ b/x
- b,/x’+
b,/x3),
(AlO)
where x = (2/3)z”* and a, = O-069444, a, = O-037133, a3 = 0.037993, b, = O-097222, b2 = 0.043885 and b, = 0.042463. W and W’ are determined from equations (A7)-(AlO) with W(z)=
For negative U(-z) U’( -z)
values
one similarly
U(z)+iV(z).
obtains
= zz”4 cos (x + 7r/4)(1-
(using
(All) the same definitions
UJX’) + zz”J sin (x+ 7r/4)(a,/x
= z”4 sin (x + 7r/4)( 1 + b2/x2) + z”’ cos (x + r/4)(
V(-z)
= zP1’4 sin (x + r/4)(
V’( -z)
= - z’/4 cos (~+~/4)(1+b~/x’)+z”~sin
of CJ and x) - a,/~‘),
(A12)
b,/x - b,/x3),
(A13)
1 - a2/x2) - 2-“4 cos (x + ~/4)(a,/x
- a,/~~),
(x+n/4)(b,/x-b,/x’).
(A14) (A15)