An isogeometric BEM using PB-spline for 3-D linear elasticity problem

Engineering Analysis with Boundary Elements 56 (2015) 154–161 Contents lists available at ScienceDirect Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound An isogeometric BEM using PB-spline for 3-D linear elasticity problem Jinliang Gu a,n, Jianming Zhang b, Lifeng Chen a, Zhihua Cai a a School of Electromechanical Engineering, Hunan University of Science and Technology, Xiangtan 411201, China State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, College of Mechanical and Vehicle Engineering, Hunan University, Changsha 410082, China b art ic l e i nf o a b s t r a c t Article history: Received 19 June 2014 Received in revised form 23 February 2015 Accepted 23 February 2015 Available online 16 March 2015 PB-spline is introduced in isogeometric BEM schema for solution of 3-D linear elasticity problems. In this schema, bivariate PB-spline is employed not only to construct the exact geometric model but also to approximate the boundary variables. PB-spline as an unstructured spline technology can largely reduce the computation cost for analysis of some special geometric models, such as a sphere, where lots of nearly singular and singular integrals will appear. Numerical tests for a sphere and a tube, show that the new method has good performance. & 2015 Elsevier Ltd. All rights reserved. Keywords: Isogeometric BEM IGABEM NURBS PB-spline BEM Linear elasticity 1. Introduction The conception of isogeometric analysis (IGA) combined with Finite Element Method(i.e., IGAFEM) was first proposed by Hughes et al. [1] and further developed by Cottrell et al. [2,3], Bazilevs et al. [4,5], and Zhang et al. [6]. Its goals are to generalize and improve on FEM in following aspects: (1) To provide exact geometry for analysis no matter how coarse the discretization. (2) To simplify mesh refinement by eliminating the need for communication with the CAD geometry once the initial mesh is constructed. (3) To more tightly weave the mesh generation process within CAD. Within the concept, basis functions generated from non-uniform rational B-spline (NURBS) play a key role in offering exact geometry representation, simplification of design optimization and tighter integration of modeling-analysis. NURBS is a tool in the most widespread use for implementation of IGA. It may be owing to the two facts, one is that NURBS is the standard approach for representation of free form curves and sculptured surfaces in CAD, and can represent elementary shapes such as sphere, cylinders, and torus exactly. This makes it possible for NURBS to express all of CAD models in a uniform geometric representation. Another is that the basis functions originated in the field of approximation theory inherently have some useful mathematic and geometric properties, such as the local support, nonnegative and partition of unit. These are attractive properties with regard to n Corresponding author. E-mail address: gujinliang1983@163.com (J. Gu). http://dx.doi.org/10.1016/j.enganabound.2015.02.013 0955-7997/& 2015 Elsevier Ltd. All rights reserved. numerical stability. Areas of application of NURBS-based or, more recently, T-splines-based IGAFEM are varied and continue to grow rapidly. However, one challenging but far from trivial task needed for resolution in IGAFEM is how to generate a trivariate NURBS or T-spline representation preserving a given NURBS or T-spline surface representation [7]. It precludes the application of the IGAFEM to complex modules and assemblages. The IGA coupled with boundary element method (BEM), i.e. IGABEM, has achieved rapid development recently [8,9], and has been successfully applied to deal with engineering problems [10– 12]. Compared with traditional BEM, 3-D bodies for analysis in IGABEM are exact geometrical representation no matter how coarse the discretization of these bodies is, it ensures that no geometrical errors produce in the analysis process. Representation of these bodies only needs to be encompassed by their boundary surfaces based on Boundary representation (B-rep) [13,14], due to the reduction of problem dimensionality by one in the BEM. This is a distinct advantage over IGAFEM, in which construction of 3-D bodies requires construction of their domains, not merely their boundary surfaces. The challenging problem of IGAFEM is avoided naturally here. It is clear that the surface modeling descriptions provided by CAD lend themselves naturally to BEM and this link represents a significant advance for both the mathematical and engineering communities. It is worthy of mention that the attempt for using exact geometric information in BEM has been performed early. The tool for implementation of it is called as Boundary Face Method (BFM) that first proposed by J.M. Zhang et al. [15]. In the scheme, exact 3-D geometric models for analysis are based on B-rep, which are constructed by CAD software, such as UG-NX. J. Gu et al. / Engineering Analysis with Boundary Elements 56 (2015) 154–161 Solution of BIEs is performed in the parametric space of geometric model that is similar with the IGA. But tools for constructing CAD model and for solving the BIEs are so different that mesh refinement needs iterative communication with CAD geometry. Thus, the BFM is not a seamless method as IGA. Nearly singular and singular integrals are key factors to deteriorate the exactness of numerical results in BEM. However, parametric mesh for NURBS surfaces must be plotted as a rectangular grid because of tensor product of a global definition of NURBS in 2-D parametric plane. This means that each row lines or column lines of this mesh traverses in their corresponding parameter directions. This mesh fashion can result in mesh singularities on some body surfaces, such as singular points on the two poles of a sphere. These singularities leads to many nearly singular and singular integrals require for evaluation in BEM. PB-spline is a point-based spline developed from NURBS, it allows a more flexible mesh fashion than NURBS in 2D or 3D parametric space. Here, we use a special case of bivariate PB-spline for performance of isogeometric BEM to solve 3-D linear elasticity problem. It combines part of the flexibility of PB-spline with the topology and structure of NURBS. It also capable of being locally refined in the same manner as T-spline keeping the original geometry and parameterization unchanged [16,17]. Its parametric mesh can allow a variable number of knots subdivision in one parametric direction. The mesh may contain a lot of T-junctions. Thus, it also can be considered as a special case of T-spline. The performance of isogeometric BEM on the special PB-spline mesh can largely alleviate the computation cost for nearly singular and singular integrals in some cases. This paper is organized as follows. In Section 2, there are expressions of B-spline and NURBS, and PB-spline functions. In Section 3, the implementation of IGA in BIEs of 3-D linear elasticity is performed. Numerical examples for 3-D linear elasticity problems are given in Section 4. Finally, we present the conclusions for our work. 2. B-spline, NURBS and PB-spline PðξÞ ¼ n X j¼1   Bi;k ξ w i pi n   P Bi;k ξ wi 155 ð4Þ j¼1 wi is the weight value, if all weights wi are set to 1, the NURBS curve (Eq. (4)) degenerates into B-spline curve (Eq. (3)). Given a control net fpij g and two series of basis functions:  fBi;k ðξÞg and fBj;l ðηÞg built from two knot vectors Ξ 1 ¼ ξ1 ; ξ2 ; ⋯; ξm þ k þ 1 Š and Ξ 2 ¼ η1 ; η2 ; ⋯; ηn þ l þ 1 , (i ¼ 1; 2; ⋯m, j ¼ 1; 2; ⋯n), bivariate B-spline and NURBS can be expressed as Pðξ; ηÞ ¼ m X n X i¼1j¼1 Pðξ; ηÞ ¼ m X n X i¼1j¼1     Bi;k ξ Bj;l η pij     Bi;k ξ Bj;l η wij pij m P n     P Bi;k ξ Bj;l η wij ð5Þ ð6Þ i¼1j¼1 If all weights wij are set to 1, the NURBS surface (Eq. (6)) degenerates into the B-spline surface (Eq. (5)). Note that this process performs on a global perspective of knot vectors as the construction of B-spline or NURBS needs all of the control points and basis functions defined in the entire parametric domain. Obviously, mesh for bivariate B-spline or NURBS must be kept as an uniform rectangular grid in 2D parametric domain because all tensor products Bi;k Bj;l are defined in global knot vector (we call it as global tensor product), such as mesh for global tensor product described in Fig. 2. 2.2. Blending function and PB-spline   If we are only interested in a basis functionBi;k ξ contained in   interval ξi ; ξi þ k þ 1 , we have  no need for the global knot vector. Only the local knot vector ξi ; ξi þ 1 ; ⋯; ξi þ k þ 1 contributes to the   definition of Bi;k ξ [16,17]. For example, an intact two degree basis function defined in a local knot vector ½ξi ; ξi þ 1 ; ξi þ 2 ; ξi þ 3 Š can be 2.1. Basis functions, B-spline and NURBS B-spline, NURBS are built from B-spline basis functions. These basis functions are defined recursively for zero degree, ( 1; ξi r ξ r ξi þ 1 Bi;k ðξÞ ¼ k¼0 ð1Þ 0; otherwise Fig. 1. Basis functions defined in an open knot vector [0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1]. and for non-zero degrees,   ξ Bi;k ξ ¼ ξi B ξi þ k ξi i;k 1   ξ þ ξi þ k þ 1 ξ B ξi þ k þ 1 ξi þ 1 i þ 1;k 1   ξ; k40 ð2Þ   Assuming that 0=0 : ¼ 0 and Ξ ¼ ξ1 ; ξ2 ; ⋯; ξn þ k þ 1 , where Ξ is a knot vector; ξi is the ith knot value; k and n are respectively degree and total number of these basis functions. In this paper, we use open knot vector with these knot values given by: ξ1 ¼ ξ2 ¼⋯ ¼ ξk þ 1 ¼ 0,ξn þ 1 ¼ ξn þ 2 ¼ ⋯ ¼ ξn þ k þ 1 ¼ 1, and ξi ¼ i ðk þ 1Þ =ðn kÞfor i ¼ k þ 2; k þ 3; ⋯n. Using Eq. (2), it is easy to deduct high degree basis functions defined successively in a global interval [0, 1], such as three degree basis functions illustrated in Fig. 1. Given a set of control points fpi g corresponding to these basis functions, univariate B-spline and NURBS can be constructed as a global fashion. PðξÞ ¼ n X j¼1   Bi;k ξ pi ð3Þ Fig. 2. A parametric mesh for a bivariate B-spline or a NURBS. 156 J. Gu et al. / Engineering Analysis with Boundary Elements 56 (2015) 154–161 expressed as 8 2 ð ξ ξi Þ > > > ξ ξ ξ > ð i þ 2 i Þð i þ 1 ξi Þ; ξi r ξ o ξi þ 1 > > < ðξ ξ Þðξ ξÞ ðξi þ 3 ξÞðξ ξi þ 1 Þ i iþ2 Bi;2 ¼ ðξi þ 2 ξi Þðξi þ 2 ξi þ 1 Þ þ ðξi þ 3 ξi þ 1 Þðξi þ 2 ξi þ 1 Þ; > > > 2 > > ð ξi þ 3 ξÞ > ; ξi þ 2 r ξ o ξi þ 3 : ξ ð i þ 3 ξi þ 1 Þ ð ξi þ 3 ξi þ 2 Þ ξi þ 1 r ξ o ξi þ 2 ð7Þ   In this case, the Bi;k ξ can be called blending function as no space is defined. In order to create a curve from several blending functions, we should make an assumption that each blending function is associated with a local knot vector and a control point, and these local knot vectors are kept some specified topology with each other. For a one dimension (1-D) case in Fig. 3, four segments extracted respectively from four different blending functionsBi;3 , ð4Þ ð2Þ ð1Þ , Bð3Þ Bi þ 1;3 , Bi þ 2;3 , Bi þ 3;3 are denoted as Bi;3 i þ 1;3 , Bi þ 2;3 , Bi þ 3;3 , defined in the subinterval [0.4, 0.6] ( In fact, theBi;3 , Bi þ 1;3 , Bi þ 2;3 , Bi þ 3;3 are respectively equal to the B3;3 , B4;3 , B5;3 , B6;3 of the Fig. 1.), spline curve defined in this subinterval can be constructed by SðξÞ ¼ iþk X Bkj;kþ 1 þ i j¼i j  ξ pj ð8Þ ðξÞ is the lth segment of the Bj;k , pj is the given control where BðlÞ j;k point corresponding to the Bj;k , k is the degree of the Bj;k . The construction of this curve has no need for the global knot vector Ξ ¼ ½0; 0; 0; 0; 0:2; 0:4; 0:6; 0:8; 1; 1; 1; 1Š, and only requires the set of local knot vectors: Ξ i ¼ ½0; 0; 0:2; 0:4; 0:6Š, Ξ i þ 1 ¼ ½0; 0:2; 0:4; 0:6; 0:8Š, Ξ i þ 2 ¼ ½0:2; 0:4; 0:6; 0:8; 1Š, Ξ i þ 3 ¼ ½0:4; 0:6; 0:8; 1; 1Š, to define the blending functions associated with four control points corresponding to them. This curve can be considered as a special case of a univariate PB-spline as the blending functions satisfy automatically a partition of unity. Generally, expression (8) is further defined as SðξÞ ¼ iþk X j¼i Bkj;kþ 1 þ i j iþk  .X kþ1þi ξ pj Bj;k j¼i j ξ  ð9Þ to ensure a partition of unity. In the bivariate PB-spline case, a “control cloud” is utilized for description of the set of control points as no organization or topology is used for their configuration. Several of the nice properties of NURBS persist in PB-splines, but several undesirable features also emerge to engineering application. Firstly, control points completely at random will result in degenerate geometries. Secondly, as each blending function from its local knot vector is constructed without regard for any other function, there is no clearly identifiable element in the classical sense, and it is difficult to speak of traditional refinement. Finally, there is no clear way of assessing the approximability of discretization spaces comprised of PB-splines. Here, we use a special case of bivariate PB-spline for performance of IGA. Its parametric mesh can allow a variable number of knots subdivision in one parametric direction. Fig. 4 shows a example of this mesh, if we wish to insert one knot into arbitrary row-line of the initial mesh, it only creates guilty of “Violation 1” Fig. 4. A parametric mesh for the special case of bivariate PB-spline. presented in the T-spline local refinement algorithm (for more details see [18]i.e., several blending functions miss the inserted knot in their local knot vectors for the current mesh). So it can be refined as h-refinement of univariate B-spline or NURBS. If we can build a relationship between local knot vectors and control points for this topology of the parametric mesh as the Fig. 4 shows, the special bivariate PB-spline can be constructed. In fact, the T-spline has a series of strategy to deal with it, such as defining the anchors in the parametric mesh to infer local knot vectors and blending functions, and building the one-to-one relationship between each anchor and each control point. Thus, this can be defined by r þk X sþl   X kþ1þr S ξ; η ¼ Bi;k i¼rj¼s i ξ Blj;lþ 1 þ s  j  η pij ð10Þ     where Bki;kþ 1 þ i r ξ and Blj;lþ 1 þ j l η are blending functions built respectively from local knot vectors h  ½ξi ;hξi þ 1 ; ⋯; ξi þ k Š and  ½ηj ; ηj þ 1; ⋯; ηj þ l Š, and ξ A ξr þ k 1 ξr þ k , η A ηs þ k 1 ηs þ k . For full generality, we assume that each control point has an associated weight wij , the expression can be further written as r þk X sþl   X   S ξ; η ¼ N ij ξ; η pij ð9Þ i¼rj¼s   where N ij ξ; η is defined by:   N ij ξ; η ¼ Bki;kþ 1 þ r ξ Bj;ll þ 1 þ s i  j  η wij rþk X sþl .X i¼rj¼s kþ1þr Bi;k ξ Bj;ll þ 1 þ s i  j ð10Þ If the 2-D parametric domain is meshed by an uniform rectangular grid as the Fig. 2, the bivariate PB-spline defined in this domain degenerates into B-spline or NURBS. 3. Implementation of IGA in BIE The implementation of IGA in our scheme is using the PBspline to construct exact geometric model in CAD and approximate physical variables in CAE in a same parametric space. The exact geometric information for analysis is easy to be obtained from the PB-spline funcition. Here, we use it for the solution of 3-D linear elasticity BIEs. 3.1. BIE for 3-D linear elasticity problem Fig. 3. Four blending functionsBi;3 , Bi þ 1;3 , Bi þ 2;3 �Bi þ 3;3 located in the subinterval [0.4, 0.6]. In order to perform the IGABEM by our method, we directly introduce the 3-D regularized displacement BIEs of linear elasticity, which is deducted from the PDEs of linear elasticity (for more details see [19,20]). The BIEs is written as Z Z   0 ¼ t j ðsÞU ij ðs; yÞdΓ ðsÞ uj ðsÞ uj ðyÞ T ij ðs; yÞdΓ ðsÞ ð11Þ Γ Γ  η wij J. Gu et al. / Engineering Analysis with Boundary Elements 56 (2015) 154–161 ¼ X i 157   Ni ξ; η p^ i ð16Þ   ^ where  N i ξ; η and pi are respectively the vector elements of ½N i Š and p^ i , which correspond one-to-one with the determinant   elements   of ½Nrs Š and ½prs Š, respectively. In Fig. 5, Nr0 ξ; η ¼ Nr1 ξ; η ¼ ⋯ ¼ N rl ξ; η . If the number of elements in every rows     are different from    each other, we have N r0 ξ; η a N r1 ξ; η a ⋯ a N rl ξ; η . N rs ξ; η does not satisfy the Kronecker delta property, but it is constricted to interpolate the node values p~ ij , i.e.,  X  p~ j ¼ S ξj ; ηj ¼ N i ξj ; ηj p^ i ð17Þ i Fig. 5. The parametric mesh of an integral boundary B-spline surface. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article). where uj ðsÞ and t j ðsÞ are respectively the displacement and stress function on the boundary Γ , which can be approximated by bivariate PB-spline. Where s is the field point (integration point), and y is the source point (collocation point). The kernels U ij ðs; yÞ and T ij ðs; yÞ are the Kelvin's fundamental solutions, which are given as 1  U ij ðs; yÞ ¼ 16π G 1 T ij ðs; yÞ ¼   3 μr 1  8π 1  μ r2  4μ δij þr ;i r ;j  ∂r  1 ∂nðsÞ i ¼ ð12Þ  1  2μ r ;i nj ðsÞ r ;j ni ðsÞ   ð13Þ j XX i ¼ j X j    2μ δij þ 3r ;i r ;j where p~ j is an arbitrary vector element corresponding to the node values p~ ij . Substituting Eq. (17) into Eq. (16) with inverse transformation, we have    X  X ½N Šij 1 p~ j S ξ; η ¼ N i ξ; η   N i ξ; η ½N Šij 1~ pj   N~ j ξ; η p~ j ð18Þ   P   where N~ j ξ; η ¼ N i ξ; η ½NŠij 1 i To use the Eq. (18), it is easy to construct the PB-spline approximation functions of boundary variables: uðsÞ and qðsÞ, respectively, such as   X     N~ i ξ; η u~ ξi ; ηi u ξ; η ¼ ð19Þ i 3.2. PB-spline approximation In our IGABEM schema, bivariate PB-spine defined in a parametric domain is used to represent integral boundary Γ of a body and to approximate the boundary values: uj ðsÞ and t j ðsÞ. Therefore, discretization for the Eq. (11) can be performed in the parametric domain. Fig. 5 shows a parametric mesh which maps to a 3-D B-rep model. The red dots denote the parametric coordinates of nodes. The purple dashed mesh-grids are integral background elements which correspond one-to-one to the nodes, for example the labeled domain of a integral element Γ uv . The black meshgrids is parametric subdivision according to knot vectors   of a quadratic bivariate B-spline: Ξ ξ ¼ 0; 0; 0; 1=3; 2=3; 1; 1; 1 and  Ξ η ¼ 0; 0; 0; 1=2 ; 1; 1; 1Š. The virtual values as the control points, are located on the corner points of each black grids, denoted by p^ ij . The node values are denoted by p~ ij . Given boundary condition is imposed on these red dots. As the Fig. 5 mesh fashion, the Eq. (11) can be discreted as XZ XZ   0¼ t j ðsÞU ij ðs; yÞdΓ ðsÞ uj ðsÞ uj ðyÞ T ij ðs; yÞdΓ ðsÞ mn Γ mn mn Γ mn ð14Þ Subsequently, the PB-spline approximation function can also be constructed with the p^ ij as control points, which is expressed as   S ξ; η ¼ k X l X r ¼0s¼0   N rs ξ; η prs ð15Þ whereprs ¼ p^ ðir rÞðjs sÞ , ir A k þ 1; ⋯; n andjs A l þ 1; ⋯; m . In Fig. 4, n ¼ 5, m ¼ 4 and k ¼ l ¼ 2. The Eq. (15) can be expanded, which then is expressed as another form, such as k X l X     S ξ; η ¼ N rs ξ; η prs r ¼0s¼0 ¼ k X r¼0 k k X X       N r0 ξ; η pr0 þ Nr1 ξ; η pr1 þ ⋯ þ N rl ξ; η prl r¼1 r¼l   X     N~ i ξ; η q~ ξi ; ηi q ξ; η ¼ ð20Þ i     where u ξ; η and q ξ; η are the 2-D parametric expressions of uðsÞ and qðsÞ, respectively. Substituting Eq. (19) and Eq. (20) in the Eq. (14), we have X ð1Þ X ð2Þ I mn þ I mn ¼ 0 ð21Þ mn mn ð1Þ and I ð2Þ where I mn mn are expressed as ! Z Z  1 1 X   ~ j ðξ; ηÞ N~ j ðξ ; η Þ T ij ðξ; ηÞJðξ; ηÞdξdη u~ j N I ð1Þ ¼ 0 0 mn l l 1 1 ð22Þ I ð2Þ mn ¼ X Z l 1 1 Z 1 1 !   ~ N rs ðξ; ηÞU ij ðs; yÞJðξ; ηÞdξdη t~ j l ð23Þ         where u~ j l ¼ u~ j ξl ; ηl and t~ j l ¼ t~ j ξl ; ηl . Then, the Eq. (14) can be written as a matrix form H u~ Gq~ ¼ 0 ð24Þ 4. Numerical examples Two 3-D geometric models are applied in the solution of linear elasticity BIEs to illustrate the accuracy and convergence of the IGABEM. These models are based on B-rep, i.e., these bodies are only encompassed by boundary surfaces. It is a distinct advantage over IGAFEM where 3-D models must be constructed by trivariate B-spline functions. To evaluate the error of numerical results, a ‘global’ L2 norm error, normalized by jvjmax , is defined by [21] vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N 2 X ðeÞ 1 u t1 e¼ v vði nÞ ð25Þ jvjmax N i ¼ 1 i 158 J. Gu et al. / Engineering Analysis with Boundary Elements 56 (2015) 154–161 Fig. 6. A sphere shell and its main size: (a) sphere shell; (b) planform of sphere shell and (c) side elevation of sphere shell. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article). Fig. 7. Two different mesh models for the sphere shell: (a) mesh according to B-spline definition; (b) mesh according to PB-spline definition. Table 1 The L2 errors of stress components from two fashions of meshes. Analytical field Degree u ¼ Linear Quadratic Cubic u ¼ Quadratic Quadratic Cubic Err_ Err_ Err_ Err_ Err_ Err_ Err_ Err_ Err_ Err_ Err_ Err_ Sxx (%) Syy (%) Szz (%) Sxx (%) Syy (%) Szz (%) Sxx (%) Syy (%) Szz (%) Sxx (%) Syy (%) Szz (%) where jvjmax is the maximum value of displacement u or stress t over N sample points. The superscripts ðeÞ and ðnÞ refer to the exact and numerical solutions, respectively. To assess the accuracy of numerical results in the IGABEM, we use the following three analytical fields, which are taken from reference [21]. (i) Linear solution: ux ¼ x þ 0:5y þ0:5z 254 nodes (PB-spline) 248 nodes (B-spline) 0.09043 0.08236 0.1192 0.04038 0.07122 0.07106 0.2936 0.3445 0.3951 0.3206 0.2437 0.1924 0.1605 0.1643 0.1943 0.07222 0.07838 0.06552 0.6451 0.6182 0.704 0.6878 0.688 0.8276 uz ¼ 3x2 þ 3y2 2z2 ð27Þ In all cases, the 3-D regularized displacement BIEs for linear elasticity is solved, combined with reasonable prescribed boundary conditions corresponding to the above analytical solutions. 4.1. Linear elasticity problem for a sphere shell uy ¼ 0:5x þ yþ 0:5z uz ¼ 0:5x þ 0:5y þ z (ii) Quadratic solution: ux ¼ 2x2 þ 3y2 þ 3z2 uy ¼ 3x2 2y2 þ3z2 ð26Þ A sphere shell is firstly used for discussion, which is constructed by three quadratic NURBS surfaces: two semisphere surfaces and one disc surface. The shell with its main size is shown in Fig. 6. The Fig. 6(b) and (c) are respectively the planform and side elevation of the body. Red dashed line in the Fig. 6(b) and (c) is used to denote the location of reference numerical solutions J. Gu et al. / Engineering Analysis with Boundary Elements 56 (2015) 154–161 4 2 Stress Sxx 0 -2 -4 Quadratic exact solution 94 nodes 254 nodes 414 nodes -6 -8 0.0 0.2 0.4 0.6 0.8 1.0 Parametric coordinates Fig. 8. Variation of the stress Sxx, along the red circle in Fig. 6(b), solved by cubic PB-spline. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article). -2 -3 -4 159 further results in many nearly singular and singular integrals need to be estimated in the BIEs. This is a challenging problem in BEM, which deteriorates the final numerical results. However, the PB-spline approach shown in Fig. 7(b) allows a variable number of nodes in each meridian, which allows for fewer nodes near the poles. Table 1 shows the numerical errors from the two fashion of mesh in Fig. 7. L2 error of the nodal stress values: Sxx, Syy and Szz, are denoted respectively by Err_Sxx, Err_Syy and Err_Szz, which are evaluated by Eq. (25). From the table, we can clearly see that the nodal values evaluated from 254 nodes by using the PB-spline are more accurate than those from 248 nodes by using the B-spline. In fact, much time for solving the BIEs also has been saved using PBspline compared with B-spline. This case indicates that PB-spline is more flexible for performance of IGABEM than B-spline. In Fig. 8 and Fig. 9, cubic PB-spline as the approximation function is used for evaluation of the convergence to the exact solutions. The stress components sampled from the red dashed line in Fig. 6(b) and (c) are considered. Three groups of nodes: 94 nodes; 254 nodes and 414 nodes, are used to discretize the semisphere shell surfaces. In Fig. 8, quadratic boundary condition corresponding to Eq. (27) is imposed on the body surfaces. And the numerical results obtained from 254 nodes and 414 nodes are good agreement with the exact solutions, even if the first data from 94 nodes has little fluctuation. In Fig. 9, quadratic boundary condition corresponding to Eq. (7) is imposed on the body. And the numerical results finally are convergent to the exact solutions by few of nodes. The two figures are testify that the PB-spline is a stable approximation method to approximate boundary variables in elasticity BIE. Stress Szz -5 -6 4.2. Linear elasticity problem for a tube -7 Quadratic exact solution 94 nodes 254 nodes 414 nodes -8 -9 -10 -11 0.0 0.2 0.4 0.6 0.8 1.0 Parametric coordinates Fig. 9. Variation of the stress Szz, along the red circle in Fig. 6(c), solved by cubic PB-spline. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article). Fig. 10. A tube and its main size. in the following illustration. In this case, Dirichlet boundary conditions corresponding to the exact solutions (Eqs. (26,27)) are imposed on the shell surface. The B-spline and PB-spline are applied for approximation of boundary variables. The surface in Fig. 7(a) meshed according to the parametric definition of bivariate B-spline. As a result, many of nodes huddle in the poles of the two semisphere surfaces that In this case, a tube is used for analysis, which is encompassed by four spline surfaces: two toroidal and two disc surfaces. The tube and its main size is shown in Fig. 10. The Dirichlet boundary conditions corresponding to the exact solutions (Eqs. 26 and 27) are imposed on the tube surface. There are three sets of nodes: (a) 276 nodes; (b) 526 nodes and (c) 892 nodes; to be used for discretization of the tube surfaces, these mesh models are presented in Fig. 11. The bivariate quadratic and cubic PB-spline are applied for approximation of boundary variables. The nodal stress components are denoted as Sxx, Syy, Sxy and Szz. Their L2 errors evaluated by Eq. (25) are denoted as Err_Sxx, Err_Syy and Err_Szz. Table 2 lists the L2 errors of nodal stress obtained by solving the linear elasticity BIEs. The table indicates that the numerical results computed with quadratic and cubic PB-spline all have high accuracy for linear displacement boundary condition (Eq. 26). It also testify that the IGABEM based on PB-spline is feasible. In order to investigate the stability of the PB-spline algorithm implemented in IGABEM with quadratic displacement boundary condition (Eq. 27), the numerical values of stress components Sxy and Szz on the toroidal and disc surfaces are observed. Fig. 12 shows a group of numerical Szz along the red dashed line on the toroidal surface of the tube. The numerical results are finally convergent to the quadratic exact solutions. Fig. 13 shows a group of numerical Sxy sampled from the red dashed line on the tube end surface. Note that the size of the end surface is much smaller than the tube size. The radial size of the face is only 0.25, but the entire size of the tube is more than 7. So the end surface here can be considered as a small size feature. To evaluate the stress values accurately, we need to use lots of elements to discretize the small size feature. Unfortunately, there are many singular and nearly singular integrals which need evaluation. However, the numerical results of Sxy are convergent to the quadratic exact solutions depicted by the Fig. 13. Especially for 892 nodes discretization, the group of numerical Sxy are a good agreement with the exact solutions. This numerical example testifies again that the 160 J. Gu et al. / Engineering Analysis with Boundary Elements 56 (2015) 154–161 Fig. 11. Three mesh models of the tube: (a) 276 nodes; (b) 526 nodes and (c) 892 nodes. -7.0 Table 2 The L2 errors of stress components on the tube Degree u ¼ Linear Quadratic Cubic Err_ Sxx (%) Err_ Syy (%) Err_ Szz (%) Err_ Sxx (%) Err_ Syy (%) Err_ Szz (%) 892nodes 526 nodes 276nodes 0.1633 0.1623 0.09876 0.1509 0.1504 0.07763 0.5291 0.5276 0.4362 0.3502 0.3508 0.2291 0.2549 0.2549 0.1601 0.2116 0.2115 0.1051 Quadratic exact solution 276 nodes 526 nodes 892 nodes -8.0 -8.5 -9.0 Stress Sxy Analytical field -7.5 -9.5 -10.0 -10.5 -11.0 7.25 -11.5 7.00 -12.0 6.75 -12.5 -0.1 6.50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 Parametric coordinates 6.25 Stress Szz 0.0 Fig. 13. Variation of the Stress Sxy,along red circle on the tube end surface,solved by quadratic PB-spline. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article). 6.00 5.75 Quadratic exact solution 276 nodes 526 nodes 892 nodes 5.50 5.25 5.00 4.75 4.50 0.0 0.2 0.4 0.6 0.8 1.0 Parametric coordinates Fig. 12. Variation of the Stress Szz, along red circle on the tube toroidal surface, solved by quadratic PB-spline. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article). numerical results converge to exact solution, and PB-spline has stronger potential in approximation with high accuracy. 5. Conclusions The conception of the IGA has been successfully implemented in the 3-D linear elasticity BIEs taking use of PB-spline. It should be a meaningful attempt to seamlessly integrate CAD modeling and BEM and to simplify the communication of mesh refinement with CAD model. The first goal of the IGA has been achieved in the linear elasticity BIEs, which is to be geometrically exact no matter how coarse the discretization is. PB-spline has been applied for CAD modeling and analysis of 3D linear elasticity BIEs. For some bodies, their NURBS mesh model produce many of extraordinary points restricted by the gobal tensor product of NURBS. The PB-spline has been proposed to avoid the computation of many nearly singular and singular integrals on these points. In approximation of boundary variables, the bivariate PB-spline functions are fitting type functions, i.e. they lack the Kronecker delta property, so an inverse transformation is performed to convert them into ones of interpolation type. Numerical tests have demonstrated that the IGABEM for solution of the 3-D linear elasticity BIEs is feasible. The satisfactory accuracy and convergence obtained from the tests indicate that the IGABEM has the stronger potential to be expanded into other areas of BEM. Because the coefficient matrix of the BIEs in IGABEM is full, numerical solution of the BIEs for complicated geometric body is a large-scale computation problem. To deal with the problem, the Fast Multipole Method (FMM) [22–24] can be combined with the IGABEM for reduction of the computation expense. And this is planned for our further research. References [1] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput. Methods Appl. Mech. Eng. 2005;194:4135–95. [2] Cottrell JA, Hughes TJR, Reali A. Studies of refinement and continuity in isogeometric analysis. Comput. Methods Appl. Mech. Eng. 2007;196:4160–83. [3] Cottrell JA, Reali A, Bazilevs Y, Hughes TJR. Isogeometric analysis of structural vibrations. Comput. Methods Appl. Mech. Eng. 2006;195:5257–96. [4] Bazilevs Y, Calo VM, Zhang Y, Hughes TJR. Isogeometric fluid-structure interaction analysis with applications to arterial blood flow. Comput. Methods Appl. Mech. Eng. 2006;38:310–22. J. Gu et al. / Engineering Analysis with Boundary Elements 56 (2015) 154–161 [5] Bazilevs Y, de Veiga L Beirao, Cottrell JA, Hughes TJR, Sangalli G. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes. Math. Model. Methods Appl. Sci. 2006;16:1031–90. [6] Zhang Y, Bazilevs Y, Goswami S, Bajaj C, Hughes TJR. Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow. Comput. Methods Appl. Mech. Eng. 2007;196:2943–59. [7] Sederberg TW, Zheng J, Bakenov A, Nasri A. T-splines and T-NURCCSs. ACM Trans. Graph. 2003;22:477–84. [8] Politis C, Ginnis AI, Kaklis PD, Belibassakis K, Feurer. C. An isogeometric BEM for exterior potential-flow problems in the plane. 2009 SIAM/ACM. Proceedings of the joint conference on geometric and physical modeling, SPM’09 2009:349–54. [9] Li K, Qian X. Isogeometric analysis and shape optimization via boundary integral. Comput. Aided Des. 2011;43(11):1427–37. [10] Simpson RN, Bordas SPA, Trevelyan J, Rabczk. T. A two-dimensional isogeometric boundary element method for elastostatic analysis. Comput. Methods Appl. Mech. Eng. 2012;209–212:87–100. [11] Scott MA, Simpson RN, Evans JA, Lipton S, Bordas SPA, Hughes TJR, et al. Isogeometric boundary element analysis using unstructured T-splines. Comput. Methods Appl. Mech. Eng. 2013;254:197–221. [12] Toru Takahashi, Toshiro Matsumoto An application of fast multipole method to isogeometric boundary element method for Laplace equation in two dimensions. Eng. Anal. Bound. Elem. 2012;36:1766–75. [13] Gavankar P, Henderson MR. Graph-based extraction of protrusions and depressions from boundary representation. Comput. Aided Des. 1990;22 (7):442–50. 161 [14] Falcidieno B, Giannini F. Automatic recognition and representation of shapebased features in a geometric modeling system. Comput. Vis., Graph., Image Process. 1989;48:93–123. [15] Zhang JM, Qin XY, Han X, Li GY. A boundary face method for potential problems in three dimensions. Int. J. Numer. Methods Eng. 2009;80:320–37. [16] Bazilevs Y, Calo VM, Cottrell JA, Evans JA, Hughes TJR, Lipton S, et al. Isogeometric analysis using T-splines. Comput. Methods Appl. Mech. Eng. 2010;199:229–63. [17] Sederberg TW, Zheng J, Bakenov A, Nasri A. T-splines and T-NURCCSs. ACM Trans. Graph. 2003;22(3):477–84. [18] Sederberg TW, Cardon DL, Finnigan GT, North NS, Zheng J, Lyche T. T-spline simplification and local refinement. ACM Trans. Graph 2004;23:276–83. [19] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques-theory and applications in engineering. Berlin: Springer; 1984. [20] Hong Hong-Ki. Jeng-Tzong chen. J. Eng. Mech. 1988;114:1028–44. [21] Jianming Zhang Masataka Tanaka. Toshiro Matsumoto. Meshless analysis of potential problems in three dimensions with the hybrid boundary node method. Int. J. Numer. Methods Eng. 2004;59:1147–66. [22] Zhang JM, Tanaka M, Endo M. The hybrid boundary node method accelerated by fast multipole method for 3D potential problems. Int. J. Numer. Methods Eng. 2005;63:660–80. [23] Zhang JM, Tanaka M. Fast HdBNM for large-scale thermal analysis of CNTreinforced composites. Comput. Mech. 2008;41:777–87. [24] Zhang JM, Tanaka M. Adaptive spatial decomposition in fast multipole method. J. Comput. Phys. 2007;226:17–28.