- Email: [email protected]
The locally corrected Nyström method applied to 3D scalar SIE in acoustic cavities using curvilinear coordinates
Engineering Analysis with Boundary Elements 79 (2017) 110–118
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
The locally corrected Nyström method applied to 3D scalar SIE in acoustic cavities using curvilinear coordinates
MARK
⁎
F. Villa-Villaa, , H. Pérez-Aguilarb, A. Mendoza-Suárezb a b
Centro de Investigaciones en Óptica, Loma del Bosque 115, Lomas del Campestre 37150, León, Guanajuato, Mexico Facultad de Ciencias Físico Matemáticas, Universidad Michoacana de San Nicolás de Hidalgo, Av. Francisco J. Mújica S/N 58030, Morelia, Mich., Mexico
A R T I C L E I N F O
A BS T RAC T
Keywords: Integral equations Helmholtz equation Acoustic scattering
In this work we make a comparison of the results obtained by the Locally Corrected Nyström method and a parametric method previously proposed by the authors to solve the scalar surface integral equation in order to determine the resonant modes in 3D acoustic cavities. Both methods are set up in basis of a parametric representation in curvilinear orthogonal coordinates. So, any geometry that involves quadrics or higher order surfaces can be considered. Different cavities whose modes can be analytically determined for the purpose of comparison of results are considered. It is found that for a spherical cavity it gives the proper results locating the modes as predicted analytically while for the cubical cavity, the Locally Corrected Nyström method gives some imprecision in the calculation of the modes despite using high order numerical quadrature, probably due to the presence of edges where derivatives of the basis vectors are not defined.
1. Introduction In the study of scattering of acoustic and electromagnetic waves by objects of complex geometry, the method of moments (MoM) is known to be a good resource to solve the surface integral equations (SIE) involved in these problems [1–4]. So, for many years people has devoted some effort to develop this method specially using basis functions defined on surface elements that involve triangular facets. Even though almost arbitrary and complex objects can be represented by triangle meshing, it is known that as the size of the objects increase, some error problems can arise probably due to the representation of curved surfaces by flat triangular surface elements. When a surface is discretized into triangles, this meshing process presents itself a complex problem since MoM requires to have a precise account of the triangles of the partition and how they are coupled in the basis functions. In this regard specialized commercial and free software has been implemented but the intervention of the user is always necessary to solve details on meshing and interfacing the MoM's algorithms to the output of this software. An additional problem with MoM is the handling of the complex and lengthy mathematical expressions, resulting from a sampling process that involves nested double surface integrals which must be solved numerically and/or analytically to implement the numerical algorithm [2]. In more recent years, the Locally Corrected Nyström method (LCN) has been proposed as an alternative to MoM [5,6]. In this method
⁎
neither basis functions nor sampling process are necessary, since the integrals over the surface elements are directly calculated by using Gauss-Legendre quadrature rules when the argument of the Green's function and its derivatives are not in the vicinity of their singularities. The algorithm of LCN is far more simple to manage and implement, despite the integrals with singular arguments deserve a special treatment. However, a similar situation holds in this respect with MoM or any other method. For some years the standard procedure to treat integrals with singular kernels in LCN has consisted in designing highly specialized quadrature rules [6,7]. This procedure has the limitation that it is not possible to find a quadrature of integrals whose kernels present high order singularities like those appearing in the dyadic Green's function of the electromagnetic case. To lower the order of the singularities in the integrals, a procedure resembling the integration by parts is applied taking advantage of its convolutional double argument [6]. This procedure introduces an external differential operator to the integrals, complicating the numerical calculations. Since most of the research has been focused on the use of plane triangles (LCN and MoM) and low order basis functions in electromagnetics and acoustics as well, now people are concentrating on developing methods based on the use of curvilinear patches and high order basis functions [8]. In the case of MoM this problem is cumbersome [8], however, LCN seems to be more suited to curvilinear formulation.
Corresponding author. E-mail address: [email protected] (F. Villa-Villa).
http://dx.doi.org/10.1016/j.enganabound.2017.04.003 Received 21 September 2016; Received in revised form 25 January 2017; Accepted 15 April 2017 0955-7997/ © 2017 Elsevier Ltd. All rights reserved.
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
In this work we propose the solution of the SIE in acoustics by our own version of LCN using a parametric representation of the arguments in the surface integrals in a curvilinear orthogonal system. This l , tu , t v , formulation involves in a natural way the orthonormal basis n composed by the normal and tangent vectors in the coordinate directions (u , v ) of a parametric surface. To solve the integrals when the argument is in the vicinity of the singularities of the Green's function and its normal derivative, we propose the use of a Taylor series expansion to treat the diagonal elements by giving an explicit procedure that allows us to obtain analytic matrix elements. So, we do not require any interpolation process or the derivation of highly specialized quadrature rules that depend on the order of the singularities. The Taylor series expansion, besides the basis vectors involves the first derivatives of the tangent vectors ∂t u/∂u , ∂t v/∂v , which are related to the local curvatures of the surface. The detailed outline of this work reads as follows: In Section 2 we give a brief description of the integral equation related to the Helmholtz equation and the boundary conditions in acoustics. In Section 3 the Locally corrected Nyström method applied to curvilinear parametric surfaces by using the singularity extraction to treat the singular integrals is outlined. In Section 4 to validate the method, the application of LCN to determine resonant and Bloch modes in different physical systems like cavities and three-dimensional phononic crystal (3DPC) is considered, and finally, we conclude our work in Section 5 with an overview on the analytical and numerical results of this work.
Fig. 1. System with a curvilinear geometry composed of two regions (interior and exterior). The observer's position is indicated by r , while r′ denotes the source position on the surface for the integral surface equation. This vector spans the surface under integration.
ψ (r′) = f (r′), r′ ∈ S′, (Dirichlet ) ∂ψ = g (r′), r′ ∈ S′, (Neumann ) ∂n′
2. Theory As a starting point let us consider the SIE for the potential function ψ (r′) of the displacement vector (or the pressure as a function of position) in an acoustic medium
⎡
∫S′ ⎢⎣G (r, r′) ∂ψ∂n(r′ ′)
− ψ (r′)
∂G (r, r′) ⎤ ′ ⎥ dS = ψ (r) θ (r) − ψ inc (r), ∂n′ ⎦
are given, where f (r′) and g (r′) are arbitrary regular functions that define the corresponding boundary condition. 3. The Nyström method
(1)
where
⎧1 if r ∈ V1 θ (r) = ⎨ ⎩ 0 if r ∉ V1,
Let S′ be a closed surface that limits the volume V1, with no incident exciting field ψ inc (r) = 0 as shown in Fig. 1. By assuming that the normal points outward to S′ and considering the convention that the surface integral is positive if the normal points outward and negative in the opposite case. It is customary to validate a numerical method by calculating the scattered field by simple objects and compare with analytical results basically for a sphere. Since analytical results are also available for resonant modes in cavities with simple geometry, we shall be interested in the fields and their characteristic modes inside S′ given the Dirichlet or Neumann boundary conditions with f (r′) = 0 or g (r′) = 0 to test our method. As we will calculate the integrals by applying a limit when the observer approaches from the exterior of S′ along the normal, this implies that θ (r) = 0 . If we divide S′ into small surface elements Si like quadrilateral curvilinear patches, the integral equation for this region will be
(2)
is the Heaviside step function, and
G (r, r′) =
e qR , 4πR
(3)
represents Green's function for the 3D space. We have called q=ik, to further simply expressions where it appears repeatedley. ψ inc (r) represents an incident field for this exterior problem, R = r − r′ being r the observer's position and r′ the integration vector that spans the surface S′ (see Fig. 1), and
∂G (r, r′) l ·∇′G (r, r′) =n ∂n′
(4)
l·R =γ (r, r′) n
(5)
l is the outward is the normal derivative of the Green's function, n normal to the surface at r′, and
N
∑∫ i =1
γ (r, r′) = a (R−3 − qR−2 ) e qR .
(6)
1 ∂ψ (1) 1 ∂ψ (2) = , ρ1 ∂n′ ρ2 ∂n′
Si
G (r, r′)
∂ψ (r′) ′ dS − ∂n′
N
∑∫ i =1
Si
ψ (r′)
∂G (r, r′) ′ dS = 0. ∂n′
(9)
Although not explicitly indicated all the functions in the integrand depend on the properties of the acoustic medium of density ρ1, and wave number k1 = ω / v1 where v1 is the velocity of propagation of the perturbation. In this equation the Green's function and its normal derivative are known, and the potential ψ (r′) together with its normal derivative ∂ψ (r′)/∂n′ constitute the functions to be determined. Let us now suppose that we want to calculate the integrals over each patch Sn by using Gauss-Legendre quadrature
The potential and its normal derivative satisfy the boundary conditions at an interface given by
ψ (1) = ψ (2),
(8)
(7)
where ρ1 and ρ2 represent the densities of the corresponding media. Sometimes the conditions [9–12].
111
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al. N
Pi
∑ ∑ wj G (r, r′ j ) i =1 j =1
∂ψ (r′ j ) ∂n′
N
dS′ j −
Pi
∑ ∑ wj ψ (r′ j )
∂G (r, r′ j ) ∂n′
i =1 j =1
dS′ j = 0, (10)
where j = 1, 2, …, Pi are the points considered for the quadrature along the i-th surface element. These Pi limits of the summations can be different. Then, we get an equation with Q = 2N (P1 + P2+⋯+PN ) terms and the same number of unknown quantities. This equation has to be evaluated at the same points on the surface for the observer's position r = rm . Here m = 1, 2, …, N (P1 + P2+⋯+PN ), in this way we have a linear system of Q/2 equations with Q unknowns N
Pi
N
Pi
∑ ∑ wj G (rm, r′ j ) Φj dS′ j − ∑ ∑ wj ψj i =1 j =1
i =1 j =1
∂G (rm, r′ j ) ∂n′
dS′ j = 0. (11)
To simplify notation the double summations can be substituted by a simple one as Q /2
Q /2
∑ wk G (rm, r′n ) Φn dS′n − ∑ wn ψn n =1
n =1
∂G (rm, r′n ) ′ dS n = 0. ∂n′
(12) Fig. 2. Schematic description of the i-th curvilinear patch Si. The symmetric sub-patches σj are centered at the j-th point determined by the Gauss-Legendre quadrature rule in question.
If r → r′ the involved integrals in Eq. (11) can not be calculated by Gaussian quadrature. As we will state below, to regularize the problem when having functions with singularities of different orders, we use the well known technique of singularity extraction to separate the functions in one smooth, and one analytically integrable part that despite of being singular can be treated by taking a limit to determine their Cauchy principal value. This procedure has demonstrated to give good results when the surface is discretized by using triangle elements [13,14]. In this work we will demonstrate that the use of higher order surface elements is possible by means of curvilinear quadrilateral patches.
In analogy with the smooth and analytical parts of G (r, r′ j ), the functions γs (r, r′ j ) and γa (r, r′) stand for the smooth and analytical parts of the γ (r, r′ j ) function given in Eq. (6). These functions shall be obtained explicitly below when considering the singularity extraction technique. At this point it worth remembering that r′ j is the position vector of point j at the i-th patch. We have added a subindex s to the smooth parts and a to denote the analytical parts of the functions, and in the analytical integrals we have taken out the source functions
3.1. Local correction of the Nyström method
ψj = ψ (r′) r′= r j ,
When the argument of the Green's function and its normal derivative approaches to zero, that is when r → r′ the GaussLegendre quadrature no longer applies. So, a process of regularization is required to calculate the surface integrals by using a technique called singularity extraction. This method consists on separating the integrands in two parts: one smooth part Gs (r, r′) that can be integrated numerically by Gaussian quadrature and another that we call Ga (r, r′) that can be integrated analytically. Thus, the integrals over a curvilinear patch shall be
∫S
G (r, r′)
i
∂ψ (r′) ′ dS = ∂n′
∫S
i
∂ψ (r′) Gs (r, r′) dS′ + ∂n′
∫S
i
Pi
Φj =
∑
j =1 Pi ⎡ =∑ Φj ⎢wj Gs (r, r′ j ) dS′ j + ⎣ j =1
∫σ
j
∫σ
Ga (r, r′) dS′ j
⎤ Ga (r, r′) dS′⎥ ⎦ (13)
∫S
i
ψ (r′)
∂G (r, r′) ′ dS = ∂n′
∫S
l ·RdS′ ψ (r′) γs (r, r′) n
+
∫S
(16)
3.2. Integrands in the vicinity of the singularity r → r′: analytic functions To treat the integrals with analytic integrands, let us assume that the surface S′ is a parametric curvilinear surface dependent on two parameters (u , v ), where u, and v represent two orthogonal directions of the curvilinear system that represents the surface in question. If S′ is divided into Q small quadrilateral surface elements σj with centroid at (uj , vj ), so dS′ = Hdudv as shown in Fig. 2. In the expression of the differential of area, H represents the Jacobian of the transformation. The integrals over each curvilinear sub-patch σj can be expressed as
∂ψ (r′) Ga (r, r′) dS′ ∂n′ Φj
j =1
, r′= r j
by assuming that they are small enough to guarantee that these functions are almost constant in the integration sub-domains or subpatches σj. These sub-elements of area are centered at the points r′ j as indicated in Fig. 2.
Pi
≈∑ wj Φj Gs (r, r′ j ) dS′ j +
∂ψ (r′) ∂n′
(15)
i
∫σ
l ·RdS′ ψ (r′) γa (r, r′) n
i
Ga (r, r′) dS′ = j
vj +
Δvj
uj +
Δuj
∫v − Δv2 ∫u − Δu2 j
j
j
2
Ga (r, r′) Hdudv,
j
2
(17)
Pi
lj ·Rj dS′ j ≈∑ wj ψj γs (r, r′ j ) n
∫σ
j =1 Pi
+
∑ j =1
ψj
∫σ
l ·RdS′ γa (r, r′) n
⎡ lj ·Rj dS′ j + =∑ ψj ⎢wj γs (r, r′ j ) n ⎣ j =1
∫σ
j
j
vj +
Δvj
uj +
Δuj
∫v − Δv2 ∫u − Δu2 j
j
j
2
j
l ·RHdudv. γa (r, r′) n
2
(18)
If we consider the change of variables η = u − uj , and ξ = v − vj in Eqs. (17)–(18), these transform to
j
Pi
l ·RdS′ = γa (r, r′) n
⎤ l ·RdS′⎥ γa (r, r′) n ⎦
Δvj
∫σ
(14) 112
Ga (r, r′) dS′ = j
Δuj
∫− Δ2v ∫− Δ2u j
2
j
2
Ga [r, r′(η + uj , ξ + vj )] Hdηdξ
(19)
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al. Δuj
Δvj
∫S
l ·RdS′ = ψ (r′) γa (r, r′) n
n
∫− Δ2v ∫− Δ2u j
j
2
r′ ≈ r′ j + t uj η + t vj ξ + …
l γa [r, r′(η + uj , ξ + vj )] n
To calculate the integrals by taking their Cauchy principal value, let us define the observer's position as
2
(η + uj , ξ + vj )·R (η + uj , ξ + vj ) Hdηdξ.
(20)
r = r′ j + ζ nj ,
If we just recall that the position vector r′ denotes a point on the surface S′, by expanding into a Taylor's series around the point r′ j = r′(uj , vj ) r′(η + uj , ξ + vj ) = r′ j + 1 ∂ 2r′ 2 ∂ξ 2
∂r′ ∂η
η+ η =0, ξ =0
ξ2 η =0, ξ =0
∂r′ ∂ξ
ξ+ η =0, ξ =0
1 ∂ 2r′ + 2 ∂η∂ξ
1 ∂ 2r′ 2 ∂η2
η2 +
R = r − r′ ≈ ζ nj − ηt uj − ξ t vj .
(21)
∂r′ ∂η
∂r′ t vj = ∂ξ
= t u u = uj = v = vj
η =0, ξ =0
= t u u = uj v = vj
η =0, ξ =0
∂r′ ∂u
∂r′ = ∂v
u = uj v = vj
u = uj v = vj
1
R ≈ (ζ 2nj2 + tu2j η2 + tv2j ξ 2 ) 2 .
t′vj =
∂ 2r′ ∂u2 ∂ 2r′ ∂v 2
u = uj v = vj
u = uj v = vj
=
=
∂ 2r′ ∂η2
(22)
, (23)
∂ 2r′ ∂ξ 2
(24)
t′uvj
t′vuj =
u = uj v = vj
∂ 2r′ ∂v∂u
u = uj v = vj
η =0, ξ =0
∂ 2r′ = ∂η∂ξ
=
(32)
t v = t vj + t′vj ξ + t′vuj η ,
(33)
l = tu × t v n
(34)
lj + (t uj × t uvj + t′uj × t vj) η + (t uj × t′vj + t uvj × t vj) ξ ≈n
(35)
+(t′uj × t uvj) η2 + (t uvj × t′vj ) ξ 2,
(36)
. the Jacobian can be approximated as
(25)
H ≈ nn [1 + 0.5(t uj ·t′uj + t uvj ·t vj) η + 0.5(t vj ·t′vj + t uvj ·t uj) ξ ].
These vectors are illustrated in Fig. 3 for a spherical surface. We will define for short hand notation, the vectors
∂ 2r′ = ∂u∂v
t u = t uj + t′uj η + t′uvj ξ,
and an approach for the normal vector can be obtained to a first order
, η =0, ξ =0
(31)
If we now consider the expansion of the tangent vectors around rn up to first order,
,
and the curvature vectors, respectively
t′uj =
∂ 2r′ ∂ξ∂η
l ·R) H ≈ nj [ζnj2 − τnu η2 − τnv ξ 2], (n (26)
(38)
that appears in one of the integrands. In the last equation we have lj ·t′uj , and τnv = n lj ·t′vj to further simplify expressions. We called τnu = n also have omitted terms containing odd powers of η and ξ, in virtue that when integrated over symmetric limits the integrals are null.
. η =0, ξ =0
(37)
With these results, we can also calculate the product
, η =0, ξ =0
(30)
lj constitute a right where the tangent vectors t uj , t vj , and the normal n hand set orthonormal basis. With these considerations the magnitude of R , shall be
ξη +…, η =0, ξ =0
we can identify from the differential geometry the tangent vectors to the surface S′ at r′ j given by
t uj =
(29)
then, after integrating; the principal value will be that of the limit when ζ → 0+ [17]. In this way taking the difference of the last two equations we obtain
η =0, ξ =0
1 ∂ 2r′ ηξ + 2 ∂ξ ∂η η =0, ξ =0
(28)
(27)
Then, we can express Eq. (21) to a first order approximation 3.3. Singularity extraction The Taylor's expansion of the Green's function
G (R ) = a
exp(qR ) 1 = R 4π
∞
∑ m =0
q m m −1 R , m!
(39)
can be separated in two parts: one smooth, and one analytic [2].
G (R ) = Gs (R ) + Ga (R ),
(40)
where the smooth part is given by ∞
Gs (R ) = a
∑ m =1 m ≠2
q m m −1 R , m!
(41)
and the analytic part by
⎛1 q2 ⎞ Ga (R ) = a ⎜ + R⎟ . ⎝R 2! ⎠
(42)
In a similar fashion,
γ (r, r′) = γs (r, r′) + γa (r, r′), ∞
γs (r, r′) = a
∑ m =3 m ≠4
Fig. 3. Basis and curvature vectors on a spherical parametric surface.
113
1 − m m m −3 q R , m!
(43)
(44)
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
⎡1 q2 1 q4 ⎤ γa (r, r′) = a ⎢ 3 − − R⎥. ⎣R 2 R 8 ⎦
Iv−1 =
(45)
Multiplying the analytic functions Ga (r, r′), and γa (r, r′) by the other factors in the integrands
⎛1 q2 ⎞ Ga (r, r′) H = a ⎜ + R⎟ , ⎝R 2! ⎠
Iu1 =
Iv1 =
Ga (r, r′) HdA ≈ a I−1 + a
∫σ
i
q2 2
I1
4. Resonant modes in cubical and spherical cavities To test the proposed method, we shall compare its results with those obtained analytically for two highly symmetric cavities under the Dirichlet and Neumann boundary conditions [15]. One condition is related to perfectly soft surfaces and the other to perfectly hard surfaces. As stated in the previous sections, the complete system of SIE over a closed surface S′ will be
(49)
∫S′ G (r, r′) ∂ψ∂n(r′ ′) dS′ − ∫S′ ψ (r′) ∂G ∂(rn,′ r′) dS′ = 0,
where we have defined, the elementary surface integrals, Δuj 2 Δuj − 2
Ip = lim+
∫
Iup = lim+
∫− Δ2v ∫− Δ2u
ζ →0
∫
Ivp = lim+ ζ →0
j
∫
+
j
p tv2j ξ 2 ) 2 dηdξ,
Δuj 2 Δuj − 2
∫
Q /2
lim ζI−3 =
I−1 =
I1 =
Q /2
∑ L mn Φn − ∑ Nmn ψn = 0. n =1
(51)
(63)
n =1
In this expression we have called
(ζ 2nj2
+
tu2j η2
+
p tv2j ξ 2 ) 2 ξ 2dηdξ,
⎧ wn G (rm, r′n ) dS′n L mn ≈ ⎨ ⎩ wn Gs (rm, r′n ) dS′n + a [I−1 + (q 2 /2) I1]
(52)
⎪
for p = −3, −1, 1. From the Iup , and Ivp , it is enough to calculate one of them and the other can be obtained from interchanging the roles of η , and ξ in the final result when the limits are symmetrical. As mentioned above, it can be shown that when we have symmetric limits, the integrals of odd powers of η and ξ are zero. This limits the Gaussian-Legendre quadrature to odd orders, since it is not possible to define symmetric limits for even orders. It can be demonstrated that lim ζ →0+ζIp = 0 for p = −1, 1 except for the integral with p = −3 where lim ζ →0+ζI−3 = 2π / tuj tvj . Calculating the iterated double integrals given in Eqs. (50)–(52), calling Δsu = tuj Δuj , Δsv = tvj Δvj , Ω = Δsu2 + Δsv2 , Δu = (Ω + Δsu )/(Ω − Δsu ), and Δv = (Ω + Δsv )/(Ω − Δsv ), we have ζ →0+
(62)
which can be expressed as a linear system of algebraic equations by using the LCN with Q/2 equations and Q unknown variables as
(50)
p
(ζ 2nj2 + tu2j η2 + tv2j ξ 2 ) 2 η2dηdξ,
2
2
Δvj 2 Δvj − 2
+
tu2j η2
Δuj
Δvj ζ →0
(ζ 2nj2
1 1 [2(6Δsv2 + Δsu2 )Δsu Δsv Ω − Δsu5 ln Δv + 4Δsv5 ln Δu ]. tuj tv3j 640
All the terms in these analytic integrals have a contribution that is dependent on the dimension of the patches and the magnitudes of the tangential vectors.
(48)
⎛ q2 q4 ⎞ l ·RHdA ≈ 0.5 + aτnu ⎜Iu−3 − γa (r, r′) n Iu − Iu ⎟ ⎝ 2 −1 8 1⎠ ⎛ q2 q4 ⎞ +aτnv ⎜Iv−3 − Iv − Iv ⎟ , ⎝ 2 −1 8 1⎠
Δvj 2 Δvj − 2
1 1 [2(6Δsu2 + Δsv2 )Δsu Δsv Ω − Δsv5 ln Δu + 4Δsu5 ln Δv ], tu3j tvj 640
(61)
With these results, our integrals that will constitute the matrix elements can be expressed approximately as
i
(59)
(60) (46)
⎛ζ q2 ζ q4 ⎞ ⎛ 1 q2 1 q4 ⎞ l ·RH = a ⎜ 3 − γa (r, r′) n − ζR⎟ +a ⎜ 3 − − R⎟ τnu η2 ⎝R 2 R 8 ⎠ ⎝R 2 R 8 ⎠ ⎛1 q2 1 q4 ⎞ +a ⎜ 3 − − R⎟ τnv ξ 2. ⎝R 2 R 8 ⎠ (47)
∫σ
⎞ Δs 3 1 1 ⎛ ⎜Δsv Δsu Ω + Δsv3 ln Δu − u ln Δv ⎟ , 3 2 tuj tvj 12 ⎝ ⎠
2π , tuj tvj
1 (Δsu ln Δv + Δsv ln Δu ), tuj tvj
1 1 (4Δsu Δsv Ω + Δsu3 ln Δv + Δsv3 ln Δu ), tuj tvj 24
Nmn
⎧ wn γ (rm, r′n ) n ln ·RmndS′n if ⎪ ⎪ wn γs (rm, r′n ) n ln·RmndS′n + 0.5 ≈⎨ ⎪+ a τnu [Iu−3 − (q 2 /2) Iu−1 − (q 4 /8) Iu1] if ⎪+ a τ [I − (q 2 /2) I − (q 4 /8) I ] nv v−3 v−1 v1 ⎩
if if
Rmn > ε Rmn ≤ ε
(64)
Rmn > ε Rmn ≤ ε (65)
where the dependence on rm , and rn is implicit in the Rmn vector. The parameter ε represents a quantity that determines the radius of a spherical vicinity with the singularity at its center. In our case we have found that ε ≈ D /25 in the case of a cubical cavity of side D or ra /25 in the case of a spherical cavity with radius ra gives good results. We shall determine the modes in a cubical and a spherical cavity independently by applying both boundary conditions, given the high difficulty found in setting up the proper matrix elements, L mn and Nmn, when we have to decide the number of terms necessary to obtain the most accurate results with the Taylor's expansions. For a perfectly soft surface the Dirichlet boundary condition ψ (r′) r′∈ S′ = 0 cancels the second integral term of Eq. (63). This allows us to analyze and set up directly just the L mn matrix elements involved in the homogeneous linear system of equations
(53)
(54)
(55)
Q /2
1 Iu−3 = 3 Δsv ln Δu , tuj tvj Iv−3 =
1 Δsu ln Δv , tuj tv3j
⎞ Δs 3 1 1 ⎛ ⎜Δsu Δsv Ω + Δsu3 ln Δv − v ln Δu ⎟ , Iu−1 = 3 2 tuj tvj 12 ⎝ ⎠
∑ L mn Φn = 0. n =1
(56)
(66)
It has been demonstrated elsewhere [16] that the modes of the cavity can be determined by finding the minimal of the determinant defined as
(57)
Δ(ω ) = ln ( det L mn ).
(67)
In the case of a perfectly hard surface, the Neumann boundary condition (Φ (r′) r′∈ S′ = 0) implies the suppression of the first integral
(58) 114
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
term of Eq. (63),
Table 1 Eigen frequencies of cubic and spherical cavities calculated analytically and numerically expressed in reduced units ω .
Q /2
∑ Nmn ψn = 0.
(68)
n =1
Cubical Cavity
The modes under such conditions will be given by the minimal of the determinant
ψ (r′) r′∈ S′ = 0 ωa ω Lmn
(69)
Δ(ω ) = ln ( det Nmn ).
The reasons to consider the cubical and spherical cavities is that in these cases the Helmholtz equation can be solved analytically, and each kind of surface (with or without curvature) involves different terms of the expansions facilitating the process of setting up and programing the algorithm. The frequencies of the resonant modes for the Helmholtz equation for a cubic cavity with edges of length D, and placed with one of its corners at the origin of a Cartesian system of coordinates, can be obtained from the analytical relation
k=
1 ω π = (l 2 + m 2 + n 2 ) 2 , cl D
ωa ω Nmn
ψ (x, 0, z ) = ψ (x, D, z ) = 0,
(72)
ψ (x, y, 0) = ψ (x, y, D ) = 0.
(73)
ωa ω Lmn
∂ψ ∂n′
∂ψ ∂n′
= (0, y, z )
= (x,0, z )
= (x, y,0)
∂ψ ∂n′
(D , y , z )
∂ψ ∂n′
(x , D , z )
∂ψ ∂n′
(x , y , D )
ωa ω Nmn
= 0, (74)
= 0, (75)
= 0,
1 1 2 (l + m 2 + n 2 ) 2 . 2
(76)
(77)
It is important to note that in both boundary conditions, the position of the modes in the reduced frequency scheme are independent of the size of the cavity. On the other hand, when we have a spherical cavity, the analytic relation to determine the modes is [15].
ω =
αli . 2π
1.732 1.704
1.871 1.888
0.000 0.027
0.500 0.497
0.707 0.691
0.866 0.861
1.000 0.994
1.118 1.103
0.500 0.505
0.715 0.717
0.917 0.918
1.000 1.003
1.112 1.114
1.230 1.229
0.331 0.331
0.532 0.532
0.715 0.718
0.718 –
0.899 0.899
0.945 0.944
Similar situation holds when we consider a cubical cavity and the L mn elements are involved. However LCN presents some instability when the Nmn matrix elements enter in the calculus since the frequencies of the modes in this cubical cavity do not match very well with those obtained analytically; this is evident from the values in Table 1. In this case the partition parameters that give the best result were Δu = D /3 and a quadrature of order 3, which results in a matrix of dimensions 486×486. It is important to mention that the results do not improve by making a finer partition or increasing the order of the quadrature. This problem may be related to the presence of edges in the cubical cavity. Although in our present formulation of LCN we are being very careful of avoiding points along the edges in view that similar procedures have given good results in the case of two-dimensional complex systems like photonic crystals [17]. It is worth mentioning that some reports have been published mentioning this problem of LCN in systems involving edges in its geometry [8]. It is interesting to note that we compared LCN with another recently proposed method where we use Regular Quadrilateral Partitions (RQP) [16]. In order to obtain similar results for the cubical cavity a matrix of 384×384 was necessary. This denotes a modest advantage of RQP over LCN for the case of the cubical cavity. However, in the case of a spherical cavity when using RQP, a partition resulting in a matrix of 384×384 was sufficient to calculate the modes with a good precision. This means that RQP demands less computer resources than LCN in this case. We show the graphics of the determinants obtained by both numerical methods in Fig. 4. It is evident that the curve obtained by LCN seems to be not clearly converged.
the parameters will be l , m, n = 0, 1, 2, …. Expressing the frequency in reduced units ω = ω 2πcl / D , using D as a normalization factor, the reduced frequency spectrum for a cubic cavity will be
ω =
1.658 1.676
Φ (r′) r′∈ S′ = 0
In contrast, with the Neumann boundary conditions
∂ψ ∂n′
1.500 1.506
Spherical Cavity ψ (r′) r′∈ S′ = 0
where k is the magnitude of the wave vector as before, and l, m, and n are positive integer numbers l , m, n = 1, 2, 3, …, when considering the Dirichlet boundary conditions (71)
1.225 1.237
Φ (r′) r′∈ S′ = 0
(70)
ψ (0, y, z ) = ψ (D, y, z ) = 0,
0.866 0.87
(78) 4.1. 3D phononic crystal band structure
In this case we have taken as a normalization parameter the radius of the sphere ra , and αli represent the i-th zeros of the spherical Bessel's functions of order l when the boundary conditions are those stated in Eq. (71). If the derivative of the scalar field is involved, the boundary conditions are those given in Eq. (75). In that case the modes are given in terms of the zeros of the first derivative of the spherical Bessel's functions [15]. In Table 1 we compare the position of modes for the cubical and spherical cavities, obtained analytically and numerically by using the LCN for L mn and Nmn respectively. These results are in good agreement in the case of the spherical cavity when the partition to define the curvilinear patches is Δu = ra /3 and the order of the Gauss-Legendre quadrature is 3. This results in a matrix of dimensions 1854×1854.
It has been demonstrated elsewhere [17] that is possible to calculate the band structure of two-dimensional photonic crystals and simulate the propagation of light within these complex systems with the integral method. So, the immediate question that arises when we are aware of this fact, is if it could be possible to determine the band structure of a 3DPC with the acoustic version of LCN proposed in present work? To answer this question, let us consider a simple cubic empty unit cell (See Fig. 5(a)). Since we have symmetry along the three independent directions of this Cartesian system, we have to apply the Bloch periodicity conditions to the fields on the boundaries of the unit cell by 115
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
To determine the band structure by using LCN we consider an integral equation with sets of terms that represent each face of the cube numbered from 1 to 6 for a system analogous to Eq. (63) given by p
∑
p (1) (1) L mn (1) Φn (1) −
n =1
∑
p (1) (1) Nmn (1) ψn (1) +
p (1) (1) L mn (3) Φn (3) −
(1) (1) Nmn (2) ψn (2)
(1) (1) Nmn (3) ψn (3) +
∑
(1) (1) L mn (4) Φn (4)
n =1 p
p
p
∑
(1) (1) (1) (1) Nmn (4) ψn (4) + ∑ L mn (5) Φn (5) −
n =1 p
n =1 p
(1) (1) ∑ Nmn (5) ψn (5) n =1
(1) (1) (1) (1) ∑ Lmn (6) Φn (6) − ∑ Nmn (6) ψn (6) = 0. n =1
+
p
∑ n =1
n =1
∑ n =1
p
∑
+
(1) (1) L mn (2) Φn −
n =1
n =1
−
p
∑
(79)
n =1
To distinguish the points corresponding to each face of the cube in the field functions and matrix elements we have added the dependence of n(j), j = 1, 2, …6 . We have also added a superscript (1) to the matrix elements and the fields to indicate that they are dependent on the properties of the medium filling the unit cell. We are assuming that the normals point outward on the surface everywhere. By applying the Bloch periodicity conditions
Fig. 4. Comparison of results of RQP and LCN in the case of a cubical cavity. It is evident that the peaks at ω = 0.866 and ω = 1.225 are not clearly defined in the case of LCN.
(1) ik Bx D ψn(1) ,Φn(1)(2) = −Φn(1)(1) eik Bx D , (2) = ψn (1) e
(80)
(1) ik By D ψn(1) ,Φn(1)(4) = −Φn(1)(3) eik By D , (4) = ψn (3) e
(81)
(1) ik Bz D ψn(1) ,Φn(1)(6) = −Φn(1)(5) eik Bz D , (6) = ψn (5) e
(82)
Eq. (79) transforms to p
p
∑
(1) (1) ik Bx D ) Φ (1) − (L mn (1) − L mn (2) e n (1)
n =1
pairs of parallel faces of this cube. The band structure of a cubical empty cell can be calculated analytically starting with the wave vector in a periodic medium
k=
n =1 p
∑ n =1 p
∑
ω = kB + G , v
n =1
⎛ ⎞ ik Bz D ⎟ Φ (1) − (1) (1) ⎜⎜L mn ⎟ n (5) (5) − L mn (6) e ⎝ ⎠
l, m, n
n =1 p
⎛
(1) ∑ ⎜⎜Nmn (5) n =1
⎝
⎞ (1) e ik Bz D ⎟ ψ (1) = 0. + Nmn ⎟ n (5) (6) ⎠
(83)
D kB + G . 2π
2π ^ l), (l i + m^j + n k D
⎛ ⎞ ik By D ⎟ ψ (1) + (1) (1) ⎜⎜Nmn ⎟ n (3) (3) + Nmn (4) e ⎝ ⎠
p
∑
In this way the band structure shall be determined by seeking the minimal of the determinant of the system in very much the same as we did for the simple cubical cavity. It is found that the band structure calculated analytically and by LCN present almost good agreement since some of the modes does not coincide perfectly as it is evident in Fig. 5(b). Trying to go further with LCN, let us consider now a simple cubic 3DPC with a spherical inclusion of radius ra centered at the origin inside the cube, as shown in Fig. 6(a). The system of linear equations for LCN in this case is
In this case D is the lattice parameter of the unit cell, kB is the Bloch vector that scans the straight segments between the points: R (π / D, π / D, π / D )–M (π / D, π / D, 0)–X (0, 0, 0)–M. Finally G represents the vectors of the reciprocal lattice
G=
⎛ ⎞ ik Bx D ⎟ ψ (1) + (1) (1) ⎜⎜Nmn ⎟ n (1) (1) + Nmn (2) e ⎝ ⎠
⎛ ⎞ ik By D ⎟ Φ (1) − (1) (1) ⎜⎜L mn ⎟ n (3) (3) − L mn (4) e ⎝ ⎠
which when expressed in terms of the reduced frequency takes the form
ω =
∑
integers.
(b) (a)
0.5
z
0
−0.5 0.5 0.5 0
0
y
−0.5
−0.5
x
Fig. 5. (a) Cubic empty cell. The points indicated by stars indicate the centroid of each patch, and those by dots are the Nyström points of the corresponding partition. (b) Band structure of a simple cubic 3DPC with empty cells. The curves corresponding to dots were calculated analytically, and those corresponding to triangles were determined by using LCN.
116
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
(a)
(b)
0.5
z
0
−0.5 0.5
0.5
0
0 −0.5
y
−0.5
x
Fig. 6. Simple cubic cell containing a centered sphere. The points indicated by stars and dots indicate the Nyström points of the partition of the sphere and the cube respectively. (b) Band structure of a 3DPC with a simple cubic lattice and a filling fraction f=0.1. The curves with triangles represent the modes determined by LCN and the those traced with dots indicate the modes calculated by PWE.
p
∑
⎛ ⎞ ik Bx D ⎟ ψ (1) (1) (1) ⎜⎜Nmn ⎟ n (1) (1) + Nmn (2) e ⎠ n =1 ⎝ p ⎛ ⎞ ik By D ⎟ ψ (1) (1) (1) − ∑ ⎜⎜Nmn ⎟ n (3) (3) + Nmn (4) e ⎠ n =1 ⎝ p ⎛ (1) − ∑ ⎜⎜Nmn (5) n =1 ⎝ p
(1) (1) ik Bx D ) Φ (1) − (L mn (1) − L mn (2) e n (1)
n =1
p ⎛ ⎞ ik By D ⎟ Φ (1) (1) (1) + ∑ ⎜⎜L mn ⎟ n (3) (3) − L mn (4) e ⎠ n =1 ⎝ p ⎛ ⎞ ik Bz D ⎟ Φ (1) (1) (1) + ∑ ⎜⎜L mn ⎟ n (5) (5) − L mn (6) e ⎠ n =1 ⎝ r ⎞ (1) eik Bz D ⎟ ψ (1) + ∑ L (1) Φ (1) − + Nmn ⎟ n (5) mn (0) n (0) (6) ⎠ n =1
ρ2 ρ1
r
r (1) ψ (1) = 0, ∑ Nmn (0) n (0) n =1
(84)
r
(2) (2) (2) (1) ∑ Lmn (0) Φn (0) + ∑ (δ mn − Nmn (0) ) ψn (0) = 0, n =1
are not well defined could be the possible cause of this fail. When analyzing the band structure of a 3DPC with a simple cubic empty lattice LCN provides approximately the proper results. Trying to go beyond with LCN, if we consider a 3DPC with a simple cubic lattice of spheres, the results are catastrophic. Since the method gives completely wrong results in the calculation of the band structure, and even introduces some spurious modes. We can conclude that despite there are a number of papers where the potentialities of the method are mentioned when analyzing the scattering of electromagnetic or acoustic waves by a Mie sphere, LCN has some limitations when considering 3D geometries.
∑
n =1
Acknowledgments
(85)
A. Mendoza-Suárez and H. Pérez-Aguilar express their gratitude to the Coordinación de la Investigación Científica de la Universidad Michoacana de San Nicolás de Hidalgo for the partial financial support granted for the development of this research project. We are also indebted to Bernardo Mendoza and the authorities of Centro de Investigaciones en Óptica for the partial support to this work by allowing us the access to the Medusa cluster.
where we have applied the boundary conditions given in Eq. (7) to the surface of the sphere indicated by a n (0) and the medium filling the sphere by (2) in the corresponding functions. The presence of the Kronecker's delta in the last equation is due to the fact that the normal points inward to the surface of the sphere; so, θ (rm) = 1 in Eq. (1) and θ (rm) ψm(2)(0) (rm) = δmn ψn(2) (0) . To examine the performance of LCN in this case that we do not have an analytical expression to determine the band structure, let us compare the results of LCN with those obtained by the well known Plane Wave Expansion (PWE) method for a filling fraction f = D 3 /[(4/3) πra3] = 0.35 of water balloons (ρ1 = 998 kg/m3, v1 = 1496.7 m/s) in mercury (ρ2 = 13, 500 kg/m3, v2 = 1450 m/s) [18]. In this case LCN clearly fails to give satisfactory results, as we can appreciate from the band structures given in Fig. 6 (b). Thus, we can conclude that LCN itself presents some intrinsic limitations with 3D structures.
References [1] Gibson W. The method of moments in electromagnetics. New York: CRC Press; 2015. [2] Hanninen I, Taskinen M, Sarvas J. Singularity subtraction integral formulae for surface integral equations with RWG rooftop and hybrid basis functions. Prog Electromagn Res 2006;63:243–78. [3] Rao SM, Raju PK, Sun SP. Application of the method of moments to acoustic scattering from fluid-filled bodies of arbitrary shape. Num Methods Biomed Eng 1992;8:117–28. [4] Chandrasekhar B, Rao SM. Acoustic scattering from complex shaped three dimensional structures. Comp Model Eng Sci 2005;8:105–17. [5] Canino LF, Ottusch JJ, Stalzer MA, Visher JL, Wandzura SM. Numerical solution of the Helmholtz equation in 2D and 3D using a high-order Nyström discretization. J Comp Phys 1998;146:627–63. [6] Gedney S. On deriving a locally corrected Nyström Scheme from a quadrature sampled moment method. IEEE Trans Ant Propag 2003;51:2402–12. [7] Kolm P, Rokhlin V. Numerical quadratures for singular and hypersingular integrals. Comp Math Appl 2001;41:327–52. [8] Peterson AF. Mapped vector basis functions for electromagnetic integral equations. USA: Morgan & Calpool; 2005. [9] Burton AJ, Miller GF. The application of integral equation methods to the numerical solution of some exterior boundary-value problems. Proc Soc Lond 1971;323:201–10. [10] Piscoya R, Ochmann M. Acoustical boundary elements: theory and virtual experiments. Arch Acoust 2014;39:453–65. [11] Zaman SI. A comprehensive review of boundary integral formulations of acoustic scattering problems. Sci Tech Spec Rev 2000:281–310.
5. Conclusions We have formulated a Locally Corrected Nyström method based on a curvilinear orthogonal system of coordinates that involves the use of the curvilinear basis vectors and their first derivatives for the study of acoustic systems for interior problems and band structure calculations in 3DPC. Although we have not treated explicitly the scattering of acoustic waves by objects of arbitrary geometry it is completely possible. The idea of analyzing resonant modes in closed cavities is just to set up and test the matrix elements under different conditions. It was found that LCN presents some trouble in finding the proper modes in the case of a cubical cavity and contrary to common sense, it gives excellent results in the case of a spherical cavity. We could think that the presence of edges where the basis vectors and their derivatives
117
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
[16] Guel-Tapia JA, Villa-Villa F, Mendoza-Suarez A, Perez-Aguilar H. Acoustic scattering of 3D complex systems having random rough surfaces by scalar integral equations. Arch Acoust 2016;41:461–72. [17] Mendoza-Suárez A, Villa-Villa F, Gaspar-Armenta JA. Numerical method based on the solution of integral equations for the calculation of the band structure and reflectance of one and two-dimensional photonic crystals. J Opt Soc Am B 2006;23:2249–56. [18] Kushwaha MS, Djafari-Rouhani B. Complete acoustic stop bands for cubic arrays of spherical liquid balloons. J Appl Phys 1996;80:3191–5.
[12] Tadeu AJB, Goginho L, Santos P. Performance of the BEM solution in 3D acoustic wave scattering. Adv Eng Soft 2001;32:629–39. [13] Tong MS, Chew WC. Super-hypersingularity treatment for solving 3D electric field integral equations. Microw Opt Tech Lett 2007;49:1383–8. [14] Tong MS, Chew WC. A novel approach for evaluating hypersingular and strongly singular surface integrals in electromagnetics. IEEE Trans Ant Propag 2010;58, [3953–3601]. [15] Arfken GB, Weber HJ, Harris FE. Mathematical methods for physicists. USA: Academic Press; 2013.
118
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
The locally corrected Nyström method applied to 3D scalar SIE in acoustic cavities using curvilinear coordinates
MARK
⁎
F. Villa-Villaa, , H. Pérez-Aguilarb, A. Mendoza-Suárezb a b
Centro de Investigaciones en Óptica, Loma del Bosque 115, Lomas del Campestre 37150, León, Guanajuato, Mexico Facultad de Ciencias Físico Matemáticas, Universidad Michoacana de San Nicolás de Hidalgo, Av. Francisco J. Mújica S/N 58030, Morelia, Mich., Mexico
A R T I C L E I N F O
A BS T RAC T
Keywords: Integral equations Helmholtz equation Acoustic scattering
In this work we make a comparison of the results obtained by the Locally Corrected Nyström method and a parametric method previously proposed by the authors to solve the scalar surface integral equation in order to determine the resonant modes in 3D acoustic cavities. Both methods are set up in basis of a parametric representation in curvilinear orthogonal coordinates. So, any geometry that involves quadrics or higher order surfaces can be considered. Different cavities whose modes can be analytically determined for the purpose of comparison of results are considered. It is found that for a spherical cavity it gives the proper results locating the modes as predicted analytically while for the cubical cavity, the Locally Corrected Nyström method gives some imprecision in the calculation of the modes despite using high order numerical quadrature, probably due to the presence of edges where derivatives of the basis vectors are not defined.
1. Introduction In the study of scattering of acoustic and electromagnetic waves by objects of complex geometry, the method of moments (MoM) is known to be a good resource to solve the surface integral equations (SIE) involved in these problems [1–4]. So, for many years people has devoted some effort to develop this method specially using basis functions defined on surface elements that involve triangular facets. Even though almost arbitrary and complex objects can be represented by triangle meshing, it is known that as the size of the objects increase, some error problems can arise probably due to the representation of curved surfaces by flat triangular surface elements. When a surface is discretized into triangles, this meshing process presents itself a complex problem since MoM requires to have a precise account of the triangles of the partition and how they are coupled in the basis functions. In this regard specialized commercial and free software has been implemented but the intervention of the user is always necessary to solve details on meshing and interfacing the MoM's algorithms to the output of this software. An additional problem with MoM is the handling of the complex and lengthy mathematical expressions, resulting from a sampling process that involves nested double surface integrals which must be solved numerically and/or analytically to implement the numerical algorithm [2]. In more recent years, the Locally Corrected Nyström method (LCN) has been proposed as an alternative to MoM [5,6]. In this method
⁎
neither basis functions nor sampling process are necessary, since the integrals over the surface elements are directly calculated by using Gauss-Legendre quadrature rules when the argument of the Green's function and its derivatives are not in the vicinity of their singularities. The algorithm of LCN is far more simple to manage and implement, despite the integrals with singular arguments deserve a special treatment. However, a similar situation holds in this respect with MoM or any other method. For some years the standard procedure to treat integrals with singular kernels in LCN has consisted in designing highly specialized quadrature rules [6,7]. This procedure has the limitation that it is not possible to find a quadrature of integrals whose kernels present high order singularities like those appearing in the dyadic Green's function of the electromagnetic case. To lower the order of the singularities in the integrals, a procedure resembling the integration by parts is applied taking advantage of its convolutional double argument [6]. This procedure introduces an external differential operator to the integrals, complicating the numerical calculations. Since most of the research has been focused on the use of plane triangles (LCN and MoM) and low order basis functions in electromagnetics and acoustics as well, now people are concentrating on developing methods based on the use of curvilinear patches and high order basis functions [8]. In the case of MoM this problem is cumbersome [8], however, LCN seems to be more suited to curvilinear formulation.
Corresponding author. E-mail address: [email protected] (F. Villa-Villa).
http://dx.doi.org/10.1016/j.enganabound.2017.04.003 Received 21 September 2016; Received in revised form 25 January 2017; Accepted 15 April 2017 0955-7997/ © 2017 Elsevier Ltd. All rights reserved.
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
In this work we propose the solution of the SIE in acoustics by our own version of LCN using a parametric representation of the arguments in the surface integrals in a curvilinear orthogonal system. This l , tu , t v , formulation involves in a natural way the orthonormal basis n composed by the normal and tangent vectors in the coordinate directions (u , v ) of a parametric surface. To solve the integrals when the argument is in the vicinity of the singularities of the Green's function and its normal derivative, we propose the use of a Taylor series expansion to treat the diagonal elements by giving an explicit procedure that allows us to obtain analytic matrix elements. So, we do not require any interpolation process or the derivation of highly specialized quadrature rules that depend on the order of the singularities. The Taylor series expansion, besides the basis vectors involves the first derivatives of the tangent vectors ∂t u/∂u , ∂t v/∂v , which are related to the local curvatures of the surface. The detailed outline of this work reads as follows: In Section 2 we give a brief description of the integral equation related to the Helmholtz equation and the boundary conditions in acoustics. In Section 3 the Locally corrected Nyström method applied to curvilinear parametric surfaces by using the singularity extraction to treat the singular integrals is outlined. In Section 4 to validate the method, the application of LCN to determine resonant and Bloch modes in different physical systems like cavities and three-dimensional phononic crystal (3DPC) is considered, and finally, we conclude our work in Section 5 with an overview on the analytical and numerical results of this work.
Fig. 1. System with a curvilinear geometry composed of two regions (interior and exterior). The observer's position is indicated by r , while r′ denotes the source position on the surface for the integral surface equation. This vector spans the surface under integration.
ψ (r′) = f (r′), r′ ∈ S′, (Dirichlet ) ∂ψ = g (r′), r′ ∈ S′, (Neumann ) ∂n′
2. Theory As a starting point let us consider the SIE for the potential function ψ (r′) of the displacement vector (or the pressure as a function of position) in an acoustic medium
⎡
∫S′ ⎢⎣G (r, r′) ∂ψ∂n(r′ ′)
− ψ (r′)
∂G (r, r′) ⎤ ′ ⎥ dS = ψ (r) θ (r) − ψ inc (r), ∂n′ ⎦
are given, where f (r′) and g (r′) are arbitrary regular functions that define the corresponding boundary condition. 3. The Nyström method
(1)
where
⎧1 if r ∈ V1 θ (r) = ⎨ ⎩ 0 if r ∉ V1,
Let S′ be a closed surface that limits the volume V1, with no incident exciting field ψ inc (r) = 0 as shown in Fig. 1. By assuming that the normal points outward to S′ and considering the convention that the surface integral is positive if the normal points outward and negative in the opposite case. It is customary to validate a numerical method by calculating the scattered field by simple objects and compare with analytical results basically for a sphere. Since analytical results are also available for resonant modes in cavities with simple geometry, we shall be interested in the fields and their characteristic modes inside S′ given the Dirichlet or Neumann boundary conditions with f (r′) = 0 or g (r′) = 0 to test our method. As we will calculate the integrals by applying a limit when the observer approaches from the exterior of S′ along the normal, this implies that θ (r) = 0 . If we divide S′ into small surface elements Si like quadrilateral curvilinear patches, the integral equation for this region will be
(2)
is the Heaviside step function, and
G (r, r′) =
e qR , 4πR
(3)
represents Green's function for the 3D space. We have called q=ik, to further simply expressions where it appears repeatedley. ψ inc (r) represents an incident field for this exterior problem, R = r − r′ being r the observer's position and r′ the integration vector that spans the surface S′ (see Fig. 1), and
∂G (r, r′) l ·∇′G (r, r′) =n ∂n′
(4)
l·R =γ (r, r′) n
(5)
l is the outward is the normal derivative of the Green's function, n normal to the surface at r′, and
N
∑∫ i =1
γ (r, r′) = a (R−3 − qR−2 ) e qR .
(6)
1 ∂ψ (1) 1 ∂ψ (2) = , ρ1 ∂n′ ρ2 ∂n′
Si
G (r, r′)
∂ψ (r′) ′ dS − ∂n′
N
∑∫ i =1
Si
ψ (r′)
∂G (r, r′) ′ dS = 0. ∂n′
(9)
Although not explicitly indicated all the functions in the integrand depend on the properties of the acoustic medium of density ρ1, and wave number k1 = ω / v1 where v1 is the velocity of propagation of the perturbation. In this equation the Green's function and its normal derivative are known, and the potential ψ (r′) together with its normal derivative ∂ψ (r′)/∂n′ constitute the functions to be determined. Let us now suppose that we want to calculate the integrals over each patch Sn by using Gauss-Legendre quadrature
The potential and its normal derivative satisfy the boundary conditions at an interface given by
ψ (1) = ψ (2),
(8)
(7)
where ρ1 and ρ2 represent the densities of the corresponding media. Sometimes the conditions [9–12].
111
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al. N
Pi
∑ ∑ wj G (r, r′ j ) i =1 j =1
∂ψ (r′ j ) ∂n′
N
dS′ j −
Pi
∑ ∑ wj ψ (r′ j )
∂G (r, r′ j ) ∂n′
i =1 j =1
dS′ j = 0, (10)
where j = 1, 2, …, Pi are the points considered for the quadrature along the i-th surface element. These Pi limits of the summations can be different. Then, we get an equation with Q = 2N (P1 + P2+⋯+PN ) terms and the same number of unknown quantities. This equation has to be evaluated at the same points on the surface for the observer's position r = rm . Here m = 1, 2, …, N (P1 + P2+⋯+PN ), in this way we have a linear system of Q/2 equations with Q unknowns N
Pi
N
Pi
∑ ∑ wj G (rm, r′ j ) Φj dS′ j − ∑ ∑ wj ψj i =1 j =1
i =1 j =1
∂G (rm, r′ j ) ∂n′
dS′ j = 0. (11)
To simplify notation the double summations can be substituted by a simple one as Q /2
Q /2
∑ wk G (rm, r′n ) Φn dS′n − ∑ wn ψn n =1
n =1
∂G (rm, r′n ) ′ dS n = 0. ∂n′
(12) Fig. 2. Schematic description of the i-th curvilinear patch Si. The symmetric sub-patches σj are centered at the j-th point determined by the Gauss-Legendre quadrature rule in question.
If r → r′ the involved integrals in Eq. (11) can not be calculated by Gaussian quadrature. As we will state below, to regularize the problem when having functions with singularities of different orders, we use the well known technique of singularity extraction to separate the functions in one smooth, and one analytically integrable part that despite of being singular can be treated by taking a limit to determine their Cauchy principal value. This procedure has demonstrated to give good results when the surface is discretized by using triangle elements [13,14]. In this work we will demonstrate that the use of higher order surface elements is possible by means of curvilinear quadrilateral patches.
In analogy with the smooth and analytical parts of G (r, r′ j ), the functions γs (r, r′ j ) and γa (r, r′) stand for the smooth and analytical parts of the γ (r, r′ j ) function given in Eq. (6). These functions shall be obtained explicitly below when considering the singularity extraction technique. At this point it worth remembering that r′ j is the position vector of point j at the i-th patch. We have added a subindex s to the smooth parts and a to denote the analytical parts of the functions, and in the analytical integrals we have taken out the source functions
3.1. Local correction of the Nyström method
ψj = ψ (r′) r′= r j ,
When the argument of the Green's function and its normal derivative approaches to zero, that is when r → r′ the GaussLegendre quadrature no longer applies. So, a process of regularization is required to calculate the surface integrals by using a technique called singularity extraction. This method consists on separating the integrands in two parts: one smooth part Gs (r, r′) that can be integrated numerically by Gaussian quadrature and another that we call Ga (r, r′) that can be integrated analytically. Thus, the integrals over a curvilinear patch shall be
∫S
G (r, r′)
i
∂ψ (r′) ′ dS = ∂n′
∫S
i
∂ψ (r′) Gs (r, r′) dS′ + ∂n′
∫S
i
Pi
Φj =
∑
j =1 Pi ⎡ =∑ Φj ⎢wj Gs (r, r′ j ) dS′ j + ⎣ j =1
∫σ
j
∫σ
Ga (r, r′) dS′ j
⎤ Ga (r, r′) dS′⎥ ⎦ (13)
∫S
i
ψ (r′)
∂G (r, r′) ′ dS = ∂n′
∫S
l ·RdS′ ψ (r′) γs (r, r′) n
+
∫S
(16)
3.2. Integrands in the vicinity of the singularity r → r′: analytic functions To treat the integrals with analytic integrands, let us assume that the surface S′ is a parametric curvilinear surface dependent on two parameters (u , v ), where u, and v represent two orthogonal directions of the curvilinear system that represents the surface in question. If S′ is divided into Q small quadrilateral surface elements σj with centroid at (uj , vj ), so dS′ = Hdudv as shown in Fig. 2. In the expression of the differential of area, H represents the Jacobian of the transformation. The integrals over each curvilinear sub-patch σj can be expressed as
∂ψ (r′) Ga (r, r′) dS′ ∂n′ Φj
j =1
, r′= r j
by assuming that they are small enough to guarantee that these functions are almost constant in the integration sub-domains or subpatches σj. These sub-elements of area are centered at the points r′ j as indicated in Fig. 2.
Pi
≈∑ wj Φj Gs (r, r′ j ) dS′ j +
∂ψ (r′) ∂n′
(15)
i
∫σ
l ·RdS′ ψ (r′) γa (r, r′) n
i
Ga (r, r′) dS′ = j
vj +
Δvj
uj +
Δuj
∫v − Δv2 ∫u − Δu2 j
j
j
2
Ga (r, r′) Hdudv,
j
2
(17)
Pi
lj ·Rj dS′ j ≈∑ wj ψj γs (r, r′ j ) n
∫σ
j =1 Pi
+
∑ j =1
ψj
∫σ
l ·RdS′ γa (r, r′) n
⎡ lj ·Rj dS′ j + =∑ ψj ⎢wj γs (r, r′ j ) n ⎣ j =1
∫σ
j
j
vj +
Δvj
uj +
Δuj
∫v − Δv2 ∫u − Δu2 j
j
j
2
j
l ·RHdudv. γa (r, r′) n
2
(18)
If we consider the change of variables η = u − uj , and ξ = v − vj in Eqs. (17)–(18), these transform to
j
Pi
l ·RdS′ = γa (r, r′) n
⎤ l ·RdS′⎥ γa (r, r′) n ⎦
Δvj
∫σ
(14) 112
Ga (r, r′) dS′ = j
Δuj
∫− Δ2v ∫− Δ2u j
2
j
2
Ga [r, r′(η + uj , ξ + vj )] Hdηdξ
(19)
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al. Δuj
Δvj
∫S
l ·RdS′ = ψ (r′) γa (r, r′) n
n
∫− Δ2v ∫− Δ2u j
j
2
r′ ≈ r′ j + t uj η + t vj ξ + …
l γa [r, r′(η + uj , ξ + vj )] n
To calculate the integrals by taking their Cauchy principal value, let us define the observer's position as
2
(η + uj , ξ + vj )·R (η + uj , ξ + vj ) Hdηdξ.
(20)
r = r′ j + ζ nj ,
If we just recall that the position vector r′ denotes a point on the surface S′, by expanding into a Taylor's series around the point r′ j = r′(uj , vj ) r′(η + uj , ξ + vj ) = r′ j + 1 ∂ 2r′ 2 ∂ξ 2
∂r′ ∂η
η+ η =0, ξ =0
ξ2 η =0, ξ =0
∂r′ ∂ξ
ξ+ η =0, ξ =0
1 ∂ 2r′ + 2 ∂η∂ξ
1 ∂ 2r′ 2 ∂η2
η2 +
R = r − r′ ≈ ζ nj − ηt uj − ξ t vj .
(21)
∂r′ ∂η
∂r′ t vj = ∂ξ
= t u u = uj = v = vj
η =0, ξ =0
= t u u = uj v = vj
η =0, ξ =0
∂r′ ∂u
∂r′ = ∂v
u = uj v = vj
u = uj v = vj
1
R ≈ (ζ 2nj2 + tu2j η2 + tv2j ξ 2 ) 2 .
t′vj =
∂ 2r′ ∂u2 ∂ 2r′ ∂v 2
u = uj v = vj
u = uj v = vj
=
=
∂ 2r′ ∂η2
(22)
, (23)
∂ 2r′ ∂ξ 2
(24)
t′uvj
t′vuj =
u = uj v = vj
∂ 2r′ ∂v∂u
u = uj v = vj
η =0, ξ =0
∂ 2r′ = ∂η∂ξ
=
(32)
t v = t vj + t′vj ξ + t′vuj η ,
(33)
l = tu × t v n
(34)
lj + (t uj × t uvj + t′uj × t vj) η + (t uj × t′vj + t uvj × t vj) ξ ≈n
(35)
+(t′uj × t uvj) η2 + (t uvj × t′vj ) ξ 2,
(36)
. the Jacobian can be approximated as
(25)
H ≈ nn [1 + 0.5(t uj ·t′uj + t uvj ·t vj) η + 0.5(t vj ·t′vj + t uvj ·t uj) ξ ].
These vectors are illustrated in Fig. 3 for a spherical surface. We will define for short hand notation, the vectors
∂ 2r′ = ∂u∂v
t u = t uj + t′uj η + t′uvj ξ,
and an approach for the normal vector can be obtained to a first order
, η =0, ξ =0
(31)
If we now consider the expansion of the tangent vectors around rn up to first order,
,
and the curvature vectors, respectively
t′uj =
∂ 2r′ ∂ξ∂η
l ·R) H ≈ nj [ζnj2 − τnu η2 − τnv ξ 2], (n (26)
(38)
that appears in one of the integrands. In the last equation we have lj ·t′uj , and τnv = n lj ·t′vj to further simplify expressions. We called τnu = n also have omitted terms containing odd powers of η and ξ, in virtue that when integrated over symmetric limits the integrals are null.
. η =0, ξ =0
(37)
With these results, we can also calculate the product
, η =0, ξ =0
(30)
lj constitute a right where the tangent vectors t uj , t vj , and the normal n hand set orthonormal basis. With these considerations the magnitude of R , shall be
ξη +…, η =0, ξ =0
we can identify from the differential geometry the tangent vectors to the surface S′ at r′ j given by
t uj =
(29)
then, after integrating; the principal value will be that of the limit when ζ → 0+ [17]. In this way taking the difference of the last two equations we obtain
η =0, ξ =0
1 ∂ 2r′ ηξ + 2 ∂ξ ∂η η =0, ξ =0
(28)
(27)
Then, we can express Eq. (21) to a first order approximation 3.3. Singularity extraction The Taylor's expansion of the Green's function
G (R ) = a
exp(qR ) 1 = R 4π
∞
∑ m =0
q m m −1 R , m!
(39)
can be separated in two parts: one smooth, and one analytic [2].
G (R ) = Gs (R ) + Ga (R ),
(40)
where the smooth part is given by ∞
Gs (R ) = a
∑ m =1 m ≠2
q m m −1 R , m!
(41)
and the analytic part by
⎛1 q2 ⎞ Ga (R ) = a ⎜ + R⎟ . ⎝R 2! ⎠
(42)
In a similar fashion,
γ (r, r′) = γs (r, r′) + γa (r, r′), ∞
γs (r, r′) = a
∑ m =3 m ≠4
Fig. 3. Basis and curvature vectors on a spherical parametric surface.
113
1 − m m m −3 q R , m!
(43)
(44)
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
⎡1 q2 1 q4 ⎤ γa (r, r′) = a ⎢ 3 − − R⎥. ⎣R 2 R 8 ⎦
Iv−1 =
(45)
Multiplying the analytic functions Ga (r, r′), and γa (r, r′) by the other factors in the integrands
⎛1 q2 ⎞ Ga (r, r′) H = a ⎜ + R⎟ , ⎝R 2! ⎠
Iu1 =
Iv1 =
Ga (r, r′) HdA ≈ a I−1 + a
∫σ
i
q2 2
I1
4. Resonant modes in cubical and spherical cavities To test the proposed method, we shall compare its results with those obtained analytically for two highly symmetric cavities under the Dirichlet and Neumann boundary conditions [15]. One condition is related to perfectly soft surfaces and the other to perfectly hard surfaces. As stated in the previous sections, the complete system of SIE over a closed surface S′ will be
(49)
∫S′ G (r, r′) ∂ψ∂n(r′ ′) dS′ − ∫S′ ψ (r′) ∂G ∂(rn,′ r′) dS′ = 0,
where we have defined, the elementary surface integrals, Δuj 2 Δuj − 2
Ip = lim+
∫
Iup = lim+
∫− Δ2v ∫− Δ2u
ζ →0
∫
Ivp = lim+ ζ →0
j
∫
+
j
p tv2j ξ 2 ) 2 dηdξ,
Δuj 2 Δuj − 2
∫
Q /2
lim ζI−3 =
I−1 =
I1 =
Q /2
∑ L mn Φn − ∑ Nmn ψn = 0. n =1
(51)
(63)
n =1
In this expression we have called
(ζ 2nj2
+
tu2j η2
+
p tv2j ξ 2 ) 2 ξ 2dηdξ,
⎧ wn G (rm, r′n ) dS′n L mn ≈ ⎨ ⎩ wn Gs (rm, r′n ) dS′n + a [I−1 + (q 2 /2) I1]
(52)
⎪
for p = −3, −1, 1. From the Iup , and Ivp , it is enough to calculate one of them and the other can be obtained from interchanging the roles of η , and ξ in the final result when the limits are symmetrical. As mentioned above, it can be shown that when we have symmetric limits, the integrals of odd powers of η and ξ are zero. This limits the Gaussian-Legendre quadrature to odd orders, since it is not possible to define symmetric limits for even orders. It can be demonstrated that lim ζ →0+ζIp = 0 for p = −1, 1 except for the integral with p = −3 where lim ζ →0+ζI−3 = 2π / tuj tvj . Calculating the iterated double integrals given in Eqs. (50)–(52), calling Δsu = tuj Δuj , Δsv = tvj Δvj , Ω = Δsu2 + Δsv2 , Δu = (Ω + Δsu )/(Ω − Δsu ), and Δv = (Ω + Δsv )/(Ω − Δsv ), we have ζ →0+
(62)
which can be expressed as a linear system of algebraic equations by using the LCN with Q/2 equations and Q unknown variables as
(50)
p
(ζ 2nj2 + tu2j η2 + tv2j ξ 2 ) 2 η2dηdξ,
2
2
Δvj 2 Δvj − 2
+
tu2j η2
Δuj
Δvj ζ →0
(ζ 2nj2
1 1 [2(6Δsv2 + Δsu2 )Δsu Δsv Ω − Δsu5 ln Δv + 4Δsv5 ln Δu ]. tuj tv3j 640
All the terms in these analytic integrals have a contribution that is dependent on the dimension of the patches and the magnitudes of the tangential vectors.
(48)
⎛ q2 q4 ⎞ l ·RHdA ≈ 0.5 + aτnu ⎜Iu−3 − γa (r, r′) n Iu − Iu ⎟ ⎝ 2 −1 8 1⎠ ⎛ q2 q4 ⎞ +aτnv ⎜Iv−3 − Iv − Iv ⎟ , ⎝ 2 −1 8 1⎠
Δvj 2 Δvj − 2
1 1 [2(6Δsu2 + Δsv2 )Δsu Δsv Ω − Δsv5 ln Δu + 4Δsu5 ln Δv ], tu3j tvj 640
(61)
With these results, our integrals that will constitute the matrix elements can be expressed approximately as
i
(59)
(60) (46)
⎛ζ q2 ζ q4 ⎞ ⎛ 1 q2 1 q4 ⎞ l ·RH = a ⎜ 3 − γa (r, r′) n − ζR⎟ +a ⎜ 3 − − R⎟ τnu η2 ⎝R 2 R 8 ⎠ ⎝R 2 R 8 ⎠ ⎛1 q2 1 q4 ⎞ +a ⎜ 3 − − R⎟ τnv ξ 2. ⎝R 2 R 8 ⎠ (47)
∫σ
⎞ Δs 3 1 1 ⎛ ⎜Δsv Δsu Ω + Δsv3 ln Δu − u ln Δv ⎟ , 3 2 tuj tvj 12 ⎝ ⎠
2π , tuj tvj
1 (Δsu ln Δv + Δsv ln Δu ), tuj tvj
1 1 (4Δsu Δsv Ω + Δsu3 ln Δv + Δsv3 ln Δu ), tuj tvj 24
Nmn
⎧ wn γ (rm, r′n ) n ln ·RmndS′n if ⎪ ⎪ wn γs (rm, r′n ) n ln·RmndS′n + 0.5 ≈⎨ ⎪+ a τnu [Iu−3 − (q 2 /2) Iu−1 − (q 4 /8) Iu1] if ⎪+ a τ [I − (q 2 /2) I − (q 4 /8) I ] nv v−3 v−1 v1 ⎩
if if
Rmn > ε Rmn ≤ ε
(64)
Rmn > ε Rmn ≤ ε (65)
where the dependence on rm , and rn is implicit in the Rmn vector. The parameter ε represents a quantity that determines the radius of a spherical vicinity with the singularity at its center. In our case we have found that ε ≈ D /25 in the case of a cubical cavity of side D or ra /25 in the case of a spherical cavity with radius ra gives good results. We shall determine the modes in a cubical and a spherical cavity independently by applying both boundary conditions, given the high difficulty found in setting up the proper matrix elements, L mn and Nmn, when we have to decide the number of terms necessary to obtain the most accurate results with the Taylor's expansions. For a perfectly soft surface the Dirichlet boundary condition ψ (r′) r′∈ S′ = 0 cancels the second integral term of Eq. (63). This allows us to analyze and set up directly just the L mn matrix elements involved in the homogeneous linear system of equations
(53)
(54)
(55)
Q /2
1 Iu−3 = 3 Δsv ln Δu , tuj tvj Iv−3 =
1 Δsu ln Δv , tuj tv3j
⎞ Δs 3 1 1 ⎛ ⎜Δsu Δsv Ω + Δsu3 ln Δv − v ln Δu ⎟ , Iu−1 = 3 2 tuj tvj 12 ⎝ ⎠
∑ L mn Φn = 0. n =1
(56)
(66)
It has been demonstrated elsewhere [16] that the modes of the cavity can be determined by finding the minimal of the determinant defined as
(57)
Δ(ω ) = ln ( det L mn ).
(67)
In the case of a perfectly hard surface, the Neumann boundary condition (Φ (r′) r′∈ S′ = 0) implies the suppression of the first integral
(58) 114
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
term of Eq. (63),
Table 1 Eigen frequencies of cubic and spherical cavities calculated analytically and numerically expressed in reduced units ω .
Q /2
∑ Nmn ψn = 0.
(68)
n =1
Cubical Cavity
The modes under such conditions will be given by the minimal of the determinant
ψ (r′) r′∈ S′ = 0 ωa ω Lmn
(69)
Δ(ω ) = ln ( det Nmn ).
The reasons to consider the cubical and spherical cavities is that in these cases the Helmholtz equation can be solved analytically, and each kind of surface (with or without curvature) involves different terms of the expansions facilitating the process of setting up and programing the algorithm. The frequencies of the resonant modes for the Helmholtz equation for a cubic cavity with edges of length D, and placed with one of its corners at the origin of a Cartesian system of coordinates, can be obtained from the analytical relation
k=
1 ω π = (l 2 + m 2 + n 2 ) 2 , cl D
ωa ω Nmn
ψ (x, 0, z ) = ψ (x, D, z ) = 0,
(72)
ψ (x, y, 0) = ψ (x, y, D ) = 0.
(73)
ωa ω Lmn
∂ψ ∂n′
∂ψ ∂n′
= (0, y, z )
= (x,0, z )
= (x, y,0)
∂ψ ∂n′
(D , y , z )
∂ψ ∂n′
(x , D , z )
∂ψ ∂n′
(x , y , D )
ωa ω Nmn
= 0, (74)
= 0, (75)
= 0,
1 1 2 (l + m 2 + n 2 ) 2 . 2
(76)
(77)
It is important to note that in both boundary conditions, the position of the modes in the reduced frequency scheme are independent of the size of the cavity. On the other hand, when we have a spherical cavity, the analytic relation to determine the modes is [15].
ω =
αli . 2π
1.732 1.704
1.871 1.888
0.000 0.027
0.500 0.497
0.707 0.691
0.866 0.861
1.000 0.994
1.118 1.103
0.500 0.505
0.715 0.717
0.917 0.918
1.000 1.003
1.112 1.114
1.230 1.229
0.331 0.331
0.532 0.532
0.715 0.718
0.718 –
0.899 0.899
0.945 0.944
Similar situation holds when we consider a cubical cavity and the L mn elements are involved. However LCN presents some instability when the Nmn matrix elements enter in the calculus since the frequencies of the modes in this cubical cavity do not match very well with those obtained analytically; this is evident from the values in Table 1. In this case the partition parameters that give the best result were Δu = D /3 and a quadrature of order 3, which results in a matrix of dimensions 486×486. It is important to mention that the results do not improve by making a finer partition or increasing the order of the quadrature. This problem may be related to the presence of edges in the cubical cavity. Although in our present formulation of LCN we are being very careful of avoiding points along the edges in view that similar procedures have given good results in the case of two-dimensional complex systems like photonic crystals [17]. It is worth mentioning that some reports have been published mentioning this problem of LCN in systems involving edges in its geometry [8]. It is interesting to note that we compared LCN with another recently proposed method where we use Regular Quadrilateral Partitions (RQP) [16]. In order to obtain similar results for the cubical cavity a matrix of 384×384 was necessary. This denotes a modest advantage of RQP over LCN for the case of the cubical cavity. However, in the case of a spherical cavity when using RQP, a partition resulting in a matrix of 384×384 was sufficient to calculate the modes with a good precision. This means that RQP demands less computer resources than LCN in this case. We show the graphics of the determinants obtained by both numerical methods in Fig. 4. It is evident that the curve obtained by LCN seems to be not clearly converged.
the parameters will be l , m, n = 0, 1, 2, …. Expressing the frequency in reduced units ω = ω 2πcl / D , using D as a normalization factor, the reduced frequency spectrum for a cubic cavity will be
ω =
1.658 1.676
Φ (r′) r′∈ S′ = 0
In contrast, with the Neumann boundary conditions
∂ψ ∂n′
1.500 1.506
Spherical Cavity ψ (r′) r′∈ S′ = 0
where k is the magnitude of the wave vector as before, and l, m, and n are positive integer numbers l , m, n = 1, 2, 3, …, when considering the Dirichlet boundary conditions (71)
1.225 1.237
Φ (r′) r′∈ S′ = 0
(70)
ψ (0, y, z ) = ψ (D, y, z ) = 0,
0.866 0.87
(78) 4.1. 3D phononic crystal band structure
In this case we have taken as a normalization parameter the radius of the sphere ra , and αli represent the i-th zeros of the spherical Bessel's functions of order l when the boundary conditions are those stated in Eq. (71). If the derivative of the scalar field is involved, the boundary conditions are those given in Eq. (75). In that case the modes are given in terms of the zeros of the first derivative of the spherical Bessel's functions [15]. In Table 1 we compare the position of modes for the cubical and spherical cavities, obtained analytically and numerically by using the LCN for L mn and Nmn respectively. These results are in good agreement in the case of the spherical cavity when the partition to define the curvilinear patches is Δu = ra /3 and the order of the Gauss-Legendre quadrature is 3. This results in a matrix of dimensions 1854×1854.
It has been demonstrated elsewhere [17] that is possible to calculate the band structure of two-dimensional photonic crystals and simulate the propagation of light within these complex systems with the integral method. So, the immediate question that arises when we are aware of this fact, is if it could be possible to determine the band structure of a 3DPC with the acoustic version of LCN proposed in present work? To answer this question, let us consider a simple cubic empty unit cell (See Fig. 5(a)). Since we have symmetry along the three independent directions of this Cartesian system, we have to apply the Bloch periodicity conditions to the fields on the boundaries of the unit cell by 115
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
To determine the band structure by using LCN we consider an integral equation with sets of terms that represent each face of the cube numbered from 1 to 6 for a system analogous to Eq. (63) given by p
∑
p (1) (1) L mn (1) Φn (1) −
n =1
∑
p (1) (1) Nmn (1) ψn (1) +
p (1) (1) L mn (3) Φn (3) −
(1) (1) Nmn (2) ψn (2)
(1) (1) Nmn (3) ψn (3) +
∑
(1) (1) L mn (4) Φn (4)
n =1 p
p
p
∑
(1) (1) (1) (1) Nmn (4) ψn (4) + ∑ L mn (5) Φn (5) −
n =1 p
n =1 p
(1) (1) ∑ Nmn (5) ψn (5) n =1
(1) (1) (1) (1) ∑ Lmn (6) Φn (6) − ∑ Nmn (6) ψn (6) = 0. n =1
+
p
∑ n =1
n =1
∑ n =1
p
∑
+
(1) (1) L mn (2) Φn −
n =1
n =1
−
p
∑
(79)
n =1
To distinguish the points corresponding to each face of the cube in the field functions and matrix elements we have added the dependence of n(j), j = 1, 2, …6 . We have also added a superscript (1) to the matrix elements and the fields to indicate that they are dependent on the properties of the medium filling the unit cell. We are assuming that the normals point outward on the surface everywhere. By applying the Bloch periodicity conditions
Fig. 4. Comparison of results of RQP and LCN in the case of a cubical cavity. It is evident that the peaks at ω = 0.866 and ω = 1.225 are not clearly defined in the case of LCN.
(1) ik Bx D ψn(1) ,Φn(1)(2) = −Φn(1)(1) eik Bx D , (2) = ψn (1) e
(80)
(1) ik By D ψn(1) ,Φn(1)(4) = −Φn(1)(3) eik By D , (4) = ψn (3) e
(81)
(1) ik Bz D ψn(1) ,Φn(1)(6) = −Φn(1)(5) eik Bz D , (6) = ψn (5) e
(82)
Eq. (79) transforms to p
p
∑
(1) (1) ik Bx D ) Φ (1) − (L mn (1) − L mn (2) e n (1)
n =1
pairs of parallel faces of this cube. The band structure of a cubical empty cell can be calculated analytically starting with the wave vector in a periodic medium
k=
n =1 p
∑ n =1 p
∑
ω = kB + G , v
n =1
⎛ ⎞ ik Bz D ⎟ Φ (1) − (1) (1) ⎜⎜L mn ⎟ n (5) (5) − L mn (6) e ⎝ ⎠
l, m, n
n =1 p
⎛
(1) ∑ ⎜⎜Nmn (5) n =1
⎝
⎞ (1) e ik Bz D ⎟ ψ (1) = 0. + Nmn ⎟ n (5) (6) ⎠
(83)
D kB + G . 2π
2π ^ l), (l i + m^j + n k D
⎛ ⎞ ik By D ⎟ ψ (1) + (1) (1) ⎜⎜Nmn ⎟ n (3) (3) + Nmn (4) e ⎝ ⎠
p
∑
In this way the band structure shall be determined by seeking the minimal of the determinant of the system in very much the same as we did for the simple cubical cavity. It is found that the band structure calculated analytically and by LCN present almost good agreement since some of the modes does not coincide perfectly as it is evident in Fig. 5(b). Trying to go further with LCN, let us consider now a simple cubic 3DPC with a spherical inclusion of radius ra centered at the origin inside the cube, as shown in Fig. 6(a). The system of linear equations for LCN in this case is
In this case D is the lattice parameter of the unit cell, kB is the Bloch vector that scans the straight segments between the points: R (π / D, π / D, π / D )–M (π / D, π / D, 0)–X (0, 0, 0)–M. Finally G represents the vectors of the reciprocal lattice
G=
⎛ ⎞ ik Bx D ⎟ ψ (1) + (1) (1) ⎜⎜Nmn ⎟ n (1) (1) + Nmn (2) e ⎝ ⎠
⎛ ⎞ ik By D ⎟ Φ (1) − (1) (1) ⎜⎜L mn ⎟ n (3) (3) − L mn (4) e ⎝ ⎠
which when expressed in terms of the reduced frequency takes the form
ω =
∑
integers.
(b) (a)
0.5
z
0
−0.5 0.5 0.5 0
0
y
−0.5
−0.5
x
Fig. 5. (a) Cubic empty cell. The points indicated by stars indicate the centroid of each patch, and those by dots are the Nyström points of the corresponding partition. (b) Band structure of a simple cubic 3DPC with empty cells. The curves corresponding to dots were calculated analytically, and those corresponding to triangles were determined by using LCN.
116
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
(a)
(b)
0.5
z
0
−0.5 0.5
0.5
0
0 −0.5
y
−0.5
x
Fig. 6. Simple cubic cell containing a centered sphere. The points indicated by stars and dots indicate the Nyström points of the partition of the sphere and the cube respectively. (b) Band structure of a 3DPC with a simple cubic lattice and a filling fraction f=0.1. The curves with triangles represent the modes determined by LCN and the those traced with dots indicate the modes calculated by PWE.
p
∑
⎛ ⎞ ik Bx D ⎟ ψ (1) (1) (1) ⎜⎜Nmn ⎟ n (1) (1) + Nmn (2) e ⎠ n =1 ⎝ p ⎛ ⎞ ik By D ⎟ ψ (1) (1) (1) − ∑ ⎜⎜Nmn ⎟ n (3) (3) + Nmn (4) e ⎠ n =1 ⎝ p ⎛ (1) − ∑ ⎜⎜Nmn (5) n =1 ⎝ p
(1) (1) ik Bx D ) Φ (1) − (L mn (1) − L mn (2) e n (1)
n =1
p ⎛ ⎞ ik By D ⎟ Φ (1) (1) (1) + ∑ ⎜⎜L mn ⎟ n (3) (3) − L mn (4) e ⎠ n =1 ⎝ p ⎛ ⎞ ik Bz D ⎟ Φ (1) (1) (1) + ∑ ⎜⎜L mn ⎟ n (5) (5) − L mn (6) e ⎠ n =1 ⎝ r ⎞ (1) eik Bz D ⎟ ψ (1) + ∑ L (1) Φ (1) − + Nmn ⎟ n (5) mn (0) n (0) (6) ⎠ n =1
ρ2 ρ1
r
r (1) ψ (1) = 0, ∑ Nmn (0) n (0) n =1
(84)
r
(2) (2) (2) (1) ∑ Lmn (0) Φn (0) + ∑ (δ mn − Nmn (0) ) ψn (0) = 0, n =1
are not well defined could be the possible cause of this fail. When analyzing the band structure of a 3DPC with a simple cubic empty lattice LCN provides approximately the proper results. Trying to go beyond with LCN, if we consider a 3DPC with a simple cubic lattice of spheres, the results are catastrophic. Since the method gives completely wrong results in the calculation of the band structure, and even introduces some spurious modes. We can conclude that despite there are a number of papers where the potentialities of the method are mentioned when analyzing the scattering of electromagnetic or acoustic waves by a Mie sphere, LCN has some limitations when considering 3D geometries.
∑
n =1
Acknowledgments
(85)
A. Mendoza-Suárez and H. Pérez-Aguilar express their gratitude to the Coordinación de la Investigación Científica de la Universidad Michoacana de San Nicolás de Hidalgo for the partial financial support granted for the development of this research project. We are also indebted to Bernardo Mendoza and the authorities of Centro de Investigaciones en Óptica for the partial support to this work by allowing us the access to the Medusa cluster.
where we have applied the boundary conditions given in Eq. (7) to the surface of the sphere indicated by a n (0) and the medium filling the sphere by (2) in the corresponding functions. The presence of the Kronecker's delta in the last equation is due to the fact that the normal points inward to the surface of the sphere; so, θ (rm) = 1 in Eq. (1) and θ (rm) ψm(2)(0) (rm) = δmn ψn(2) (0) . To examine the performance of LCN in this case that we do not have an analytical expression to determine the band structure, let us compare the results of LCN with those obtained by the well known Plane Wave Expansion (PWE) method for a filling fraction f = D 3 /[(4/3) πra3] = 0.35 of water balloons (ρ1 = 998 kg/m3, v1 = 1496.7 m/s) in mercury (ρ2 = 13, 500 kg/m3, v2 = 1450 m/s) [18]. In this case LCN clearly fails to give satisfactory results, as we can appreciate from the band structures given in Fig. 6 (b). Thus, we can conclude that LCN itself presents some intrinsic limitations with 3D structures.
References [1] Gibson W. The method of moments in electromagnetics. New York: CRC Press; 2015. [2] Hanninen I, Taskinen M, Sarvas J. Singularity subtraction integral formulae for surface integral equations with RWG rooftop and hybrid basis functions. Prog Electromagn Res 2006;63:243–78. [3] Rao SM, Raju PK, Sun SP. Application of the method of moments to acoustic scattering from fluid-filled bodies of arbitrary shape. Num Methods Biomed Eng 1992;8:117–28. [4] Chandrasekhar B, Rao SM. Acoustic scattering from complex shaped three dimensional structures. Comp Model Eng Sci 2005;8:105–17. [5] Canino LF, Ottusch JJ, Stalzer MA, Visher JL, Wandzura SM. Numerical solution of the Helmholtz equation in 2D and 3D using a high-order Nyström discretization. J Comp Phys 1998;146:627–63. [6] Gedney S. On deriving a locally corrected Nyström Scheme from a quadrature sampled moment method. IEEE Trans Ant Propag 2003;51:2402–12. [7] Kolm P, Rokhlin V. Numerical quadratures for singular and hypersingular integrals. Comp Math Appl 2001;41:327–52. [8] Peterson AF. Mapped vector basis functions for electromagnetic integral equations. USA: Morgan & Calpool; 2005. [9] Burton AJ, Miller GF. The application of integral equation methods to the numerical solution of some exterior boundary-value problems. Proc Soc Lond 1971;323:201–10. [10] Piscoya R, Ochmann M. Acoustical boundary elements: theory and virtual experiments. Arch Acoust 2014;39:453–65. [11] Zaman SI. A comprehensive review of boundary integral formulations of acoustic scattering problems. Sci Tech Spec Rev 2000:281–310.
5. Conclusions We have formulated a Locally Corrected Nyström method based on a curvilinear orthogonal system of coordinates that involves the use of the curvilinear basis vectors and their first derivatives for the study of acoustic systems for interior problems and band structure calculations in 3DPC. Although we have not treated explicitly the scattering of acoustic waves by objects of arbitrary geometry it is completely possible. The idea of analyzing resonant modes in closed cavities is just to set up and test the matrix elements under different conditions. It was found that LCN presents some trouble in finding the proper modes in the case of a cubical cavity and contrary to common sense, it gives excellent results in the case of a spherical cavity. We could think that the presence of edges where the basis vectors and their derivatives
117
Engineering Analysis with Boundary Elements 79 (2017) 110–118
F. Villa-Villa et al.
[16] Guel-Tapia JA, Villa-Villa F, Mendoza-Suarez A, Perez-Aguilar H. Acoustic scattering of 3D complex systems having random rough surfaces by scalar integral equations. Arch Acoust 2016;41:461–72. [17] Mendoza-Suárez A, Villa-Villa F, Gaspar-Armenta JA. Numerical method based on the solution of integral equations for the calculation of the band structure and reflectance of one and two-dimensional photonic crystals. J Opt Soc Am B 2006;23:2249–56. [18] Kushwaha MS, Djafari-Rouhani B. Complete acoustic stop bands for cubic arrays of spherical liquid balloons. J Appl Phys 1996;80:3191–5.
[12] Tadeu AJB, Goginho L, Santos P. Performance of the BEM solution in 3D acoustic wave scattering. Adv Eng Soft 2001;32:629–39. [13] Tong MS, Chew WC. Super-hypersingularity treatment for solving 3D electric field integral equations. Microw Opt Tech Lett 2007;49:1383–8. [14] Tong MS, Chew WC. A novel approach for evaluating hypersingular and strongly singular surface integrals in electromagnetics. IEEE Trans Ant Propag 2010;58, [3953–3601]. [15] Arfken GB, Weber HJ, Harris FE. Mathematical methods for physicists. USA: Academic Press; 2013.
118