The interior resonance problem of time domain integral equations (TDIEs) formulated to analyze acoustic field interactions on penetrable objects is investigated. Two types of TDIEs are considered: The first equation, which is termed the time domain potential integral equation (TDPIE), suffers from the interior resonance problem, i.e., its solution is replete with spurious modes that are excited at the resonance frequencies of the acoustic cavity in the shape of the scatterer. Numerical experiments demonstrate that, unlike the frequency-domain integral equations, the amplitude of these modes in the time domain could be suppressed to a level that does not significantly affect the solution. This is achieved by increasing the numerical solution accuracy through the use of a higher-order discretization in space and the band limited approximate prolate spheroidal wave function with high interpolation accuracy as basis function in time. The second equation is obtained by linearly combining TDPIE with its normal derivative. The solution of this equation, which is termed the time domain combined potential integral equation (TDCPIE), does not involve any spurious interior resonance modes but it is not as accurate as the TDPIE solution at non-resonance frequencies. In addition, TDCPIE's discretization calls for treatment of hypersingular integrals.

Many applications in engineering and physical sciences call for simulations of acoustic scattering from penetrable objects, i.e., scatterers that internally support nonzero velocity potential and pressure field.1–11 An acoustic scattering problem involving penetrable objects is also known as an acoustic transmission problem. The time-harmonic (frequency-domain) acoustic transmission problem can be analyzed by solving a set of integral equations enforced on the surface of the scatterer.12 This set of equations is obtained by using the Kirchhoff-Helmholtz theorem to express the scattered fields of the exterior and interior problems in terms of (unknown) velocity potential on the surface of the scatterer and its normal derivative. The exterior and interior problems involve the unbounded domains with the material properties (density and wave speed) of the background medium and the scatterer, respectively. Numerical schemes developed to solve these integral equations discretize only the surface of the scatterer and implicitly enforce the radiation condition at infinity,12 offering advantages over finite element and finite difference methods that directly solve the Helmholtz equation, and require a volumetric discretization of the whole computation domain and use absorbing boundary conditions to approximate the radiation condition.

On the other hand, traditional integral equation formulations suffer from so-called “interior resonance” problem.13–21 This problem is observed when the excitation frequency approaches any one of the resonance frequencies of the acoustic cavity that is in the shape of the scatterer and has the density and the wave speed of the background medium. At these frequencies, the surface integral operator has a null space and the corresponding equation does not have a unique solution. Several approaches have been proposed to address the interior resonance problem of the frequency-domain integral equations. Examples of these include the combined Helmholtz integral equation formulation13 and the Burton-Miller scheme.14 

Even though interior resonance problem is well-studied for the frequency-domain integral equations, there are only a couple of studies that investigate the spurious interior resonance modes in the solution of time domain integral equations (TDIEs) of acoustics.22–25 In Ref. 22, a spurious resonance-free Burton-Miller-type time domain combined field integral equation is formulated to analyze acoustic scattering from sound-rigid bodies. In Ref. 26, interior resonance modes observed in the solution of the time domain electric field integral equation (of electromagnetics) that is enforced on perfect electrically conducting scatterers are investigated. Theoretically, TDIEs should not admit any interior resonance modes since their solution is obtained under zero initial condition and the interior resonance modes do not satisfy this initial condition.26,27 However, the interior resonance modes are still observed in the time domain solutions. It is discussed in Ref. 26 that this is because of the numerical errors introduced due to discretization and matrix inversions carried out during time marching.

In this work, the interior resonance problem of two different TDIEs formulated to analyze the transient acoustic transmission problem is investigated: The first equation, which is termed time domain potential integral equation (TDPIE) (in unknowns velocity potential and its normal derivative) here, is the time-domain equivalent of the frequency-domain integral equation that is traditionally used in the literature to solve the acoustic transmission problem.6,9,10,28 This equation suffers from the interior resonance problem. Its solution is replete with spurious modes that oscillate (without any decay) with the resonance frequencies of the acoustic cavity that is in the shape of the scatterer and has the density and the wave speed of the background medium.15 These modes are excited when the band of excitation includes their resonance frequency. In this work, it is demonstrated that unlike the frequency-domain integral equations, the amplitude of these modes in the time domain could be suppressed to a level, which does not significantly affect the solution, by increasing the accuracy of the numerical solution. This is achieved by using a higher-order Nyström discretization29–33 in space and the band limited approximate prolate spheroidal wave (APSW) function34–36 with high interpolation accuracy as basis function in time. The APSW function is a “two-sided” interpolator, which means that the resulting TDIE solver is non-causal. To address this and enable time-marching, a highly-accurate extrapolation scheme37,38 is used to “restore” the causality.

On the other hand, the second equation investigated in this work, which is termed time domain combined potential integral equation (TDCPIE), completely eliminates the interior resonance problem. The frequency-domain counterpart of TDCPIE has been introduced in Ref. 8 and it has been theoretically shown that its frequency-domain solution is unique at the interior resonance frequencies. Numerical results in this work verify that TDCPIE does not admit any interior resonance modes. It should be noted here that TDCPIE is obtained by linearly combining TDPIE with its normal derivative. Coupling parameters of this combination are carefully selected to enable the computation of the strongly-singular and hypersingular integrals that appear in the expressions of the matrix elements resulting from the Nyström discretization in space.29–33 In addition to this more-complicated singularity treatment method, the TDCPIE solution is, in general, less accurate than the TDPIE solution (except at the interior resonance frequencies of TDPIE).

Note that a preliminary version of this work has been described in a conference contribution.39 This paper includes a significantly more detailed formulation of TDPIE and TDCPIE and their discretization as well as a systematic investigation of the effects of solution accuracy on the spurious resonance modes associated with TDPIE. The remainder of the paper is organized as follows. In Sec. II, TDPIE and TDCPIE are derived. Section III describes the spatial and temporal discretization schemes and the marching-on-in-time method that is used to solve the resulting matrix system. In Sec. IV, numerical results are presented to validate the accuracy of TDPIE and TDCPIE solutions and demonstrate the relationship between the numerical errors and the interior resonance modes observed in the solution of the TDPIE. Section V concludes the paper with a short summary.

Let Ω2 and Ω1 denote the supports of an acoustically penetrable scatterer and the unbounded background medium in which this scatterer resides, respectively (Fig. 1). Let S represent the surface that separates these two domains, i.e., the surface of the scatterer. The wave speed and the density in Ωk,k{1,2}, are ck and ρk, respectively. An acoustic field with velocity potential φi(r,t) is incident on S. It is assumed that φi(r,t) is band limited to maximum frequency fmax and vanishingly small for t0 on rS. In response to this excitation, scattered fields with velocity potentials φks(r,t) are generated in Ωk. Total velocity potentials in Ω1 and Ω2 are expressed as φ1(r,t)=φi(r,t)+φ1s(r,t) and φ2(r,t)=φ2s(r,t), respectively.

FIG. 1.

(Color online) Description of the acoustic scattering problem.

FIG. 1.

(Color online) Description of the acoustic scattering problem.

Close modal

Using the Kirchhoff-Helmholtz theorem,40tφ1(r,t) and tφ2(r,t) are expressed as28 

tφ1(r,t)=tφi(r,t)tS1[nφ1](r,t)+tD1[φ1](r,t),rΩ1,
(1)
tφ2(r,t)=tS2[nφ2](r,t)tD2[φ2](r,t),rΩ2.
(2)

Here, t denotes the temporal derivative, n=n̂(r)·,n̂(r) is the outward pointing unit normal at point rS, and the spatiotemporal integral operators Sk[x](r,t) and Dk[x](r,t) are given by

Sk[x](r,t)=SGk(|rr|,t)*x(r,t)ds,Dk[x](r,t)=SnGk(|rr|,t)*x(r,t)ds,

where “*” denotes temporal convolution and

Gk(|rr|,t)=δ(t|rr|/ck)4π|rr|

is the time domain Green function. Note that since the Green function is in the form of a Dirac delta function δ(tt0), the temporal convolutions in operators Sk[x](r,t) and Dk[x](r,t) reduce to retarded time integrals.

On S, the acoustic pressure field and the normal component of the velocity field are continuous, i.e., the velocity potential satisfies the following boundary conditions:28 

ρ1tφ1(r,t)=ρ2tφ2(r,t),rS,
(3)
nφ1(r,t)=nφ2(r,t),rS.
(4)

Taking the limit of Eqs. (1) and (2) as r approaches S from Ωk,k{1,2} and inserting Eqs. (3) and (4) into the resulting equations yield TDPIE as28 

12tφ1(r,t)=tφi(r,t)tS1[nφ1](r,t)+tD̃1[φ1](r,t),rS,
(5)
ρ12ρ2tφ1(r,t)=tS2[nφ1](r,t)ρ1ρ2tD̃2[φ1](r,t),rS.
(6)

Here, “∼” on top of D̃k means that the space integral in Dk is evaluated in the principal value sense.41 

Note that the time-derivative form of TDIEs is used in this paper because this form of TDPIE is the one originally formulated in Ref. 28 to analyze the transient acoustic transmission problem. Ideally, one could remove the time derivative and still solve the resulting equation. However, as explained in Refs. 42–45, for the solution of a TDIE to be stable, basis and testing functions should be selected from appropriate Sobolev spaces. More specifically, the basis function should reside in the domain space of the integral operator and the testing function should reside in the space dual to the range of the integral operator. Following the analysis carried out in Refs. 42 and 43, the testing function for TDPIE without the time derivative should be selected as the time-derivative of the Dirac delta function. This means that a marching method, which relies on point-testing in time (i.e., testing with Dirac delta function) as described in Sec. III of this paper, does not yield a stable solution when it is used to solve TDPIE without the time derivative.

Taking the normal derivative of Eqs. (1) and (2) yields

tnφ1(r,t)=tnφi(r,t)tD1[nφ1](r,t)+tN1[φ1](r,t),rΩ1,
(7)
tnφ2(r,t)=tD2[nφ2](r,t)tN2[φ2](r,t),rΩ2,
(8)

where the spatiotemporal integral operators Dk[x](r,t) and Nk[x](r,t) are given by

Dk[x](r,t)=SnGk(|rr|,t)*x(r,t)ds,Nk[x](r,t)=Snn2Gk(|rr|,t)*x(r,t)ds,

and nn2=nn denotes the double normal derivative. Taking the limit of Eqs. (7) and (8) as r approaches S from Ωk,k{1,2} and inserting Eqs. (3) and (4) into the resulting equations yield the normal derivative of TDPIE in Eqs. (5) and (6) as

12tnφ1(r,t)=tnφi(r,t)tD̃1[nφ1](r,t)+tN1[φ1](r,t),rS,
(9)
12tnφ1(r,t)=tD̃2[nφ1](r,t)ρ1ρ2tN2[φ1](r,t),rS.
(10)

Linearly combining Eqs. (5) and (6) and (9) and (10) as α1 [Eq. (5)] +α2ρ2/ρ1 [Eq. (6)] and β1[Eq. (9)] +β2[Eq. (10)] yields TDCPIE as

α1+α22tφ1(r,t)α1tD̃1[φ1](r,t)+α2tD̃2[φ1](r,t)+α1tS1[nφ1](r,t)α2ρ2ρ1tS2[nφ1](r,t)=α1tφi(r,t),rS,
(11)
β1tN1[φ1](r,t)+β2ρ1ρ2tN2[φ1](r,t)+β1+β22tnφ1(r,t)+β1tD̃1[nφ1](r,t)β2tD̃2[nφ1](r,t)=β1tnφi(r,t),rS.
(12)

Here, αk and βk are real constants. They are selected as described in Sec. III C to cancel out the strongly-singular and hypersingular terms in Eqs. (11) and (12) and enable accurate computation of the matrix entries arising from their discretization.

To solve Eqs. (5) and (6) and (11) and (12) numerically, first, S is discretized into a mesh of curvilinear triangles and surface unknowns φ1(r,t) and nφ1(r,t) are expanded in space and time as

φ1(r,t)=i=1Ntq=1Npn=1NnIi1|qnϑ(r)qn(r)T(tiΔt),
(13)
nφ1(r,t)=i=1Ntq=1Npn=1NnIi2|qnϑ(r)qn(r)T(tiΔt).
(14)

In Eqs. (13) and (14), Nt is the number of time steps, Np is the number of curvilinear triangles, Nn is the number of interpolation nodes on each triangle, qn(r) is the Lagrange interpolation function defined at rqn (node n on triangle q),29,ϑ(r) is the inverse of the Jacobian of the coordinate transformation between the unit right triangle and the Cartesian coordinate system, T(t) is the temporal basis function, which is constructed using the APSW function,34–36 Δt is the time step size, and Ii1|qn and Ii2|qn are the unknown expansion coefficients to be solved for.

Substituting Eqs. (13) and (14) into Eqs. (5) and (6) and (11) and (12) and point-testing the resulting equations in space at rpm,p=1,,Np,m=1,,Nn (i.e., Nyström discretization in space), and in time at t=jΔt yield the following system of matrix equations:

Z0Ij=Vji=1j1ZjiIii=j+1j+NhwZjiIi,j=1,,Nt.
(15)

Here, Nhw is the half-width of the APSW function T(t), Ij=[Ij1Ij2]T, where Ij1 and Ij2 store the unknown expansion coefficients in Eqs. (13) and (14), Vj=[Vj1Vj2]T store the tested excitation vectors at time t=jΔt, and

Zji=[Zji11Zji12Zji21Zji22]

store the discretized retarded time integrals between nodes of mesh elements. Expressions of elements of Vj and Zji for TDPIE and TDCPIE are provided in Secs. III A and III B, respectively. Note that the system of matrix equations in Eq. (15) is not causal, i.e., Ij can not be solved for without knowing “future” unknowns Ij+1,Ij+2,…, Ij+Nhw [see the second summation on the right-hand side of Eq. (15)]. An extrapolation scheme34,37,38 can be used to represent these future unknowns Ij+1,Ij+2,…, Ij+Nhw in terms of “past/current” unknowns IjN+1,…, Ij1,Ij:

Ij+nm=1N0{pn}m+NIj+m,n=1,,Nhw.

Here, N is the number of samples used in the extrapolation and pn stores the extrapolation coefficients. Inserting this expression into Eq. (15) converts it into a causal form as

Z¯0Ij=Vji=1j1Z¯jiIi,j=1,,Nt,
(16)

where the modified matrices Z¯ji are expressed as

Z¯ji={Zji+n=1Nhw{pn}Nj+iZn,0ji<NZji,jiN.

The extrapolation coefficients pn are computed using the method described in Refs. 37 and 38. One can also use the scheme described in Ref. 34, but this comes with a slight decrease in accuracy.

The system of matrix equations in Eq. (16) is now in a form that can be recursively solved for the unknown coefficient vectors Ij,j=1,,Ntvia time marching as briefly described next. For j = 1, I1 is found by solving Eq. (16) with the right-hand side V1 (contribution from the summation is zero at the first time step). For j = 2, the right-hand side of Eq. (16) is computed by subtracting Z¯1I1 (only term coming from the summation) from V2. Then, I2 is found by solving Eq. (16) with this right-hand side. For j=3, the right-hand side of Eq. (16) is computed by subtracting Z¯2I1+Z¯1I2 from V3. Then, Eq. (16) is solved for I3. This recursive time marching algorithm is continued until all Ij,j=1,,Nt are obtained.

The marching-on-in-time method described above, thanks to the higher-order discretization in space and the use of the APSW function as basis function in time, yields a highly-accurate numerical solution. Indeed, results presented in Sec. IV show that, compared to the traditional marching-on-in-time methods, which uses polynomial temporal basis functions, it significantly suppresses the amplitude of the spurious non-decaying oscillations due to the interior resonance modes.

In the next two sections, Secs. III A and III B, the expressions of the elements of Vj and Zji in Eq. (15) are provided for TDPIE and TDCPIE, respectively.

The elements of Vj1 and Vj2 are given as Vj1|pm=tφi(rpm,t)|t=jΔt and Vj2|pm=0, respectively. The elements of Zji are expressed as

Zji11|pm,qn=12ϑ(rpm)tT(t)|t=(ji)ΔtδpqδmntSqnG1(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt,Zji12|pm,qn=tSqG1(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt,Zji21|pm,qn=12ϑ(rpm)tT(t)|t=(ji)Δtδpqδmn+tSqnG2(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt,Zji22|pm,qn=ρ2ρ1tSqG2(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt.
(17)

Here, R=|rpmr|,Sq is the surface of the curvilinear triangle q, and δpq=1 for p = q, and δpq=0 for pq.

The elements of Vj1 and Vj2 are given as Vj1|pm=α1tφi(rpm,t)|t=jΔt and Vj2|pm=β1tnφi(rpm,t)|t=jΔt, respectively. The elements of Zji are expressed as

Zji11|pm,qn=α1+α22ϑ(rpm)tT(t)|t=(ji)Δtδpqδmnα1tSqnG1(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt+α2tSqnG2(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt,Zji12|pm,qn=α1tSqG1(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔtα2ρ2ρ1tSqG2(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt,Zji21|pm,qn=β1tSqnn2G1(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt+β2ρ1ρ2tSqnn2G2(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt,Zji22|pm,qn=β1+β22ϑ(rpm)tT(t)|t=(ji)Δtδpqδmn+β1tSqnG1(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔtβ2tSqnG2(R,t)*T(tiΔt)ϑ(r)qn(r)ds|t=jΔt.
(18)

Computation of the matrix elements in Eqs. (17) and (18) calls for treatment of the singularities of the spatial integrals. This singularity treatment is described in Sec. III C.

There are four kinds of integrands in Eqs. (17) and (18): Gk(R,t)*T(tiΔt),nGk(R,t)*T(tiΔt),nGk(R,t)*T(tiΔt), and nn2Gk(R,t)*T(tiΔt). When p = q, and as r approaches rpm, the integrals of Gk(R,t)*T(tiΔt) and nGk(R,t)*T(tiΔt) become weakly-singular, the integral of nGk(R,t)*T(tiΔt) becomes strongly-singular, and the integral of nn2Gk(R,t)*T(tiΔt) becomes hypersingular.46 The weakly-singular integrals are computed using the Duffy transformation.47 The strongly-singular integrals in Eq. (17) are computed using the approach described in Ref. 30. The strongly-singular and hypersingular integrals in Eq. (18) are treated by extending the singularity subtraction and cancellation methods48–51 to account for the temporal discretization, more specifically the temporal basis function. This extension is briefly described next.

T(tR/ck) is expanded using the Taylor series around R=0 as

T(tR/ck)=u=0(1)utuT(t)u!ckuRu.
(19)

Using this expansion in n{T(tR/ck)/(4πR)} and nn2{T(tR/ck)/(4πR)} yields

nT(tR/ck)4πR=(n̂·R̂)4π[tT(tR/ck)ckR+T(tR/ck)R2]=(n̂·R̂)4π[u=0(1)utu+1T(t)u!cku+1Ru1+u=0(1)utuT(t)u!ckuRu2]=(n̂·R̂)4π[tT(t)ckR1+T(t)R2+(1)tT(t)ckR1+O(Ru0)]=(n̂·R̂)T(t)4πR2+O(Ru0),
(20)
nn2T(tR/ck)4πR=(n̂·n̂)4π[tT(tR/ck)ckR2+T(tR/ck)R3](n̂·R̂)(n̂·R̂)4π[t2T(tR/ck)ck2R+3tT(tR/ck)ckR2+3T(tR/ck)R3]=(n̂·n̂)4π[u=0(1)utu+1T(t)u!cku+1Ru2+u=0(1)utuT(t)u!ckuRu3](n̂·R̂)(n̂·R̂)4π[u=0(1)utu+2T(t)u!cku+2Ru1+3u=0(1)utu+1T(t)u!cku+1Ru2+3u=0(1)utuT(t)u!ckuRu3]=(n̂·n̂)4π[tT(t)ckR2+(1)t2T(t)ck2R1+T(t)R3+(1)tT(t)ckR2+t2T(t)2ck2R1+O(Ru0)](n̂·R̂)(n̂·R̂)4π[t2T(t)ck2R1+3tT(t)ckR2+3(1)t2T(t)ck2R1+3T(t)R3+3(1)tT(t)ckR2+3t2T(t)2ck2R1+O(Ru0)]=[(n̂·n̂)3(n̂·R̂)(n̂·R̂)]T(t)4πR3+[(n̂·n̂)+(n̂·R̂)(n̂·R̂)]t2T(t)8πck2R+O(Ru0).
(21)

Here, R̂=(rpmr)/R and O(Ru0) represents the higher-order terms that are not singular as R0. The singular terms on the right hands of Eqs. (20) and (21) are subtracted from nGk(R,t)*T(tiΔt) and nn2Gk(R,t)*T(tiΔt), respectively. Then, the integrals of these terms are added back to yield the final expressions for Zji11|pm,qn and Zji21|pm,qn in Eq. (18) as

Zji11|pm,qn=α1+α22ϑ(rpm)tT(t)|t=(ji)Δtδpqδmnα1tSq[nT(tR/c1)4πR(n̂·R̂)T(t)4πR2]ϑ(r)qn(r)ds|t=(ji)Δtα1Sq(n̂·R̂)tT(t)4πR2ϑ(r)qn(r)ds|t=(ji)Δt+α2tSq[nT(tR/c2)4πR(n̂·R̂)T(t)4πR2]ϑ(r)qn(r)ds|t=(ji)Δt+α2Sq(n̂·R̂)tT(t)4πR2ϑ(r)qn(r)ds|t=(ji)Δt,
(22)
Zji21|pm,qn=β1tSq{nn2T(tR/c1)4πR[(n̂·n̂)3(n̂·R̂)(n̂·R̂) ]T(t)4πR3[(n̂·n̂)+(n̂·R̂)(n̂·R̂) ]t2T(t)8πc12R }ϑ(r)qn(r)ds|t=(ji)Δtβ1Sq[(n̂·n̂)3(n̂·R̂)(n̂·R̂) ]tT(t)4πR3ϑ(r)qn(r)ds|t=(ji)Δtβ1Sq[(n̂·n̂)+(n̂·R̂)(n̂·R̂) ]t3T(t)8πc12Rϑ(r)qn(r)ds|t=(ji)Δt+β2ρ1ρ2tSq{nn2T(tR/c2)4πR[(n̂·n̂)3(n̂·R̂)(n̂·R̂) ]T(t)4πR3[(n̂·n̂)+(n̂·R̂)(n̂·R̂) ]t2T(t)8πc22R }ϑ(r)qn(r)ds|t=(ji)Δt+β2ρ1ρ2Sq[(n̂·n̂)3(n̂·R̂)(n̂·R̂) ]tT(t)4πR3ϑ(r)qn(r)ds|t=(ji)Δt+β2ρ1ρ2Sq[(n̂·n̂)+(n̂·R̂)(n̂·R̂) ]t3T(t)8πc22Rϑ(r)qn(r)ds|t=(ji)Δt.
(23)

The first and the third integrals on the right-hand side of Eq. (22) and the first and the fourth integrals on the right-hand side of Eq. (23) are the integrals of the integrands remaining after the singularity subtraction. Since these integrands are “smooth” (i.e., not singular as R0), their integrals can be computed accurately using a Gaussian quadrature rule.52 The third and the sixth integrals on the right-hand side of Eq. (23) are subtracted weakly-singular integrals and they are computed using the Duffy transformation.47 The second and the fourth integrals on the right-hand side of Eq. (22) are subtracted strongly-singular integrals and they cancel out each other for α1=α2. Similarly, the second and the fifth integrals on the right-hand side of Eq. (23) are subtracted hypersingular integrals and they cancel out each other for β1ρ2=β2ρ1.

In this section, numerical results, which demonstrate the relationship between numerical errors and interior resonance modes, are presented. TDPIE and TDCPIE are used to analyze acoustic scattering from a penetrable unit sphere that resides in an unbounded backgroud medium. It is assumed that the sphere is centered at the origin. The wave speed in the background medium and inside the sphere is c1=300 and c2=200m/s, respectively. The ratio of the densities in these two media is ρ1/ρ2=1.5. In all simulations, the excitation is a plane wave with velocity potential

φi(r,t)=φ0Gmod(tk̂i·r/c1),
(24)

where φ0 is the amplitude, k̂i is the unit vector along the direction of propagation, and Gmod(t)=cos[2πf0(ttp)]exp[(ttp)2/(2σ2)] is a modulated Gaussian pulse with center frequency f0, time delay tp, and duration σ. The excitation parameters are selected as φ0=1m2/s,k̂i=ẑ,f0=120Hz,tp=10σ, and σ=3/(2πfbw), where the effective bandwidth fbw=80Hz. Note that this definition of σ ensures that 99.997% of the wave energy is within the frequency band [fmin,fmax] with fmin=f0fbw and fmax=f0+fbw.53 Also, this specific selection of f0=120Hz and fbw=80Hz ensures that the frequency of the lowest interior resonance mode (the first cavity mode of the Dirichlet problem), 150Hz, is within the frequency band [fmin,fmax], i.e., this resonance mode could possibly be excited using the Gaussian pulse described above.15 The time step size is chosen as Δt=1/(2γfmax) with oversampling factor γ. Unless stated otherwise, the APSW function with an half-width of Nhw=7 is used to construct T(t). The surface of the sphere is discretized using Np=396 curvilinear triangles with Nn=6 interpolation nodes on each triangle. For TDCPIE, the linear combination coefficients in Eqs. (11) and (12) are chosen as α1=α2=β1=1 and β2=2/3. Lower-upper (LU) decomposition is used to solve the matrix system in Eq. (16) (at every time step) to ensure that the error in the matrix solution is at the machine precision level.54 

For the first set of simulations, the oversampling factor is selected as γ=6 resulting in Δt=0.42ms. Sixteen-point Gaussian and 9-point Gauss-Legendre quadrature rules52 are used to compute the two-dimensional (2D) surface integral with “smooth” integrand and the one-dimensional (1D) line integral needed for the Duffy transformation,47 respectively. First, three different TDPIE solutions are compared: The first solution is obtained using the proposed method where the temporal basis function T(t) is constructed using the APSW function. The second and the third solutions are obtained using the “traditional” marching-on-in-time method, where T(t) is constructed using third- and fourth-order Lagrange polynomial (LP) functions, respectively.55,56 It should be noted here that all solutions are obtained using the same discretization in space. Also, note that the traditional scheme with LP functions does not require an extrapolation scheme (as described in Sec. III) to ensure the causality of the time marching. Figure 2 plots the magnitude of the expansion coefficient Ij1|qn, q = 3, n = 3 (corresponding to rqn=(0.78,0.58,0.23)m), j=1,,Nt computed by solving TDPIE using the APSW, the third-order LP (termed “LP3”), and the fourth-order LP (termed “LP4”) functions. The figure clearly shows that all three solutions are corrupted by non-decaying oscillations that are especially visible in the late time. These oscillations correspond to the spurious interior resonance mode with frequency 150 Hz. The figure also shows that the amplitude of this resonance mode is significantly smaller in the TDPIE solution obtained using the APSW function. This is observed because the APSW function leads to a more accurate solution than the LP functions and this observation supports the claim that spurious resonance modes are induced in the TDPIE solution because of the numerical errors. It should be noted here that the solution obtained using the APSW function is more accurate because it is band limited (with a short duration time) and has excellent interpolation properties and continuous derivatives at t=jΔt while the LP functions have discontinuous derivatives and a wideband spectrum.34,35,37

FIG. 2.

(Color online) Magnitude of the expansion coefficient Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m], j=1,,Nt computed by solving TDPIE with the APSW, the LP3, and the LP4 functions and solving TDCPIE with the APSW function.

FIG. 2.

(Color online) Magnitude of the expansion coefficient Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m], j=1,,Nt computed by solving TDPIE with the APSW, the LP3, and the LP4 functions and solving TDCPIE with the APSW function.

Close modal

Next, the solutions of TDPIE and TDCPIE, both of which are obtained using the APSW function, are compared. Figure 2 plots the magnitude of the expansion coefficient Ij1|qn, q = 3, n = 3 (corresponding to rqn=(0.78,0.58,0.23)m), j=1,,Nt computed by solving these two equations. Clearly, the solution of TDPIE is corrupted by non-decaying oscillations while the solution of TDCPIE is free from any resonances. Figures 3(a), 3(b), and 3(c) compare the normalized Fourier transform of φ1(r,t) [i.e., Fourier transform of φ1(r,t) divided by the Fourier transform of Gmod(t)], which is computed after solving TDPIE and TDCPIE, in the time range t[0,1.67s],t[0,0.27s], and t[0,0.17s] to the frequency-domain (time-harmonic) total velocity potential computed at r=(0.78,0.58,0.23)m using the Mie series solution,57 respectively. Note that the normalization is required to ensure that the time-harmonic response (with equal excitation amplitude at each frequency) is obtained from the solutions of TDPIE and TDCPIE. Figure 3(a) clearly shows that both simulations are accurate within the effective band of the excitation, except at 150Hz where the solution of TDPIE is corrupted by the interior resonance mode. Figure 3(a) also shows that TDPIE is more accurate than TDCPIE at other frequencies. This might be explained by the fact that a second-kind surface integral equation (e.g., TDCPIE) is usually less accurate than its first-kind counterpart (e.g., TDPIE).58 The results in these three figures show that the visibility of the error due to the spurious resonance mode (in the solution of TDPIE) at 150 Hz decreases with the decreasing length of the Fourier transform window. More specifically, in Fig. 3(c), which presents the results with the shortest Fourier transform window, the error at 150 Hz cannot be distinguished from the numerical errors at other frequencies. This might be interpreted as that the amplitude of the spurious mode is already small enough not to be visible in the main scattering response. But when the Fourier transform window is increased to include the late time response, the presence of the spurious mode becomes more visible. Note that all three figures, as expected show that the solution of TDCPIE is free from any interior resonance modes.

FIG. 3.

(Color online)Normalized Fourier transform of φ1(r,t) computed after solving TDPIE and TDCPIE in the time range (a) t[0,1.67s], (b) t[0,0.27s], and (c) t[0,0.17s] and the frequency-domain total velocity potential computed at r=(0.78,0.58,0.23)m using the Mie series.

FIG. 3.

(Color online)Normalized Fourier transform of φ1(r,t) computed after solving TDPIE and TDCPIE in the time range (a) t[0,1.67s], (b) t[0,0.27s], and (c) t[0,0.17s] and the frequency-domain total velocity potential computed at r=(0.78,0.58,0.23)m using the Mie series.

Close modal

The effect of interior resonances on the scattered velocity potential is investigated by comparing the scattering cross section (SCS) of the sphere computed using the normalized Fourier transformed solutions of TDPIE and TDCPIE to SCS computed using the Mie series solution.57 Figures 4(a) and 4(b) plot SCS versus θ for ϕ=0° at 120 and 150 Hz, respectively. As shown in Fig. 4(a), SCS computed using TDPIE and TDCPIE solutions at 120 Hz shows good agreement with the Mie results. On the other hand, as shown in Fig. 4(b), the resonance mode at 150Hz dramatically changes SCS computed using the TDPIE solution.

FIG. 4.

(Color online)SCS computed using the Fourier transformed solutions of TDCPIE and TDPIE and the Mie series solution versus θ for ϕ=0° at (a) 120Hz and (b) 150Hz.

FIG. 4.

(Color online)SCS computed using the Fourier transformed solutions of TDCPIE and TDPIE and the Mie series solution versus θ for ϕ=0° at (a) 120Hz and (b) 150Hz.

Close modal

In this set of simulations, the effect of the computation accuracy of the integrals in Eqs. (17) and (18) on the amplitude of the interior resonance modes is investigated. Three sets of computation accuracy are considered by using different number of quadrature points to compute the 2D surface integral with “smooth” integrand and the 1D line integral needed for the Duffy transformation:47 

  • 16-point Gaussian and 9-point Gauss-Legendre quadrature rules

  • 7-point Gaussian and 5-point Gauss-Legendre quadrature rules

  • 4-point Gaussian and 3-point Gauss-Legendre quadrature rules

In all simulations, the oversampling factor is selected as γ=6 resulting in Δt=0.42ms. To clearly identify the interior resonance mode, Fourier transforms of the late-time data of Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m] computed after solving TDCPIE and TDPIE with integration accuracy sets (i), (ii), and (iii) in the time range t[0.8s,1s] are plotted in Fig. 5. The figure shows that the solutions of TDPIE with all three sets exhibit spurious interior resonance mode at 150Hz. However, as expected, no spurious resonance mode is present in the solutions of TDCPIE. Furthermore, the amplitude of the interior resonance mode observed in the solution of TDPIE with set (iii) is stronger than those with sets (i) and (ii). This is because set (iii) yields larger numerical errors than sets (i) and (ii).

FIG. 5.

(Color online) Fourier transforms of Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m] computed after solving TDCPIE and TDPIE with the integration accuracy sets (i), (ii), and (iii) in the time range t[0.8s,1s].

FIG. 5.

(Color online) Fourier transforms of Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m] computed after solving TDCPIE and TDPIE with the integration accuracy sets (i), (ii), and (iii) in the time range t[0.8s,1s].

Close modal

To compare the accuracy of the solutions of TDCPIE and TDPIE with sets (i), (ii), and (iii), the L2-norm error in the normalized Fourier transform of φ1(r,t) (over all 2376 interpolation nodes on the sphere surface) and the L2-norm error in SCS (for θ=[0°,180°] with 361 sample points and ϕ=0°) are plotted versus frequency in Figs. 6(a) and 6(b), respectively. Note that the reference data used in the computation of the L2-norm errors is obtained using the Mie series solution. Figure 6 shows that the accuracy of the TDPIE solution in the vicinity of 150Hz is significantly affected by the interior resonance mode. As expected, larger numerical errors increase the amplitude of the resonance mode. Figure 6 also shows that the spurious interior resonance mode is not observed in the TDCPIE solution regardless of the integral computation accuracy.

FIG. 6.

(Color online)L2-norm error (a) in the normalized Fourier transform of φ1(r,t) (at all the interpolation nodes on the sphere surface) and (b) in SCS for θ=[0°,180°] and ϕ=0° (at 361 sample points) computed after solving TDCPIE and TDPIE with the integration accuracy sets (i), (ii), and (iii).

FIG. 6.

(Color online)L2-norm error (a) in the normalized Fourier transform of φ1(r,t) (at all the interpolation nodes on the sphere surface) and (b) in SCS for θ=[0°,180°] and ϕ=0° (at 361 sample points) computed after solving TDCPIE and TDPIE with the integration accuracy sets (i), (ii), and (iii).

Close modal

To visualize the interior resonance mode at 150Hz, the normalized Fourier transform of φ1(r,t) (at all interpolation nodes on the sphere surface) computed after solving TDCPIE and TDPIE with set (iii) are presented in Figs. 7(a) and 7(b), respectively. Figure 7(c) presents the difference in the normalized Fourier transform of φ1(r,t) obtained from the TDCPIE and TDPIE solutions. As expected, the pattern of the difference follows the amplitude of the interior resonance mode.59 

FIG. 7.

(Color online)Patterns of the normalized Fourier transform of φ1(r,t) over the sphere surface at 150Hz. (a) TDCPIE solution. (b) TDPIE solution. (c) Difference between TDCPIE and TDPIE solutions.

FIG. 7.

(Color online)Patterns of the normalized Fourier transform of φ1(r,t) over the sphere surface at 150Hz. (a) TDCPIE solution. (b) TDPIE solution. (c) Difference between TDCPIE and TDPIE solutions.

Close modal

In this set of simulations, the effect of time step size Δt (i.e., temporal discretization density) on the amplitude of the interior resonance modes is demonstrated. Three different oversampling factors are used: γ{6,7.5,10} resulting in Δt={0.42,0.33,0.25}ms, respectively. In all simulations, the integration accuracy set (i) described in Sec. IV B is used to ensure that the numerical error resulting from the computation of the space integrals in matrix elements is suppressed to be sufficiently small. Figure 8 compares the magnitude of the expansion coefficient Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m], j=1,,Nt computed by solving TDPIE with γ=6,γ=7.5, and γ=10. The figure shows that the amplitude of interior resonance mode as identified by the non-decaying oscillations in the late time decreases by reducing the time step size. Figure 9 plots the Fourier transforms of Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m] computed after solving TDCPIE and TDPIE with γ=6,γ=7.5, and γ=10 in the time range t[0.8s,1s]. As expected, interior resonance mode is observed in the TDPIE solution at 150Hz. Furthermore, the figure shows that using a larger γ (or smaller Δt) reduces the amplitude of the resonance mode. As expected, no interior resonance mode is observed in the TDCPIE solution regardless of the value of γ used. Figures 10(a) and 10(b) plot the L2-norm error in the normalized Fourier transform of φ1(r,t) (over all 2376 interpolation nodes on the sphere surface) and the L2-norm error in SCS (for θ=[0°,180°] at 361 sample points and ϕ=0°) computed after solving TDCPIE and TDPIE with γ=6,γ=7.5, and γ=10, respectively. Note that the reference data used in the computation of the L2-norm errors is obtained using the Mie series solution. Figure 10 shows that interior resonance mode is observed in all of the TDPIE solutions with different γ but its amplitude could be significantly suppressed by increasing γ. Another point to note here is that, even though the TDCPIE solution does not admit any interior resonance modes, it is usually less accurate than the TDPIE solution within the whole effective band of the excitation except in the vicinity of 150Hz.

FIG. 8.

(Color online) Magnitude of the expansion coefficient Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m], j=1,,Nt computed by solving TDPIE with γ=6,γ=7.5, and γ=10.

FIG. 8.

(Color online) Magnitude of the expansion coefficient Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m], j=1,,Nt computed by solving TDPIE with γ=6,γ=7.5, and γ=10.

Close modal
FIG. 9.

(Color online) Fourier transforms of Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m] computed after solving TDCPIE and TDPIE with γ=6,γ=7.5, and γ=10 in the time range t[0.8s,1s].

FIG. 9.

(Color online) Fourier transforms of Ij1|qn, q = 3, n = 3 [corresponding to rqn=(0.78,0.58,0.23)m] computed after solving TDCPIE and TDPIE with γ=6,γ=7.5, and γ=10 in the time range t[0.8s,1s].

Close modal
FIG. 10.

(Color online) L2-norm error (a) in the normalized Fourier transform of φ1(r,t) (at all the interpolation nodes on the sphere surface) and (b) in SCS for θ=[0°,180°] and ϕ=0° (at 361 sample points) computed after solving TDCPIE and TDPIE with γ=6,γ=7.5, and γ=10.

FIG. 10.

(Color online) L2-norm error (a) in the normalized Fourier transform of φ1(r,t) (at all the interpolation nodes on the sphere surface) and (b) in SCS for θ=[0°,180°] and ϕ=0° (at 361 sample points) computed after solving TDCPIE and TDPIE with γ=6,γ=7.5, and γ=10.

Close modal

The interior resonance problem of TDPIE and TDCPIE that are formulated to analyze the time domain acoustic field interactions on penetrable scatterers is investigated. Numerical results demonstrate that the solution of TDPIE is corrupted by the spurious interior resonance modes that oscillate (without any decay) with the resonance frequencies of the acoustic cavity in the shape of the scatterer and has the density and the wave speed of the background medium. However, unlike the frequency-domain integral equations, the amplitude of these modes in the time domain can be suppressed by reducing the numerical error. This is achieved by using a higher-order Nyström discretization in space and the band limited APSW function with high interpolation accuracy as basis function in time.

The solution of TDCPIE, which is obtained by linearly combining TDPIE with its normal derivative, is free from spurious interior resonance modes but in general is less accurate than the solution of TDPIE (except at the interior resonance frequencies of TDPIE). In addition, the discretization of TDCPIE calls for treatment of strongly-singular and hypersingular integrals.

This publication is based upon work supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No. 2019-CRG8-4056. The authors would like to thank the King Abdullah University of Science and Technology Supercomputing Laboratory (KSL) for providing the required computational resources.

1.
C.
Yeh
, “
Scattering of acoustic waves by a penetrable prolate spheroid. I. Liquid prolate spheroid
,”
J. Acoust. Soc. Am.
42
(
2
),
518
521
(
1967
).
2.
D.
Brill
and
G. C.
Gaunaurd
, “
Acoustic resonance scattering by a penetrable cylinder
,”
J. Acoust. Soc. Am.
73
(
5
),
1448
1455
(
1983
).
3.
B.
Wang
,
K.
Zhao
,
J.
Fan
, and
G.
Zheng
, “
Modeling method on acoustic scattering from penetrable objects using a hybrid Kirchhoff/ray approach
,”
J. Acoust. Soc. Am.
141
(
5
),
3531
3531
(
2017
).
4.
G. T.
Schuster
, “
Solution of the acoustic transmission problem by a perturbed Born series
,”
J. Acoust. Soc. Am.
77
(
3
),
880
886
(
1985
).
5.
N. B.
Kakogiannos
and
J. A.
Roumeliotis
, “
Acoustic scattering from a sphere of small radius coated by a penetrable one
,”
J. Acoust. Soc. Am.
98
(
6
),
3508
3515
(
1995
).
6.
R.
Kittappa
and
R.
Kleinman
, “
Acoustic scattering by penetrable homogeneous objects
,”
J. Math. Phys.
16
(
2
),
421
432
(
1975
).
7.
D. C.
Thomas
,
K. L.
Gee
, and
R. S.
Turley
, “
A balloon lens: Acoustic scattering from a penetrable sphere
,”
Am. J. Phys.
77
(
3
),
197
203
(
2009
).
8.
R. E.
Kleinman
and
P. A.
Martin
, “
On single integral equations for the transmission problem of acoustics
,”
SIAM J. Appl. Math.
48
(
2
),
307
325
(
1988
).
9.
M.
Costabel
and
E.
Stephan
, “
A direct boundary integral equation method for transmission problems
,”
J. Math. Anal. Appl
106
(
2
),
367
413
(
1985
).
10.
R.
Kress
and
G.
Roach
, “
Transmission problems for the Helmholtz equation
,”
J. Math. Phys.
19
(
6
),
1433
1437
(
1978
).
11.
H.
Wu
,
Y.
Liu
, and
W.
Jiang
, “
A fast multipole boundary element method for 3D multi-domain acoustic scattering problems based on the Burton–Miller formulation
,”
Eng. Anal. Boundary Elem.
36
(
5
),
779
788
(
2012
).
12.
D.
Colton
and
R.
Kress
,
Integral Equation Methods in Scattering Theory
(
SIAM
,
New York
,
2013
). pp.
1
271
.
13.
H. A.
Schenck
, “
Improved integral formulation for acoustic radiation problems
,”
J. Acoust. Soc. Am.
44
(
1
),
41
58
(
1968
).
14.
A. J.
Burton
and
G. F.
Miller
, “
The application of integral equation methods to the numerical solution of some exterior boundary-value problems
,”
Proc. R. Soc. Lond. Ser. A
323
(
1553
),
201
210
(
1971
).
15.
C.-J.
Zheng
,
H.-B.
Chen
,
H.-F.
Gao
, and
L.
Du
, “
Is the Burton–Miller formulation really free of fictitious eigenfrequencies?
,”
Eng. Anal. Boundary Elem.
59
,
43
51
(
2015
).
16.
E.
Schulz
and
R.
Hiptmair
, “
Spurious resonances in coupled domain-boundary variational formulations of transmission problems in electromagnetism and acoustics
,” arXiv:2003.14357 (
2020
).
17.
A.
Buffa
and
R.
Hiptmair
, “
Regularized combined field integral equations
,”
Numer. Math.
100
(
1
),
1
19
(
2005
).
18.
Y.
Boubendir
,
V.
Dominguez
,
D.
Levadoux
, and
C.
Turc
, “
Regularized combined field integral equations for acoustic transmission problems
,”
SIAM J. Appl. Math.
75
(
3
),
929
952
(
2015
).
19.
D.
Jones
, “
Integral equations for the exterior acoustic problem
,”
Q. J. Mech. Appl. Math
27
(
1
),
129
142
(
1974
).
20.
J.-Y.
Hwang
and
S.-C.
Chang
, “
A retracted boundary integral equation for exterior acoustic problem with unique solution for all wave numbers
,”
J. Acoust. Soc. Am.
90
(
2
),
1167
1180
(
1991
).
21.
Z. Y.
Qian
,
Z. D.
Han
,
P.
Ufimtsev
, and
S. N.
Atluri
, “
Non-hyper-singular boundary integral equations for acoustic problems, implemented by the collocation-based boundary element method
,”
Comput. Model. Eng. Sci.
6
,
133
144
(
2004
).
22.
A.
Ergin
,
B.
Shanker
, and
E.
Michielssen
, “
Analysis of transient wave scattering from rigid bodies using a Burton–Miller approach
,”
J. Acoust. Soc. Am.
106
(
5
),
2396
2404
(
1999
).
23.
D. J.
Chappell
,
P. J.
Harris
,
D.
Henwood
, and
R.
Chakrabarti
, “
A stable boundary element method for modeling transient acoustic radiation
,”
J. Acoust. Soc. Am.
120
(
1
),
74
80
(
2006
).
24.
H.-W.
Jang
and
J.-G.
Ih
, “
Stabilization of time domain acoustic boundary element method for the exterior problem avoiding the nonuniqueness
,”
J. Acoust. Soc. Am.
133
(
3
),
1237
1244
(
2013
).
25.
J. H.
Kao
, “
A time-shifting algorithm for alleviating convergence difficulties at interior acoustic resonance frequencies
,”
Appl. Sci
11
(
6
),
2701.1–2701.17
(
2021
).
26.
Y.
Shi
,
H.
Bagci
, and
M.
Lu
, “
On the internal resonant modes in marching-on-in-time solution of the time domain electric field integral equation
,”
IEEE Trans. Antennas Propag.
61
(
8
),
4389
4392
(
2013
).
27.
B.
Shanker
,
A. A.
Ergin
,
K.
Aygun
, and
E.
Michielssen
, “
Analysis of transient electromagnetic scattering from closed surfaces using a combined field integral equation
,”
IEEE Trans. Antennas Propag.
48
(
7
),
1064
1074
(
2000
).
28.
J.
Li
,
D.
Dault
, and
B.
Shanker
, “
A quasianalytical time domain solution for scattering from a homogeneous sphere
,”
J. Acoust. Soc. Am.
135
(
4
),
1676
1685
(
2014
).
29.
G.
Kang
,
J.
Song
,
W. C.
Chew
,
K. C.
Donepudi
, and
J.-M.
Jin
, “
A novel grid-robust higher order vector basis function for the method of moments
,”
IEEE Trans. Antennas Propag.
49
(
6
),
908
915
(
2001
).
30.
R.
Chen
,
S. B.
Sayed
,
N.
Alharthi
,
D.
Keyes
, and
H.
Bagci
, “
An explicit marching-on-in-time scheme for solving the time domain Kirchhoff integral equation
,”
J. Acoust. Soc. Am.
146
(
3
),
2068
2079
(
2019
).
31.
R.
Chen
,
S. B.
Sayed
,
H. A.
Ulku
, and
H.
Bagci
, “
An explicit time marching scheme for efficient solution of the magnetic field integral equation at low frequencies
,”
IEEE Trans. Antennas Propag.
69
(
2
),
1213
1218
(
2021
).
32.
R.
Chen
and
H.
Bagci
, “
On higher-order Nyström discretization of scalar potential integral equation for penetrable scatterers
,” in
Proceedings of the Applied Computational Electromagnetics Society Symposium
, Miami, FL (April 14–18,
2019
).
33.
R.
Chen
and
H.
Bagci
, “
Explicit solution of time domain scalar potential surface integral equations for penetrable scatterers
,” in
Proceedings of the IEEE International Symposium on Antennas and Propagation
, Montreal, Canada (July 5–10,
2020
), pp.
1001
1002
.
34.
D. S.
Weile
,
G.
Pisharody
,
N.-W.
Chen
,
B.
Shanker
, and
E.
Michielssen
, “
A novel scheme for the solution of the time-domain integral equations of electromagnetics
,”
IEEE Trans. Antennas Propag.
52
(
1
),
283
295
(
2004
).
35.
R. A.
Wildman
,
G.
Pisharody
,
D. S.
Weile
,
S.
Balasubramaniam
, and
E.
Michielssen
, “
An accurate scheme for the solution of the time-domain integral equations of electromagnetics using higher order vector bases and bandlimited extrapolation
,”
IEEE Trans. Antennas Propag.
52
(
11
),
2973
2984
(
2004
).
36.
J.
Knab
, “
Interpolation of band-limited functions using the approximate prolate series (Corresp.)
,”
IEEE Trans. Inform. Theory
25
(
6
),
717
720
(
1979
).
37.
S. B.
Sayed
,
H. A.
Ulku
, and
H.
Bagci
, “
A stable marching on-in-time scheme for solving the time-domain electric field volume integral equation on high-contrast scatterers
,”
IEEE Trans. Antennas Propag.
63
(
7
),
3098
3110
(
2015
).
38.
A.
Glaser
and
V.
Rokhlin
, “
A new class of highly accurate solvers for ordinary differential equations
,”
J. Sci. Comput.
38
(
3
),
368
399
(
2009
).
39.
R.
Chen
and
H.
Bagci
, “
On the internal resonance modes of time domain surface integral equations for acoustic transmission problems
,” in
Proceedings of the URSI General Assembly and Scientific Symposium
, Rome, Italy (August 28–September 4,
2021
), pp.
1
4
.
40.
A. D.
Pierce
,
Acoustics: An Introduction to Its Physical Principles and Applications
(
Springer
,
New York
,
2019
), pp.
1
768
.
41.
A.
Ishimaru
, “
Electromagnetic wave propagation
,” in
Radiation, and Scattering: From Fundamentals to Applications
(
Wiley-IEEE
,
New York
,
2017
), pp.
401
442
.
42.
T. E.
Roth
and
W. C.
Chew
, “
Stability analysis and discretization of A-ϕ time domain integral equations for multiscale electromagnetics
,”
J. Comput. Phys.
408
,
109102
(
2020
).
43.
T. E.
Roth
and
W. C.
Chew
, “
Lorenz gauge potential-based time domain integral equations for analyzing subwavelength penetrable regions
,”
IEEE J. Multiscale Multiphys. Comput. Technol.
6
,
24
34
(
2021
).
44.
T.
Ha-Duong
, “
On retarded potential boundary integral equations and their discretisation
,” in
Topics in Computational Wave Propagation
(
Springer
,
Berlin, Germany
,
2003
), pp.
301
336
.
45.
T.
Ha-Duong
,
B.
Ludwig
, and
I.
Terrasse
, “
A Galerkin BEM for transient acoustic scattering by an absorbing obstacle
,”
Int. J. Numer. Methods Eng.
57
(
13
),
1845
1882
(
2003
).
46.
Y.
Liu
and
F. J.
Rizzo
, “
A weakly singular form of the hypersingular boundary integral equation applied to 3-D acoustic wave problems
,”
Comput. Methods Appl. Mech. Eng.
96
(
2
),
271
287
(
1992
).
47.
M. G.
Duffy
, “
Quadrature over a pyramid or cube of integrands with a singularity at a vertex
,”
SIAM J. Numer. Anal.
19
(
6
),
1260
1262
(
1982
).
48.
S.
Jarvenpaa
,
M.
Taskinen
, and
P.
Yla-Oijala
, “
Singularity subtraction technique for high-order polynomial vector basis functions on planar triangles
,”
IEEE Trans. Antennas Propag.
54
(
1
),
42
49
(
2006
).
49.
J.
Lin
,
W.
Chen
, and
C. S.
Chen
, “
Numerical treatment of acoustic problems with boundary singularities by the singular boundary method
,”
J. Sound Vib.
333
(
14
),
3177
3188
(
2014
).
50.
S.
Li
and
Q.
Huang
, “
An improved form of the hypersingular boundary integral equation for exterior acoustic problems
,”
Eng. Anal. Boundary Element.
34
(
3
),
189
195
(
2010
).
51.
J.
Li
,
X.
Fu
, and
B.
Shanker
, “
Decoupled potential integral equations for electromagnetic scattering from dielectric objects
,”
IEEE Trans. Antennas Propag.
67
(
3
),
1729
1739
(
2019
).
52.
J.-M.
Jin
,
Theory and Computation of Electromagnetic Fields
(
Wiley
,
Hoboken, NJ
,
2010
), pp.
1
572
.
53.
H.
Bagci
,
A. E.
Yilmaz
,
J.-M.
Jin
, and
E.
Michielssen
, “
Fast and rigorous analysis of EMC/EMI phenomena on electrically large and complex cable-loaded structures
,”
IEEE Trans. Electromagn. Compat.
49
(
2
),
361
381
(
2007
).
54.
E.
Anderson
,
Z.
Bai
,
C.
Bischof
,
L. S.
Blackford
,
J.
Demmel
,
J.
Dongarra
,
J.
Du Croz
,
A.
Greenbaum
,
S.
Hammarling
,
A.
McKenney
, and
D.
Sorensen
,
LAPACK Users' Guide
(
SIAM
,
Warrendale, PA
,
1999
), pp.
1
407
.
55.
K.
Aygun
,
B.
Shanker
,
A. A.
Ergin
, and
E.
Michielssen
, “
A two-level plane wave time-domain algorithm for fast analysis of EMC/EMI problems
,”
IEEE Trans. Electromagn. Compat.
44
(
1
),
152
164
(
2002
).
56.
H.
Bagci
,
A. E.
Yilmaz
,
V.
Lomakin
, and
E.
Michielssen
, “
Fast solution of mixed-potential time-domain integral equations for half-space environments
,”
IEEE Trans. Geosci. Remote Sens.
43
(
2
),
269
279
(
2005
).
57.
S.
Turley
, “
Acoustic scattering from a sphere
,” in
Class Notes
(
Department of Physics and Astronomy, Brigham Young University
,
Provo, UT
,
2006
).
58.
S.
Yan
,
J.-M.
Jin
, and
Z.
Nie
, “
Accuracy improvement of the second-kind integral equations for generally shaped objects
,”
IEEE Trans. Antennas Propag.
61
(
2
),
788
797
(
2013
).
59.
H.
Bagci
,
F. P.
Andriulli
,
K.
Cools
,
F.
Olyslager
, and
E.
Michielssen
, “
A Calderón multiplicative preconditioner for the combined field integral equation
,”
IEEE Trans. Antennas Propag.
57
(
10
),
3387
3392
(
2009
).