Derivation and implementation of the boundary integral formula for the convective acoustic wave equation in time domain Yong Woo Lee and Duck Joo Leea) Division of Aerospace Engineering, KAIST, Daejeon, Republic of Korea
(Received 3 June 2014; revised 21 September 2014; accepted 6 October 2014) Kirchhoff’s formula for the convective wave equation is derived using the generalized function theory. The generalized convective wave equation for a stationary surface is obtained, and the integral formulation, the convective Kirchhoff’s formula, is derived. The formula has a similar form to the classical Kirchhoff’s formula, but an additional term appears due to a moving medium effect. For convenience, the additional term is manipulated to a final form as the classical Kirchhoff’s formula. The frequency domain boundary integral can be obtained from the current time domain boundary integral form. The derived formula is verified by comparison with the analytic solution of source in the uniform flow. The formula is also utilized as a boundary integral equation. Time domain boundary element method (BEM) analysis using the boundary integral equation is conducted, and the results show good agreement with the analytical solution. The formula derived here can be useful for sound radiation and scattering by arbitrary C 2014 Acoustical Society of America. bodies in a moving medium in the time domain. V [http://dx.doi.org/10.1121/1.4898427] PACS number(s): 43.20.Rz, 43.20.Bi [JDM] I. INTRODUCTION
The boundary integral formula is useful for sound radiation problems. The well-known Helmholtz integral formula is a boundary integral formula for the Helmholtz equation in frequency domain. The acoustic property at the observer position can be obtained through boundary integration. Another boundary integral formula is the Kirchhoff’s formula, which is the formula for wave equation in the time domain. The two formulas in the preceding text are based on the wave equation with stationary medium. For moving medium, Blokhintzev1 derived an integral formula for the convective wave equation using coordinate transformation. Blokhintzev transformed convective wave equation using a new system of coordinates and derived the final form using Green’s theorem in time domain. Blokhintzev’s integral formula contains surface integration in transformed coordinates. Wu et al.2 derived the convective integral formula in the frequency domain using the weighted residual integration and Green’s theorem. The integration of the formula of Wu et al. is in physical coordinates, and it is convenient to use. Guo3 also obtained an integral formula in the frequency domain using weighted residual integration and Green’s theorem. Another method for the derivation of the boundary integral formula is the generalized function theory.4 Ffowcs-Williams and Hawking5 extended Lighthills6 acoustic analogy theory for a moving body using this generalized function theory. The generalized function theory makes it possible to handle discontinuous areas, such as body surface, and it simplifies the derivation process. Recently, the convective acoustic analogy formula has been derived7 using the generalized function theory. Farassat et al. utilized the generalized function theory for the derivation a)
Author to whom correspondence should be addressed. Electronic mail: [email protected]
J. Acoust. Soc. Am. 136 (6), December 2014
Pages: 2959–2967 of Kirchhoff’s formula.8,9 Farassat et al. set a moving surface surrounding the moving acoustic sources and used the generalized potential, which is discontinuous at the surface. Kirchhoff’s formula offers an advantage in which it can be utilized as a boundary integral equation itself to solve the radiation and scattering problem for the moving surface.10,11 Thus the scattered field can be obtained by solving the boundary integral equation in the time domain. The time domain boundary element method (BEM) analysis has stability problems.12,13 An averaging scheme is used to suppress the instability, but it increases the computational cost.14–16 The high-order approximation for space and time is also used to reduce the instability problem.17–19 For the exterior problem, the fictitious internal resonance in the body also can cause the instability problem.20 To solve this internal resonance problem, the time domain Burton–Miller formulation is introduced.21 But the Burton–Miller formulation requires an additional equation for the normal derivative of the variable and it increases computational cost. Recently, the vector filtering technique was introduced to suppress the instability by filtering the instability wave, and it shows good results.22 In this paper, we derived the Kirchhoff’s formula for the convective wave equation in physical coordinates. The obtained formula contains an additional term compared to the classical Kirchhoff’s formula due to a moving medium effect. The additional term was modified to have the same form as in the classical Kirchhoff’s formula for easy implementation. The generalized function theory was employed for the derivation of the formula. The derived formula is verified by comparison with an analytic solution of simple source. In addition, time domain BEM analysis is conducted to obtain the utility of the formula as a boundary integral equation, which can handle the radiation and scattering problems for arbitrary bodies in a moving medium. We used a simple implicit time marching scheme to avoid the longtime instability problem.23 If the frequency of the acoustic
0001-4966/2014/136(6)/2959/9/$30.00
C 2014 Acoustical Society of America V
2959
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
wave does not match with the internal resonance frequency, the instability due to the internal resonance does not occur. Thus the implicit time marching scheme is appropriate to see the utility of the formula. The numerical technique for the convective time domain BEM using implicit time marching scheme is demonstrated, and a simple test case is solved using the time domain BEM technique without instability problem.
ðfÞ @H @f @f ¼ H0 ð f Þ ¼ dð f Þ ; @x @x @x
II. DERIVATION OF GENERALIZED CONVECTIVE ACOUSTIC WAVE EQUATION
To obtain Kirchhoff’s formula for the convective acoustic wave equation, the generalized convective wave equation was derived first. For simplification, we assumed the medium moves in the positive x direction with constant speed U0 . Then the convective wave equation can be written as 1 @2/ U0 @ 2 / U2 @ 2 / þ2 2 þ 20 2 r2x / ¼ 0; 2 2 a0 @t a0 @t@x1 a0 @x1
(1)
where / is the velocity potential, a0 is the speed of sound, x1 is the component of the position vector ~ x ¼ ðx1 ; x2 ; x3 Þ, and x . If rx denotes gradient for the observer position vector ~ there are bodies, Eq. (1) is valid in the region outside the body surface. The body surface can be any kind of surface, including a rigid and permeable surface. The surface can be represented by f ð~ x ; tÞ ¼ 0 such that f ð~ x ; tÞ < 0 for the region inside the surface and f ð~ x ; tÞ > 0 for the region outside the surface. Function f has set as jrx f j ¼ 1, which simplifies the derivation process.9 Consequently, rx f is the surface normal unit vector n^ pointing out from the surface. We assumed that Kirchhoff’s surface is stationary so f is a function of ~ x only (Fig. 1). Then / can be extended inside the surface by assuming that / ¼ 0 inside the surface using the function f. The generalized potential is written using the Heaviside step function as9 ~ ¼ Hðf Þ / ¼ /
/ for f > 0 0 for f < 0:
Notice that the generalized variable is discontinuous at f ¼ 0, and it is not differentiable. However, the generalized function theory gives derivatives to this discontinuous variable, the so-called the generalized derivative.8,9 The Heaviside step function is the generalized function and its derivative is the Dirac-delta function. Thus the generalized derivative of the Heaviside step function can be written as
(2)
(3)
where the overbar denotes the generalized differentiation to distinguish from the ordinary differentiation. By utilizing the above generalized differentiation, the convective wave equation can be changed to the generalized convective wave equation that governs the whole region. By using the generalized derivatives, the generalized convective wave equation can be rewritten as follows: 2~ 2~ 2~ 1 @ / U0 @ / U02 @ / ~ 2/ þ 2 þ r x a20 @t2 a20 @t@x1 a20 @x21 U0 @/ U 2 @/ n1 dð f Þ þ 20 n1 dð f Þ ¼2 2 a0 @t a0 @x1 @ þ ½/n1 dð f Þ @x1
@/ x ½/^ dð f Þ r n dð f Þ: @n
(4)
Equation (4) is the convective acoustic wave equation governing the whole region. Equation (4) has source terms on the right-hand side and they are the surface sources. III. DERIVATION OF KIRCHHOFF’S FORMULA FOR THE CONVECTIVE WAVE EQUATION
The generalized convective wave equation has source terms on the right-hand side, and the solution can be obtained by using the Greens function. Greens function for the convective wave equation is G¼
dðs t þ R=a0 Þ ; 4pR
(5)
where R¼
M0 ðx1 y1 Þ þ R
;
(6)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ ðx1 y1 Þ2 þ b2 ½ðx2 y2 Þ2 þ ðx3 y3 Þ2 ;
(7)
b2
b¼
FIG. 1. The schematic of a stationary surface in a moving medium. 2960
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 M02 :
(8)
Here s is source time, t is observer time, x is observer position, y is source position, and M0 ¼ U0 =a0 . The difference of the Greens function for the convective equation is that the radius is changed to the amplitude radius R and the phase radius R. The solution of Eq. (4) can be obtained using the Greens function in the preceding text as Y. W. Lee and D. J. Lee: Convective boundary integral formulation
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
ð ð ðt ð dð gÞ dð gÞ dðgÞ @/ U02 t @/ n1 dð f Þ dVds þ 2 n1 dð f Þ dVds rx / n^dð f Þ dVds 3 3 3 4pR 4pR 4pR @s @y a 1 1 R 1 R 0 1 R ðt ð ð ð dðgÞ dðgÞ U 2 @ t x /n1 dð f Þ dVds r /^ n dð f Þ dVds; þ 20 4pR 4pR 3 a0 @x1 1 R3 1 R
2U0 ~ ð~ / x ; tÞ ¼ 2 a0
ðt
ð
(9)
where g ¼ s t þ R=a0 . All integrals of Eq. (9) have Dirac-delta functions, and these restrict the integration range to zero of the function variable. Thus the integration for space reduces to the surface, f ¼ 0, and the integration for time reduces to retarded time which is the solution of g ¼ 0, sr ¼ t
R : a0
(10)
This means that the integrand has to be evaluated at retarded time sr for desired observer time. It should be noted that the retarded time is defined with the phase radius R in Eq. (6), which is not the geometric radius. When the medium moves, the propagation time R=a0 of the acoustic wave is different from that of a stationary medium. The last two integrals have spatial derivatives outside the integral and the spatial derivative can be changed to a temporal derivative using the property of the Dirac-delta function.7 Hence, Eq. (9) can be changed to ð ð ð @/ 1 U02 @/ 1 @/ 1 U02 @ R~1 n1 dS þ 2 n1 dS dS 3 /n1 dS 4pR ret 4pR ret a0 f ¼0 @y1 4pR ret a0 @t f ¼0 f ¼0 @s f ¼0 @n 4pR ret ð ð ð U02 1 @ rx R rx R R~1 /n1 dSþ /^ n dSþ /^ n dVds; 2 a0 @t f ¼0 4pR ret a0 f ¼0 4pR2 ret 4pR2 ret f ¼0
U0 ~ ð~ / x ; tÞ ¼ 2 2 a0
ð
where R~1 ¼ @R=@x1 , and R~1 ¼ @R =@x1 . The bracket ½ ret means that the integrands have to be evaluated at retarded time sr . If Eq. (11) is used directly, we need to differentiate the integral at the observer position with respect to the time due to the temporal derivative outside the integral. It is convenient to use if we move the temporal derivative from outside the integral to inside the integral by utilizing the relation between observer time and retarded time sr t þ R=a0 ¼ 0. Because R and a0 are constant, the observer time derivative can be replaced by the source time derivative directly. Then Eq. (11) can be arranged as ð @/ @/ @/ ~ þ B/ þ C þD A n1 dS; 4p/ ð~ x ; tÞ ¼ @s @n @y1 S ret (12)
similar form to the classical Kirchhoff’s formula except for the coefficients and @/=@y1 term. Equation (12) requires additional information at surface @/=@y1 , which does not appear in the classical Kirchhoff’s formula. If we can obtain @/=@y1 at the surface, Eq. (12) can be employed directly, but it usually requires additional computation. Wu et al.2 and Guo3 manipulated this term in the frequency domain through decomposition of gradient to the surface normal and tangential component; this resulted in @/=@y1 changing to a combination of @/=@n and /. Using the preceding approach in this time domain, the second form of Kirchhof’s formula for the convective wave equation was obtained, ð @/ @/ ~ ð~ þ B0 / þ C0 A0 dS; (17) 4p/ x ; tÞ ¼ @s @n ret S where
where 2M0 n1 M02 n1 R~1 þ n^ rx R ; A¼ a0 R B¼
(11)
M02 n1 R~1 þ n^ rx R ; R2
1 C¼ ; R M2 D ¼ 0 : R
(13) (14)
(15) (16)
Equation (12) is Kirchhoff’s formula for the convective acoustic wave equation in physical coordinates, and it has a J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
2M0 n1 M02 n1 R~1 þ n^ rx R ; a0 R M02 n1 R~1 þ n^ rx R Mn ~ 0 Ms ; rs B ¼ R R2 A0 ¼
C0 ¼
1 Mn2 ; R
(18) (19) (20)
~0 n^, M ~s ¼ M ~0 Mn n^, and rs ¼ r @=@n. where Mn ¼ M Equation (17) now has the same form as the classical Kirchhoff’s formula. Thus Eq. (17) can be utilized in the same way as the classical Kirchhoff’s formula with the coefficients in Eqs. (18) to (20). Y. W. Lee and D. J. Lee: Convective boundary integral formulation
2961
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
IV. BOUNDARY INTEGRAL EQUATION FORM OF THE CONVECTIVE KIRCHHOFF’S FORMULA
If the observer is on the surface, the formula changes to a boundary integral equation. Equation (17) is valid for the observer position in the region outside the surface. To obtain the boundary integral equation, we consider that the observer approaches the surface from the region outside the surface. The domain is augmented by a hemi-ellipsoid, and its semiprincipal axes are in the positive x1 ; x2 ; x3 direction. The length of semi-principal axes is =b in the positive x2 and x3 direction and in the positive x1 direction. If tends to zero, the observer approaches the surface, and the boundary integral equation will be obtained. Kirchhoff’s formula for the augmented domain is ð 0 @/ 0 0 @/ þB/þC A dS 4p/ð~ x ; tÞ ¼ @s @n ret S ð @/ @/ þ B0 / þ C0 þ A0 dS: (21) @s @n ret S We consider the region outside the surface and tilde over the velocity potential is removed. Because tends to zero, it is clear that all terms in the second integral converges to zero except the / term. The second integral of Eq. (21) can be written as ð M02 n1 R~1 þ n^ rx R / dS: (22) R2 ret S To evaluate Eq. (22) more conveniently, we apply the Prandtl–Glauert transformation x1 ¼ x1 ; x2 ¼ bx2 ; x3 ¼ bx3 :
(23)
By applying the Prandtl–Glauert transformation, the hemiellipsoid changes to a hemi-sphere with radius . Equation (22) then changes to ð n rx R / dS; (24) 2 S R ret where n ¼ ry f is the surface normal vector in transformed x 2 y2 Þ2 þ ð x 3 y3 Þ2 g0:5 domain, and R ¼ fð x 1 y1 Þ2 þ ð is the distance between the observer and source in transformed domain. If tends to zero, Eq. (24) will be 2p/ð~ x ; tÞ. Thus the boundary integral formula can be arranged as ð 0 @/ 0 0 @/ ð Þ þB/þC A dS; (25) C~ x /ð~ x ; tÞ ¼ @s @n ret S where Cð~ x Þ is
8 < 4p region outside the surface Cð~ x Þ ¼ 2p on the smooth surface : 0 region inside the surface:
(26)
V. IMPLEMENTATION OF THE KIRCHHOFF CONVECTIVE BOUNDARY INTEGRAL EQUATION
The formula has been changed into the form of the classical Kirchhoff’s formula, but the surface gradient term of Eq. (19) is still not easy to handle. Guo3 solved this problem for constant element in the frequency domain BEM. Because each element of Kirchhoff’s surface is constant, / can be placed outside the integral, and the integrand will be a surface gradient term only. This is utilized in the time domain problem. This can be changed to a line integral by divergence theorem, ð Mn ~ Mn ~ rs M s dS ¼ M s n~dl; R S l R
ð
(27)
where n~ is a line normal unit vector pointing outward from the element and parallel to the element. After discretizing the surface into N elements, Eq. (25) can be written in discretized form as Cð~ x Þ /ð~ x ; tÞ ¼
N X j¼1
@/j @/j þ Bj 0 /j þ Cj 0 Aj @s @n 0
DSj : ret
(28) If we know the information at the surface, Eq. (28) can be used easily to obtain the velocity potential at the observer position. However, if we do not know the information at the surface, it can be obtained using Eq. (28) by positioning the observer on the surface with the time domain BEM technique. Generally, the surface normal derivative term is given as a boundary condition, and the unknown is the velocity potential. To obtain the velocity potential at the surface using Eq. (28), the time has to be discretized too. The velocity potential at a certain time step can be obtained by summing up the velocity potential and its spatial and temporal derivative terms at the corresponding retarded time. But generally the retarded time does not coincide with the discretized time steps although the observer time step is discretized because the propagation time R=a0 is not an integer order. Thus the values at the retarded time, located between two integer time steps, have to be evaluated using interpolation with the corresponding integer time steps. An interpolation method using the shape functions is used to evaluate the values at the retarded time. In this study, we used linear shape functions. The final form of discretized boundary integral equation is
" # Nj X 2 @/ ~ y j ; ðl þ mÞ Dt @/ ~ y j ; ðl þ mÞ Dt X sm DSj A0ij y j ; ðl þ mÞ Dt þ C0ij þ B0ij / ~ ; 2p/ð~ x i ; kDtÞ ¼ @s @n j¼1 m¼1
(29)
where ~ x i is collocation point, k is designated time step, sm is the mth shape function, and l is discretized retarded time index, 2962
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
kDt Rij =a0 Rij =a0 sr 1 ¼ int k 1: l ¼ int 1 ¼ int Dt Dt Dt
(30)
Rij is the phase radius in Eq. (6) between ith and jth elements and defined as
M0 ðxi1 yj1 Þ þ Rij ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n o 2 ðxi1 yj1 Þ2 þ b ðxi2 yj2 Þ2 þ ðxi3 yj3 Þ2 b2
As stated before, the retarded time has to be evaluated using the phase radius, Rij , not the geometrical radius, in a moving medium. Thus the size of boundary elements has to be reduced as the medium Mach number increases. To obtain the velocity potential at the surface at the designated time step using Eq. (29), the time marching procedure is required. If a0 Dt is taken shorter than the shortest Rij in Eq. (30), l þ m in Eq. (29) is less than k (m ¼ 2 for the linear interpolation shape function) for all elements. This means the velocity potential at the designated observer time step can be obtained by summing up the past information in the right-hand side of Eq. (29). This is the explicit time marching scheme. The explicit time marching scheme does not requires matrix solving process, but many time steps are required because of the small time step Dt mentioned in the preceding text. Further it is reported that the explicit time marching scheme causes the instability.14,15 The error included in numerical calculation of the time derivative term causes this instability. This problem can be solved by using the implicit time marching scheme.19 Further, the implicit time marching scheme allows to select large time step size. For the implicit time marching scheme,
:
(31)
the time step has to be chosen to satisfy Courant–Friedrichs–Lewy (CFL) condition as CFL ¼ a0
Dt > 1; minðRij Þ
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
(32)
where Dt is the time step, minðRij Þ is the shortest phase radius between the elements in Eq. (31). If CFL > 1, the discretized retarded time index l in Eq. (30) is l ¼ k 2 and l þ m in the right-hand side of Eq. (29) becomes k for the linear shape functions (m ¼ 2). Then the right-hand side contains the information at time step k. Thus Eq. (29) has to be solved implicitly. The time marching scheme can be explicit if Rij becomes larger to make the CFL number less than 1 although the size of time step is same. The temporal derivative term in Eq. (29) is evaluated using the finite difference scheme. The central difference scheme is used to increase the order of the scheme. But if CFL > 1, the time marching scheme becomes implicit, and the central difference scheme requires the information at time step k þ 1. This violates the causality. Thus the backward difference is used for the temporal derivative at time step k when CFL > 1. Now, Eq. (29) can be arranged to a linear system ½a f/k g ¼ ½c:
FIG. 2. A sphere used for the test.
the
(33)
FIG. 3. The schematic of the test. Y. W. Lee and D. J. Lee: Convective boundary integral formulation
2963
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
By solving the linear system at each time step, the velocity potential at the designated time step k can be obtained. The coefficient matrix ½a is a sparse matrix, and the sparsity depends on the CFL number. As the CFL number increases, the right-hand side of Eq. (29) contains more information at the designated time step, and the coefficient matrix ½a becomes dense. And the sparsity is dependent on the medium Mach number too because the CFL number increases as the minimum phase radius Rij decreases. Thus selection of an appropriate time step with consideration of the Mach number is needed.
Cð~ x Þ /ð~ xÞ e
ixt
VI. VERIFICATION A. The frequency domain formula
Equation (12) is derived in the time domain. If the variable is harmonic, then the formula can be transformed to the frequency domain. Wu et al. derived a formula in the frequency domain. Thus we verified our result by comparison with the result of Wu et al. The velocity potential is assumed to be harmonic with angular frequency x. Consequently, Eq. (12) can be written as follows,
ð
@ ixs @ ixs @ ixs þ B/eixs þ C þD n1 /e /e ¼ A /e @t @n @y1 S ð @/ @/ ¼ þ Deixs ixA/eixs þ B/eixs þ Ceixs n1 dS: @n @y1 S
dS ret
(34)
s is not a function of space, so it can be placed outside of the normal derivative. By using the relation between the observer time and the retarded time, Eq. (12) can be rearranged,
Cð~ x Þ /ð~ xÞ ¼
ð ixA/e
ikR
þ B/e
ikR
S
þ Ce
ikR
@/ ikR @/ þ De n1 dS: @n @y1
For direct comparison, it is necessary to express Eq. (35) with Green’s function. Green’s function for the convective Helmholtz equation is G¼
eikR : 4pR
(36)
Because the integration is for the surface, the derivative of G has to be for ~ y . The spatial derivative of G is eikR ikeikR ry R ry G ¼ 2 ry R 4pR 4pR eikR ikeikR ¼ rx R: 2 rx R þ 4pR 4pR Thus Eq. (35) can be arranged as follows:
(37)
Cð~ x Þ /ð~ xÞ ¼
(35)
ð @/ @G / 2ikM0 n1 G/ G @n @n S @/ @G / n1 Þ dS: þ M02 G @y1 @y1
(38)
Equation (38) is the Kirchhoff’s formula in the frequency domain and it is the same as the formula of Wu et al.2 B. A point monopole in a moving medium
A test using Eq. (25) is conducted. The test case is a point monopole in a moving medium. The monopole has an amplitude of 1 m3/s and a frequency of 50 Hz, and it is placed at the origin. The medium Mach number is 0.7 in the positive x1 direction, and the corresponding b is 0.714. Because there is no real surface, an artificial surface is generated. The Kirchhoff’s formula has no restriction on the type
FIG. 4. Velocity potential by a point monopole at (a) ð2; 0; 0Þ (b) ð0; 2; 0Þ. 2964
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
FIG. 5. Velocity potential by a point monopole at surface point ð1; 0; 0Þ.
of surface, and the artificial permeable surface can be used as a surface of integration. The surface is a sphere, and the radius of the sphere is 1 m. The sphere is modeled by 1298 triangular elements. For easy numerical implementation, the property at each element is assumed to be constant. Figure 2 shows the sphere used for the test. Figure 3 shows the schematic of the problem. The acoustic property at the surface is obtained from the exact solution for the validation of Eq. (25). The acoustic velocity potential of a point monopole is /ð~ x ; tÞ ¼
1 expðixðt R=a0 ÞÞ: 4pR
(39)
Figure 4 shows good agreement between the result obtained by Eq. (25) and the exact solution of Eq. (39). The velocity potential at ð0; 2; 0Þ shows higher amplitude than ð2; 0; 0Þ. This is the moving medium effect. The amplitude radius R in Eq. (7) changes for the observer position although the geometrical radius is identical. Because the amplitude radius for ð2; 0; 0Þ is b times that for ð0; 2; 0Þ, the amplitude of velocity potential at ð0; 2; 0Þ is 1=b times higher than that at ð2; 0; 0Þ. Also the phase of two points is different. This comes from the difference of the phase radius. Due to the different phase radius, propagation time is different for each observer position.
C. Time domain BEM analysis
Kirchhoff’s formula for a moving medium can also be utilized as a boundary integral equation. A time domain BEM analysis is conducted using the boundary integral equation to see the utility of the equation. The test case is the same as the previous case except that the velocity potential at the surface is unknown. It is obtained by the current time domain BEM analysis derived in this paper. For a simple test, the linear temporal shape function is employed. The time step is 6:25 104 s, and the corresponding CFL number is 3. The velocity potential is obtained by solving the linear system. The GMRes technique is used to accelerate the linear system solving.24 Figure 5 shows the result at surface point ð1; 0; 0Þ. Figure 5 shows good agreement between the surface velocity potential obtained from the time domain BEM and the analytic solution. The instability is not observed and the solution does not diverge. D. Acoustic scattering in a moving medium by a wing-like body
The time domain BEM technique is applied to a scattering problem. Because there is no exact solution for the scattering in a moving medium, we compared our results with previous numerical results by Myers et al.,11 that are
FIG. 6. Geometry of a wing-like body.
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
2965
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
FIG. 7. Acoustic pressure at ð0; 0; 105LÞ: (a) M0 ¼ 0; (b) M0 ¼ 0:5. Dashed line: Incident pressure; dashed-dotted line: Scattered pressure; solid line: Total pressure.
obtained using the frequency domain BEM. A wing-like body is in a moving medium, and a point monopole source is placed above the body. The body has a rectangular wing shape, and the chord length of the wing is 2L, and the span length is 6L. The maximum thickness ratio to the chord length of the wing is 10%, and the location of the maximum thickness is the mid-chord. The center of the body is set to the origin of the coordinate. Figure 6 shows the geometry of the body. A point monopole source is located above the body and its position is ð0; 0; 10LÞ. The velocity potential of a point monopole in a moving medium is written in Eq. (39). The non-dimensionalized wave number kL is set to 1. The medium is moving in the positive x1 direction with Mach number of 0 and 0.5. To compare the present results with those of Myers et al., we obtained the acoustic pressure using the velocity potential as @/ @/ 0 : (40) þ U0 p ¼ q0 @t @x1 The body surface is discretized into 1450 triangular elements for M0 ¼ 0 and 5730 triangular elements for M0 ¼ 0:5. Due to the moving medium effect, the wave length is
shortened in the upstream direction, and more elements are used. The size of time step is 3:36 104 s, and the corresponding CFL number is 3.57 for M0 ¼ 0 and 10.43 for M0 ¼ 0:5. Thus the implicit time marching scheme is used. The GMRes technique is used to accelerate the linear system solving also. Figure 7 shows the signal of acoustic pressure at ð0; 0; 105LÞ for M0 ¼ 0 and M0 ¼ 0:5. The phase difference between the incident pressure and the scattered pressure is observed, and the total pressure contains the phase interference. The phase interference is more clearly observed when M0 ¼ 0:5. Because of the phase interference, amplitude of the total pressure is lower than the incident pressure. Figure 8 shows directivity patterns of the total pressure. The directivity is plotted in the x1 x3 plane, which is the midspan plane and the radius is 105L. h ¼ 0 is the positive x1 direction. To guarantee the validity of the result, the result obtained using frequency domain BEM analysis is plotted. Myers et al. introduced a pressure directivity function D and the directivity is plotted with the function, D¼
rob p0rms ; L p00
(41)
FIG. 8. Acoustic pressure at ð0; 0; 105LÞ: (a) M0 ¼ 0; (b) M0 ¼ 0:5. Diamond line: Myers and Hausmann (Ref. 11); dashed-dotted line: Frequency domain BEM; solid line: Present(time domain BEM). 2966
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
where p00 is the RMS incident pressure at the origin The directivity pattern obtained by present method shows similar pattern to the results of Myers et al. Spikes are observed in the directivity pattern as a result of phase interference and the moving medium effect broke the symmetry. The amplitude of the directivity function is higher in the upstream direction than the downstream direction when the medium is moving and this is because of the moving medium effect. VII. CONCLUSION
Kirchhoff’s formula for the convective acoustic wave equation is derived in the physical coordinates using generalized function theory. An additional term appears due to moving medium, and it is manipulated by decomposing it to the surface normal component and tangential component for easy numerical implementation. The final formula is rearranged in the form of the classical Kirchhoff’s formula. If the variable is harmonic, the derived formula is changed to the frequency domain formula. The transformed formula is exactly the same as the previously derived frequency domain formula. The formula is also verified by comparing with a simple exact solution. A point monopole in a moving medium case is tested, and the result shows a good agreement with the exact solution. If we set the observer position on the surface, the formula changes to a boundary integral equation. The velocity potential on the left-hand side of the formula becomes half when the observer is placed on the surface, and it is the same as the classical Kirchhoff’s boundary integral equation. We set a simple numerical technique to solve the boundary integral equation in time domain. The same technique as the classical time domain BEM technique is introduced, but the medium Mach number is considered. The phase radius related to the retarded time is dependent on the medium Mach number, and it influences to the size of the time step and the element. A time domain BEM analysis is conducted to verify the convective Kirchhoff’s integral equation. An implicit time marching scheme is used, and the result is obtained by solving a linear system. The CFL number for a moving medium is higher than a stationary medium because it requires more propagation time in the upstream direction. This means that a relatively shorter time step is required in the moving medium analysis for the same CFL number. The result obtained by the time domain BEM analysis for a point monopole in a moving medium showed a good agreement with the exact solution. The acoustic scattering by a wing-like body was solved using the time domain BEM analysis, and the result showed good agreement also. Thus the derived time domain boundary integral formula can be used for sound radiation and scattering in a moving medium.
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
1
D. Blokhintzev, “The propagation of sound in an inhomogeneous and moving medium I,” J. Acoust. Soc. Am. 18(2), 322–328 (1946). 2 T. Wu and L. Lee, “A direct boundary integral formulation for acoustic radiation in a subsonic uniform flow,” J. Sound Vib. 175(1), 51–63 (1994). 3 Y. Guo, “Computation of sound propagation by boundary element method,” NASA Contract Report No. NAS1-00086-A003 (2005) available at http://ntrs.nasa.gov (Last viewed 15 October 2014). 4 M. J. Lighthill, An Introduction to Fourier Analysis and Generalised Functions (Cambridge University Press, Cambridge, UK, 1958), pp. 10–27. 5 J. F. Williams and D. L. Hawkings, “Sound generation by turbulence and surfaces in arbitrary motion,” Philos. Trans. R. Soc. Ser. A. 264(1151), 321–342 (1969). 6 M. J. Lighthill, “On sound generated aerodynamically. I. General theory,” Proc. R. Soc. London Ser. A. 211(1107), 564–587 (1952). 7 A. Najafi-Yazdi, G. A. Bres, and L. Mongeau, “An acoustic analogy formulation for moving sources in uniformly moving media,” Proc. R. Soc. A Math. Phys. 467(2125), 144–165 (2011). 8 F. Farassat, “Discontinuities in aerodynamics and aeroacoustics: The concept and applications of generalized derivatives,” J. Sound Vib. 55(2), 165–193 (1977). 9 F. Farassat and M. Myers, “Extension of Kirchhoff’s formula to radiation from moving surfaces,” J. Sound Vib. 123(3), 451–460 (1988). 10 M. Myers and J. Hausmann, “On the application of the Kirchhoff formula for moving surfaces,” J. Sound Vib. 139(1), 174–178 (1990). 11 M. Myers and J. Hausmann, “Computation of acoustic scattering from a moving rigid surface,” J. Acoust. Soc. Am. 91(5), 2594–2605 (1992). 12 B. Rynne, “Stability and convergence of time marching methods in scattering problems,” IMA J. Appl. Math. 35(3), 297–310 (1985). 13 B. Rynne, “Instabilities in time marching methods for scattering problems,” Electromagnetics 6(2), 129–144 (1986). 14 P. Smith, “Instabilities in time marching methods for scattering: Cause and rectification,” Electromagnetics 10(4), 439–451 (1990). 15 B. Rynne and P. Smith, “Stability of time marching algorithms for the electric field integral equation,” J. Electromagnet. Wave. 4(12), 1181–1205 (1990). 16 P. J. Davies and D. B. Duncan, “Averaging techniques for time–marching schemes for retarded potential integral equations,” Appl. Num. Math. 23(3), 291–310 (1997). 17 D. Chappell, P. Harris, D. Henwood, and R. Chakrabarti, “A stable boundary element method for modeling transient acoustic radiation,” J. Acoust. Soc. Am. 120(1), 74–80 (2006). 18 G. Manara, A. Monorchio, and R. Reggiannini, “A space-time discretization criterion for a stable time-marching solution of the electric field integral equation,” IEEE T. Antenna Propag. 45(3), 527–532 (1997). 19 M. Bluck and S. Walker, “Analysis of three-dimensional transient acoustic wave propagation using the boundary integral equation method,” Int. J. Num. Methods Eng. 39(8), 1419–1431 (1996). 20 H. Wang, D. Henwood, P. Harris, and R. Chakrabarti, “Concerning the cause of instability in time-stepping boundary element methods applied to the exterior acoustic problem,” J. Sound Vib. 305(1), 289–297 (2007). 21 A. Ergin, B. Shanker, and E. Michielssen, “Analysis of transient wave scattering from rigid bodies using a Burton-Miller approach,” J. Acoust. Soc. Am. 106(5), 2396–2404 (1999). 22 H.-W. Jang and J.-G. Ih, “Stabilization of time domain acoustic boundary element method for the interior problem with impedance boundary conditions,” J. Acoust. Soc. Am. 131(4), 2742–2752 (2012). 23 M. Bluck and S. Walker, “Time-domain BIE analysis of large threedimensional electromagnetic scattering problems,” IEEE T. Antenna Propag. 45(5), 894–901 (1997). 24 V. Fraysse, L. Giraud, S. Gratton, and J. Langou, “Algorithm 842: A set of GMRes routines for real and complex arithmetics on high performance computers,” ACM T. Math. Software 31(2), 228–238 (2005).
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
2967
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
(Received 3 June 2014; revised 21 September 2014; accepted 6 October 2014) Kirchhoff’s formula for the convective wave equation is derived using the generalized function theory. The generalized convective wave equation for a stationary surface is obtained, and the integral formulation, the convective Kirchhoff’s formula, is derived. The formula has a similar form to the classical Kirchhoff’s formula, but an additional term appears due to a moving medium effect. For convenience, the additional term is manipulated to a final form as the classical Kirchhoff’s formula. The frequency domain boundary integral can be obtained from the current time domain boundary integral form. The derived formula is verified by comparison with the analytic solution of source in the uniform flow. The formula is also utilized as a boundary integral equation. Time domain boundary element method (BEM) analysis using the boundary integral equation is conducted, and the results show good agreement with the analytical solution. The formula derived here can be useful for sound radiation and scattering by arbitrary C 2014 Acoustical Society of America. bodies in a moving medium in the time domain. V [http://dx.doi.org/10.1121/1.4898427] PACS number(s): 43.20.Rz, 43.20.Bi [JDM] I. INTRODUCTION
The boundary integral formula is useful for sound radiation problems. The well-known Helmholtz integral formula is a boundary integral formula for the Helmholtz equation in frequency domain. The acoustic property at the observer position can be obtained through boundary integration. Another boundary integral formula is the Kirchhoff’s formula, which is the formula for wave equation in the time domain. The two formulas in the preceding text are based on the wave equation with stationary medium. For moving medium, Blokhintzev1 derived an integral formula for the convective wave equation using coordinate transformation. Blokhintzev transformed convective wave equation using a new system of coordinates and derived the final form using Green’s theorem in time domain. Blokhintzev’s integral formula contains surface integration in transformed coordinates. Wu et al.2 derived the convective integral formula in the frequency domain using the weighted residual integration and Green’s theorem. The integration of the formula of Wu et al. is in physical coordinates, and it is convenient to use. Guo3 also obtained an integral formula in the frequency domain using weighted residual integration and Green’s theorem. Another method for the derivation of the boundary integral formula is the generalized function theory.4 Ffowcs-Williams and Hawking5 extended Lighthills6 acoustic analogy theory for a moving body using this generalized function theory. The generalized function theory makes it possible to handle discontinuous areas, such as body surface, and it simplifies the derivation process. Recently, the convective acoustic analogy formula has been derived7 using the generalized function theory. Farassat et al. utilized the generalized function theory for the derivation a)
Author to whom correspondence should be addressed. Electronic mail: [email protected]
J. Acoust. Soc. Am. 136 (6), December 2014
Pages: 2959–2967 of Kirchhoff’s formula.8,9 Farassat et al. set a moving surface surrounding the moving acoustic sources and used the generalized potential, which is discontinuous at the surface. Kirchhoff’s formula offers an advantage in which it can be utilized as a boundary integral equation itself to solve the radiation and scattering problem for the moving surface.10,11 Thus the scattered field can be obtained by solving the boundary integral equation in the time domain. The time domain boundary element method (BEM) analysis has stability problems.12,13 An averaging scheme is used to suppress the instability, but it increases the computational cost.14–16 The high-order approximation for space and time is also used to reduce the instability problem.17–19 For the exterior problem, the fictitious internal resonance in the body also can cause the instability problem.20 To solve this internal resonance problem, the time domain Burton–Miller formulation is introduced.21 But the Burton–Miller formulation requires an additional equation for the normal derivative of the variable and it increases computational cost. Recently, the vector filtering technique was introduced to suppress the instability by filtering the instability wave, and it shows good results.22 In this paper, we derived the Kirchhoff’s formula for the convective wave equation in physical coordinates. The obtained formula contains an additional term compared to the classical Kirchhoff’s formula due to a moving medium effect. The additional term was modified to have the same form as in the classical Kirchhoff’s formula for easy implementation. The generalized function theory was employed for the derivation of the formula. The derived formula is verified by comparison with an analytic solution of simple source. In addition, time domain BEM analysis is conducted to obtain the utility of the formula as a boundary integral equation, which can handle the radiation and scattering problems for arbitrary bodies in a moving medium. We used a simple implicit time marching scheme to avoid the longtime instability problem.23 If the frequency of the acoustic
0001-4966/2014/136(6)/2959/9/$30.00
C 2014 Acoustical Society of America V
2959
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
wave does not match with the internal resonance frequency, the instability due to the internal resonance does not occur. Thus the implicit time marching scheme is appropriate to see the utility of the formula. The numerical technique for the convective time domain BEM using implicit time marching scheme is demonstrated, and a simple test case is solved using the time domain BEM technique without instability problem.
ðfÞ @H @f @f ¼ H0 ð f Þ ¼ dð f Þ ; @x @x @x
II. DERIVATION OF GENERALIZED CONVECTIVE ACOUSTIC WAVE EQUATION
To obtain Kirchhoff’s formula for the convective acoustic wave equation, the generalized convective wave equation was derived first. For simplification, we assumed the medium moves in the positive x direction with constant speed U0 . Then the convective wave equation can be written as 1 @2/ U0 @ 2 / U2 @ 2 / þ2 2 þ 20 2 r2x / ¼ 0; 2 2 a0 @t a0 @t@x1 a0 @x1
(1)
where / is the velocity potential, a0 is the speed of sound, x1 is the component of the position vector ~ x ¼ ðx1 ; x2 ; x3 Þ, and x . If rx denotes gradient for the observer position vector ~ there are bodies, Eq. (1) is valid in the region outside the body surface. The body surface can be any kind of surface, including a rigid and permeable surface. The surface can be represented by f ð~ x ; tÞ ¼ 0 such that f ð~ x ; tÞ < 0 for the region inside the surface and f ð~ x ; tÞ > 0 for the region outside the surface. Function f has set as jrx f j ¼ 1, which simplifies the derivation process.9 Consequently, rx f is the surface normal unit vector n^ pointing out from the surface. We assumed that Kirchhoff’s surface is stationary so f is a function of ~ x only (Fig. 1). Then / can be extended inside the surface by assuming that / ¼ 0 inside the surface using the function f. The generalized potential is written using the Heaviside step function as9 ~ ¼ Hðf Þ / ¼ /
/ for f > 0 0 for f < 0:
Notice that the generalized variable is discontinuous at f ¼ 0, and it is not differentiable. However, the generalized function theory gives derivatives to this discontinuous variable, the so-called the generalized derivative.8,9 The Heaviside step function is the generalized function and its derivative is the Dirac-delta function. Thus the generalized derivative of the Heaviside step function can be written as
(2)
(3)
where the overbar denotes the generalized differentiation to distinguish from the ordinary differentiation. By utilizing the above generalized differentiation, the convective wave equation can be changed to the generalized convective wave equation that governs the whole region. By using the generalized derivatives, the generalized convective wave equation can be rewritten as follows: 2~ 2~ 2~ 1 @ / U0 @ / U02 @ / ~ 2/ þ 2 þ r x a20 @t2 a20 @t@x1 a20 @x21 U0 @/ U 2 @/ n1 dð f Þ þ 20 n1 dð f Þ ¼2 2 a0 @t a0 @x1 @ þ ½/n1 dð f Þ @x1
@/ x ½/^ dð f Þ r n dð f Þ: @n
(4)
Equation (4) is the convective acoustic wave equation governing the whole region. Equation (4) has source terms on the right-hand side and they are the surface sources. III. DERIVATION OF KIRCHHOFF’S FORMULA FOR THE CONVECTIVE WAVE EQUATION
The generalized convective wave equation has source terms on the right-hand side, and the solution can be obtained by using the Greens function. Greens function for the convective wave equation is G¼
dðs t þ R=a0 Þ ; 4pR
(5)
where R¼
M0 ðx1 y1 Þ þ R
;
(6)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ ðx1 y1 Þ2 þ b2 ½ðx2 y2 Þ2 þ ðx3 y3 Þ2 ;
(7)
b2
b¼
FIG. 1. The schematic of a stationary surface in a moving medium. 2960
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 M02 :
(8)
Here s is source time, t is observer time, x is observer position, y is source position, and M0 ¼ U0 =a0 . The difference of the Greens function for the convective equation is that the radius is changed to the amplitude radius R and the phase radius R. The solution of Eq. (4) can be obtained using the Greens function in the preceding text as Y. W. Lee and D. J. Lee: Convective boundary integral formulation
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
ð ð ðt ð dð gÞ dð gÞ dðgÞ @/ U02 t @/ n1 dð f Þ dVds þ 2 n1 dð f Þ dVds rx / n^dð f Þ dVds 3 3 3 4pR 4pR 4pR @s @y a 1 1 R 1 R 0 1 R ðt ð ð ð dðgÞ dðgÞ U 2 @ t x /n1 dð f Þ dVds r /^ n dð f Þ dVds; þ 20 4pR 4pR 3 a0 @x1 1 R3 1 R
2U0 ~ ð~ / x ; tÞ ¼ 2 a0
ðt
ð
(9)
where g ¼ s t þ R=a0 . All integrals of Eq. (9) have Dirac-delta functions, and these restrict the integration range to zero of the function variable. Thus the integration for space reduces to the surface, f ¼ 0, and the integration for time reduces to retarded time which is the solution of g ¼ 0, sr ¼ t
R : a0
(10)
This means that the integrand has to be evaluated at retarded time sr for desired observer time. It should be noted that the retarded time is defined with the phase radius R in Eq. (6), which is not the geometric radius. When the medium moves, the propagation time R=a0 of the acoustic wave is different from that of a stationary medium. The last two integrals have spatial derivatives outside the integral and the spatial derivative can be changed to a temporal derivative using the property of the Dirac-delta function.7 Hence, Eq. (9) can be changed to ð ð ð @/ 1 U02 @/ 1 @/ 1 U02 @ R~1 n1 dS þ 2 n1 dS dS 3 /n1 dS 4pR ret 4pR ret a0 f ¼0 @y1 4pR ret a0 @t f ¼0 f ¼0 @s f ¼0 @n 4pR ret ð ð ð U02 1 @ rx R rx R R~1 /n1 dSþ /^ n dSþ /^ n dVds; 2 a0 @t f ¼0 4pR ret a0 f ¼0 4pR2 ret 4pR2 ret f ¼0
U0 ~ ð~ / x ; tÞ ¼ 2 2 a0
ð
where R~1 ¼ @R=@x1 , and R~1 ¼ @R =@x1 . The bracket ½ ret means that the integrands have to be evaluated at retarded time sr . If Eq. (11) is used directly, we need to differentiate the integral at the observer position with respect to the time due to the temporal derivative outside the integral. It is convenient to use if we move the temporal derivative from outside the integral to inside the integral by utilizing the relation between observer time and retarded time sr t þ R=a0 ¼ 0. Because R and a0 are constant, the observer time derivative can be replaced by the source time derivative directly. Then Eq. (11) can be arranged as ð @/ @/ @/ ~ þ B/ þ C þD A n1 dS; 4p/ ð~ x ; tÞ ¼ @s @n @y1 S ret (12)
similar form to the classical Kirchhoff’s formula except for the coefficients and @/=@y1 term. Equation (12) requires additional information at surface @/=@y1 , which does not appear in the classical Kirchhoff’s formula. If we can obtain @/=@y1 at the surface, Eq. (12) can be employed directly, but it usually requires additional computation. Wu et al.2 and Guo3 manipulated this term in the frequency domain through decomposition of gradient to the surface normal and tangential component; this resulted in @/=@y1 changing to a combination of @/=@n and /. Using the preceding approach in this time domain, the second form of Kirchhof’s formula for the convective wave equation was obtained, ð @/ @/ ~ ð~ þ B0 / þ C0 A0 dS; (17) 4p/ x ; tÞ ¼ @s @n ret S where
where 2M0 n1 M02 n1 R~1 þ n^ rx R ; A¼ a0 R B¼
(11)
M02 n1 R~1 þ n^ rx R ; R2
1 C¼ ; R M2 D ¼ 0 : R
(13) (14)
(15) (16)
Equation (12) is Kirchhoff’s formula for the convective acoustic wave equation in physical coordinates, and it has a J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
2M0 n1 M02 n1 R~1 þ n^ rx R ; a0 R M02 n1 R~1 þ n^ rx R Mn ~ 0 Ms ; rs B ¼ R R2 A0 ¼
C0 ¼
1 Mn2 ; R
(18) (19) (20)
~0 n^, M ~s ¼ M ~0 Mn n^, and rs ¼ r @=@n. where Mn ¼ M Equation (17) now has the same form as the classical Kirchhoff’s formula. Thus Eq. (17) can be utilized in the same way as the classical Kirchhoff’s formula with the coefficients in Eqs. (18) to (20). Y. W. Lee and D. J. Lee: Convective boundary integral formulation
2961
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
IV. BOUNDARY INTEGRAL EQUATION FORM OF THE CONVECTIVE KIRCHHOFF’S FORMULA
If the observer is on the surface, the formula changes to a boundary integral equation. Equation (17) is valid for the observer position in the region outside the surface. To obtain the boundary integral equation, we consider that the observer approaches the surface from the region outside the surface. The domain is augmented by a hemi-ellipsoid, and its semiprincipal axes are in the positive x1 ; x2 ; x3 direction. The length of semi-principal axes is =b in the positive x2 and x3 direction and in the positive x1 direction. If tends to zero, the observer approaches the surface, and the boundary integral equation will be obtained. Kirchhoff’s formula for the augmented domain is ð 0 @/ 0 0 @/ þB/þC A dS 4p/ð~ x ; tÞ ¼ @s @n ret S ð @/ @/ þ B0 / þ C0 þ A0 dS: (21) @s @n ret S We consider the region outside the surface and tilde over the velocity potential is removed. Because tends to zero, it is clear that all terms in the second integral converges to zero except the / term. The second integral of Eq. (21) can be written as ð M02 n1 R~1 þ n^ rx R / dS: (22) R2 ret S To evaluate Eq. (22) more conveniently, we apply the Prandtl–Glauert transformation x1 ¼ x1 ; x2 ¼ bx2 ; x3 ¼ bx3 :
(23)
By applying the Prandtl–Glauert transformation, the hemiellipsoid changes to a hemi-sphere with radius . Equation (22) then changes to ð n rx R / dS; (24) 2 S R ret where n ¼ ry f is the surface normal vector in transformed x 2 y2 Þ2 þ ð x 3 y3 Þ2 g0:5 domain, and R ¼ fð x 1 y1 Þ2 þ ð is the distance between the observer and source in transformed domain. If tends to zero, Eq. (24) will be 2p/ð~ x ; tÞ. Thus the boundary integral formula can be arranged as ð 0 @/ 0 0 @/ ð Þ þB/þC A dS; (25) C~ x /ð~ x ; tÞ ¼ @s @n ret S where Cð~ x Þ is
8 < 4p region outside the surface Cð~ x Þ ¼ 2p on the smooth surface : 0 region inside the surface:
(26)
V. IMPLEMENTATION OF THE KIRCHHOFF CONVECTIVE BOUNDARY INTEGRAL EQUATION
The formula has been changed into the form of the classical Kirchhoff’s formula, but the surface gradient term of Eq. (19) is still not easy to handle. Guo3 solved this problem for constant element in the frequency domain BEM. Because each element of Kirchhoff’s surface is constant, / can be placed outside the integral, and the integrand will be a surface gradient term only. This is utilized in the time domain problem. This can be changed to a line integral by divergence theorem, ð Mn ~ Mn ~ rs M s dS ¼ M s n~dl; R S l R
ð
(27)
where n~ is a line normal unit vector pointing outward from the element and parallel to the element. After discretizing the surface into N elements, Eq. (25) can be written in discretized form as Cð~ x Þ /ð~ x ; tÞ ¼
N X j¼1
@/j @/j þ Bj 0 /j þ Cj 0 Aj @s @n 0
DSj : ret
(28) If we know the information at the surface, Eq. (28) can be used easily to obtain the velocity potential at the observer position. However, if we do not know the information at the surface, it can be obtained using Eq. (28) by positioning the observer on the surface with the time domain BEM technique. Generally, the surface normal derivative term is given as a boundary condition, and the unknown is the velocity potential. To obtain the velocity potential at the surface using Eq. (28), the time has to be discretized too. The velocity potential at a certain time step can be obtained by summing up the velocity potential and its spatial and temporal derivative terms at the corresponding retarded time. But generally the retarded time does not coincide with the discretized time steps although the observer time step is discretized because the propagation time R=a0 is not an integer order. Thus the values at the retarded time, located between two integer time steps, have to be evaluated using interpolation with the corresponding integer time steps. An interpolation method using the shape functions is used to evaluate the values at the retarded time. In this study, we used linear shape functions. The final form of discretized boundary integral equation is
" # Nj X 2 @/ ~ y j ; ðl þ mÞ Dt @/ ~ y j ; ðl þ mÞ Dt X sm DSj A0ij y j ; ðl þ mÞ Dt þ C0ij þ B0ij / ~ ; 2p/ð~ x i ; kDtÞ ¼ @s @n j¼1 m¼1
(29)
where ~ x i is collocation point, k is designated time step, sm is the mth shape function, and l is discretized retarded time index, 2962
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
kDt Rij =a0 Rij =a0 sr 1 ¼ int k 1: l ¼ int 1 ¼ int Dt Dt Dt
(30)
Rij is the phase radius in Eq. (6) between ith and jth elements and defined as
M0 ðxi1 yj1 Þ þ Rij ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n o 2 ðxi1 yj1 Þ2 þ b ðxi2 yj2 Þ2 þ ðxi3 yj3 Þ2 b2
As stated before, the retarded time has to be evaluated using the phase radius, Rij , not the geometrical radius, in a moving medium. Thus the size of boundary elements has to be reduced as the medium Mach number increases. To obtain the velocity potential at the surface at the designated time step using Eq. (29), the time marching procedure is required. If a0 Dt is taken shorter than the shortest Rij in Eq. (30), l þ m in Eq. (29) is less than k (m ¼ 2 for the linear interpolation shape function) for all elements. This means the velocity potential at the designated observer time step can be obtained by summing up the past information in the right-hand side of Eq. (29). This is the explicit time marching scheme. The explicit time marching scheme does not requires matrix solving process, but many time steps are required because of the small time step Dt mentioned in the preceding text. Further it is reported that the explicit time marching scheme causes the instability.14,15 The error included in numerical calculation of the time derivative term causes this instability. This problem can be solved by using the implicit time marching scheme.19 Further, the implicit time marching scheme allows to select large time step size. For the implicit time marching scheme,
:
(31)
the time step has to be chosen to satisfy Courant–Friedrichs–Lewy (CFL) condition as CFL ¼ a0
Dt > 1; minðRij Þ
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
(32)
where Dt is the time step, minðRij Þ is the shortest phase radius between the elements in Eq. (31). If CFL > 1, the discretized retarded time index l in Eq. (30) is l ¼ k 2 and l þ m in the right-hand side of Eq. (29) becomes k for the linear shape functions (m ¼ 2). Then the right-hand side contains the information at time step k. Thus Eq. (29) has to be solved implicitly. The time marching scheme can be explicit if Rij becomes larger to make the CFL number less than 1 although the size of time step is same. The temporal derivative term in Eq. (29) is evaluated using the finite difference scheme. The central difference scheme is used to increase the order of the scheme. But if CFL > 1, the time marching scheme becomes implicit, and the central difference scheme requires the information at time step k þ 1. This violates the causality. Thus the backward difference is used for the temporal derivative at time step k when CFL > 1. Now, Eq. (29) can be arranged to a linear system ½a f/k g ¼ ½c:
FIG. 2. A sphere used for the test.
the
(33)
FIG. 3. The schematic of the test. Y. W. Lee and D. J. Lee: Convective boundary integral formulation
2963
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
By solving the linear system at each time step, the velocity potential at the designated time step k can be obtained. The coefficient matrix ½a is a sparse matrix, and the sparsity depends on the CFL number. As the CFL number increases, the right-hand side of Eq. (29) contains more information at the designated time step, and the coefficient matrix ½a becomes dense. And the sparsity is dependent on the medium Mach number too because the CFL number increases as the minimum phase radius Rij decreases. Thus selection of an appropriate time step with consideration of the Mach number is needed.
Cð~ x Þ /ð~ xÞ e
ixt
VI. VERIFICATION A. The frequency domain formula
Equation (12) is derived in the time domain. If the variable is harmonic, then the formula can be transformed to the frequency domain. Wu et al. derived a formula in the frequency domain. Thus we verified our result by comparison with the result of Wu et al. The velocity potential is assumed to be harmonic with angular frequency x. Consequently, Eq. (12) can be written as follows,
ð
@ ixs @ ixs @ ixs þ B/eixs þ C þD n1 /e /e ¼ A /e @t @n @y1 S ð @/ @/ ¼ þ Deixs ixA/eixs þ B/eixs þ Ceixs n1 dS: @n @y1 S
dS ret
(34)
s is not a function of space, so it can be placed outside of the normal derivative. By using the relation between the observer time and the retarded time, Eq. (12) can be rearranged,
Cð~ x Þ /ð~ xÞ ¼
ð ixA/e
ikR
þ B/e
ikR
S
þ Ce
ikR
@/ ikR @/ þ De n1 dS: @n @y1
For direct comparison, it is necessary to express Eq. (35) with Green’s function. Green’s function for the convective Helmholtz equation is G¼
eikR : 4pR
(36)
Because the integration is for the surface, the derivative of G has to be for ~ y . The spatial derivative of G is eikR ikeikR ry R ry G ¼ 2 ry R 4pR 4pR eikR ikeikR ¼ rx R: 2 rx R þ 4pR 4pR Thus Eq. (35) can be arranged as follows:
(37)
Cð~ x Þ /ð~ xÞ ¼
(35)
ð @/ @G / 2ikM0 n1 G/ G @n @n S @/ @G / n1 Þ dS: þ M02 G @y1 @y1
(38)
Equation (38) is the Kirchhoff’s formula in the frequency domain and it is the same as the formula of Wu et al.2 B. A point monopole in a moving medium
A test using Eq. (25) is conducted. The test case is a point monopole in a moving medium. The monopole has an amplitude of 1 m3/s and a frequency of 50 Hz, and it is placed at the origin. The medium Mach number is 0.7 in the positive x1 direction, and the corresponding b is 0.714. Because there is no real surface, an artificial surface is generated. The Kirchhoff’s formula has no restriction on the type
FIG. 4. Velocity potential by a point monopole at (a) ð2; 0; 0Þ (b) ð0; 2; 0Þ. 2964
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
FIG. 5. Velocity potential by a point monopole at surface point ð1; 0; 0Þ.
of surface, and the artificial permeable surface can be used as a surface of integration. The surface is a sphere, and the radius of the sphere is 1 m. The sphere is modeled by 1298 triangular elements. For easy numerical implementation, the property at each element is assumed to be constant. Figure 2 shows the sphere used for the test. Figure 3 shows the schematic of the problem. The acoustic property at the surface is obtained from the exact solution for the validation of Eq. (25). The acoustic velocity potential of a point monopole is /ð~ x ; tÞ ¼
1 expðixðt R=a0 ÞÞ: 4pR
(39)
Figure 4 shows good agreement between the result obtained by Eq. (25) and the exact solution of Eq. (39). The velocity potential at ð0; 2; 0Þ shows higher amplitude than ð2; 0; 0Þ. This is the moving medium effect. The amplitude radius R in Eq. (7) changes for the observer position although the geometrical radius is identical. Because the amplitude radius for ð2; 0; 0Þ is b times that for ð0; 2; 0Þ, the amplitude of velocity potential at ð0; 2; 0Þ is 1=b times higher than that at ð2; 0; 0Þ. Also the phase of two points is different. This comes from the difference of the phase radius. Due to the different phase radius, propagation time is different for each observer position.
C. Time domain BEM analysis
Kirchhoff’s formula for a moving medium can also be utilized as a boundary integral equation. A time domain BEM analysis is conducted using the boundary integral equation to see the utility of the equation. The test case is the same as the previous case except that the velocity potential at the surface is unknown. It is obtained by the current time domain BEM analysis derived in this paper. For a simple test, the linear temporal shape function is employed. The time step is 6:25 104 s, and the corresponding CFL number is 3. The velocity potential is obtained by solving the linear system. The GMRes technique is used to accelerate the linear system solving.24 Figure 5 shows the result at surface point ð1; 0; 0Þ. Figure 5 shows good agreement between the surface velocity potential obtained from the time domain BEM and the analytic solution. The instability is not observed and the solution does not diverge. D. Acoustic scattering in a moving medium by a wing-like body
The time domain BEM technique is applied to a scattering problem. Because there is no exact solution for the scattering in a moving medium, we compared our results with previous numerical results by Myers et al.,11 that are
FIG. 6. Geometry of a wing-like body.
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
2965
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
FIG. 7. Acoustic pressure at ð0; 0; 105LÞ: (a) M0 ¼ 0; (b) M0 ¼ 0:5. Dashed line: Incident pressure; dashed-dotted line: Scattered pressure; solid line: Total pressure.
obtained using the frequency domain BEM. A wing-like body is in a moving medium, and a point monopole source is placed above the body. The body has a rectangular wing shape, and the chord length of the wing is 2L, and the span length is 6L. The maximum thickness ratio to the chord length of the wing is 10%, and the location of the maximum thickness is the mid-chord. The center of the body is set to the origin of the coordinate. Figure 6 shows the geometry of the body. A point monopole source is located above the body and its position is ð0; 0; 10LÞ. The velocity potential of a point monopole in a moving medium is written in Eq. (39). The non-dimensionalized wave number kL is set to 1. The medium is moving in the positive x1 direction with Mach number of 0 and 0.5. To compare the present results with those of Myers et al., we obtained the acoustic pressure using the velocity potential as @/ @/ 0 : (40) þ U0 p ¼ q0 @t @x1 The body surface is discretized into 1450 triangular elements for M0 ¼ 0 and 5730 triangular elements for M0 ¼ 0:5. Due to the moving medium effect, the wave length is
shortened in the upstream direction, and more elements are used. The size of time step is 3:36 104 s, and the corresponding CFL number is 3.57 for M0 ¼ 0 and 10.43 for M0 ¼ 0:5. Thus the implicit time marching scheme is used. The GMRes technique is used to accelerate the linear system solving also. Figure 7 shows the signal of acoustic pressure at ð0; 0; 105LÞ for M0 ¼ 0 and M0 ¼ 0:5. The phase difference between the incident pressure and the scattered pressure is observed, and the total pressure contains the phase interference. The phase interference is more clearly observed when M0 ¼ 0:5. Because of the phase interference, amplitude of the total pressure is lower than the incident pressure. Figure 8 shows directivity patterns of the total pressure. The directivity is plotted in the x1 x3 plane, which is the midspan plane and the radius is 105L. h ¼ 0 is the positive x1 direction. To guarantee the validity of the result, the result obtained using frequency domain BEM analysis is plotted. Myers et al. introduced a pressure directivity function D and the directivity is plotted with the function, D¼
rob p0rms ; L p00
(41)
FIG. 8. Acoustic pressure at ð0; 0; 105LÞ: (a) M0 ¼ 0; (b) M0 ¼ 0:5. Diamond line: Myers and Hausmann (Ref. 11); dashed-dotted line: Frequency domain BEM; solid line: Present(time domain BEM). 2966
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
where p00 is the RMS incident pressure at the origin The directivity pattern obtained by present method shows similar pattern to the results of Myers et al. Spikes are observed in the directivity pattern as a result of phase interference and the moving medium effect broke the symmetry. The amplitude of the directivity function is higher in the upstream direction than the downstream direction when the medium is moving and this is because of the moving medium effect. VII. CONCLUSION
Kirchhoff’s formula for the convective acoustic wave equation is derived in the physical coordinates using generalized function theory. An additional term appears due to moving medium, and it is manipulated by decomposing it to the surface normal component and tangential component for easy numerical implementation. The final formula is rearranged in the form of the classical Kirchhoff’s formula. If the variable is harmonic, the derived formula is changed to the frequency domain formula. The transformed formula is exactly the same as the previously derived frequency domain formula. The formula is also verified by comparing with a simple exact solution. A point monopole in a moving medium case is tested, and the result shows a good agreement with the exact solution. If we set the observer position on the surface, the formula changes to a boundary integral equation. The velocity potential on the left-hand side of the formula becomes half when the observer is placed on the surface, and it is the same as the classical Kirchhoff’s boundary integral equation. We set a simple numerical technique to solve the boundary integral equation in time domain. The same technique as the classical time domain BEM technique is introduced, but the medium Mach number is considered. The phase radius related to the retarded time is dependent on the medium Mach number, and it influences to the size of the time step and the element. A time domain BEM analysis is conducted to verify the convective Kirchhoff’s integral equation. An implicit time marching scheme is used, and the result is obtained by solving a linear system. The CFL number for a moving medium is higher than a stationary medium because it requires more propagation time in the upstream direction. This means that a relatively shorter time step is required in the moving medium analysis for the same CFL number. The result obtained by the time domain BEM analysis for a point monopole in a moving medium showed a good agreement with the exact solution. The acoustic scattering by a wing-like body was solved using the time domain BEM analysis, and the result showed good agreement also. Thus the derived time domain boundary integral formula can be used for sound radiation and scattering in a moving medium.
J. Acoust. Soc. Am., Vol. 136, No. 6, December 2014
1
D. Blokhintzev, “The propagation of sound in an inhomogeneous and moving medium I,” J. Acoust. Soc. Am. 18(2), 322–328 (1946). 2 T. Wu and L. Lee, “A direct boundary integral formulation for acoustic radiation in a subsonic uniform flow,” J. Sound Vib. 175(1), 51–63 (1994). 3 Y. Guo, “Computation of sound propagation by boundary element method,” NASA Contract Report No. NAS1-00086-A003 (2005) available at http://ntrs.nasa.gov (Last viewed 15 October 2014). 4 M. J. Lighthill, An Introduction to Fourier Analysis and Generalised Functions (Cambridge University Press, Cambridge, UK, 1958), pp. 10–27. 5 J. F. Williams and D. L. Hawkings, “Sound generation by turbulence and surfaces in arbitrary motion,” Philos. Trans. R. Soc. Ser. A. 264(1151), 321–342 (1969). 6 M. J. Lighthill, “On sound generated aerodynamically. I. General theory,” Proc. R. Soc. London Ser. A. 211(1107), 564–587 (1952). 7 A. Najafi-Yazdi, G. A. Bres, and L. Mongeau, “An acoustic analogy formulation for moving sources in uniformly moving media,” Proc. R. Soc. A Math. Phys. 467(2125), 144–165 (2011). 8 F. Farassat, “Discontinuities in aerodynamics and aeroacoustics: The concept and applications of generalized derivatives,” J. Sound Vib. 55(2), 165–193 (1977). 9 F. Farassat and M. Myers, “Extension of Kirchhoff’s formula to radiation from moving surfaces,” J. Sound Vib. 123(3), 451–460 (1988). 10 M. Myers and J. Hausmann, “On the application of the Kirchhoff formula for moving surfaces,” J. Sound Vib. 139(1), 174–178 (1990). 11 M. Myers and J. Hausmann, “Computation of acoustic scattering from a moving rigid surface,” J. Acoust. Soc. Am. 91(5), 2594–2605 (1992). 12 B. Rynne, “Stability and convergence of time marching methods in scattering problems,” IMA J. Appl. Math. 35(3), 297–310 (1985). 13 B. Rynne, “Instabilities in time marching methods for scattering problems,” Electromagnetics 6(2), 129–144 (1986). 14 P. Smith, “Instabilities in time marching methods for scattering: Cause and rectification,” Electromagnetics 10(4), 439–451 (1990). 15 B. Rynne and P. Smith, “Stability of time marching algorithms for the electric field integral equation,” J. Electromagnet. Wave. 4(12), 1181–1205 (1990). 16 P. J. Davies and D. B. Duncan, “Averaging techniques for time–marching schemes for retarded potential integral equations,” Appl. Num. Math. 23(3), 291–310 (1997). 17 D. Chappell, P. Harris, D. Henwood, and R. Chakrabarti, “A stable boundary element method for modeling transient acoustic radiation,” J. Acoust. Soc. Am. 120(1), 74–80 (2006). 18 G. Manara, A. Monorchio, and R. Reggiannini, “A space-time discretization criterion for a stable time-marching solution of the electric field integral equation,” IEEE T. Antenna Propag. 45(3), 527–532 (1997). 19 M. Bluck and S. Walker, “Analysis of three-dimensional transient acoustic wave propagation using the boundary integral equation method,” Int. J. Num. Methods Eng. 39(8), 1419–1431 (1996). 20 H. Wang, D. Henwood, P. Harris, and R. Chakrabarti, “Concerning the cause of instability in time-stepping boundary element methods applied to the exterior acoustic problem,” J. Sound Vib. 305(1), 289–297 (2007). 21 A. Ergin, B. Shanker, and E. Michielssen, “Analysis of transient wave scattering from rigid bodies using a Burton-Miller approach,” J. Acoust. Soc. Am. 106(5), 2396–2404 (1999). 22 H.-W. Jang and J.-G. Ih, “Stabilization of time domain acoustic boundary element method for the interior problem with impedance boundary conditions,” J. Acoust. Soc. Am. 131(4), 2742–2752 (2012). 23 M. Bluck and S. Walker, “Time-domain BIE analysis of large threedimensional electromagnetic scattering problems,” IEEE T. Antenna Propag. 45(5), 894–901 (1997). 24 V. Fraysse, L. Giraud, S. Gratton, and J. Langou, “Algorithm 842: A set of GMRes routines for real and complex arithmetics on high performance computers,” ACM T. Math. Software 31(2), 228–238 (2005).
Y. W. Lee and D. J. Lee: Convective boundary integral formulation
2967
Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 134.225.1.226 On: Fri, 26 Dec 2014 13:48:11
