- Email: [email protected]
A two-dimensional Isogeometric Boundary Element Method for elastostatic analysis
Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
Contents lists available at SciVerse ScienceDirect
Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma
A two-dimensional Isogeometric Boundary Element Method for elastostatic analysis R.N. Simpson a, S.P.A. Bordas a,⇑, J. Trevelyan b, T. Rabczuk c a
School of Engineering, Institute of Mechanics and Advanced Materials, Cardiff University, Queen’s Buildings, The Parade, Cardiff CF24 3AA, United Kingdom School of Engineering and Computing Sciences, Durham University, South Road, Durham DH1 3LE, United Kingdom c Institute of Structural Mechanics, Bauhaus-University Weimar, Marienstraße 15, 99423 Weimar, United Kingdom b
a r t i c l e
i n f o
Article history: Received 11 April 2011 Received in revised form 19 July 2011 Accepted 9 August 2011 Available online 19 August 2011 Keywords: Isogeometric analysis Boundary Element Method NURBS
a b s t r a c t The concept of isogeometric analysis, where functions that are used to describe geometry in CAD software are used to approximate the unknown fields in numerical simulations, has received great attention in recent years. The method has the potential to have profound impact on engineering design, since the task of meshing, which in some cases can add significant overhead, has been circumvented. Much of the research effort has been focused on finite element implementations of the isogeometric concept, but at present, little has been seen on the application to the Boundary Element Method. The current paper proposes an Isogeometric Boundary Element Method (BEM), which we term IGABEM, applied to two-dimensional elastostatic problems using Non-Uniform Rational B-Splines (NURBS). We find it is a natural fit with the isogeometric concept since both the NURBS approximation and BEM deal with quantities entirely on the boundary. The method is verified against analytical solutions where it is seen that superior accuracies are achieved over a conventional quadratic isoparametric BEM implementation. 2011 Elsevier B.V. All rights reserved.
1. Introduction In the conventional engineering design process, geometries from CAD software are converted into suitable meshes for numerical analysis in a procedure commonly referred to as pre-processing. In some cases this task may take up the majority of the analysis time, and, for this reason, any improvements or speedups in the meshing procedure are of great industrial interest. A major step towards achieving this goal was presented in the seminal paper by Hughes et al. [9] introducing the idea of isogeometric analysis. The concept stems from the mismatch between CAD software, where geometries are most commonly described using functions known as Non-Uniform Rational B-Splines (NURBS) while in numerical analysis software, most often polynomial functions are used to describe the geometry and unknown fields. In contrast, isogeometric analysis uses the same functions that are used in CAD software to approximate the unknown fields in the numerical analysis stage. Of course, this comes at the expense of having to evaluate functions of greater complexity than simple polynomial shape functions, but efficient algorithms exist (see [14]) and efficiency is preserved. More recently, different functions other than NURBS have been used as an approximation advancing the efficiency of the method even further. These include T-splines [4,23,22] and PHT-splines ⇑ Corresponding author. Tel.: +44 0 7956520270. E-mail addresses: [email protected] (R.N. Simpson), stephane.bordas@ alumni.northwestern.edu (S.P.A. Bordas), [email protected] (J. Trevelyan). 0045-7825/$ - see front matter 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2011.08.008
[6,13] which have the attractive property that refinement can be achieved locally compared to the expensive global refinement achieved with NURBS due to the tensor product nature of the approximation. In addition, these approximations provide a solution to the formation of gaps between NURBS patches which often require ‘repairing’ by designers. Much of the recent literature on isogeometric analysis focuses on application to the Finite Element Method (FEM) with little mention given to other methods, with the exception of isogeometric collocation methods [2]. One method which has received little attention is the Boundary Element Method. In this paper we introduce an isogeometric BEM (IGABEM) in which NURBS are used to approximate the geometry along with the displacement and traction fields around the boundary. First, an overview of NURBS functions is outlined in which emphasis is given to their flexibility in modelling complex geometries. The isogeometric method and how it can be applied to BEM is then shown where the differences between a conventional BEM implementation are clearly stated. Finally, to verify the IGABEM implementation, convergence studies are carried out on two-dimensional problems with analytical solutions while comparisons with the conventional BEM are made throughout. 2. Problem definition All examples in the paper are based on linear elastostatic analysis, and to avoid any ambiguity, the problem is defined here in a formal manner. The equations of equilibrium are stated as:
88
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
rij;j þ bi ¼ 0; i; j ¼ 1; 2; 3;
ð1Þ
where rij is the Cauchy stress tensor and bi are body-force components. Assuming small displacements, strains eij are defined by
1 ui;j þ uj;i ; 2
eij ¼
i; j ¼ 1; 2; 3;
ð2Þ
where ui represents a displacement component and a comma implies differentiation. Assuming that the material is elastic, homogenous and isotropic, Eqs. (1) and (2) are related by Lamé’s equation:
rij ¼ kekk dij þ 2leij ; i; j; k ¼ 1; 2; 3;
ð3Þ
where a repeated index implies summation, the Kronecker-delta function dij is defined by
dij ¼
0
i–j;
1 i ¼ j;
ð4Þ
the Lamé constant k is given by
k¼
2ml ; ð1 2mÞ
ð5Þ
and l is the shear modulus:
l¼
E ; 2ð1 þ mÞ
ð6Þ
where E denotes Young’s modulus of elasticity and m denotes Poisson’s ratio. In addition to Eq. (2), strains have to satisfy the following compatibility equations:
oejk oeik oeij oeij o ¼ 0; þ þ oxj oxk oxi oxi oxj oxk
ð7Þ
where i, j, k = 1, 2, 3 and i – j – k. We can now define our problem. Let X denote an open set with boundary C as shown in Fig. 1 subject to the following boundary conditions:
i on Cu ; ui ¼ u t i ¼ t i on Ct ;
ð8Þ
Before any consideration is given towards isogeometric BEM and its implementation, some discussion on the construction of B-splines and NURBS must be made, since any isogeometric method is based heavily on such approximations. It is important to understand that B-splines and NURBS are parametric functions of the form F(n) in which the parameter n lies in the parameter space. This is analogous to the parent coordinate system that is often used to construct shape functions in FEM and BEM. We also note that in fact, NURBS are based on B-spline approximations, and as we shall see later, they differ only by applying a ‘‘projection’’ allowing complex geometries to be formed. For this reason, we first outline how B-spline functions are formed while making any definitions or terminology clear along the way. Then, NURBS are described in a similar fashion where the benefits of using such approximations over B-splines are made clear. 3.1. B-spline interpolation One of the most important elements of B-splines and NURBS is the knot vector, which is defined as a non-decreasing sequence of coordinates defined in the parameter space [14]. For example, the following is a valid knot vector: N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}. In addition, this knot vector is known as an open knot vector in which knots are repeated p + 1 times at the start and end of the vector where p is defined as the order of the spline (p = 0, 1, 2, . . . for constant, linear, quadratic functions etc.). Open knot vectors are ubiquitous in CAD software and as such, it can be assumed that in all future examples open knot vectors are used. Therefore, it is found that the knot vector is made up of a total of n + p + 1 components, where n is defined as the number of basis functions which make up the B-spline. These functions are denoted by Na,p with 1 6 a 6 n and are defined as follows for p = 0:
Na;0 ¼
ð9Þ
1 if na 6 n < naþ1 ; 0
otherwise;
ð13Þ
and for p = 1, 2, 3 . . .:
where
C ¼ Cu [ Ct ;
3. B-splines and NURBS
Cu \ Ct ¼ ;:
ð10Þ
Defining the closed set X as:
X ¼ X [ C;
ð11Þ
i ; t i and bi, find ui ðyÞ; y 2 X, subject to the the problem is: given u equilibrium condition of (1). Further to this, in all examples in the present paper plane strain is assumed and therefore:
ezz ¼ exz ¼ eyz ¼ 0:
Na;p ¼
naþpþ1 n n na Na;p1 ðnÞ þ Naþ1;p1 ðnÞ; naþp na naþpþ1 naþ1
ð14Þ
By inspecting Eq. (14), it can be seen that this is a recursive relation requiring the evaluation of lower order basis functions. This comes at a computational cost, but efficient algorithms exist [14,1]. Using these basis functions, it is possible to interpolate a set of coordinates Pa, commonly referred to as control points, in the following manner:
ð12Þ CðnÞ ¼
n X
Na;p ðnÞPa :
ð15Þ
a¼1
This interpolation is heavily dependent on the knot vector and to illustrate this, an example B-spline along with its associated control points is shown in Fig. 2 with p = 2 and N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}. The basis functions for this curve are also illustrated in Fig. 3. The curve is interpolatory at the ends of the curve and also at n = 3, which can be explained by the following feature of B-spline curves: If a knot if repeated k times then the continuity of the B-spline curve at that point will be Cpk. If the curve is C0 continuous at a point, then the curve is interpolatory at that point.
Fig. 1. Definition of boundary value problem.
Therefore, the curve is C1 continuous at n = 0, 5 and C0 continuous at n = 3. In addition, it is interpolatory at each of these locations. We find that at the corresponding locations in the physical
89
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
the first derivatives of the basis functions shown in Fig. 3 are plotted in Fig. 4. 3.3. NURBS interpolation
Fig. 2. B-spline curve and N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}.
associated
control
points
using
knot
vector
We now introduce NURBS where, before they are formally defined, we consider their benefits over B-splines. An excellent example is the case of a quarter-circle which presents a challenge for B-splines with few control points as shown in Fig. 5. Using three control points, a large discrepancy is seen between the interpolated curve and the exact. If the circle is now approximation using NURBS, we find that the circle can be interpolated exactly using three control points as shown in Fig. 6). This is achieved through the use of certain ‘weightings’ which, in the case of a circle can be calculated analytically [17] and can be considered an additional coordinate. Therefore, denoting ds as the number of space dimensions, NURBS operate in Rds þ1 space. To interpolate a set of control points Pa using NURBS, the following relation is used:
CðnÞ ¼
n X
Ra;p ðnÞPa ;
ð18Þ
a¼1
where the basis functions Ra,p(n) are defined by projecting the Bspline basis functions given by (13) and (14) using a set of weights wa as:
Na;p ðnÞwa Ra;p ðnÞ ¼ Pn : ^;p ðnÞwa ^ ^¼1 N a a
Fig. 3. B-spline basis N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}.
functions
for
p=2
and
knot
vector
space, the control points lie on the curve, but for any other control point this is not necessarily true. This is in contrast to traditional polynomial approximations used in BEM and FEM routines where the approximation is interpolatory at every nodal point. The lack of the Kronecker-delta property is also commonly seen in meshless methods (see [12] for a recent review).
ð19Þ
In the case that all weights are equal to one (wa = 1 " a), the basis functions given by (19) degenerate to the B-spline basis functions as expressed by (14). Therefore B-splines can be considered as a subset of NURBS. In keeping with the previous section on B-splines, it is instructive to plot the basis functions corresponding to a general NURBS curve. Using the same control points as in Fig. 2, but applying weightings of w4 = 3 and w6 = 3, with all other points given by wa = 1, a NURBS curve is defined as shown in Fig. 7. It can be seen that by applying stronger weights to points four and six the curve is drawn towards the associated control points. Inspecting the basis functions (Fig. 8) and comparing to the equivalent B-spline basis functions (3), the additional weighting given to R4,2 and R6,2 is evident.
3.2. B-spline derivatives In any FEM or BEM implementation the derivatives of the basis functions used in the approximations of the unknown fields are required. It can be shown that the first derivative of (14) is given by
d p p Na;p ðnÞ ¼ Na;p1 ðnÞ Naþ1;p1 ðnÞ; dn naþp na naþpþ1 naþ1
ð16Þ
with higher order derivatives given by the general formula: ðk1Þ
ðkÞ Na;p
¼p
Na;p1 naþp na
ðk1Þ
Naþ1;p1 naþpþ1 naþ1
! ;
ð17Þ
where NðkÞ a;p denotes the kth derivative of the basis function. Once again, it should be noted that these formulae are recursive in nature with an associated computational cost. To illustrate these formulae,
Fig. 4. B-spline basis function N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}.
derivatives
for
p=2
and
knot
vector
90
Fig. 5. Approximation N = {0, 0, 0, 1, 1, 1}.
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
of
quarter-circle
using
B-spline.
Knot
vector
Fig. 8. NURBS basis w = {1, 1, 1, 3, 1, 3, 1, 1}.
functions
for
p = 2,
N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}
and
dependence on the previously stated derivative expressions of (16) and (17). It can be shown that the first derivative of Eq. (19) is given by
WðnÞN0a;p W 0 ðnÞNa;p ðnÞ d Ra;p ðnÞ ¼ wa ; dn WðnÞ2
ð20Þ
with
WðnÞ ¼
n X
Na^;p ðnÞwa ;
ð21Þ
^¼1 a
N0a;p
d Na;p ; dn
ð22Þ
and
W 0 ðnÞ ¼
n X
N0a^;p ðnÞwa :
ð23Þ
^¼1 a
Fig. 6. Approximation of quarter-circle using NURBS. Knot vector N = {0, 0, 0, 1, 1, 1} pffiffiffi with weights w1 ¼ 1; w2 ¼ 1= 2 and w1 = 1.
Any higher order derivatives are determined using the expression shown in Appendix A. Using Eq. (20), the first order derivatives of the basis functions seen in Fig. 8 can be determined and are illustrated in Fig. 9. Comparing to Fig. 4, the effect of the weighting has an obvious effect and consequently, any numerical integration routine which evaluates the integral of these functions must take account of this. Hughes et al. [10] presented efficient rules to eval-
Fig. 7. B-spline basis function derivatives for p = 2, N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5} and w = {1, 1, 1, 3, 1, 3, 1, 1}.
3.4. NURBS derivatives Just as derivatives of B-splines are required for any isogeometric implementation, the same is also true for NURBS where we find a
Fig. 9. NURBS basis function derivatives for p = 2, N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5} and w = {1, 1, 1, 3, 1, 3, 1, 1}.
91
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
uate such integrands, although in most scenarios, Gauss–Legendre (GL) quadrature is sufficiently accurate. 4. Isogeometric BEM We now turn our attention towards how the unknown fields in a boundary value problem in elastostatics can be approximated using the NURBS approximation given by (18). This section therefore develops the main idea of the paper: applying the isogeometric concept within a BEM framework. We emphasise that this is a natural combination, in that the NURBS-based boundary representation (or ‘‘b-rep’’) of an object in a CAD system provides directly both the boundary geometry upon which the BEM operates and the approximation space in which it seeks its solution. No extra approximation is made in the interior of the object, either in the geometry, the approximation space, or the physics of the problem, as long as a fundamental solution is available (See Fig. 10). We organise the section as follows: first, the displacement boundary integral equation (DBIE) which forms the basis of the BEM is defined along with commonly used BEM notation. Next, the NURBS approximation of both the displacement and traction fields is given, the collocation positions are defined, and the methods used to evaluate all boundary integrals are described. Finally, an overview of the solution procedure is described to allow for numerical implementation.
Fig. 11. Definition of problem domain and boundary with source and field points.
Z C
t i ui dC þ
Z X
bi ui dX ¼
Z C
t i ui dC þ
Z X
bi ui dX:
ð25Þ
For simplicity in this article we will neglect body forces, i.e. bi = 0 to remove one of the volume integrals.1 The other may be removed by an appropriate selection of the (still arbitrary) stress state rij ; eij . It is most common in the BEM development for elasticity to take this state to result from an infinite point force injected through the body force term, i.e.
bi ¼ Dðx X0 Þei ;
ð26Þ 0
4.1. Boundary integral equation To allow understanding of later sections, a few key concepts of BEM are given presently where more comprehensive derivations of the method are given by several other authors e.g. [5,15]. We begin by defining the domain of the boundary value problem as X with a boundary C = oX. In addition, we define two points x0 2 C and x 2 C known as the ‘source point’ and ‘field point’ respectively, and the distance between them r: = kx x0 k. These are all depicted in Fig. 11. One of the most common ways of deriving the BEM is to utilise Betti’s reciprocal theorem which states that, given two states rij, eij and rij ; eij , then these can be related (assuming linear elasticity) as:
Z X
ij ij d
r e X¼
Z X
ij ij d
r e X;
where D is the Dirac delta function, X is a point within the domain (X0 2 X, X0 R C) and ei are the unit vector components. The solution for such a problem (of a point force within an infinite domain) is commonly referred to as Kelvin’s solution and solves for the components ui and ti as:
ui ¼ U ij ej
ð27Þ
where Uij and Tij are called ‘fundamental solutions’. Appendix B gives explicit expressions for the fundamental solutions for 2D elasticity problems. If these expressions are now substituted into (25) along with the property of the Dirac delta function given by R 0 0 X Dðx X Þuj dX ¼ uj ðX Þ, the following boundary integral equation can be written:
uj ðX0 Þ þ
Z
T ij ðX0 ; xÞuj ðxÞ dCðxÞ ¼
C
ð24Þ
Applying the divergence theorem and using equilibrium equations of (1), the displacement–strain relation of (2), and the definition of tractions given by ti = rijnj, Eq. (24) can be expressed as:
t i ¼ T ij ej ;
Z
U ij ðX0 ; xÞt j ðxÞ dCðxÞ;
ð28Þ
C
where the dependence on the source and field points has been made explicit. This equation consists almost entirely of quantities specified on the boundary except that the point X0 is restricted to the interior of the domain. Considering the limit as this point is taken toward the boundary, a jump term Cij appears and the first integral in (28) must be evaluated in Cauchy Principal Value (CPV) limiting sense. This is written as:
Z C ij ðx0 Þuj ðx0 Þ þ -- T ij ðx0 ; xÞuj ðxÞ dCðxÞ C Z 0 ¼ U ij ðx ; xÞt j ðxÞ dCðxÞ; i; j ¼ 1; 2;
ð29Þ
C
R where the integral sign -- denotes that the integrand is a CPV integral and all terms are now specified on the boundary. This equation defines the displacement boundary integral equation. To allow for numerical implementation, we discretise the boundary into a non-overlapping set of Ne elements giving:
C¼
Ne [
Ce ;
Ci \ Cj ¼ 0;
i–j:
ð30Þ
e¼1
Fig. 10. NURBS pffiffiffi basis functions corresponding to knot vector N = {0, 0, 0, 1, 1, 1} and w ¼ f1; 1= 2; 1g for quarter circle geometry.
1 This restriction is made for simplicity in this article, but there are various formulations to allow for body forces in BEM implementations. Body forces may easily be included in the isogeometric BEM framework.
92
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
This allows (29) to be rewritten:
C ij ðx0 Þuj ðx0 Þ þ
Ne X m Z X e¼1 l¼1
¼
Ne X m Z X e¼1 l¼1
þ1
1
þ1 1
T ij x0 ; x ^n Nl ^n J e ^n d^n uelj
U ij x0 ; x ^n Nl ^n J e ^n d^n t elj ;
ð31Þ
where l is a local node number, on element e, that varies from 1 to m = 2, 3,. . . for linear, quadratic elements etc., Nl is the local shape function for node l; Je ¼ dC=d^ n is the Jacobian of transformation, ^ n 2 ½1; 1 is a local coordinate system and uelj and telj are displacements and tractions, respectively, at local node l on element e. Fig. 12 illustrates a discretised domain and an example quadratic boundary element. To arrive at a system of equations, the conventional procedure is to collocate at a series of points around the boundary. This corresponds to placing the point x0 at a discrete set of points; most often, for classical (i.e. unenriched) BEM formulations the nodal positions are chosen since this automatically provides a square set of equations. In this way, by grouping the first two terms on the left-hand side of (31) the system of equations can be written in matrix notation as:
½Hfug ¼ ½Gftg;
Fig. 13. Example IGABEM domain created using NURBS with knot vector N = {0, 0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 9, 9}.
ð32Þ
where H is a square matrix containing a combination of the integrals of the Tij kernel and the jump terms, u is vector of displacement components, G is a rectangular matrix of Uij kernel integrals and t is a vector of traction components. Rearranging this set of equations by placing all unknown components on the left-hand side and known components on the right, the following can be written:
½Afxg ¼ fbg;
ð33Þ
where the vector x contains all unknown displacement and traction components. Eq. (33) may then be solved using an appropriate solver. Fig. 14. NURBS basis functions for IGABEM example.
4.2. Isogeometric approximation We now introduce the isoparametric concept to the displacement boundary integral equation where for clarity, a simple domain is used to illustrate each stage of the method. This is shown in Fig. 13 where NURBS were used to create the boundary and the first and last points are coincident to form a closed boundary. The knot vector used to create the curve is N = {0, 0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 9, 9} where the associated basis functions are illustrated in Fig. 14. To define the boundary elements, we take the unique values contained in the knot vector. Therefore, the element boundaries for this example are defined as {0, 1, 2, 3, 4, 5, 6, 7, 8, 9} in the parameter space (Fig. 14) which, when mapped onto the physical space, are depicted in Fig. 15. In
the isoparametric fashion, displacements and tractions are interpolated using the NURBS approximation (18) as:
uj ðnÞ ¼ tj ðnÞ ¼
n X a¼1 n X
Ra;p ðnÞdj ;
a
ð34Þ
Ra;p ðnÞqaj :
ð35Þ
a¼1 a
where dj and qaj represent displacement and traction coefficients respectively. Each of these is associated with a particular control point. These approximations are not used directly however, since the basis functions Ra,p(n) have a local support and should only be evaluated on those parts of the boundary where they are non-zero.
Fig. 12. BEM discretisation.
93
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
Fig. 16. Element containing collocation point in both global and local coordinate systems.
dC dn : Jð^nÞ ¼ dn d^n Fig. 15. Element definitions for IGABEM domain.
For example, in the case e = 1 with the element domain n = [0, 1], we can see from Fig. 14 that only the basis functions corresponding to a = 1, 2, 3 are non-zero. Therefore, it seems sensible to define elemental approximations of uj and tj as:
uj ðnÞ ¼
pþ1 X
Nel ðnÞdj ;
le
ð36Þ
Nel ðnÞqlej ;
ð37Þ
l¼1
t j ðnÞ ¼
pþ1 X l¼1
ð41Þ
Appendix C gives explicit expressions to evaluate this quantity. Some mention should be made of the first term in Eq. (40) relating to the jump term Cij which is distributed over the element e using the basis functions Nel . This is done because the local collocation 0 ^ coordinate n may not necessarily lie at an interpolatory point in the parameter space. This is in contrast to the conventional BEM formulation using polynomial shape functions where, by collocating at nodal coordinates, the basis functions are found to be interpolatory. To take account of this, we follow the same technique as used in [21,18] to distribute the jump term.
4.3. Collocation
where
Nel ðnÞ Ra;p ðnÞ;
ð38Þ
and the global basis function a is related to the local number l and element number e through a connectivity matrix:
a ¼ connðe; lÞ:
ð39Þ
The connectivity matrix for this example is shown in Table 1. Substituting Eqs. (36) and (37) into (29), the discretised isogeometric DBIE is given by 0
C ij ðx Þ
pþ1 X
Nel
pþ1 Z Ne X X ^n0 dle þ j e¼1 l¼1
l¼1
¼
pþ1 Z Ne X X e¼1 l¼1
þ1
U ij
1
þ1
1
le T ij x0 ; x ^n Nel ^n J ^n d^n dj
x0 ; x ^n Nel ^n J ^n d^n qlej ;
ð40Þ
where the local coordinate system ^ n ¼ ½1; 1 is used, e denotes the element containing the collocation point x0 ; ^ n0 denotes the local coordinate of the collocation point (see Fig. 16) and the Jacobian of transformation, Jð^ nÞ, is given by
Using the collocation method to arrive at a system of equations, a definition must be made to obtain the position of collocation points for isogeometric BEM. The solution is not immediately obvious, since if we follow the traditional procedure of collocating at nodal points (which correspond to control points in IGABEM), collocation will occur at points outside the boundary (see Fig. 13). In this study, we choose to collocate at the Greville abscissae [7,11] where the collocation points in parameter space are given by
n0a ¼
naþ1 þ naþ2 þ þ naþp ; p
a ¼ 1; 2; . . . ; n 1:
ð42Þ
Note that the index a varies from 1 to n 1 since we do not wish to collocate at the position corresponding to the start/end point twice. Applying this relation to the knot vector N = {0, 0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 9, 9} used in the previous section, we find the location of the collocation points in both parameter space and physical space are as shown in Figs. 17 and 18 respectively.
Table 1 Connectivity matrix a = conn (e, l) for IGABEM example. l
e
1 2 3 4 5 6 7 8 9
1
2
3
1 3 4 5 6 7 8 9 11
2 4 5 6 7 8 9 10 12
3 5 6 7 8 9 10 11 13
Fig. 17. Location of collocation points in parameter space for IGABEM example.
94
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
Fig. 18. Location of collocation points in physical space for IGABEM example.
4.4. Integration As is well-known in boundary element techniques, integration of singular integrands plays a key role in the implementation of any such method. The current work is no exception, and a few key points should be raised. Firstly, it is not possible to use rigid body motion [24,16] to evaluate the strongly singular integrals containing the Tij kernels and to see why this is the case, some further explanation is required. 4.4.1. Preclusion of rigid-body motion The system of equations given by (32) is constructed by considering the collocation point x0 to lie at each nodal point in turn, forming each row of the system of equations. In most cases we find that the collocation point does not lie in the same element we are integrating over (commonly referred to as the ‘field’ element), and the resulting Tij integrand is non-singular. But in the case when the collocation point does lie in the field element, we must use appropriate techniques to evaluate the resulting strongly singular Tij integrals. In the rigid-body technique all non-singular integrals for a particular collocation point are evaluated first, and rigid body displacements are then applied to evaluate the final remaining integrals. This is under the assumption that the number of singular integrals remaining is equal to the number of rigid body translations (e.g. x, y, z). Using polynomial shape functions and collocating at nodal points (as is the case in standard BEM), this is known to be true. This is because for a particular field element and collocation point, only the case where the collocation point and associated basis function coincide is singular – other integrands are found to be non-singular (see Fig. 19). In the case of NURBS however, we find that we cannot guarantee that the number of remaining singular integrals is equal to the number of rigid body translations. To see how this can arise, we consider a typical element that might be found in an isogeometric implementation. Fig. 20 illustrates an element with p = 3 and collocation points positioned according to Eq. (42). Considering collocation at the point n = 0 we find that only R1,3 is non-zero. But we can also show that the basis functions R2,3, R3,3 and R4,3 are of OðrÞ and therefore, as in the case of conventional BEM, only one set of singular integrals remain and rigid body motion can be applied. In the case of the second and third collocation point (n = 1/3, 2/3) we can see that all basis functions are non-zero, but more importantly, they are all of Oð1Þ. In this scenario we find that the four integrands corresponding to each of the basis functions are singular and rigid-body motion is insufficient. Explicit singular integration routines are therefore necessary. In the current implementation we choose to use a subtraction of singularity scheme as proposed by Guiggiani and Casalini [8]. This
Fig. 19. Example boundary element for case where source point lies in the field element but resulting integrand is non-singular.
Fig. 20. Example NURBS element with collocation at point where basis functions are found to be non-zero and of Oð1Þ.
creates a non-singular integrand that can be easily computed using GL quadrature and accounts for the singular part through an analytical integration. Weakly singular integrals are seen in the Uij integrands, but these do not present issues for implementation. One of the most convenient techniques to evaluate integrals of this type is the transformation as proposed by Telles [20] which removes the singularity by ensuring that the Jacobian of the transformation approaches zero at the singular point. The technique is simple to implement and efficient, and it is adopted in the implementation of IGABEM in the present study. 4.5. Solution procedure Using the previously outlined techniques for the formation of elements, collocation points and integration over boundary elements, it is possible to perform an IGABEM analysis as illustrated in Fig. 21. This is very similar to a conventional BEM analysis but with the following main differences:
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
95
Fig. 23. Mesh and basis functions for coarse mesh, p = 2.
Fig. 21. Solution procedure for IGABEM.
Fig. 22. Hole within an infinite plate - geometry, boundary conditions and material properties.
The mesh is now completely generated from CAD data. The definition of elements is given by unique knot vector values. NURBS basis functions are substituted in place of polynomial shape functions for both the geometry and unknown fields. Collocation occurs at points defined by the Greville abscissae. Explicit singular integration routines are used to calculate the Tij integrands – rigid body motion is no longer applicable.
Fig. 24. Mesh and basis functions for coarse mesh, p = 3.
5. Examples To illustrate the accuracy and flexibility of IGABEM, several 2D elastostatic examples are now presented. The first two problems of a hole within an infinite plate and an L-shaped wedge have
96
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
Fig. 25. Comparison of deformed profile with exact solution for three elements per line, p = 2.
Fig. 26. L2 relative error norm plot for plate with hole.
Fig. 28. IGABEM mesh and basis functions for L-plate problem for one element per line and p = 2.
Fig. 27. Geometry, boundary conditions and material properties for L-plate problem.
closed-form solutions and therefore a convergence study carried out on each. In addition, a comparison is made with an isoparametric quadratic BEM to assess the accuracy of IGABEM. Finally, the problem of an open spanner is illustrated to show the ability of the method to handle arbitrary CAD geometries. 5.1. Hole within an infinite plate The problem of a hole within an infinite plate subject to a traction Tx at infinity has an exact solution given by e.g. [3]. The problem is defined in Fig. 22 where symmetry is applied to the bottom and right edges and exact tractions are applied along the top and
left edges. Using the polar coordinate system (q, h) defined in Fig. 22, these tractions can be found using:
! ! Tx R2 Tx R2 R4 rqq ðq; hÞ ¼ 1 2 þ 1 4 2 þ 3 4 cos 2h; 2 q 2 q q ! ! T R2 T R4 rhh ðq; hÞ ¼ x 1 þ 2 x 1 þ 3 4 cos 2h; 2 q 2 q ! 2 4 T R R rqh ðq; hÞ ¼ x 1 þ 2 2 3 4 sin 2h; 2 q q
ð43Þ ð44Þ ð45Þ
where an appropriate transformation to cartesian coordinates and the relation:
ti ¼ rij nj
ð46Þ
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
97
is required. Choosing the order of the NURBS basis function to be quadratic (p = 2) in the first instance, the minimum number of control points required to represent the geometry of this problem exactly is n = 11, as shown in Fig. 23(a), where the first and last points are coincident. The knot vector in this case is found to be:
N ¼ f0; 0; 0; 1=9; 1=9; 3=9; 3=9; 6=9; 6=9; 7=9; 7=9; 1; 1; 1g;
ð47Þ
with the weights:
pffiffiffi wa ¼ f1; 1; 1; 1= 2; 1; 1; 1; 1; 1; 1; 1g a ¼ 1; 2; ; n:
ð48Þ
Using the definition of elements and collocation points as described in Sections 4.2 and 4.3, the IGABEM mesh for this set of control points and knot vector is shown in Fig. 23(b). The NURBS basis functions for p = 2 in parameter space are shown in Fig. 23(c). In
Fig. 30. Element and collocation point definitions for three elements per line and deformed profile.
Fig. 31. L2 relative error norm plot for L-plate problem.
Fig. 29. IGABEM mesh and basis functions for L-plate problem for one element per line and p = 3.
addition, to evaluate the effect of higher order NURBS basis functions, the geometry and unknown fields were described by setting p = 3 with the mesh and basis functions illustrated in Fig. 24. Fig. 29 illustrates the corresponding control points, elements, collocation points and resulting basis functions for the plate with a hole problem. The knot insertion algorithm (equivalent to h-refinement) as described in [9,14] was used to carry out a convergence study with
98
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
uniform refinement applied around the boundary for both p = 2 and p = 3. In addition, the problem was analysed using quadratic isoparametric boundary elements using an equivalent mesh refinement strategy. Exactly the same number of Gauss points were used to evaluate each of the boundary integrals given by the second and third terms in Eq. (31) for both the IGABEM and BEM analyses. Fig. 25 illustrates an IGABEM mesh with three elements per line and the deformed IGABEM profile. Excellent agreement with the analytical solution is seen. Using the following definition to calculate the relative L2 error norm in displacements around the boundary:
eL2 ¼
ku uex kL2 ; kuex kL2
ð49Þ
where
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uZ dp u X kukL2 ¼ t ðui Þ2 dC;
ð50Þ
C i¼1
a comparison can be made between IGABEM and BEM (Fig. 26). In the case of IGABEM with p = 2 and quadratic BEM, both methods converge at the same rate but importantly, IGABEM demonstrates a consistently lower error for all meshes. For IGABEM with p = 3,
we see, as expected, that a higher convergence rate is obtained with lower errors than the equivalent second order mesh. 5.2. L-shaped wedge The next problem which was considered was the L-shaped wedge which exhibits a singularity at the wedge apex. The analytical solution to this problem is given by Szabó and Babuška [19] where a wedge angle of 2a = 3p/2 was used. Considering only the mode 1 loading case, exact tractions were applied along all faces with appropriate displacement constraints as shown in Fig. 27. Material properties E = 1e5 and m = 0.3 were used under plane strain conditions. The problem was solved using four different methods: quadratic BEM with uniform h-refinement, IGABEM with p = 2 and uniform h-refinement, IGABEM with p = 3 and uniform h-refinement and finally IGABEM with p = 2 and graded h-refinement towards the wedge apex. For the case of one element per line and p = 2, the control points are shown in Fig. 28(a) with collocation points and elements shown in Fig. 28(b). The knot vector for this example is given by
N ¼ f0; 0; 0; 1=6; 1=6; 2=6; 2=6; 3=6; 3=6; 4=6; 4=6; 5=6; 5=6; 1; 1; 1g;
Fig. 32. Open spanner problem.
ð51Þ
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
99
Fig. 33. Von Mises stress for open spanner problem.
with wa = 1 " a and p = 2. The resulting NURBS basis functions (which in this case are B-splines) are shown in Fig. 28(c). Likewise, the control points, elements, collocation points and basis functions for p = 3 are shown in Fig. 29. Fig. 30(a) shows a mesh of three elements per line with the deformed IGABEM profile shown against the exact solution. As can be seen, good agreement can be seen, but the effect of the singularity at the wedge apex is clearly evident (cf. Fig. 25). Plotting the same deformed profile but using the graded mesh (with p = 2), we can visually verify that the apex singularity is captured more accurately. Most importantly though, Fig. 31 shows that comparing IGABEM against BEM, consistently higher accuracies are achieved for an equivalent number of degrees of freedom for both p = 2 and p = 3, with the graded mesh achieving the most accurate results of all the analyses. 5.3. Open spanner Finally, to illustrate the ability of IGABEM to handle arbitrary geometries which may be taken directly from CAD, the problem of an open spanner as illustrated in Fig. 32(a) was analysed using IGABEM. A uniform traction of ty = 10 is applied at the right edge with all other faces traction-free. The displacement boundary conditions are also shown in Fig. 32(a). The control points for this example are shown in Fig. 32(b) while Fig. 32(c) illustrates the collocation points and element edges. To allow for comparison, the problem was also solved using FEM. Fig. 32(c) shows the deformed mesh profile obtained using both FEM and IGABEM where good agreement is obtained. Von Mises stresses obtained from both methods are shown in Fig. 33 where good agreement is also seen. This example illustrates an important concept: through IGABEM, it is possible to integrate CAD and BEM into a seamless
design/analysis. The method only requires knot vectors and control points to completely define the mesh and these are taken directly from CAD. Therefore, it is possible to analyse an extremely wide variety of geometries (only restricted by the capabilities of NURBS) using the present method.
6. Conclusions A formulation for isogeometric boundary element analysis has been presented. The problem geometry and the approximation spaces for the displacement and traction solutions are all provided, using the isoparametric concept, from a set of NURBS basis functions. This offers a number of key improvements that together offer a significant advantage over conventional piecewise polynomial formulations. Firstly, the NURBS basis permits exact representation of commonly used geometries including circular arcs. Secondly, a meshing process is no longer required because the control points defining the surfaces play the role of nodes in the IGABEM analysis. Thirdly, the use of the NURBS basis has been shown to offer improved accuracy, for a given model size, in comparison with conventional piecewise quadratic BEM schemes. The combination of the isogeometric concept and the BEM is a natural and productive one, since both the NURBS boundary representation and the BEM deal with quantities entirely on the problem boundary. The authors expect this development to have a major significance in the future of computational mechanics with the deeper integration of stress analysis into CAD. Future developments include the extension of the method to three-dimensional problems and the application of the method to contact problems where it is expected that the use of NURBS will demonstrate benefits over conventional BEM implementations.
100
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
which allows the desired expression to be stated as:
Appendix A. Higher order NURBS derivatives Any higher order derivatives of the NURBS basis functions are given by the following general formula:
AaðkÞ ðnÞ
k W ðjÞ ðnÞRðkjÞ a;p ðnÞ j¼1 j ; WðnÞ
Pk
ðkÞ Ra;p ðnÞ ¼
where the superscript function and
(k)
ð52Þ
denotes the kth derivative of relevant
ðkÞ AðkÞ a ðnÞ ¼ wa N a;p ðnÞ;
ð53Þ
(no summation over index a) with
k j
¼
k! : j!ðk jÞ!
ð54Þ
Appendix B. 2D Fundamental solutions The fundamental solutions for Tij and Uij in plane strain are given by
U ij ðx0 ; xÞ ¼
1 1 ð3 4mÞ ln dij þ r ;i r ;j ; 8plð1 mÞ r
T ij ðx0 ; xÞ ¼
ð55aÞ
1 or ð1 2mÞ dij þ 2r ;i r;j ð1 2mÞ r;i nj r;j ni ; 4pð1 mÞr on ð55bÞ
with
r ¼ kx0 xk;
ð56Þ
lis the shear modulus, m is Poisson’s ratio, dij is the Kronecker-delta function, ni is the i-component of the outward-pointing normal at x and r,i or/oxi. For plane stress problems the fundamental solutions (55a) and (55b) may be used with an appropriately modified set of elastic constants. Appendix C. Transformation of coordinate system To transform the integrals into the local coordinate system ^ n 2 ½1; 1, the Jacobian given by (41) is required. This consists of two parts: dC/dn and dn=d^ n. The first is given by
dC ¼ dn
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 dx dy þ ; dn dn
ð57Þ
which can be determined using the derivatives of the NURBS basis functions (Section 3.4). To evaluate the second expression, we first assume that the element is defined by the interval in the parameter space as n 2 [n1, n2]. Then the following linear relationship between n and ^ n exists:
h n¼
i ðn2 n1 Þ^n þ ðn2 þ n1 Þ 2
;
ð58Þ
dn n2 n1 ¼ : 2 d^n
ð59Þ
References [1], 2011. [2] F. Auricchio, L.B. Da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli, Isogeometric collocation methods, Math. Models Methods Appl. Sci. 20 (11) (2010) 2075– 2107. [3] J.R. Barber, Elasticity, Kluwer Academic Publishers, 2002. [4] Y. Bazilevs, V.M. Calo, J.A. Cottrell, J.A. Evans, T.J.R. Hughes, S. Lipton, M.A. Scott, T.W. Sederberg, Isogeometric analysis using T-splines, Comput. Methods Appl. Mech. Engrg. 199 (5–8) (2010) 229–263. [5] T. Cruse, Numerical solutions in three dimensional elastostatics, Int. J. Solids Struct. 5 (1969) 1259–1274. [6] J. Deng, F. Chen, X. Li, C. Hu, W. Tong, Z. Yang, Y. Feng, Polynomial splines over hierarchical T-meshes, Graph. Models 70 (2008) 76–86. [7] T. Greville, Numerical procedures for interpolation by spline functions, J. Soc. Ind. Appl. Math. Ser. B. Numer. Anal. (1964). [8] M. Guiggiani, P. Casalini, Direct computation of cauchy principal value integrals in advanced boundary elements, Int. J. Numer. Methods Engrg. 24 (1987) 1711–1720. [9] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (39–41) (2005) 4135–4195. [10] T.J.R. Hughes, A. Reali, G. Sangalli, Efficient quadrature for NURBS-based isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 199 (5–8) (2010) 301–313. [11] R. Johnson, Higher order B-spline collocation at the Greville abscissae, Appl. Numer. Math. 52 (2005) 63–75. [12] V.P. Nguyen, T. Rabczuk, S. Bordas, M. Duflot, Meshless methods: a review and computer implementation aspects, Math. Comput. Simul. 79 (3) (2008) 763– 813. [13] N. Nguyen-Thanh, H. Nguyen-Xuan, S.P.A. Bordas, T. Rabczuk, Isogeometric analysis using polynomial splines over hierarchical T-meshes for twodimensional elastic solids, Comput. Methods Appl. Mech. Engrg. (2011) 1–37. [14] L. Piegl, W. Tiller, The NURBS Book, Springer, 1995. [15] F.J. Rizzo, An integral equation approach to boundary value problems of classical elastostatics, Q. Appl. Math. 25 (1) (1967) 83–95. [16] F.J. Rizzo, D.J. Shippy, An advanced boundary integral equation method for three-dimensional thermoelasticity, Int. J. Numer. Methods Engrg. 11 (1977) 1753–1768. [17] D.F. Rogers, An Introduction to NURBS: with Historical Perspective, Elsevier, 2001. [18] R. Simpson, J. Trevelyan, A partition of unity enriched dual Boundary Element Method for accurate computations in fracture mechanics, Comput. Methods Appl. Mech. Engrg. 200 (1–4) (2010) 1–10. [19] B. Szabó, I. Babuška, Finite Element Analysis, John Wiley and Sons, 1991. [20] J.C.F. Telles, A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary integrals, Int. J. Numer. Methods Engrg. 24 (1987) 959–973. [21] J. Trevelyan, G. Coates, On adaptive definition of the plane wave basis for wave boundary elements in acoustic scattering: the 2D case, Comput. Model. Engrg. Sci. 55 (2010) 147–170. [22] T.-K. Uhm, K.-S. Kim, Y.-D. Seo, S.-K. Youn, A locally refinable T-spline finite element method for CAD/CAE integration, Struct. Engrg. Mech. 30 (2) (2008) 225–245. [23] T.-K. Uhm, S.-K. Youn, T-spline finite element method for the analysis of shell structures, Int. J. Numer. Methods Engrg. 80 (4) (2009) 507–536. [24] J.O. Watson, Advanced Implementation of the Boundary Element Method for Two- and Three-Dimensional Elastostatics, Applied Science Publishers, London, 1979.
Contents lists available at SciVerse ScienceDirect
Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma
A two-dimensional Isogeometric Boundary Element Method for elastostatic analysis R.N. Simpson a, S.P.A. Bordas a,⇑, J. Trevelyan b, T. Rabczuk c a
School of Engineering, Institute of Mechanics and Advanced Materials, Cardiff University, Queen’s Buildings, The Parade, Cardiff CF24 3AA, United Kingdom School of Engineering and Computing Sciences, Durham University, South Road, Durham DH1 3LE, United Kingdom c Institute of Structural Mechanics, Bauhaus-University Weimar, Marienstraße 15, 99423 Weimar, United Kingdom b
a r t i c l e
i n f o
Article history: Received 11 April 2011 Received in revised form 19 July 2011 Accepted 9 August 2011 Available online 19 August 2011 Keywords: Isogeometric analysis Boundary Element Method NURBS
a b s t r a c t The concept of isogeometric analysis, where functions that are used to describe geometry in CAD software are used to approximate the unknown fields in numerical simulations, has received great attention in recent years. The method has the potential to have profound impact on engineering design, since the task of meshing, which in some cases can add significant overhead, has been circumvented. Much of the research effort has been focused on finite element implementations of the isogeometric concept, but at present, little has been seen on the application to the Boundary Element Method. The current paper proposes an Isogeometric Boundary Element Method (BEM), which we term IGABEM, applied to two-dimensional elastostatic problems using Non-Uniform Rational B-Splines (NURBS). We find it is a natural fit with the isogeometric concept since both the NURBS approximation and BEM deal with quantities entirely on the boundary. The method is verified against analytical solutions where it is seen that superior accuracies are achieved over a conventional quadratic isoparametric BEM implementation. 2011 Elsevier B.V. All rights reserved.
1. Introduction In the conventional engineering design process, geometries from CAD software are converted into suitable meshes for numerical analysis in a procedure commonly referred to as pre-processing. In some cases this task may take up the majority of the analysis time, and, for this reason, any improvements or speedups in the meshing procedure are of great industrial interest. A major step towards achieving this goal was presented in the seminal paper by Hughes et al. [9] introducing the idea of isogeometric analysis. The concept stems from the mismatch between CAD software, where geometries are most commonly described using functions known as Non-Uniform Rational B-Splines (NURBS) while in numerical analysis software, most often polynomial functions are used to describe the geometry and unknown fields. In contrast, isogeometric analysis uses the same functions that are used in CAD software to approximate the unknown fields in the numerical analysis stage. Of course, this comes at the expense of having to evaluate functions of greater complexity than simple polynomial shape functions, but efficient algorithms exist (see [14]) and efficiency is preserved. More recently, different functions other than NURBS have been used as an approximation advancing the efficiency of the method even further. These include T-splines [4,23,22] and PHT-splines ⇑ Corresponding author. Tel.: +44 0 7956520270. E-mail addresses: [email protected] (R.N. Simpson), stephane.bordas@ alumni.northwestern.edu (S.P.A. Bordas), [email protected] (J. Trevelyan). 0045-7825/$ - see front matter 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2011.08.008
[6,13] which have the attractive property that refinement can be achieved locally compared to the expensive global refinement achieved with NURBS due to the tensor product nature of the approximation. In addition, these approximations provide a solution to the formation of gaps between NURBS patches which often require ‘repairing’ by designers. Much of the recent literature on isogeometric analysis focuses on application to the Finite Element Method (FEM) with little mention given to other methods, with the exception of isogeometric collocation methods [2]. One method which has received little attention is the Boundary Element Method. In this paper we introduce an isogeometric BEM (IGABEM) in which NURBS are used to approximate the geometry along with the displacement and traction fields around the boundary. First, an overview of NURBS functions is outlined in which emphasis is given to their flexibility in modelling complex geometries. The isogeometric method and how it can be applied to BEM is then shown where the differences between a conventional BEM implementation are clearly stated. Finally, to verify the IGABEM implementation, convergence studies are carried out on two-dimensional problems with analytical solutions while comparisons with the conventional BEM are made throughout. 2. Problem definition All examples in the paper are based on linear elastostatic analysis, and to avoid any ambiguity, the problem is defined here in a formal manner. The equations of equilibrium are stated as:
88
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
rij;j þ bi ¼ 0; i; j ¼ 1; 2; 3;
ð1Þ
where rij is the Cauchy stress tensor and bi are body-force components. Assuming small displacements, strains eij are defined by
1 ui;j þ uj;i ; 2
eij ¼
i; j ¼ 1; 2; 3;
ð2Þ
where ui represents a displacement component and a comma implies differentiation. Assuming that the material is elastic, homogenous and isotropic, Eqs. (1) and (2) are related by Lamé’s equation:
rij ¼ kekk dij þ 2leij ; i; j; k ¼ 1; 2; 3;
ð3Þ
where a repeated index implies summation, the Kronecker-delta function dij is defined by
dij ¼
0
i–j;
1 i ¼ j;
ð4Þ
the Lamé constant k is given by
k¼
2ml ; ð1 2mÞ
ð5Þ
and l is the shear modulus:
l¼
E ; 2ð1 þ mÞ
ð6Þ
where E denotes Young’s modulus of elasticity and m denotes Poisson’s ratio. In addition to Eq. (2), strains have to satisfy the following compatibility equations:
oejk oeik oeij oeij o ¼ 0; þ þ oxj oxk oxi oxi oxj oxk
ð7Þ
where i, j, k = 1, 2, 3 and i – j – k. We can now define our problem. Let X denote an open set with boundary C as shown in Fig. 1 subject to the following boundary conditions:
i on Cu ; ui ¼ u t i ¼ t i on Ct ;
ð8Þ
Before any consideration is given towards isogeometric BEM and its implementation, some discussion on the construction of B-splines and NURBS must be made, since any isogeometric method is based heavily on such approximations. It is important to understand that B-splines and NURBS are parametric functions of the form F(n) in which the parameter n lies in the parameter space. This is analogous to the parent coordinate system that is often used to construct shape functions in FEM and BEM. We also note that in fact, NURBS are based on B-spline approximations, and as we shall see later, they differ only by applying a ‘‘projection’’ allowing complex geometries to be formed. For this reason, we first outline how B-spline functions are formed while making any definitions or terminology clear along the way. Then, NURBS are described in a similar fashion where the benefits of using such approximations over B-splines are made clear. 3.1. B-spline interpolation One of the most important elements of B-splines and NURBS is the knot vector, which is defined as a non-decreasing sequence of coordinates defined in the parameter space [14]. For example, the following is a valid knot vector: N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}. In addition, this knot vector is known as an open knot vector in which knots are repeated p + 1 times at the start and end of the vector where p is defined as the order of the spline (p = 0, 1, 2, . . . for constant, linear, quadratic functions etc.). Open knot vectors are ubiquitous in CAD software and as such, it can be assumed that in all future examples open knot vectors are used. Therefore, it is found that the knot vector is made up of a total of n + p + 1 components, where n is defined as the number of basis functions which make up the B-spline. These functions are denoted by Na,p with 1 6 a 6 n and are defined as follows for p = 0:
Na;0 ¼
ð9Þ
1 if na 6 n < naþ1 ; 0
otherwise;
ð13Þ
and for p = 1, 2, 3 . . .:
where
C ¼ Cu [ Ct ;
3. B-splines and NURBS
Cu \ Ct ¼ ;:
ð10Þ
Defining the closed set X as:
X ¼ X [ C;
ð11Þ
i ; t i and bi, find ui ðyÞ; y 2 X, subject to the the problem is: given u equilibrium condition of (1). Further to this, in all examples in the present paper plane strain is assumed and therefore:
ezz ¼ exz ¼ eyz ¼ 0:
Na;p ¼
naþpþ1 n n na Na;p1 ðnÞ þ Naþ1;p1 ðnÞ; naþp na naþpþ1 naþ1
ð14Þ
By inspecting Eq. (14), it can be seen that this is a recursive relation requiring the evaluation of lower order basis functions. This comes at a computational cost, but efficient algorithms exist [14,1]. Using these basis functions, it is possible to interpolate a set of coordinates Pa, commonly referred to as control points, in the following manner:
ð12Þ CðnÞ ¼
n X
Na;p ðnÞPa :
ð15Þ
a¼1
This interpolation is heavily dependent on the knot vector and to illustrate this, an example B-spline along with its associated control points is shown in Fig. 2 with p = 2 and N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}. The basis functions for this curve are also illustrated in Fig. 3. The curve is interpolatory at the ends of the curve and also at n = 3, which can be explained by the following feature of B-spline curves: If a knot if repeated k times then the continuity of the B-spline curve at that point will be Cpk. If the curve is C0 continuous at a point, then the curve is interpolatory at that point.
Fig. 1. Definition of boundary value problem.
Therefore, the curve is C1 continuous at n = 0, 5 and C0 continuous at n = 3. In addition, it is interpolatory at each of these locations. We find that at the corresponding locations in the physical
89
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
the first derivatives of the basis functions shown in Fig. 3 are plotted in Fig. 4. 3.3. NURBS interpolation
Fig. 2. B-spline curve and N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}.
associated
control
points
using
knot
vector
We now introduce NURBS where, before they are formally defined, we consider their benefits over B-splines. An excellent example is the case of a quarter-circle which presents a challenge for B-splines with few control points as shown in Fig. 5. Using three control points, a large discrepancy is seen between the interpolated curve and the exact. If the circle is now approximation using NURBS, we find that the circle can be interpolated exactly using three control points as shown in Fig. 6). This is achieved through the use of certain ‘weightings’ which, in the case of a circle can be calculated analytically [17] and can be considered an additional coordinate. Therefore, denoting ds as the number of space dimensions, NURBS operate in Rds þ1 space. To interpolate a set of control points Pa using NURBS, the following relation is used:
CðnÞ ¼
n X
Ra;p ðnÞPa ;
ð18Þ
a¼1
where the basis functions Ra,p(n) are defined by projecting the Bspline basis functions given by (13) and (14) using a set of weights wa as:
Na;p ðnÞwa Ra;p ðnÞ ¼ Pn : ^;p ðnÞwa ^ ^¼1 N a a
Fig. 3. B-spline basis N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}.
functions
for
p=2
and
knot
vector
space, the control points lie on the curve, but for any other control point this is not necessarily true. This is in contrast to traditional polynomial approximations used in BEM and FEM routines where the approximation is interpolatory at every nodal point. The lack of the Kronecker-delta property is also commonly seen in meshless methods (see [12] for a recent review).
ð19Þ
In the case that all weights are equal to one (wa = 1 " a), the basis functions given by (19) degenerate to the B-spline basis functions as expressed by (14). Therefore B-splines can be considered as a subset of NURBS. In keeping with the previous section on B-splines, it is instructive to plot the basis functions corresponding to a general NURBS curve. Using the same control points as in Fig. 2, but applying weightings of w4 = 3 and w6 = 3, with all other points given by wa = 1, a NURBS curve is defined as shown in Fig. 7. It can be seen that by applying stronger weights to points four and six the curve is drawn towards the associated control points. Inspecting the basis functions (Fig. 8) and comparing to the equivalent B-spline basis functions (3), the additional weighting given to R4,2 and R6,2 is evident.
3.2. B-spline derivatives In any FEM or BEM implementation the derivatives of the basis functions used in the approximations of the unknown fields are required. It can be shown that the first derivative of (14) is given by
d p p Na;p ðnÞ ¼ Na;p1 ðnÞ Naþ1;p1 ðnÞ; dn naþp na naþpþ1 naþ1
ð16Þ
with higher order derivatives given by the general formula: ðk1Þ
ðkÞ Na;p
¼p
Na;p1 naþp na
ðk1Þ
Naþ1;p1 naþpþ1 naþ1
! ;
ð17Þ
where NðkÞ a;p denotes the kth derivative of the basis function. Once again, it should be noted that these formulae are recursive in nature with an associated computational cost. To illustrate these formulae,
Fig. 4. B-spline basis function N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}.
derivatives
for
p=2
and
knot
vector
90
Fig. 5. Approximation N = {0, 0, 0, 1, 1, 1}.
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
of
quarter-circle
using
B-spline.
Knot
vector
Fig. 8. NURBS basis w = {1, 1, 1, 3, 1, 3, 1, 1}.
functions
for
p = 2,
N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5}
and
dependence on the previously stated derivative expressions of (16) and (17). It can be shown that the first derivative of Eq. (19) is given by
WðnÞN0a;p W 0 ðnÞNa;p ðnÞ d Ra;p ðnÞ ¼ wa ; dn WðnÞ2
ð20Þ
with
WðnÞ ¼
n X
Na^;p ðnÞwa ;
ð21Þ
^¼1 a
N0a;p
d Na;p ; dn
ð22Þ
and
W 0 ðnÞ ¼
n X
N0a^;p ðnÞwa :
ð23Þ
^¼1 a
Fig. 6. Approximation of quarter-circle using NURBS. Knot vector N = {0, 0, 0, 1, 1, 1} pffiffiffi with weights w1 ¼ 1; w2 ¼ 1= 2 and w1 = 1.
Any higher order derivatives are determined using the expression shown in Appendix A. Using Eq. (20), the first order derivatives of the basis functions seen in Fig. 8 can be determined and are illustrated in Fig. 9. Comparing to Fig. 4, the effect of the weighting has an obvious effect and consequently, any numerical integration routine which evaluates the integral of these functions must take account of this. Hughes et al. [10] presented efficient rules to eval-
Fig. 7. B-spline basis function derivatives for p = 2, N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5} and w = {1, 1, 1, 3, 1, 3, 1, 1}.
3.4. NURBS derivatives Just as derivatives of B-splines are required for any isogeometric implementation, the same is also true for NURBS where we find a
Fig. 9. NURBS basis function derivatives for p = 2, N = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5, 5} and w = {1, 1, 1, 3, 1, 3, 1, 1}.
91
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
uate such integrands, although in most scenarios, Gauss–Legendre (GL) quadrature is sufficiently accurate. 4. Isogeometric BEM We now turn our attention towards how the unknown fields in a boundary value problem in elastostatics can be approximated using the NURBS approximation given by (18). This section therefore develops the main idea of the paper: applying the isogeometric concept within a BEM framework. We emphasise that this is a natural combination, in that the NURBS-based boundary representation (or ‘‘b-rep’’) of an object in a CAD system provides directly both the boundary geometry upon which the BEM operates and the approximation space in which it seeks its solution. No extra approximation is made in the interior of the object, either in the geometry, the approximation space, or the physics of the problem, as long as a fundamental solution is available (See Fig. 10). We organise the section as follows: first, the displacement boundary integral equation (DBIE) which forms the basis of the BEM is defined along with commonly used BEM notation. Next, the NURBS approximation of both the displacement and traction fields is given, the collocation positions are defined, and the methods used to evaluate all boundary integrals are described. Finally, an overview of the solution procedure is described to allow for numerical implementation.
Fig. 11. Definition of problem domain and boundary with source and field points.
Z C
t i ui dC þ
Z X
bi ui dX ¼
Z C
t i ui dC þ
Z X
bi ui dX:
ð25Þ
For simplicity in this article we will neglect body forces, i.e. bi = 0 to remove one of the volume integrals.1 The other may be removed by an appropriate selection of the (still arbitrary) stress state rij ; eij . It is most common in the BEM development for elasticity to take this state to result from an infinite point force injected through the body force term, i.e.
bi ¼ Dðx X0 Þei ;
ð26Þ 0
4.1. Boundary integral equation To allow understanding of later sections, a few key concepts of BEM are given presently where more comprehensive derivations of the method are given by several other authors e.g. [5,15]. We begin by defining the domain of the boundary value problem as X with a boundary C = oX. In addition, we define two points x0 2 C and x 2 C known as the ‘source point’ and ‘field point’ respectively, and the distance between them r: = kx x0 k. These are all depicted in Fig. 11. One of the most common ways of deriving the BEM is to utilise Betti’s reciprocal theorem which states that, given two states rij, eij and rij ; eij , then these can be related (assuming linear elasticity) as:
Z X
ij ij d
r e X¼
Z X
ij ij d
r e X;
where D is the Dirac delta function, X is a point within the domain (X0 2 X, X0 R C) and ei are the unit vector components. The solution for such a problem (of a point force within an infinite domain) is commonly referred to as Kelvin’s solution and solves for the components ui and ti as:
ui ¼ U ij ej
ð27Þ
where Uij and Tij are called ‘fundamental solutions’. Appendix B gives explicit expressions for the fundamental solutions for 2D elasticity problems. If these expressions are now substituted into (25) along with the property of the Dirac delta function given by R 0 0 X Dðx X Þuj dX ¼ uj ðX Þ, the following boundary integral equation can be written:
uj ðX0 Þ þ
Z
T ij ðX0 ; xÞuj ðxÞ dCðxÞ ¼
C
ð24Þ
Applying the divergence theorem and using equilibrium equations of (1), the displacement–strain relation of (2), and the definition of tractions given by ti = rijnj, Eq. (24) can be expressed as:
t i ¼ T ij ej ;
Z
U ij ðX0 ; xÞt j ðxÞ dCðxÞ;
ð28Þ
C
where the dependence on the source and field points has been made explicit. This equation consists almost entirely of quantities specified on the boundary except that the point X0 is restricted to the interior of the domain. Considering the limit as this point is taken toward the boundary, a jump term Cij appears and the first integral in (28) must be evaluated in Cauchy Principal Value (CPV) limiting sense. This is written as:
Z C ij ðx0 Þuj ðx0 Þ þ -- T ij ðx0 ; xÞuj ðxÞ dCðxÞ C Z 0 ¼ U ij ðx ; xÞt j ðxÞ dCðxÞ; i; j ¼ 1; 2;
ð29Þ
C
R where the integral sign -- denotes that the integrand is a CPV integral and all terms are now specified on the boundary. This equation defines the displacement boundary integral equation. To allow for numerical implementation, we discretise the boundary into a non-overlapping set of Ne elements giving:
C¼
Ne [
Ce ;
Ci \ Cj ¼ 0;
i–j:
ð30Þ
e¼1
Fig. 10. NURBS pffiffiffi basis functions corresponding to knot vector N = {0, 0, 0, 1, 1, 1} and w ¼ f1; 1= 2; 1g for quarter circle geometry.
1 This restriction is made for simplicity in this article, but there are various formulations to allow for body forces in BEM implementations. Body forces may easily be included in the isogeometric BEM framework.
92
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
This allows (29) to be rewritten:
C ij ðx0 Þuj ðx0 Þ þ
Ne X m Z X e¼1 l¼1
¼
Ne X m Z X e¼1 l¼1
þ1
1
þ1 1
T ij x0 ; x ^n Nl ^n J e ^n d^n uelj
U ij x0 ; x ^n Nl ^n J e ^n d^n t elj ;
ð31Þ
where l is a local node number, on element e, that varies from 1 to m = 2, 3,. . . for linear, quadratic elements etc., Nl is the local shape function for node l; Je ¼ dC=d^ n is the Jacobian of transformation, ^ n 2 ½1; 1 is a local coordinate system and uelj and telj are displacements and tractions, respectively, at local node l on element e. Fig. 12 illustrates a discretised domain and an example quadratic boundary element. To arrive at a system of equations, the conventional procedure is to collocate at a series of points around the boundary. This corresponds to placing the point x0 at a discrete set of points; most often, for classical (i.e. unenriched) BEM formulations the nodal positions are chosen since this automatically provides a square set of equations. In this way, by grouping the first two terms on the left-hand side of (31) the system of equations can be written in matrix notation as:
½Hfug ¼ ½Gftg;
Fig. 13. Example IGABEM domain created using NURBS with knot vector N = {0, 0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 9, 9}.
ð32Þ
where H is a square matrix containing a combination of the integrals of the Tij kernel and the jump terms, u is vector of displacement components, G is a rectangular matrix of Uij kernel integrals and t is a vector of traction components. Rearranging this set of equations by placing all unknown components on the left-hand side and known components on the right, the following can be written:
½Afxg ¼ fbg;
ð33Þ
where the vector x contains all unknown displacement and traction components. Eq. (33) may then be solved using an appropriate solver. Fig. 14. NURBS basis functions for IGABEM example.
4.2. Isogeometric approximation We now introduce the isoparametric concept to the displacement boundary integral equation where for clarity, a simple domain is used to illustrate each stage of the method. This is shown in Fig. 13 where NURBS were used to create the boundary and the first and last points are coincident to form a closed boundary. The knot vector used to create the curve is N = {0, 0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 9, 9} where the associated basis functions are illustrated in Fig. 14. To define the boundary elements, we take the unique values contained in the knot vector. Therefore, the element boundaries for this example are defined as {0, 1, 2, 3, 4, 5, 6, 7, 8, 9} in the parameter space (Fig. 14) which, when mapped onto the physical space, are depicted in Fig. 15. In
the isoparametric fashion, displacements and tractions are interpolated using the NURBS approximation (18) as:
uj ðnÞ ¼ tj ðnÞ ¼
n X a¼1 n X
Ra;p ðnÞdj ;
a
ð34Þ
Ra;p ðnÞqaj :
ð35Þ
a¼1 a
where dj and qaj represent displacement and traction coefficients respectively. Each of these is associated with a particular control point. These approximations are not used directly however, since the basis functions Ra,p(n) have a local support and should only be evaluated on those parts of the boundary where they are non-zero.
Fig. 12. BEM discretisation.
93
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
Fig. 16. Element containing collocation point in both global and local coordinate systems.
dC dn : Jð^nÞ ¼ dn d^n Fig. 15. Element definitions for IGABEM domain.
For example, in the case e = 1 with the element domain n = [0, 1], we can see from Fig. 14 that only the basis functions corresponding to a = 1, 2, 3 are non-zero. Therefore, it seems sensible to define elemental approximations of uj and tj as:
uj ðnÞ ¼
pþ1 X
Nel ðnÞdj ;
le
ð36Þ
Nel ðnÞqlej ;
ð37Þ
l¼1
t j ðnÞ ¼
pþ1 X l¼1
ð41Þ
Appendix C gives explicit expressions to evaluate this quantity. Some mention should be made of the first term in Eq. (40) relating to the jump term Cij which is distributed over the element e using the basis functions Nel . This is done because the local collocation 0 ^ coordinate n may not necessarily lie at an interpolatory point in the parameter space. This is in contrast to the conventional BEM formulation using polynomial shape functions where, by collocating at nodal coordinates, the basis functions are found to be interpolatory. To take account of this, we follow the same technique as used in [21,18] to distribute the jump term.
4.3. Collocation
where
Nel ðnÞ Ra;p ðnÞ;
ð38Þ
and the global basis function a is related to the local number l and element number e through a connectivity matrix:
a ¼ connðe; lÞ:
ð39Þ
The connectivity matrix for this example is shown in Table 1. Substituting Eqs. (36) and (37) into (29), the discretised isogeometric DBIE is given by 0
C ij ðx Þ
pþ1 X
Nel
pþ1 Z Ne X X ^n0 dle þ j e¼1 l¼1
l¼1
¼
pþ1 Z Ne X X e¼1 l¼1
þ1
U ij
1
þ1
1
le T ij x0 ; x ^n Nel ^n J ^n d^n dj
x0 ; x ^n Nel ^n J ^n d^n qlej ;
ð40Þ
where the local coordinate system ^ n ¼ ½1; 1 is used, e denotes the element containing the collocation point x0 ; ^ n0 denotes the local coordinate of the collocation point (see Fig. 16) and the Jacobian of transformation, Jð^ nÞ, is given by
Using the collocation method to arrive at a system of equations, a definition must be made to obtain the position of collocation points for isogeometric BEM. The solution is not immediately obvious, since if we follow the traditional procedure of collocating at nodal points (which correspond to control points in IGABEM), collocation will occur at points outside the boundary (see Fig. 13). In this study, we choose to collocate at the Greville abscissae [7,11] where the collocation points in parameter space are given by
n0a ¼
naþ1 þ naþ2 þ þ naþp ; p
a ¼ 1; 2; . . . ; n 1:
ð42Þ
Note that the index a varies from 1 to n 1 since we do not wish to collocate at the position corresponding to the start/end point twice. Applying this relation to the knot vector N = {0, 0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 9, 9} used in the previous section, we find the location of the collocation points in both parameter space and physical space are as shown in Figs. 17 and 18 respectively.
Table 1 Connectivity matrix a = conn (e, l) for IGABEM example. l
e
1 2 3 4 5 6 7 8 9
1
2
3
1 3 4 5 6 7 8 9 11
2 4 5 6 7 8 9 10 12
3 5 6 7 8 9 10 11 13
Fig. 17. Location of collocation points in parameter space for IGABEM example.
94
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
Fig. 18. Location of collocation points in physical space for IGABEM example.
4.4. Integration As is well-known in boundary element techniques, integration of singular integrands plays a key role in the implementation of any such method. The current work is no exception, and a few key points should be raised. Firstly, it is not possible to use rigid body motion [24,16] to evaluate the strongly singular integrals containing the Tij kernels and to see why this is the case, some further explanation is required. 4.4.1. Preclusion of rigid-body motion The system of equations given by (32) is constructed by considering the collocation point x0 to lie at each nodal point in turn, forming each row of the system of equations. In most cases we find that the collocation point does not lie in the same element we are integrating over (commonly referred to as the ‘field’ element), and the resulting Tij integrand is non-singular. But in the case when the collocation point does lie in the field element, we must use appropriate techniques to evaluate the resulting strongly singular Tij integrals. In the rigid-body technique all non-singular integrals for a particular collocation point are evaluated first, and rigid body displacements are then applied to evaluate the final remaining integrals. This is under the assumption that the number of singular integrals remaining is equal to the number of rigid body translations (e.g. x, y, z). Using polynomial shape functions and collocating at nodal points (as is the case in standard BEM), this is known to be true. This is because for a particular field element and collocation point, only the case where the collocation point and associated basis function coincide is singular – other integrands are found to be non-singular (see Fig. 19). In the case of NURBS however, we find that we cannot guarantee that the number of remaining singular integrals is equal to the number of rigid body translations. To see how this can arise, we consider a typical element that might be found in an isogeometric implementation. Fig. 20 illustrates an element with p = 3 and collocation points positioned according to Eq. (42). Considering collocation at the point n = 0 we find that only R1,3 is non-zero. But we can also show that the basis functions R2,3, R3,3 and R4,3 are of OðrÞ and therefore, as in the case of conventional BEM, only one set of singular integrals remain and rigid body motion can be applied. In the case of the second and third collocation point (n = 1/3, 2/3) we can see that all basis functions are non-zero, but more importantly, they are all of Oð1Þ. In this scenario we find that the four integrands corresponding to each of the basis functions are singular and rigid-body motion is insufficient. Explicit singular integration routines are therefore necessary. In the current implementation we choose to use a subtraction of singularity scheme as proposed by Guiggiani and Casalini [8]. This
Fig. 19. Example boundary element for case where source point lies in the field element but resulting integrand is non-singular.
Fig. 20. Example NURBS element with collocation at point where basis functions are found to be non-zero and of Oð1Þ.
creates a non-singular integrand that can be easily computed using GL quadrature and accounts for the singular part through an analytical integration. Weakly singular integrals are seen in the Uij integrands, but these do not present issues for implementation. One of the most convenient techniques to evaluate integrals of this type is the transformation as proposed by Telles [20] which removes the singularity by ensuring that the Jacobian of the transformation approaches zero at the singular point. The technique is simple to implement and efficient, and it is adopted in the implementation of IGABEM in the present study. 4.5. Solution procedure Using the previously outlined techniques for the formation of elements, collocation points and integration over boundary elements, it is possible to perform an IGABEM analysis as illustrated in Fig. 21. This is very similar to a conventional BEM analysis but with the following main differences:
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
95
Fig. 23. Mesh and basis functions for coarse mesh, p = 2.
Fig. 21. Solution procedure for IGABEM.
Fig. 22. Hole within an infinite plate - geometry, boundary conditions and material properties.
The mesh is now completely generated from CAD data. The definition of elements is given by unique knot vector values. NURBS basis functions are substituted in place of polynomial shape functions for both the geometry and unknown fields. Collocation occurs at points defined by the Greville abscissae. Explicit singular integration routines are used to calculate the Tij integrands – rigid body motion is no longer applicable.
Fig. 24. Mesh and basis functions for coarse mesh, p = 3.
5. Examples To illustrate the accuracy and flexibility of IGABEM, several 2D elastostatic examples are now presented. The first two problems of a hole within an infinite plate and an L-shaped wedge have
96
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
Fig. 25. Comparison of deformed profile with exact solution for three elements per line, p = 2.
Fig. 26. L2 relative error norm plot for plate with hole.
Fig. 28. IGABEM mesh and basis functions for L-plate problem for one element per line and p = 2.
Fig. 27. Geometry, boundary conditions and material properties for L-plate problem.
closed-form solutions and therefore a convergence study carried out on each. In addition, a comparison is made with an isoparametric quadratic BEM to assess the accuracy of IGABEM. Finally, the problem of an open spanner is illustrated to show the ability of the method to handle arbitrary CAD geometries. 5.1. Hole within an infinite plate The problem of a hole within an infinite plate subject to a traction Tx at infinity has an exact solution given by e.g. [3]. The problem is defined in Fig. 22 where symmetry is applied to the bottom and right edges and exact tractions are applied along the top and
left edges. Using the polar coordinate system (q, h) defined in Fig. 22, these tractions can be found using:
! ! Tx R2 Tx R2 R4 rqq ðq; hÞ ¼ 1 2 þ 1 4 2 þ 3 4 cos 2h; 2 q 2 q q ! ! T R2 T R4 rhh ðq; hÞ ¼ x 1 þ 2 x 1 þ 3 4 cos 2h; 2 q 2 q ! 2 4 T R R rqh ðq; hÞ ¼ x 1 þ 2 2 3 4 sin 2h; 2 q q
ð43Þ ð44Þ ð45Þ
where an appropriate transformation to cartesian coordinates and the relation:
ti ¼ rij nj
ð46Þ
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
97
is required. Choosing the order of the NURBS basis function to be quadratic (p = 2) in the first instance, the minimum number of control points required to represent the geometry of this problem exactly is n = 11, as shown in Fig. 23(a), where the first and last points are coincident. The knot vector in this case is found to be:
N ¼ f0; 0; 0; 1=9; 1=9; 3=9; 3=9; 6=9; 6=9; 7=9; 7=9; 1; 1; 1g;
ð47Þ
with the weights:
pffiffiffi wa ¼ f1; 1; 1; 1= 2; 1; 1; 1; 1; 1; 1; 1g a ¼ 1; 2; ; n:
ð48Þ
Using the definition of elements and collocation points as described in Sections 4.2 and 4.3, the IGABEM mesh for this set of control points and knot vector is shown in Fig. 23(b). The NURBS basis functions for p = 2 in parameter space are shown in Fig. 23(c). In
Fig. 30. Element and collocation point definitions for three elements per line and deformed profile.
Fig. 31. L2 relative error norm plot for L-plate problem.
Fig. 29. IGABEM mesh and basis functions for L-plate problem for one element per line and p = 3.
addition, to evaluate the effect of higher order NURBS basis functions, the geometry and unknown fields were described by setting p = 3 with the mesh and basis functions illustrated in Fig. 24. Fig. 29 illustrates the corresponding control points, elements, collocation points and resulting basis functions for the plate with a hole problem. The knot insertion algorithm (equivalent to h-refinement) as described in [9,14] was used to carry out a convergence study with
98
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
uniform refinement applied around the boundary for both p = 2 and p = 3. In addition, the problem was analysed using quadratic isoparametric boundary elements using an equivalent mesh refinement strategy. Exactly the same number of Gauss points were used to evaluate each of the boundary integrals given by the second and third terms in Eq. (31) for both the IGABEM and BEM analyses. Fig. 25 illustrates an IGABEM mesh with three elements per line and the deformed IGABEM profile. Excellent agreement with the analytical solution is seen. Using the following definition to calculate the relative L2 error norm in displacements around the boundary:
eL2 ¼
ku uex kL2 ; kuex kL2
ð49Þ
where
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uZ dp u X kukL2 ¼ t ðui Þ2 dC;
ð50Þ
C i¼1
a comparison can be made between IGABEM and BEM (Fig. 26). In the case of IGABEM with p = 2 and quadratic BEM, both methods converge at the same rate but importantly, IGABEM demonstrates a consistently lower error for all meshes. For IGABEM with p = 3,
we see, as expected, that a higher convergence rate is obtained with lower errors than the equivalent second order mesh. 5.2. L-shaped wedge The next problem which was considered was the L-shaped wedge which exhibits a singularity at the wedge apex. The analytical solution to this problem is given by Szabó and Babuška [19] where a wedge angle of 2a = 3p/2 was used. Considering only the mode 1 loading case, exact tractions were applied along all faces with appropriate displacement constraints as shown in Fig. 27. Material properties E = 1e5 and m = 0.3 were used under plane strain conditions. The problem was solved using four different methods: quadratic BEM with uniform h-refinement, IGABEM with p = 2 and uniform h-refinement, IGABEM with p = 3 and uniform h-refinement and finally IGABEM with p = 2 and graded h-refinement towards the wedge apex. For the case of one element per line and p = 2, the control points are shown in Fig. 28(a) with collocation points and elements shown in Fig. 28(b). The knot vector for this example is given by
N ¼ f0; 0; 0; 1=6; 1=6; 2=6; 2=6; 3=6; 3=6; 4=6; 4=6; 5=6; 5=6; 1; 1; 1g;
Fig. 32. Open spanner problem.
ð51Þ
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
99
Fig. 33. Von Mises stress for open spanner problem.
with wa = 1 " a and p = 2. The resulting NURBS basis functions (which in this case are B-splines) are shown in Fig. 28(c). Likewise, the control points, elements, collocation points and basis functions for p = 3 are shown in Fig. 29. Fig. 30(a) shows a mesh of three elements per line with the deformed IGABEM profile shown against the exact solution. As can be seen, good agreement can be seen, but the effect of the singularity at the wedge apex is clearly evident (cf. Fig. 25). Plotting the same deformed profile but using the graded mesh (with p = 2), we can visually verify that the apex singularity is captured more accurately. Most importantly though, Fig. 31 shows that comparing IGABEM against BEM, consistently higher accuracies are achieved for an equivalent number of degrees of freedom for both p = 2 and p = 3, with the graded mesh achieving the most accurate results of all the analyses. 5.3. Open spanner Finally, to illustrate the ability of IGABEM to handle arbitrary geometries which may be taken directly from CAD, the problem of an open spanner as illustrated in Fig. 32(a) was analysed using IGABEM. A uniform traction of ty = 10 is applied at the right edge with all other faces traction-free. The displacement boundary conditions are also shown in Fig. 32(a). The control points for this example are shown in Fig. 32(b) while Fig. 32(c) illustrates the collocation points and element edges. To allow for comparison, the problem was also solved using FEM. Fig. 32(c) shows the deformed mesh profile obtained using both FEM and IGABEM where good agreement is obtained. Von Mises stresses obtained from both methods are shown in Fig. 33 where good agreement is also seen. This example illustrates an important concept: through IGABEM, it is possible to integrate CAD and BEM into a seamless
design/analysis. The method only requires knot vectors and control points to completely define the mesh and these are taken directly from CAD. Therefore, it is possible to analyse an extremely wide variety of geometries (only restricted by the capabilities of NURBS) using the present method.
6. Conclusions A formulation for isogeometric boundary element analysis has been presented. The problem geometry and the approximation spaces for the displacement and traction solutions are all provided, using the isoparametric concept, from a set of NURBS basis functions. This offers a number of key improvements that together offer a significant advantage over conventional piecewise polynomial formulations. Firstly, the NURBS basis permits exact representation of commonly used geometries including circular arcs. Secondly, a meshing process is no longer required because the control points defining the surfaces play the role of nodes in the IGABEM analysis. Thirdly, the use of the NURBS basis has been shown to offer improved accuracy, for a given model size, in comparison with conventional piecewise quadratic BEM schemes. The combination of the isogeometric concept and the BEM is a natural and productive one, since both the NURBS boundary representation and the BEM deal with quantities entirely on the problem boundary. The authors expect this development to have a major significance in the future of computational mechanics with the deeper integration of stress analysis into CAD. Future developments include the extension of the method to three-dimensional problems and the application of the method to contact problems where it is expected that the use of NURBS will demonstrate benefits over conventional BEM implementations.
100
R.N. Simpson et al. / Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100
which allows the desired expression to be stated as:
Appendix A. Higher order NURBS derivatives Any higher order derivatives of the NURBS basis functions are given by the following general formula:
AaðkÞ ðnÞ
k W ðjÞ ðnÞRðkjÞ a;p ðnÞ j¼1 j ; WðnÞ
Pk
ðkÞ Ra;p ðnÞ ¼
where the superscript function and
(k)
ð52Þ
denotes the kth derivative of relevant
ðkÞ AðkÞ a ðnÞ ¼ wa N a;p ðnÞ;
ð53Þ
(no summation over index a) with
k j
¼
k! : j!ðk jÞ!
ð54Þ
Appendix B. 2D Fundamental solutions The fundamental solutions for Tij and Uij in plane strain are given by
U ij ðx0 ; xÞ ¼
1 1 ð3 4mÞ ln dij þ r ;i r ;j ; 8plð1 mÞ r
T ij ðx0 ; xÞ ¼
ð55aÞ
1 or ð1 2mÞ dij þ 2r ;i r;j ð1 2mÞ r;i nj r;j ni ; 4pð1 mÞr on ð55bÞ
with
r ¼ kx0 xk;
ð56Þ
lis the shear modulus, m is Poisson’s ratio, dij is the Kronecker-delta function, ni is the i-component of the outward-pointing normal at x and r,i or/oxi. For plane stress problems the fundamental solutions (55a) and (55b) may be used with an appropriately modified set of elastic constants. Appendix C. Transformation of coordinate system To transform the integrals into the local coordinate system ^ n 2 ½1; 1, the Jacobian given by (41) is required. This consists of two parts: dC/dn and dn=d^ n. The first is given by
dC ¼ dn
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 dx dy þ ; dn dn
ð57Þ
which can be determined using the derivatives of the NURBS basis functions (Section 3.4). To evaluate the second expression, we first assume that the element is defined by the interval in the parameter space as n 2 [n1, n2]. Then the following linear relationship between n and ^ n exists:
h n¼
i ðn2 n1 Þ^n þ ðn2 þ n1 Þ 2
;
ð58Þ
dn n2 n1 ¼ : 2 d^n
ð59Þ
References [1]