1 Introduction

In the conventional boundary integral equation (BIE) for the interior problem of plane elasticity, the solution for the displacement along the boundary is generally not unique in the Neumann boundary value problem (BVP). Clearly, the relevant integral operator acting upon the displacement is not invertible. Alternatively speaking, after discretization, the relevant influence matrix is singular. People have been trying to find all possible ways for solving the mentioned problem with a unique solution.

Blazquez et al.^[1] theoretically and numerically studied the removal of rigid body motions in the solution of elastostatic problems by use of the direct boundary element method (BEM). Lutz et al.^[2] studied the regularization approach for the BEM which usually led to non-symmetric matrices. Vodicka and Mantic^[3] deduced an explicit equation for the evaluation of critical scales for a given boundary, when the single-layer operator failed to be invertible. It was proved that there were either two simple critical scales or one double critical scale for any domain boundary. Vodicka et al.^[4] studied the removal of the non-uniqueness in the solution of elastostatic problems. Chen et al.^[5] pointed out that there were two kinds of rank-deficiency problems in the BEM. For the interior problem of plane elasticity, one rank-deficiency problem happened in the Neumann BVP. In this case, the additional rigid body motions in the displacement solution were the result of the rank-deficiency in the influence matrix. For the exterior problem, the other rank-deficiency problem happened in the Dirichlet BVP. In this case, the non-unique solutions of the traction along the boundary were the result of the rank-deficiency in the influence matrix when the used scale coincided with the degenerate scale. Chen et al.^[6] added a free constant and an extra constraint to the traditional method of the fundamental solution (MFS), and clearly explained why the free constant and extra constraints were both required by use of the degenerate kernel for the closed-form fundamental solution.

Vodicka and Petrik^[7] pointed out that the degenerate scale was usually encountered from a size effect which caused non-unique solutions for a certain type of BVPs. There were three degenerate scales for a general anisotropic material. If the standard plane strain isotropic elasticity was considered, two degenerate scales would be obtained. Based on a complex variable boundary integral equation (CVBIE), Chen^[8] provided a numerical solution for the elastic inclusion problem. The solution depended on the formulation of the interior problem and the exterior problem. The properties for the operators derived from the interior problem and the exterior problem were studied in detail. Chen and Wang^[9] studied the properties of the integral operators and the solutions for the CVBIE in plane elasticity for a multiply-connected region. Chen et al.^[10-11] revisited the direct and indirect BEMs to examine the uniqueness of the solution by introducing the Fichera idea and the self-regularized technique. The boundary integral equation method in conjunction with the degenerate kernel, the direct searching technique (singular value decomposition), and the only two-trials technique were analytically and numerically used to find the degenerate scales, respectively.

Based on the complex variable, a conventional CVBIE in plane elasticity is provided in this paper. We use the Somigliana identity with a particular source point in the studied interior domain of the fundamental stress field, and assume that the displacements vanish at this point in the physical stress field. After taking this step, we will obtain an additional integral equality. By adding both sides of this integral equality to both sides of the conventional CVBIE, we will obtain the amended BIE. In the amended BIE, the structures of both the integral operators or the operator acting upon the displacement and the operator acting upon the traction are changed. The method based on the discretization of the amended BIE is called the amended influence matrix method. With this method, for the Neumann BVP of an interior region, we obtain a unique solution for the displacement. The rigid body motion solution for the displacement is removed in the suggested method. Several numerical examples are provided to prove the efficiency of the suggested method.

2 Analysis

Based on the complex variable, a CVBIE in plane elasticity is provided in this paper. After using the Somigliana identity between a particular fundamental stress field and a physical stress field, an additional integral equality is obtained. By adding both sides of this integral equality to both sides of the conventional CVBIE, we obtain the amended BIE. The amended BIE is useful in the present study.

2.1 Complex variable BIE for interior BVPs

In recent years, the CVBIE in plane elasticity has been suggested to solve different BVPs by Chen^[8] and Chen and Wang^[9]. The CVBIE is equivalent to the formulation based on real variables^[12-13]. The suggested CVBIE is more straightforward, since all involved arguments are expressed in an explicit form.

In the suggested CVBIE in plane elasticity, the limit process for t₀ (t₀ ∈ S⁺, a domain point) approaching t₀ (t₀ ∈ Γ, a boundary point) is easy to see (see Fig. 1)^[8-9]. This process can be derived exactly from the behavior of the Cauchy type integral. However, this behavior is not easy to derive in the conventional BIE.

Particularly, it has been proven that if the load applied on the contour in the exterior BVP is not in equilibrium, the obtained results from the conventional BIE are inaccurate^[8-9].

Moreover, the usage of the CVBIE or the conventional BIE will give the same result for the studied problem in this paper.

After some manipulation, the BIE for the interior problem takes the following form^[8] (see Fig. 1):

(1)

where t₀ ∈ Γ for the interior problem, and Γ denotes the boundary of the interior region S⁺ (see Fig. 1). In Eq. (1), U(t) and Q(t) denote the displacement and traction along the boundary Γ, which are defined by

(2)

Two elastic constants and two kernels are defined by

(3)

(4)

where G is the shear modulus of elasticity, υ is the Poisson ratio, and

In this paper, the plane strain condition and υ =0.3 are assumed. In Eq. (1), the increment "dt" is going forward in an anti-clockwise direction.

For the physical field, we assume that the traction Q(t)=σ_N (t)+iσ_NT (t) is applied along the contour Γ (see Fig. 1). The applied loading must satisfy

(5)

In the analysis, we define the following three rigid body motions:

(6)

where t=x_t +iy_t.

Clearly, the relevant stresses vanish for the three rigid body motions.

From the left-hand term in Eq. (1), we have the following property^[8-9]:

(7)

where U^(j) (j=1, 2, 3) denote the three rigid body motions.

After discretization, at many discrete points (t₀ =t_0j, j=1, 2, …, N), the integral equation (1) can be reduced to

(8)

where the matrix [M_U] is the matrix acting upon the displacement vector, and the matrix [M_Q] is the matrix acting upon the traction vector.

2.2 Formulation of the amended influence matrix method

Now, we introduce the additional fundamental stress field with the concentrated forces p_x and p_y applied at the point z=t_p(t_p ∈ S⁺) (see Fig. 2). For the physical stress field, assume the following condition:

(9)

Since the point z=t_p(t_p∈ S⁺) is one definite point by t_p∈ S⁺, the assumed condition shown by Eq. (9) is reasonable.

After using the Betti reciprocal theorem between the fundamental stress field and the physical stress field shown in Figs. 2 and 3, we find

(10)

where t_p ∈ S⁺.

From the left-hand term in Eq. (10), we have the following property^[8-9]:

(11)

where j=1, 2, 3, and t_p ∈ S⁺.

Clearly, Eq. (11) is equivalent to

(12)

where j=1, 2, 3, and t_p ∈ S⁺.

In Eqs. (11) and (12), U^(j) (j=1, 2, 3) denote the three rigid body motions.

After using the condition (9), from Eq. (10), we find

(13)

where t_p ∈ S⁺.

After discretization, the equality shown by Eq. (13) can be reduced to

(14)

Note that the point t_p is a assumed point (t_p ∈ S⁺), and the two matrices [M_UC] and [M_QC] depend on the position of t_p(t_p ∈ S⁺).

Obviously, Eq. (14) can be realized as an equality rather than the discrete form of some BIE. The detailed explanation for the formulation of the equality shown by Eq. (14) can be referred to Appendix A.

By multiplying a factor -γ(0 < γ < 1) to both sides of Eq. (14) and adding them to both sides of Eq. (8), we have

(15)

where

(16)

It is easy to see that, the amended portion in the matrix [M_UA] is -γ [M_UC] (0 < γ < 1), and the amended portion in the matrix [M_QA] is -γ [M_QC] (0 < γ < 1).

Instead of Eq. (8), the governing equation in the formulation is changed into Eq. (15). Therefore, the suggested method is called the amended influence matrix method.

The purpose for introducing the value of γ is to examine whether the suggested method is valid for some range of γ or not. Physically, the option -γ =-1 or γ =1 corresponds to one unit negative force applied at the point t_p (see Figs. 2 and 3).

Since the matrix [M_U] in Eq. (8) is modified into [M_UA] in Eq. (15) (note that [M_UA]=[M_U]-γ [M_UC]), the invertible properties for the matrices [M_U] and [M_UA] are quite different. This can be seen from the following results. From Eq. (7), we have

(17)

From Eqs. (7) and (12), we have

(18)

where j=1, 2, 3.

From Eq. (18), we find that {U}={U^(j)} (j=1, 2, 3) or the three rigid body motion displacements and {Q}={0} are no longer the solution for [M_UA]{U}=[M_QA]{Q}. Moreover, for [M_UA]{U}={0}, we only have the trivial solution {U}={0}. Therefore, it is expected that the inverse matrix [M_UA]^-1 exists. Alternatively speaking, for the Neumann BVP, we have a unique solution for {U}, which is evaluated from [M_UA]{U}=[M_QA]{Q}.

It is known that there are three degrees of freedom or three rigid body motions for the displacement solution in the Neumann BVP for the interior problem^[1]. To remove the mentioned three rigid body motions, the point support method is generally used. In the point support method, it is an effective way to let the values of the three diagonal elements in the matrix [M_U] be extremely huge. The modified matrix for [M_U] is denoted by [M_U^*]. Therefore, the discrete form of the BIE or [M_U]{U}=[M_Q]{Q} shown by Eq. (8) will be changed into [M_U^* ]{U}=[M_Q]{Q}.

From the aforementioned derivation, we see that there is no direct relation between the suggested method and the point support method. Alternatively speaking, it is not easy to compare the two methods within the framework of the general BIE formulation.

In the usual way for avoiding the rigid body motion by use of the support constraint technique, we should assume that the values for the three diagonal elements in the influence matrices are huge. Therefore, the calculated results for the displacements will be affected by the assumed three diagonal elements significantly.

2.3 Numerical examination for the amended influence matrix method

As seen from the above-mentioned analysis, the governing equation in the amended influence matrix method has been shown by Eq. (15) or [M_UA]{U}=[M_QA]{Q}.

To examine the validity of the suggested method, we must provide an exact solution beforehand. In the exact solution, the displacement and traction vector assumed along the boundary Γ are denoted by {U}_ext and {Q}_ext, respectively.

First, we consider the Neumann BVP for an elliptic plate in plane elasticity (see Fig. 1). It is known that, the integral operator on the left-hand side of Eq. (8) is not invertible, and the inverse matrix [M_U]^-1 does not exist. However, after making modification to the influence matrix, the inverse matrix [M_UA]^-1 exists.

In the Neumann BVP, if {Q}_ext is substituted into the right-hand side of Eq. (15), we have

(19)

where {f₁}=[M_QA]{Q}_ext.

From Eq. (19), we can get a solution for the displacement vector as follows:

(20)

It is known that the relevant solution for the displacement vectors may be different from each other by the three rigid body motions. To verify the accuracy of the suggested method, a technique is suggested below.

To evaluate the stress component σ_T (see Fig. 1), we suggest the following technique^[8]. In fact, in the plane strain case, the strain component ε_T (in the T-direction) can be expressed by Hooke's law (see Fig. 1) as follows:

(21)

or

(22)

In Eq. (22), the component σ_N is from the input data {Q}_ext, and ε_T is the strain in the T-direction, which can be evaluated from the numerical solution of the displacements along the boundary or from {U}_num. In fact, the elongation of a boundary element can be found from the displacement solution, and the strain ε_T can be evaluated accordingly. Thus, the values of σ_T at many discrete points along the boundary can be evaluated. The stress component vector σ_T evaluated from Eq. (22) is denoted by {σ_T }_num. After comparing this vector {σ_T }_num with the vector {σ_T }_ext, which comes from the input data, we can judge the achieved accuracy accordingly.

Secondly, we consider the Dirichlet BVP for an elliptic plate in plane elasticity (see Fig. 1). If the used size does not coincide with the degenerate scale, we can obtain the solution for the Dirichlet BVP in the interior region directly from Eq. (1) or Eq. (8) in the discrete form. Therefore, there is no need to use the amended influence matrix method in the Dirichlet BVP. In the Dirichlet BVP, if {U}_ext is substituted into the right-hand side of Eq. (15), we have

(23)

where {f₂ }=[M_UA]{U}_ext.

From Eq. (23), we can get a solution for the displacement vector

(24)

Similarly, we can evaluate {σ_T}_num from Eq. (24). However, the component {ε_T } is obtained from the input vector {U}_ext, and the component {σ_N } is obtained from the computed vector {Q}_num. Finally, we can compare {σ_T }_num with {σ_T }_ext.

3 Numerical examples

Two numerical examples are provided to prove the efficiency of the suggested method. In all numerical examples, the plane strain condition and υ =0.3 are assumed.

Example 1 In the first example, we provide a numerical solution for the Neumann BVP for a finite elliptic plate under the conditions: (ⅰ) b/a=0.5; (ⅱ) γ =0.25 or γ =0.75; (ⅲ) t_{p, j} =((acos α_j)/2, (bsin α_j)/2), where α_j =(j-1)π /6 (j=1, 2, 3, 4); (ⅳ) σ_y =p, σ_x=0, and σ_xy =0; (ⅴ) 2Gu=-υ p(x-x_p), and 2Gv=(1-υ)p(y-y_p), where t_p =x_p +iy_p (see Fig. 4).

The additional source points t_{p, j} =((acos α_j)/2, (bsin α_j)/2) (j=1, 2, 3, 4) are indicated in Fig. 4. The technique for the discretization of the BIE suggested by Chen^[8] is used to the present example. A total of 96 divisions for the elliptic boundary are used in the discretization.

After solving the Neumann BVP, the calculated results for σ_T are expressed as follows:

(25)

For the following cases: (ⅰ) b/a=0.5; (ⅱ) σ_y =p, σ_x =0, and σ_xy =0; (ⅲ) γ =0.25 or γ =0.75; and (ⅳ) t_{p, j} =((acos α_j)/2, (bsin α_j)/2), where α_j =(j-1)π /6 (j=1, 2, 3, 4), the calculated results of f(θ) (=σ_T/p) for γ =0.25 and γ =0.75 are plotted in Figs. 5 and 6, respectively. The exact solutions are also plotted in those figures.

Clearly, the accuracy of the computed results for f(θ) (=σ_T /p) depends on two factors: (i) γ =0.25 or γ =0.75; (ⅱ) t_{p, j} =((acos α_j)/2, (bsin α_j)/2), where α_j =(j-1)π/6 (j=1, 2, 3, 4. In all the eight combinations for γ and t_{p, j} shown in Fig. 4, we find that the relative maximum error is less than 0.6%.

Example 2 In fact, it is not necessary to use the amended influence method for the Dirichlet BVP for a finite elliptic plate. However, in order to examine the validity of the suggested method, we present the solution for the Dirichlet BVP in the second example.

In Example 2, we provide a numerical solution for the Dirichlet BVP for a finite elliptic plate (see Figs. 1 and 4). The same conditions used in the first example are used in this example.

After solving the Dirichlet BVP, we obtain

(26)

For the γ =0.25 or 0.75 case, the calculated results for g(θ) (=σ_T/p) are plotted in Figs. 7 and 8, respectively. The exact solutions are also plotted in those figures.

Clearly, the accuracy of the calculated results for g(θ) (=σ_T/p) depends on two factors: (ⅰ) γ =0.25 or γ =0.75 and (ⅱ) t_{p, j} =((acos α_j)/2, (bsin α_j)/2), where α_j =(j-1)π /6 (j=1, 2, 3, 4). In all the eight combinations for γ and t_{p, j}, we find that the relative maximum error is less than 0.2%.

4 Conclusions

In the amended influence matrix method, the matrix [M_U] acting upon the displacement is modified to [M_UA]([M_UA]=[M_U ]-γ [M_UC]), where -γ [M_UC] is the amended portion. It has proven that the inverse matrix [M_UA]^-1 exists. Therefore, we can evaluate a unique solution for the displacement along the boundary immediately in the Neumann BVP. It is seen that the error in the calculation mainly depends on two factors: (ⅰ) the value of γ, e.g., γ =0.25 or 0.75, and (ⅱ) the positions of t_p (see Fig. 4). The numerical examples prove that the error in the computation is extremely small.

A real calculation has been proven that a sufficient accurate result has been achieved for a wider range with the four different t_p points in Fig. 4 and the two different values of γ (0.25 or 0.75) in Figs. 5, 6, 7, and 8, where a total of 8 combinations are used in the examination. A rather complicated calculation from calculation from different t_p points and γ values with a total of 8 combinations provides the accurate result. This is an alternative way to show the efficiency of the proposed method.

In the author comment, the suggested method can be used to the interior BVP for 2D Laplace equations and 3D BIEs of elasticity.

Appendix A

Derivation for the equality [M_UC]{U}=[M_QC]{Q} shown by Eq. (14) (see Fig. 2)

After using the following condition:

(A1)

from a discretization of Eq. (13), we obtain

(A2)

Obviously, Eq. (A2) can be realized as an equality rather than the discrete form of some BIE.

The derivation for Eq. (A2) needs to have a detailed explanation.

In the first step, for a definite point z=t_p, where t_p denotes one definite point in the domain, the equality shown by Eq. (13) can be written in the discretization form as follows:

(A3)

where the displacement and traction vectors {U} and {Q} are defined by

(A4)

In Eq. (A4), (u_j, v_j) represents the displacement for the jth element at the boundary, and (σ_{N, j}}, σ_{NT, j}}) represents the traction for the jth element at the boundary. Note that we use N elements in the derivation. In addition, one complex value equals two values in the real expression.

In the second step, one may extend the matrix relation shown by Eq. (A3) into a wider range as follows:

(A5)

Finally, Eq. (A5) can be rewritten as follows:

(A6)

From Eq. (A5) and Eq. (14), we find the following properties for the elements in the matrices [M_UC] and [M_QC]:

(A7)

where j=2, …, N, and k=1, 2, …, 2N.

Since the matrices [M_U] and [M_UC] (or [M_Q] and [M_QC]) have the same dimension, the linear combination shown by the following equations

(A8)

are reasonable.