A Subdivision Transformation Method for Weakly Singular Boundary Integrals in Thin-structural Problem

source points, the shape information of the elements and the nearest distance from the source point to the integral element, a subdivision technology is constructed at first. With this subdivision technology, the integral element can be divided into several integral blocks with good shapes. And then, a simpler coordinate transformation method is constructed to remove the weak singularities of the integral blocks obtained by the subdivision technology. Compared with the conventional polar coordinate transformation method, the present transformation method does not need to calculate their integral interval, which is simpler and more effective to implement. Finally, the paper gives three numerical examples to verify the accuracy and validity of the present method.


Introduction
Thin structures have very good bearing performance, high strength and stiffness, small weight and can bear a considerable load with a small thickness, which are widely used in engineering structures [1][2][3]. Such as components for automobile and ship structures that improve loadbearing capacity and reduce weight, coating structures that improve heat and wear resistance of metal surfaces, as a buffer and energy absorption component used in aerospace, rail transportation fields, etc. To ensure the force and material strength of thin structure, it is necessary to carry out numerical analysis of mechanical properties according to its geometric parameters.
The boundary element method (BEM) is a more accurate and effective method developed after finite element method. Different from the basic idea of the finite element method (FEM), BEM only divides the elements on the boundary of the domain and approximate the boundary conditions with the functions satisfying the governing equation [4][5][6][7]. Therefore, compared with the FEM [8,9], BEM has the advantages of dimensionality reduction, fewer elements and simpler data preparation, which is suitable for solving stress concentration, infinite domain and semi-infinite domain problems [10,11]. However, as the singularities of the kernel function in the boundary integral equation (BIE), the computational accuracy of weak singular integral directly affects the final results, so additional processing schemes must be given to remove these singularities to ensure the successful solution of BEM. Especially for thin-structural problems, poorly shaped elements (such as elements with large angles or narrow lengths) will appear in the discrete geometric model of thin structure [12][13][14], which will seriously affect the accuracy of singular and near singular integrals. Many works have been demonstrated that the near singular integrals can be accurately evaluated by analytical and semi-analytical methods [15,16], Sinh and other nonlinear transformation methods [4,[17][18][19], etc. Therefore, accurate and effective calculation of weakly singular integral in BIE is the key to the implementation of boundary element analysis [20][21][22][23].
At present, there are many methods have been developed to handle the weak singular integrals, such as integral simplification method [24], element subdivision methods [25][26] and polar coordinate transformation methods [20,[27][28][29][30][31]. When the position of the source point is close to the boundaries of the element, the computational accuracy cannot be guaranteed only by using the coordinate transformation method [32]. To obtain high calculation accuracy, the integral element needs to be subdivided into several integral blocks with good shape. And there are many element subdivision methods, such as the singular points directly connected to the element vertexes [12], interval block method [12], tree structure methods [25], etc. When these subdivision methods are employed, some integral blocks with poor shape (such as the integral blocks with large angle and large aspect ratio) will be subdivided, which increases the difficulty in dealing with singular integral and decrease the computational accuracy. In addition, the polar coordinate transformation is the key to removing the weak singularities of the integral of the kernel. It converts the surface integral into the double integral of circumferential and radial direction, and then transforms the integral of radial direction to eliminate the singularities. Therefore, when this transformational method is employed, it needs to recalculate the integral interval of the integral block in the radial and circumferential directions every time, which is more troublesome to implement [32].
In this paper, the weakly singular boundary integrals of the displacement kernel function in the BIE of elastic thin-structural problems are considered. When the BIE is discrete by several elements, for the numerical integration of a discrete element, firstly, the subdivision techniques of the quadrilateral elements are developed according to the location of the source point, the element shape and the nearest distance from the source point to the element; and then, based on this subdivision technology, a simpler coordinate transformation method is constructed to remove the weak singularities of the singular integral blocks obtained by the proposed subdivision technology. Compared with the conventional polar coordinate trans-formation method, the presented method does not need to calculate their integral interval and is simpler and more effective. Finally, the subdivision technology and coordinate transformation method are applied to the BEM analysis for the elastic structure to realize an accurate and effective solution to the thin-structural problems.
The outlines of the paper are as follows: the BIE and its numerical discretization of the thin-structural problems are described in section 2. In section 3, the treatment of the weakly singular integrals is presented in detail. The verification results of numerical examples are given in Section 4. At last, the conclusions are given.

The boundary integral equations for thin-structural problems
The BIE of 3D elastic analysis for thin-structural problems is shown in Eq. (1) [33]. P and Q are the source and field points; cij(P) is a constant associated with the boundary type; * ( , ) ij u P Q and * ( , ) ij t P Q are the kernel function of displacement and traction, respectively, and their expressions are shown in Eq. (2) [34]: , where: G and v denote the shear modulus and Poisson's ratio; r is the distance from the source point to field point; uj(Q) and tj(Q) represent the variables of displacement and traction on the boundary Γ; ni and ni are the components of the external normal vector in i and j direction, respectively. By discretizing the integral equation Eq. (1) with N elements, Eq. (1) can be expressed as: in which, n is the number of nodes on each element, and Nα(Q) is the shape function of the αth node on the element. It can be seen from Eq. (2), the integrals of the displacement kernel function have weak singularities, while the integrals of the traction kernel function have strong singularities. The numerical results have shown that strongly singular integrals can be accurately calculated by the approximate expansion method and regularization method [35,36]. In this paper, the weakly singular integrals of the displacement kernel function are mainly considered.

The processing method of weak singular boundary integral
To construct the weak singularity elimination technique proposed in this paper, one discrete element Γe is considered at first, and the integral of the displacement kernel function is simplified into the following form (Eq. (4)).
Where f(P, Q) is a smooth function without singularity; ϕ(Q) is the corresponding shape function: We can see in Eq. (4), the weak singularity will arise when r is equal to zero. Additional schemes for handling weak singular integrals need to be considered. To achieve the purpose of solving the weakly singular integral accurately, the element subdivision technique associated with the weak singular boundary integral needs to be considered at first.

The element subdivision technique for solving singular integrals
As the uncertainty of the position of the source point (which may fall on the interior, edges and vertices of the element), the subdivision method will be different. When the position of the source point is given, the nearest distance d from the source point to the edges of the element can be calculated. And then a very small quadrilateral with a length of 0.2d is constructed near the source point and the source point is included in it.
In which, the selection of 0.2d is an empirical value obtained from a numerical test to avoid mutual interference between integral blocks. The specific operations are as follows: Firstly, to ensure that the shape of the singular integral block is good, we need to choose a regular quadrilateral that includes the source points. The choice of side length of the quadrilateral is very important. Too large selection of the side length will lead to a too large singular integral region, which requires too many Gaussian points and reduces the computational efficiency. In addition, the interference between integral blocks should be avoided. Therefore, this paper starts from 0.5d. In the selection, it is found that if the side length of the constructed quadrilateral is large, especially for the element whose source point S falls near the center, small integral blocks will be divided near the edge, and the density of integral blocks cannot be guaranteed, which is the reason why 0.2d is selected as the side length of the small quadrilateral in this paper.
By extending the four sides of the small quadrilateral and intersecting the integral element, several quadrilateral integral blocks can be obtained. Take the quadrilateral element as an example, the specific subdivision schemes are shown in Figs. 1-5. Fig. 1 is the result of the first step subdivision when the source point falls on the vertex of the integral element.  In the second step, the vertices of the small quadrilateral containing the source point are connected with the source point, and several triangular integral blocks containing singular points S can be obtained (as show in Fig. 4). For the sake of integral, the quadrilateral R1-R8 without source point in Figs. 1 -3 should also be further subdivided. The subdivision method is as follows Fig. 5): If the value of l /d is less than 1, the integral block can be seen as a regular integral element. If l /d is greater than 1, the integral block is equally divided into four sub-elements. Repeat for subelements until all sub-elements are regular integral elements. Then the normal Gauss integral can be applied to the above regular integral elements.
In addition, if the shape of the integral element is narrow, the element can be subdivided according to the proportion of the longest (Lmax) and shortest edge (Lmin) of the element (as shown in Fig. 6, if Lmax/Lmin < 2, the element cannot be subdivided (Fig. 6, a). If 2 < Lmax/Lmin < 3, the element needs to be divided into two integral sub-elements We can see from the above element subdivision method, the source points are first contained in a small quadrilateral, and then the vertices of the small quadrilateral are connected to the singular points, and four singular integral blocks are obtained at last. The integral block without source points can be solved by conventional Gaussian integral. However, due to the singularity of the triangular integral block containing the source point, a simple Hammer integral cannot eliminate the singularity in the integrand, so it needs special treatment.
Taking Eq. (4) as the integral form, to remove the weak singularities of the integral block containing the source point, the (α, β) local coordinate system is employed at first [37], and the corresponding mapping system is shown in Fig. 7. The specific process of the coordinate transformation method is as follows: From Eqs. (5) -(9), we can see the Jacobian obtained by this coordinate transformation method contains quasi-zero factor α, which just offsets the quasi-zero factor in the denominator of integral equation Eq. (4), and the singularity can be eliminated. Moreover, this coordinate transformation method can transform the integral interval of α and β to [0,1] directly, which does not need to calculate the integral interval of each integral block. Compared with the traditional polar coordinate transformation, it is simpler and more effective.

Numerical examples
In this part, three numerical examples are given to verify the accuracy and effectiveness of the proposed method. Firstly, the accuracy of the proposed method on one element is verified, where different source positions and element shape types are considered. Finally, the present method is employed to analyse the elastic analysis of the thin-structural problem. The computational formula of relative error is shown in Eq.   Table 1. Table 1 The results of singular integrals of quadrilateral elements  Table 2.
When a = 5, the four vertices of the quadrilateral element are given as  Table 4. Table 2 The numerical results of singular integrals of quadrilateral elements (a=2)

Source points
Gaussian points  Table 3 The numerical results of singular integrals of quadrilateral elements (a=5)

Source points
Gaussian points  Table 4 The numerical results of singular integrals of quadrilateral elements (a=10)  Table 1 - Table 4 and Fig. 6, it can be seen when using the same Gaussian integral points, the accuracy of the proposed method has remained high with the increase of the aspect ratio of the integral element (even if the integral element becomes narrow and long). With the proposed method, the weakly singular integrals in the BIE can be accurately evaluated, at the same time, the influence of weak singular integrals can be eliminated in the solution of thin structural problems.

The thin plate structure
In the second example, the proposed method is employed for the elastic analysis of thin plate structure. As shown in Fig. 10, the length, width and height of the thin plate are l=10 mm and h=1 mm. Young's modulus and Poisson's ratio are E=1 MPa and v=0.25, respectively. To facilitate comparison, the quadratic displacement fields with an analytical solution are imposed on the boundary of the thin plate and the formulas are expressed as follows: The geometrical model of the plate is discretized into 140 quadrilateral elements (Fig. 11). The flanks of the plate are discretized by several slender quadrilateral elements. A series of sample points are selected on line x=10

The thin plate with a hole
To further verify the computational accuracy and generality of the proposed method, a thin plate with a hole is considered in the third example. The geometric parameters of the thin plate are as shown in Fig. 16. The Young's modulus and Poisson's ratio are E = 1 MPa, v = 0.25. Eq. (13) is treated as the boundary conditions which are imposed on the external boundary of the thin plate. 10 1 10 x z y Unit：mm Fig. 16 The model of the plate with a hole Fig. 17 The discrete model of the thin plate with a hole 80 quadrilateral elements with different shapes areused in all boundaries of the geometric model ( Fig. 17 is the corresponding discrete model). The flanks of the model are discretized by several narrow elements. The sample points are selected on the centreline of the right side, and the numerical results obtained by using the proposed method and the traditional method are shown in Figs. 18 and 19. It can be seen from Fig. 18 and Fig. 19, compared with the traditional method, the numerical results obtained Fig. 18 The results of the traction in x and z direction by the proposed method are in better agreement with the analytical solution and have higher computational accuracy, while the numerical results of the traditional method have diverged, which further verify the accuracy of the proposed method. Fig. 19 The results of the Von Mises stress and the traction in y direction

The cylindrical revolution structures
To verify the validity and universality of the proposed method, a cylindrical revolution structure is considered in the last example. As shown in Fig.20

Conclusions
Based on the properties of the weak singular integrals in the boundary integral equations of elastic problems, an element subdivision technique combined with αβ coordinate transformation method is proposed in this paper. The proposed method is simpler to implement, no matter where the source point is in the element, no matter what the element shape is, accurate and stable computational results can be obtained. Then, the proposed method is employed to analyze the thin-structural problem, and the numerical computational results show that the proposed method can be used to evaluate thin-structural problems accurately and efficiently.