Shape Sensitivity Analysis of Acoustic Scattering with Series Expansion Boundary Element Methods

1  Introduction

In modern urban environments, traffic noise pollution has become a critical environmental issue, significantly impacting daily life, work efficiency, and public health [1,2]. With the rapid advancement of urbanization, traffic volumes on highways, railways, and city roads have surged, exacerbating noise pollution and diminishing residents’ quality of life [3,4]. Prolonged exposure to high noise levels has been linked to serious health risks, including hearing impairment, sleep disturbances, and cardiovascular diseases. Addressing this growing concern necessitates effective noise mitigation strategies and advanced acoustic modeling techniques [5,6].

Among the traffic noise control measures, sound barriers are widely used along highways, railways, and urban roads because of their significant noise reduction effect and convenient construction. Sound barriers mainly block the propagation path of noise and reduce the noise level in the far field area, thereby improving the acoustic environment in the noise-affected area [7]. However, the noise reduction performance of traditional sound barriers is limited by material properties, geometric shapes, and surrounding environmental factors, and it is difficult to maintain a good attenuation effect under different noise frequencies and complex scenarios. In addition, the impact of ground reflection on noise propagation cannot be ignored. The design of optimized sound barriers needs to comprehensively consider ground effects, boundary conditions, and noise propagation characteristics to improve the overall noise reduction capability.

The boundary element method (BEM) has long been recognized as a powerful tool for acoustic analysis, particularly in solving unbounded field problems. Unlike finite element methods (FEM), BEM discretizes only the boundary of the domain, drastically reducing computational effort while maintaining high fidelity in capturing wave phenomena such as diffraction, reflection, and scattering. This feature makes BEM particularly advantageous for road noise prediction and noise barrier optimization. Pioneering work by Rizzo [8] laid the foundation for BEM in elastodynamics, which later expanded into acoustics, as demonstrated by Hothersall et al. [9] in their studies on barrier shape effects and ground reflections. Subsequent advancements, including fast multipole methods (FMM), adaptive integration techniques, and hybrid collocation schemes for coupled vibro-acoustic problems, further enhanced BEM’s applicability to large-scale engineering problems [10,11].

Despite its strengths, conventional BEM faces critical challenges in broadband acoustic computations. The frequency-dependent kernel function necessitates repetitive reconstruction of system matrices across frequency bands, leading to prohibitive computational costs and memory demands for large-scale problems [12–14]. Many approaches have been proposed to address this limitation [15–19], recent innovations integrate Taylor series expansions to decouple the kernel function into frequency-dependent and frequency-independent components. By precomputing and reusing frequency-independent terms, this approach significantly reduces computational overhead while preserving accuracy, thereby extending BEM’s utility to wideband noise barrier optimization [20,21]. Notably, the modified singular boundary method (SBM) introduced by Wu et al. [22] incorporates a combined Helmholtz integral equation formulation and self-regularization technique, which improves the accuracy and efficiency of numerical solutions, especially in external acoustics.

This study integrates Taylor-expanded BEM with reduced-order computation for sensitivity analysis of noise barrier structures. By employing Taylor series expansion, the frequency-dependent kernel is decomposed into separable components, enabling efficient reuse of precomputed terms across broadband frequencies [23–26]. Furthermore, a reduced-order model is adopted to alleviate the computational burden associated with large-scale matrices [27–29]. To overcome singularities inherent in boundary integral equations, Cauchy principal value integration and Hadamard finite-part integration are systematically applied, ensuring numerical stability and accuracy. The proposed isogeometric sensitivity boundary integral equation, based on the direct differentiation method, facilitates efficient sensitivity evaluation with respect to control point coordinates. This methodology not only enhances the computational efficiency of broadband acoustic optimization but also provides a theoretical foundation for acoustic optimization calculations.

Specifically, the research content of this section includes the following aspects:

1. Acoustic modeling of ground reflection: Based on the image source method, the sound field integral equation including the ground reflection effect is established, and its influence on the acoustic response is analyzed.

2. Derivation of acoustic sensitivity formula: Based on the formula in Section 3, the direct differential method derivation considering the ground reflection term is extended to construct a new sensitivity analysis formula.

3. Broadband calculation and order reduction implementation: Taylor expansion is performed on the sensitivity formula under ground reflection conditions, and the system order reduction is completed in combination with the SOAR method to reduce the calculation scale and memory requirements [30,31].

4. Numerical verification and analysis: For typical acoustic problems, numerical simulations including ground reflection conditions are performed to verify the accuracy of the derived formula and the effectiveness of the order reduction method.

2  Acoustic State Analysis Using IGABEM

Isogeometric Boundary Element Method (IGABEM) is an advanced computational approach for acoustic state analysis, leveraging the smooth basis functions of Non-Uniform Rational B-Splines to enhance numerical accuracy and geometric representation [32]. By integrating geometric modeling with numerical simulation, IGABEM reduces discretization errors and improves solution continuity, making it particularly effective for solving Helmholtz equations in unbounded acoustic domains [33–36]. This method enables precise analysis of sound propagation, diffraction, and scattering while maintaining computational efficiency. It is widely applied in acoustic state analysis to evaluate sound field characteristics in complex engineering scenarios, such as noise control and structural acoustics.

2.1 Construction of Model

By combining Isogeometric Analysis (IGA) [37–39] and the Boundary Element Method (BEM), IGABEM utilizes Non-Uniform Rational B-Splines (NURBS) as basis functions for both geometry representation and physical field approximation [40,41]. NURBS are widely used in Computer-Aided Design (CAD) due to their ability to accurately model complex curves and surfaces with fewer control points compared to traditional polynomial-based methods, eliminating the geometric discretization errors that often occur in traditional methods [42]. The rational form of NURBS allows exact representation of conic sections (e.g., circles, ellipses) that cannot be captured by standard polynomials, ensuring geometric fidelity in acoustic domains with curved boundaries. Such preservation is particularly critical for acoustic simulations, where boundary conditions (e.g., impedance or pressure distributions) are highly sensitive to geometric accuracy. In the constructed curve, the coordinates of any arbitrary point can be represented as follows:

x(ξ)=∑i=0nRi,pg(ξ)Pi,(1)

in which the coordinates of the control point are denoted by Pi, and the NURBS basis function is represented by Ri,pg, defined as follows:

Ri,pg(ξ)=Ni,pg(ξ)wiW(ξ),(2)

where Ni,pg(ξ) is the i-th B-spline basis function of degree pg, wi is the weight associated with the control point Pi, W(ξ) is the weighted summation of the B-spline basis functions:

W(ξ)=∑i=0nwiNi,pg(ξ).(3)

For boundary integral calculation in IGABEM, Gauss-Legendre quadrature is needed to handle singular and nearly singular kernels appearing in the Helmholtz integral equation.

2.2 IGABEM Formulations for Acoustic State Analysis with NURBS

This subsection provides a detailed explanation of the isogeometric boundary element formulation used to calculate the radiated and scattered acoustic fields generated by incident waves. Initially, the commonly employed Burton-Miller boundary integral equation is utilized to accurately evaluate the acoustic pressure field at all frequencies, as follows [43,44]:

C(x)[ϕ(x)+αq(x)]+∫S[F(x,y)+αH(x,y)]ϕ(y)dS(y)=∫S[G(x,y)+αK(x,y)]q(y)dS(y)+ϕinc(x)+α∂ϕinc(x)∂n(x),(4)

in which x and y denote the source and field points, respectively. The term C(x) corresponds to the jump condition, which is equal to 1/2 under conditions of smooth boundaries. The parameter α acts as the coupling coefficient, defined by α=i/k when k>1 and by α=i otherwise [45]. The sound pressure flux is represented by q(x)=∂ϕ(x)/∂n(x). ϕinc stands for the sound pressure of the incident wave at the point x. The values for G(x,z) and its normal derivative in the context of 2D full-space acoustics are presented as follows:

{G(x,y)=i4H0(1)(kr)F(x,y)=∂G(x,y)∂n(y)=−ik4H1(1)(kr)∂r∂n(y)K(x,y)=∂G(x,y)∂n(x)=−ik4H1(1)(kr)∂r∂n(x)H(x,y)=∂2G(x,y)∂n(x)∂n(y)=ik4rH1(1)(kr)nj(x)nj(y)−ik24H2(1)(kr)∂r∂n(x)∂r∂n(y),(5)

where r=|x−y|, which represents the distance between the source point and the field point. The term nj is a Cartesian component of n(x) or n(y), and ∂r∂n=r,jnj.

The assumed ideal conditions of the 2D full-space acoustics are often not met in actual engineering. In many application scenarios such as urban noise control, sound barrier design, and architectural acoustic analysis, the sound field is often affected by complex boundary conditions. Ground reflection, as a common boundary effect, will significantly change the distribution characteristics and sensitivity of the sound field. Therefore, ignoring ground reflection can lead to deviations in sensitivity calculations, which can affect the reliability of design and optimization.

In order to be closer to actual engineering needs, we derive the acoustic formula of a two-dimensional half-space, which refers to a scenario that introduces ground reflection conditions in acoustic analysis. Its core feature is that sound wave propagation includes not only the direct path, but also the impact of the ground reflection path on the sound field distribution. In the analysis of urban traffic noise or building noise, the modeling of two-dimensional half-space and ground reflection is equivalent. The study of this model is essentially to accurately describe the sound field distribution and sound pressure changes under complex boundary conditions. The ground reflection effect is particularly significant in the near-surface area. The sound pressure superposition effect caused by it will change the propagation law of noise, which is of great significance to the design of noise barriers and the optimization of the sound environment around buildings. In order to effectively deal with the impact of ground reflection on the sound field distribution, this chapter introduces the mirror point processing method.

The core idea of the mirror point method is to regard the ground as the symmetric boundary of the sound wave, and to construct the mirror point of the sound source below the ground, so that the ground reflection is equivalent to the mirror sound source in the whole space. The position of the mirror point is determined by the symmetry relationship of the original sound source relative to the ground, and its acoustic characteristics are consistent with the original sound source, but the corresponding correction coefficient may be introduced due to the sound absorption characteristics of the ground. This method greatly simplifies the sound field modeling under ground reflection conditions, so that it can be directly applied to the sound field superposition calculation. The fundamental solution for half-space problems, referred to as G~(x,y), can be represented in the following form [45,46]:

G~(x,y)=G(x,y)+γG(x~,y)(6)

As shown in Fig. 1, in Eq. (6), x~ represents the mirror image of x relative to the horizontal plane. γ denotes the reflection coefficient, where γ=1 when the boundary is a rigid plane, and γ=−1 when the boundary is flexible. By replacing G(x,y) in Eq. (4) with x~, it can be rewritten as follows:

C(x)[ϕ(x)+αq(x)]+∫S[F~(x,y)+αH~(x,y)]ϕ(y)dS(y)=∫S[G~(x,y)+αK~(x,y)]q(y)dS(y)+ϕ~inc(x)+α∂ϕ~inc(x)∂n(x),(7)

where

{G~(x,y)=i4H0(1)(kr)+γi4H0(1)(kr~)F~(x,y)=∂G~(x,y)∂n(y)=−ik4H1(1)(kr)∂r∂n(y)−γik4H1(1)(kr~)∂r~∂n(y)K~(x,y)=∂G~(x,y)∂n(x)=−ik4H1(1)(kr)∂r∂n(x)−γik4H1(1)(kr~)∂r~∂n(x)H~(x,y)=∂2G~(x,y)∂n(x)∂n(y)=ik4rH1(1)(kr)nj(x)nj(y)−ik24H2(1)(kr)∂r∂n(x)∂r∂n(y)+γik4r~H1(1)(kr~)nj(x)nj(y)−γik24H2(1)(kr~)∂r~∂n(x)∂r~∂n(y).(8)

Figure 1: The propagation of sound waves emitted by the source (considering ground reflection).

In the framework of IGABEM, the geometry is constructed by discretizing the physical field using NURBS basis functions.

ϕ(ξ)=∑i=0nRi,pg(ξ)ϕ^iq(ξ)=∑i=0nRi,pg(ξ)q^i(9)

Since NURBS fundamental functions don’t have the Kronecker-delta characteristic, ϕ^i and q^i do not represent the actual sound pressure and flux values at the collocation points but rather the coefficients of sound pressure and sound pressure flux.

Substituting Eqs. (1) and (9) into Eq. (7) and evaluating the equation at the discrete collocation point x, the following expression can be obtained:

C(x(ξ¯i))[ϕ(x(ξ¯i))+αq(x(ξ¯i))]+∑e=1Ne∑j=0n{∫ee+1[F~(x(ξ¯i),y(ξ))+αH~(x(ξ¯i),y(ξ))]Rj,pg(ξ)J(ξ)dξ}ϕ^j=∑e=1Ne∑j=0n{∫ee+1[G~(x(ξ¯i),y(ξ))+αK~(x(ξ¯i),y(ξ))]Rj,pg(ξ)J(ξ)dξ}q^j+ϕ~inc(x(ξ¯i))+α∂ϕ~inc(x(ξ¯i))∂n(x(ξ¯i)),(10)

in which subscript i denotes the number of control points. Control points are projected onto the curve to derive the collocation points. Eq. (10) can be rewritten in the following matrix form:

HΦ=Gq+Φinc,(11)

where the matrices H and G are dense, non-symmetric, and frequency-dependent. The vectors Φ and q represent acoustic pressure and acoustic flux, respectively. Φinc denotes the incident wave vector. In the paper, the impedance boundary condition q(x)=ikβ(x)ϕ(x) is introduced, allowing Eq. (11) to be rewritten as follows:

[H−GB]Φ=Φinc,(12)

where

B=ik[β1⋯0⋮⋱⋮0⋯βn],(13)

where the i-th element has a normalized admittance denoted by βi.

From Eq. (8), it can be inferred that a wave number k exists in the kernel function, resulting in Eq. (12) being frequency-dependent. This frequency dependence leads to high computational costs when solving broadband problems, as the coefficient matrix must be recalculated at each frequency point. To address this computational burden, a Taylor series expansion is employed in this study to decouple the kernel function into frequency-dependent and frequency-independent components, as described in the following process.

3  Acoustic Shape Sensitivity Analysis Using IGABEM

Shape sensitivity analysis is a method for evaluating the impact of geometric variations on system performance or objective functions. This approach enables engineers to understand how geometric modifications affect performance metrics, thereby facilitating design optimization and enhancing product performance [47,48]. In noise barrier design, shape sensitivity analysis is employed to evaluate the impact of geometric variations on noise reduction performance, thereby guiding the optimization process to enhance acoustic efficiency. By computing the sensitivities of acoustic wave reflection, diffraction, and absorption with respect to key geometric parameters, such as height, angle, and curvature, engineers can identify the most effective shape features for noise attenuation [49]. It not only improves the performance of noise barriers but also reduces material usage and costs, achieving both economic and environmental benefits.

By directly differentiating Eq. (7), the expression for acoustic shape sensitivity considering ground reflections can be derived:

C(x)[ϕ˙(x)+αq˙(x)]+∫S[F~˙(x,y)+αH~˙(x,y)]ϕ(y)dS(y)+∫S[F~(x,y)+αH~(x,y)]ϕ˙(y)dS(y) +∫S[F~(x,y)+αH~(x,y)]ϕ(y)dS˙(y)=∫S[G~˙(x,y)+αK~˙(x,y)]ϕ(y)dS(y) +∫S[G~(x,y)+αK~(x,y)]ϕ˙(y)dS(y)+∫S[G~(x,y)+αK~(x,y)]ϕ(y)dS˙(y) + ϕ~˙inc(x)+α∂ϕ~˙inc(x)∂n(x)(14)

The fundamental solution and its derivatives are governed by the coordinates of both the field and source points. As the shape undergoes continuous modifications, their values may be affected by variations in the shape design variables. Consequently, the sensitivities of the coordinates, represented by G~˙, F~˙, K~˙ and H~˙, can be expressed as follows:

{G~˙(x,y)=−ik4H1(1)(kr)r˙−γik4H1(1)(kr~)r~˙F~˙(x,y)=−ik4H1(1)(kr)[r˙r∂r∂n(y)+(∂r∂n(y))]−ik24H2(1)(kr)r˙∂r∂n(y)−γik4H1(1)(kr~)[r~˙r~∂r~∂n(y)+(∂r~∂n(y))]−γik24H2(1)(kr~)r~˙∂r~∂n(y)K~˙(x,y)=−ik4H1(1)(kr)[r˙r∂r∂n(x)+(∂r∂n(x))]−ik24H2(1)(kr)r˙∂r∂n(x)−γik4H1(1)(kr~)[r~˙r~∂r~∂n(x)+(∂r~∂n(x))]−γik24H2(1)(kr~)r~˙∂r~∂n(x)H~˙(x,y)=ik4rH1(1)(kr)[n˙j(x)nj(y)+nj(x)n˙j(y)]+ik34H3(1)(kr)r˙r,jnj(x)r,lnl(y)−ik24H2(1)(kr)nj(x)nj(y)r˙r−ik24H2(1)(kr)r,jnj(x)[r˙,lnl(y)+r,ln˙l(y)]−ik24H2(1)(kr)r,lnl(y)[r˙,lnj(x)+r,ln˙j(x)]−ik24H2(1)(kr)2r˙r,jnj(x)r,lnl(y)r+γik4r~H1(1)(kr~)[n˙j(x)nj(y)+nj(x)n˙j(y)]+γik34H3(1)(kr~)r~˙r~,jnj(x)r~,lnl(y)−γik24H2(1)(kr~)nj(x)nj(y)r~˙r~−γik24H2(1)(kr~)r~,jnj(x)[r~˙,lnl(y)+r~,ln˙l(y)]−γik24H2(1)(kr~)r~,lnl(y)[r~˙,lnj(x)+r~,ln˙j(x)]−γik24H2(1)(kr)2r~˙r~,jnj(x)r~,lnl(y)r~(15)

in which

r˙=r,j(y˙j−x˙j)(16)

r~˙=r~,i(y˙i−x~˙i)(17)

r~˙,i=(y˙i−x~˙i)r~−(yi−x~)r~˙r~2(18)

(∂r~∂n)⋅=r~,ini=r~˙,ini+r~,in˙i(19)

r˙,lnl(y)+r,ln˙l(y)=(y˙l−x˙l)nl(y)r−r˙r,lnl(y)r+r,jn˙l(y)(20)

where ()˙ represents the differential of the design variable. nl˙(y) and dS˙(y) are expressed as follows [50]:

n˙l(y)=−y˙j,lnj(y)+−y˙j,mnj(y)nm(y)nl(y)(21)

dS˙(y)=[y˙l,l−y˙l,jnl(y)nj(y)]dS(y)(22)

The sensitivity of boundary physical quantities is expressed using NURBS function interpolation as follows:

ϕ˙(ξ)=∑i=0nRi,pg(ξ)ϕ^˙(ξ)q˙(ξ)=∑i=0nRi,pg(ξ)q^˙(ξ)(23)

It is important to note that ϕ^˙(ξ) and q^˙(ξ) in the Eq. (23) do not represent the actual sensitivity values of the physical quantities on the boundary. Instead, they are referred to as the coefficients of boundary physical quantity sensitivity. Meanwhile, dS(y) can be expressed as follows:

dS(y)=J(ξ)dξdS˙(y)=J˙(ξ)dξ(24)

By substituting Eqs. (1), (9), (23), and (24) into Eq. (14), the isogeometric discrete form of the sensitivity boundary integral equation is obtained as follows:

C(x)ϕ˙(x)+αC(x)q˙(x)+∑e=1Ne∑j=0n∫S[F~˙(x,y)+αH~˙(x,y)]Rj,pg(x)ϕ(x)dS(y)+∑e=1Ne∑j=0n∫S[F~(x,y)+αH~(x,y)]Rj,pg(x)ϕ˙(x)dS(y)+∑e=1Ne∑j=0n∫S[F~(x,y)+αH~(x,y)]Rj,pg(x)ϕ(x)dS˙(y)=∑e=1Ne∑j=0n∫S[G~˙(x,y)+αK~˙(x,y)]Rj,pg(x)ϕ(x)dS(y)+∑e=1Ne∑j=0n∫S[G~(x,y)+αK~(x,y)]Rj,pg(x)ϕ˙(x)dS(y)+∑e=1Ne∑j=0n∫S[G~(x,y)+αK~(x,y)]Rj,pg(x)ϕ(x)dS˙(y)+ϕ~˙inc(x)+α∂ϕ~˙inc(x)∂n(x)(25)

To determine all unknown boundary state values required for shape sensitivity analysis, solve Eq. (12) and use all values on the boundary and known boundary sensitivity values to solve Eq. (25). According to equation Eq. (14), let α=0,C(x)=1, the value of any point x in the problem domain can be calculated. Finally, the sensitivity between the objective function and the shape design variables of certain calculation points in the acoustic scattering domain can be calculated. Finally, the sensitivity of the objective function of certain calculation points in the acoustic scattering domain to the shape design variables can be calculated.

4  Rapid Frequency Sweep Analysis of Sound Pressure

The kernel function contains the wave number k, which causes Eq. (11) to be frequency-dependent. This frequency dependence significantly increases the computational cost when solving broadband problems because the coefficient matrix needs to be recalculated at each frequency [51]. To address this computational burden, this chapter uses the Taylor series expansion method to decouple the kernel function into frequency-dependent and frequency-independent terms; and uses the SOAR method to construct a ROM that retains the key features of the original FOM to accelerate the solution process, thereby significantly improving computational efficiency [52–54].

4.1 Frequency Sweep Analysis for Acoustic Sensitivity

The first-order n-th Hankel function can be expressed as a Taylor series expansion at the fixed frequency expansion point k0r [55].

Hn(1)(kr)=∑m=0∞(kr−k0r)mm![Hn(1)(kr)]kr=k0r,(m)(26)

in which k0 represents the fixed frequency expansion point. As shown in Fig. 2, the midpoint of the frequency range is used as the fixed expansion point for this interval.

Figure 2: Calculate the frequency band.

By substituting Eq. (26) and Eq. (15) into Eq. (14), the integral form for kernel function at the fixed frequency expansion point k0 is obtained:

∫SG~˙(x,y)q(y)dS(y)=∑m=0∞(k−k0)mm!k(Sg~1m+Sg~2m)∫SG~(x,y)q˙(y)dS(y)=∑m=0∞(k−k0)mm!(Sg~3m+Sg~4m)∫SG~(x,y)q(y)dS~˙(y)=∑m=0∞(k−k0)mm!(Sg~5m+Sg~6m)(27)

∫SF~˙(x,y)ϕ(y)dS(y)=∑m=0∞(k−k0)mm![k(Sf~1m+Sf~2m)+k2(Sf~3m+Sf~4m)]∫SF~(x,y)ϕ˙(y)dS(y)=∑m=0∞(k−k0)mm!(Sf~5m+Sf~6m)∫SF~(x,y)ϕ(y)dS˙(y)=∑m=0∞(k−k0)mm!(Sf~7m+Sf~8m)(28)

∫SαK˙(x,y)q(y)dS(y)=∑m=0∞(k−k0)mm![k(Sk~1m+Sk~2m)+k2(Sk~3m+Sk~4m)]∫SαK(x,y)q˙(y)dS(y)=∑m=0∞(k−k0)mm!k(Sk~5m+Sk~6m)∫SαK(x,y)q(y)dS˙(y)=∑m=0∞(k−k0)mm!k(Sk~7m+Sk~8m)(29)

∫SαH~˙(x,y)ϕ(y)dS(y)=∑m=0∞(k−k0)mm

where

Sg~1m=−∫Sirm4[H1(1)(kr)]kr=k0r(m)r˙q(y)dS(y)Sg~2m=−∫Sγir~m4[H1(1)(kr~)]kr~=k0r~(m)r~˙q(y)dS(y)Sg~3m=∫Sirm4[H0(1)(kr)]kr=k0r(m)q˙(y)dS(y)Sg~4m=∫Sγir~m4[H0(1)(kr~)]kr~=k0r~(m)q˙(y)dS(y)Sg~5m=∫Sirm4[H0(1)(kr)]kr=k0r(m)q(y)dS˙y)Sg~6m=∫Sγir~m4[H0(1)(kr~)]kr~=k0r~(m)q(y)dS˙(y)(31)

Sf~1m=−∫Sirm4[H1(1)(kr)]kr=k0r(m)[(y˙j−x˙j)nj(y)r+r,jn˙j(y)]ϕ(y)dS(y)Sf~2m=−∫Sγir~m4[H1(1)(kr~)]kr~=k0r~(m)[(y˙j−x~˙j)nj(y)r~+r~,jn˙j(y)]ϕ(y)dS(y)Sf~3m=∫Sirm4[H2(1)(kr)]kr=k0r(m)rr˙,jnj(y)ϕ(y)dS(y)Sf~4m=∫Sγ1˙r~m4[H2(1)(kr~)]kr~=k0r~(m)r~˙r~,jnj(y)ϕ(y)dS(y)Sf~5m=−∫Sirm−14[(kr)H1(1)(kr)]kr=k0r(m)∂r∂n(y)ϕ˙(y)dS(y)Sf~6m=−∫Sγir~m−14[(kr~)H1(1)(kr~)]kr~=k0r~(m)∂r~∂n(y)ϕ˙(y)dS(y)Sf~7m=−∫Sirm−14[(kr)H1(1)(kr)]kr=k0r(m)∂r∂n(y)ϕ(y)dS˙(y)Sf~8m=−∫Sγir~m−14[(kr~)H1(1)(kr~)]kr~=k0r~(m)∂r~∂n(y)ϕ(y)dS˙(y)(32)

Sk~1m=−∫Sαirm4[H1(1)(kr)]kr=k0r(m)[(y˙j−x˙j)nj(x)r+r,jn˙j(x)]q(y)dS(y)Sk~2m=−∫Sγαir~m4[H1(1)(kr~)]kr~=k0r~(m)[(y˙j−x~˙j)nj(x)r~+r~,jn˙j(x)]q(y)dS(y)Sk~3m=∫Sαirm4[H2(1)(kr)]kr=k0r(m)r˙r,jnj(x)q(y)dS(y)Sk~4m=∫Sγαir~m4[H2(1)(kr~)]kr~=k0r~(m)r~˙r~,jnj(x)q(y)dS(y)Sk~5m=−∫Sαirm4[H1(1)(kr)]kr=k0r(m)∂r∂n(x)q˙(y)dS(y)Sk~6m=−∫Sγαir~m4[H1(1)(kr~)]kr~=k0r~(m)∂r~∂n(x)q˙(y)dS(y)Sf~7m=−∫Sαirm4[H1(1)(kr)]kr=k0r(m)∂r∂n(x)q(y)dS˙(y)Sf~8m=−∫Sγαir~m4[H1(1)(kr~)]kr~=k0r~(m)∂r~∂n(x)q(y)dS˙(y)(33)