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 higherorder 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 nonconforming elements, which are also called discontinuous elements. In the nonconforming element, collocation nodes that represent and do not coincide with geometric nodes, but they do so in the conforming element. The nonconforming 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 nonconforming 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 threedimensional 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 nonconforming 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 leastsquare 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 nonuniqueness 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 leastsquare method since the set of simultaneous equations becomes an overdetermined equation. The other approach is referred to as the BurtonMiller method BurtonMiller:1971 ; Benthien:1997 ; Diwan:2013 ; Langrenne:2015 . Burton and Miller BurtonMiller: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 nonuniqueness problem with comparisons of other methods, including the CHIEF and BurtonMiller 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 zerothickness 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 dualBEM.
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 dualBEM.
In the application of HBIEs, the regularization of hypersingularity is a key issue. The authors developed an analytical regularization of the twodimensional 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 leastsquare 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 partialHBIE 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 subnodes for a double node are then shown; and lastly, the rank deficiency conditions are discussed based on comparisons between two equations for the two subnodes.
2.1 Discretization of CBIEs
A complexvalued time harmonic scalar wave satisfies the following Helmholtz equation with an assumed time factor , where is an imaginary unit and is the angular frequency:
(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
(2) 
where the differential operator operates only on , but not on . This equation can be solved analytically. In twodimensional problems, the outward propagating wave that obeys Eq. (2) is
(3) 
where the function is a second kind 0th 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),
where denotes a boundary surrounding , is the position of the source point on the boundary, is the position of the field point,
is the outwardpointing 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:(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 2dimensional problems. In addition, when we introduce a fundamental solution to a Laplace equation, , which satisfies
(6) 
the coefficient can be expressed by a boundary integral:
(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:
(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
(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:(10) 
where denotes the number of nodes in a boundary element; e.g., for the linear element, for the secondorder element. A global node number can be mapped from the local node numbers by a permutation matrix :
(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
(12) 
where and are defined as
(13)  
(14) 
The integral on the righthand 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
(15)  
(16)  
(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 subnodes at the same position; . The node at is connected to both and a zerosized element between and , and vice versa. The normal direction of the subnodes 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 zerosized element is identically zero regardless of the normal direction, we do not need to define the normal direction for the zerosized 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 subnode in the same way as ordinary nodes, and the boundary condition is imposed for each subnode. 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 subnodes 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 zerosized elements is zero, the righthand sides of Eq. (7) for and are the same. Therefore, the coefficients and must have the same value:
(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 :
(19) 
Therefore, summarizing the discussions in the last paragraph in Sec. 2.1 and this result, we obtain the relations for the subnodes as
(20)  
(21) 
2.3 Rank deficiency conditions in CBIEs
By using the double node technique, the boundary nodes including subnodes 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
(22)  
(23)  
(24) 
where the overbars 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 multimedia 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
(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 subnodes 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 subnodes belong to or are presented.
The combinations of the boundary conditions where the two subnodes, 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 subnode equations of Eq. (8) for and can be arranged so that the terms including known values are moved to the righthand 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 lefthand sides and the righthand sides are written as
(28)  
(29) 
Comparing the coefficients of the terms on the lefthand side of the two subnode 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 lefthand 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 righthand sides, which are not related to rank reduction. Therefore, the lefthand sides of the two equations are identical, and the rank is always reduced by these two subnode equations; i.e., the coefficient matrix becomes a singular matrix.
2.3.2 Case of two Neumann conditions
As in the previous subsection, two subnode equations of Eq. (8) in the case of and can be arranged as
(30)  
(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 subnode equations is not reduced in the case of and . In addition, the righthand sides are the same since they do not include the Kronecker’s delta. This means that the results of these two subnode equations involve the relation .
2.3.3 Case of coupled Dirichlet and Neumann conditions
In the case of and , two subnode equations are
(32)  
(33) 
Since the second terms on the lefthand sides of the above equations are different, the rank of the coefficient matrix is not reduced in the case of and .
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 subnodes. Eliminating from Eq. (8) using Eq. (24), and arranging the equation, we obtain:
(34) 
When , the term with the factor for unknown remains in the third terms of the lefthand side. This satisfies the condition described at the end of Sec. 2.3.3, which ensures that the equations of the two subnodes 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 subnodes at the corner; two subnodes labeled and do not belong to the interface boundary, and the other two subnodes labeled and belong to the interface boundary. There are eight quantities related to these four subnodes. Of the eight quantities, two quantities at the subnodes 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
(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 subsections, the coefficients unrelated to the subnodes do not affect rank reduction; therefore, we only consider the nature of the submatrix 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 subnode equations for , , , and , are written as
(36) 
By applying Eqs. (18), (19) and (20) to the submatrix to eliminate coefficients with the subscript and , the submatrix is rewritten as
(37) 
Furthermore, by applying elementary matrix operations, the submatrix is transformed to
(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 nonzero components in both the second and the fourth row are only found in the fourth column; therefore, the rank of this submatrix is reduced by one.
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 subnodes at the corner node. After applying Eqs. (18), (19) and (20), and eliminating the quantities related to and , the terms related to the subnodes in the four subnode equations are expressed as
(41) 
By applying the elementary matrix operations to the coefficient matrix, the matrix is transformed to
(42) 
We can understand that the rank of this matrix is reduced by one since the determinant of a 2x2 submatrix composed of nonzero elements in the second and the fourth row in the above submatrix 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 subnode 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 subnode equations of a double node for CBIEs. The rank of the coefficient matrix is reduced in the following cases:

both subnodes of a double node are imposed by Dirichlet conditions,

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

subnodes 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 subnode equations or the submatrix associated with the subnodes 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 :
(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 .
3.1 Regularized gradient field
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.
The gradient is expressed by a 2x2 dyadic tensor
and a vector as(44) 
The second term can be evaluated by two types of boundary integrals as
(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 lefthand side of Eq. (44) and . The components of , and are written as
(46)  
(47)  
(48)  
(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 nonsingular integral in Eq. (47) is expressed by the shape function as
(50)  
Comments
There are no comments yet.