THE UN IVER S I TY OF MICHIGAN COLLEGE OF ENGINEERING Department of Nuclear Engineering Laboratory for Fluid Flow and Heat Transfer Phenomena Technical Report GEOMETRICAL ACOUSTIC PROPAGATION WITH FREQUENCY DEPENDENCE L. P. Solomon ORA Project 01357 supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. GK-1889 WASHINGTON, D. C. administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR April 1968

ABSTRACT One characteristic of geometrical acoustic solutions is their frequency independence. Geometrical acoustics stems from the solution of the eikonal equation which is a form of the wave equation in the limit of infinite frequency. Since the effects of the medium and the boundary conditions on the acoustic fields are definitely frequency dependent, it is desirable to obtain approximate solutions that reflect this dependence on frequency. Although data are available on propagation losses for signals of different frequencies, geometrical acoustics does not offer a means of comparison. An asymptotic expansion in terms of inverse powers of the frequency is presented. The first term in the expansion obeys the eikonal equation; all the ray-path solutions obtained from geometrical acoustics are applicable. Here, the first term of the expansion, which is the geometrical acoustic field, is obtained in closed form using a simplified profile in which the first derivative is continuous everywhere. The calculation of the higher-order terms then involve operators that have known variable coefficients which are functions of the lower-order solutions. ii

NOMENCLATURE a = constant s b = constant s c(z) = speed of sound c = reference speed of sound o c = minimum speed of sound f,g = functions - see Eq. (24) i=r k = wave number = w/c o k = constant s L,L = differential operators ac 1/2 m(y) = bv /f(y) n(z) = c0/C(Z) p = constant p,q,rs,uw = summation indices P,Q,R = functions - see Eq, (33) t = arbitrary integer x,y,z = spatial coordinates x,z = translated coordinates v = amplitude function v = amplitude functions P iii

a= constant s = parameter 5 = profile parameter o t f.,. = arbitrary constants J J 7 = speed of sound profile constant V = del operator = Laplacian operator 0 = angle of ray path to horizontal i = x - x s T = time c = phase function X amplitude function = solution of reduced wave equation 4 = solution of wave equation Q= slope of speed of sound profile = - radian frequency () = source condition iv

1. INTRODUCTION In a preceding paper [herein referred to as (I)] it was pointed out that one property of all ray path solutions obtained using geometrical acoustics is that they are independent of frequency. Since geometrical acoustics is actually a solution of the wave equation in the limit of infinite frequency, it seems that one method for obtaining frequency dependent solutions would be to obtain an asymptotic expansion in terms of inverse powers of the frequency. It was further pointed out in (I), that one of the practical difficulties in dealing with asymptotic expansions is the calculation of higher order terms. The difficulty in calculating the higher order terms is directly related to the complexity of the leading term. This observation about asymptotic expansions in general prompted an investigation to find simple expressions for the first term in the asymptotic expansion of the wave equation. This first term, being Frequency independent, would be the geometrical acoustic solution (i.e., the ray paths); closed form expressions for the ray paths are given in (I). This paper presents a general asymptotic expansion for the wave equation, 2 following the method given by Friedlander and Keller. The form of the expansion allows for the possibility of exponential decay factors, for inverse fractional powers of the wave number k, and for dispersive effects. The ex3 4 pansions of Luneberg and Kline, which are a special case of the expansions given here, do not include these effects. 1

2. GENERAL THEORY Let V be some function which satisfies the wave equation, 1 = o (1) C aT where T is the time, c is the speed of propagation and will be considered to be a function of depth z only, and V2 is the Laplacian operator. Further, assume that $ = l e (2) where w is the frequency in radians/second, and where' is a function of the spatial variables only. The substitution of Eq. (2) into Eq. (1) gives rise to the reduced wave equation ~2 +?2 =* = 0 (5) c (z) Now let ik,-k a(x,y,z) 4(x,y,z) = v(x,y,z)ekk (4) 2

where k = -- (5) c 0 Substitution of Eq. (4) into Eq. (3) gives k2v[n2 - (V)2 + k v[-2i Vv.Vx + k[2i Vv.V0 + iv V2Q] + k [-2Vv.X - vV2x] v +k )2 +2v = 0 (6) where 2 2 0 n - = (7) o -(z) Finally let v(x,y,z) = 0 k/ v (x,y,z) (8) p=O p where t is an integer to be determined. It is positive, and dependent upon a. Substitution of Eq. (8) into Eq. (6) gives k2-/t v [n2 - (Vk)2] + Z kl-/t v [-2i V.VX] p=00 p q=0 q 00 1-r/t 00 -s/t + Z k / [2i Vv V + iv V20] + Z k [-2W VX - v V2] r=0 r r s=0 s s 00 k2au/t 00 W/0 + k2 -u/ v (vx)2 + k- vv v = (9) u=O u w=O w Certain comments may now be made about the range of a. If a is negative or zero, then the exponential decay which was anticipated for large values of k will not be obtained. In fact, as Friedlander and Keller point out, if oa is negative or zero, then the Taylor expansion of exp(-k X} in powers of k X 5

proceeds in nonincreasing powers of k. If this factor is expanded, and the resulting series multiplied by the expansion given in Eq. (8), then the product so obtained will be of the same form. This suggests that the v may be redeP fined, and the factor exp(-k X] merely be eliminated from Eq. (4). This expansion, as well as for the case of a = 1 coincides with the Luneberg-Kline case. Thus for new forms of the expansion, only positive values of a will be considered. Let us now group terms in Eq. (9). This requires that all exponents of k be equal, i.e., q = p + t( - 1) r = p - t s = p + t(a 2) u = p + 2t(a - 1) w = p- 2 (10) The value of t will, of course, depend upon the value of a. One will always be able to find a t such that t(a - 1) is an integer if a is assumed to be a rational number, i.e., a ratio of two integers. Note that if t(a - 1) is an integer, 2t(a - 1) will also be an integer, as well as t(a - 2) = t(a - 1) - t. With the assumption that ais a rational positive number, and using Eq. (1), Eq. (9) may then be rewritten: 4

po k2-p v[ - (vt) ] + v [t-2i V1.7X] p=O p+t(a-l) 12i V pt + ivt 2] + [-2 v'p+t(a-l) + +t(a-2) p+2t(a-l)(VX + p-2t = (1) Let us now consider various values of a. Case 1: a> 1 For p = O, it is required that the coefficients of v ) ----------- 9 t(a-l)' and v 2t ) be zero. 2t(a-i) V. Vx = (vX)2 = 0 This requires that X = constant. The next term will require that (V)2 - n2 = 0 which is the eikonal equation, and corresponds to geometrical acoustics. Let La - 2V.V + V2 Then for successive values of p the following argument may be made: p =t L v 0 a o p= t +1 L vl = 0 p = 2t - 1 L v =0 a t-1 5

This implies that v =. v j = O0, 1,... t - 1 and P = 1. J J o o p 2t L v - 2v a t o p = 2t + 1 L v - vl t+L p =- t- 1 L - 2V a 2t-1 t-l t t This implies, as above, that vt+j =P j vt j =, 1,... t - 1 and P =1. t Any v. will be repeated t times, with arbitrary constants.. Thus in principle all v may be calculated. The important result is that for a > 1, X must be a constant, and 0 satisfies the eikonal equation. Case 2: a = 1 Using exactly the same procedure as in Case 1, it can be shown that the conditions on X and $ are (V0)2 - n2 = (VX)2 VO x = 0 Equations (14) and (15) are satisfied by the condition that X = constant. This corresponds to the Luneberg-Kline expansion mentioned above. Furthermore, if it is assumed that X = constant in a solution, then the operators necessary to calculate v. are identical to those of Case 1. a Case 5: 1/2 < al < 1 Proceeding as usual, setting the coefficient of the first term (p = 0) equal to zero gives the eikonal equation, (V0)2 - n2 = 0 6

The next term which gives a contribution requires that V.Vx = 0 The next term, however, requires that (vx)2 0 Again, only X = constant will satisfy this condition. Therefore, this case reduces to the two previous cases and generates the same v;. J Case 4: a = 12 Using the same procedure as before, the first term (p = 0) again yields the eikonal equation. It is clear that t may be any even number. The next term in the expansion, corresponding to p = t/z, requires that V7Vx- = 0 Let us define L/ = 2i v.v + i + (Vx)2 1/2 Then for the next term, i.e., p = t, p = t LL1/2 = 0 p = t + 1 L1/2 V1 0 p = 3t/2 L1/2 vt/2 = 2VOV - (=X )v p = 2t L12 v = 2VX.t/2 - (X)vt/2 - v 7

This case does not require that X = constant; only that surfaces of constant < and surfaces of constant X be orthogonal to each other. Case 5: a < 1/2 Again, as in Case 4, the first term (p = 0) requires that the eikonal equation be satisfied. The next term, when p = t(l - o), will again require that surfaces of constant D be orthogonal to surfaces of constant X. The next term will be when p = t. Recalling the definition of L from Case 1, we may write p =t L v = 0 Oa o p = t + 1 L v = O p = 2t(l - O) L t(La = (X)2 p = t(2 - a) L Vt( = -(VX)2v + 2Vx.Vv - (Vx)v p = 2t L vt =-(VX)2v + 2VX VV - (X)va - etc. Case 6: c = 0 As usual, the first term (p = 0) gives the condition that ~ must satisfy the eikonal equation. The next term will occur when p = t. This gives as a condition 8

- 2(V0.Vx)v + 2V v -V0 + v V2 = 0 o 0 0 One should note that there is no condition on X. If it is assumed to be a constant, then the above condition becomes (2V -V + V)v 0 0 This is the operator L. The choice of X = constant is acceptable without loss of generality for the case of a = 0 because the solution was assumed in the form ikO-k X v(x,y,z)e If a = O, this reduces to ikl-X ikl v(x,y,z)e = v(x,y,z)e Thus, for a = 0, it is merely necessary to carry the same procedure through as before, after redefining the v.. The condition on < and X can then be summarized: 1/2 < a X = constant (V0)2 n2 = 0 O < a 1/2 VO.Vx (V) )2- n2 = 0 = 0 X = constant (V~)2 - n2 = 0 9

3. RELATIONSHIP BETWEEN THE EIKONAL EQUATION AND GEOMETRICAL ACOUSTICS The eikonal equation for 0 is (VD)2 - n2 = (12) The surfaces 0 = constant are the wave surfaces or wave fronts. The equation for the normals to these surfaces is dx dy dz ( S0/-x - a/ay 0/z The relationships between the normals to the surfaces and the index of refraction is well known.' For most cases of physical interest in underwater sound propagation, the speed of propagation-and therefore the index of refraction n-is a function of the depth only. This assumption requires that any path in (x,y,z) space lie in a plane normal to the x,y plane, and be fixed in space. Let us assume that this is the x - z plane. Figure 1 illustrates the appropriate coordinate system. It is clear, therefore, that if n(z) is specified, and the paths are known so that G(x,z) is known at each point along the path, then 0 may also be calculated by the following expressions -x n(z) cos 9(x,z) ax - = n(z) sin O(x,z) (14) az If the path is given by z = z(x), then Eq. (14) may be rewritten as do n(z) dz sin 9(z) 10

z - QI2' S Figure 1. Raypath coordinate system. Figure 1. Ray path coordinate system.

4. UTILIZATION OF THE EXPANSION TECHNIQUE IN SOME EXAMPLES Several examples have been calculated using the speed of sound variation, c(z), employed in (I). (i) c(z) = constant c o n(z) = 1 (16) Then from Eq. (14) = x cos 0 + z sin 0 (17) s s For all a > 1/2 and a = 0, X = constant, and the solution for v is u o v = v + v [x sin - z cos ] o 00 01 s S However, since X = constant, this would allow the amplitude to be unbounded as x or z is increased without limit. Thus v =0 and v = v = constant. 01 0 00 The same argument holds for a = 0. Since v constant, Vv = 0, this implies o o that all v. are constant. Therefore, J + 1 + k +.. exp k(x cos e + z sin s) - k o/t k2/t s s (18) For 0 < a < 1/2, X must satisfy V71Vx = 0 12

The solution for X can easily be found to be a constant. For a = 1/2 it is found that v. are all constants. The same result is found for 0 < < 1/2. j Thus, for all a > 0 the solution for X is given by Eq. (18). (ii) c(z) = Qz + 7 c n(z) = + (19) Oz + y Now the ray path is given by x + z2 = p2 (20) where (1 - A2 2)1/2 x = x- x + (21a).s AQ = z + (21b) c sin e s s sz x = (21c) s Q s As - 1 p A= (21d) cos E A == -s (21e) c s Utilization of Eq. (14) allows the calculation of c px c - = ~_ + 0- -in _ - + x (22) 2 2 - s 22z z where c sin e c 1 + sin 0 o s o ___1 s s ~ cos2 0 2Q cos s s 13

Again let us consider the calculation for X, which are of course dependent upon a. For all 1/2 < a and ac = 0, X = X = constant. The equation for v is (2V-V + V')v = 0 0 Performing the indicated operations, and realizing that along the path, z = z(x), v may be expressed as o v = v exp- (2 (23) where f(y) = - 2 y 2 - + 4p(p2 y2)1 (24a) g(y) = 0402 + 4p(p2 -y2)1/2 -2 (24b) and where y is a dummy variable, After considerable algebraic manipulation, it can be shown that Eq. (23) becomes v = v - (C(2(P - q) - p)/(2(p _ q) + p))2/ (25) v v z (25) O o (p2 + 4pq 2 / where -2 2 -2 q2 = p -z = x The calculation for vl follows directly from the general expansion procedure given in Section 3. As expected, vl may be expressed as z vl = const. v - v / rny) dy (26) o o 14

where V2% 0 m(y) =f() (iii () c(z) = c cosh (z - m) c 0 n(z) - cosh (z - ) (27) The ray path is given by 2 ~(l_ -)1/2 sin 3(x - x ) + cos (x - x ) (28) a s a s a s s s s where 5 = sinh 5(Z - zm) a = 1 - c) 1/2/a s s s cos e (29) s cosh B(Z - (29) and where for 0 > 0, the + sign is taken, and for 9 < O, the - sign is taken. s s Using Eq. (14) as before, the solution can be attained: k s b + (a 22)1 a + (a 1/2 s- a s + s s = k x + log - 2a log S 4P s - (2)1/2 s s S where b2 1 + a2 s S C /5 k = c/(31) (1+ a 15

For a > 1/2 and a = 0, X = X = constant. 0 The equation which v satisfies is, as before, O (2v V+ V2q)v o 0 Performing the indicated operation, and using the fact that g = t(x) along the path, v may be expressed as V f Q, (y) dy ln _ R(y) + P() (32 v 0 where (a - 2)1/2 P(O) = s(1 + 2)1/2 fs(l + a2) Q(S) = 1/2 2(1 + r2)(a - t2)1/2 R() = +f(a2 - s2)/ cos P(x - x ) - PS sin (x - x ) (33) and where S + S 1/2 ( 2 1/2 cos P(x - x) = ) (1 ) (34) s a a a s s s s2 as s 1/2 1/ sin P(x - x ) + (1 - )1/2(1 ) 1 2 s + =a -aa a a a- a s s S S s S It can be seen that due to the algebraic expression given by Eq. (34) which we required for R(S) given in Eq. (33), the task of integrating Eq. (32) exactly is extremely difficult. Nevertheless, it may be evaluated using numerical techniques. 16

The next term in the expansion, vl, can be obtained by the methods described above. Although the calculation is tedious, it poses no analytical problems. 17

BIBLIOGRAPHY 1. Solomon, L. P., Ai, D. K., and Haven, G., "Acoustic Propagation in a Continuously Refracting Medium," J. Acoustical Society of America, to be published. 2. Friedlander, F. G., and Keller, J. B., "Asymptotic Expansions of Solutions of (V2 + k2)u = 0," Communications on Pure and Applied Mathematics, Vol. III, pp. 387-394, (1955). 3. Luneberg, R. K., "The Mathematical Theory of Optics," Lecture Notes, Brown University, 1944, and "Propagation of Electromagnetic Waves," Lecture Notes, Brown University, 1948. 4. Kline, M., "An Asymptotic Solution of Maxwell's Equations," Comm. Pure and Apple Math., Vol. 4, No. 2-3, pp. 225-263. 5. Officer, C. B., "Introduction to the Theory of Sound Transmission," McGrawHill, 1958, New York, pp. 41f. 18