A family of exponentially-gradient elements for numerical computation of singular boundary value problems

Engineering Analysis with Boundary Elements 80 (2017) 184–198

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

A family of exponentially-gradient elements for numerical computation of singular boundary value problems

MARK

⁎

Morteza Eskandari-Ghadia, , Delaram Mehdizadeha, Ali Morshedifardb, Mohammad Rahimiana a b

School of Civil Engineering, College of Engineering, University of Tehran, P.O Box 11155-4563, Tehran, Iran Department of Civil and Environmental Engineering, University of California at Irvine, Irvine, CA 92617, USA

A R T I C L E I N F O

A BS T RAC T

Keywords: Exponentially gradient elements Singularity High gradient Mixed boundary value problem Boundary element method Finite element method

In this paper, a new family of single-parameter exponentially gradient elements (EG-elements) are introduced, which can be used in various numerical procedures such as boundary and finite element methods. These elements have the ability to accurately interpolate the unknown values in regions, where either high gradient or singularity of the unknown field occurs. The shape functions of two-dimensional EG-elements are high gradient at either a corner of the element or at an edge of the element. Another advantage of this element is that the regular quadratic shape functions are obtained as a special case of EG-elements by adjusting the single parameter of the element, which allows this element to be used as regular Lagrange quadratic element, where it is appropriate. Some mixed boundary value problems are solved with the use of EG-elements in a boundary element program to show the capability of these elements for capturing the solution with less number of elements and higher accuracy.

1. Introduction Many kinds of Mixed Boundary Value Problems (MBVP) are encountered in engineering mechanics [1–5] and potential theory [6,7], where the natural boundary conditions are discontinuous and thus the derivatives of the solution for the MBVP is singular at discontinuities. Contact problems, crack investigation [7], scattering problems and problems in fracture mechanics are some examples in mechanics, and electrostatic potential of plates charged to a prescribed potential or electrostatic field due to two parallel plates are examples in potential theory [6]. When an analytical solution is sought for this type of MBVP, a very precise attention should be paid, and when it is considered to be solved numerically, some difficulties are encountered both in reproducing the singularities of the problem and in numerical integrations [7,8]. Because of these issues, mathematicians and engineers are interested in this kind of problems. In addition, the numerical research group is also interested in the phenomenon. There exist a few researches for either numerical integration of singular functions or numerical solutions for singular integral equations [7–11], many of which proposed to express the unknown function in terms of linear combination of some known functions, such as Chebyshev [9] and Jacobi [12,13] series functions. When one expresses an unknown function, say f , in the form of linear combination of some other known functions in its entire domain, where both regular and singular

⁎

behavior of f are seen, it should be noticed that a large number of the known functions should be used to capture all the behavior of f . On the other hand, if f is expressed in terms of some function in a subset of its domain, then a small number of the known functions may be adequate. Polynomials from constant to higher degree such as 3rd or 4th degree are used in ordinary finite element method, called Lagrange elements. Since, these elements due to their shape functions are naturally smooth, they cannot capture the singular or high gradient behavior present in the solution of various MBVPs appropriately and a large number of elements are needed to have a poor approximation for the solution. Pak and Ashlock [14] introduced a family of twoparameter power based adaptive elements which can be adjusted to capture both regular and high-gradient variations of different functions. Although the shape functions they derived can be used for numerical analysis of some high-gradient functions, however, a combination of two real parameters should be selected for adjusting the behavior of the element. On the other hand, the complicated form of the kernel function used in deriving the shape functions, makes them less attractive for applications. In this work, we are going to introduce a new one-parameter family of shape functions with an adaptive behavior to capture both regular and singular behavior of different phenomena. Since, integration of singular functions needs some special attention, nonsingular, however high gradient shape functions are used in this new element. To this

Corresponding author. E-mail address: [email protected] (M. Eskandari-Ghadi).

http://dx.doi.org/10.1016/j.enganabound.2017.03.013 Received 4 November 2016; Received in revised form 20 January 2017; Accepted 9 March 2017 0955-7997/ © 2017 Elsevier Ltd. All rights reserved.

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

end, an exponential function is used as the kernel function of the new element. Since, both the exponential function and its derivatives are linearly dependent, some extra terms are added to the shape functions to make a more complete basis, in the sense of vector space [see [15,16]], for the new element. On the other hand, since the main part of the gradient of the new elements is expressed by the exponential function, we call them the Exponentially Gradient-elements (EGelements). The EG-elements in this paper are introduced in such a way that their behavior from regular to very high gradient can be controlled with only one parameter, and the so called Lagrange element can be derived as a special case of the EG family of elements. In this way, any EG-element introduced in this paper will be an adaptive element. Some two- and three-dimensional contact boundary value problems are numerically solved with the application of EG family of elements in a boundary element method to show the power of EGfamily line and surface elements.

Fig. 1. The EG-kernel for different values of m. s(s + 1)

function f (s; m ). We call the term 2 em(s −1) in the interpolation function f (s; m ) as the EG-kernel function, which can be localized and sharpened toward the EG end of the parent domain via the m parameter, as can be seen from Fig. 1. As is clear from the formulation of the EG-kernel function and as seen from Fig. 1, the EG-kernel function is smooth everywhere for different values of m , however, the larger the value of m the higher the gradient of the kernel function is achieved. If the usual nodal requirements for an element with three nodes are imposed, the following equations are obtained:

2. Basis for construction of EG-family elements A linear vector (function) space is a set containing some vectors (functions) that except some other properties satisfy two main axioms, which are closure under addition and closure under multiplication by real numbers [see [15,16]]. A Finite Dimensional Linear Vector Space (FDLVS) may be a subset of a linear vector (function) space that is spanned on a basis with finite members. An arbitrary vector, say u, can be exactly expressed in an FDLVS if there exists a basis for the FDLVS containing u, otherwise, u cannot be exactly expressed as a linear combination of the members of the FDLVS, and thus any linear combination of the members of FDLVS is an approximation for u. Any smooth function can be expressed in an FDLVS that spanned on a basis containing smooth functions, however, a high gradient or singular function cannot be well expressed on the domain of interest in an FDLVS that is spanned on a basis containing only smooth functions with low gradient behavior. In this case, high gradient functions are needed to make a good approximation for the singular functions. There are many cases in computational mechanics that high gradient or singular functions are encountered. Thus, an FDLVS spanned on a basis containing high gradient functions is needed to handle the singular functions. On the other hand, the functions in the basis of an FDLVS should be continuous on the domain of interest to be a place for making good approximations for continuous functions. In this paper, we are going to present an FDLVS spanned on a basis containing only three functions with an arbitrary gradient, where the gradient is adjusted with only one parameter. This FDLVS is used for making a new three-node finite element used in one-dimensional boundary value problems. With the use of standard methods, the two-dimensional high gradient elements are made. The Lagrange three-node one-dimensional element will be a special case of the proposed element.

f ( − 1) = A − B = f1 f (0) = A = f2 f (1) = A + B + C = f3

where f1, f2 and f3 are the nodal values of the variable being interpolated. A simultaneous solution of these three equations gives the constants A, B and C in terms of f1, f2 and f3 as follows:

A = f2 ,

C = f1 − 2f2 + f3

(3)

1 s(−2 + em(s −1)(s + 1)) 2 = 1 + s − em(s −1)s(s + 1) 1 = em(s −1)s(s + 1) 2

N1EG = N2EG N3EG

(4)

The EG-shape functions should be selected in such a way that they can make a basis for an FDLVS containing completely regular to very high gradient functions defined on a finite subset of real numbers. In this way, one can use linear combinations of those shape functions to describe any regular to very high gradient functions defined on the subset. Of course, the shape functions should possess the delta function property, the partition of unity property, C°-continuity, and consistency [17] conditions. The delta function property implies that the shape functions should have unit values at their home node and vanish at remote nodes of the element, while the partition of unity property ensures that the shape functions should sum to unity. The C°continuity at the subset makes sure of the continuity of any function directly expressed by the shape functions, and consistency condition requires that the interpolation function represents exactly any polynomial function up to some order. Because of analyticity of the exponential function, the shape functions used in the EG-elements are smooth with any value of m , however, their gradient is controlled by the value of m . Examples of the shape functions are shown in Fig. 2 for a wide range of the parameter m . As it can be seen from the display, they exhibit exactly the kind of localizability and adaptability that are lacking in current shape function ensembles. An important property that can be observed is that the shape functions of the one-dimensional Lagrange quadratic element can be easily reproduced by choosing

A family of single parameter Exponentially-Gradient (EG) shape functions is introduced, which can be utilized for making one- to threedimensional finite and boundary elements. The shape functions are adjusted in such a way that they can vary from regular to very high gradient behavior. To do so, first the following interpolation function is presented for derivation of the shape functions of a three-node onedimensional element:

s(s + 1) m(s −1) e 2

B = f2 − f1 ,

Substitution of Eq. (3) into Eq. (1) makes it possible to rearrange the interpolated function in the form of f (s; m )=N1(s; m )f1 + N2(s; m )f2 + N3(s; m )f3 where Ni(s ), (i = 1, 2, 3) are the EG-shape functions for a 3-node one-dimensional EG-element obtained as

3. One-dimensional EG-elements

f (s; m ) = A + Bs + C

(2)

(1)

where s is the parent coordinate in the domain [ − 1, 1] (see Fig. 1), A, B and C are constants, which will be derived from the nodal values of the variable being interpolated and m is the sole EG-parameter that can be chosen in such a way that a desired behavior is obtained from the

185

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Fig. 2. EG-shape functions for various values of m.

m = 0 in EG-shape functions, an analytic continuity and inclusiveness that is precious, as it can, for example, be utilized to facilitate smoother element class transitions and enhance solution convergence even in a uniform mesh setting. It is observed from both the kernel function (Fig. 1) and shape functions (Fig. 2) that as m increases, we have a higher gradient at the near-singular side of the element, which is assumed to be at s = 1. A wide range of problems with singularities can be handled by judiciously choosing the EG-parameter for these elements. To express the capability of the EG-shape functions for capturing different gradients, we interpolate the function q 1 p for

Fig. 3. Comparing the function

1 q

(1 − x 2) p

for

p q

1 2

= ,

1 3

and

2 3

(solid lines) with their

interpolation functions using 4 quadratic elements (dotted lines), and 3 quadratic elements plus 1 EG-element with m = 10 (dashed lines). The rightmost node is fixed at x = 0.999 .

(1 − x 2)

0 < x,

p q

< 1 with both quadratic Lagrange elements and EG-elements.

Fig. 3 illustrates a comparison among the function p q

1

1

1 q

(1 − x 2)

p

for

2

= 2 , 3 and 3 with both quadratic Lagrange elements and EGelements with m = 10 . As seen, the quadratic Lagrange elements can be poorly fitted to the functions especially near the edge x = 1, while the EG-elements are very close to the functions. It is easy to show that p the EG-elements can capture the behavior of q 1 p with q > 1 for

Fig. 4. A nine-node quadrilateral EG-edge element. The high gradient edge is adjusted at s = 1.

(1 − x 2)

0 < x < 1 very well. 4.1. Nine-node quadrilateral EG-edge element 4. Two-dimensional EG-elements

Consider the nine-node quadrilateral element in parent s − t plane illustrated in Fig. 4, where the edge s = 1 is a place, where the interpolated function is high gradient. The shape functions may be constructed from tensor products of a vector containing the three-node one-dimensional EG-shape functions introduced in the previous section in s− direction and a vector containing the one-dimensional

We shall now employ the above results in deriving two-dimensional EG-element's shape functions that can handle high-gradient phenomena near one edge or a corner of the element, which are respectively called here as EG-edge element or EG-corner element.

186

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

m in s− and n in t− direction, can be used, however, the graphical results are presented for identical values of m and n in both directions. The 9-node EG-edge and corner elements all satisfy the corresponding smoothness, interpolation, localization, and completeness conditions as outlined earlier. The above nine-node edge and corner element shape functions are shown in Figs. 5 and 6, respectively.

Lagrange quadratic interpolation functions in t− direction, which themselves are the degeneration of the shape functions given in Eq. (4) for m = 0 except that s should be replaced by t . Thus, shape functions for a nine-node quadrilateral EG-edge element are given as follows

N1EEG −9(s, t ; m ) = N1EG (s; m )N1L (t ) = N2EEG −9(s, t ; m ) = N3EG (s; m )N1L (t ) = N3EEG −9(s, t ; m ) = N3EG (s; m )N3L (t ) = N4EEG −9(s, t ; m ) = N1EG (s; m )N3L (t ) = N5EEG −9(s, t ; m ) = N2EG (s; m )N1L (t ) = N6EEG −9(s, t ; m ) = N3EG (s; m )N2L (t ) = N7EEG −9(s, t ; m ) = N2EG (s; m )N3L (t ) = N8EEG −9(s, t ; m ) = N1EG (s; m )N2L (t ) = N9EEG −9(s, t ; m ) = N2EG (s; m )N2L (t ) =

1 s(−2 + em(s −1)(s + 1))(t − 1)t 4 1 m(s −1) e s(s + 1)(t − 1)t 4 1 m(s −1) e s(s + 1)t (t + 1) 4 1 s(−2 + em(s −1)(s + 1))t (t + 1) 4 1 (1 + s − em(s −1)s(s + 1))(t − 1)t 2 1 m(s −1) e s(s + 1)(1 − t 2 ) 2 1 (1 + s − em(s −1)s(s + 1))t (t + 1) 2 1 s(−2 + em(s −1)(s + 1))(1 − t 2 ) 2 (1 + s − em(s −1)s(s + 1))(1 − t 2 )

4.3. Eight-node quadrilateral EG-edge element To obtain an eight-node serendipity-type quadrilateral EG-element that can be computationally attractive, one may start from the ninenode element and employ a node-ejection technique [14]. As an example, one may subtract (1/4)N9 from shape functions N1 through N4 and add (1/2)N9 to previous shape functions N5 through N8. Since, the sum of the nine-node shape functions introduced in the previous section equals unity, the eight-node shape functions so constructed will also sum to unity. In addition, the other properties will remain unchanged. This approach is equivalent to distributing N9 to the other eight nodes. The eight-node EG-edge element shape functions so derived are

(5)

N1EEG −8(s, t ; m ) = N1EG (s; m )N1L (t ) − 1/4N9EEG −9(s, t ; m ) 1 = (t − 1)(1 + s − em(s −1)s(s + 1) + t − st ) 4 N2EEG −8(s, t ; m ) = N3EG (s; m )N1L (t ) − 1/4N9EEG −9(s, t ; m ) 1 = − (s + 1)(−1 + em(s −1)s − t )(t − 1) 4 N3EEG −8(s, t ; m ) = N3EG (s; m )N3L (t ) − 1/4N9EEG −9(s, t ; m ) 1 = (s + 1)(t + 1)(−1 + em(s −1)s + t ) 4 N4EEG −8(s, t ; m ) = N1EG (s; m )N3L (t ) − 1/4N9EEG −9(s, t ; m ) 1 = (t + 1)(−1 + em(s −1)s(s + 1) + t − s(t + 1)) 4 N5EEG −8(s, t ; m ) = N2EG (s; m )N1L (t ) + 1/2N9EEG −9(s, t ; m ) 1 = (s + 1)(−1 + em(s −1)s )(t − 1) 2 = N3EG (s; m )N2L (t ) + 1/2N9EEG −9(s, t ; m ) N6EEG −8(s, t ) 1 = − (s + 1)(t 2 − 1) 2 N7EEG −8(s, t ; m ) = N2EG (s; m )N3L (t ) + 1/2N9EEG −9(s, t ; m ) 1 = − (s + 1)(−1 + em(s −1)s )(t + 1) 2

4.2. Nine-node quadrilateral EG-corner element Similarly, the shape functions for a nine-node quadrilateral EGcorner element with a high gradient behavior at the corner of the edges s = 1 and t = 1 may be constructed from a tensor product of two similar vectors of the three-node one-dimensional EG-shape functions. Thus, they are given as

N1CEG −9(s, t ; m, n ) = N1EG (s; m )N1EG (t ; n ) 1 = s(−2 + em(s −1)(s + 1))t (−2 + en(t −1)(t + 1)) 4 N2CEG −9(s, t ; m, n ) = N3EG (s; m )N1EG (t ; n ) 1 = em(s −1)s(s + 1)t (−2 + en(t −1)(t + 1)) 4 N3CEG −9(s, t ; m, n ) = N3EG (s; m )N3EG (t ; n ) 1 = em(s −1)+ n(t −1)s(s + 1)t (t + 1) 4 N4CEG −9(s, t ; m, n ) = N1EG (s; m )N3EG (t ; n ) 1 = en(t −1)s(−2 + em(s −1)(s + 1))t (t + 1) 4 N5CEG −9(s, t ; m, n ) = N2EG (s; m )N1EG (t ; n ) 1 = (1 + s − em(s −1)s(s + 1))t (−2 + en(t −1)(t + 1)) 2 N6CEG −9(s, t ; m, n ) = N3EG (s; m )N2EG (t ; n ) 1 = em(s −1)s(s + 1)(1 + t − en(t −1)t (t + 1)) 2 N7CEG −9(s, t ; m, n ) = N2EG (s; m )N3EG (t ; n ) 1 = en(t −1)(1 + s − em(s −1)s(s + 1))t (t + 1) 2 N8CEG −9(s, t ; m, n ) = N1EG (s; m )N2EG (t ; n ) 1 = s(−2 + em(s −1)(s + 1))(1 + t − en(t −1)t (t + 1)) 2 N9CEG −9(s, t ; m, n ) = N2EG (s; m )N2EG (t ; n ) = (1 + s − em(s −1)s(s + 1))(1 + t − en(t −1)t (t + 1))

N8EEG −8(s, t )

= N1EG (s; m )N2L (t ) + 1/2N9EEG −9(s, t ; m ) 1 = (s − 1)(t 2 − 1) 2

(7)

These shape functions are illustrated in Fig. 7 for m = 10 . It should be noted that we have the node numbering convention for eight-node quadrilateral EG-edge element in Fig. 4 by omitting node 9.

4.4. Eight-node quadrilateral EG-corner element With the same procedure of node-ejection technique, the shape functions for EG-corner elements are given as

(6)

As indicated in the relations (6), different EG-parameters, namely

187

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

N1CEG −8(s , t ; m, n ) = N1EG (s ; m )N1EG (t ; n ) − 1/4N9CEG −9 (s , t ; m, n ) 1 = (−1 − t − en(t −1)(s − 1)t (t + 1) 4

N1EEG −6(s′, t′) = (s′ + t′ − 1)(2s′ + 2t′ − 1) ⎛ ⎞ − 2mt′ N2EEG −6(s′, t′; m ) = s′⎜⎜ −1 + s′ + e s′+ t′ (s′ − t′) + t′⎟⎟ ⎝ ⎠

+ s(−1 − em(s −1)(s + 1)(t − 1) + 3t )) N2CEG −8(s ,

t ; m, n ) = =

N3CEG −8(s ,

t ; m, n ) = =

N4CEG −8(s ,

t ; m, n ) = =

N5CEG −8(s , t ; m )

N6CEG −8(s , t ; n )

N7CEG −8(s , t ; m )

N8CEG −8(s , t ; n )

N3EG (s ;

m )N1EG (t ; n ) − 1/4N9CEG −9(s , t ; m, n )

N3EEG −6(s′, t′; m ) = e

1 − (s + 1)(1 + em(s −1)s(t − 1) + t − en(t −1)t (t + 1)) 4 N3EG (s ; m )N3EG (t ; n ) − 1/4N9CEG −9(s , t ; m, n ) 1 (s + 1)(t + 1)(−1 + em(s −1)s + en(t −1)t ) 4 N1EG (s ; m )N3EG (t ; n ) − 1/4N9CEG −9 (s , t ; m, n ) ⎛ ⎞ 1 (t + 1) ⎜ (s + 1)( −1 + em(s −1)s ) − en(t −1)(s − 1)t ⎟ ⎝ ⎠ 4

= N2EG (s ; m )N1EG (t ; n ) + 1/2N9CEG −9(s , 1 = (s + 1)(−1 + em(s −1)s )(t − 1) 2 = N3EG (s ; m )N2EG (t ; n ) + 1/2N9CEG −9(s , 1 = − (s + 1)(t + 1)(−1 + en(t −1)t ) 2 = N2EG (s ; m )N3EG (t ; n ) + 1/2N9CEG −9(s , 1 = − (s + 1)(−1 + em(s −1)s )(t + 1) 2 = N1EG (s ; m )N2EG (t ; n ) + 1/2N9CEG −9(s , 1 = (s − 1)(t + 1)(−1 + en(t −1)t ) 2

− 2mt′ s′+ t′

s′(s′ − t′) − t′ − (s′ − 2t′)(s′ + t′)

N4EEG −6(s′,

t′) = − 4s′(s′ + t′ − 1) ⎛ ⎞ − 2mt′ N5EEG −6(s′, t′; m ) = 2s′⎜⎜s′ + t′ + e s′+ t′(− s′ + t′)⎟⎟ ⎝ ⎠ N6EEG −6(s′, t′) = − 4t′(s′ + t′ − 1)

(10)

These shape functions are illustrated in Fig. 10 for m = 10 .

t ; m, n )

4.6. Six-node triangular EG-corner element

t ; m, n )

Similarly, the shape functions for a six-node triangular EG-corner element with a high gradient behavior at the corner of the edges s = 1 and t = 1 may be constructed as follows

t ; m, n )

N1CEG −6(s, t ; n )

= N1CEG −8(s, t ; m, n ) + N2CEG −8(s, t ; m, n ) + N5CEG −8(s, t ; m ) 1 = t (−2 + en(t −1)(t + 1)) 2

t ; m, n )

N2CEG −6(s, t ; m, n ) = N3CEG −8(s, t ; m, n ) + N9CEG −9(s, t ; m, n )/8 1 = (s + 1)(t + 1)(−1 + em(s −1)s 8

(8)

+ en(t −1)(1 + em(s −1)s )t )

Examples of the foregoing set of eight-node corner EG-element shape functions with sharp edge gradients are shown in Fig. 8 for m = n = 10 .

N3CEG −6(s,

4.5. Six-node triangular EG-edge element

N4CEG −6(s,

t ; m, n ) = N4CEG −8(s, t ; m, n ) + N9CEG −9(s, t ; m, n )/8 1 = (t + 1)((s + 1)(−1 + em(s −1)s ) 8 + en(t −1)(1 − 3s + em(s −1)s(s + 1))t )

= N6CEG −8(s, t ; n ) 1 = − (s + 1)(t + 1)(−1 + en(t −1)t ) 2 N5CEG −6(s, t ; m, n ) = N7CEG −8(s, t ; m ) − 2N9CEG −9(s, t ; m, n )/8 1 = − (s + 1)(−1 + em(s −1)s )(t + 1)(1 + en(t −1)t ) 4 1 CEG −6 CEG −8 (s , t ; n ) = N8 (s, t ; n ) = (s − 1)(t + 1)(−1 + en(t −1)t ) N6 2 (11)

With reference to Fig. 9a, the shape functions of a six-node triangular EG-element may be obtained by collapsing one side of an eight-node square EG-element on a node and modifying the shape functions using the procedure described in [18]. The six-node EG-edge element shape functions so derived are

N1EEG −6(s, t ) = N1EEG −8(s, t ; m ) + N2EEG −8(s, t ; m ) 1 + N5EEG −8 (s, t ; m ) = t (t − 1) 2 N2EEG −6(s, t ; m ) = N3EEG −8(s, t ; m ) + N9EEG −9(s, t ; m )/8 1 = (s + 1)(t + 1)( − 1 + t + em(s−1)s(t + 1)) 8

t; n)

The shape functions in (s′, t′) natural coordinate system are presented as

N1CEG −6(s′, t′; n ) = (2s′ + 2t′ − 1)(−1 + e2n′(s′+ t′−1)(s′ + t′)) 2mt′

N3EEG −6(s, t ; m ) = N4EEG −8(s, t ; m ) + N9EEG −9(s, t ; m )/8 1 = (t + 1)( − 1 + t + s( − 1 − 3t + em(s−1)(s + 1)(t + 1))) 8 1 N4EEG −6(s, t ) = N6EEG −8(s, t ) = − (s + 1)(t 2 − 1) 2 N5EEG −6(s, t ; m ) = N7EEG −8(s, t ; m ) − 2N9EEG −9(s, t ; m )/8 1 = − (s + 1) ( − 1 + em(s−1)s )(t + 1)2 4 1 N6EEG −6(s, t ) = N8EEG −8(s, t ) = (s − 1)(t 2 − 1) (9) 2

N2CEG −6(s′, t′; m, n ) =

If we consider a map in the form of s = (s′ − t′)/(s′ + t′) and t = 2(s′ + t′) − 1 between the collapsed coordinates (s, t ) and the natural coordinates (s′, t′) shown in Fig. 9b, the shape functions can be derived in terms of (s′, t′) natural coordinate system with the high gradient edge changing from s = 1 to t′ = 0 . The shape functions so derived are

N6CEG −6(s′,

e− s′+ t′ s′(s′ − t′)(1 + e2n(−1+ s′+ t′)(2s′ + 2t′ − 1)) 2(s′ + t′) 1

+ 2 s′(−1 + e2n(s′+ t′−1)(2s′ + 2t′ − 1)) 1

N3CEG −6(s′, t′; n ) = 2 (−s′ − e2n(s′+ t′−1)(s′ − 2t′)(2s′ + 2t′ − 1)) − 2mt′

e s′+ t′

+ 2(s′ + t′) s′(s′ − t′ + e2n(s′+ t′−1)(s′ − t′)(2s′ + 2t′ − 1)) N4CEG −6 (s′, t′; n ) = − 2s′(−1 + e2n(s′+ t′−1)(2s′ + 2t′ − 1)) N5CEG −6(s′, t′; m, n ) =

e s′+ t′ s ⎛ ⎜ −s′ s′ + t′ ⎝ − 2mt′

⎞ 2mt′ + t′ + e s′+ t′ (s′ + t′)⎟ ⎠

× (1 + e2n(s′+ t′−1)(2s′ + 2t′ − 1)) t′; n ) = − 2t′(−1 + e2n(s′+ t′−1)(2s′ + 2t′ − 1)) (12) Examples of the foregoing set of six-node corner EG-element shape functions with high gradient behavior at the corner of the edges t′ = 0 and s′ = 1 − t′ are shown in Fig. 11 for m = n = 10 . 188

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Fig. 5. Nine node EG-edge element shape functions for m = 10 . N

5. Applications in boundary element method

cu +

j =1

To examine the performance of EG-elements introduced in this paper, it is useful to consider their implementation within the framework of boundary element method, where singular behavior exists. We are going to demonstrate the performance of both one- and twodimensional EG-elements by implementing them in the boundary element procedure. If we show the displacements at a field or a boundary point ξ by uj (ξ ) due to the tractions pj (x) applied at some point x , then the displacements at any point ξ ⊂ Γ may be given by the following integral equation [19]:

cij (ξ )uj (ξ ) =

∫Γ uij*(ξ; x)pj (x)dΓ (x) − ∫Γ pij*(ξ; x)uj(x)dΓ (x)

⎛

L

∑ ⎜⎜∑ ⎝ l =1

⎞ J l wl ( p* NT(s, t ))l ⎟⎟un = ⎠

N

⎛

L

∑ ⎜⎜∑ j =1

⎝ l =1

⎞ J l wl (u*NT(s, t ))l ⎟⎟pn ⎠ (14)

where L is the number of Gauss points on the surface element, J l is the determinant of the Jacobi matrix for the surface elements and wl is the weighting coefficients at Gauss points. In what follows, proposed EG-elements are applied to a set of mixed boundary value problems to evaluate their numerical performance. For the boundary value problem of a rigid footing resting on the surface of a half-space under some prescribed displacements, there are singularities for tractions at the edges of the footing [6]. The three examples considered in this paper address this problem and show how utilizing EG edge and corner elements help capture this singular behavior. The first example considers the mixed boundary value problem [6] due to contact of a rigid circular disc with an elastic half-space under vertical movement, which is an axisymmetric case. With reference to a cylindrical coordinate system (r , θ , z ) attached at the center of the disc, the only variable on z = 0 plane is the distance from the origin, r , and the boundary elements are on the one-dimensional radial line. The second example considers the plane strain problem of a rigid strip foundation resting on a half-space under some prescribed vertical displacement. The one-dimensional boundary elements in this case are along the width of the strip footing. The three-dimensional mixed boundary value problem of a rigid square footing in contact with an isotropic half-space is considered as the third example, where two-

(13)

where pij* (ξ; x) and uij*(ξ; x) are respectively the traction and the displacement Green's functions, cij (ξ ) is the smoothness function of the surface point ξ and pj (x) is the traction at the field point x . Eq. (13) is valid for both three- and two-dimensional boundary value problems. Dividing the boundary Γ into N surface elements, namely Γj ( j = 1, 2, . . . , N ), writing the integrals in terms of the parent coordinate system (s, t ), and employing the standard Gaussian quadrature for the integration, Eq. (13) changes to a set of linear algebraic equations as:

189

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Fig. 6. Nine-node EG-corner element shape functions for m = n = 10 .

dimensional boundary elements should be utilized in boundary element method. Before going to present the examples, since the shape functions are no longer polynomials, it is needed to explain the method of integration. The integrals that we have to evaluate are in the form of 3

I=

∑∫ i =1

=

−1

(15)

where 1

1

∫−1 N1(s)g(s)ds = ∫−1 12 s(−2 + em(s−1)(s + 1))g(s)ds 1 2 1 = 2 1 = 2 =

1

1

∫−1 −2sg(s)ds + 12 ∫−1 em(s−1)[s(s + 1)]g(s)ds ∫−1 −2sg(s)ds + 12 ∫−1 em(s−1)g (s)ds

(16)

1

1

∫−1 N2(s)g(s)ds = ∫−1 (1 + s − em(s−1)s(s + 1))g(s)ds 1

=

∫−1 ((1 + s)g(s) + em(s−1)[−s(s + 1)]g(s))ds

=

∫−1 (1 + s)g(s)ds + ∫−1 em(s−1)[−s(s + 1)]g(s)ds

=

∫−1 (1 + s)g(s)ds − ∫−1 em(s−1)g (s)ds

1

1

1

1

∫−1 em(s−1)[s(s + 1)]g(s)ds= 12 ∫−1 em(s−1)g (s)ds

1

1

1

1 2

(18)

n

weights for evaluation of the integral ∫ em(s −1)g (s )ds = ∑i =1 wg i (si ), −1 where g (s ) is a polynomial function of degree not greater than 5. Note that the m -values have been selected from Section 7 (Table 3) for different functions with different orders of singularity. As seen in Table 1, the larger the m-value, the closer to the singular side the Gauss points are. On the other hand, increasing the regular Gauss points causes the Gauss points to be distributed over the domain, which makes sure of having more Gauss points near the singular side. Thus, as we mentioned earlier, alternatively, one may evaluate the integrals with a sufficient number of classical (regular) Gauss points. In this paper, we have evaluated the integrals with the use of classical Gauss quadrature with 10 Gauss points to capture the rapid variations of the exponential function in EG-element shape functions. Table 2 repre-

∫−1 (−2sg(s) + em(s−1)[s(s + 1)]g(s))ds 1

⎞

In addition, g (s ) = s(s + 1)g(s ). Each of these integrals can be evaluated with the use of classic Gauss quadrature, with a sufficient number of Gauss points. On the other hand, one may evaluate these integrals with special Gauss points and Gauss weighting values when the function g (s ) is a polynomial. Table 1 presents the Gauss points and

1

Ni(s )g(s )ds

1⎛

1

∫−1 N3(s)g(s)ds = ∫−1 ⎜⎝ 12 em(s−1)s(s + 1)⎟⎠g(s)ds

1

1

sents the evaluation of the integral ∫ em(s −1)g (s )ds using (1) classic −1 Gauss quadrature with 10 Gauss points, (2) special Gauss quadrature with 3 Gauss points, and (3) analytical integration, where the accuracy

1

(17) 190

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Fig. 7. Eight-node EG-edge element shape functions for m = 10 .

Boundary integral equations for this problem can be formed by the following Green's function for a unit vertical ring load on the surface of a half-space [22]:

of the classic Gauss quadrature with 10 Gauss points may be seen. 6. Numerical examples

2(1 − ν ) r K (kl ), πμ r+r π /2 dθ K (kl ) = ∫ 0 (1 − klsin2θ )1/2

uzz(r; r ) =

6.1. Example one: contact of a rigid circular footing with an elastic half-space The mixed boundary value problem due to axisymmetric frictionless contact of a rigid circular disc of radius a undergoing a vertical displacement of δ with an isotropic half-space can be reduced to the determination of only the normal contact traction tz(r ) between the disc and the half-space (see Fig. 12). The exact solution for contact traction tz(r ) is given as [21]:

tz(r ) =

F (δ ) 2πa

1 2

a − r2

,

2μδ π (1 − ν )

1 a2 − r 2

4r r ≤ 1, (r + r )2 (21)

where uzz is the displacement component in the z− direction at the field point on the surface of the half-space when a unit ring load is applied in the z− direction at the source ring on the surface of the half-space; and K (kl ) denotes the complete elliptic integral of the first kind [22], r is the radial coordinate of the field point and r is the radius of the ring load. It should be noted that use of the half-space Green's function makes the integral on the left hand side of Eq. (14) identically equal to zero. Moreover, it can be shown by consideration of rigid-body movements in a semi-infinite body that cij = δij [19], where δij is the Kroneker delta. Thus, the BIE is reduced to:

0≤r≤a (19)

where a and r are the radius and radial coordinate of the disc, and F (δ ) is the total force corresponding to the vertical displacement δ . By integrating the stress distribution, one can find the constant settlement δ beneath the disc due to the total force F , which can be used to express the vertical traction in terms of the imposed vertical displacement as:

tz(r ) =

kl =

N

u=

⎛

L

∑ ⎜⎜∑ j =1

⎝ l =1

⎞ J l wl (u*NT(s, t ))l ⎟⎟ pn ⎠

(22)

where u* is the half-space Green's function. Owing to the use of the exact half-space Green's functions, the contact surface is the sole part of the boundary that should be discretized. To show the efficiency of the element introduced here, we first illustrate a comparison between the numerical results evaluated based on the quadratic elements (m = 0 in

,0≤r≤a (20)

where μ and ν are respectively the shear modulus and the Poisson's ratio of the half-space as the subgrade of the rigid disc. 191

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Fig. 8. Eight-node EG-corner element shape functions for m = n = 10 .

Fig. 9. Collapsing an eight-node square EG element into a six-node triangular EG element and mapping the triangular element into (s′, t′) coordinate system.

close to the exact solution. However, to recognize the best solution with the least error, we define an error function to measure the discrepancy between the analytical and numerical solutions over all elements as follows:

our EG-element), EG-elements with various values for m , and the exact solution for ν = 0.3 and E = 2μ(1 + ν ) = 30 GPa , using ten uniform elements along the radius of the disc. To this end, the one-dimensional three-node EG-element is used as the outmost element for the unknown traction, together with nine quadratic elements for the rest of the radius in a boundary element method for this example. As shown in Fig. 13, one can easily see that the solution for m = 0 (the quadratic element) fares the worst, possessing irregular fluctuations of the traction near the singularity of the edge. On the other hand, the solutions obtained with the use of EG-elements are all smooth and very

N −1

error (m ) =

0.9

1

∑e =1 ∫ ( f e (s ) − ge(x (s )))2 J e ds + ∫ ( f N (s ) − g N (x (s )))2 J N ds −1 −1 1

∫0 (g(x ))2dx (23)

× 100% e

e

where the functions f (x ) and g (x ) are the interpolated and analytical 192

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Fig. 10. Six-node EG-edge element shape functions for m = 10 .

Fig. 11. Six-node EG-corner element shape functions for m = n = 10 .

values of the surface traction over element e , and N is the number of elements. This error is usually called as relative L 2− norm error over the whole domain of the problem. The upper limit of integration is chosen thusly in order to avoid the singularity at the edge of the disc and increase accuracy. If the function error (m ) is plotted for various values of m , we get a minimum at a particular value of m . It can be seen from Fig. 14 that with a mesh of 10 uniform elements along the disc radius, the error function is minimum for a wide range of m between 11 and 16 and particularly for m = 13, and it is maximum for m = 0 (quadratic element). With m = 13 for the EG-element and a mesh of 10 non-uniform elements, we can see from Fig. 15 that the numerical

solution is very close to the exact solution not only in terms of value but also smoothness. 6.2. Example two: contact of a rigid strip footing with an elastic halfspace The next example deals with vertical displacement of an infinitely long rigid footing with a width of 2a resting on an elastic half-space (see Fig. 16). This example is utilized as a plane strain mixed boundary value problem in a Cartesian coordinate system. The half-space is defined by Ω = {(x , y ) − ∞ < x < + ∞, y > 0} and characterized by

193

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Table 1 Gauss points and weights for different values of EG-parameter.

Table 2 1 Comparing the numerical values of ∫ em(s −1)g (s )ds evaluated using classic Gauss −1 quadrature with 10 Gauss points, special Gauss quadrature with 3 Gauss points, and analytical integration.

m

Points, si

Weights, wi

6

−0.023592046 0.6236955348 0.9316585976

0.0019250580 0.0472345036 0.1175060810

8

0.2154116791 0.7136600647 0.9480999996

0.0013116378 0.0348765201 0.0888118281

0.3710866283 0.7705953281 0.9584263975

0.0010395784 0.0278551018 0.0711053195

0.5507183243 0.8361228674 0.9703018235

0.0007420907 0.0198941291 0.0507923516

0.7265241268 0.9002486800 0.9819228454

0.0004517068 0.0121094667 0.0309170874

0.8537222074 0.9466446427 0.9903308243

0.0002416106 0.0064771566 0.0165370467

10

14

23

43

F 2

π a −x

2

(2) g (s )

(3) Analytical integration

(4) Special Gauss quadrature with 3 Gauss points

(5) Classic Gauss quadrature with 10 Gauss points

(6) Error of column (5) with respect to Analytical solution (%)

6

1 s

s2

0.16667 0.13889 0.12037

0.16667 0.13889 0.12037

0.16667 0.13889 0.12037

0.00 0.00 0.00

s3

0.10648

0.10648

0.10648

0.00

s4

0.09568

0.09568

0.09568

0.00

s5 1 s

0.08694

0.08694

0.08694

0.00

s2

0.12500 0.10938 0.09766

0.12500 0.10938 0.09766

0.12500 0.10938 0.09766

0.00 0.00 0.00

s3

0.08838

0.08838

0.08838

0.00

s4

0.08081

0.08081

0.08081

0.00

s5 1 s

0.07449

0.07449

0.07449

0.00

s2

0.10000 0.09000 0.08200

0.10000 0.09000 0.08200

0.10000 0.09000 0.08200

0.00 0.00 0.00

s3

0.07540

0.07540

0.07540

0.00

s4

0.06984

0.06984

0.06984

0.00

s5 1 s

0.06508

0.06508

0.06508

0.00

s2

0.07143 0.06633 0.06195

0.07143 0.06633 0.06195

0.07143 0.06633 0.06195

0.00 0.00 0.00

s3

0.05815

0.05815

0.05815

0.01

s4

0.05481

0.05481

0.05481

0.01

s5 1 s

0.05185

0.05185

0.05184

0.02

s2

0.04348 0.04159 0.03986

0.04348 0.04159 0.03986

0.04343 0.04153 0.03978

0.11 0.15 0.19

s3

0.03828

0.03828

0.03818

0.26

s4

0.03682

0.03682

0.03670

0.33

s5 1 s

0.03547

0.03547

0.03532

0.42

s2

0.02326 0.02271 0.02220

0.02326 0.02271 0.02220

0.02216 0.02153 0.02093

4.71 5.20 5.72

s3

0.02171

0.02171

0.02035

6.27

s4

0.02124

0.02124

0.01978

6.84

s5

0.02079

0.02079

0.01924

7.44

8

10

the Young's modulus E = 2μ(1 + ν ) = 30 GPa and Poisson's ratio ν = 0.3. The width of the strip footing at y = 0 is discretized using one-dimensional quadratic elements and three-node EG-elements for the two edges. The exact contact stress distribution, which was discovered by M. Sadowsky in 1928, is given by [21]:

tz(x ) =

(1) m

14

,−a≤x≤a (24)

23

where F is the total load per unit length of the footing. From the line load solution, the displacement Green's function is [23]:

1−ν uzz*(r ) = ln(r ) πμ

(25) 43

where r = x − ξ is the distance between the field point x and the source point ξ on the surface of the half-space (see Fig. 17). We can conclude that the displacement under the rigid strip foundation is computed by: a

(1 − ν )F δ= tz(x )uzz*( x )dx = −a π 2μ

∫

a

∫−a

1 a2 − x 2

ln( x )dx (26)

which results in the vertical force, F , in terms of δ . Substituting for F in Eq. (24) gives the traction distribution in terms of the vertical rigid displacement. The analytical solutions given in Eqs. (24) and (25) are used as benchmarks in this example. In two-dimensional problems, when the collocation node lies within the integration element and the singularity of the kernel is of order ln(1/r ), weakly singular integrals appear [24]. For the integration of the integrals including ln(1/r ) as part of their kernel, one can use a modified Gauss Quadrature called Gauss-Laguerre integration [20] in the form of

∫0

1

⎛1⎞ f (s )ln⎜ ⎟ds ≈ ⎝s ⎠

L

∑ wl f (sl ) l =1

(27)

where L is the number of integration points. Note that for this integration scheme s = 0 is the singular point and the limits are s = 0 and s = 1, so a change in coordinates has to be made before Eq. (27) can be applied [20]. For the case where the source point, ξ , is located at the mid-side node, the integration has to split into two regions, one over −1 < s < 0 , and the other over 0 < s < 1. For the computation of product u*(ξ; s )NT(s ), as in (22), the intrinsic coordi-

Fig. 12. A rigid circular footing on the surface of a half-space.

194

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Fig. 13. Vertical traction under rigid frictionless disc for mesh of 10 uniform elements with EG element at the right end.

Fig. 16. A rigid strip footing on the surface of a half-space.

Fig. 14. The error function for different EG parameters. Fig. 17. Boundary value problem for two-dimensional surface load Green's functions.

discrepancy between the analytical and numerical solutions (see Fig. 18). However, as can be seen from the results, they do not strongly depend on the optimal value of the EG parameter. The results are accurate enough for a wide range of m between 11 and 16. Choosing m = 13 for the EG-element and a mesh of 10 non-uniform elements along 0 ≤ x ≤ a , Fig. 19 illustrates the vertical component of the contact traction for the surface punch problem against the exact solution and the case when only quadratic elements are used. The results again demonstrate the effectiveness of using EG boundary elements as a tool for improving the accuracy of numerical approximation.

6.3. Example three: contact of a rigid square footing with an elastic half-space Fig. 15. Vertical traction under rigid frictionless disc for mesh of 10 non-uniform elements.

To evaluate the performance of two-dimensional EG-elements, the problem of a square rigid punch in frictionless contact with an elastic half-space is considered (see Fig. 20). A surface square footing with a dimension of 2a × 2a is discretized using 8-node EG corner and EG edge quadrilateral elements for the outmost elements of the mesh, and the standard eight-node quadrilateral elements for all other traction, and geometric variations. To form the boundary integral equations of the problem, the vertical component of the displacement Green's function is required [19]:

nates for the two sub-regions are computed by s = − s and s = s for the two regions respectively [20]. With the standard quadratic element and the 3-node EG-element employed to interpolate the contact traction for the outmost element of the mesh, and standard quadratic elements for all other traction and geometric variations, numerical solutions are generated. The value of m is chosen in the same way as the preceding section, by minimizing the 195

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Fig. 21. Static stiffness of a square surface punch (frictionless contact).

Fig. 18. The error function for different EG parameters.

Fig. 22. Normal contact traction of a frictionless punch using 10×10 non-uniform mesh with m = 13 for outmost elements.

where uzz* is the displacement component in the z direction at the field point on the surface of the half space when a unit point load is applied in the z−direction at (ξ, η, ς ) = (0, 0, 0). In case of two-dimensional quadrilateral elements, one method for dealing with the singular integration of the kernel with (1/ r ) that occurs when the source point is a node of the integration element is to split the element up into triangular sub-elements. For each sub-element a local coordinate system should be chosen in such a way that the Jacobian of the transformation approaches zero at the source point. Numerical integration formulae are then applied over two or three sub-elements depending if the source point is a corner or mid-side node [20]. For comparison, the results obtained by the BEM approach for F vertical stiffness Kνν = δ , as done in this paper, are plotted beside the solution of Dempsey and Li [25] in Fig. 21. The relationship between the load per unit area, F , of a 2a × 2a square footing on an elastic halfspace and the rigid vertical displacement, δ , has been obtained by Dempsey and Li [25] by using an indirect integral equation approach as (see also [26]):

Fig. 19. Vertical traction under rigid frictionless strip footing for mesh of 10 nonuniform elements.

⎛ 4μa ⎞ F = 1.1523 ⎜ ⎟δ ⎝1 − ν ⎠

With the experience we have from the previous examples and since the numerical evaluations are not very sensitive with the m−value, we set m = 13 for this example. It is observed that when EG-elements are used, the vertical stiffness values converge much faster to the benchmark value compared to the case where only standard quadratic elements are used (Fig. 21). It should also be noted that even if we take another value such as m = 10 for this problem, there is still a significant improvement in the results (Fig. 21). The distribution of normal contact tractions is plotted in Fig. 22. The boundary element discretization used to obtain the results is non-uniform and consists of

Fig. 20. A rigid square footing on the surface of a half-space.

uzz*(x, y, z = 0; 0, 0, 0) =

1−ν 2πμ x 2 + y 2

(29)

(28)

196

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al.

Table 3 The best estimation for m -value for some singular functions with different degrees of singularity when 0 ≤ x ≤ 0.999 with the estimated error. Function

1/(1 − x 2 )1/5

1/(1 − x 2 )1/3

1/(1 − x 2 )1/2

1/(1 − x 2 )2/3

1/(1 − x 2 )

1/(1 − x 2 )3/2

m Error (%)

6 1.69

8 4.48

10 10.44

14 17.28

23 24.04

43 21.41

8. Conclusion Based on the concept of finite dimensional linear vector space, a family of one- and two-dimensional elements has been introduced in this paper, which can interpolate different functions from very regular to very high-gradient functions. The advantage of these elements is that the gradient of the shape functions can be adjusted with the use of only one parameter. The kernel of the shape functions is an exponentially gradient function, that's why we call these elements as ExponentiallyGradient elements (EG-elements). The elements can be used in both finite element and boundary element methods to simulate any regular or singular behavior. The shape functions are in such a way that the elements can be used for either point singularities as encountered in three dimensional boundary value problems or line singularities as can be seen in two-dimensional problems. Some examples have been solved to prove the capabilities of the EG-elements introduced in this paper, where excellent results are observed. Fig. 23. The error function in terms of EG-parameters for the function 1/(1 − x 2 )2/3.

Acknowledgement

10 × 10 quadrilateral elements, and the EG-parameter is 13. On top of the improvements observed in the convergence of the stiffness values, it should be mentioned that if EG-elements are used a significant improvement is resulted in the anomalies that always occur near the edges of the foundation when using quadratic elements for this problem. Fig. 22 shows a three dimensional variation of vertical contact tractions in between the rigid square footing and the half-space when EG-elements are used.

The partial support from the University of Tehran through 27840/ 1/07 to M. E.-G. during this work is gratefully acknowledged. References [1] Gladwell GML. Forced tangential and rotatory vibration of a rigid circular disc on a semi-infinite solid. Int J Eng Sci 1968;6:591–607. http://dx.doi.org/10.1016/ 0020-7225(68)90061-X. [2] Pak RYS. Asymmetric wave propagation in an elastic half-space by a method of potentials. J Appl Mech ASME 1987;54:121–6. [3] Eskandari-Ghadi M, Fallahi M, Ardeshir-Behrestaghi A. Forced vertical vibration of rigid circular disc on a transversely isotropic half-space. J Eng Mech 2010;136:913–22. http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000114. [4] Eskandari-Ghadi M, Mirzapour A, Ardeshir-Behrestaghi A. Rocking vibration of a rigid circular disc in a transversely isotropic full-space. Int J Numer Anal Methods Geomech 2011;35:1587–603. http://dx.doi.org/10.1002/nag.976. [5] Eskandari-Ghadi M, Ardeshir-Behrestaghi A. Forced vertical vibration of rigid circular disc buried in an arbitrary depth of a transversely isotropic half space. Soil Dyn Earthq Eng 2010;30:547–60. http://dx.doi.org/10.1016/j.soildyn.2010.01.011. [6] Sneddon IN. Mixed boundary value problems in potential theory. Amsterdam: North-Holland Publishing Company; 1966. [7] Erdogan F. Mixed boundary-value problems in mechanics. Mechanics today, 4. New York: Pergamon Press, Inc.; 1978. [8] Erdogan F. Mixed boundary value problems in mechanics of materials “Some Reflections on Forty Years of Solving Mixed Boundary Value Problems in Inhomogeneous Elasticity”. AIP Conference Proceedings, vol. 973, AIP; 2008, p. 772–83. doi: http://dx.doi.org/10.1063/1.2896880. [9] Erdogan F, Gupta GD. On the numerical solution of singular integral equations. Q Appl Math 1972;29:525–34. [10] Martin PA, Rizzo FJ, Cruse TA. Smoothness-relaxation strategies for singular and hypersingular integral equations. Int J Numer Methods Eng 1998;42:885–906. [11] Eskandari-Ghadi M, Mahmoodian M, Pak RYS, Ardeshir-Behrestaghi A. Analytical solution of torsion vibration of a finite cylindrical cavity in a transversely isotropic half-space. ZAMM Z Fur Angew Math Und Mech 2012;92:583–95. http:// dx.doi.org/10.1002/zamm.201100123. [12] Pak RYS, Abedzadeh F. Torsional traction on an open finite cylindrical cavity. Proc R Soc A Math Phys Eng Sci 1992;438:133–44. http://dx.doi.org/10.1098/ rspa.1992.0097. [13] Ardeshir-Behrestaghi A, Eskandari-Ghadi M, Navayi neya B, Vaseghi-Amiri J. Dynamic Reissner–Sagoci problem for a transversely isotropic half-space containing finite length cylindrical cavity. Soil Dyn Earthq Eng 2014;66:252–62. http:// dx.doi.org/10.1016/j.soildyn.2014.07.009. [14] Pak RYS, Ashlock JC. Method of adaptive-gradient elements for computational mechanics. J Eng Mech 2007;133:87–97. http://dx.doi.org/10.1061/(ASCE)07339399(2007)133:1(87). [15] Apostol TM. Calculas volume 1 one-variable Calculus, with an introduction to linear algebra, 2nd edition, 1. North-Holland Publishing Company: Wiley; 1967.

7. Estimating the best m -value for some singular functions In this section, we estimate the optimum EG-parameter, m -value, for some singular functions. By a curve-fitting procedure, we interpolate the singular functions 1/ q (1 − x 2 ) p using 10 uniform 3-node one-dimensional elements (9 quadratic elements plus 1 EG-element). To estimate the best value for m , we minimize the error function defined as

error (m ) =

2 ∼e 10 ∑e =1 ∫ (f (x ) − f e (x )) dx Γe

∫0

0.999

(f (x ))2 dx

× 100% ,

∼e f (x ) =

3

∑ Ni(x )f (xi ) i =1

(30)

∼e where the function f (x ) is the interpolated value of the function f (x ) = 1/ q (1 − x 2 ) p over each element. The best values for m -parameter derived in this manner are listed in Table 3 for 1/ q (1 − x 2 ) p . As can be seen, the value of m increases with p / q . Comparing the optimum m -value for 1/ 1 − x 2 with example one, we see that the value is 10 rather than 13. The reason is in example one, we had the interpolated solution in the right-most node at 1 from BEM, but here the right-most node is fixed at 0.999, as the function is singular at 1. However, it is clear from Fig. 14 that the error for 10 < m < 15 does not change very much. The optimum EG-parameter derived by minimizing the error function depends on the upper limit of the integrals in Eq. (30). However, any non-zero value for m improves the solution. As an example, the error function is plotted for 1/ 3 (1 − x 2 )2 in Fig. 23. 197

Engineering Analysis with Boundary Elements 80 (2017) 184–198

M. Eskandari-Ghadi et al. [16] Halmos PR. Finite-dimensional vector spaces. New York, NY: Springer New York; 1987. http://dx.doi.org/10.1007/978-1-4612-6387-6. [17] Liu GR. Meshfree Methods: Moving Beyond the Finite Element Method, Second Edition. 2009. doi: http://dx.doi.org/10.1063/1.289688010.1115/1.1553432. [18] Bathe KJ. Finite element procedures. New Jersey: Prentice Hall; 1996. [19] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Berlin Heidelberg: Springer; 1984. http://dx.doi.org/10.1007/978-3-642-48860-3. [20] Beer G, Smith I, Duenser C. The boundary element method with programming. Vienna: Springer Vienna; 2008. http://dx.doi.org/10.1007/978-3-211-71576-5. [21] Davis RO, Selvadurai APS. Elasticity and geomechanics. New York: Cambridge University Press; 1996. [22] Lubarda VA. Circular loads on the surface of a half-space: displacement and stress

[23]

[24] [25] [26]

198

discontinuities under the load. Int J Solids Struct 2013;50:1–14. http://dx.doi.org/ 10.1016/j.ijsolstr.2012.08.029. Telles JCF, Brebbia CA. Boundary element solution for half-plane problems. Int J Solids Struct 1981;17:1149–58. http://dx.doi.org/10.1016/0020-7683(81)900949. Aliabadi MH. The boundary element methodVolume 2: applications in solids and structures. Wiley; 2002. Dempsey JP, Li H. A rigid rectangular footing on an elastic layer. Géotechnique 1989;39:147–52. http://dx.doi.org/10.1680/geot.1989.39.1.147. Guzina BB, Pak RYS, Martinez-Castro AE. Singular boundary elements for threedimensional elasticity problems. Eng Anal Bound Elem 2006;30:623–39. http:// dx.doi.org/10.1016/j.enganabound.2006.02.010.