The Singular Function Boundary Integral Method for the 2-D and 3-D Laplace Equation Problems in Mechanics, with a Boundary Singularity
Miltiades C. Elliotis
Department of Mathematics and Statistics, University of Cyprus, Nicosia, Cyprus
Email address:
To cite this article:
Miltiades C. Elliotis. The Singular Function Boundary Integral Method for the 2-D and 3-D Laplace Equation Problems in Mechanics, with a Boundary Singularity. Pure and Applied Mathematics Journal. Vol. 5, No. 6, 2016, pp. 192-204. doi: 10.11648/j.pamj.20160506.14
Received: October 3, 2016; Accepted: October 15, 2016; Published: November 10, 2016
Abstract: In this study the Singular Function Boundary Integral Method (SFBIM) is implemented in the case of a planar elliptic boundary value problem in Mechanics, with a point boundary singularity. The method is also extended in the case of a typical problem of Solid Mechanics, concerning the Laplace equation problem in three dimensions, defined in a domain with a straight edge singularity on the surface boundary. In both the 2-D and 3-D cases, the general solution of the Laplace equation is approximated by the leading terms (which contain the singular functions) of the local asymptotic solution expansion. The singular functions are used to weight the governing equation in the Galerkin sense. For the 2-D Laplacian model problem of this study, which is defined over a domain with a re-entrant corner, the resulting discretized equations are reduced to boundary integrals by means of Green’s second identity. For the 3-D model problem of this work, the volume integrals of the discretized equations are reduced to surface integrals by implementing Gauss’ divergence theorem. The Dirichlet boundary conditions are then weakly enforced by means of Lagrange multipliers. The values of the latter are calculated together with the singular coefficients, in the 2-D case or the Edge Flux Intensity Functions (EFIFs), in the 3-D model problem, which appear in the local solution expansion. For the planar problem, the numerical results are favorably compared with the analytic solution. Especially for the extension of the method in three dimensions, the preliminary numerical results compare favorably with available post-processed finite element results.
Keywords: Laplace Equation, Boundary Singularity, Straight Edge Singularity, Singular Coefficients, Edge Flux Intensity Functions, Singular Function Boundary Integral Method
1. Introduction
In the past few decades, many numerical methods have been proposed for the solution of 2-D elliptic boundary value problems in Mechanics, with a boundary singularity. The lack of adequate accuracy and poor convergence, were some of the difficulties which appeared in the attempt of many researchers to solve this kind of problems. Remedies used were special mesh-refinement schemes, multigrid methods, singular elements, p/hp finite elements and many other techniques (e.g. [1-3, 5, 24, 25]). An extensive survey of the treatment of singularities in elliptic 2-D boundary value problems is provided in [8]. Of course, in the absence of a boundary singularity, there are many other efficient numerical techniques available in the literature (e.g. [1, 2, 4, 22, 29]) to be used for the solution of problems in Mechanics.
For the 2-D problems with a boundary singularity there are more sophisticated techniques in the literature, than those mentioned above, which incorporate the form of the local asymptotic expansion. For example, in [18] a singular finite element method is developed for Stoke’ s flow problems, in which special elements are employed, in the neighborhood of the singularity, with which the radial form of the local expansion is utilized, in order to resolve the convergence difficulties and improve the accuracy of the global solution. The interest in such two-dimensional Laplace equation problems with a boundary singularity [20], is motivated by the need to compute the singular coefficients of the local solution expansion which, in the area of Solid Mechanics, normally represents the Airy stress function Φ(r,θ) and expressed as follows:
(1)
In the above asymptotic expansion, μ_{j} are the eigenvalues and ^{}f_{j}(θ) are the eigenfuntions of the problem and they are known. The polar coordinates (r, θ) have as a center the singular point. Singular coefficients β_{j }are the primary unknowns in this kind of problems. They are also known as generalized stress intensity factors (GSIFs). These are determined by the boundary conditions on the boundary parts which are away from the singularity. They are of significant importance in many applications in the area of Fracture Mechanics [16]. Thus, after the calculation of the leading singular coefficients β_{j}, the approximation of Φ(x,y) is known (by converting into Cartesian coordinates) and the approximate stress state, at any position, is found from the known expressions: σ_{xx}=∂^{2}Φ/∂y^{2}, σ_{yy}=∂^{2}Φ/∂x^{2} and σ_{xy}=-∂^{2}Φ/∂x∂y.
In the past two decades, the Singular Function Boundary Integral Method (SFBIM) was developed [9-14, 17, 19], which belongs to the general group of collocation methods [26] and with which the unknown singular coefficients are calculated directly. Thus, the method gives directly the approximation of the Airy stress function (for solids) or the stream function (for fluids) in planar Laplacian and biharmonic problems. Basically, for the implementation of the SFBIM to solve 2-D problems with a boundary singularity, it is required that the solution is approximated by the leading terms of expansion (1). The Dirichlet conditions, on the boundary parts, which are away from the singularity, are enforced by means of Lagrange multipliers. The method has been tested on standard elliptic problems exhibiting exponential convergence and high accuracy with respect to the number of singular functions. This behavior of the method has been theoretically proved in [28]. The results obtained in the 2-D case have encouraged the extension of the method in 3-D problems with a straight-edge singularity [6, 13].
The 2-D model problem, treated in this study, is a Laplacian problem defined over a disk, with a boundary ∂Ω on which a quarter of a circle is missing and thus a boundary singularity at the centre, is created. This could be a heat transfer problem in two dimensions, but in this study it is a plane stress problem of a perfectly elastic solid, on which the interest is to find the unknown singular coefficients β_{j}, of expansion (1). As for the 3-D case, we consider a Laplacian boundary value problem in a 3-D domain. This is a model for an elastic cylindrical body, having a typical V-notch, made of an isotropic material which obeys to Hooke’s Law and which is subjected to certain physical conditions. The form of the local solution depends on the geometry in the neighbourhood of the straight-edge singularity and on the boundary conditions of the boundary parts which share the singularity. Also, the local solution is characterized by the presence of certain eigenpairs (arising from the 2-D problem) and the so-called Edge Flux Intensity Functions (EFIFs), which are the primary unknowns in the 3-D case [30]. The interest in such a problem is motivated by the need to compute generalized stress intensity functions for the V-notched solids loaded by static loads. For the solution of this class of 3-D problems, few methods have been proposed so far, such as the J-integral method [21], the B- and H-integral methods [23], and more recently the methods by Costabel et al. [7] and Yosibash et al. [32-35] and Zaltzman [36], in which the EFIFs are computed by means of a post-processing procedure in a p-version finite element scheme.
Based on [7] and [30], the solution u of the 3-D Laplace equation, in a domain with an edge singularity, may be written in terms of cylindrical coordinates (r, θ, z) as
(2)
In the above expansion, α_{κ} ∈ R and φ_{κ}(θ,α_{κ}) is the known eigenpair, of the two-dimensional Laplace operator. The functions φ_{κ}(θ,α_{κ}) are analytic in θ. The functions Α^{(}^{ακ)}(z) are the EFIFs. These are functions of z and for a large class of problems they are the primary unknowns. In Solid Mechanics [16], function u represents the first stress invariant (i.e. u=σ_{ii}).
In the SFBIM, the solution is approximated by the leading terms of the asymptotic expansion (1) for the 2-D case, or the asymptotic expansion (2) for the 3-D case. The leading terms of these expansions are also used to discretize the governing differential equation in the Galerkin sense. The discretized equations are reduced to boundary integrals by means of the Green’s second identity, for the 2-D case, or the Gauss divergence theorem, for the 3-D case. A particular feature of the SFBIM is that the Dirichlet conditions are weakly enforced by means of Lagrange multipliers. In two-dimensional problems the coefficients of the asymptotic expansion are constants; these are calculated directly by the SFBIM. As we have already mentioned above, in three dimensions the EFIFs are functions of the axial direction. In the present approach, these are approximated by polynomials the coefficients of which are primary unknowns and are calculated directly by the method.
The outline of the rest of this paper is as follows: in Section 2 we present both a 2-D Laplace equation problem, with a point boundary singularity and a 3-D Laplace equation problem with a straight-edge singularity. The asymptotic local solution expansion is also given for each one of these problems. In Section 3 the formulation of the SFBIM is presented both in two and three dimensions. Numerical results are given in Section 4 together with comparisons between the results obtained with the SFBIM and the exact solution. Especially for the extension of the method in three dimensions, the numerical results obtained are compared with post-processed finite element results available in the literature [30]. Finally, conclusions derived from this study are summarized in Section 5.
2. Governing Equation and Local Solution in 2-D and 3-D
2.1. The Planar Model Problem of the Laplace Equation
The 2-D model problem is a problem of a Laplace equation defined over a domain Ω. This domain is a very thin disk with radius R and with a re-entrant corner on its boundary, which has its tip on the center O of the disk (Figure 1). Thus, at point O there is a boundary singularity. Boundary parts OA and OB, which share the singularity, intersect vertically at O and each one has the Neumann condition: ∂Φ/∂n=0. Circular part AB has the Dirichlet condition: Φ=f(θ). The material of the disk is isotropic and linearly elastic. Thus, basically, with the specific geometry and boundary conditions of this problem we need to find an expression for the Airy stress function Φ(r,θ), which is the solution of the Laplace equation. The values of Young’s modulus and of Poisson’s ratio are not necessary for the solution of the problem. It must be mentioned that this model problem belongs to the general class of problems of two-dimensional plates with geometric singularities on the boundary.
The planar model problem of the Laplace equation, illustrated by Figure 1 is the following:
(3)
with the boundary conditions
(4)
For this specific model problem the Dirichlet condition on the boundary S_{C} has the following trigonometric expression:
Φ(θ)=(2/5)R^{2/3}cos(2θ/3)+(4/7)R^{4/3}cos(4θ/3)
where R is the radius of boundary S_{C}.
Using the separation of variables approach and the boundary conditions on S_{A} and S_{B} it can be easily found that the solution around the singular point is found to be Φ(r,θ)=βr^{μ}^{j}cos(μ_{j}θ). Furthermore, by employing the principle of superposition we obtain the local solution expansion, which is expressed in terms of polar coordinates (r,θ), centered at O:
(5)
Also, by considering the 2π-periodicity property of the trigonometric functions, singular coefficients β_{j} can be calculated as Fourier coefficients by using Euler’s formula. So, the exact solution of the 2-D model problem, defined by (3) and (4), is found to be
on Ω.(6)
Using the maximum principle theorem in its weak form (e.g. [15]) it can be shown that the above solution of the 2-D problem is unique.
2.2. The 3-D Model Problem of the Laplace Equation
The geometry and the boundary conditions of the 3-D model problem of the Laplace equation, are illustrated in Figure 2. The domain Ω is bounded by a cylindrical surface S_{C} around the z-axis, a circular sector S_{D} lying on the xy-plane, a flat circular sector S_{E} perpendicular to the z-axis and at a distance L from S_{D} (here L=2) and two flat boundaries S_{A} and S_{B} (one parallel to the xz plane and the other parallel to the yz plane) intersecting vertically on the z-axis and thus creating a V-notch, which is the straight-edge singularity of this problem.
This solid body is made of an isotropic linearly elastic material and has no voids inside. For the specific geometry and boundary conditions of this model problem, the solution u, which, as we have already mentioned, normally represents the first stress invariant, is independent of any combination of values of the Young modulus E and the Poisson ratio ν. As we know, these are important parameters for the calculation of stresses and strains in linearly elastic materials which obey to Hooke’s law.
This problem is a version of a more general problem which was first suggested by Yosibash et al. [30] and is as follows: Find u such that
in Ω,(7)
with the following boundary conditions
(8)
where
Naturally, the Laplace equation (7) is expressed in cylindrical co-ordinates, which are more convenient for both the analytical and numerical treatment of the problem.
In [30] it was shown that the solution u, of the three-dimensional general Laplace equation problem with a governing equation ∇^{2}_{3-D} u = 0 in Ω, such that the boundary conditions (8) hold, is obtained by augmenting the 2-D solution of the Laplace equation. Thus, according to [30], the solution in the neighborhood of the edge singularity, is given by
(9)
where (r,θ,z)∈Ω are cylindrical coordinates, with the z-axis along the boundary straight-edge singularity. Functions A^{(ak)}(z) are the EFIFs, which were described earlier. In this work they are chosen to be polynomials of degree N_{p} and with unknown coefficients a_{k,j} as follows:
(10)
The constants c_{ki} are given by
and the eigenfunctions φ_{k} and eigenvalues α_{k} are given by
(11)
Here, ω is the external angle defined by the flat boundaries ODEA and ODCB as shown in Figure 2; in the present work ω = 3π/2. Moreover, for the model problem of Figure 2, the first two EFIFs are known and they are of the form
for k=1 or 2
We must emphasize, at this stage, that the model problem of the 3-D case is selected in a way that the exact solution is known and given by a (finite) sum. Investigation of the uniqueness of this solution is similar with that of the 2-D case. Furthermore, since the only two non-zero EFIFs are polynomials, the accuracy of the EFIFs, computed by our method, can be measured by simply comparing the polynomial coefficients of the true and approximate EFIFs. In general, however, the true EFIFs are not always polynomials and the accuracy of the EFIFs computed by the numerical method, would need to be measured using some appropriate function norm.
3. The Singular Function Boundary Integral Method
3.1. The SFBIM in the 2-D Case
The first step of the method, in both the planar and the three-dimensional problems, is the approximation of the local solution expansion with its leading terms. Thus, for the 2-D problem of Figure 1, expansion (5) is approximated by its first N_{a} terms:
(12)
The above series can also be written, more simply, as
(13)
where V_{j} =r^{μj}cos(μ_{j}θ). Functions V_{j} are called singular functions and from now on, they are going to play an important role in the development of the formulation of the method. Note that these functions satisfy the governing equation and the boundary conditions along S_{A} and S_{B}.
Next step is to weight the governing equation by the singular functions V_{i}, in the Galerkin sense. Thus, we obtain the first set of discretized linear equations:
(14)
So, we have our first set of N_{a} discretized equations. Next, we apply Green’s second identity by considering that both the approximation of Φ and the singular functions V_{i} satisfy the Laplace equation and we obtain
(15)
The integrand of the first integral in (15) is equal to zero along boundary parts S_{A} and S_{B}, because there it is ∂Φ/∂n=0 and ∂V_{i}/∂n=0. Thus, the system of discretized equations (15) is simplified further to
(16)
For boundary part S_{C} the Dirichlet boundary condition is imposed by employing the Langrange multipliers function which is expanded in terms of quadratic basis functions M_{j} as follows:
(17)
In the above expansion, λ_{j} are the so-called Lagrange multipliers which are auxiliary parameters and they are additional unknowns in our problem. As we will see, they are calculated together with singular coefficients β_{j} after solving the complete system of equations. The quadratic basis functions M_{j}, in expansion (17), are defined as follows:
For j odd
(18)
and for j even
(19)
In expressions (18) and (19) h_{s} is the mesh-width, which is constant in this approach.
The normal derivative ∂V_{i}/∂n, which appears in (16), is expressed as follows
(20)
where ∇ is the normal unit vector on boundary S_{C} (see Figure 1). Furthermore, using (20) the normal derivative of V_{i} is given by
(21)
which is an easy expression to be used in calculations.
Before we proceed further, it is more convenient to change, in all the above formulation, variable s with variable θ by considering that θ=s/R (Figure 1). Now, there is one more step before we have our complete system of linear equations. So, according to the SFBIM the Dirichlet condition on S_{C} must be weighted with the basis functions M_{i}. Thus, finally, the discretized linear equations (16) take the form of the following system of N_{a}+N_{λ} discretized equations:
(22)
The integrands in (22) are not singular and integrations take place away from the singularity. If we substitute the approximation of Φ and function λ in (22), with their expressions in (13) and (17), respectively, then the system of equations (22) can further be written as
(23)
The above system of linear equations can also be expressed in matrix form as follows
(24)
or
K.b=c (25)
where K is called stiffness matrix. It is a symmetric matrix but it becomes singular when N_{a} < N_{λ}.
3.2. The SFBIM in the 3-D Case
In the SFBIM the solution of the problem is approximated by the leading terms of the local solution expansion given by (9). Thus, in (9) we employ the first N_{a} terms and we substitute the EFIFs with their polynomial expression given by (10). Then the approximate solution is written as follows:
(26)
where N (which also appears in (9)) is an additional parameter that allows us to ensure that the approximation of u satisfies the 3-D Laplace equation, by selecting it according to the restriction N_{p}<2N+1 (see Appendix). We note that, in principle, N could be taken to be infinity since after all, the sum would terminate after a finite number of terms due to the fact that we are differentiating a polynomial of degree N_{p}.
Following the notation used in previous applications of the SFBIM (e.g. [9-14]), the above expansion is written as follows:
(27)
In the above expression the coefficients a_{k,j} will be referred to as singular coefficients and they are the coefficients of the polynomial expression (10) of the EFIFs. The functions W_{k}^{(j)} are the singular functions (the same name was adopted for the 2-D approach) and in this case they have the form
(28)
It is easy to verify that ∇^{2}_{3-D} (W_{k}^{(j}^{)})=0 in domain Ω and that the W_{k}^{(j)}’s satisfy exactly the boundary conditions on S_{A} and S_{B} (see Appendix).
As in the 2-D case, we weight the governing equation by the singular functions W_{k}^{(j)} in the Galerkin sense. This gives the first N_{a}(N_{p}+1) discretized equations:
(29)
Recalling that W_{k}^{(j)} satisfies the governing equation and using Gauss’ divergence theorem we obtain
(30)
where ∂Ω=S_{A}∪S_{B}∪S_{C}∪S_{D}∪S_{E}. Since the singular functions W_{k}^{(j)} satisfy exactly the boundary conditions on S_{A} and S_{B} and considering the boundary conditions (8), the boundary integral in (30) is identically zero along S_{A} and S_{B}. The discretized equations (30) then become
(31)
As in the 2-D approach, the Dirichlet condition on the cylindrical boundary S_{C}, which is away from the singularity, is imposed by means of a Lagrange multiplier function λ, which is expanded in terms of basis functions M_{i}(θ,z):
on S_{C}, (32)
where N_{λ} is the number of the discrete Lagrange multipliers λ^{(i)} on S_{C}. The nodal values of λ appear as additional unknowns in the problem. The additional N_{λ} required equations are then obtained by weighting the Dirichlet boundary condition on S_{C} by the bilinear basis functions M_{i}(θ,z)^{}in the Galerkin sense. Thus, in the end, we obtain the following linear system of (N_{p}+1)N_{a}+N_{λ} discretized equations:
(33)
(34)
Integration in the neighborhood of the point singularities O and D must be avoided. Thus, the EFIFs should be of degree equal to one or zero, for this problem (i.e. N_{p}=1 or N_{p}=0), in order to force the surface integrals on S_{C} and S_{D} to vanish. Therefore, by choosing a value of N_{p} less or equal to 1 and by substituting the approximation of u and of the Lagrange multiplier function λ, with their expressions in (27) and (32), respectively, equations (33) and (34) take the form
(35)
(36)
It should be noted that, as in the 2-D case, the integrands in the above equations are non-singular and that all integrations are carried out far from the boundaries SA and SB and the points O and D causing the edge singularity. The surface integrals in (35) and (36) are two-dimensional integrals and are estimated using standard techniques, such as Gauss-Legendre quadrature. For example, if we use a 4-point rule the first of the integrals in equations (33) is estimated as follows:
where
Here, L_{z }^{(l)} and L_{θ }^{(l)} are the dimensions of each element (l) along directions z and θ, respectively and N_{E} is the number of boundary elements employed on the surface S_{C}.
The system of discretized linear equations, (35) and (36), can be written in block form as follows:
K.a_{λ}=C (37)
or
(38)
where the vector A contains the unknown coefficients a_{k,j} of the approximation of the EFIFs and the vector Λ contains the unknown (discrete) Lagrange multipliers λ^{(i)}. Clearly, the stiffness matrix K is symmetric. As in the case of 2-D problems, the number (Ν_{p }+1)N_{a} of the unknown coefficients a_{k,j} should be greater than or equal to the number of Lagrange multipliers N_{λ}, since the matrix Κ becomes singular when (Ν_{p }+1)N_{a}< N_{λ}.
4. Numerical Experiments
4.1. Numerical Results for the 2-D Problem
Expansion (17) indicates clearly that the Lagrange multiplier function λ, which is used to impose the Dirichlet boundary condition on S_{C}, has specific numerical values λ^{(j)} at positions j which are chosen according to a selected subdivision of S_{C} into boundary elements. For the 2-D model problem (Figure 1), boundary S_{C} is subdivided into N_{E} elements of equal size. Therefore, the number of Lagrange multipliers is N_{λ}=2N_{E}+1.
The integrals in (23) are calculated numerically by subdividing each quadratic element into 10 subintervals and using a 15-point Gauss-Legendre quadrature over each subinterval. Since the stiffness matrix is symmetric only the elements which are on or above its principal diagonal are calculated. Now, in order to obtain the optimum combination of N_{a} and N_{λ}, several series of runs have been made. Previous applications of the method in 2-D problems [9-14], indicated that N_{λ} should be large enough. This means that we must have an adequate number of elements on the boundary which will also help in having a better accuracy in numerical integration. However, N_{λ} should be smaller than or equal to N_{a} in order to avoid ill-conditioning of the stiffness matrix. In general, very high values of N_{a} are also avoided because although double precision is used in the codes, the computer accuracy cannot handle the contributions of the higher-order singular functions, which become very small for r < 1 or very large for r > 1. For the specific 2-D problem, the radius of the disk was chosen to be R=1. Hence, we do not have the phenomena described above. Thus, N_{λ} was varied from 3 to 15 and N_{a} was varied from 13 to 30. After performing calculations within this range of combinations of values for N_{a} and N_{λ}, it was observed that convergence occurs when N_{a}=N_{λ}=13 (optimum choice).
Table 1 presents the convergence in the values of the four leading coefficients, with respect to the number of Lagrange multipliers N_{λ} and for N_{a}=13. As in previous implementations of the method [9-14], one may observe that the values of singular coefficients converge rapidly with N_{λ}. In fact, in [28] a theoretical analysis of the method proved algebraic convergence in N_{λ}. Also, very accurate estimates are obtained. Clearly, since the exact values of the leading singular coefficients of expansion (5) are β_{1}=0, β_{2}=2/5=0.40000000, β_{3}=4/7≈0.571428571428571, β_{4}=0 and β_{j}=0, for j≥5, we can see, from Table 1, how fast convergence is achieved at N_{λ}=13.
The values of the leading singular coefficients calculated for N_{λ}=13 and various values of N_{a} are shown in Table 2. Extremely fast convergence with respect to N_{a} is observed (much faster than in previous implementations) and extremely accurate estimates are obtained. In fact, the leading singular coefficients converge from the very beginning. For N_{a}>13 the solution starts to deteriorate. This is because the size of the system increases and thus the number of numerical calculations is increased together with the numerical errors which, of course, are not avoidable at the calculation of the matrix elements or at the numerical inversion of the stiffness matrix of the model problem. For the optimum combination N_{a}=N_{λ}=13, the plot of λ as a function of variable θ, is shown in Figure 3. As expected, it has the form of a cosine trigonometric function because on S_{C} we have ∂V_{i}/∂r=μ_{i}r^{μi-1}cos(μ_{i}θ). In this graph the curve is smooth and free from oscillations.
In Table 3 the converged values of the singular coefficients, calculated for the optimal choices N_{a}=N_{λ}=13 and R=1, are presented. The values of this table indicate clearly that the contribution of the higher order terms, vanishes immediately for this 2-D model problem and that the approximate values are, practically, equal to the exact values.
Figure 4 shows the plot of the errors which appear in the calculated values of the leading singular coefficients β_{2} and β_{3}, as N_{λ} varies, when N_{a}=13. Here, the error is defined as the absolute value of the difference between the exact and the approximate solution of coefficients β_{j}:
│β_{j}^{(exact)} - β_{j}^{(approx)}│.
This is not the only definition of error. Any suitable norm can also be used to define the error in numerical calculations. However, in the present problem the absolute value of the difference between the exact value and the approximate value of β_{j} were much adequate in our present research.
N_{λ} | β_{1} | β_{2} | β_{3} | β_{4} |
3 | 0.0000000000000000 | 0.398200979875739 | 0.551151493230821 | 0.014748 |
5 | 0.0000000000000001 | 0.399978936129256 | 0.568905347448065 | -0.000307 |
7 | 0.0000000000000000 | 0.399998216881525 | 0.571283357528031 | 0.000000 |
9 | 0.0000000000000001 | 0.399999714617777 | 0.571400112413111 | 0.000000 |
11 | -0.0000000000000000 | 0.399999934075381 | 0.571420857864075 | 0.000000 |
13 | -0.0000000000000000 | 0.400000000000000 | 0.571428571428572 | -0.000000 |
N_{α} | β_{1} | β_{2} | β_{3} | β_{4} |
13 | -0.0000000000000000 | 0.400000000000000 | 0.571428571428572 | -0.000000 |
14 | -0.0000000000000000 | 0.399999980543728 | 0.571428571428572 | -0.000000 |
15 | 0.0000000000000000 | 0.399999980543728 | 0.571426036455448 | -0.000000 |
... | ............................. | .............................. | ............................. | .............. |
18 | 0.0000000000000000 | 0.399999980543728 | 0.571426036455448 | -0.000000 |
... | ............................. | .............................. | ............................. | .............. |
24 | 0.0000000000000000 | 0.399999980543414 | 0.571426024116464 | -0.000000 |
... | ............................... | .............................. | ............................. | ............. |
30 | 0.0000000000000003 | 0.399999980270283 | 0.571425960739584 | -0.000000 |
Both curves in Figure 4 indicate an algebraic convergence with N_{λ}. This behavior will be explained theoretically below.
j | β_{j} (approximate) | β_{j} (exact) |
1 | 0.000000000000000 | 0.000000000000000 |
2 | 0.400000000000000 | 0.400000000000000 |
3 | 0.571428571428572 | 0.571428571428571 |
4 | 0.000000000000000 | 0.000000000000000 |
5 | 0.000000000000000 | 0.000000000000000 |
... | ............................... | ............................... |
The above behavior of the present technique was also demonstrated in previous implementations of the method, in both 2-D and 3-D problems [9-14].
According to reference [28], if λ ∈ H^{k}(∂Ω_{AB}), for some k≥1 and λ_{h} is the approximation of the Lagrange multiplier function with h being always the mesh-width, then there exist positive constants C and γ∈(0,1), independent of N_{α} and h such that
(39)
where m=min{k, p+1} and H^{k}(∂Ω_{AB}), k∈N, is the usual Sobolev space on boundary part AB, which contains functions that have k generalized derivatives in the space of the square integrable functions L^{2}(Ω). Also, in [28] it was shown that
(40)
Inequalities (39) and (40) clearly indicate exponential convergence of the method with respect to the number of singular functions N_{a} and an algebraic convergence with the number of Lagrange multipliers N_{λ}. This fast convergence is visible in the results presented in Tables 1 and 2 and in the plot of Figure 4. In this study it is k=2 for the Sobolev space used. Also, quadratic basis functions were chosen (p=2). Now, if the two errors in (39) are forced to be equal to each other then we obtain
(41)
Hence using (41) and the "optimal" pair N_{a}=N_{λ}=13 we find that γ≈0.77≈0.80 ∈ (0,1). Therefore, with this value of γ and by prescribing parameter N_{a}, one can find N_{λ} from (41) and thus can have another pair of values. Using these values convergence of the method can be achieved.
4.2. Numerical Results for the 3-D Problem
In order to implement the SFBIM for the 3-D model problem, the boundary part S_{C} is subdivided into standard 3-D boundary elements (having dimensions along z and θ directions only) as shown in Figure 5. Specifically, N_{E}=N_{z}´N_{θ} elements are employed. Figure 5 shows the division of cylindrical boundary S_{C} into elements. As it is indicated by equations (35) and (36), calculations of the integrals, finally take place in a 2-D space (Figure 5). On the same figure the bilinear basis functions are also presented. The total number of Lagrange multipliers is N_{λ}=(N_{z}+1)´(N_{θ}+1). The surface integrals are estimated using a 9-point Gaussian quadrature rule over each element.
For all computations presented we take N_{p}=1 in (26), i.e. the EFIFs are approximated by polynomials of linear form. As explained earlier, this choice is made in order to force boundary integrals on S_{D} and S_{E} to vanish. In this way the integration in the vicinity of point singularities O and D is avoided. Now, with this choice of N_{p}, the parameter N in (26) is determined to be equal to 1; choosing a value N > 1 will not yield any additional terms in (26). Systematic runs have been carried out for the model problem in order to study the effects of N_{a} and N_{λ} and of the other parameters, on the numerical results.
The runs for different values of R included various combinations between N_{a} and N_{λ} in an attempt to find the "ideal’’ combination of these parameters, for which the method converges. Calculations were made with different combinations of N_{z} and N_{θ} but the pair of N_{z}=2 and N_{θ}=2 indicated (in combination with a suitable value of N_{a}) convergence of the method. Table 4 contains the values of the first two singular coefficients of A^{(a1)}(z) (i.e. of a_{1,1} and a_{1,2}) obtained with N_{λ}=9, R =0.01 and for different values of N_{a}. Note that since N_{p}+1=2, the step increment for N_{a}^{p}=(N_{p}+1)N_{a} in all the trials, is equal to 2. One may immediately observe very high rate of convergence with respect to N_{a} and great accuracy in the values obtained, even up to the 14^{th} decimal digit. Also, Table 5 shows the values of the first two EFIFs at the points z = 0.5 and 1, respectively, (i.e. of A^{(a1)}(0.5) and A^{(a2)}(1.0)) calculated for the same values of N_{λ} and R and for various values of N_{a}.
The values of z were selected so that a comparison can be made between the SFBIM and the energy projection method presented in [30]. We may observe, again, very fast convergence with respect to N_{a} and high accuracy in the values obtained. In both Tables 4 and 5 convergence corresponds to the "optimal" combination of values N_{a}=7, N_{a}^{p}=(N_{p}+1)N_{a}=14 and N_{λ} = 9. The value of the radius of the domain, for our 3-D model problem, is taken as R = 0.01. A model problem of the same geometry and the same dimensions was also adopted in [30], where a special finite element technique was implemented in order to solve it.
Table 6 contains the converged approximate values of the singular coefficients obtained for R = 0.01 and for the "optimal" combination of N_{a} and N_{λ}. One may observe that the converged values of the coefficients, obtained with the SFBIM, are practically the same with those contained in the exact form of the EFIFs which is known for this model problem as already explained above.
a_{1,1} | a_{1,2} | |
10 | 0.99999999999993294 | 0.50000000000005129 |
12 | 1.00000000000000020 | 0.50000000000000766 |
14 | 1.00000000000000004 | 0.50000000000000001 |
16 | 0.99999999999998842 | 0.50000000000014344 |
18 | 0.99999999998523833 | 0.50000000001876197 |
10 | 1.25000000000006 | 1.50000000000008 |
12 | 1.25000000000000 | 1.50000000000001 |
14 | 1.25000000000000 | 1.50000000000000 |
16 | 1.25000000000106 | 1.50000000000113 |
18 | 1.25000000001046 | 1.50000000001140 |
Finally, Table 7 compares the values obtained by the SFBIM for R=0.01 and z=1.0 with the results given by the energy projection method [30], for the same values of R and z. The exact solution, for our 3-D model problem, is also presented on the same table. Clearly, the SFBIM gives results which are comparable with the exact values.
In fact, for the accuracy considered in this work, the approximate values obtained with the SFBIM, coincide with the exact values. Also, it is obvious that the SFBIM, is significantly more accurate than the finite element technique which is implemented in [30]. This is because the SFBIM is a direct method and no post-processing is required. Also, the numerical error and the CPU time in the SFBIM, are much smaller than those observed in other numerical techniques, for this type of boundary value problems in Applied Mathematics.
i | k | a_{i,k} |
1 | 1 | 1.0000000000000004 |
1 | 2 | 0.5000000000000001 |
2 | 1 | 1.0000000000000003 |
2 | 2 | 0.5000000000000048 |
... | ... | ................................... |
... | ... | ................................... |
i | |||||
1 | 1.500 | 0.000 | 1.499 | 0.001 | 1.500 |
2 | 1.500 | 0.000 | 1.500 | 0.000 | 1.500 |
In Figure 6 there is a graph presenting with "squares" the approximation of the first EFIF A^{(a1)}(z)=A_{1}(z) obtained at convergence, for which it is N_{a}^{p} =N_{a}(N_{p}+1)=14 and N_{λ}=9. On the same graph and with a solid line, the analytic form of the first EFIF A^{(a1)}(z) is also presented. One may observe that, practically, the approximation of the polynomial A^{(a1)}(z) coincides with the exact solution. Finally, the graph of Figure 7 shows how the error changes its value as N_{a}^{p} varies. In this case, the error is defined by using the infinity norm as follows:
.
At N_{a}^{p} =N_{a}(N_{p}+1)=14, the error is minimized and it is of the order of 10^{-16}. This value of error corresponds to the optimum combination of the parameters (i.e. for N_{a}=7, N_{p}=1, N=1 and N_{λ}=9). After a significant number of runs it was found that this is the minimum of the Ε vs. N_{a}^{p} graph for all values of N_{λ}. A general conclusion derived from these graphs is that the method requires a small number of boundary elements and that it exhibits very fast convergence.
5. Conclusions
The Singular Function Boundary Integral Method (SFBIM) has been formulated for both a 2-D elliptic boundary value problem with a point boundary singularity and a 3-D Laplace equation problem with a straight-edge boundary singularity. With this method the singular coefficients and the Edge Flux Inensity Functions (EFIFs), are the primary unknowns in the 2-D and the 3-D cases, respectively. The latter are approximated by polynomials, the coefficients of which are unknowns in the formulation.
Both the singular coefficients, of the two-dimensional local solution and the coefficients of the EFIFs, in the 3-D case, are calculated directly and not by post-processing of the solution. The implementation of the SFBIM to both a 2-D and a 3-D Laplacian model problems, yielded highly accurate results for the singular coefficients and the EFIFs and exhibited fast convergence, as in other two-dimensional applications [9-12] of the method. For the planar problem, the numerical results are favorably compared with the analytic solution. The numerical results for the 3-D case compare favorably with both the exact solution and those of the energy projection method [30].
Currently a particular version of the method is being investigated for other three-dimensional Laplace equation problems. This approach will aim at the better treatment of the inner sum in Eq. (26). Also, a hybrid method is being developed in which the advantages of both the SFBIM and the finite element techniques are being exploited for problems in Mechanics with boundary singularities.
Appendix
Herein we will show that ∇^{ }^{2}(W_{ }^{j}_{k})=0 and thus that the approximation of u also satisfies the Laplace equation. So, we consider, first, expression (28). Then, the explicit form of W_{ }^{j}_{k }is as follows:
So, we have
where
and
Also,
and
Therefore, after substituting in the expression of ∇^{2}(W_{ }^{j}_{k}) the second order partial derivatives of W_{ }^{j}_{k},_{ }with the above expansions and performing several algebraic operations, we obtain
It can be easily verified that all coefficients of the derivatives ∂^{2i}_{z }(z^{j-1}), in the first brackets, cancel each other. Also, by considering the definition of parameter c_{ki}, the expression of ∇^{2}(W_{ }^{j}_{k}) takes the form
Obviously, the algebraic operations inside the parentheses give us a zero result. Thus, the expression of ∇^{2 }(W_{ }^{j}_{k}) shrinks to
which, of course, leads us to ∇^{2}(W_{ }^{j}_{k}) = 0 if j-1 £ 2N+1 and since the maximum value of j is N_{p}+1 we may say that W_{ }^{j}_{k }satisfy the 3-D Laplacian equation with the restriction N_{p }£ 2N+1.
It is also easily shown that singular functions W_{ }^{j}_{k }satisfy the boundary conditions on S_{A} and S_{B}. Indeed, we have
and
References