Desingularization of matrix equations employing hypersingular integrals in boundary element methods using double nodes

TL;DR

This paper addresses singularities in boundary element method matrix equations caused by double nodes at corners, proposing hypersingular integral equations to improve accuracy and reduce computational cost.

Contribution

It introduces a novel application of hypersingular integral equations (HBIE) to resolve singularities in boundary element methods with double nodes, comparing partial and full HBIE approaches.

Findings

Partial-HBIE improves accuracy over full-HBIE.

Both HBIE approaches resolve singularities effectively.

Partial-HBIE has lower computational cost.

Abstract

In boundary element methods, the method of using double nodes at corners is a useful approach to uniquely define the normal direction of boundary elements. However, matrix equations constructed by conventional boundary integral equations (CBIE) become singular under certain combinations of double node boundary conditions. In this paper, we analyze the singular conditions of the CBIE formulation for cases where the boundary conditions at the double node are imposed by combinations of Dirichlet, Neumann, Robin, and interface conditions. To address this singularity we propose the use of hypersingular integral equations (HBIE) for wave propagation problems that obey Helmholtz equation. To demonstrate the applicability of HBIE, we compare three types of simultaneous equations: (i) CBIE, (ii) partial-HBIE in which HBIE is only applied to the double nodes at corners while CBIE is applied to the other nodes, and (iii) full-HBIE in which HBIE is applied to all nodes. Based on our numerical results, we observe the following results. The singularity of the matrix equations for problems with any combination of boundary conditions can be resolved by both full-HBIE and partial-HBIE, and partial-HBIE exhibits better accuracy than full-HBIE. Furthermore, the computational cost of partial-HBIE is smaller than that of full-HBIE.

Equations

$${\nabla }^{2}u(\mathbf{r})+k^{2}u(\mathbf{r})=0,\mathbf{r}\in \Omega ,$$

$${\nabla }^2{u}^{\ast }(\mathbf{r},{\mathbf{r}}_i)+k^2{u}^{\ast }(\mathbf{r},{\mathbf{r}}_i)=-\delta (\mathbf{r}-{\mathbf{r}}_i),$$

$$u^{\ast }(\mathbf{r},{\mathbf{r}}_i)=\frac{1}{4\mathrm{j}}{H}_0^{(2)}(kr),r=|\mathbf{r}-{\mathbf{r}}_i|,$$

$$c({\mathbf{r}}_i)$$

$$c({\mathbf{r}}_i)u({\mathbf{r}}_i)\triangleq {\int }_{\Omega }u(\mathbf{r})\delta (\mathbf{r}-{\mathbf{r}}_i)d\Omega ={\int }_{\Omega }\delta (\mathbf{r}-{\mathbf{r}}_i)d\Omega u({\mathbf{r}}_i).$$

$${\nabla }^2{u}_L^{\ast }(\mathbf{r},{\mathbf{r}}_i)=-\delta (\mathbf{r}-{\mathbf{r}}_i),$$

$$c({\mathbf{r}}_i)={\int }_{\Omega }\delta (\mathbf{r}-{\mathbf{r}}_i)d\Omega =-{\oint }_{\Gamma }\nabla u_L^{\ast }\cdot \mathbf{n}d\Gamma ,$$

$$\sum _{j\in I}\left(c_i{\delta }_i^j+h_i^j\right)u_j-\sum _{j\in I}{g}_i^j{q}_j=0,$$

$$\sum _{j\in I}{a}_i^j{x}_j={\overline{b}}_i,{\overline{b}}_i=\sum _{j\in I}{b}_i^j{\overline{x}}_j\text{mbox}foreachi\in I,$$

$$q(\mathbf{r})=\sum _{l=1}^{N_l}{\varphi }^{(k,l)}(\mathbf{r})q_{(k,l)}\text{mbox}on{\Gamma }_k,$$

$$q_{(k,l)}=\sum _j{m}_{(k,l)}^j{q}_j,$$

$${\oint }_{\Gamma }{u}^{\ast }(\mathbf{r},{\mathbf{r}}_i)q(\mathbf{r})d\Gamma$$

$$=\sum _k{\int }_{{\Gamma }_k}{u}^{\ast }(\mathbf{r},{\mathbf{r}}_i)\left(\sum _l{\varphi }^{(k,l)}(\mathbf{r})\left(\sum _j{m}_{(k,l)}^j{q}_j\right)\right)d\Gamma$$

$$=\sum _j\left(\sum _k\sum _l{m}_{(k,l)}^j{g}_i^{(k,l)}\right)q_j=\sum _j{g}_i^j{q}_j,$$

$$g_i^{(k,l)}={\int }_{{\Gamma }_k}{u}^{\ast }(\mathbf{r},{\mathbf{r}}_i){\varphi }^{(k,l)}(\mathbf{r})d\Gamma ,$$

$$g_i^j=\sum _k\sum _l{m}_{(k,l)}^j{g}_i^{(k,l)}.$$

$${\oint }_{\Gamma }{q}^{\ast }(\mathbf{r},{\mathbf{r}}_i)u(\mathbf{r})d\Gamma =\sum _j{h}_i^j{u}_j,$$

$$h_i^{(k,l)}={\int }_{{\Gamma }_k}{q}^{\ast }(\mathbf{r},{\mathbf{r}}_i){\varphi }^{(k,l)}(\mathbf{r})d\Gamma ,$$

$$h_i^j=\sum _k\sum _l{m}_{(k,l)}^j{h}_i^{(k,l)},$$

$$c_{\alpha }=c_{\beta }.$$

$$h_{\alpha }^{\alpha }=h_{\beta }^{\alpha }=h_{\alpha }^{\beta }=h_{\beta }^{\beta }=0$$

$$g_{\alpha }^j=g_{\beta }^j$$

$$h_{\alpha }^j=h_{\beta }^j$$

$$u_j={\overline{u}}_j$$

$$q_j={\overline{q}}_j$$

$${\kappa }_u{u}_j+{\kappa }_q{q}_j={\overline{\psi }}_j$$

$$u_{j^{(2)}}=u_{j^{(1)}},q_{j^{(2)}}=-{\kappa }_{21}{q}_{j^{(1)}}$$

$$-g_{\alpha }^{\alpha }{q}_{\alpha }-g_{\alpha }^{\beta }{q}_{\beta }+\sum _{j\in I\setminus \{\alpha ,\beta \}}{a}_{\alpha }^j{x}_j={\overline{b}}_{\alpha },$$

$$-g_{\beta }^{\alpha }{q}_{\alpha }-g_{\beta }^{\beta }{q}_{\beta }+\sum _{j\in I\setminus \{\alpha ,\beta \}}{a}_{\beta }^j{x}_j={\overline{b}}_{\beta },$$

$$\sum _{j\in I\setminus \{\alpha ,\beta \}}{a}_i^j{x}_j=\sum _{j\in I_N}{h}_i^j{u}_j-\sum _{j\in I_D,j\ne \alpha ,\beta }{g}_i^j{q}_j,$$

$${\overline{b}}_i=-\sum _{j\in I_D}\left(c_i{\delta }_i^j+h_i^j\right){\overline{u}}_j+\sum _{j\in I_N}{g}_i^j{\overline{q}}_j.$$

$$\left(c_{\alpha }+h_{\alpha }^{\alpha }\right)u_{\alpha }+h_{\alpha }^{\beta }{u}_{\beta }$$

$$h_{\beta }^{\alpha }{u}_{\alpha }+\left(c_{\beta }+h_{\beta }^{\beta }\right)u_{\beta }$$

$$-g_{\alpha }^{\alpha }{q}_{\alpha }$$

$$-g_{\beta }^{\alpha }{q}_{\alpha }$$

$$\sum _{j\in I_N}\left(c_i{\delta }_i^j+h_i^j\right)u_j-\sum _{j\in I_D}{g}_i^j{q}_j-\sum _{j\in I_R}\left(\frac{{\kappa }_q}{{\kappa }_u}\left(c_i{\delta }_i^j+h_i^j\right)+g_i^j\right)q_j={\overline{b}}_i.$$

Full text

††journal: Engineering Analysis with Boundary Elements (Accepted June 11, 2019)

1 Introduction

The boundary element method (BEM), the finite difference method (FDM), and the finite element method (FEM) have been commonly used to solve boundary value problems. In the BEM, a set of simultaneous equations for determining unknown variables at nodes on the boundary is constructed in discretized boundary integral equations. The variables in simultaneous equations are nodal field values and normal derivatives only on individual boundary nodes, whereas the variables in the FDM or FEM are field values at domain nodes which are placed in the entire domain enclosed by the boundary. Therefore, the number of variables in the BEM is much smaller than that in the FDM or FEM, which is one of the advantages of BEM. Furthermore, the BEM can be easily applied to external problems, such as wave scattering problems, since it does not require the placement of nodes in a domain spreading to infinity and it does not require any other boundary conditions to represent radiation at a boundary enclosing the domain considered.

In most boundary problems, the field values, $u$, along the boundary are continuous; however, the normal derivatives, $q$, are discontinuous at corners since the normal directions at any corner point are different. By using a linear element or higher-order elements, the boundary elements share the nodes at both ends of the element with adjacent boundary elements. In this case, the normal direction, $\boldsymbol{n}$, at the corner node cannot be defined uniquely since the single node at any corner belongs to two boundary elements with different normal directions.

There are two approaches to addressing the definition problem of the normal direction. The first approach involves the use of non-conforming elements, which are also called discontinuous elements. In the non-conforming element, collocation nodes that represent $u$ and $q$ do not coincide with geometric nodes, but they do so in the conforming element. The non-conforming element has been applied to several problems; e.g., elastostatic problems [1, 2, 3, 4, 5, 6, 7], fluid flow problems [8, 9], and acoustic problems [10]. Although the accuracy between the non-conforming element and the conforming element were compared by Manolis and Banerjee [1], and Parreira [2], they arrived at different conclusions.

The second approach includes a double node technique [11] or a multiple node technique for three-dimensional problems. In the double node technique, two normal derivatives, $q_{\alpha}=\boldsymbol{n}_{\alpha}\cdot\nabla u$ and $q_{\beta}=\boldsymbol{n}_{\beta}\cdot\nabla u$, and a field, $u$, are defined at the corner node, where $\boldsymbol{n}_{\alpha}$ and $\boldsymbol{n}_{\beta}$ denote the directions normal to the two boundary elements connected to the corner node. However, a set of simultaneous equations, called a matrix equation, becomes singular under certain boundary conditions; i.e., the rank of the matrix equation is reduced since some node equations are redundant. Further details will be presented in Sec. 2.

To address the rank reduction problem caused by the double nodes, there are two categories of approaches. The first category involves the use of a local relation for each double node [12, 13, 14, 15]. By using a Taylor expansion around the corner node, this relation is described as a linear combination of the two normal derivatives and the field values at neighbor nodes of the corner node. The local relation is replaced by one of the redundant equations that reduces the rank; therefore, the matrix equation still includes a square matrix. In the second category, extra node equations are employed [16, 17, 18, 19, 20]. Mitra and Ingber [16] proposed a technique for replacing one of the redundant equations in each corner by an extra node equation with respect to an extra collocation node placed outside the domain considered. Following this study, the authors mentioned that “external collocation yields a coefficient matrix with a large number” [17], and they improved the method using the extra node equations so that the location of the extra node on the boundary elements connects to the double node [17, 18]. Subia, Ingber, and Mitra demonstrated that there are no significant differences in accuracy between the method of the extra node equation and the non-conforming method [18]. The method that uses the extra node equations was extended to the problems of interface boundaries at which two or more domains are connected [19, 20]. Using these methods, the number of variables is not increased since the field at the extra node is known. Therefore, the coefficient matrix is the square matrix, which is similar to the matrix of the first category. If we simply add the local relations shown in the first category or the extra node equations in the second category instead of replacing them, the number of equations becomes larger than the number of variables, which is referred to as an overdetermined problem. We can solve this equation by using least-square methods, but the computational time required to solve this equation is much greater than solving the general simultaneous equation; therefore, the replacements are generally applied. The replacement of the equation should be performed individually while examining the types of boundary conditions at the corner. The individual examination increases the complexity of the programming of widely applicable BEM codes that include many types of boundary conditions.

In addition to the corner node problems, there are rank deficient or large condition number problems. We focus on two problems related to hypersingular boundary integral equations. The first is a non-uniqueness or a spurious solution problem in an external field for a wave scattering problem that obeys the Helmholtz equation, in which the domain considered is outside of the boundary enclosing a scatterer. In this situation, spurious solutions can be obtained at the eigenfrequency of the scatterer. To resolve this problem, two major approaches were proposed. Schenck employed equations related to additional nodes in the scatterer region and the method is called ‘Combined Helmholtz Integral Equation Formulation (CHIEF)’ [21]. Chen et al. [22] applied the additional node equations to interior problems in which additional points are placed outside of the boundary. These methods are similar to the extra node equation methods for the corner problems shown in the previous paragraph; however, they require a solver based on the least-square method since the set of simultaneous equations becomes an overdetermined equation. The other approach is referred to as the Burton-Miller method [23, 24, 25, 26]. Burton and Miller [23] represented a boundary integral equation (BIE) for each node using a linear combination of two types of BIE called a conventional BIE (CBIE) and a hypersingular BIE (HBIE). In the CBIE, the field $u$ at the field point is denoted by integrals over the boundary on which sources are distributed. The HBIE is obtained by taking a gradient of the CBIE. The singularity of the integral in the HBIE is stronger than that of the CBIE; therefore, it is called a hypersingular integral. Bentihien and Schenck [24] reviewed the non-uniqueness problem with comparisons of other methods, including the CHIEF and Burton-Miller methods.

The other rank deficient problem is found at a degenerate boundary which appears either at a crack in an elastostatic problem [27, 28, 29, 30] or at both surfaces of a thin metal with zero-thickness in an electromagnetic problem [29]. To resolve these problems, the CBIEs are applied to one side of the crack or the thin metal, and the HBIEs are applied to the other side. These methods are referred to as dual-BEM.

As shown in the previous two paragraphs, the use of HBIEs is effective for resolving rank deficient problems. In this paper, we will demonstrate that the corner singular problem for wave propagation problems can be solved by only using node equations based on HBIEs. In addition, to suppress the rank deficiency caused by the corners, we do not need to use HBIEs for every node. We will also illustrate that only the replacement of the node equations related to the corner node by HBIEs is sufficient, which is similar to the aforementioned dual-BEM.

In the application of HBIEs, the regularization of hypersingularity is a key issue. The authors developed an analytical regularization of the two-dimensional Helmholtz equation [31]; the other regularization methods are also found in the references of that study. In the regularization of HBIEs, we include the relations between two normal derivatives at the double nodes, similar to the local relation methods described above. Therefore, the method by HBIE can be considered an extension that uses the local relation. However, the contributions to the double node from all the nodes are considered in the method by HBIE; whereas the local relation method represents the relations between local nodal quantities only.

The method proposed in this paper to overcome the corner problem caused by the double nodes does not require local relations at corners, extra node equations, or least-square methods. Our method can also be applied to any kind of boundary conditions of boundary elements that include corner nodes. In addition, the proposed method can also be applied to interface boundary conditions. From these characteristics, no additional effort is required in prepossessing to prepare the input data for solving the boundary value problem.

The outline of this paper is as follows. In Sec. 2, we demonstrate why the rank of the coefficient matrix of the CBIE is reduced in the case where the double node is employed. We also demonstrate the condition that results in rank deficiency based on the nature of the discretized node equations. In Sec. 3, we illustrate why the HBIE does not cause a rank deficiency. The numerical results and discussions for simple waveguide problems are presented in Sec. 4 to demonstrate that the HBIE is applicable to any combination of boundary conditions. The advantages of the partial-HBIE method in which the HBIE is applied to only the double nodes and the CBIE is applied to other nodes are also shown. Finally, the summary is presented in Sec. 5.

2 Rank deficiency problem in CBIEs

To illustrate rank deficient conditions for a set of discretized node equations of CBIEs, the discretization using linear elements is first shown; the double node technique that defines two sub-nodes for a double node are then shown; and lastly, the rank deficiency conditions are discussed based on comparisons between two equations for the two sub-nodes.

2.1 Discretization of CBIEs

A complex-valued time harmonic scalar wave satisfies the following Helmholtz equation with an assumed time factor ${\rm e}^{{\rm j}\omega t}$, where $\rm j$ is an imaginary unit and $\omega$ is the angular frequency:

[TABLE]

where $k$ indicates the wave number, which is a ratio of $\omega$ to the wave velocity, and $\Omega$ represents the spatial domain considered. A fundamental solution $u^{*}(\boldsymbol{r},\boldsymbol{r}_{i})$, which represents the contribution to a field point $\boldsymbol{r}_{i}$ from a unit source placed at a source point $\boldsymbol{r}$ in free space satisfies

[TABLE]

where the differential operator $\nabla$ operates only on $\boldsymbol{r}$, but not on $\boldsymbol{r}_{i}$. This equation can be solved analytically. In two-dimensional problems, the outward propagating wave that obeys Eq. (2) is

[TABLE]

where the function $H_{0}^{(2)}(kr)$ is a second kind 0-th order Hankel function.

Using Green’s second identity and some integral operations for Eqs. (1) and (2), we obtain the conventional boundary integral equation (CBIE),

[TABLE]

where $\Gamma$ denotes a boundary surrounding $\Omega$, $\boldsymbol{r}$ is the position of the source point on the boundary, $\boldsymbol{r}_{i}$ is the position of the field point, $\boldsymbol{n}$ is the outward-pointing normal unit vector; the contour integral is evaluated as a Cauchy principal value, and $c(\boldsymbol{r}_{i})$ is the result of the following evaluation of Dirac’s delta function:

[TABLE]

The coefficient $c(\boldsymbol{r}_{i})$ depends on the shape of the boundary $\Gamma$ at the field point $\boldsymbol{r}_{i}$. Because of the nature of Dirac’s delta function, when $\boldsymbol{r}_{i}$ is located inside and outside the domain, $c(\boldsymbol{r}_{i})$ evaluates to 1 and 0, respectively. In the case where $\boldsymbol{r}_{i}$ is located on $\Gamma$, $c(\boldsymbol{r}_{i})$ is equal to the ratio of the interior angle $\varDelta\theta$ to the whole angle; e.g., $\varDelta\theta/2\pi$ for 2-dimensional problems. In addition, when we introduce a fundamental solution to a Laplace equation, $u_{L}^{*}$, which satisfies

[TABLE]

the coefficient $c(\boldsymbol{r}_{i})$ can be expressed by a boundary integral:

[TABLE]

which is called an equipotential condition [32].

When $\boldsymbol{r}_{i}$ is located on $\Gamma$ and the boundary is partitioned into a number of boundary elements, Eq. (2.1) can be written as a node equation for the node $i$ as a discrete algebraic expression:

[TABLE]

where both the superscripts and the subscripts $i$ and $j$ denote the node numbers and not the element numbers; $I$ denotes a set of the node numbers, which includes $N$ members; $\delta_{i}^{j}$ denotes the Kronecker’s delta; and $q\equiv\nabla u\cdot\boldsymbol{n}$. The factors $g_{i}^{j}$ and $h_{i}^{j}$ are results of boundary integrals in Eq. (2.1) as shown below.

In Eq. (8), there are two variables, $u_{j}$ and $q_{j}$, at $\boldsymbol{r}_{j}$. One of them is specified by a boundary condition as a boundary value which is denoted by $\overline{x}_{j}$, and the other is an unknown variable denoted $x_{j}$; therefore, the number of variables is equal to that of the nodes, $N$. Since the node $i$ can be placed at every boundary node, we can obtain $N$ equations of Eq. (8) as

[TABLE]

where both $a_{i}^{j}$ and $b_{i}^{j}$ are either $h_{i}^{j}$ or $g_{i}^{j}$ when simple boundary conditions are specified. The matrix composed of $a_{i}^{j}$ is called a coefficient matrix in the following discussions.

The coefficient $g_{i}^{j}$ represents the contribution to $u_{i}$ from the nodal quantity $q_{j}$ at the node $j$, and $h_{i}^{j}$ represents the contribution to $u_{i}$ from $u_{j}$. These coefficients are evaluated by boundary integrals, which depend on a method that interpolates $u(\boldsymbol{r})$ and $q(\boldsymbol{r})$ along the boundary element. To examine the nature of $g_{i}^{j}$ and $h_{i}^{j}$, we present evaluations using a shape function to interpolate them. Let us consider the discretization of the boundary integral of $u^{*}(\boldsymbol{r},\boldsymbol{r}_{i})\,q(\boldsymbol{r})$. By denoting the $k$-th boundary element as $\Gamma_{k}$ and a local node number as $l$, $q(\boldsymbol{r})$ on $\Gamma_{k}$ is represented by a linear combination of shape functions $\phi^{(k,l)}(\boldsymbol{r})$ and $q$ at the boundary nodes:

[TABLE]

where $N_{l}$ denotes the number of nodes in a boundary element; e.g., $N_{l}=2$ for the linear element, $N_{l}=3$ for the second-order element. A global node number $j$ can be mapped from the local node numbers $(k,l)$ by a permutation matrix $m_{(k,l)}^{j}$:

[TABLE]

where $m^{j}_{(k,l)}$ has the value 1 if $(k,l)$ and $j$ are associated, and 0 otherwise. Although the notation of $m^{j}_{(k,l)}$ is used here, it provides a symbolic meaning, and it is treated as a mapping function such as $j=m(k,l)$ in actual coding to avoid summation procedures with respect to $k$. By using these definitions, the boundary integral is discretized as

[TABLE]

where $g_{i}^{(k,l)}$ and $g_{i}^{j}$ are defined as

[TABLE]

The integral on the right-hand side in Eq. (13) is evaluated by a numerical integral such as Gauss quadrature or by an analytical integral in the case of $\boldsymbol{r}_{i}\in\Gamma_{k}$, which is called a singular integral. Similarly, the other integral in Eq. (2.1) is evaluated as

[TABLE]

where $q^{*}(\boldsymbol{r},\boldsymbol{r}_{i})=\nabla u^{*}(\boldsymbol{r},\boldsymbol{r}_{i})\cdot\boldsymbol{n}(\boldsymbol{r})$, which is not the derivative at $\boldsymbol{r}_{i}$; i.e., $q^{*}(\boldsymbol{r},\boldsymbol{r}_{i})$ does not depend on $\boldsymbol{n}(\boldsymbol{r}_{i})$.

To discuss the rank of the matrix equation shown in Eq. (9), the dependencies of both $g_{i}^{j}$ and $h_{i}^{j}$ on $\boldsymbol{n}(\boldsymbol{r}_{i})$ which is the normal unit vector at $\boldsymbol{r}_{i}$ are important. First, $g_{i}^{j}$ is independent of $\boldsymbol{n}(\boldsymbol{r}_{i})$ since Eqs. (13) and (14) do not include $\boldsymbol{n}(\boldsymbol{r}_{i})$ regardless of $j$ or $k$. In contrast, the case of $h_{i}^{j}$ is different from that of $g_{i}^{j}$. Although $h_{i}^{j}$ in Eqs. (16) and (17) looks independent of $\boldsymbol{n}(\boldsymbol{r}_{i})$ even when $q^{*}(\boldsymbol{r},\boldsymbol{r}_{i})$ is expressed using $\boldsymbol{n}(\boldsymbol{r})$, there is an exception at $i=j$ in which $\boldsymbol{n}(\boldsymbol{r})$ becomes $\boldsymbol{n}(\boldsymbol{r}_{i})$. The dependency of this exceptional case is presented in the next section.

2.2 Double node technique

In the node equation for node $i$ shown in Eq. (8), both $u_{j}$ and $q_{j}$ are values of the nodes, which are located at the two ends of the boundary element in the case of the linear element. When the node $\boldsymbol{r}_{i}$ is located at a corner, the normal derivative $q_{i}$ of the linear element cannot be defined uniquely since the node belongs to two elements with different normal directions. Using a double node technique is one of the solutions.

Figure 1 illustrates a configuration of a double node. The node $\boldsymbol{r}_{i}$ is represented by two sub-nodes at the same position; $\boldsymbol{r}_{\alpha}=\boldsymbol{r}_{\beta}$. The node at $\boldsymbol{r}_{\alpha}$ is connected to both $\Gamma_{{\rm A}}$ and a zero-sized element between $\boldsymbol{r}_{\alpha}$ and $\boldsymbol{r}_{\beta}$, and vice versa. The normal direction of the sub-nodes $\alpha$ and $\beta$ can be determined from the directions normal to $\Gamma_{{\rm A}}$ and $\Gamma_{{\rm B}}$, respectively; therefore, the direction of the derivatives of $q_{\alpha}$ and $q_{\beta}$ can be defined individually. Since the integral along the zero-sized element is identically zero regardless of the normal direction, we do not need to define the normal direction for the zero-sized element. Although the node $\alpha$ does not belong to $\Gamma_{\rm B}$, contributions of $g^{({\rm B},l)}_{\alpha}$ and $h^{({\rm B},l)}_{\alpha}$, which are evaluated by analytical integrals as singular integrals for accurate evaluation, are similar to $g^{({\rm A},l)}_{\alpha}$ and $h^{({\rm A},l)}_{\alpha}$.

Both $u$ and $q$ are defined at each sub-node in the same way as ordinary nodes, and the boundary condition is imposed for each sub-node. Either $u_{\alpha}$ or $u_{\beta}$ is unnecessary since $u_{\alpha}=u_{\beta}$; however, the same representation as the ordinary nodes simplifies the programming effort. In this case, the number of unknown variables increases by the number of double nodes, which equals the number of corners. Since one double node $i$ is replaced by the two sub-nodes $\alpha$ and $\beta$ for each corner, the number of node equations is also increased by the number of corners. Consequently, a set of CBIEs can be expressed by a matrix equation with a square coefficient matrix even when we apply double nodes.

Since the integral along the zero-sized elements is zero, the right-hand sides of Eq. (7) for $i=\alpha$ and $i=\beta$ are the same. Therefore, the coefficients $c_{\alpha}$ and $c_{\beta}$ must have the same value:

[TABLE]

In the case of linear elements, the singular integral of $h_{i}^{j}$ for the corners ($i,j=\alpha,\beta$) becomes zero since the $\boldsymbol{r}-\boldsymbol{r}_{i}$ for $i=\alpha$ or $\beta$ is perpendicular to $\boldsymbol{n}_{j}$ for $j=\alpha$ or $\beta$:

[TABLE]

Therefore, summarizing the discussions in the last paragraph in Sec. 2.1 and this result, we obtain the relations for the sub-nodes as

[TABLE]

2.3 Rank deficiency conditions in CBIEs

By using the double node technique, the boundary nodes including sub-nodes are defined at the ends of boundary elements, and the boundary conditions are defined at each node. In problems with a single medium, there are three common types of boundary conditions; Dirichlet condition, Neumann condition, and Robin condition. These conditions are given, respectively, as

[TABLE]

where the over-bars denote values imposed by the boundary conditions, $\kappa_{u}$ and $\kappa_{q}$ are the given constants determined by the problem considered, and the sets of the node numbers with corresponding boundary conditions are denoted by $I_{D}$, $I_{N}$, and $I_{R}$, respectively.

When we consider a multi-media problem, there exists a condition at the interface between the media. The interface condition is expressed by two continuous conditions for $u$ and $q$. The continuous conditions are given as

[TABLE]

where media numbers are denoted by (1) and (2); the nodes $j^{(1)}$ and $j^{(2)}$ are located at the same positions; $\kappa_{21}$ is determined by media constants of $\Omega^{(1)}$ and $\Omega^{(2)}$; and $I_{I}$ denotes the set of nodes with interface conditions.

In this section, we first present three simple examples when two sub-nodes of a double node at a corner point belong to $I_{D}$ or $I_{N}$, and discuss why the rank deficient problem arises. Then, the case in which the sub-nodes belong to $I_{R}$ or $I_{I}$ are presented.

The combinations of the boundary conditions where the two sub-nodes, $\alpha$ and $\beta$, belong to $I_{D}$ or $I_{N}$ are classified into the following three cases:

both include the Dirichlet conditions: $\alpha\in I_{D}$ and $\beta\in I_{D}$,

• 2.

both include the Neumann conditions: $\alpha\in I_{N}$ and $\beta\in I_{N}$,

• 3.

one includes the Dirichlet condition and the other includes the Neumann condition:

 $\alpha\in I_{D}$ and $\beta\in I_{N}$, or $\alpha\in I_{N}$ and $\beta\in I_{D}$.

2.3.1 Case of two Dirichlet conditions

Under these conditions ($\alpha\in I_{D}$ and $\beta\in I_{D}$), the sub-node equations of Eq. (8) for $i={\alpha}$ and ${\beta}$ can be arranged so that the terms including known values are moved to the right-hand side and the terms related to $\alpha$ and $\beta$ are moved out from the summation as

[TABLE]

where $I=I_{D}\cup I_{N}$, i.e., all nodes, and $I{\setminus\{\alpha,\beta\}}$ denotes a set of all node numbers except $\alpha$ and $\beta$; and the detailed descriptions of the third terms on the left-hand sides and the right-hand sides are written as

[TABLE]

Comparing the coefficients of the terms on the left-hand side of the two sub-node equations shown in Eqs. (26) and (27), we can evaluate whether the rank of the coefficient matrix is reduced or not. From Eqs. (20) and (21), the contributions $h_{i}^{j}$ and $g_{i}^{j}$ ($i=\alpha,\beta$) to the node $j$ are the same. The coefficients $c_{\alpha}$ and $c_{\beta}$ share the same value from Eq. (18); however, the multiplied terms, $u_{j}$, are different since $u_{j}$ is also multiplied by the Kronecker’s delta in Eq. (8). Therefore, we can examine rank reduction by examining which terms simply include the Kronecker’s delta. The Kronecker’s delta is not found in the first and the second terms on the left-hand sides in Eqs. (26) and (27). For the third terms, node $j$ does not include $\alpha$ and $\beta$; therefore, the $\delta_{\alpha}^{j}$ and $\delta_{\beta}^{j}$ are not included. They only appear in the right-hand sides, which are not related to rank reduction. Therefore, the left-hand sides of the two equations are identical, and the rank is always reduced by these two sub-node equations; i.e., the coefficient matrix becomes a singular matrix.

2.3.2 Case of two Neumann conditions

As in the previous sub-section, two sub-node equations of Eq. (8) in the case of $\alpha\in I_{N}$ and $\beta\in I_{N}$ can be arranged as

[TABLE]

Although the definitions of $a_{i}^{j}$ and $\overline{b}_{i}$ are different from Eqs. (28) and (29), the third terms of Eqs. (30) and (31) are identical since they do not include the Kronecker’s delta. In contrast, the first and second terms are different. By applying Eqs. (18) and (19), the coefficients associated with $u_{\alpha}$ and $u_{\beta}$ in Eq. (30) are $c_{\alpha}$ and 0, respectively; while those in Eq. (31) are 0 and $c_{\alpha}$, respectively. Therefore, the rank of the matrix that includes the sub-node equations is not reduced in the case of $\alpha\in I_{N}$ and $\beta\in I_{N}$. In addition, the right-hand sides are the same since they do not include the Kronecker’s delta. This means that the results of these two sub-node equations involve the relation $u_{\alpha}=u_{\beta}$.

2.3.3 Case of coupled Dirichlet and Neumann conditions

In the case of $\alpha\in I_{D}$ and $\beta\in I_{N}$, two sub-node equations are

[TABLE]

Since the second terms on the left-hand sides of the above equations are different, the rank of the coefficient matrix is not reduced in the case of $\alpha\in I_{D}$ and $\beta\in I_{N}$.

According to the examples in Secs. 2.3.2 and  2.3.3, we can understand that the two sub-node equations for $\alpha$ and $\beta$ are different when either $c_{\alpha}$ or $c_{\beta}$ is included in the coefficients of unknown variables in the two sub-node equations.

2.3.4 Case of Robin condition

We consider the case in which the Robin conditions are imposed on at least one of the two sub-nodes. Eliminating $u_{j}$ from Eq. (8) using Eq. (24), and arranging the equation, we obtain:

[TABLE]

When $i=\alpha\in I_{R}$, the term with the factor $c_{\alpha}$ for unknown $q_{\alpha}$ remains in the third terms of the left-hand side. This satisfies the condition described at the end of Sec. 2.3.3, which ensures that the equations of the two sub-nodes are independent and not identical. It does not depend on the type of boundary condition of the node $\beta$.

2.3.5 Case of interface condition

Figure 2 illustrates a configuration around a double node where two domains are in contact. The corner node labeled P is represented by the four sub-nodes at the corner; two sub-nodes labeled $\alpha^{(1)}$ and $\alpha^{(2)}$ do not belong to the interface boundary, and the other two sub-nodes labeled $\beta^{(1)}$ and $\beta^{(2)}$ belong to the interface boundary. There are eight quantities related to these four sub-nodes. Of the eight quantities, two quantities at the sub-nodes $\alpha^{(1)}$ and $\alpha^{(2)}$ are considered ordinary boundary conditions, and two variables at $\beta^{(1)}$ and $\beta^{(2)}$ can be eliminated by the interface conditions shown in Eq. (25), which, for $\beta^{(1)}$ and $\beta^{(2)}$, are rewritten as

[TABLE]

where the double suffixes, such as $u_{\beta^{(1)}}$, are denoted by $u_{{\beta(1)}}$ for simplicity. Consequently, there are four unknown variables and four equations with respect to the double node when we apply these conditions before constructing the simultaneous equations.

As mentioned in the previous sub-sections, the coefficients unrelated to the sub-nodes do not affect rank reduction; therefore, we only consider the nature of the sub-matrix composed of the coefficients for the four unknown variables.

In the case of $\alpha^{(1)}\in I_{D}$ and $\alpha^{(2)}\in I_{D}$, the given values are $\overline{u}_{{\alpha(1)}}$ and $\overline{u}_{{\alpha(2)}}$, and the unknown variables are $q_{{\alpha(1)}}$, $u_{{\beta(1)}}$, $q_{{\beta(1)}}$, $q_{{\alpha(2)}}$, $u_{{\beta(2)}}$, and $q_{{\beta(2)}}$. The independent variables are reduced by using Eq. (35) such as $\{q_{{\alpha(1)}},q_{{\beta(1)}},q_{{\alpha(2)}},u_{{\beta(1)}}\}$. The terms associated with these independent variables in the individual sub-node equations for $i=\alpha^{(1)}$, $\beta^{(1)}$, $\alpha^{(2)}$, and $\beta^{(2)}$, are written as

[TABLE]

By applying Eqs. (18), (19) and (20) to the sub-matrix to eliminate coefficients with the subscript $\beta^{(1)}$ and $\beta^{(2)}$, the sub-matrix is rewritten as

[TABLE]

Furthermore, by applying elementary matrix operations, the sub-matrix is transformed to

[TABLE]

where the second row is the result of replacement with the difference between the first and the second row of the matrix in Eq. (49), and the fourth row is obtained by similar operations. The non-zero components in both the second and the fourth row are only found in the fourth column; therefore, the rank of this sub-matrix is reduced by one.

In other cases, the product of the sub-matrices and the unknown vector after applying Eqs. (18), (19) and (20) are as follows. In the case of $\alpha^{(1)}\in I_{N}$ and $\alpha^{(2)}\in I_{N}$,

[TABLE]

In the case of $\alpha^{(1)}\in I_{D}$ and $\alpha^{(2)}\in N_{N}$,

[TABLE]

Since $\kappa_{21}>0$ generally, Eqs. (63) and (72) do not become singular.

Furthermore, in the case of $\alpha^{(1)}\in I_{R}$ or $\beta^{(2)}\in I_{R}$, the coefficient matrix becomes more complex than Eq. (63) or Eq. (72); therefore, the rank is also not reduced.

2.3.6 Case of internal interface condition

In the case where the node is located at a corner of two interface boundary elements as shown in Fig.​  3(a), there are four sub-nodes at the corner node. After applying Eqs. (18), (19) and (20), and eliminating the quantities related to $\alpha^{(2)}$ and $\beta^{(2)}$, the terms related to the sub-nodes in the four sub-node equations are expressed as

[TABLE]

By applying the elementary matrix operations to the coefficient matrix, the matrix is transformed to

[TABLE]

We can understand that the rank of this matrix is reduced by one since the determinant of a 2x2 sub-matrix composed of non-zero elements in the second and the fourth row in the above sub-matrix is zero. When we increase the number of the domains as shown in Fig.​  3(b), we can illustrate the singularity, although the details are not shown here. Therefore, the sub-node equations become singular when the node is located at an intersection of the interface boundaries which partitions the original single domain into two or more domains.

2.4 Summary of rank deficiency conditions for CBIEs

In this section, we analyze the rank deficiency conditions for the sub-node equations of a double node for CBIEs. The rank of the coefficient matrix is reduced in the following cases:

both sub-nodes of a double node are imposed by Dirichlet conditions,

• 2.

sub-nodes are located at an intersection of an interface boundary element and two boundary elements imposed by Dirichlet conditions,

• 3.

sub-nodes are located at intersections of interface boundary elements with different normal directions and they do not belong to any boundary with ordinary boundary conditions.

The original sub-node equations or the sub-matrix associated with the sub-nodes for these cases include a pair of Eqs. (26) and (27) for the first case, Eq. (44) for the second case, and Eq. (81) for the last case. These equations have common characteristics; all the equations include $q_{\alpha}$ and $q_{\beta}$ (or $q_{\alpha(1)}$ and $q_{\beta(1)}$ for the interface boundary) as the unknown variables. The coefficients associated with these unknowns are $g_{i}^{j}$ ($i,j\in\{\alpha,\beta\}$), which possess the characteristics shown in Eq. (20). The terms with $q_{\alpha}$ and $q_{\beta}$ are canceled because of this characteristic during the transformations of the equations, and the ranks have been reduced.

If the coefficients $g_{i}^{j}$ did not satisfy Eq. (20); in other words, if they included information on $\boldsymbol{n}_{i}$, these singularities could be avoided.

3 Regularization of coefficient matrix using HBIEs

As described in Sec. 2.4, the reason for a rank deficient matrix in CBIE is the coefficients $g_{i}^{j}$ do not carry information concerning the normal direction of the node $i$. To provide this information, we employ a hypersingular boundary integral equation (HBIE). The HBIE can be derived by taking a gradient of the CBIE shown in Eq. (2.1) with respect to $\boldsymbol{r}_{i}$:

[TABLE]

where $\nabla\!_{i}$ denotes the gradient with respect to $\boldsymbol{r}_{i}$. The gradient, $\nabla\!_{i}\nabla u^{*}$, has a stronger singularity than the CBIE, and this is known as hypersingularity. However, the integral of the hypersingular function can be regularized. In this study, we employ the regularization for the Helmholtz equation which uses the fundamental solution of Laplace’s equation [31].

3.1 Regularized gradient field

In this section, after presenting a regularization of the gradient at $\boldsymbol{r}_{i}$ according to the method shown in Ref. [31], a discretized form of the normal derivative, $q$, is derived.

The gradient is expressed by a 2x2 dyadic tensor $\overleftrightarrow{{\bf{C_{\mathit{i}}}}}$ and a vector $\boldsymbol{J}_{i}(u,q)$ as

[TABLE]

The second term can be evaluated by two types of boundary integrals as

[TABLE]

where $\boldsymbol{J}_{i}^{\gamma,\rm sing}(u,q)$ is associated to the boundary integral along the boundary elements that include the node $i\in\{\alpha,\beta\}$ and $\gamma\in\{{\rm A},{\rm B}\}$ (Fig.​  1); whereas $\boldsymbol{J}_{i}^{k,\rm reg}$ corresponds to integrals of the boundary elements not including the node $i$. Similar to the factor $c(\boldsymbol{r}_{i})$ in CBIEs, the singular integral included in Eq. (3) contributes to $\overleftrightarrow{{\bf{C_{\mathit{i}}}}}$ on the left-hand side of Eq. (88) and $\boldsymbol{J}_{i}^{\gamma,\rm sing}(u,q)$. The components of $\boldsymbol{J}_{i}^{\gamma,\rm sing}(u,q)$, $\boldsymbol{J}_{i}^{k,\rm reg}(u,q)$ and $\overleftrightarrow{{\bf{C_{\mathit{i}}}}}$ are written as

[TABLE]

where $u_{L}^{*}$ denotes the fundamental solution to the Laplace equation that satisfies Eq. (6), $q_{L}^{*}$ is its normal derivative with respect to $\boldsymbol{n}(\boldsymbol{r})$ at $\boldsymbol{r}$; and $\boldsymbol{\tau}_{\gamma}$, $\theta_{\gamma}$, and $L_{\gamma}$ are the tangential unit vectors along $\Gamma_{\gamma}$, the azimuth angle of $\Gamma_{\gamma}$ from $\boldsymbol{r}_{i}$, and the length of $\Gamma_{\gamma}$, respectively.

Similar to CBIEs, the non-singular integral in Eq. (91) is expressed by the shape function as

[TABLE]

which are evaluated using numerical integrals. By using the permutation matrix, the second term on the left-hand side in Eq. (89) is denoted by a form of linear combinations as

[TABLE]

The two singular integrals on the right-hand side of Eq. (3.1) can be analytically obtained by using Taylor series expansions around $\boldsymbol{r}_{i}$ for the two second-order derivatives. The first term on the right-hand side of Eq. (89) is also expressed as linear combinations with respect to the nodes neighboring the sub-nodes $\alpha$ and $\beta$ as

[TABLE]

By defining

[TABLE]

the vector $\boldsymbol{J}_{i}(u,q)$ is expressed as the following form with a linear combination:

[TABLE]

By substituting this equation into Eq. (88) and by multiplying both sides by $\overleftrightarrow{{\bf{C_{\mathit{i}}^{\rm{-1}}}}}$, the gradient field at $\boldsymbol{r}_{i}$ is expressed as

[TABLE]

3.2 Rank deficiency conditions in HBIEs

Let us consider the relations of $\boldsymbol{s}_{i}^{j}$ and $\boldsymbol{t}_{i}^{j}$ when $\boldsymbol{r}_{i}$ is a double node; i.e., $\boldsymbol{r}_{i}=\boldsymbol{r}_{\alpha}=\boldsymbol{r}_{\beta}$. If both $u_{j}$ and $q_{j}$ at all nodes are determined as known values, the gradient in Eq. (103) must be expressed uniquely. Since the dyadic tensor $\overleftrightarrow{{\bf{C_{\mathit{i}}^{\rm{-1}}}}}$ on the right-hand side of Eq. (103) is determined by a geometric configuration of the boundary elements connected to $\boldsymbol{r}_{i}$ as shown in Eq. (94), $\overleftrightarrow{{\bf{C_{\mathit{\alpha}}}}}=\overleftrightarrow{{\bf{C_{\mathit{\beta}}}}}$ is satisfied. Therefore, to uniquely determine the gradient, i.e., $\nabla u(\boldsymbol{r}_{\alpha})=\nabla u(\boldsymbol{r}_{\beta})$, the coefficient vectors for $\alpha$ and $\beta$ must be the same. Strictly speaking, there is an exception in $\boldsymbol{s}_{i}^{j}$ for $j=\alpha,\beta$. In this case, the condition $\boldsymbol{s}_{\alpha}^{j}\neq\boldsymbol{s}_{\beta}^{j}$ is permitted when $\boldsymbol{s}_{\alpha}^{\alpha}+\boldsymbol{s}_{\alpha}^{\beta}=\boldsymbol{s}_{\beta}^{\alpha}+\boldsymbol{s}_{\beta}^{\beta}$ is satisfied since $u_{\alpha}=u_{\beta}$. Therefore,

[TABLE]

Since $q_{i}=\boldsymbol{n}(\boldsymbol{r}_{i})\cdot\nabla u(\boldsymbol{r}_{i})$, by taking an inner product of $\boldsymbol{n}(\boldsymbol{r}_{i})$ to Eq. (103), the discretized equation for the HBIE is obtained as

[TABLE]

where

[TABLE]

When $\boldsymbol{n}_{\alpha}\neq\boldsymbol{n}_{\beta}$, we can derive the following conditions using Eqs. (105) and (106):

[TABLE]

These relations mean that the two sub-node equations of Eq. (107) for $i=\alpha$ and $i=\beta$ are different. Even in the case of $\boldsymbol{n}_{\alpha}=\boldsymbol{n}_{\beta}$, the two sub-node equations are also different since the coefficients for $u_{\alpha}$ of the two equations are always different because of Eq. (104) and the coefficients for $q_{\alpha}$ are always different because of the nature of $\delta_{i}^{j}$.

The aforementioned difference in the two equations is not a sufficient condition to not reduce the rank since we have not proved that the determinant of the sub-matrix associated with the sub-nodes is not equal to zero for any problem. However, the HBIE is highly applicable to the problem in which the node equations constructed by the CBIE become singular.

4 Numerical results and discussions

To demonstrate the validity of the method in regularizing the coefficient matrix using HBIEs, we analyzed simple models as shown in Fig.​  4. The actual models correspond to electromagnetic wave propagation problems in a waveguide where the wall parallel to the propagation direction is made of metal. When we consider $u$ as the $z$-component of the electric field, $u=0$ at the side-walls of the waveguide since the component of the electric field parallel to the metal is zero. In the case of model-(A), the metal walls are located at $\Gamma_{D}$, which is called the TE10 mode for a waveguide with rectangular cross-section and width $W$. In the case of model-(B), virtual walls $\Gamma_{N}$ with $q=0$ are located at $x=\pm W/2$, which are equivalent to the placement of physical metal walls with $u=0$ at $x=\pm W$, which is called the TE20 mode. These incident conditions are given as functions of $x$ as

[TABLE]

In these modes, a propagation field to $\pm y$-directions are proportional to $\exp({\mp{\rm j}k_{y}})$, respectively, where

[TABLE]

By using this characteristic, the incident boundary condition on $\Gamma^{\rm inc}$ can be rewritten as a Robin condition:

[TABLE]

On the boundary $\Gamma^{\rm T}$ at $y=L$, we analyzed the following three conditions for termination:

[TABLE]

where $R$ denotes the reflection coefficient of the electric field on $\Gamma_{T}$; and the non-reflective condition in (c) is equivalent to the incident condition with $\overline{u^{\rm inc}}=0$ in Eq. (116). Among the above three, termination-(a) (short type) has two double nodes at $(\pm W/2,L)$, where both sub-nodes have two Dirichlet conditions; therefore, the rank of the coefficient matrix of CBIEs will be reduced as mentioned in Sec. 2.3.1. The exact solutions to these models are given by

[TABLE]

To simplify, we used $\overline{u^{\rm inc}_{0}}=1$, $W=2$, $L=2$, and $\lambda=2\pi/k=1$; under this condition, $\max|u|=2$ for the cases (a) and (b), and $\max|u|=1$ for the case (c), $\lambda_{x}=2\pi/k_{x}=4$ and $\lambda_{y}=2\pi/k_{y}\simeq 1.033$. As examples of $\hat{u}(x,y)$, the real part of $\hat{u}(x,y)$ for the models (A-b) and (B-b) are shown in Fig.​  5.

We analyzed three types of the following simultaneous equations:

1. (i)

CBIE: BIEs for all nodes including sub-nodes are obtained from the CBIE, 2. (ii)

partial-HBIE: BIEs for only sub-nodes related to double nodes are obtained from the HBIE, and BIEs for the other nodes are obtained from the CBIE. 3. (iii)

full-HBIE: BIEs for all nodes including sub-nodes are obtained from the HBIE,

To solve the complex-valued simultaneous equations, we employed two subroutines based on an LU decomposition provided in Lapack [33], in which the subroutine names are ‘zgetrf’ and ‘zgetrs.’ After determining both $u$ and $q$ at all boundary nodes by solving the simultaneous equations, the internal field $u(\boldsymbol{r}_{i\,^{\prime}})$ can be evaluated by

[TABLE]

which is derived from the discretized CBIE shown in Eq. (2.1). This process is the same for all types of simultaneous equations.

Figure 6 presents the distributions of errors based on the three types of simultaneous equations (i), (ii), and (iii) for the model-(A-b) and model-(B-b). In these models with termination-(b) (open type), the coefficient matrices, even in the case of the CBIE shown in Fig.​  6(i), are not singular as mentioned in Sec. 2.4.

Based on the comparison between the sub-figures (A-b-i) and (A-b-iii) or between (B-b-i) and (B-b-iii) in Fig.​  6, we can observe that the error of the full-HBIE is several times larger than that the CBIE. This difference can be explained by two reasons. The first reason is the difference in singularity to evaluate the coefficients. The coefficients $w_{i}^{j}$ and $v_{i}^{j}$ in the HBIE are evaluated from $\boldsymbol{s}_{i}^{j}$ and $\boldsymbol{t}_{i}^{j}$ with a multiplication of dyadic $\overleftrightarrow{{\bf{C_{\mathit{i}}}}}$ as shown in Eqs. (108) and (109), respectively. The vectors $\boldsymbol{s}_{i}^{j}$ and $\boldsymbol{t}_{i}^{j}$ are evaluated by the boundary integral shown in Eqs. (96), (97) and (98). Their integrands include the first or second order derivatives of the fundamental solution. The strongest singularity in the HBIE is $O(r^{-2})$, while the strongest singularity in the CBIE is $O(r^{-1})$. This error emerges significantly in the contributions between two nodes with short distances. The second reason is the multiplication by $\overleftrightarrow{{\bf{C_{\mathit{i}}}}}$. In the worst case, the errors of $v_{i}^{j}$ and $w_{i}^{j}$ are multiplied by the maximum norm,       $\overleftrightarrow{{\bf{C_{\mathit{i}}^{\rm{-1}}}}}$      , and the errors of $\boldsymbol{s}_{i}^{j}$ and $\boldsymbol{t}_{i}^{j}$, respectively. This amplification affects all coefficients regardless of the distances between nodes. According to Ref. [31], \,\hbox{ \vrule height=3.66875pt,depth=4.0pt\,\vrule height=3.66875pt,depth=4.0pt\,\hbox{\overleftrightarrow{{\bf{C_{\mathit{i}}^{\rm{-1}}}}}}\,\vrule height=3.66875pt,depth=4.0pt\,\vrule height=3.66875pt,depth=4.0pt\,}\,\leq 4\pi/(\pi-2)\simeq 11 in the case of $\Delta\theta=\pi/2$ and $L_{{\rm A}}=L_{{\rm B}}$.

In contrast, the error in the case of the partial-HBIE shown in Fig.​  6(A-b-ii) is similar to the case of the CBIE in (A-b-i); and the relation between (B-b-ii) and (B-b-i) is also similar. The number of unknowns in the analyses presented in Fig.​  6 was 164 including four double nodes; i.e., the number of HBIEs was only 8 and the number of CBIEs was 156 in the partial-HBIE. Since the number of HBIEs with a larger error is sufficiently smaller than that of CBIEs, the total error of partial-HBIEs does not increase so much as that of CBIEs.

Based on the comparisons between model-(A) and model-(B) in Fig.​  6, we can observe that the error of model-(B) is larger than that of model-(A) for each equation type. This reason is the same as the first reason for the difference between CBIEs and full-HBIEs. The coefficients $g_{i\,^{\prime}}^{j}$ and $h_{i\,^{\prime}}^{j}$ in Eq. (119) are results of the boundary integrals in which integrands include $u^{*}$ and $q^{*}$, respectively. Since the singularity of $q^{*}$ is stronger than that of $u^{*}$, the contribution of $h_{i\,^{\prime}}^{j}$ to $u(\boldsymbol{r}_{i\,^{\prime}})$ is larger than that of $g_{i\,^{\prime}}^{j}$, especially in the case where the distance between the field point and boundary elements is smaller than a wavelength. The error caused by $h_{i\,^{\prime}}^{j}$ is also larger than $g_{i\,^{\prime}}^{j}$. In model-(A), there are many nodes with $\overline{u}=0$. Therefore, the larger error caused by $h_{i\,^{\prime}}^{j}$ does not appear, and the smaller error caused by $g_{i\,^{\prime}}^{j}$ becomes dominant. In contrast, in model-(B), since there are no nodes with $\overline{u}=0$, the larger error caused by $h_{i\,^{\prime}}^{j}$ remains.

We analyzed the errors of different boundary element sizes. The error is evaluated based on average sampling points as follows:

[TABLE]

where $N_{s}$ is the number of sampling points that are intersections of the grid in Fig.​  5 ($N_{s}=41^{2}$), and $N_{s}$ is unchanged for all results regardless of the element size. The errors are not normalized since the averaged intensities of the exact solutions, $\langle|\hat{u}|\rangle$, have almost the same order of magnitude; $\simeq$0.8 for the termination types (a) and (b), and $\simeq$0.6 for (c). Figure 7 illustrates the dependence between the error and the number of elements in a wavelength, $N_{\lambda}=\lambda/\Delta_{\rm Ele}$ where $\Delta_{\rm Ele}$ denotes the boundary element size. In the case of (a), there are no plots for the short termination type (i) since the simultaneous equations become singular. In all cases of (a), (b) and (c), the errors decrease as $N_{\lambda}$ increases with a decay proportional to $1/N_{\lambda}^{2}$, except for the points of $N_{\lambda}\gtrsim 200$ in (b). This property is reasonable if $h_{i}^{j}$, $g_{i}^{j}$, $v_{i}^{j}$, and $w_{i}^{j}$ have accuracies of $O(\Delta_{\rm Ele})$ in the case of the linear element, and the truncated error is proportional to $O(\Delta_{\rm Ele}^{2})$. However, there is an error which does not show this characteristic [31]; the order of the error of $\boldsymbol{s}_{i}^{j,{\rm reg}}$ for a short distance between the nodes $i$ and $j$ obeys $O(\Delta_{\rm Ele})$. This error arises when $N_{\lambda}\gtrsim 200$. In the cases (b) and (c), the error of the full-HBIE (iii) is several times larger than the CBIE (i), and that of the partial-HBIE (ii) is almost the same as (i), which is similar to the result previously shown in Fig.​  6.

Figure 8 presents the computational time, which does not include the CPU time of the file input and output processes. The computation consists of several major steps; the computation of the components of the coefficient matrix $a_{i}^{j}$ (through $h_{i}^{j}$ and $g_{i}^{j}$ in Eqs. (17) and (14), respectively; or $v_{i}^{j}$ and $w_{i}^{j}$ in Eqs. (108) and (109), respectively); solving the matrix equation; and the evaluation of the internal field using Eq. (119); for which individual computational costs are proportional to $N_{\lambda}^{2}$, $N_{\lambda}^{3}$, and $N_{\lambda}$, respectively. Higher-order terms appear with increasing $N_{\lambda}$. In the case where $N_{\lambda}\lesssim 5$, most of the computation time is exhausted in minor common steps such as initializing tables for the Hankel functions. The computational cost for the full-HBIE (iii) is larger than that for the others. This is because of the difference between the evaluation time of $v_{i}^{j}$ and $w_{i}^{j}$ for the HBIE and that of $h_{i}^{j}$ and $g_{i}^{j}$ for the CBIE. In the HBIE, the cost of evaluating $v_{i}^{j}$ and $w_{i}^{j}$ is mainly exhausted in the numerical integrals of non-singular elements for the three vectors in Eqs. (96), (97) and (98). In contrast, two scalar integrals in Eqs. (16) and (13) are dominant in the CBIE. The cost of evaluating a coefficient with a vector is twice that of a scalar in two-dimensional problems, and the number of components in the HBIE is 3/2 times greater than the CBIE. Moreover, the operator $\boldsymbol{n}\cdot\nabla_{i}\nabla$ in Eqs. (96) and (97) has two vector components $\boldsymbol{e}_{r}\boldsymbol{e}_{r}\cdot\boldsymbol{n}$ and $\boldsymbol{n}$. Because some of the terms have common factors, the sum of costs was reduced from these estimations; however, the cost of evaluating the coefficient matrix component in the HBIE is almost four times larger than that in the CBIE. Even when the simultaneous equations are singular, it can be solved as a minimal-norm solution of underdetermined equations by using a solver based on a singular value decomposition (SVD). The details are not included in this paper because the authors do not understand whether the minimal-norm solution is always correct or not. By limiting the examples shown here, the accuracy of the CBIE by using a solver based on SVD was almost the same as in the case of the partial-HBIE. However, the computational time of SVD was much larger than in the case of LU decomposition; e.g., the time to solve 6,900 s for the CBIE using SVD called ‘zgelss’ in Lapack, 220 s for the partial-HBIE using the LU decomposition in the case of $N_{\lambda}=1,000$.

The above results can be summarized as follows. First, as mentioned in Sec. 2.3, the set of simultaneous equations constructed by the CBIE for all nodes becomes singular when both boundary conditions of the double node are imposed by Dirichlet conditions. In contrast, when the node equations are constructed by the HBIE for all or a part of the nodes, the set of simultaneous equations does not become singular. Next, the accuracy of the HBIE is unfortunately worse than the CBIE when the set of equations is regular; however, in the case of the partial-HBIE where only the equations of the sub-nodes belonging to the double nodes are given by the HBIE and the others are given by the CBIE, the reduction in accuracy is negligibly small. Finally, more computational cost is required to compute the component of the coefficient matrix by the HBIE than the CBIE. The rise in computational cost can be suppressed by applying the HBIE to sub-nodes only. Therefore, we can conclude that the replacement of the CBIE by the HBIE only for the sub-nodes (partial-HBIE) is the best solution from the viewpoint of singularity, accuracy, and computational cost.

To demonstrate the applicability of the HBIE in the cases with interface boundaries, we evaluated the models shown in Fig.​  9. In these models, the original model shown in Fig.​  4(A) or (B) is partitioned into two or four sub-domains by interface boundaries, and the exact solutions are the same as Eqs. (117) and (118). The number of sub-nodes for each double node is four at two intersections in the model-(A:

$\,|\,$

) and at four corners in (A:

$\times$

) and (B:

$\times$

). In addition, at the intersections of the four interface boundaries in (A:

$\times$

) and (B:

$\times$

), the number of sub-nodes for each double node is eight. The combinations of four types of boundary conditions (Dirichlet, Neumann, Robin, and interface conditions), can be examined by these models.

Figure 10 presents comparisons of the errors. As predicted in Sec. 2.3.5, when the simultaneous equations are constructed by the CBIE, the set of equations for each model including the interface boundary conditions is always singular. Similar to the above discussions on the single region problems, the error in the partial-HBIE (ii) is less than that in the full-HBIE (iii) in the case of multi-regions problems. Based on the comparison between model-(A) and either (A:

$\,|\,$

) or (A:

$\times$

), we can observe that the errors in the multi-media problems with the interface boundaries (models-(A:

$\,|\,$

) or (A:

$\times$

)) are larger than that in the single region problem (model-(A)). One of the reasons could be the same reason for which the model, including the boundary with $\overline{q}=0$ such as model-(B), has a larger error than that with $\overline{u}=0$ such as model-(A), and this is discussed in the description of Fig.​  6. The internal field is evaluated by Eq. (119) as the boundary integral where the boundary encloses the domain considered. In the case of the multi-regions partitioned by the interface boundaries, the boundary enclosing a single region must include the continuous boundary where $u\neq 0$. Therefore, the error in $h_{i^{\prime}}^{j}$ contributing from the interface boundary in Eq. (119) is added to the total error in the multi-region problems; whereas it is not added from the boundary with $\overline{u}=0$ in model-(A). The other reason is a difference in the distances between the internal points and their nearest boundary; an average of distances in the multi-region problem is shorter than that of the single region problem. Since the contributions $h_{i^{\prime}}^{j}$ and $g_{i^{\prime}}^{j}$ increase with decreasing distance, the error in the multi-region problems is larger than that in the single-region problem. The difference in error between the multi- and single-region problem can also be found in the comparison between model-(B) and (B:

$\times$

).

Consequently, we can demonstrate that the formulation based on the HBIE is applicable without rank deficiency even in the cases involving the interface boundaries, which is similar to the corner nodes in the single region problems.

One may question whether the partial-HBIE method can avoid spurious solutions of an external problem shown in Sec. 1 since the partial-HBIE method uses both CBIEs and HBIEs like the Burton-Miller method does to avoid spurious solutions. In the Burton-Miller method, a linear combination of CBIEs and HBIEs with an appropriate combination factor is used to prevent spurious solutions from arising when some of the sub-matrices become singular. Whereas, in the partial-HBIE method, the matrix equation consists of two sets of equations, CBIEs and HBIEs, without any modifications. Each set of equations has sub-matrices which may potentially produce spurious solutions. The partial-HBIE method, therefore, cannot avoid spurious solutions. To avoid them, we should use other methods, such as the Burton-Miller method [23], CHIEF [21], or a virtual boundary method [34, 35] which divides the external region into multiple regions by virtual boundaries to stop the external region from surrounding the internal region. If we use a virtual boundary method, the issue of multiple-duplicated nodes has to be solved; that is, however, not very difficult when using the method proposed in this paper.

5 Conclusion

The method of using double nodes at corners is a useful approach to uniquely define the normal direction. However, a set of simultaneous equations in CBIE formulation produces rank deficient problems in the following cases: both sub-nodes belonging to any double node are imposed by Dirichlet conditions; an intersection of the interface boundary located between different media is not connected to the boundary imposed by ordinary boundary conditions; and an interface boundary is connected to the two boundaries imposed by Dirichlet conditions. This means that the applicable problem that uses the double nodes are limited in the CBIE formulation.

In contrast, when the coefficient matrix is constructed by HBIEs, the rank is not reduced for any combination of boundary conditions, including interface conditions. However, the contribution coefficients between nodes in HBIEs are less accurate than those in CBIEs for the problem without rank deficiency because of two reasons; a HBIE exhibits a stronger singularity of the integrand than a CBIE, and most of the coefficients are multiplied by the dyadic tensor with a large norm. Furthermore, the computational cost of evaluating the coefficients of HBIEs is higher than that of CBIEs.

To address the rank deficiency problem in CBIEs and the drawbacks in HBIEs, the coupling approach presented in this paper called partial-HBIE is the best choice. In the partial-HBIE, most node equations are constructed by CBIEs, and only the sub-node equations related to corners are constructed by HBIEs.

The method that uses HBIEs demonstrates the following advantages compared to other methods: it does not require any additional local relation between nodal points around double nodes, any extra boundary integral equation, and it does not require a least-square method, which can be computationally time-consuming. Furthermore, the partial-HBIE can be applied by only switching the sub-node equation for the double node from a CBIE to a HBIE; therefore, we can be relieved from the efforts involved in preparing input data and complex coding.

Acknowledgments

This work was supported by JSPS KAKENHI Grant Number 18K04158.