On the instability of time-domain acoustic boundary element method due to the static mode in interior problems

Journal of Sound and Vibration 332 (2013) 6463–6471

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

On the instability of time-domain acoustic boundary element method due to the static mode in interior problems Hae-Won Jang, Jeong-Guon Ih n Center for Noise and Vibration Control, Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Daejeon 305-701, Republic of Korea

a r t i c l e i n f o

abstract

Article history: Received 30 July 2012 Received in revised form 23 June 2013 Accepted 21 July 2013 Handling Editor: J. Astley Available online 15 August 2013

In the analysis of interior acoustic problems, the time domain boundary element method (TBEM) suffers the monotonically increasing instability when using the direct Kirchhoff integral. This instability is related to the non-oscillatory static acoustic mode describing the constant spatial response in an enclosure. In this work, nonphysical natures of non-oscillatory static mode influencing the instability of TBEM calculation are investigated, and a method for stabilization is studied. In TBEM calculation, the static mode is represented by two nonoscillatory eigenmodes with different eigenvalues. The monotonically increasing instability is caused by the unstable poles of non-oscillatory eigenmodes as well as very small, very low frequency noise of an input signal. Interior problems with impedance boundary condition also exhibit the monotonically increasing instability stemming from its pseudo nonoscillatory static mode due to the lack of dissipation at very low frequencies. Calculation of transient sound fields within rigid and lined boxes provides numerical evidences. It is noted that the stabilization effort by modifying the coefficient matrix based on the spectral decomposition can be used only for correcting the unstable pole. The filtering method based on the eigen-analysis must be additionally used to avoid the remaining instability caused by very low frequency noise of input signal. & 2013 Elsevier Ltd. All rights reserved.

1. Introduction The time-domain boundary element method (TBEM) in acoustics employs the time stepping algorithm to solve the Kirchhoff integral [1,2]. Notwithstanding that this method has a great potential to solve various transient wave propagation problems, it has not been used widely due to the exponentially diverging instability. Such instability is mostly related to the oscillatory eigenmodes, which are the fictitious internal modes of an interior space of a radiator for exterior problems [3,4] and the actual internal modes of an enclosure for interior problems [5,6]. Exponentially diverging instability happens usually in dealing with high order modes due to an increase of numerical approximation errors with an increase of their modal frequencies. It is known that TBEM calculation also involves with a monotonically increasing instability in time with a very slow rate [5–7], which is spatially constant. This instability is associated with the non-oscillatory eigenmode, which is a spatially constant internal response in an enclosure. The eigenvalues of non-oscillatory eigenmodes are usually very close to 1 due to both very small numerical error and physically scarce dissipation at very low frequencies. Because of this characteristic,

n

Corresponding author. Tel.: +82 42 350 3035; fax: +82 42 350 8220. E-mail addresses: [email protected], [email protected] (J.-G. Ih).

0022-460X/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2013.07.018

6464

H.-W. Jang, J.-G. Ih / Journal of Sound and Vibration 332 (2013) 6463–6471

TBEM response bears the monotonically increasing instability at the late times. Although the actual non-oscillatory static mode is physically negligible due to the presence of various dissipation mechanisms, most of interior problems with impedance boundary conditions suffer the instability due to the pseudo non-oscillatory eigenmodes [6]. It is because most impedance conditions can be regarded as a nearly rigid condition at very low frequencies due to its finite but very large magnitude. In this paper, the causes of monotonically increasing instability in TBEM calculation for the interior problems are investigated. The nonphysical nature of non-oscillatory eigenmodes are explained by using the analytically calculated solution of the linear wave equation for a sound field when the point source is given in an interior space. One can also ascribe the unstable behavior of non-oscillatory eigenmode to both large magnitude of involved eigenvalues and very low frequency noise in an input signal. Such two investigated causes of the instability are corrected by combined use of filtering method [6] and spectral decomposition method [7]. It should be noted that, if one adopts only one method from the two previously suggested methods, the instability problem cannot be solved without the degradation of accuracy. Transient sound fields within a rigid box and a lined box are used as demonstaration examples. 2. Basic theory 2.1. Kirchhoff integral and its discretization The acoustic TBEM is based on the following Kirchhoff integral as [1]:    Z  1 ðr0 rs Þ U n Rs ∂p 1 ∂p ðrs ; t ret Þ  ðrs ; t ret Þ dS; cðr0 Þpðr0 ; t 0 Þ ¼ pðrs ; t ret Þ þ 2 Rs Rs ∂n c0 ∂t s Rs

(1)

where c(r0) denotes the solid angle, p(r, t) the pressure at location r and at time t, Rs the distance between the field point at r0 and the surface point at rs (Rs ¼ |r0  rs|), t0 the current time for calculation, tret the retarded time (tret ¼t0  Rs/c0), c0 the speed of sound, S the domain boundary, n the unit normal vector toward the acoustic domain of concern, and (∂/∂n) the normal derivative on the surface. When all surface points are taken as field points, the discrete Kirchhoff integral becomes m

m

i¼0

i¼0

∑ ψi Pni ¼ ∑ φi Uni ;

(2)

where Pn  i and Un  i are the surface pressure and velocity vectors at time (n  i)Δt, Δt is a time step size, and m is the maximum retarded time steps related to the longest distance between boundary nodes as well as the numerical time delay due to both temporal discretization and interpolation. The surface velocity vector Un  i comes from ∂p/∂n in Eq. (1), which is replaced by a time derivative of surface normal velocity using the Euler equation. The coefficient matrices, ψi and φi, represent the contributions from dipole and monopole sources on the surface. For a problem involved in all types of boundary conditions, including an impedance boundary condition, the surface variable vectors can be restructured according to the type of boundary conditions [6]. The time domain impedance boundary condition is described by using the following infinite impulse response (IIR) digital filter expression: na

nb

l¼0

l¼0

∑ al uðnlÞ ¼ ∑ bl pðnlÞ :

(3)

Here, al and bl are the filter coefficients, and na and nb denote the filter orders. There are a great variety of impedance functions in the IIR filter type depending on the properties of wall materials, baking conditions, and so on. For example, when the same impedance data is assigned to all impedance boundary nodes, the final TBEM formulation can be expressed as N

na þm

i¼0

i¼0

∑ Ai Y ni ¼ ∑ Bi Xni :

(4)

Here, Xn i and Yn i are the restructured input and output vector, Ai and Bi are the coefficient matrices, and N is sum of the maximum retard time steps m and the highest filter order of the IIR filter type impedance: N¼max(na,nb)+m. Eq. (4) shows that TBEM calculation adapts the recursive, discrete convolution sum of a multi-input multi-output (MIMO) IIR system [6]. Although such a system calculates a time response directly, the responses often suffers the exponentially diverging instability as well as the monotonically increasing instability, which is related to the non-oscillatory static mode describing a spatially constant response in an enclosure. 2.2. Eigen-analysis The final MIMO IIR formulation of Eq. (4) can be rewritten as [8] 〈Y〉n ¼ M〈Y〉n1 þ 〈X〉n ;

(5)

H.-W. Jang, J.-G. Ih / Journal of Sound and Vibration 332 (2013) 6463–6471

6465

where 2

3 Yn 6 Y n1 7 6 7 6 7 7; ⋮ 〈Y〉n ¼ 6 6 7 6 7 ⋮ 4 5 YnðN1Þ

2 6 6 6 M ¼ 6 6 6 4

A0 1 A1

A0 1 A2

I

0

0 ⋮

I 0

0

⋯

A0 1 AN

::: ⋯

0

0 :::

⋯ 0

⋮ ⋮

0

I

0

3

2

7 7 7 7; 7 7 5

6 6 6 〈X〉n ¼ 6 6 6 4

1

⋯ 0 ⋮ ⋮ 0

1

T 3 7 7 n þm 7 a 7 ∑ A0 1 Bi Xni : 7 7 i¼0 5

(6a–c)

The output response 〈Y〉n represents the sequence of the surface pressure vectors during the maximum time delay, and the input 〈X〉n is given by boundary condition. The eigensolutions of the single iterative matrix M determine the characteristics of natural radiation modes. The eigenvector describes the time-marching radiation mode shape in complex domain, and the eigenvalue includes the information of both modal frequency and modal dissipation [4,6]. When some of eigenmodes have their eigenvalues outside a unit circle due to uncertain numerical approximation errors, TBEM calculation suffers the exponentially diverging instability. Based on the results of eigen-analysis, the sequential time response can be expanded into the eigenspace and its complementary space as [6] " # 〈Y〉n ¼

Ne

∑ Di ðnΔt Þ U〈u〉i þ ðI〈U〉〈U〉þ Þ〈Y〉n :

(7)

i¼1

Here, Ne is the number of eigensolutions, which can be calculated by multiplying the number of time steps of the maximum time delay and the number of boundary nodes, Di (nΔt) denotes the ith modal coefficient at a discrete time nΔh, (…)+ represents a Moore–Penrose pseudo inverse, and 〈u〉i is the ith eigenvector. 3. Non-oscillatory static mode 3.1. Analytic time-domain modal solution To investigate the characteristics of the non-oscillatory static mode, the interior sound field calculated by solving the following non-homogenous linear wave equation is examined: ∇2 pðr; tÞ

1 ∂2 pðr; tÞ ¼ f ðtÞδðjrqjÞ: c0 2 ∂t 2

(8)

Because the Kirchhoff integral equation is an integral form of the linear wave equation, the non-oscillatory eigenmodes of TBEM calculation exhibit the same characteristics with that of the linear wave equation. Here, q denotes a location of the point source and f(t) is an input function on itself. Based on the modal expansion theory, an interior sound field can be presented in time domain as [9] pðr; tÞ ¼ ∑i Di ðtÞΨ i ðrÞ:

(9)

Here, Di(t) denotes the ith modal coefficient at the continuous time t, of which discretized version is Di(nΔt) in Eq. (7). The mode shapes Ψi(r) satisfies the given boundary conditions and the homogeneous Helmholtz equation, which are generally orthogonal to each other. Because of a perfect modeling of the sound radiation in the continuous domain, the complementary components to the natural modes are not included. Meanwhile, the exterior sound field cannot be expressed by such modal expansion theory in the physical sense because the volume of concern is infinite so that there is no refection wave. However, mathematically, the Kirchhoff integral as well as Helmholtz integral has fictitious internal modes, which cause the well-known non-uniqueness problem. Substituting Eq. (9) into Eq. (8) and decomposing into the orthogonal mode set, one can obtain ordinary differential equations as ( ) 2 1 d Di ðtÞ Ψ ðqÞ 2 ki Di ðtÞ þ 2 ¼R i f ðtÞ; (10) Ψ c0 dt 2 V i Ψ i dV where ki is the wavenumber of the ith mode. Then, by solving the above equations, the modal coefficients of the nonoscillatory mode as well as oscillatory ones can be calculated. For the interior sound field with rigid boundaries, the modal coefficients can be expressed as Z t 2 _ i ð0Þt þ Rc0 Ψ i ðqÞ τf ðtτÞdτ for ki ¼ 0; (11a) Di ðtÞ ¼ Di ð0Þ þ D V Ψ i Ψ i dV 0 Di ðtÞ ¼

sin ðωi tÞ _ c0 2 Ψ i ðqÞ Di ð0Þ þ cos ðωi tÞDi ð0Þ þ R ωi ωi V Ψ i Ψ i dV

Z 0

t



 sin ðωi τÞf ðt; zÞ dτ

for

ki a0;

(11b)

_ i ð0Þ are an initial modal component and its initial time derivative, where ωi is the ith modal frequency, and Di(0) and D respectively. Because a sound field is perfectly reactive, imaginary parts of ki and ωi are zero. Eq. (11a) shows the nonoscillatory modal coefficient, which is determined by convoluting an input signal from a point source with a monotonically

6466

H.-W. Jang, J.-G. Ih / Journal of Sound and Vibration 332 (2013) 6463–6471

varying function τ. This linear function is the modal coefficient of an impulse response function of the sound field produced by the Dirac-delta function on a point source. When a sinusoidal function is given as an input signal, the part of convolution can be expressed as Z t t sin ðωq tÞ τ sin ðωq ðtτÞÞdτ ¼  ; (12) ωq ω2q 0 where ωq is a frequency of a sinusoidal source function. It shows that the monotonically increasing rate becomes steeper as the excited frequency becomes lower. After a point source excitation drops off, the non-oscillatory modal component does not stop to increase, but increases in a constant rate as Z t 0 Z t0  2 _ i ð0Þt þ Rc0 Ψ i ðqÞ Di ðtÞ ¼ Di ð0Þ þ D f ðτÞdτ t τf ðτÞdτ for f ðτÞ ¼ 0 when τ 4 t 0 : (13) 0 0 V Ψ i Ψ i dV This constant rate is a mean amplitude of an input signal, i.e., the component at 0 Hz. The last term in RHS represents DC offset, which is not contributed from 0 Hz component, but very low frequency input component. Such DC offset becomes bigger as the dominant frequency becomes lower as mentioned in Eq. (12). Eqs. (12) and (13) reveal that the non-oscillatory static mode is very sensitive to very low frequency and 0 Hz components of an input signal. On the other hand, the oscillatory modes hardly influence on the monotonically varying response because its impulse response is a sinusoidal function as shown in Eq. (11b). The interior sound field enclosed by dissipative boundaries calculated from the linear wave equation also includes the non-oscillatory mode. It is because the acoustic surface impedance at 0 Hz should have an infinite magnitude regardless of the type of boundary condition. This means that the surface velocity should always be zero at 0 Hz. Otherwise, an acoustic particle on the surface is forced to have a rigid body motion. One should note that such a phenomenon is definitely improbable. Thus, the monotonically varying instability can occur in TBEM calculation even for interior problems with dissipative boundary conditions. Meanwhile, in the numerical view point, the time domain impedance expression adopted in this study, i.e., IIR digital filter in Eq. (3), hardly realizes an infinite magnitude of impedance at 0 Hz. For such finite but very large magnitude of impedance at 0 Hz as well as at very low frequencies, the analytic solution of linear wave equation can include the pseudo non-oscillatory mode, of which modal coefficient can be written as Z t 2 sinhðαi τÞ _ i ð0Þ sinhðαi tÞ þ Rc0 Ψ i ðqÞ Di ðtÞ ¼ Di ð0Þcoshðαi tÞ þ D f ðtτÞdτ for αi ⪡1; ωi ¼ 0: (14) αi αi V Ψ i Ψ i dV 0 The pseudo non-oscillatory mode has very small attenuation constant αi only, and its modal coefficient of an impulse response is sinh(αiτ)/αi. This hyperbolic sine function is the solution of Eq. (10) for purely imaginary wavenumber ki ¼j  αi/c0 and the Dirac-delta function δ(t) on a point source. The modal coefficient of an impulse response in Eq. (14) can be approximated as a linear function τ in the early times, i.e., αiτ⪡1. Thus, in the same manner with the non-oscillatory mode for rigid boundaries, the pseudo non-oscillatory mode is also very sensitive to the input excitation at very low frequencies and 0 Hz. However, the response diverges exponentially in the late times due to the exponentially increasing behavior of the hyperbolic sine function. The aforementioned time domain analytic solutions of linear wave equation, in this paper, are used as references to evaluate the performance of stabilized TBEM calculation. However, their monotonically increasing instability due to 0 Hz input component is definitely unphysical. The non-oscillatory modal component should be stopped to increase after an impulse source is turned off. Such unphysical behavior disobeys the perturbation theory, which is the basic assumption of linear or pseudo-linear waves. The monotonically increasing response is not small enough compared to the static quantity. When a point source with 0 Hz is considered, this point source is gradually expanding with a constant velocity and eventually becomes an infinite size. Such monotonically increasing instability at 0 Hz input is physically implausible, but it might happen in the theoretical calculation due to the mathematical nature in computing Kirchhoff integral and the linear wave equation as well. In the same vein, very low frequency input components on the point source, which causes DC offset as mentioned in Eqs. (12) and (13), represent very large displacement on the surface of an infinitesimal point source due to the 2nd order temporal integration procedures. Thus, very large DC offset is also against the assumption of the perturbation theory and the infinitesimal size of point source. On the other hand, the sufficiently small DC offset satisfying the aforementioned assumptions can be regarded as an increase of potential energy in the interior space. Mathematically, the monotonically increasing instability can be avoided by simply excluding 0 Hz component from the input signal. However, in the discrete TBEM calculation, 0 Hz input noise cannot perfectly be eliminated due to the presence of uncertain, very small truncation error. Furthermore, a finite number of eigenmodes including two non-oscillatory ones are not enough to accurately describe the extremely low frequency components because the sound radiation modeling in discrete domain is not perfect. From this reason, in numerical TBEM calculation, very low frequency input noise cannot cause DC offset but it causes the monotonically increasing instability. Therefore, even if the source excitation is stopped and the involved eigenvalues of non-oscillatory eigenmodes are less than 1, TBEM calculation can reveal the monotonically increasing instability. The study subject dealt with in this work is originated purely from the mathematical treatment in calculating the interior sound field, which includes a constantly growing response related to the non-oscillatory mode as mentioned in Eqs. (11a) and (14). The purpose of this study, therefore, is to investigate the nonphysical nature and to find a method to alleviate such mathematical phenomenon.

H.-W. Jang, J.-G. Ih / Journal of Sound and Vibration 332 (2013) 6463–6471

6467

3.2. Non-oscillatory static eigenmodes in TBEM calculation TBEM calculation uses two non-oscillatory eigenmodes to describe the arbitrary variation of non-oscillatory modal component. Two non-oscillatory eigenmodes have linearly independent eigenvectors, although they have the same spatially constant mode shapes. This is because their eigenvectors represent time marching mode shapes with slopes of different eigenvalues. When the input vector 〈X〉n are bounded, the coefficient of non-oscillatory static mode Dnon(nΔt) can be expressed as Dnon ðnΔtÞ ¼ D1;0 λ1 n þ D2;0 λ2 n :

(15)

Here, λ1 and λ2 denote the eigenvalues of non-oscillatory eigenmodes, and D1,0 and D2,0 are their modal coefficients at the initial time. Because of scarce dissipation and very small numerical approximation error in very low frequency range, the eigenvalue of two non-oscillatory eigenmodes are very close to 1. Thus, Eq. (15) can be approximated as the following linear function:

 fD1;0 lnλ1 þ D2;0 lnλ2 g Dnon ðnΔtÞ ¼ D1;0 ð1 þ nlnλ1 Þ þ D2;0 ð1 þ nlnλ2 Þ ¼ D1;0 þ D2;0 þ ðnΔt Þ: Δt

(16)

The temporal variation of non-oscillatory modal component is approximated as a linear function with an arbitrary slope and an arbitrary DC level. It also shows that after the input signal drops off, the non-oscillatory modal component Dnon(nΔt) varies monotonically with a constant rate as implied in Eq. (13). In a physical sense, the non-oscillatory static modes cannot persist in the presence of dissipation mechanisms, e.g., an impedance boundary condition, geometric divergence in an exterior problem, etc. However, the calculation of the normal derivative integral for an exterior problem includes the fictitious non-oscillatory static modes, which is the non-trivial solution of normal derivative integral [7]. Besides, most of interior problems with impedance boundary condition has two pseudo non-oscillatory static modes. It is because most impedance conditions can be regarded as a nearly rigid boundary condition at very low frequencies due to its finite but very large resistence. 3.3. Causes of monotonically increasing instability and their stabilization methods TBEM calculation frequently suffers the erroneously monotonic variation of non-oscillatory static mode. When the magnitudes of their eigenvalues are larger than 1 in a very small scale, the response shows an exponential growth with a very slow rate. This exponential divergence stems from the numerical error in calculating the boundary integral and is often called the monotonically increasing instability due to its very slow exponential slope and non-oscillatory behavior. For the same reason, this instability becomes significant at the late time steps in most cases. Another important cause of the monotonically increasing instability is the input noise at 0 Hz and very low frequencies, which should be differentiated from the exponential divergence of non-oscillatory modes caused by the numerical errors in calculating the boundary integral. Such monotonically increasing instability is originated from the mathematical nature of both Kirchhoff integral and linear wave equation as manifested in Eqs. (12) and (13). Due to this fact, even if the involved eigenvalues are less than 1, TBEM response can exhibit an erroneous monotonic variation at the late times. To differentiate between aforementioned two important causes, for calculating the time domain analytic solutions of linear wave equation in the next chapter, unwanted 0 Hz input noise is included although the instability is expected. On the aforementioned possible causes of the monotonically varying instability, the previously suggested stabilization methods need to be considered to be applied for each case. For a simple example, a high pass filter can be used to avoid the monotonically varying instability by suppressing all low frequency components including 0 Hz [5]. However, because a high pass filter suppresses even the components of non-oscillatory static mode, it should be used only when the physically monotonic variation related to very low frequencies is not involved in. After projecting the time response onto the eigenspace as written in Eq. (7), the coefficients of non-oscillatory eigenmodes can be filtered as follows [6]: " # Ne

〈Y〉n ¼ ∑ ðF i Di ðnΔt ÞÞ〈u〉i þ I〈U〉〈U〉þ 〈Y〉n : (17) i¼1

Here, Fi denotes a filter coefficient, which plays a role as the damping applied to the ith eigenmode. Then, the sequential response at the next time step becomes " # Ne

〈Y〉nþ1 ¼ ∑ ðF i λi Di ðnΔt ÞÞ〈u〉i þ M I〈U〉〈U〉þ 〈Y〉n : (18) i¼1

In comparison with Eq. (7), the magnitude of the ith eigenvalue is replaced by |Fi∙λi|, which represents the corrected magnitude change rate per single time step. By setting all the replaced eigenvalues |F∙λ| to be less than 1, the exponential divergence due to unstable poles can be eliminated. To stabilize the monotonically increasing instability, the replaced eigenvalues of two non-oscillatory eigenmodes can be set as 1. In this case, the rate of monotonic variation of non-oscillatory static mode becomes zero. This method does not influence the other oscillatory modes at all, and it can also solve all possible causes of the monotonically increasing instability as like the high pass filter. However, because this method filters

6468

H.-W. Jang, J.-G. Ih / Journal of Sound and Vibration 332 (2013) 6463–6471

even the non-oscillatory static mode, it should be used when the non-oscillatory static modes are not physically involved in or their components become constant at the late times. To suppress the exponential divergence of non-oscillatory mode, a stabilization method based on the spectral decomposition can be adopted [7]. Because the non-oscillatory static mode describes an ideally constant sound field in time and space, it is the non-trivial solution of the following homogeneous matrix equation, which can be derived from Eq. (4): ! N

∑ Ai Y ni ¼ ½MYni ¼ 0:

(19)

i¼0

Theoretically, the non-oscillatory static mode cannot be included in the eigensolutions of the matrix [M]. However, due to numerical approximation errors, [M] has a few numbers of eigensolutions, of which eigenvectors are nearly constant and eigenvalues are very small compared to the remaining eigenvalues. Such eigensolutions are directly related to the unstable poles of non-oscillatory eigenmodes. The real symmetric matrix [M] can be spectrally decomposed into the orthogonal projections of its eigenvectors. To obtain a stable solution, the projected non-oscillatory eigenvectors should be excluded. In the recursive TBEM calculation, this exclusion is applied to the coefficient matrix at the current time step, A0. This method has little influence on the oscillatory modes, and it can delay the onset time of the exponential growth of non-oscillatory modes because their eigenvalues become closer to a unit circle. Also, the monotonic variation of an inherent response can be described even after using this method, because the modified eigenvalues of non-oscillatory eigenmodes are still different. 4. Examples and discussions Transient sound fields within a rigid box and a lined box are calculated as demonstration examples for the discussion on non-oscillatory eigenmodes and effect of possible stabilization methods. The geometry of the box is 0.5  0.7  0.8 m3, for which boundary element model is composed of 264 nodes and 524 linear triangular elements. The maximum element size is 0.141 m, which limits the high frequency for reliable calculation as 410 Hz under λ/6 criterion. The time step size Δt is selected as 0.2 ms (fs ¼5 kHz, c0Δt/Δh¼0.69). The maximum retarded time of the longest diagonal of the box is 17∙Δt. An octave-band impulse, which is the inverse Fourier transform of the acoustic excitation of the octave band pass filter centered at 63 Hz, is given at a point at (0.10, 0.15, 0.20) m. Table 1 summarizes the characteristics of the lowest three internal modes including non-oscillatory static mode for a sound field in a rigid box. In this table, the index (Nx,Ny,Nz) represents the mode shape, and two (0, 0, 0) modes denote the non-oscillatory eigenmodes. Here, f and D are the modal frequency and decay rate, respectively, which can be calculated from the eigenvalue λ as [6]. D ¼ ð20=ΔtÞlog 10 ðjλjÞ;

f ¼ ð1=2πÞðargðλÞ=ΔtÞ

(20a,b)

In the rigid box example, the unstable oscillatory eigenmodes related to the exponentially diverging instability in very high frequency range are not included. In spite of the perfectly reactive characteristics, the low order modes are dissipative in a certain extent due to increasing numerical approximation errors with frequency. Two non-oscillatory eigenmodes are scarcely dissipative, and they have monotonically increasing and decreasing characteristics with the same absolute rates. In spite of their scarce dissipation, one of them has the unstable pole outside the unit circle on a very small scale, which causes the exponential divergence with a very slow rate. By using the spectral decomposition method [7], the unstable non-oscillatory eigen mode can be nearly stabilized, as the decay rates of non-oscillatory eigenmodes become 2.1  10  10 dB/s. Moreover, the modal characteristics of other oscillatory modes are hardly influenced. Two non-oscillatory eigen-modes become the complex conjugate to each other, nearly at 0 Hz with the same magnitudes of eigenvalues. Although the arbitrary variation of non-oscillatory static mode can be described mathematically by using two such eigen-modes, very small ratio of the imaginary part to its real part can be a cause of the numerical errors in the increase rate and the DC value of non-oscillatory modal components. Fig. 1 compares the calculated time responses with the analytically calculated response by using modal expansion method in Section 3.1. Comparison is made at a point (0.40, 0.10, 0.60) m in a rigid box excited by an octave-band impulse centered at 63 Hz. For an accurate analytic solution, the lowest 1000 modes below 3.3 kHz are used in the calculation, because all interior modes of an enclosure with rigid boundaries are physically non-dissipative. One can find that the

Table 1 Characteristics of the lowest 3 eigenmodes for a sound field in a rigid box. Mode (Nx,Ny,Nz)

Analytic f (Hz)

TBEM D (dB/s)

f (Hz)

(0,0,0)

0.0

0

0.0

(0,0,1) (0,1,0)

214.7 245.4

0 0

215.6 246.4

TBEM using the spectral decomposition method [7] |λ|

D (dB/s) 5

1+3.2  10 1  3.2  10  5 0.9994 0.9992

 1.4 +1.4 26.8 33.3

|λ|

f (Hz) 6

3.2  10  3.2  10  6 215.6 246.4

D (dB/s)  15

1+4.8  10 1+4.8  10  15 0.9994 0.9992

 2.1  10  10  2.1  10  10 26.8 33.4

H.-W. Jang, J.-G. Ih / Journal of Sound and Vibration 332 (2013) 6463–6471

6469

Fig. 1. A comparison of the analytic solution with the calculated interior field pressure in a rigid box excited by an octave band impulse centered at 63 Hz: ————, analytic solution; ——, TBEM before using stabilization algorithms; —□—, TBEM stabilized by the spectral decomposition method [7]; —●—, TBEM stabilized by using the filtering method [6] in addition to the spectral decomposition method [7]. (a) Early times and (b) late times, 0.2  20 s.

Fig. 2. A comparison of the analytic solution with the frequency spectra of the calculated early time response below 0.2 s in a rigid box: ————, analytic solution; ——, TBEM before using stabilization algorithms; —□—, TBEM stabilized by the spectral decomposition method [7]; —●—, TBEM stabilized by using the filtering method [6] in addition to the spectral decomposition method [7].

numerically calculated responses in early time steps (tr0.2 s) agree well with the analytic response, of which the average difference norm is 0.7 percent. Frequency spectra are also compared in the Fig. 2, in which their dominant components below 0.3 kHz also exhibit very small difference. Although discrepancy becomes large with an increase of frequency, its contribution to the total difference is quite small. However, in the late time steps, the exponential divergence with a very slow increase rate is shown in the result if the stabilization algorithms are not employed. The analytic solution also shows the monotonic variation with a constant rate in the late times due to non-zero mean amplitude of the input signal as mentioned in Eq. (13), while its energy ratio of non-zero mean amplitude to the total input energy is  75.5 dB. It should be noted that such very small 0 Hz input component causes the instability mathematically as mentioned in Section 3.1. One can find that the exponential divergence problem with a very slow rate is solved by applying the stabilization method based on the spectral decomposition [7]. However, the monotonic variation due to non-zero mean amplitude of an input signal persists yet, and its increase rate is not accurately predicted due to very small ratio of the imaginary part of non-oscillatory eigensolution to its real part. As the effective remedy against the monotonic variation due to unwanted 0 Hz input noise, one should additionally use the high pass filter [5] or the filtering method based on the eigen-analysis [6] as shown in Fig. 1. For the example of a lined box, an absorptive liner with a rigid backing is assumed on the left side at z ¼0.8 m, while a rigid condition is assumed on the other sides. An empirical model [10] is used to obtain the frequency domain surface impedance of a 2 cm-thick absorptive layer, of which the flow resistivity is set to 20,000 Rayls/m. The frequency domain impedance is fitted to the form of IIR digital filter in Eq. (3) with filter coefficients as a0 ¼4380, a1 ¼  7362, a2 ¼3707, a3 ¼  586.1, b0 ¼ 1, b1 ¼  2.346, b2 ¼1.800, and b3 ¼  0.453. In Fig. 3, the given frequency domain impedance and its fitted time domain impedance are presented. The impedance at 0 Hz, which is the ratio of the sum of filter coefficients bl to that of filter coefficients am, is 1.389  105 Ns/m3, and the corresponding reflection coefficient is 0.994; one should note that the impedance should have an infinite magnitude physically.

6470

H.-W. Jang, J.-G. Ih / Journal of Sound and Vibration 332 (2013) 6463–6471

Fig. 3. A comparison of time-domain surface impedance for an absorptive liner with a rigid backing condition: ——, IIR digital filter; ———— the frequency domain surface impedance from an empirical model [10]. (a) Resistance and (b) reactance.

Table 2 Characteristics of the lowest 3 eigenmodes for a sound field in a lined box. Mode (Nx,Ny,Nz)

Analytic f (Hz)

TBEM D (dB/s)

f (Hz)

(0,0,0)

0

1.4

0.0

(0,0,1) (0,1,0)

212.5 244.5

80.7 45.8

214.2 245.8

TBEM using the spectral decomposition method [7] |λ|

D (dB/s)

f (Hz)

|λ|

D (dB/s)

1+2.2  10  5 1  4.5  10  5 0.9981 0.9987

 1.0 +1.9 81.1 58.2

0.0

1+1.5  10  10 1 2.2  10  05 0.9981 0.9987

 6.5  10  6 +1.0 81.1 58.2

214.2 245.8

Fig. 4. A comparison of the analytic solution with the calculated interior field pressure in a lined box by TBEM: ————, analytic solution; ——, TBEM before using stabilization algorithms; —□—, TBEM stabilized by the spectral decomposition method [7]; —●—, TBEM stabilized by using the filtering method [6] in addition to the spectral decomposition method [7]. (a) Early times and (b) late times, 0.2  20 s.

The characteristics of the lowest three modes for the sound field in a lined box are summarized in Table 2. Because of very large but finite impedance at 0 Hz corresponding to a decay rate of 1.4 dB/s, the rate of monotonically decreasing nonoscillatory static mode is a little steeper than that of monotonically increasing mode. Here, the frequency response for the fitted IIR filter type impedance is 1.07  106 N/m3 s at 0 Hz. After using the spectral decomposition method [7], the exponentially increase rate of the non-oscillatory eigenmode becomes nearly zero as its decay rate is  6.5  10  6 dB/s. Also, the remaining oscillatory modes are little influenced similar to the rigid box case. On the other hand, monotonically decreasing mode is still dissipative due to the finite impedance at 0 Hz. Fig. 4 compares the analytic solution with the calculated responses before and after employing stabilization algorithms for TBEM. The calculated response in the early time steps (tr0.2 s) agrees well with analytical solution, of which the average difference norm is 2.0 percent. However, it diverges exponentially in a very slow rate at the late time steps due to the presence of monotonically increasing pseudo non-oscillatory static mode. In the calculation of analytic solution, the lowest 9 modes below 0.5 kHz are used because high order modes are much less influential on the low frequency response than those of rigid enclosure. By employing the spectral decomposition method [7], the unstable non-oscillatory eigenmode

H.-W. Jang, J.-G. Ih / Journal of Sound and Vibration 332 (2013) 6463–6471

6471

Fig. 5. A comparison of the analytic solution with the frequency spectrum of the calculated early time response below 0.2 s in a lined box: ————, analytic solution; ——, TBEM before using stabilization algorithms; —□—, TBEM stabilized by the spectral decomposition method [7]; —●—, TBEM stabilized by using the filtering method [6] in addition to the spectral decomposition method [7].

having its eigenvalue a little larger than 1 can be nearly stabilized. However, the monotonically varying instability due to the low frequency input noise persists similar to the rigid box example. In particular, the calculated response does not converge to zero even though the energy dissipation on the impedance boundary condition is included. By using the filtering method [6], this remaining instability can be fully stabilized as shown in Fig. 4. Also, their frequency spectra in Fig. 5 show good agreements with each other. 5. Conclusions In this paper, nonphysical characteristics of the non-oscillatory static mode and its monotonically increasing instability in TBEM calculation for the interior acoustic problem are studied. In TBEM calculation, the arbitrary variation of non-oscillatory modal component is approximated as a linear function by using two non-oscillatory eigenmodes with different eigenvalues. It is asserted that the monotonically increasing instability is caused by both unstable poles and very low frequency input noise. Calculation of interior problems with impedance boundary condition also reveals a similar monotonically increasing instability due to the presence of pseudo non-oscillatory eigenmodes. In the reconsideration of the previously suggested stabilization algorithms, it is noted that the spectral decomposition method [7] can stabilize only the exponentially diverging instability caused by the unstable poles. It is concluded that, for the complete stabilization of the computed result, the filtering method [6] should be used in additional to the spectral decomposition method. Acknowledgments This work was partially supported by BK21Project and the NRF Grants (nos.2012R1A2A2A01009874 and 2010-0028680). References [1] P.H.L. Groenenboom, C.A. Brebbia (Ed.), Wave propagation phenomena, Vol.2, Progress in Boundary Element Methods, Pentech, London1983, pp. 24–52. [2] W.J. Mansur, C.A. Brebbia, Formulation of boundary element method for transient problems governed by the scalar wave equation, Applied Mathematical Modelling 6 (1982) 307–311. [3] P.D. Smith, Instabilities in time marching methods for scattering: cause and rectification, Electromagnetics 10 (1990) 439–451. [4] H. Wang, D.J. Henwood, P.J. Harris, R. Chakrabarti, Concerning the cause of instability in time stepping boundary element methods applied to the exterior acoustic problem, Journal of Sound and Vibration 305 (2007) 289–297. [5] M. Stutz, M. Moser, M. Ochmann, Instability problems using the time domain BEM for impulse response calculations, Proceedings of the 6th Forum Acusticum, Aalborg, Denmark, 2011, pp. 253–257. [6] H.-W. Jang, J.-G. Ih, Stabilization of time-domain acoustic boundary element method for the interior problem with impedance boundary conditions, Journal of Acoustical Society of America 131 (2012) 2742–2752. [7] M. Parot, C. Thirard, C. Puillet, Elimination of a non-oscillatory instability in a retarded potential integral equation, Engineering Analysis with Boundary Elements 31 (2007) 133–151. [8] S.J. Dodson, S.P. Walker, M.J. Bluck, Implicitness and stability of time domain integral equation scattering analyses, Applied Computational Electromagnetics 13 (1998) 291–301. [9] R. Alimenti, P. Mezzanotte, L. Roselli, R. Sorrentino, Efficient analysis of waveguide components by FDTD combined with time domain modal expansion, IEEE Microwave and Guided Wave Letters 5 (1995) 351–353. [10] Y. Miki, Acoustical properties of porous materials—modifications of Delany–Bazley models, Journal of Acoustical Society of Japan 11 (1990) 19–24.