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

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.

## Authors

• 1 publication
• 1 publication
• 1 publication
• 1 publication
04/28/2020

### Boundary element methods for Helmholtz problems with weakly imposed boundary conditions

We consider boundary element methods where the Calderón projector is use...
01/16/2020

### About Three Dimensional Double-Sided Dirichlet and Neumann Boundary Value Problems for the Laplacian

The orthogonality of Hilbert spaces whose elements can be represented as...
05/05/2021

### Finite Element Methods for Isotropic Isaacs Equations with Viscosity and Strong Dirichlet Boundary Conditions

We study monotone P1 finite element methods on unstructured meshes for f...
01/19/2022

### Models for information propagation on graphs

In this work we propose and unify classes of different models for inform...
07/06/2020

### Corrected Trapezoidal Rules for Boundary Integral Equations in Three Dimensions

The manuscript describes a quadrature rule that is designed for the high...
10/16/2019

### Computational homogenization of time-harmonic Maxwell's equations

In this paper we consider a numerical homogenization technique for curl-...
07/03/2020

### Analytic solution of system of singular nonlinear differential equations with Neumann-Robin boundary conditions arising in astrophysics

In this paper, we propose a new approach for the approximate analytic so...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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, , along the boundary are continuous; however, the normal derivatives, , 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, , 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 and 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 Manolis:1986 ; Parreira:1988 ; Olukoko:1993 ; Blazquez:1994 ; Huesmann:1994 ; Paris:1995 ; Blazquez:1998 , fluid flow problems Patterson:1982 ; Dyka:1989 , and acoustic problems Silva:1993 . Although the accuracy between the non-conforming element and the conforming element were compared by Manolis and Banerjee Manolis:1986 , and Parreira Parreira:1988 , they arrived at different conclusions.

The second approach includes a double node technique Brebbia:1984:Sec_5_2 or a multiple node technique for three-dimensional problems. In the double node technique, two normal derivatives, and , and a field, , are defined at the corner node, where and 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 Walker:1989 ; Yan:1994 ; Kassab:1994 ; Gao:2000 . 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 Mitra:1987 ; Mitra:1993 ; Subia:1995 ; Deng:2013 ; Zheng:2018 . Mitra and Ingber Mitra:1987 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” Mitra:1993 , 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 Mitra:1993 ; Subia:1995 . 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 Subia:1995 . The method that uses the extra node equations was extended to the problems of interface boundaries at which two or more domains are connected Deng:2013 ; Zheng:2018 . 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)’ Schenck:1968 . Chen et al. IL_Chen:2001 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 Burton-Miller:1971 ; Benthien:1997 ; Diwan:2013 ; Langrenne:2015 . Burton and Miller Burton-Miller:1971 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 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 Benthien:1997 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 Portela:1992 ; JT_Chen:1994 ; Chyuan:2003 ; Lu:2010 or at both surfaces of a thin metal with zero-thickness in an electromagnetic problem Chyuan:2003 . 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 Tomioka:2010 ; 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 , where is an imaginary unit and is the angular frequency:

 ∇2u(r)+k2u(r)=0,r∈Ω, (1)

where indicates the wave number, which is a ratio of to the wave velocity, and represents the spatial domain considered. A fundamental solution , which represents the contribution to a field point from a unit source placed at a source point in free space satisfies

 ∇2u∗(r,ri)+k2u∗(r,ri)=−δ(r−ri), (2)

where the differential operator operates only on , but not on . This equation can be solved analytically. In two-dimensional problems, the outward propagating wave that obeys Eq. (2) is

 u∗(r,ri)=14jH(2)0(kr),r=|r−ri|, (3)

where the function 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),

 c(ri) u(ri)=∮Γ[u∗(r,ri)∇u(r)⋅n−u(r)∇u∗(r,ri)⋅n]dΓ,

where denotes a boundary surrounding , is the position of the source point on the boundary, is the position of the field point,

is the outward-pointing normal unit vector; the contour integral is evaluated as a Cauchy principal value, and

is the result of the following evaluation of Dirac’s delta function:

 c(ri)u(ri)≜∫Ωu(r)δ(r−ri)dΩ=∫Ωδ(r−ri)dΩu(ri). (5)

The coefficient depends on the shape of the boundary at the field point . Because of the nature of Dirac’s delta function, when is located inside and outside the domain, evaluates to 1 and 0, respectively. In the case where is located on , is equal to the ratio of the interior angle to the whole angle; e.g., for 2-dimensional problems. In addition, when we introduce a fundamental solution to a Laplace equation, , which satisfies

 ∇2u∗L(r,ri)=−δ(r−ri), (6)

the coefficient can be expressed by a boundary integral:

 c(ri)=∫Ωδ(r−ri)dΩ=−∮Γ∇u∗L⋅ndΓ, (7)

which is called an equipotential condition Brebbia:1980 .

When is located on and the boundary is partitioned into a number of boundary elements, Eq. (2.1) can be written as a node equation for the node as a discrete algebraic expression:

 ∑j∈I(ciδji+hji)uj−∑j∈Igjiqj=0, (8)

where both the superscripts and the subscripts and denote the node numbers and not the element numbers; denotes a set of the node numbers, which includes members; denotes the Kronecker’s delta; and . The factors and are results of boundary integrals in Eq. (2.1) as shown below.

In Eq. (8), there are two variables, and , at . One of them is specified by a boundary condition as a boundary value which is denoted by , and the other is an unknown variable denoted ; therefore, the number of variables is equal to that of the nodes, . Since the node can be placed at every boundary node, we can obtain equations of Eq. (8) as

 ∑j∈Iajixj=¯¯bi,¯¯bi=∑j∈Ibji¯¯¯xj for each i∈I, (9)

where both and are either or when simple boundary conditions are specified. The matrix composed of is called a coefficient matrix in the following discussions.

The coefficient represents the contribution to from the nodal quantity at the node , and represents the contribution to from

. These coefficients are evaluated by boundary integrals, which depend on a method that interpolates

and along the boundary element. To examine the nature of and , we present evaluations using a shape function to interpolate them. Let us consider the discretization of the boundary integral of . By denoting the -th boundary element as and a local node number as , on is represented by a linear combination of shape functions and at the boundary nodes:

 q(r)=Nl∑l=1ϕ(k,l)(r)q(k,l) on\ Γk, (10)

where denotes the number of nodes in a boundary element; e.g., for the linear element, for the second-order element. A global node number can be mapped from the local node numbers by a permutation matrix :

 q(k,l)=∑jmj(k,l)qj, (11)

where has the value 1 if and are associated, and 0 otherwise. Although the notation of is used here, it provides a symbolic meaning, and it is treated as a mapping function such as in actual coding to avoid summation procedures with respect to . By using these definitions, the boundary integral is discretized as

 ∮Γu∗(r,ri)q(r)dΓ =∑k∫Γku∗(r,ri)(∑lϕ(k,l)(r)(∑jmj(k,l)qj))dΓ =∑j(∑k∑lmj(k,l)g(k,l)i)qj=∑jgjiqj, (12)

where and are defined as

 g(k,l)i=∫Γku∗(r,ri)ϕ(k,l)(r)dΓ, (13) gji=∑k∑lmj(k,l)g(k,l)i. (14)

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 , which is called a singular integral. Similarly, the other integral in Eq. (2.1) is evaluated as

 ∮Γq∗(r,ri)u(r)dΓ=∑jhjiuj, (15) h(k,l)i=∫Γkq∗(r,ri)ϕ(k,l)(r)dΓ, (16) hji=∑k∑lmj(k,l)h(k,l)i, (17)

where , which is not the derivative at ; i.e., does not depend on .

To discuss the rank of the matrix equation shown in Eq. (9), the dependencies of both and on which is the normal unit vector at are important. First, is independent of since Eqs. (13) and (14) do not include regardless of or . In contrast, the case of is different from that of . Although in Eqs. (16) and (17) looks independent of even when is expressed using , there is an exception at in which becomes . The dependency of this exceptional case is presented in the next section.

### 2.2 Double node technique

In the node equation for node shown in Eq. (8), both and 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 is located at a corner, the normal derivative 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 is represented by two sub-nodes at the same position; . The node at is connected to both and a zero-sized element between and , and vice versa. The normal direction of the sub-nodes and can be determined from the directions normal to and , respectively; therefore, the direction of the derivatives of and 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 does not belong to , contributions of and , which are evaluated by analytical integrals as singular integrals for accurate evaluation, are similar to and .

Both and are defined at each sub-node in the same way as ordinary nodes, and the boundary condition is imposed for each sub-node. Either or is unnecessary since ; 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 is replaced by the two sub-nodes and 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 and are the same. Therefore, the coefficients and must have the same value:

 cα=cβ. (18)

In the case of linear elements, the singular integral of for the corners () becomes zero since the for or is perpendicular to for or :

 hαα=hαβ=hβα=hββ=0 for linear elements. (19)

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

 gjα=gjβ for any j, (20) hjα=hjβ for any j. (21)

### 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

 uj=¯¯¯uj for j∈ID, (22) qj=¯¯¯qj for j∈IN, (23) κuuj+κqqj=¯¯¯¯ψj for j∈IR, (24)

where the over-bars denote values imposed by the boundary conditions, and are the given constants determined by the problem considered, and the sets of the node numbers with corresponding boundary conditions are denoted by , , and , 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 and . The continuous conditions are given as

 uj(2)=uj(1),qj(2)=−κ21qj(1) for j(1),j(2)∈II, (25)

where media numbers are denoted by (1) and (2); the nodes and are located at the same positions; is determined by media constants of and ; and 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 or , and discuss why the rank deficient problem arises. Then, the case in which the sub-nodes belong to or are presented.

The combinations of the boundary conditions where the two sub-nodes, and , belong to or

are classified into the following three cases:

• both include the Dirichlet conditions: and ,

• both include the Neumann conditions: and ,

• one includes the Dirichlet condition and the other includes the Neumann condition:
and , or and .

#### 2.3.1 Case of two Dirichlet conditions

Under these conditions ( and ), the sub-node equations of Eq. (8) for and can be arranged so that the terms including known values are moved to the right-hand side and the terms related to and are moved out from the summation as

 (26) (27)

where , i.e., all nodes, and denotes a set of all node numbers except and ; and the detailed descriptions of the third terms on the left-hand sides and the right-hand sides are written as

 ∑j∈I∖{α,β}ajixj=∑j∈INhjiuj−∑j∈ID,j≠α,βgjiqj, (28) ¯¯bi=−∑j∈ID(ciδji+hji)¯¯¯uj+∑j∈INgji¯¯¯qj. (29)

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 and () to the node are the same. The coefficients and share the same value from Eq. (18); however, the multiplied terms, , are different since 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 does not include and ; therefore, the and 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 and can be arranged as

 (cα+hαα)uα+hβαuβ +∑j∈I∖{α,β}ajαxj=¯¯bα, (30) hαβuα+(cβ+hββ)uβ +∑j∈I∖{α,β}ajβxj=¯¯bβ. (31)

Although the definitions of and 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 and in Eq. (30) are and 0, respectively; while those in Eq. (31) are 0 and , respectively. Therefore, the rank of the matrix that includes the sub-node equations is not reduced in the case of and . 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 .

#### 2.3.3 Case of coupled Dirichlet and Neumann conditions

In the case of and , two sub-node equations are

 −gααqα +hβαuβ+∑j∈I∖{α,β}ajαxj=¯¯bα, (32) −gαβqα +(cβ+hββ)uβ+∑j∈I∖{α,β}ajβxj=¯bβ. (33)

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 and .

According to the examples in Secs. 2.3.2 and  2.3.3, we can understand that the two sub-node equations for and are different when either or 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 from Eq. (8) using Eq. (24), and arranging the equation, we obtain:

 ∑j∈IN(ciδji+hji)uj−∑j∈IDgjiqj−∑j∈IR(κqκu(ciδji+hji)+gji)qj=¯¯bi. (34)

When , the term with the factor for unknown 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 .

#### 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 and do not belong to the interface boundary, and the other two sub-nodes labeled and 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 and are considered ordinary boundary conditions, and two variables at and can be eliminated by the interface conditions shown in Eq. (25), which, for and , are rewritten as

 uβ(2)=uβ(1),qβ(2)=−κ21qβ(1), (35)

where the double suffixes, such as , are denoted by 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 and , the given values are and , and the unknown variables are , , , , , and . The independent variables are reduced by using Eq. (35) such as . The terms associated with these independent variables in the individual sub-node equations for , , , and , are written as

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−gα(1)α(1)−gβ(1)α(1)0hβ(1)α(1)−gα(1)β(1)−gβ(1)β(1)0cβ(1)+hβ(1)β(1)0κ21gβ(2)α(2)−gα(2)α(2)hβ(1)α(2)0κ21gβ(2)β(2)−gα(2)β(2)cβ(2)+hβ(1)β(2)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜⎝qα(1)qβ(1)qα(2)uβ(1)⎞⎟ ⎟ ⎟ ⎟ ⎟⎠. (36)

By applying Eqs. (18), (19) and (20) to the sub-matrix to eliminate coefficients with the subscript and , the sub-matrix is rewritten as

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−gα(1)α(1)−gβ(1)α(1)00−gα(1)α(1)−gβ(1)α(1)0cα(1)0κ21gβ(2)α(2)−gα(2)α(2)00κ21gβ(2)α(2)−gα(2)α(2)cα(2)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (37)

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

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−gα(1)α(1)−gβ(1)α(1)00000cα(1)0κ21gβ(2)α(2)−gα(2)α(2)0000cα(2)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (38)

where the second row is the result of replacement with the difference between the first and the second row of the matrix in Eq. (37), 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 and ,

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝cα(1)0−gβ(1)α(1)00cα(1)−gβ(1)α(1)000κ21gβ(2)α(2)cα(2)0cα(2)κ21gβ(2)α(2)0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜⎝uα(1)uβ(1)qβ(1)uα(2)⎞⎟ ⎟ ⎟ ⎟ ⎟⎠. (39)

In the case of and ,

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−gα(1)α(1)−gβ(1)α(1)00−gα(1)α(1)−gβ(1)α(1)cα(1)00κ21gβ(2)α(2)0cα(2)0κ21gβ(2)α(2)cα(2)0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜⎝qα(1)qβ(1)uβ(1)uα(2)⎞⎟ ⎟ ⎟ ⎟ ⎟⎠. (40)

Since generally, Eqs. (39) and (40) do not become singular.

Furthermore, in the case of or , the coefficient matrix becomes more complex than Eq. (39) or Eq. (40); 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 and , the terms related to the sub-nodes in the four sub-node equations are expressed as

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−gα(1)α(1)−gβ(1)α(1)cα(1)0−gα(1)α(1)−gβ(1)α(1)0cα(1)κ21gα(2)α(2)κ21gβ(2)α(2)cα(2)0κ21gα(2)α(2)κ21gβ(2)α(2)0cα(2)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜⎝qα(1)qβ(1)uα(1)uβ(1)⎞⎟ ⎟ ⎟ ⎟ ⎟⎠. (41)

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

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−gα(1)α(1)−gβ(1)α(1)cα(1)000cα(1)−cα(1)κ21gα(2)α(2)κ21gβ(2)α(2)cα(2)000cα(2)−cα(2)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (42)

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,

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

• 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. (36) for the second case, and Eq. (41) for the last case. These equations have common characteristics; all the equations include and (or and for the interface boundary) as the unknown variables. The coefficients associated with these unknowns are (), which possess the characteristics shown in Eq. (20). The terms with and are canceled because of this characteristic during the transformations of the equations, and the ranks have been reduced.

If the coefficients did not satisfy Eq. (20); in other words, if they included information on , 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 do not carry information concerning the normal direction of the node . 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 :

 ∇i[c(ri)u(ri)]=∮Γ [(∇iu∗(r,ri))(∇u(r))⋅n −u(r)(∇i∇u∗(r,ri))⋅n]dΓ, (43)

where denotes the gradient with respect to . The gradient, , 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 Tomioka:2010 .

In this section, after presenting a regularization of the gradient at according to the method shown in Ref. Tomioka:2010 , a discretized form of the normal derivative, , is derived.

and a vector as

 ←→Ci⋅∇u(ri)+Ji(u,q)=0. (44)

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

 Ji(u,q)=∑γ=A,BJγ,singi(u,q)+∑k≠A,BJk,regi(u,q), (45)

where is associated to the boundary integral along the boundary elements that include the node and (Fig.  1); whereas corresponds to integrals of the boundary elements not including the node . Similar to the factor in CBIEs, the singular integral included in Eq. (3) contributes to on the left-hand side of Eq. (44) and . The components of , and are written as

 Jγ,singi(u,q)=−∫Lγ01r(∂u∗∂r−∂u∗L∂r)dr⋅u(ri)nγ +∫Lγ0r∂u∗∂rdr⎛⎝−12nγ∂2u∂τ2γ∣∣ ∣∣ri+τγ∂2u∂τγ∂nγ∣∣∣ri⎞⎠, (46) Jk,regi(u,q)=∫Γk{(∇iq∗)u−(∇iq∗L)u(ri)−(∇iu∗)q}dΓ, (47) ←→Ci=⎛⎜ ⎜ ⎜⎝Diffγ:A−B[(2θγ+sin2θγ4π]Diffγ:A−B[(−cos2θγ4π+u∗(Lγ)]Diffγ:A−B[(−cos2θγ4π−u∗(Lγ)]Diffγ:A−B[(2θγ−sin2θγ4π]⎞⎟ ⎟ ⎟⎠, (48) Diffγ:A−B[(fγ]≜fA−fB, (49)

where denotes the fundamental solution to the Laplace equation that satisfies Eq. (6), is its normal derivative with respect to at ; and , , and are the tangential unit vectors along , the azimuth angle of from , and the length of , respectively.

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

 s(k,l),regi=∫Γk∇i∇u∗(r,ri)ϕ(k,l)(r)⋅ndΓ, (50) s(k,l),regL,i=∫Γk∇i∇u∗L(r,r