NURBS-enhanced boundary element method based on independent geometry and field approximation for 2D potential problems

Engineering Analysis with Boundary Elements 83 (2017) 158–166

Contents lists available at ScienceDirect

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

NURBS-enhanced boundary element method based on independent geometry and field approximation for 2D potential problems Wei Zhou a,b, Biao Liu a,b, Qiao Wang a,b,∗, Yonggang Cheng a,b, Gang Ma a,b, Xiaolin Chang a,b, Xudong Chen c a

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China School of Water Resources and Hydropower Engineering, Wuhan University, Wuhan 430072, China c College of Civil and Transportation Engineering, Hohai University, Nanjing 210098, China b

a r t i c l e

i n f o

Keywords: Non-uniform rational B-spline Isogeometric analysis Boundary element method Potential problems

a b s t r a c t Non-uniform rational B-spline (NURBS) in Isogeometric analysis (IGA) is coupled with the boundary element method (BEM) for 2D potential problems in this paper. The geometry and field are usually approximated by the same basis functions in IGA, such as the B-spline or the NURBS basis functions. In the proposed method, these two kinds of approximation are performed independently, i.e. the geometry is reproduced by the NURBS basis functions while the field is approximated by the traditional Lagrangian basis functions which are used in the conventional BEM. The proposed method has the advantage that the geometry can be reproduced exactly at all stages in IGA methods. Actually, one can use the computer aided design (CAD) software or NURBS library to perform the operations related to the geometry. The field approximation is performed in parameter space and separated from the geometry. Thus, it can be implemented easily as the conventional BEM since most algorithms for BEM can be applied directly, such as the methods for treatment of the singular integrals, and the boundary conditions can be imposed directly. Numerical examples have demonstrated the accuracy of the proposed method. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Isogeometric analysis (IGA) [1,2] has achieved many attentions in recent years and has been implemented in the finite element method (FEM) successfully. In IGA, the analysis can be performed directly from the model obtained by computer aided design (CAD) software, i.e. the geometry data used to construct the model can be used as the input file for numerical analysis. The basic idea is using the same parametric functions in CAD for geometry modeling to approximate the fields in numerical analysis, and usually the non-uniform rational B-spline (NURBS) is applied. Based on this idea, IGA has at least two advantages than conventional methods: (i) the geometry of the computed domain can be reproduced exactly at all the stages; (ii) the refinement becomes much simpler since it can be done in parameter space based on the known geometry data, and the refinement for NURBS can be performed easily by knot insertion or order elevation. In most CAD software, the boundary representation data is used to model the solid and the boundary surfaces or boundary curves are represented by parametric functions, such as NURBS basis functions. Thus the numerical methods based on the boundary integral equation (BIE), such as the boundary element method (BEM), seem to be ideally suited for

∗

IGA. BEM [3–6] has been applied in many areas since only the boundary of the domain is needed to be discretized into elements and high accuracy can be obtained. Coupling the BEM with IGA is a natural idea and it has already attracted many researchers’ attentions. Simpson et al. [7] proposed an isogeometric boundary element method for 2D elastostatic analysis. An extended isogeometric BEM (XIBEM) which uses partition of unity enriched NURBS functions was proposed by Peake et al. [8] for 2D Helmholtz problems. Gu et al. [9] introduced the IGA in BIE for 3D potential problem and the local bivariate B-spline function was deducted in the work. The application of isogeometric boundary element method for 3D tunnel simulation was investigated by Marussig et al. [10]. A nonsingular isogeometric boundary element method for Stokes flows in 3D case was presented by Heltai et al. [11] and the convergence was tested numerically. Nguyen et al. [12] applied the IGA in the weakly singular symmetric Galerkin boundary element method (SGBEM) for two-dimensional crack problems. To further reduce the computational cost of the isogeometric boundary element method, Takahashi and Matsumoto [13] applied the fast multipole method (FMM) [14–18] to accelerate the method for Laplace equation in 2D case. In most of the isogeometric boundary element methods proposed, the geometry and the fields are approximated by the same basis

Corresponding author at: State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China. E-mail address: [email protected] (Q. Wang).

http://dx.doi.org/10.1016/j.enganabound.2017.07.013 Received 11 November 2016; Received in revised form 7 June 2017; Accepted 17 July 2017 0955-7997/© 2017 Elsevier Ltd. All rights reserved.

W. Zhou et al.

Engineering Analysis with Boundary Elements 83 (2017) 158–166

6

functions, such as the B-spline basis functions or NURBS basis functions [19,20], which is known as the isoparametric fashion and much higher accuracy can be obtained [21]. However, if using the NURBS basis functions to approximate the fields, the collocation points need to be redefined since the control points may not lie on the boundary and the basis functions may not satisfy the property of Kronecker-delta. In addition, the computation of singular integrals may also become more difficult and the rigid-body motion method can no longer be used to compute the strongly singular integrals [7]. In this paper, a NURBS-enhanced boundary element method (NEBEM) for 2D potential problems is proposed. In NEBEM, the geometry and fields are approximated independently, namely, the geometry is approximated by NURBS basis functions while the fields are approximated by Lagrangian basis functions in parameter space. By doing this, the geometry of the model can also be reproduced exactly at all stages like the isogeometric boundary element method, while the fields can be approximated easily like the conventional BEM. Some other NURBS-enhanced methods have been proposed, such as the NURBSenhanced finite element method [22], the NURBS-enhanced reproducing kernel (RK) approximation [23], the NURBS-enhanced maximumentropy schemes [24] and the IGA-LME meshfree approximants [25,26]. Compared with BEM, there are at least three advantages for the proposed method: (i) the discretization error of geometry is avoid, which results in higher accuracy. (ii) the mesh procedure is much simpler since the knot vector is used to create the elements, and the refinement of NURBS is very simple. Some researchers even classified some IGA methods into meshfree methods based on this point. (iii) the proposed method can be integrated with Computer Aided Design (CAD) like many other IGA methods. Compared with IGA BEM, the proposed method have two advantages: (i) the filed approximation is much simpler and the shape functions have the delta function property, which results in easy code programming and easy implementation of boundary condition. (ii) the treatment of singular integrals is simpler since many methods in BEM can be applied in the propose method directly, such as the rigid body motion method can still be applied to compute the strongly singular integrals. This paper is organized as follows. The B-spline and NURBS are introduced briefly in the second section and the NEBEM for 2D potential problems is proposed in the third section. Finally, numerical examples are given in the fourth section to demonstrate the accuracy of the proposed method.

5 4 3

B-spline

x2

2

0 -1 -2 -3 -5

-2

-1

0

1

2

3

4

5

2.1. B-spline basis functions and B-spline curves The B-spline curves and surfaces are based on the B-spline basis functions. Let Ξ = {t0 ,⋅⋅⋅, tm′ } be a nondecreasing sequence of real numbers, i.e., ti ≤ ti + 1 , i is from 0 to m′ − 1. Then the B-spline basis functions can be defined as { 1, 𝑡 ∈ [𝑡𝑖 , 𝑡𝑖+1 ) 𝑁𝑖0 (𝑡) = (3) 0, elsewise 𝑁𝑖𝑝 (𝑡) =

𝑡𝑖+𝑝+1 − 𝑡 𝑡 − 𝑡𝑖 𝑁 𝑝−1 (𝑡) + 𝑁 𝑝−1 (𝑡), 𝑝 ≥ 1 𝑡𝑖+𝑝 − 𝑡𝑖 𝑖 𝑡𝑖+𝑝+1 − 𝑡𝑖+1 𝑖+1

(4)

where p is the degree of the basis functions, which denoted by 𝑁𝑖𝑝 (𝑡). ti are called knots and Ξ is the knot vector. If 𝑡𝑗−1 < 𝑡𝑗 = 𝑡𝑗+1 = ⋯ 𝑡𝑗+𝑚𝑗 −1 < 𝑡𝑗+𝑚𝑗 , then knot tj is repeated mj times. The knot vector is called open knot vector if knots are repeated p + 1 times at the start and end of the vector. Then a pth-degree B-spline curve can be defined as 𝑛 ∑ 𝑖=0

𝑁𝑖𝑝 (𝑡)𝐏𝑖

(5)

where Pi are the control points and C(t) is a piecewise polynomial curve. The degree p, the number of control points n + 1, and the number of knots m′ + 1 are related by

as

𝑚′ = 𝑛 + 𝑝 + 1

(1)

(6)

Fig. 1 shows an example of geometry based on two B-spline curves with their associate control points. The control points for the left curve are P0 = {0, 5}, P1 = { − 2, 5}, P2 = { − 4, 4}, P3 = { − 2, −2} and P4 = {0, 0} while the control points of the right curve are P0 = {0, 0}, P1 = {2, −2}, P2 = {4, 4}, P3 = {2, 5} and P4 = {0, 5}. The degree of the curves are 2 and the knot vector for both curves are Ξ = {0, 0, 0, 1/3, 2/3, 1, 1, 1}. The B-spline basis functions for the knot vector are also illustrated in Fig. 2.

where t is an independent parameter, x{x1 ,x2 } are the coordinates of a point on the curve, which can be represented separately as an explicit function of t by Eq. (1). C(t) is a vector-valued function denoting the coordinates of points on the curve by parameter t. For example, the parametric functions of a circle centered at point {0, 0} can be written as 0 ≤ 𝑡 ≤ 2𝜋

-3

Fig. 1. B-spline curves and control points.

In 2D space, the parametric form of a general curve can be written

𝑥2 = 𝑟 sin(𝑡),

-4

x1

2. B-spline and non-uniform rational B-spline (NURBS)

𝑥1 = 𝑟 cos(𝑡),

Control point

1

𝐂(𝑡) =

𝐂(𝑡) = (𝑥1 (𝑡), 𝑥2 (𝑡))

B-spline

(2)

where r is the radius of the circle. In engineering software, it may be hard to construct analytical functions like those in Eq. (2) to represent a very complex curve or surface. To overcome this problem, most CAD packages using piecewise polynomial functions to construct the curves and the non-uniform rational B-spline (NURBS) basis functions are one of the most investigated methods. In this section, the B-spline and NURBS is introduced briefly and more details can be found in book [20]. The B-spline is outlined first since the NURBS is based on the B-spline approximations.

2.2. NURBS basis functions and NURBS curves B-spline can be considered as a special case of NURBS and the basis functions of NURBS can be defined by projecting the B-spline basis functions using a set of weights as 𝑅𝑝𝑖 (𝑡) =

𝑤𝑖 𝑁𝑖𝑝 (𝑡)

𝑛 ∑

𝑖′ =0

159

𝑤𝑖′ 𝑁𝑖𝑝′ (𝑡)

(7)

W. Zhou et al.

Engineering Analysis with Boundary Elements 83 (2017) 158–166

2

N0

1.0

2

R0

1.0

2

N1

2

R1

2

N2

0.8

2

R2

0.8

2

N3

2

R3

2

N4

2

R4

0.6

Ri (t)

2

2

Ni (t)

0.6

0.4

0.4

0.2

0.2

0.0

0.0 0.0

0.2

0.4

0.6

0.8

0.0

1.0

0.2

0.4

0.6

0.8

1.0

t

t Fig. 2. B-spline basis functions.

Fig. 4. NURBS basis functions.

6

-1

of isogeometric boundary element method is using the same parametric functions used by CAD to approximate the geometry and field. By doing this, the geometry can be reproduced exactly and the analysis can be performed from the data in CAD directly. In this section, a NURBS-enhanced boundary element method (NEBEM) is proposed for 2D potential problems. Unlike the IGA BEM, the parametric functions are only used to reproduce the geometry of the curves in NEBEM and the conventional Lagrangian basis functions are applied to approximate the field in parameter space. The control equation for potential problems is Laplace equation and the well-known regularized boundary integral equation can be written as

-2

𝟎=

5 4 3

B-spline B-spline

2

Control point

x2

NURBS

1

NURBS

0

-3 -5

-4

-3

-2

-1

0

1

2

3

4

Fig. 3. NURBS curves and control points.

𝑛 ∑ 𝑖=0

𝑅𝑝𝑖 (𝑡)𝐏𝑖

∫Γ

𝑞 ∗ (𝐱, 𝐲)[𝑢(𝐲) − 𝑢(𝐱)]dΓ(𝐲)

(9)

𝑢∗ (𝐱 , 𝐲 ) =

1 1 ) ln( 2𝜋 𝑟(𝐱 , 𝐲 )

(10)

𝑞 ∗ (𝐱 , 𝐲 ) =

𝜕 𝑢∗ (𝐱 , 𝐲 ) 𝜕𝐧(𝐲)

(11)

𝑁𝑖𝑝 (𝑡)

where are the NURBS basis functions, wi are the weights and are the B-spline basis functions defined in Eqs. (3) and (4). Then a pthdegree NURBS curve can be defined as 𝐂(𝑡) =

𝑢∗ (𝐱, 𝐲)𝑞(𝐲)dΓ(𝐲)−

where Γ is the boundary of the bounded domain Ω, x and y are points on the boundary, u and q are the potential and normal flux on the boundary, u∗ (x, y) and q∗ (x, y) are fundamental solutions of Laplace equation and in 2D case can be written as

5

x1

𝑅𝑝𝑖 (𝑡)

∫Γ

where r(x, y) = ‖x − y‖ is the distance between x and y, and n is the unit outward normal at the boundary point. The discretization of Eq. (9) includes two parts, one part is partitioning the boundary into cells for the convenience of integration and the other part is approximating the field by discretizing the problem into finite degrees of freedom (DOFs).

(8)

From Eq. (7) it can be indicated that the NURBS basis functions will reduce to B-spline basis functions if all the weights are set to 1.0 since 𝑛 ∑ 𝑁𝑖𝑝′ (𝑡) = 1, which is known as the partition of unity property of B𝑖′ =0

spline and NURBS basis functions. NURBS can reproduce many geometries exactly while B-spline only approximate such shapes, such as circular arcs and spheres. Using the same control points and knot vector as for B-spline curve shown in Fig. 1, Fig. 3 shows the NURBS curves with weights w = {1, 1, 2, 2, 1} for the left curve and w = {1, 2, 2, 1, 1} for the right curve. The relative NURBS basis functions for the left curve are plotted in Fig. 4.

3.1. Geometry approximation In NEBEM, the geometry is approximated by parametric functions and in this paper they are NURBS introduced in Section 2. The integral cells, which can be called isogeometric boundary elements, are defined in parameter space and the unique knot vector can be used directly as the boundary elements. The unique knot vector is a subset of knot vector Ξ and all the knots are repeated just once. For the NURBS curves in Fig. 3, the unique knot vector Λ = {0, 1/3, 2/3, 1} and the isogeometric elements can be seen in Fig. 5. Suppose the unique knot vector is defined as Λ = {𝜉 0 ,⋅⋅⋅, 𝜉 m } with 𝜉 j < 𝜉 j + 1 , j is from 0 to m − 1. Then the NURBS based boundary can be

3. NURBS-enhanced boundary element method In conventional BEM, the geometry and field are approximated by Lagrangian basis functions. The discretization error for geometry approximation cannot be avoided for complex geometries. The key idea 160

W. Zhou et al.

Engineering Analysis with Boundary Elements 83 (2017) 158–166

6

3.2. Field approximation

5

Rather than the field is approximated by NURBS basis functions in most IGA, in NEBEM it is approximated by Lagrangian basis functions as used in conventional BEM. The approximation is performed on each isogeometric element in parameter space and the collocation nodes can be constructed on each element very easily as conventional BEM. If each element contains two nodes at the ends, i.e. using the unique knot vector Λ as the set of collocation nodes, one can define a new local coordinate transformation at element Γj [𝜉 j ,𝜉 j + 1 ] as

4 3

Control point NURBS Element edge

x2

2 1 0

𝑡(𝑠) =

-1

𝜉𝑗+1 − 𝜉𝑗 2

𝑠+

𝜉𝑗+1 + 𝜉𝑗

(16)

2

where s ∈ [ − 1, 1] is the local coordinate and 𝑠0𝑗 = −1, 𝑠1𝑗 = 1 are nodes on element Γj , and 𝐱𝑗+𝑖 = 𝐱𝑗𝑖 = 𝐱(𝑠𝑖𝑗 ), where i = 0, 1 (see Fig. 7). Then the field u(y) and q(y) can be approximated by

-2 -3 -5

-4

-3

-2

-1

0

1

2

3

4

5

𝑛

x1

𝑢(𝐲) = 𝑢(𝐲(𝑡(𝑠))) =

𝑗 ∑

𝑖=0

Fig. 5. Elements for NURBS curves.

𝜓𝑗𝑖 (𝑠)𝑢𝑖𝑗

(17)

𝜓𝑗𝑖 (𝑠)𝑞𝑗𝑖

(18)

𝑛

𝑞(𝐲) = 𝑞(𝐲(𝑡(𝑠))) =

𝑗 ∑

𝑖=0

where 𝑢𝑖𝑗 and 𝑞𝑗𝑖 are the displacements and tractions at point 𝑠𝑖𝑗 . nj + 1 is the number of nodes in each element and in this case nj = 1. 𝜓𝑗𝑖 (𝑠) are the shape functions defined on Γj as

Fig. 6. Isogeometric boundary elements.

𝑚 −1 ∑ 𝑗=0

Γ𝑗 [𝜉𝑗 , 𝜉𝑗+1 ]

𝑛 ∑ 𝑖=0

𝑅𝑝𝑖 (𝑡)𝐏𝑖

𝑗=0

=

∫Γ𝑗

𝑡(𝑠) = 𝜉𝑗1 −

𝑗=0

∫ Γ𝑗

𝑞 ∗ (𝐱, 𝐲(𝑡))[𝑢(𝐲(𝑡)) − 𝑢(𝐱)]𝐽 d𝑡

𝜉𝑗2 − 𝜉𝑗0 2

𝑠+

𝜉𝑗2 + 𝜉𝑗0 − 2𝜉𝑗1 2

𝑠2

(21)

and the shape functions are defined as

(13)

𝜓𝑗0 (𝑠) = 0.5𝑠(𝑠 − 1)

(22)

𝜓𝑗1 (𝑠) = (1 − 𝑠)(1 + 𝑠)

(23)

𝜓𝑗2 (𝑠) = 0.5𝑠(𝑠 + 1)

(24)

In Eq. (21), 𝜉𝑗0 , 𝜉𝑗1 and 𝜉𝑗2 can be three knots, i.e. one can also use the knot vector to create quadratic elements. However, in this paper, we only use two knots to define a quadratic element by defining 𝜉𝑗0 =

𝜉𝑗 , 𝜉𝑗2 = 𝜉𝑗+1 and 𝜉𝑗1 = (𝜉𝑗 + 𝜉𝑗+1 )∕2, then Eq. (21) becomes the same as Eq. (16), which results in an easier computation. Higher order elements can be obtained in the same ways. Substituting Eqs. (17) and (18) into Eq. (14), one can have

𝑢∗ (𝐱, 𝐲(𝑡))𝑞(𝐲(𝑡))𝐽 d𝑡

𝑚 −1 ∑

(20)

𝑠2𝑗 = 1. The corresponding local coordinate transformation for the quadratic element Γj is

(12)

where y are vectors of the geometric coordinates. Eq. (13) is a global approximation and it seems to be very time consuming. However, the NURBS basis functions have the property of local support: 𝑅𝑝𝑖 (𝑡) = 0 for t ∉ [ti ,ti + p + 1 ]. In addition, at most p + 1 of the basis functions are nonzero in any given knot span. Thus, the approximation in Eq. (13) can be computed efficiently with local basis functions and efficient algorithm can be found in Ref. [20]. By partitioning the boundary into isogeometric elements, Eq. (9) can be written as 𝑚 −1 ∑

𝜓𝑗1 (𝑠) = 0.5 + 0.5𝑠

and are zero for s ∉ Γj . The above element is linear element, and one can define a quadratic element with three nodes as shown in Fig. 8, where 𝑠0𝑗 = −1, 𝑠1𝑗 = 0 and

where Γj [𝜉 j ,𝜉 j + 1 ] is the ith element with range [𝜉 j ,𝜉 j + 1 ] in parameter space (see Fig. 6). The approximation for the geometry can be stated as follows by Eq. (8): 𝐲(𝑡) =

(19)

𝜓𝑗𝑖 (𝑠)

subdivided into isogeometric boundary elements as Γ=

𝜓𝑗0 (𝑠) = 0.5 − 0.5𝑠

𝑛

𝑚 −1 ∑

(14)

𝑗=0

where t are parameters of points y and J is the Jacobi coefficient computed by √ ( ) ( ) 𝜕 𝑦1 2 𝜕 𝑦2 2 𝐽 = + (15) 𝜕𝑡 𝜕𝑡

=

∫Γ𝑗

𝑢∗ (𝐱, 𝐲(𝑡(𝑠)))

𝑖=0

𝜓𝑗𝑖 (𝑠)𝑞𝑗𝑖 𝐽 𝐽1 ds 𝑛

𝑚 −1 ∑ 𝑗=0

𝑗 ∑

∫Γ𝑗

𝑞 ∗ (𝐱, 𝐲(𝑡(𝑠)))[

𝑗 ∑

𝑖=0

𝑛

𝜓𝑗𝑖 (𝑠)𝑢𝑖𝑗 −

𝑗 ∑

𝑖=0

𝜓𝑗𝑖 (𝑣)𝑢𝑖𝑗 ]𝐽 𝐽1 ds

(25)

where v is the local coordinate of point x in the relative element and J1 = (𝜉 j + 1 − 𝜉 j )/2. 161

W. Zhou et al.

Engineering Analysis with Boundary Elements 83 (2017) 158–166

Fig. 7. Local coordinate for linear elements.

Fig. 8. Local coordinate for quadratic element.

Suppose xl , where l = 0, 1, ⋅⋅⋅m − 1, are collocation nodes on the boundary, and l is the global index of a node. Then for each node xl , Eq. (25) can be written as

∫Γ𝑗

𝑢∗ (𝐱𝑙 , 𝐲(𝑡(𝑠)))

𝑚 −1 ∑

𝑖=0

𝑞 ∗ (𝐱𝑙 , 𝐲(𝑡(𝑠)))

∫Γ𝑗

𝑗=0

𝑗 ∑

4

𝜓𝑗𝑖 (𝑠)𝑞𝑗𝑖 𝐽 𝐽1 ds

[ 𝑛𝑗 ∑ 𝑖=0

3 𝑛

𝜓𝑗𝑖 (𝑠)𝑢𝑖𝑗 −

𝑗 ∑

𝑖=0

] 𝜓𝑗𝑖 (𝑣𝑙 )𝑢𝑖𝑗 𝐽 𝐽1 ds

𝑗=0

=

∫Γ𝑗

𝑢 (𝐱𝑙 , 𝐲(𝑡(𝑠)))

𝑗 ∑

-2

𝑖=0

𝜓𝑗𝑖 (𝑠)𝑞𝑗𝑖 𝐽 𝐽1 ds

-3 -5

𝑞 ∗ (𝐱𝑙 , 𝐲(𝑡(𝑠)))

∫Γ𝑗

𝑚 −1 ∑

𝑗 ∑

𝑖=0

-3

𝐇𝐮 = 𝐆𝐪

∫ Γ𝑗

𝑞 ∗ (𝐱𝑖 , 𝐲(𝑡(𝑠)))𝜓𝑘 (𝑠)𝐽 𝐽1 d𝑠

∑

−𝛿𝑖𝑘

Γ𝑗 ∈K𝑖

∑ Γ𝑗 ∈K𝑘

∫Γ𝑗

∫Γ𝑗

𝑞 ∗ (𝐱𝑖 , 𝐲(𝑡(𝑠)))𝐽 𝐽1 d𝑠

𝑢∗ (𝐱𝑖 , 𝐲(𝑡(𝑠)))𝜓𝑘 (𝑠))𝐽 𝐽1 d𝑠

[ 𝐮 = 𝑢0 ,

𝑢1 ,

...,

𝑢𝑚

[ 𝐪 = 𝑞0 ,

𝑞1 ,

...,

𝑞𝑚

1

2

3

4

5

node xi . From Eqs. (28)–(32), it can be indicated that the formulations for NEBEM are very similar to those in conventional BEM, and the methods for treatment of the singular integrals in BEM can also be applied in NEBEM. In addition, the boundary conditions can also be imposed directly as in BEM.

where

Γ𝑗 ∈K𝑘

0

(27)

(28)

∑

-1

Fig. 9. Elements and control points after h-refinement.

Finally, Eq. (27) can be written in matrix form as

𝐻𝑖𝑘 =

-2

x1

𝜓𝑗𝑖 (𝑠)𝑢𝑖𝑗 𝐽 𝐽1 ds

𝑞 ∗ (𝐱𝑙 , 𝐲(𝑡(𝑠))𝑢𝑘𝑗 𝐽 𝐽1 ds

∫Γ𝑗

𝑗=0

𝐺𝑖𝑘 =

-4

𝑛

𝑚 −1 ∑ 𝑗=0

−

-1

𝑛

∗

1 0

where vl is the local coordinate of point xl computed from element Γj and 𝜓𝑗𝑖 (𝑣𝑙 ) ≡ 0 if vl ∉ Γj . Since the shape function 𝜓𝑗𝑖 (𝑣𝑙 ) in Eq. (26) has the interpolating property, Eq. (26) can be written as 𝑚 −1 ∑

Control point NURBS Element edge

2

(26)

x2

=

5

𝑛

𝑚 −1 ∑ 𝑗=0

6

]T

3.3. Refinement methods (29)

The geometry approximation can capture the curves exactly, but the field approximation may still lead to errors. To reduce the errors from field approximation, refinement of the model is still needed and it will not change the shape of the geometry at all stages. There are three popular methods for refinement:

(30)

(32)

(i) h-refinement, which is known as knot insertion; (ii) p-refinement, which is known as order elevation; (iii) k-refinement, which is a combination of h-refinement and prefinement.

where xi ,i = 0, 1, ⋅⋅⋅m are collocation nodes, 𝜓 k (s),k = 0, 1, ⋅⋅⋅m are the global shape functions for node xk . One can note that 𝜓 k (s) are nonzero only if s ∈ Γj . In Eqs. (29) and (30), Kk are the set of elements which have contribution to 𝜓 k (s) and Ki are the set of elements containing

In this paper the h-refinement is applied for all the examples. Fig. 9 shows the elements and control points for the NURBS curves in Fig. 3 after h-refinement.

]T

(31)

162

W. Zhou et al.

Engineering Analysis with Boundary Elements 83 (2017) 158–166

3.4. Treatment of singular integrals The accuracy of the method is highly dependent on the accuracy of singular integrals. Many methods have been proposed for singular integrals in BEM [27–31] and most of these methods can be applied in the proposed method directly. In this paper, the strongly singular integrals are computed indirectly by the rigid body motion method, while the self-adaptive coordinate transformation method [32] is applied to evaluate the weakly singular integrals. Suppose the one-dimensional integral is 𝐼=

1

∫−1

𝑓 (𝜂)d𝜂

(33)

where f(𝜂) is singular at a point 𝜂. ̄ A transformation is chosen as 𝜂(𝛾) = 𝑎𝛾 3 + 𝑏𝛾 2 + 𝑐𝛾 + 𝑑

(34)

such that the following requirements are met: d2 𝜂 |𝜂̄ = 0 d𝛾 2

(35)

d𝜂 | =0 d𝛾 𝜂̄

(36)

𝜂(1) = 1

(37)

𝜂(−1) = −1

(38)

Fig. 10. Relative error of linear solution.

The unknown parameters a, b, c and d can be obtained by solving Eqs. (35), (36), (37) and (38). And Eq. (33) can be computed by 𝐼=

1

∫−1

𝑓 (𝜂(𝛾))𝐽2 d𝛾

where 𝐽2 =

(39)

𝜕𝜂 . 𝜕𝛾

4. Numerical examples Three numerical examples have been studied in this section. For the purpose of error estimation, a formula is defined as √( )/( 𝑁 ) √ 𝑁 √ ∑ ∑ (𝑒) 2 (𝑒) (𝑛) 2 √ 𝑒= (𝑢𝑖 − 𝑢𝑖 ) (𝑢𝑖 ) (40) 𝑖=1

where tively.

𝑢(𝑖𝑒)

and

𝑖=1

𝑢(𝑖𝑛)

Fig. 11. Relative error of quadratic solution.

refer to the exact and numerical solutions, respec-

4.1. Dirichlet problems Dirichlet problems for the model with NURBS curves shown in Fig. 3 are investigated and three types of analytical solutions are employed as: (i) Linear solution: 𝑢 = 𝑥1 + 𝑥2

(41)

(ii) Quadratic solution 𝑢 = 𝑥21 − 𝑥22 + 𝑥1 𝑥2

(42)

(iii) Cubic solution 𝑢 = 𝑥31 + 𝑥32 − 3𝑥21 𝑥2 − 3𝑥1 𝑥22

(43)

Seven refinements are applied to discretize the model and the relative errors of u, q1 = 𝜕 u/𝜕 x1 and q2 = 𝜕 u/𝜕 x2 for each solution are plotted in Figs. 10, 11 and 12, respectively. It can be indicated that the proposed method has good convergence rates and the accuracy is very high.

Fig. 12. Relative error of cubic solution.

163

W. Zhou et al.

Engineering Analysis with Boundary Elements 83 (2017) 158–166 Table 1 Knot vectors and control points for model in example 3. Curve ID

Knot vector

Control points

Weights

1 (p = 2) 2 (p = 2) 3 (p = 2) 4 (p = 2) 5 (p = 2)

0.0, 0.0, 0.0, 1.0, 1.0, 1.0 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 0.0, 0.0, 0.0, 1.0, 1.0, 1.0

{0.0, 2.0}; {0.0, 1.0}; {0.0, 0.0} {0.0, 0.0}; {0.0, 1.5}; {3.0, 0.0} {3.0, 0.0}; {3.0, 1.0}; {4.0, 1.0} {4.0, 1.0}; {4.0, 1.5}; {4.0, 2.0} {4.0, 2.0}; {2.0, 2.0}; {0.0, 2.0}

1.0, 1.0, 1.0, 1.0, 1.0, 1.0 1.0, 0.707106781, 1.0 1.0, 1.0, 1.0 1.0, 1.0, 1.0

Fig. 13. Model of flow around a cylinder.

4.2. Potential flow A potential flow around a cylinder of radius a = 1 in an infinite domain is considered in this example. Due to the symmetry of the problem, only a part of the upper left quadrant of the field is modeled as shown in Fig. 13, where L = 4 and W = 2 in the figure. The exact solution for this problem is given by [33] 𝑢 = 𝑥2 [1 −

𝑥22

𝑎2 ] + (𝑥1 − 𝐿)2

Fig. 14. Results of potential for example 2.

(44)

where u represents the stream function. The model is constructed by five separated NURBS curves shown in Fig. 13 and the data of the curves can be found in Table 1, where p is the degree of the NURBS curve. The boundary conditions of this problem are as follows: Curve 1: 𝑢 = 𝑥2 (𝑥22 + 15)∕(𝑥22 + 16); Curve 2: u = 0; Curve 3: u = 0; Curve 4: 𝜕 u/𝜕 n = 0; Curve 5: u = [6 + 2(x1 − 4)2 ]/[4 + (x1 − 4)2 ]. The results of the potential on the diagonal line, i.e. from point {0, 0} to point {4, 2} are plotted in Fig. 14. The results obtained from the conventional BEM with linear elements are also plotted for the purpose of comparison. Both the results are obtained by 15 elements and one can observe that the accuracy of the NEBEM is a little higher than the results obtained by BEM. To further show the high accuracy of NEBEM, more elements are constructed and the relative errors of NEBEM and BEM are shown in Fig. 15. From Fig. 15, it can be observed that high rates of convergence are achieved by both NEBEM and BEM, and the accuracy has been improved by NEBEM. However, the gap from the results of the proposed method and BEM tends to be unchanged and difference between the proposed method and BEM decreases when the number of elements increases. One of the reasons is that the discretization error of geometry becomes smaller and smaller when the number of elements increases for traditional BEM.

Fig. 15. Relative error for example 2.

(i) Solution 1: 𝑢 = 𝑥41 +𝑥42 −6𝑥21 𝑥22

(45)

(ii) Solution 2: 𝑢 = sinh(𝑥1 ∕10) sin(𝑥2 ∕10)

4.3. Mixed problems

(46)

Both linear and quadratic elements proposed in Section 3 are applied to approximate the fields. The relative errors of u, q1 and q2 computed from 90 inner points with respective to the number of nodes (degrees of freedom) are shown in Figs. 17 and 18. It can be indicated that the error becomes smaller as the number of nodes increases for both linear and quadratic elements. In addition, the relative errors of u, q1 and q2 for both solutions obtained by quadratic elements are always less than those

A model shown in Fig. 16, which constructed in a CAD software with mixed boundary conditions is studied in this example. Neumann boundary conditions are imposed on lines 1, 2 and 3 in Fig. 16, and Dirichlet boundary conditions are imposed on the rest boundaries. Two kinds of analytical solutions are used for this example as 164

W. Zhou et al.

Engineering Analysis with Boundary Elements 83 (2017) 158–166

Fig. 18. Relative error of solution 2 for example 3.

Fig. 16. Model constructed in CAD software.

Acknowledgments Financial support for the project from the National Natural Science Foundation of China (No. 51609181) is acknowledged. References [1] Hughes TJ, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput. Methods Appl Mech Eng 2005;194(39):4135–95. [2] Cottrell JA, Hughes TJ, Bazilevs Y. Isogeometric analysis: toward integration of CAD and FEA. John Wiley & Sons; 2009. [3] Jaswon MA, Symm GT. Integral equation methods in potential theory and elastostatics. United Kingdom: Oxford Univ Press; 1977. [4] Yao Z-H, Xu J-D, Wang H-T, Zheng X-P. Simulation of CNT composites using fast multipole BEM. J Mar Sci Technol 2009;17(3):194–202. [5] Yao Z, Wang H. Some benchmark problems and basic ideas on the accuracy of boundary element analysis. Eng Anal Boundary Elem 2013;37(12):1674–92. [6] Tan F, Jiao Y-Y. The combination of the boundary element method and the numerical manifold method for potential problems. Eng Anal Boundary Elem 2017;74:19–23. [7] Simpson RN, Bordas SP, Trevelyan J, Rabczuk T. A two-dimensional isogeometric boundary element method for elastostatic analysis. Comput Methods Appl Mech Eng 2012;209:87–100. [8] Peake MJ, Trevelyan J, Coates G. Extended isogeometric boundary element method (XIBEM) for two-dimensional Helmholtz problems. Comput Methods Appl Mech Eng 2013;259:93–102. [9] Gu J, Zhang J, Li G. Isogeometric analysis in BIE for 3-D potential problem. Eng Anal Boundary Elem 2012;36(5):858–65. [10] Marussig B, Beer G, Duenser C. Isogeometric boundary element method for the simulation in tunneling. Appl Mech Mater 2014;553:495–500. [11] Heltai L, Arroyo M, DeSimone A. Nonsingular isogeometric boundary element method for Stokes flows in 3D. Comput Methods Appl Mech Eng 2014;268:514–39. [12] Nguyen B, Tran H, Anitescu C, Zhuang X, Rabczuk T. An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems. Comput Methods Appl Mech Eng 2016;306:252–75. [13] Takahashi T, Matsumoto T. An application of fast multipole method to isogeometric boundary element method for Laplace equation in two dimensions. Eng Anal Boundary Elem 2012;36(12):1766–75. [14] Wang H, Yao Z. A new fast multipole boundary element method for large scale analysis of mechanical properties in 3D particle-reinforced composites. Comput Model Eng Sci 2005;7(1):85–95. [15] Liu Y. A new fast multipole boundary element method for solving large‐scale two‐dimensional elastostatic problems. Int J Numer Methods Eng 2006;65(6):863–81. [16] Liu Y, Nishimura N. The fast multipole boundary element method for potential problems: a tutorial. Eng Anal Boundary Elem 2006;30(5):371–81. [17] Wang Q, Zhou W, Cheng Y, Ma G, Chang X. A line integration method for the treatment of 3D domain integrals and accelerated by the fast multipole method in the BEM. Comput Mech 2017;59(4):611–24. [18] Wang Q, Zhou W, Cheng Y, Ma G, Chang X, Huang Q. The boundary element method with a fast multipole accelerated integration technique for 3D elastostatic problems with arbitrary body forces. J Sci Comput 2017;71(3):1238–64. [19] Rogers DF. An introduction to NURBS: with historical perspective. Elsevier; 2000. [20] Piegl L, Tiller W. The NURBS book. Springer Science & Business Media; 2012. [21] Li K, Qian X. Isogeometric analysis and shape optimization via boundary integral. Comput-Aided Des 2011;43(11):1427–37.

Fig. 17. Relative error of solution 1 for example 3.

obtained by linear elements when the number of nodes increases, and the accuracy of quadratic elements is much higher than the accuracy of linear elements.

5. Conclusions In this paper, an NURBS-enhanced boundary element method (NEBEM) is proposed for 2D potential problems. The geometry is reproduced by the NURBS basis functions, which are commonly used in CAD software, while the field is approximated by traditional Lagrangian basis functions like in BEM. Thus, the method can be considered as a combination of isogeometric boundary element and conventional BEM. It has the advantage that the geometry can be reproduced exactly at all the stages as in IGA. In addition, the methods for treatment of the singular integrals in BEM can be applied directly in NEBEM and the boundary conditions can also be imposed directly as in BEM. The method will be time-consuming for large scale computation like BEM. Coupling it with the FMM [34–36] is under research. The domain integrals will also exist for nonhomogeneous problems in the proposed method as in BEM, and using the dual reciprocity method (DRM) [37,38] or other methods [39–42] to overcome this problem is possible. Multi-domain problems with interface [43] will also be investigated in the future work. 165

W. Zhou et al.

Engineering Analysis with Boundary Elements 83 (2017) 158–166

[22] Sevilla R, Fernández‐Méndez S, Huerta A. NURBS‐enhanced finite element method (NEFEM). Int J Numer Methods Eng 2008;76(1):56–83. [23] Chi S-W, Lin S-P. Meshfree analysis with the aid of NURBS boundary. Comput Mech 2016;58(3):371–89. [24] Greco F, Coox L, Maurin F, Desmet W. NURBS-enhanced maximum-entropy schemes. Comput Methods Appl Mech Eng 2017;317:580–97. [25] Rosolen A, Arroyo M. Blending isogeometric analysis and local maximum entropy meshfree approximants. Comput Methods Appl Mech Eng 2013;264:95–107. [26] Greco F, Coox L, Desmet W. Maximum-entropy methods for time-harmonic acoustics. Comput Methods Appl Mech Eng 2016;306:1–18. [27] Lv J-H, Feng X-T, Chen B-R, Jiang Q, Guo H-S. The CPCT based CBIE and HBIE for potential problems in three dimensions. Eng Anal Boundary Elem 2016;67:53–62. [28] Lv J-H, Feng X-T, Li S-J, Jiang Q, Guo H-S. Solid analysis of micron-sized thin structures with BEM for steady-state heat conduction. Eng Anal Boundary Elem 2016;71:11–19. [29] Gao X-W, Zhang J-B, Zheng B-J, Zhang C. Element-subdivision method for evaluation of singular integrals over narrow strip boundary elements of super thin and slender structures. Eng Anal Boundary Elem 2016;66:145–54. [30] Gao X-W, Feng W-Z, Yang K, Cui M. Projection plane method for evaluation of arbitrary high order singular boundary integrals. Eng Anal Boundary Elem 2015;50:265–74. [31] Rong J, Wen L, Xiao J. Efficiency improvement of the polar coordinate transformation for evaluating BEM singular integrals on curved elements. Eng Anal Boundary Elem 2014;38:83–93. [32] Telles J. A self‐adaptive co‐ordinate transformation for efficient numerical evaluation of general boundary element integrals. Int J Numer Methods Eng 1987;24(5):959–73. [33] Zhang J, Yao Z, Li H. A hybrid boundary node method. Int J Numer Methods Eng 2002;53(4):751–63.

[34] Wang Q, Miao Y, Zheng J. The hybrid boundary node method accelerated by fast multipole expansion technique for 3D elasticity. Comput Model Eng Sci 2010;70(2):123–51. [35] Wang Q, Miao Y, Zhu H, Zhang C. An O (N) fast multipole hybrid boundary node method for 3D elasticity. Comput Mater Continua 2012;28(1):1–25. [36] Wang Q, Miao Y, Zhu H. A fast multipole hybrid boundary node method for composite materials. Comput Mech 2013;51(6):885–97. [37] Zhou F, Zhang J, Sheng X, Li G. Shape variable radial basis function and its application in dual reciprocity boundary face method. Eng Anal Boundary Elem 2011;35(2):244–52. [38] Zhou F, Zhang J, Sheng X, Li G. A dual reciprocity boundary face method for 3D non-homogeneous elasticity problems. Eng Anal Boundary Elem 2012;36(9):1301–10. [39] Zhou W, Wang Q, Cheng Y, Ma G. A fast multipole method accelerated adaptive background cell-based domain integration method for evaluation of domain integrals in 3D boundary element method. Eng Anal Boundary Elem 2016;67:1–12. [40] Gao X-W. The radial integration method for evaluation of domain integrals with boundary-only discretization. Eng Anal Boundary Elem 2002;26(10):905–16. [41] Wang Q, Zhou W, Cheng Y, Ma G, Chang X. Line integration method for treatment of domain integrals in 3D boundary element method for potential and elasticity problems. Eng Anal Boundary Elem 2017;75:1–11. [42] Wang Q, Zhou W, Cheng Y, Ma G, Chang X, Huang Q. An adaptive cell-based domain integration method for treatment of domain integrals in 3D boundary element method for potential and elasticity problems. Acta Mech Solida Sin 2017;30(1):99–111. [43] Gao XW, Feng WZ, Zheng BJ, Yang K. An interface integral equation method for solving general multi‐medium mechanics problems. Int J Numer Methods Eng 2016;107:696–720.

166