Document Sample

Simulation of transient sound radiation using the Time Domain Boundary Element Method M. St¨tz1 , M. Ochmann2 u 1 u Beuth Hochschule f¨r Technik Berlin, 13353 Berlin, Germany, Email: stuetz@tfh-berlin.de 2 u Beuth Hochschule f¨r Technik Berlin, 13353 Berlin, Germany, Email: ochmann@tfh-berlin.de Introduction retarded time. The boundary surface Γ is divided into N planar elements with a uniform spatial pressure and ﬂux The boundary element method (BEM) is a widely used distribution. For q a constant approach and for p a linear numerical tool. Calculations may be performed in the approach in time are chosen. A more detailed description frequency domain (FD) or in the time domain (TD). The can be found in [10] and [13]. With these approaches one fundamentals of the TD-BEM were developed in 1960 can derive the following discretized integral equation from Friedman and Shaw [1] and later in 1968 from Cruse and Rizzo [2]. The FD-BEM is limited to the simulation N i of stationary processes. Instationary processes like 1 n − 2πpi (x0 ) = q Ψ(tri )dΓn moving sources, engine acceleration, impulses etc. must n=1 m=1 Γn r m be simulated in the time domain. Mansur [3] developed ∂r 1 one of the ﬁrst boundary element formulations in the time + (i − m + 1)pn − (i − m)pn m m−1 Ψ(tri )dΓ domain for the scalar wave equation and for elastody- Γn ∂n r2 namics with zero initial conditions. The extension of this (3) formulation to non-zero initial conditions was presented a by Antes [5]. Later J¨ger [6] and Antes and Baaran [7] extended the time domain formulations to analyze 3D 1 ; tri ∈ [tm−1 , tm ) noise radiation caused by moving sources. Since then this Ψ(tri ) = (4) 0 ; else method is subject to continuous enhancement, whereby however the frequency range algorithms exhibit a clear Using the Collocation method the following correspond- lead in development. ing linear system of equations can be derived The main reason that the TD-BEM-Method is not being commercially used till now is the instable behavior that µmax may appear. Ergin et al [8], Chappell et al [9] and ˆ − 2π pi = ˆ Gµ qi−µ+1 u St¨tz and Ochmann [10] showed numerical evidence that µ=1 the instable behavior is caused by the eigenfrequencies µmax of the structure. Ergin and Chappell used the Burton- + p p Hµ [µˆi−µ+1 − (µ − 1)ˆi−µ ] . (5) Miller-Method [11] to prevent instabilities. In this paper µ=1 a diﬀerent approach, the use of the so-called CHIEF Method, will be presented. Gµ and Hµ are square matrices describing the system at a relative time step µ = i − m + 1. The matrix entries Numerical model are calculated in the following way A detailed derivation of the boundary integral equation µ 1 used here can be found e.g. in Meise [4]. gab = Ψ(trµ )dΓ (6) rab Γb q ∂r p 1 ∂p ∂rab 1 4πd(x0 ) p(x0 , t0 ) = − + + dΓ hµ = r ∂n r2 cr ∂t ab 2 Ψ(trµ )dΓ . ∂n rab (7) ret Γ Γb (1) p denotes the sound pressure, r the distance between Compared to the FD-BEM the resulting matrices the observation point and source point and c the speed are very sparse, because for each matrix with index of sound. Γ is the boundary surface and q is the ﬂux, µ the integration must only be performed within following the relationship derived from Eulers equation the integration limits of an inner radius (µ − 1)c∆t ∂vn and an outer radius µc∆t. A diﬃculty that arises q=− 0 (2) from that condition is that one must not integrate ∂t over the whole element but over the part that lays The index ret in Eq 1 states that all variables are within the integration limits. The sparsity of the taken with respect to the retarded time tret = t − r .c matrices is explained due to the fact that most Discretizing time itself in equidistant time steps ti = elements are not within the integration limits and i∆t with (i = 1, 2, ...) we get tri = i∆t − r for the c therefore have a zero entry. Eq. 5 can be written as integration area and element is becoming much smaller p1 ˆ G1 0 ··· ··· 0 ˆ q1 than the element itself. The integrals are evaluated ˆ p2 G2 G1 0 ··· 0 q2 ˆ numerically using a Gauss Quadratur. By increasing the . . . − 2π . . = G3 G2 G1 . . . . + numbers of Gauss points the accuracy can be increased which will decrease the inﬂuence of the eigenfrequencies . . .. . . . . . . . . to a certain degree, as can be seen in Figure 3. The ˆ pi 0 ··· Gµmax · · · G1 ˆ qi inﬂuence does not vanish but the aberration from the H1 0 ··· ··· 0 ˆ p1 analytic solution does not exceed 1 dB. 2H2 H1 0 ··· 0 p2 ˆ . . . . 3H3 − H2 2H2 H1 . . 90 90 1 1 . . .. . 120 60 120 60 . 0.8 0.8 . . . 0.6 0.6 150 30 150 30 0 ··· −(µmax − 1)Hµmax · · · H1 ˆ pi 0.4 0.4 0.2 0.2 (8) 180 0 180 0 to solve the time response at once. This results in 2 210 330 210 330 Block-Toeplitz-Matrices that are not easy to handle because of the relatively big size of (number of elements 240 270 300 240 270 300 x time steps)2 . Thats why the time response is solved in an iterative way, starting with solving the ﬁrst row of (a) ∆t=1.1 ms β=1.9 (b) ∆t=0.38 ms β=0.66 ˆ ˆ ˆ Eq. 8 (−2π p1 = G1 q1 + H1 p1 ) and proceeding till the last row. This approach is numerically much easier to Figure 1: Example of eigenvalues within unit circle of system handle. operator T for diﬀerent time step sizes (a) and (b) Accuracy and stability 95 A common problem of the TD-BEM reported in the literature is the occurrence of spurious instabilities. In 90 case of instability the results start to oscillate at high p [dB] frequency with exponentially rising amplitude. To look 0.45 at the stability of the method the iterativ solving process 0.61 β 85 0.82 can be written as 1.03 1.29 pi = T pi−1 + φi (9) analytical 80 with T being the system operator. If all eigenvalues of T 85 170 255 340 425 510 595 680 f [Hz] lay within the unit circle the system can be considered as stable, because p can not grow with proceeding time Figure 2: Inﬂuence of eigenfrequencies for diﬀerent β on if φi = 0. sound radiation of test sphere ) We are now looking at the simple example of a sphere 95 with r = 1m consisting of 384 ﬂat rectangular ele- ments. Figure 1 shows that for this given example the eigenvalues lay all within the unit circle. But with 90 decreasing time step the absolut values of the eigenvalues p [dB] are increasing, making the system less stable. 3 Using the sphere to simulate a monopol the results are 85 6 12 stable as shown in Figure 2. But apparently the sound 24 radiated by the sphere is inﬂuenced by the eigenfrequen- analytisch cies of the structure. This connection however is still 80 85 170 255 340 425 510 595 680 mathematically unproven. In Figures 2, 3 and 5 the f [Hz] eigenfrequencies are marked by the vertical dotted lines. Figure 3: Inﬂuence of eigenfrequencies depends on accuracy The inﬂuence of the eigenfrequencies depends strongly of integration, results for diﬀerent numbers of Gauss Points on β, which describes the ratio of time step size ∆t and are shown element size h. β = c∆t/h (10) CHIEF Method Small changes of β can change the inﬂuence of the The problem of the natural frequencies is well known eigenfrequencies substantially. Nevertheless in none in frequency domain BEM calculation, also called the of the examined cases an unstable behavior, like an nonuniqueness problem. Various methods have been de- exponential rise of the sound pressure, could be observed. rived to ﬁnd a modiﬁed equation with a valid solution for With decreasing β the evaluation of the matrix entries all frequencies. Well known approaches are the combined is getting more complicated, because the intersection of ﬁeld integral equation by Burton and Miller [11] and the combined Helmholtz integral equation formulation (CHIEF) by Schenk [12]. For time domain calculations ﬁrst Ergin [8] in 1999 and later Chappell et al [9] in 2006 used a Burton-Miller approach to avoid the above collocation point chief point described problem. The principle of the CHIEF method is very simple. As shown in section “Numerical model” a linear system of equations can be derived using the collocation method. Figure 4: CHIEF points are located inside of the structure The collocation points are located on the boundary and collocation points are located on the surface surface. Placing some additional points inside of the structure we can derive a second linear system of equa- 95 tions. But in this case we set the pressure in this points to zero to force the inner sound ﬁeld to zero. 90 cmax p [dB] 0= Gc qi−µ+1 µˆ 3 6 µ=s 85 12 cmax 24 c + p p Hµ [µˆi−µ+1 − (µ − 1)ˆi−µ ] (11) analytical CHIEF µ=s 80 85 170 255 340 425 510 595 680 765 850 935 f [Hz] If we combine both linear systems Eq. 5 and Eq. 11 we get an overdetermined system. Figure 5: The use of CHIEF points prevents the inﬂuence As an example here the resulting matrix of eigenfrequencies. c containing all Hµ and Hµ matrices is shown The analytical far ﬁeld approximation for a circular piston in an inﬁnite wall can be found in the literature: H1 0 ··· ··· 0 2H2 H1 0 ··· 0 e−ikr J1 (kRsinα) . . p(r, α) = iω υn R2 · (13) 3H3 − H2 2H2 H1 . r kRsinα . .. . . . Figure 6 shows the test results for 2 immission points. 0 · · · −(µmax − 1)Hµmax · · · H1 ........................................................... Both points are located in 7.5m distance, P1 with c hight h=1.2m and P2 with h=3m as can be seen in sHs 0 ··· ··· 0 (s + 1)Hs+1 − (s − 1)Hs sHs c c c 0 ··· 0 Figure 7. The results show a very good agreement with . . .. . the analytical solution. Also the directivity which can . 0 ··· c −(se − 1)Hse c · · · sHs be observed in point P2 is reproduced quite well. The (12) diﬀerence in the low frequency region is expected. Short circuit eﬀects of the numerical model without wall are The matrix containing the Gµ and Gc matrices is set up µ not present in the analytical model with an inﬁnite wall. in a similar way. To the authors knowledge the use of the A second test case examined the sound radiation of a CHIEF method in time domain BEM calculations was railway wheel. The impulse response function in the time never shown before. This may have diﬀerent reasons. domain was calculated by a FEM calculation. Based One reason is that the CHIEF method does not work on these data the sound radiation of the wheel could using the iterativ solution process. Applying the method be calculated. A comparison to a FD-BEM calculation to the non-iterativ one-step solution makes it necessary performed with Virtual Lab, shows a good agreement to solve the huge rectangular matrix (Eq. 12). The (Figure 8). There is no measurement data available, so no solution process is not an easy task because of the conclusion about the simulation accuracy can be made. immense size of the matrices involved. To deal with this problem, the special symmetry of the matrices must be used. Using an iterativ solver that can handle rectangular matrices, like the LSQR method, only one Conclusion line of the matrix has to be build at a time. The results The TD-BEM is a promising method for the computation in Figure 5 show that the CHIEF method successfully of the transient sound radiation. Numerical test cases regularizes the solution at the inner eigenfrequencies of show the reliability and stability of this method. The the structure. accurate evaluation of the boundary integrals plays a fun- damental role, because errors in the integration process Test cases - circular piston and seem to stimulate the inﬂuence of the eigenfrequencies of the structure. In order to ensure accurate simulation railway wheel results, it was shown that the CHIEF method can be A good test case to see if this method can simulate the used to regularize the eigenfrequencies. More research sound radiation of railway wheels is a circular piston. regarding the eﬃciency and the limits of the CHIEF 80 numerisch P1 TD−BEM analytisch P1 0 70 FD−BEM numerisch P2 p [dB] 60 analytisch P2 −20 50 −40 p [dB] 40 0 500 1000 1500 2000 2500 3000 3500 30 f [Hz] sound power density 20 0 10 −20 0 −40 0 500 1000 1500 2000 2500 3000 f [Hz] −60 −80 0 500 1000 1500 2000 2500 3000 3500 Figure 6: Sound radiated by a circular piston (in free ﬁeld), f [Hz] compared to far ﬁeld approximation (in inﬁnite wall). Figure 8: Sound radiation of an railway wheel computed with TD-BEM and V-Lab (FD-BEM). sungswesen der TU Carolo-Wilhelmina zu Braun- schweig, Germany; 1993 [7] Antes H, Baaran J. Noise radiation from moving surfaces. Engineering Analysis with Boundary Elements; Volume 25, Issue 9, October 2001, Pages 725-740 Combustion [8] A.A.Ergin, B. Shanker, E. Michielssen: Analysis of Figure 7: Setup of circular piston test case with 2 immission transient wave scattering from rigid bodies using a points located in 7.5m distance, P1 with hight h=1.2m and Burton-Miller approach, JASA, 106, 1999, pp. 2396- P2 with h=3m. 2404 [9] D.J.Chappell, PJ. Harris, D. Henwood, R. method in TD-BEM calculations has to be done. The Chakrabarti: A stable boundary element method results will be presented in one of the next publications. for modeling transient acoustic radiation,JASA,120, 2006, pp. 74-80 References [10] M.St¨tz, M.Ochmann, Stability behaviour and u [1] M. Friedmann, R. Shaw: Diﬀraction of pulses by results of a transient boundary element method cylindical Obstacles of Arbitrary cross Section, J. for exterior radiation problems, Acoustics’08-Paris, Appl. Mech.,Vol. 29 (1962), S.40-46 2008 [2] T.A. Cruse, F.J. Rizzo: A direct formulation [11] A. J. Burton, G. F. Miller: The Application of and numerical solution of the general transient Integral Equation Methods to the Numerical So- elastodynamic problem I, J. Math. Analysis and lution of Some Exterior Boundary-Value Problems, Applic., Vol.22 (1968), S.244-259 Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, Vol. 323, [3] Mansur, W. J.: A Time-Stepping technique No. 1553, A Discussion on Numerical Analysis of to solve wave propagation problems using the Partial Diﬀerential Equations (Jun. 8, 1971), pp. boundary Element Method, PhD thesis, University 201-210 of Southampton, 1983 [12] Harry A. Schenck: Improved Integral Formulation [4] Meise T.: BEM calculation of scalar wave for Acoustic Radiation Problems, JASA , July 1968, propagation in 3-D frequency and time domain (in Vol 44, Issue 1, pp. 41-58 German). PhD Thesis, Technical Reports Nr.90-6, a u Fakult¨t f¨r Bauingenieurwesen, Ruhr-Universit¨t, a u [13] M.St¨tz, M.Ochmann : Calculation of the acoustic Bochum, Germany; 1990 radiation of an open turbulent ﬂame with a transient boundary element method, International Congress [5] H. Antes, A Boundary Element Procedure for on Acoustics in Madrid (ICA), 2007 Transient Wave Propagations in Two-dimensional Isotropic Elastic Media, Finite Elements in Analysis and Design 1 (1985) 313-322 a [6] J¨ger M. Entwicklung eines eﬃzienten Randele- u mentverfahrens f¨r bewegte Schallquellen. PhD u Thesis. Fachbereich f¨r Bauingenieur- und Vermes-

DOCUMENT INFO

Shared By:

Categories:

Tags:

Stats:

views: | 5 |

posted: | 11/25/2011 |

language: | English |

pages: | 4 |

OTHER DOCS BY chenmeixiu

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.