Use of plane wave basis for the numerical solutions of acoustic wave problems using element based methods has become an attractive approach for extending the allowable frequency range for simulations beyond that available using piecewise polynomial elements. The non-uniqueness of the solution at characteristic frequencies resulting from the use of the conventional boundary integral equation is well known. The standard methods of overcoming this problem are the so-called CHIEF method and that of Burton and Miller. The latter method introduces a hypersingular integral which can be treated in several ways. In this paper we present results for Partition of Unity BEM (PUBEM) for Helmholtz problem and compare the performance of CHIEF against a Burton-Miller formulation regularised using the approach of Chen et al

Coates, G.

Diwan, G.C.

Trevelyan, J.

English

Durham Research Online

Proceedings of the 9th UK Conference on Boundary Integral Methods, University of Aberdeen,UK, 8-9th July 2013COMPARISON OF CHIEF AND BURTON-MILLER APPROACHES IN COL-LOCATION PARTITION OF UNITY BEM FOR HELMHOLTZ PROBLEMSG.C. DIWAN, J. TREVELYAN and G. COATESSchool of Engineering and Computing Sciences, Durham University, Durham DH1 4JR, UKe-mail: g.c.diwan@durham.ac.ukAbstract. Use of plane wave basis for the numerical solutions of acoustic wave problems usingelement based methods has become an attractive approach for extending the allowable fre-quency range for simulations beyond that available using piecewise polynomial elements. Thenon-uniqueness of the solution at characteristic frequencies resulting from the use of the con-ventional boundary integral equation is well known. The standard methods of overcoming thisproblem are the so-called CHIEF method and that of Burton and Miller. The latter methodintroduces a hypersingular integral which can be treated in several ways. In this paper wepresent results for Partition of Unity BEM (PUBEM) for Helmholtz problem and compare theperformance of CHIEF against a Burton-Miller formulation regularised using the approach ofChen et al.1. INTRODUCTIONThe theory of the boundary element method (BEM) for solving boundary integral equations(BIE) is well established. It is known that the Conventional BIE (CBIE) based on Green’sfunction representations for an exterior acoustic problem results in a non-unique solution atcharacteristic frequencies for the corresponding interior problem and that this is a purely math-ematical phenomenon. There are several methods to handle the non-uniqueness. One of theextensively used method is the so called Combined Helmholtz Integral Equation Formulation(CHIEF) due to Schenck [1], where some extra Helmholtz integral equations evaluated at inte-rior points are added in the original system matrix. Although this results in an over-determinedsystem, the CHIEF method ensures a unique solution at the characteristic frequency. How-ever, one needs to choose the interior points such that they do not lie on the nodal lines of theinterior modes of Helmholtz problem thus failing to provide any necessary constraint for unique-ness of the solution. This can therefore introduce uncertainties for complicated geometries athigh frequencies as the nodal lines become densely packed in the interior making it difficultto find suitable locations for the placement of interior points. Another method to avoid thenon-uniqueness problem is due to Burton and Miller [2]. It was shown in their work that theintegral equation resulting from linear combination of the CBIE and its normal derivative atthe collocation point always results in a unique solution. The main problem with this method isthe evaluation of the hypersingular integral which arises as a result of the differentiation of theCBIE. In the present study, we compare the CHIEF method with one regularized form of theBurton-Miller formulation for acoustic scattering from hard cylinders in two dimensions usingPUBEM. The two methods are compared for their accuracy of the solution and efficiency.2. BOUNDARY INTEGRAL EQUATIONThe mathematical formulation for deriving the CBIE from the Helmholtz equation is well estab-lished [3]. The CBIE for an acoustic scattering (or radiation) problem governed by the Helmholtzdifferential equation is given byc(p)ϕ(p) +∫Γ∂G∂nqϕ(q)dΓ(q) =∫ΓG∂ϕ(q)∂nqdΓ(q) + ϕi(p) (1)where p is the collocation or source point, q the field point, G the free space Green’s function forthe Helmholtz problem, nq and np the normals respectively at points q and p pointing away fromacoustic domain Ω, ϕ(q) the unknown acoustic potential and ϕi(p) the known incident acousticwave. c(p) is the free coefficient which depends on the local geometry of Γ at p. In this studywe assume Γ is smooth and take c(p) = 12 . The normal derivative of (1) at the collocation pointp is given byc(p)∂ϕ(p)∂np+∫Γ∂2G∂np∂nqϕ(q)dΓ(q) =∫Γ∂G∂np∂ϕ(q)∂nqdΓ(q) +∂ϕi(p)∂np(2)and the Combined Hypersingular BIE (CHBIE) due to Burton and Miller [2] isc(p)ϕ(p) + αc(p)∂ϕ(p)∂np+∫Γ∂G∂nqϕ(q)dΓ(q) + α∫Γ∂2G∂np∂nqϕ(q)dΓ(q) =∫ΓG∂ϕ(q)∂nqdΓ(q) + α∫Γ∂G∂np∂ϕ(q)∂nqdΓ(q) + ϕi(p) + α∂ϕi(p)∂np(3)where α is a coupling constant most commonly taken as i/k. In the present study, we analysethe acoustic scattering from sound hard cylinders. A sound hard surface is where the normalderivative of the total acoustic potential vanishes. Therefore, all the terms involving the normalderivative of acoustic potential vanish. As mentioned earlier, the main drawback of (3) is thenumerical treatment of the hypersingular integral, i.e. the last integral on the left hand side.Chen et al [4] give the following weakly singular form of the hypersingular integral∫Γ∂2G∂np∂nqϕ(q)dΓ(q) =∫Γ[∂2G∂np∂nq− ∂2G0∂np∂nq]ϕ(q)Γ(q) +∫Γ[ϕ(q)− ϕ(p)−∇ϕ(p) · (q − p)] ∂2G0∂np∂nqdΓ(q)(4)+∫Γ∇ϕ(p) · nq∂G0∂npdΓ(q)− 12∇ϕ(p) · npwhere G0 is the free space Green’s function for the Laplace equation. For the present case ofa hard boundary, the last term in the right hand side of (4) vanishes. Consequently, the finalequation for this case of a hard boundary can be expanded asc(p)ϕ(p) +∫Γ∂G∂nqϕ(q)dΓ(q) + α∫Γ[∂2G∂np∂nq− ∂2G0∂np∂nq]ϕ(q)Γ(q)+ (5)α∫Γ[ϕ(q)− ϕ(p)−∇ϕ(p) · (q − p)] ∂2G0∂np∂nqdΓ(q) + α∫Γ∇ϕ(p) · nq∂G0∂npdΓ(q) = ϕi(p) + α∂ϕi(p)∂npThe acoustic potential at a point x on the boundary Γ using plane wave basis can be approxi-mated asϕ(x) =3∑j=1NjMj∑mAjmeikdjm·x x ∈ Γ (6)where Nj is the jth shape function, Ajm the unknown which can be thought of as the amplitudeof the mth plane wave with wave number k associated with node j. The direction of the mthplane wave at node j is given by unit vector djm and x is the location of the point where thepotential ϕ is sought. We consider a three noded element with Mj as the number of plane wavesassociated with the jth node. In the context of the BEM, the plane wave basis defined in (6) canbe used to express the unknown acoustic potential on the boundary Γ. Choosing appropriatelocations on the boundary Γ as collocation point p yields the following set of linear equationsfrom (5)[H]{a} = {b} (7)where the vector a contains the amplitudes of plane waves, Ajm, which can be used to quicklyrecover the acoustic potential on the boundary Γ using (6).3. NUMERICAL EXAMPLESThis section presents the error analyses for two problems viz. i) plane wave scattering from asingle sound hard cylinder and ii) from an array of four cylinders. At the outset, it is convenientto define a parameter τ which gives the number of degrees of freedom per wavelength for a givenproblem, i.e.,τ =Tka(8)where T is the total number of degrees of freedom in the system for one cylinder and a is theradius of the cylinder. Thus for the problem of scattering from a single cylinder with unit radius,τ = T/k where T will be simply the multiplication of the total number of nodes on the scattererboundary and number of plane waves per node. It may be noted that the introduction of planewaves into the basis makes the boundary integrals oscillatory in nature and it is very importantthat these integrals be computed as accurately as possible. We follow the strategy described in[5], namely element subdivision, for evaluating these integrals.For all the results presented here the parameter τ ≈ 3.0 unless otherwise mentioned. Thisvalue has been found to be sufficient to recover solutions with acceptable engineering accuracyof 1% and moderate condition numbers which can be efficiently handled with the Singular ValueDecomposition (SVD) algorithm, see [5]. All the results are obtained with 30 integration (Gauss)points per wavelength unless otherwise mentioned. For both the single cylinder and four cylinderexamples, we use two 3-noded continuous elements per cylinder along with the trigonometricshape functions presented by Peake et al [6]. For all computations the integration points areplaced analytically on the scatterer boundary. We now define the relative L2 error for the totalacoustic potential ϕ on the boundary Γ, E2(ϕ) asE2(ϕ) =∥ϕ− ϕ̃∥∥ϕ̃∥(9)where ϕ is the numerically computed solution and ϕ̃ the analytical solution computed using theinfinite or approximate series for a given scattering problem. The 2-norm condition number forthe matrix H, κ(H) may be defined asκ(H) =σmax(H)σmin(H)(10)where σmax(H) and σmin(H) are respectively the maximum and minimum singular values of thematrix H computed using the SVD algorithm. As discussed earlier, the placement of interiorcollocation points for the CHIEF method might become an issue. For the numerical examplespresented in this study, the interior points are placed completely randomly in the interior of thecylinder(s). The number of interior points used here is 20% of the total number of equationsin (7) since this has been found to give stable results for the CHIEF method. Also the CHIEFpoints in the interior of the cylinder(s) are placed such that they are sufficiently away from theboundary.40 60 80 100 120 140 16010−1210−1010−810−610−410−2kaE2 (φ)  BM with unsplit Telles transformationBM with split Telles transformationCHIEF with split Telles transformationFigure 1: E2(ϕ) for the single cylinder problem3.1 Scattering from a single sound hard cylinderWe first investigate the performance of CHIEF and Burton-Miller methods for the classicalproblem of plane wave scattering from an acoustically hard cylinder of infinite extent. We usethe infinite series in [7] to compute the analytical solution for the scattered potential, ϕs, on thesurface of a hard cylinder. The total acoustic potential ϕ can be computed by simply performinga complex addition of incident wave to the scattered potential obtained using the infinite seriesfrom [7], i.e., ϕ = ϕi+ϕs. The relative L2 error for the total acoustic potential is then computedusing (9). Fig. 1 shows the relative L2 error, E2(ϕ) for CHIEF and Burton-Miller methods.As seen from Fig.1, CHIEF provides better accuracy compared to Burton-Miller results. Notethat when the weak singularity in (5) is handled with the Telles scheme without splitting theinterval containing the singularity, the Burton-Miller formulation gives poorer results. Despitethe regularization, the integrals in the (5) are slowly converging. As is evident, to achieve a com-parable accuracy to that of the CHIEF method, the integrals in the regularized Burton-Millerformulation needs to use Telles transformation with splitting the interval towards left and theright of the singularity.3.2 Scattering from an array of four cylindersThe scattering from a multi-cylinder array presents a more challenging case as it involves multiplereflections from individual cylinders which ultimately forms the total acoustic field. The recursivemultiple reflections make this problem an ideal candidate to test the efficacy of PUBEM toobtain an accurate solution. We consider a setting of four unit radius sound hard cylinders ofinfinite extent with their centres placed at (-2,-2), (2,-2), (2,2) and (-2,2) in a two dimensionalhomogeneous unbounded acoustic medium (air). Let this array then be impinged by a unitamplitude plane wave with wavenumber k at an angle of θI = 45◦ with the horizontal. We usethe formula proposed by Linton and Evans [8] (eq. 2.15) to compare our PUBEM solution forthe total acoustic potential on the surface of each cylinder. The formula proposed by Lintonand Evans is based on the addition theorem that combines the separable solutions of Helmholtzequation, see [9] for details. The addition theorems can be efficiently used to compute thesolution but the infinite series has to be truncated in practice. Theoretically of course, aninfinite sum should result in a converged solution. However, when solving even the truncatedsystem of linear equations, the addition of extra terms in the series can make the matrix formed0 20 40 60 80 10010101020103010401050Number of terms in the Linton−Evans seriesCondition number of the system matrix  k = 2.4048k = 10k = 20Figure 2: Stability of Linton-Evans series, eq. (2.15) in [8].Cyl. 1 Cyl. 2 Cyl. 3 Cyl. 4CHIEF BM CHIEF BM CHIEF BM CHIEF BMτ ≈ 3.0 5.98E-04 1.03 2.15E-04 2.71 2.78E-04 6.26 1.67E-04 2.01τ ≈ 3.5 8.67E-06 4.06E-05 1.28E-05 9.39E-05 7.6E-06 3.29E-05 516E-06 3.82E-05τ ≈ 3.9 2.01E-07 4.07E-06 1.87E-07 4.01E-06 3.46E-07 6.32E-06 1.92E-07 5.17E-06Table 1: PUBEM results - E2(ϕ) for scattering from four cylinder array for k = 36.9171 andθI = 45◦, 100 terms in Linton-Evans seriesusing (2.15) in [8] highly ill-conditioned. Fig.2 shows the dependence of the condition numberof the system matrix formed from (2.15) in [8] on the number of terms included in the series.Note that k = 2.4048 is an irregular wavenumber (first zero of the first kind Bessel function, J0).Clearly the reason for such significantly high condition numbers is the wide spread of eigenvalueswith the growing number of terms in the series.We useM terms in the Linton-Evans series, thus, a system of linear equations of sizeNc(2M+1) is formed where Nc is the number of cylinders (4 in the present case). We use a linear leastsquares solver with QR factorisation to solve this system of linear equations using suitableroutines from the LAPACK library and obtain the total acoustic potential on each cylindersurface. This solution is considered as the reference solution and used to compute the relativeL2 error (see (9)) for our PUBEM solution with the CHIEF and Burton-Miller methods. For theerror analysis of the four cylinder problem, we consider two cases of the wavenumber, namely,k = 36.9171 and k = 150. It may be noted that both k = 36.9171 and k = 150 are irregularwavenumbers. The L2 error results shown in Tables 1-2 are obtained using two continuouselements per cylinder with trigonometric shape functions as before. All the results are obtainedwith 30 integration points per wavelength.We have used M = 100 for k = 36.9171 and M =200 for k = 150, in the Linton-Evans series.The condition number of the coefficient matrix for the Linton-Evans series for k = 36.9171 with100 terms was 14.28 and that for k = 150 with 200 terms was 16.29. It can be noted fromTables 1-2 that the accuracy of both CHIEF and regularised Burton-Miller methods improveswith more plane waves per node i.e. by increasing the value of the parameter τ .Cyl. 1 Cyl. 2 Cyl. 3 Cyl. 4CHIEF BM CHIEF BM CHIEF BM CHIEF BMτ ≈ 2.2 8.26 15.26 12.53 27.67 42.43 98.36 15.68 37.90τ ≈ 2.6 6.73E-05 3.0E-03 7.30E-05 2.6E-03 7.23E-05 4.8E-03 7.98E-05 7.0E-03τ ≈ 3.0 6.48E-05 6.47E-05 6.40E-05 6.39E-05 6.68E-05 6.70E-05 6.40E-05 6.46E-05Table 2: PUBEM results - E2(ϕ) for scattering from four cylinder array for k = 150 and θI = 45◦,200 terms in Linton-Evans series.4 CONCLUSIONWe have presented a plane wave enriched BEM formulation of the regularised Burton-Millerequations for the exterior acoustic scattering problem in two dimensions. The error analysespresented for the classical single and the multiple scattering problems show that the CHIEFmethod outperforms Burton-Miller method by at least 1 order of magnitude for the problemsconsidered in this paper. The Burton-Miller method can prove competitive despite the difficultand slowly converging integrals if suitable coordinate transformations are implemented. TheCHIEF method may be preferred over the Burton-Miller formulation, at least for simpler ge-ometries and moderate wavenumbers (k < 200) as the former does not have the problem ofhypersingular integrals and provided that a sufficient number of interior collocation points arechosen that ensure the linear independence of the coefficient matrix H.AcknowledgementsThe first author gratefully acknowledges generous financial support from Durham University.REFERENCES1. H.A. Schenck (1968) Improved Integral Formulation for Acoustic Radiation Problems,Journal of the Acoustical Society of America, 44, pp.41-58.2. A.J. Burton and G.F. Miller (1971) The application of integral equation methods to thenumerical solution of some exterior boundary-value problems, Proceedings of the RoyalSociety of London. Series A, Mathematical and Physical Sciences, 323, pp.201-210.3. R.D. Ciskowski and C.A. Brebbia (1991) Boundary Element Methods in Acoustics, CMPand Elsevier Applied Science.4. K. Chen, J. Cheng and P.J. Harris (2009) A new study of the Burton and Miller methodfor the solution of a 3D Helmholtz problem, IMA Journal of Applied Mathematics, 74(2),pp.163-177.5. E. Perrry-Debain, O. Laghrouche, P. Bettess and J. Trevelyan (2004) Plane-wave basisfinite elements and boundary elements for three-dimensional wave scattering, Philosoph-ical Transactions of the Royal Society of London. Series A: Mathematical, Physical andEngineering Sciences, 362, pp.561-577.6. M.J. Peake, J. Trevelyan and G. Coates (2012) Novel basis functions for the partitionof unity boundary element method for Helmholtz problems, International Journal forNumerical Methods in Engineering, 93, pp.905-918.7. J.J. Bowman, T.B.A. Senior and P.L.E. Uslenghi (1970) Electromagnetic and AcousticScattering by Simple Shapes, North-Holland Pub. Co.8. C.M. Linton and D.V. Evans (1990) The interaction of waves with arrays of vertical circularcylinders,Journal of Fluid Mechanics, 215, pp.549-569.9. P.A. Martin (1985) Integral-equation methods for multiple-scattering problems I. Acous-tics, The Quarterly Journal of Mechanics and Applied Mathematics, 38, pp.105-118.

Comparison of CHIEF and Burton-Miller approaches in collocation Partition of Unity BEM for Helmholtz problems

https://durham-repository.worktribe.com/file/1152730/1/Accepted%20Conference%20Proceeding

Comparison of CHIEF and Burton-Miller approaches in collocation Partition of Unity BEM for Helmholtz problems

Abstract

Similar works

Full text

Available Versions

Durham Research Online