CAD-based collocation eigenanalysis of 2-D elastic structures

Computers and Structures 182 (2017) 55–73

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

CAD-based collocation eigenanalysis of 2-D elastic structures C.G. Provatidis National Technical University of Athens, School of Mechanical Engineering, Mechanical Design and Control Systems Division, Heroon Polytechniou 9, 15780 Zografou, Greece

a r t i c l e

i n f o

Article history: Received 16 November 2015 Accepted 17 November 2016

Keywords: Coons-Gordon Bézier B-splines CAD/CAE Global collocation Eigenvalues

a b s t r a c t This paper investigates the performance of the global collocation method for the numerical eigenfrequency extraction of 2-D elastic structures. The method is applied to CAD-based macroelements, starting from the older blending function Coons-Gordon interpolation (based on Lagrange polynomials) and extending to tensor product Bézier and B-splines. Numerical findings show equivalence between Lagrangian and Bézierian macroelements, while a mass lumping procedure is proposed for the former ones. Concerning B-splines, the influence of multiplicity of inner knots and the position of collocation points is thoroughly investigated. The theory is supported by 2-D numerical examples on rectangular and curvilinear structures of simple shape under plane stress conditions, in which the approximate solution rapidly converges towards the exact solution faster than that of the conventional finite element of similar mesh density. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction The use of CAD-based global approximation for the numerical solution of partial differential equations (CAE: computeraided-engineering) is as old as the theory of computeraided-design (CAD) itself. It is well known that an industrial team (at General Motors) under Gordon’s leadership, in early 1970s, used blending function methods, based on the ideas put forward in [1], to produce some interesting element families [2,3]. Although this team presented the mathematical background for the common description between the geometric model and the unknown variable (CAD/CAE integration), unknown reasons (perhaps the high computational cost) prevented further dissemination of this excellent idea. One decade later, E1-Zafrany and Cookson [4,5] also used Coons’ and Barnhill’s ideas for quadrilateral and triangular patches, respectively, whereas Zhaobei and Zhiqiang proposed the use of Coons’ interpolation for the analysis of plates and shells [6]. Nevertheless computational results concerning CAD-based isoparametric macroelements (occupying a Coons patch ABCD) were presented for the first time by Kanarachos, Deriziotis and Provatidis [7,8] in 2D potential and elasticity (static and dynamic) problems, where the so-called ‘‘C”-elements were successfully compared with conventional finite elements and boundary elements of similar mesh density. For a detailed review (of over 160 E-mail address: [email protected] URL: http://users.ntua.gr/cprovat http://dx.doi.org/10.1016/j.compstruc.2016.11.007 0045-7949/Ó 2016 Elsevier Ltd. All rights reserved.

references) on the use of CAD-based macroelements the reader is referred to [9]. Summarizing some of the most important previous findings concerning macroelements that occupy a 2D quadrilateral patch ABCD or a 3D hexahedral block ABCDEFGH, in chronological correspondence with the progress in CAD-theory (Coons, Gordon, Bézier, B-splines and NURBS) (see, for instance, [10]), it has been reported that: (i) (Boundary-only) Coons interpolation is capable of creating a broad family of arbitrary-noded elements that may be equivalent to that of Serendipity type. For example, the conventional 4- up to 8-noded 2D elements, as well as the 8- and 20-noded 3D elements can be directly derived applying Coons interpolation [11–14]. (ii) Gordon-Coons (transfinite blending function) interpolation when applied to a structured macroelement of which the boundary and internal nodal points lay at the same normalized (n, g) positions, degenerates to the classical Lagrangian type finite element [14]. (iii) Coons-Gordon interpolation allows for dealing with a (relatively) unstructured mesh of internal nodes. Using a single quadrilateral macroelement, not only simple shapes such as a rectangular or a circle can be treated, but also it was possible to perform analysis until the complexity of a pi-shaped domain (see, [14–16], among others). For more complex shapes, domain decomposition using large macroelements becomes necessary.

56

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

(iv) For a given number of subdivisions per direction (x and y, or n and g) along the sides of the patch ABCD, Coons-Gordon interpolation allows for the construction of a large number of alternative macroelements. The two main influencing factors are: a. The choice of blending functions and b. The univariate interpolation along the four external sides of the macroelement as well as along its so-called ‘‘interboundaries” [17]. For the completeness of the paper, some basics are summarized in Appendices A and B. (v) Bézier tensor product representation of surface geometry allows also for a global interpolation of the variable, u, within a patch. Since Bézier (Bernstein) polynomials share the same functional space with Lagrange polynomials, i.e. the classes fnn g and fgn g, a CAD-based isogeometric solution is identical with an isoparametric solution based on Lagrangian type polynomials [18]. (vi) B-splines tensor product representation of surfaces is an ‘‘analogue” to the Lagrange tensor product interpolation, in the sense that univariate Lagrange polynomials are merely substituted by B-splines in both directions. Its applicability in conjunction with Galerkin-Ritz formulation has been summarized by Höllig [19]. An extension of this theory is the dominating NURBS-based isogeometric analysis (shortly IGA) [20]. (vii) Barnhill interpolation within triangular patches is capable of producing relevant macroelements. For example, conventional three - and six-noded triangular elements can be easily derived applying Barnhill interpolation (see, [9]), whereas this interpolation has been applied in conjunction with Bézier patches [21].

Second, a Bézier tensor product is introduced for both the geometry and displacement representations. Investigation is performed to determine whether the numerical solution using Bézier formulation is identical with that which is obtained using a uniform mesh of the abovementioned Lagrange-polynomial tensor product interpolation. Third, B-splines tensor products are tested. Particular attention is paid on the handling of corners. The role of multiplicity of inner knots in conjunction with several alternative sets of collocation points is investigated. Fourth, the above results based on the global collocation are compared with the conventional bilinear FEM as well as with the Galerkin-Ritz formulation based on the same global shape functions with those used in the global collocation schemes, for the same multiplicity of inner knots. 2. One-dimensional shape functions This section refers to either of the Cartesian (x- and y-) directions, or even the curvilinear n- and g-directions (along the sides AB and DA, respectively, of a patch ABCD). Below, n corresponds to either of the subdivisions nn or ng , whereas the corresponding normalized domain is ½0; 1. Given the breakpoints x0 ; . . . ; xn , where x0  0 and xn  1, among several alternative interpolations, in this paper we focus on the following. I. Lagrange polynomials

Lj;n ðxÞ ¼

ðx  x0 Þ . . . ðx  xj1 Þðx  xjþ1 Þ . . . ðx  xn Þ ðxj  x0 Þ . . . ðxj  xj1 Þðxj  xjþ1 Þ . . . ðxj  xn Þ

ð1Þ

II. Bezier (Bernstein) polynomials As a general remark, despite the elegant formulation of the above CAD-based interpolations, even in B-splines (which are characterized by compact support) the computer effort is sometimes comparable or even higher than what the conventional FEA requires. For this reason, in 2005 it was proposed to preserve the CAD-based approximation but use a global collocation method instead of the time consuming Galerkin-Ritz formulation [22, p. 6704]. Relevant works in 2D structures concern Poisson equation [23,24], elastostatics [25], and acoustics [26]. Early studies in 1-D elastodynamics are [27,28]. In the framework of Lagrange polynomials applied to the aforementioned 1-D problems, it was shown that the lumped mass can be easily set equal to the identity matrix when the nodal points are put at the non-uniform position of Gauss points [28]. Within this context, this paper continues extending previously published ideas, from 1-D to 2-D elastodynamics. A preliminary report in 2010 by Filippatos [29] concerning circular and rectangular structures under plane stress conditions, for both Dirichlet and traction-free boundary conditions, shows that, global collocation in conjunction with Lagrange polynomials leads to acceptable results. This paper substantially extends the latter report and is structured as follows: First, uniform interpolation based on Lagrange-polynomial tensor products is used to interpolate either of displacement components (ux or uy) as well as the geometry [x(n,g) or y(n,g)]. Using Dirichlet and Neumann (traction-free) boundary conditions, several versions of global collocation (nodal collocation, orthogonal collocation as well as images of Demko’s and Greville’s abscissae) are thoroughly tested. Orthogonal collocation refers to collocation points located at several positions such as those used in Gaussian quadrature, also the roots of Chebyshev polynomials (of first and second kind). Furthermore, a previous technique for mass lumping is successfully transferred from 1-D [28] to 2-D Lagrange-based collocation problems.

Bi;n ðxÞ ¼


 n i n! xi ð1  xÞni x ð1  xÞni ¼ i!ðn  iÞ! i

ð2Þ

III. B-Splines We start with the abovementioned breakpoints,

fxb g ¼ ½x0 ; . . . ; xn ;

ð3Þ

and then we introduce the knot vector fVg:

fVg ¼ ½v0 ; . . . ; vm ;

ð4Þ

which highly depends on the chosen multiplicity k of internal knots (usually single, double, etc.), as follows:  Multiplicity k = 1:

2

3

6 7 fVgk¼1 ¼ 4x0 ; . . . ; x0 ; x1 ; x2 ; . . . ; xn1 ; xn ; . . . ; xn 5 |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} pþ1

ð5Þ

pþ1

 Multiplicity k = 2: 2

fVgk¼2

3

6 7 ¼ 4x0 ; . . . ; x0 ; x1 ; x1 ; x2 ; x2 ; . . . ; xn1 ; xn1 ; xn ; . . . ; xn 5; |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflffl{zfflffl} |fflffl{zfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} pþ1

2

2

2

ð6Þ

pþ1

and so on. It is noted that:  the maximum allowable multiplicity is k = p  1.  In all cases, C pk -continuity is ensured. Therefore, Eqs. (5) and (6) and so on, lead to the unified relationship:

m ¼ 2ðp þ 1Þ þ kðn  1Þ  1;

k ¼ 1; 2; . . .

ð7Þ

57

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

Based on the abovementioned computed knot vector fVg, the vector of control points is denoted by:

where the number of control points (nc þ 1) is related to the number ðm þ 1Þ of elements in the knot vector V, as follows:

where P n , P g and Png are projections explained elsewhere [2,3,7,8,15]. In a series of papers, it has been previously shown that Eq. (16) is the ‘‘vehicle” to produce cardinal global shape functions, where the latter operate within the entire patch [14,22]. As previously mentioned, these shape functions depend on the choice of univari-

m ¼ nc þ p þ 1

BC CD ate interpolation, i.e. the trial functions BAB j ðnÞ, Bj ðgÞ, Bj ðnÞ and

fPg ¼ ½P0 ; . . . ; Pnc ;

ð8Þ

ð9Þ

Moreover, it is widely known, for example [30,31, p. 50], that the ith B-spline function of p-degree, denoted by N i;p ðxÞ, is defined as



Ni;0 ðxÞ ¼ Ni;p ðxÞ ¼

1; if vi 6 x < viþ1 0;

otherwise

x  vi viþpþ1  x Ni;p1 ðxÞ þ Niþ1;p1 ðxÞ viþpþ1  viþ1 viþp  vi

ð10Þ ð11Þ

Also, the first derivative of N i;p ðxÞ can be effectively calculated by the recursion [31, p.59]:

N0i;p ðxÞ ¼

p p Ni;p1 ðxÞ  Niþ1;p1 ðxÞ viþp  vi viþpþ1  viþ1

ð12Þ

In both interpolations, i.e. Bézier and B-splines, considering an arbitrary 1D domain ½0; L, then, for every position x 2 ½0; L, with normalized coordinate n ¼ x=L 2 ½0; 1, we can determine the values of nc þ 1 basis functions, /i ðxÞ or /i ðnÞ, i ¼ 0; . . . ; nc (where /i stands for either of Bi;n and N i;p ), so as: (i) The Cartesian coordinate is approximated as

xðnÞ ¼

n X /i ðnÞ  xi

ð13Þ

i¼0

n X /i ðnÞ  ai

In this paper we distinguish three broad categories (here, NURBS is not dealt): 3.2.1. Lagrange polynomials Based on the abovementioned nn and ng segments, the computational mesh consists of q ¼ ð1 þ nn Þ  ð1 þ ng Þ nodal points, of which nb ¼ 2ðnn þ ng Þ belong to the boundary. Each of these nodal points (designated by ‘I’) lies at the intersection of the lines n ¼ ni and g ¼ gj , and is associated with one global shape function given by: 

I ¼ 1; . . . ; q  ð1 þ nn Þ  ð1 þ ng Þ i ¼ 1; . . . ; ð1 þ nn Þ;

j ¼ 1; . . . ; ð1 þ ng Þ

Remaining in the framework of Bézier and B-splines, it is worthy to mention that the coefficients ai in Eq. (14) are generally different than the nodal values ui associated to the breakpoints, except of the ends where a0 ¼ u0 and an ¼ un . In contrast, when working in conjunction with Lagrange polynomials, the coefficients ai in Eq. (14) refer to nodal values. Remark. In case of Bézier polynomials there are no internal breakpoints but only the so-called control points, which are produced by the knot vector:

8 9 > > < = V ¼ x0 ; . . . ; x0 ; xn ; . . . ; xn ¼ fv0 ; v1 ; . . . ; v2pþ1 g; > :|fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl}> ;

ð15Þ

with p denoting the polynomial degree (p ¼ nn or p ¼ ng ). It is well known that the 2ðp þ 1Þ elements in the vector V give ðp þ 1Þ control points [31]. In other words, in case of Bézier (Bernstein polynomial) representation, the number of coefficients is identical with those involved in a Taylor series i.e. a full polynomial of degree p.

3. Two-dimensional shape functions The two-dimensional formulas can be categorized in two broad families: 3.1. Bivariate Coons-Gordon interpolation Both the geometry, x ¼ ½x; yt , and the displacement, u ¼ ½u; vt , within the patch ABCD are approximated by:

uðn; gÞ ¼ Pn þ Pg  Png ;

where the indexes in Lagrange polynomials are n  nn and m  ng . It is noted that Eq. (16) degenerates to (17) when using Lagrange polynomials for both the univariate (along AB, BC, CD and DA) and the blending functions as well (see, [14]; also Appendix B). 3.2.2. Bézier-Bernstein polynomials Defining the knot vector along each side of a quadrilateral ABCD, the basis functions Bi;n ðnÞ and Bj;m ðgÞ (with n ¼ nn and m ¼ ng ) are determined. Then, the aforementioned onedimensional control points produce q ¼ ð1 þ nn Þ  ð1 þ ng Þ control points. Each of the aforementioned control points is associated with one global shape function given by:  /I;Bezier ðn; gÞ ¼ Bi;n ðnÞ  Bj;m ðgÞ;

pþ1

ð16Þ

;

ð17Þ

ð14Þ

i¼0

pþ1

3.2. Tensor products

/I;Lagrange ðn; gÞ ¼ Li;n ðnÞ  Lj;m ðgÞ;

(ii) The variable is approximated as

uðnÞ ¼

BDA j ðgÞ between the nodal values along the four sides (AB, BC, CD and DA, respectively) of the patch, as well as on the blending functions (see, for example, [17]). When there is no internal node, the blending functions are of first degree (i.e. E0 ðwÞ ¼ 1  w, E1 ðwÞ ¼ w, with w denoting either of n; g), whereas their number and degree increases with the number of internal nodes per direction [14,17]. Details are presented in Appendices A and B.

I ¼ 1; . . . ; q  ð1 þ nn Þ  ð1 þ ng Þ i ¼ 1; . . . ; ð1 þ nn Þ;

j ¼ 1; . . . ; ð1 þ ng Þ

:

ð18Þ

When the indexes in Eq. (18) obtain the next four values, i.e. i ¼ 1, ð1 þ nn Þ or j ¼ 1, ð1 þ ng Þ, then rather improperly we refer to ‘‘boundary” (literally ‘outer’ or ‘extreme’), otherwise to ‘‘inner” control points of the Bézier patch. 3.2.3. B-Splines interpolation Based on the basis functions N i;pn ðnÞ and N j;pm ðgÞ, the aforementioned one-dimensional control points produce q ¼ ð1 þ nc;n Þ ð1 þ nc;g Þ control points (cf. Eqs. (7) and (9)). Each of the aforementioned control points is associated with one global shape function given by: ( /I;BSplines ðn; gÞ ¼ Ni;pn ðnÞ  Nj;pm ðgÞ;

I ¼ 1; . . . ; q  ð1 þ nc;n Þ  ð1 þ nc;g Þ i ¼ 1; . . . ; ð1 þ nc;n Þ;

j ¼ 1; . . . ; ð1 þ nc;g Þ

:

ð19Þ

58

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

Similarly, when the indexes in Eq. (19) obtain the next four values, i.e. i ¼ 1, ð1 þ nc;n Þ or j ¼ 1, ð1 þ nc;g Þ, then we schematically refer to ‘‘boundary” (outer or extreme), otherwise to ‘‘inner” control points of the B-splines patch. Then, for every position ðx; yÞ within the patch ABCD, we can determine the values of q basis functions, /i ðx; yÞ or /i ðn; gÞ; i ¼ 1; . . . ; q, so as: (i) The Cartesian coordinate is approximated as

9 q X > > /i ðn; gÞ  xi > > =

xðn; gÞ ¼

i¼0

ð20Þ

q > X > > yðn; gÞ ¼ /i ðn; gÞ  yi > ; i¼0

Let

us

consider

the

field

of

displacements

uðx; yÞ ¼

T f uðx; yÞ vðx; yÞ g on the domain X, which is required to fulfil the governing equation (26). The boundary conditions that correspond to this problem are of two types as follows:

 ðtÞ on C1 (Dirichlet type). (a) Essential conditions, such as u ¼ u (b) Natural conditions of the type t ¼ tðtÞ on C2 (Neumann type), with t denoting the traction, the components of which are related with the stress tensor as follows:

 x þ sxy n y t x ¼ rx n  x þ ry n y ; t y ¼ syx n

ð28Þ

 x and n  y are the components of the outward normal vector where n on the boundary C2 . The total boundary is C ¼ C1 þ C2 .

(ii) The variable is approximated as

uðn; gÞ ¼

4.2. Boundary conditions

q X /i ðn; gÞ  ai

ð21Þ

4.3. Boundary and domain discretisation

i¼0

As previously occurred in one-dimensional approximation, the coefficients ai in Eq. (21) are generally different than the nodal values ðui ; vi Þt associated to the breakpoints, except those related to the boundary. 4. Solution of two-dimensional problem 4.1. Governing equations In the case of an elastic, isotropic, homogeneous 2D structure, the governing equation in the domain X is given by

€ in X; LT r þ b ¼ qu where r ¼ f rx

ð22Þ T



ry sxy g is the stress tensor, b ¼ bx by



@=@x body force, L ¼ 0 T

0 @=@y

T

is the @=@y is the stress operator, q is the @=@x

€ ¼ fu € v € gT is the acceleration vector. mass density and u Considering the Hookean elasticity matrix in 2-D problems (plane stress or plane strain)

0

d11 B E ¼ @ d21 0

d12 d22 0

0

1

C 0 A; d33

ð23Þ

the relationship between the stress, r and the strain  T is

e ¼ ex ey cxy

r ¼ Ee:

ð24Þ

In addition, the relationship between the strain and the displacement vector, u ¼ f u

e ¼ Lu:

v gT is:

ð25Þ

Therefore, the final governing equation (Naviér- Kirchhoff) is written as follows:

€; Du þ b ¼ qu

ð26Þ

where

D ¼ LT EL

4.3.1. Lagrange polynomials When working in conjunction with Lagrange polynomials, the boundary is actually discretised into nb nodes or, equivalently into 2nb degrees of freedom (DOF) of which n1 DOF belong to C1 and n2 DOF to C2 , so as n1 þ n2 ¼ 2nb . While the number of unknown displacement components along C2 is always n2 , the number of prescribed tractions may be slightly different, particularly when a certain number of c corners appears, because two independent normal tractions per node are possible. Therefore, the total number ^ 2 ¼ n2 þ c. Moreover, of traction components increases from n2 to n ni internal nodes (2ni DOF) exist within the domain X of the patch ABCD. 4.3.2. Bézier and B-splines In contrast, when working in conjunction with either Bézier or B-splines approximations, the unknown generalized variables, a, do not necessarily correspond to boundary displacements, but they are associated to q control points (the latter may also lay within or even outside the domain X; sometimes, mostly along straight line segments, these may also lie along the true boundary C). Within this context, the q control points (cf. Eq. (18) or Eq. (19)) are divided in two categories, i.e. the nc;in ‘‘internal” and the nc;b ‘‘boundary” control points (for the definition of the terms ‘‘internal” and ‘‘boundary” we refer to Section 3.2), so as q ¼ nc;in þ nc;b . In more detail, if a side of the quadrilateral patch ABCD (e.g. AB) is straight, the corresponding control points lie on this side (AB). In contrast, if the side is curved, then only the extreme control points (P0 and Pn ) will belong to the boundary and even they coincide with the corners (e.g., A, B), whereas the rest will be either inside or outside the domain in accordance to the curvature of the curve AB. Unifying the two approaches (Lagrange polynomials on one hand versus Bézier or B-splines on the other), nb will always represent the number of boundary nodes (or control points, accordingly), whereas ni the number of internal nodes (or control points, accordingly). Also, q will represent the total number of nodes (or control points, accordingly), whereas the total number of DOF in the macroelement will be equal to 2q. 4.4. Global collocation

ð27Þ

is a quadratic operator. In this way, stress equilibrium equations (Newton’s second law) take the standard form of a partial differential equation.

4.4.1. Using Lagrange polynomials When approximating the displacement components through Lagrange polynomials (Eq. (17)), the procedure is performed in two steps as follows.

59

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

STEP-1: Fulfilment of the governing equation In the absence of body force (f = 0), the equation of motion (26) becomes:

2 4

2

2

2

@ @ ðd11 @x 2 þ d33 @y2 Þ

ðd12 þ d33 Þ

@ ðd12 þ d33 Þ @x@y

2

2

@ @x@y

@ ðd11 @y 2

þ

2

@ d33 @x 2Þ

3 5



uðx; yÞ vðx; yÞ



 ¼q

€ ðx; yÞ u : €ðx; yÞ v ð29Þ

By virtue of either of Eqs. (17)–(19), Eq. (29) becomes:

q

q X

€j ¼ /j u

j¼1

q

q X

€j ¼ /j v

j¼1

# @ 2 /j @ 2 /j þ d uj 33 @x2 @y2 j¼1 " # q X @ 2 /j vj þ ðd12 þ d33 Þ @x@y j¼1 q X

q X

"

"

2

ðd12 þ d33 Þ

q X j¼1

Muu

0

0

Mvv



"

#

@ /j uj @x@y

# @ 2 /j @ 2 /j d11 2 þ d33 2 vj @y @x

ð31Þ

€ ðtÞ u v€ ðtÞ



 þ

Kuu

Kuv

Kvu

Kvv



uðtÞ

vðtÞ

¼



0 ; 0

ð32Þ

½Muu ij ¼ ½M vv ij ¼ q/j ðxi Þ;

ð33Þ

while those of the stiffness matrix by:

@ 2 /j ðxi Þ @ 2 /j ðxi Þ þ d ; 33 @x2 @y2

½K uv ij ¼ ðd12 þ d33 Þ

@ 2 /j ðxi Þ @x@y

½K vu ij ¼ ðd12 þ d33 Þ

@ 2 /j ðxi Þ ; @x@y

@ 2 /j ðxi Þ @ 2 /j ðxi Þ þ d ½K vv ij ¼ d11 33 @y2 @x2

ð34aÞ

q
X

 @/j @/  y d33 j uj þn @x @y j¼1
 q X @/ @/  x d12 j þ n  y d33 j vj ; n þ @y @x j¼1

ð34bÞ

 x d11 n

q
X

 @/j @/  x d33 j uj þn @x @y j¼1
 q X @/ @/  x d33 j þ n  y d11 j vj þ n @x @y j¼1

A22

8 9  > < U1 > = A2i  U2 ¼ ft2 g > : > ; Ui

ð35Þ

ð37aÞ

ð38Þ

Since the value of Dirichlet b.c.’s does not play any role in the eigen 1 ¼ 0 and also t2 ¼ 0. value problem, we take that U In the framework of this study several alternative procedures have been tested with respect to handling the Neumann type boundary conditions, as follows. One possibility is to select so many collocation points as the number of the unrestrained nodes. This fits well when the governing equation is fulfilled at the boundary C2 as well. Then the corresponding traction-free boundary condition has to be incorporated through a Lagrange multiplier technique. This approach has been adopted by Auricchio et al. [40]. However in the case of Lagrange polynomials in conjunction with collocation points located at the Gauss points plus the Neumann type boundary, it was found that fictitious complex eigenvalues may appear (see Appendix D). Due to this fact, and since proper collocation points should be determined in order to avoid the aforementioned complex values, this approach is not adopted in this study. Another possibility is to select so many collocation points as the number of internal points, as proposed in the past [24–29]. In this case, the Neumann boundary condition is an extra equation, and can be handled by eliminating the DOF that appear along the Neumann type boundary, namely C2 . Then, solving Eq. (38) in U2 we obtain:

ð39Þ

Before we continue, concerning the smoothness of the Neumanntype boundary C2 to which Eq. (39) is applied, two cases are distinguished as follows: (1) Smooth boundary In the case of a smooth boundary C2 , a couple of tractions (tx, ty) and also a couple of displacement components (ux, uy) appear at each boundary node. As a result, the number of traction components equals to the number of displacement components thus the matrix A22 is perfectly square, and therefore Eq. (39) implies

U2 ¼ ðA22 Þ1  A2i  Ui

ð40Þ

Substituting Eq. (40) into Eq. (37a), in conjunction with  1 ¼ 0 and t2 ¼ 0, the latter leads to the final matrix form: U

€ i þ KUi ¼ 0; MU

and similarly

ty ¼

Kc2

8 9  > < U1 > = Kci  U2 ¼ f0g > : > ; Ui

A22 U2 ¼ A2i Ui

STEP-2: Fulfilment of the boundary conditions For the boundary part C2 the tractions (given by Eq. (28)) are expressed first in terms of the stress tensor (by virtue of Eq. (24)) and then in terms of the nodal displacements (using Eq. (25)) and corresponding shape functions (Eq. (17)), as follows:

tx ¼

8 9 € > > 1= :€ > ; Ui

First, we deal with the general case where the boundary consists of both parts C1 (Dirichlet type) and C2 (Neumann type). Since nodal points along C2 exist, the system of Eqs. (35) and (36) leads to the following matrix form:

½ A21

where the elements of the mass matrix are given by:

½K uu ij ¼ d11

Mc2

ð30Þ

Collocating successively Eqs. (30) and (31) at a certain number of points ðxc ; yc Þ, we derive a set of equations that can be written in matrix form as follows:



½ Mc1

d11

j¼1

þ

Splitting the columns of mass and stiffness matrices in a similar manner, in free vibrational motion Eq. (32) becomes

ð41Þ

where

 y d12 n

M ¼ Mci  Mc2 ðA22 Þ1 A2i ; ð36Þ

In the sequence, we split the displacement vector in three discrete  1 ), boundary ones sets, i.e. boundary displacements of given value (U of given traction (U2 ), as well as internal displacements (Ui ).

K ¼ Kci  Kc2 ðA22 Þ1 A2i

ð42Þ

(2) Irregular boundary. In the appearance of a corner point C along C2 , three independent stress and traction components exist at C. As a result, A22 becomes now a non-square matrix ^ 2  ðn2 þ 2ni Þ whereas U2 continues to be a vector of size ½n of size ½ðn2 þ 2ni Þ  1. Therefore, in this case the Neumann

60

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

type boundary conditions could be applied in the leastsquares sense. In order to eliminate U2 , we can, for example,

½ Mcol;1

Mcol;2

multiply both sides of Eq. (39) by its transpose AT22 thus receiving: h i 1

U2 ¼ ðAT22 A22 Þ

 AT22  A2i Ui

ð43Þ

Substituting Eq. (43) into Eq. (37a), the latter leads again to the expression of Eq. (41), where now the mass and stiffness matrices change into: 1

M ¼ Mci  Mc2 ðAT22 A22 Þ AT22 A2i ; 1

K ¼ Kci  Kc2 ðAT22 A22 Þ AT22 A2i

ð44Þ

Due to the complexity in the matrices involved in Eq. (44), below a remedy will be proposed in the ‘Remarks’ (see remark No.4). Second, let us now consider the particular case where the boundary consists of only the part C1 (Dirichlet type). Since  1 ¼ 0 whereas U2 disappears, Eq. (37a) now becomes: U

€ i þ Kci Ui ¼ f0g; Mci U

ð45Þ

where Mci and Kci are square matrices of size ½ð2ni Þ  ð2ni Þ. In conclusion, irrespectively to the type of boundary conditions (Dirichlet type only, or mixed Dirichlet and Neumann type), in the above proposed formulations only the internal DOFs participate in the final system of equations (cf. Eqs. (42), (44) and (45)). Remarks.. (1) The final matrices, M and K described by Eqs. (42), (44) and (45), are generally nonsymmetric. (2) In contrast to B-splines collocation which is characterized by banded matrices, the collocation methods using (global) Lagrange polynomials lead to fully-populated global system matrices, M and K. (3) When corners appear along the boundary C2 , the leastsquares idea can be used to condense the Neumann boundary equations from the global system (cf. Eq. (44)). The multiplication by AT22 is a global operation that can be computed fairly easily with MATLAB. However, difficulties may arise when dealing with large scale problems. (4) As a practical remedy for the abovementioned corner effect, we can form expressions for the traction vectors t0 ¼ ½t0x ; t0y 

T

and t00 ¼ ½t00x ; t00y  along both edges passing through the same corner (with outward normal vectors n0 and n00 , respectively) and then force only the corresponding normal tractions to vanish, i.e. t 0n ¼ tx n0x þ t y n0y ¼ 0 and t00n ¼ tx n00x þ t y n00y ¼ 0. In this way the shear tractions are ignored but then the matrix A22 becomes square and therefore the abovementioned least-squares scheme is by-passed. Not only this scheme is trivial and at no cost but it also was found to lead to slightly better results than the least-squares scheme. (5) Concerning the computation of global derivatives involved in Eqs.(34a,b) as well as Eqs. (35) and (36), the reader is referred to Appendix C. T

4.4.2. Using generalized coefficients In case of applying Bézier (Eq. (18)) or B-splines (Eq. (19)) approximation, we have to deal with ‘‘nodeless” generalized coefficients, a, the latter corresponding one-by-one to the control points. Following the same procedure as in Section 4.4.1, the collocation leads to the following general matrix equations system:

8 9 €1 > > :€ > ; aI

Kcol;2

8 9 > < a1 > = Kcol;I   a2 ¼ f0g > : > ; aI ð37bÞ

where a1 and a2 are vectors of coefficients that refer to the abovementioned nc;1 and nc;2 DOF related to ‘‘boundary” control points under Dirichlet and Neumann conditions, respectively, while aI refers to the associated nc;in ‘‘internal” control points (i.e. 2nc;in DOF). Again, we obtain Eq. (41) where Ui is now replaced by aI, and Eqs. (42), (44), and (45). In other words, the only significant difference is the use of different global shape functions (Eq. (18) or Eq. (19)). 4.5. The well-known Galerkin/Ritz approach Applying the Galerkin method to Eq. (26), for the free vibration problem (with u ¼ Ua) one obtains the well known matrix formulation [32]:

€ðtÞg þ ½KfaðtÞg ¼ f0g ½Mfa

ð46Þ

where ½M and ½K are the mass and stiffness matrices, respectively, which are given by:

Z

M¼q

X

UT  U dX;

Z

K¼

X

ðLUÞT  E  ðLUÞ dX

ð47Þ

Eq. (46) has been previously applied, in frequency and transient elastodynamics, in conjunction with Lagrange polynomials [8,12,16,22]. Moreover, B-splines implementation of Galerkin-Ritz [Eqs. (46) and (47)] has been previously presented in a textbook ([19] and papers therein). For the sake of completeness, we mention that: (1) Dirichlet boundary conditions are easily implemented deleting both the nc;1 rows and columns which correspond to the restrained coefficients. (2) Neumann boundary conditions treat the nc;2 ‘‘boundary” DOF in the same way with the 2nc;in unrestrained (internal) ones. For example, in the particular case of a free-free problem, no deletion of rows or columns is required. 5. Numerical implementation 5.1. Global collocation method Several forms of the global collocation method were implemented as follows. 1. Tensor product Lagrange polynomials were applied in conjunction with a global arrangement of the collocation points at (i) the internal points (nodal collocation), (ii) the Gaussian points, (iii) the roots of Chebyshev polynomials of first and second kind, (iv) images of the Demko’s abscissae, and (v) images of b n and b the Grevilles’ abscissae ( T S nþ2 : for definitions and examples the reader is referred to [34]). For a structured discretization that consists of nx  ny subdivisions [i.e. ðnx þ 1Þðny þ 1Þ nodes], the involved polynomial degrees are px ¼ nx and py ¼ ny , so the number of collocation points was taken equal to ðnx  1Þ  ðny  1Þ, globally distributed within the entire patch. After imposing the b.c.’s, this choice ensures the same number of collocation points as the number of the unknowns. By the nature of Lagrange polynomials, the aforementioned models are C p1 -continuous, whereas the meaning of multiplicity is not here applicable (each nodal point is merely associated with two corresponding DOFs).

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

2. Tensor product Bézier (Bernstein) polynomials (of polynomial degrees px  py ) were also applied in full agreement with the abovementioned Lagrange polynomials. By their nature, Bernstein polynomials do not contain any inner breakpoints, thus the meaning of multiplicity is again not applicable. The global collocation scheme consists of ðpx  1Þ  ðpy  1Þ collocation points for the entire domain. 3. Tensor product B-splines, based on nx  ny cells made of the breakpoints, with polynomial degrees px and py towards the x- and y-directions respectively, were treated in two different broad ways (using single and multiple knots) as follows. a. The case of single multiplicity (k ¼ 1) was handled using a set of ½ðpx  1Þ þ ðnx  1Þ  ½ðpy  1Þ þ ðny  1Þ collocation points. In order to ensure a global arrangement within the entire patch, these collocation points were taken at the images of Demko’s or Greville’s abscissae [30, p. 192,34]. b. The case of multiple inner knots was achieved taking the values kx ¼ px  1 and ky ¼ py  1. Then, ðpx  1Þ  ðpy  1Þ collocation points within each of the nx ny cells lead to so many linear equations as exactly the number of unknowns. Actually, this choice gives ncol ¼ nx ny ðpx  1Þ  ðpy  1Þ collocation points. Among the nc ¼ ½ðpx þ 1Þ þ kx ðnx  1Þ ½ðpy þ 1Þ þ ky ðny  1Þ control points, a number of control points nb ¼ 2½px þ kx ðnx  1Þ þ py þ ky ðny  1Þ belong to the boundary. Taking, for instance Dirichlet b. c.’s for the case of 2D elasticity (2 DOF per control point), after imposing the b.c.’s the number of remaining equations is neq ¼ 2ðnc  nb Þ, which is easily shown to coincide with twice the number of ‘‘inner” control points (2ncol ) [33]. In addition, concerning the abovementioned B-splines approximation, for single knots (k = 1), we can use more collocation points than those mentioned in the above case 3a, for example, the same as those for kx ¼ px  1 and ky ¼ py  1, i.e. ðpx  1Þ  ðpy  1Þ per cell. In such case, after imposing the boundary conditions, the matrices Mcol;I and Kcol;I become non-square. This shortcoming is easily resolved applying the least-squares technique, of which the academic implementation is to left-multiply Eq. (41) by the transpose of the matrix Mcol;I , thus leading to

 a  €I þ K M |{z} |{z} aI ¼ 0 nI nI

ð48Þ

nI nI

where

 ¼ ðMcol;I ÞT  ðMcol;I Þ; M

 ¼ ðMcol;I ÞT  ðKcol;I Þ K

ð49Þ

neq ¼ 2ni ¼ 2ðnx  1Þ  ðny  1Þ, while a Neumann problem leads to matrices of dimensions neq ¼ 4ðnx þ 1Þ  ðny þ 1Þ, a fact that is entirely different than the abovementioned collocation technique. (ii) Bézierian type elements, for which the abovementioned ðpx þ 1Þ  ðpy þ 1Þ numerical integration scheme was used. (iii) B-splines elements, for which a scheme using ðpx þ 1Þ  ðpy þ 1Þ Gaussian quadrature, for each of the nx ny cells (made of breakpoints), was applied. 6. Numerical examples In order to assess the accuracy of the proposed procedure in 2-D problems, three plane stress examples will be worked out using the proposed Global Collocation Method (GCM), which will be compared with the Galerkin-Ritz approach. The latter concerns the conventional 4-node (bilinear) elements (FEM), as well as the high-order ones (denoted as CPM) in either of Lagrange polynomials or B-splines formulation for two alternative multiplicities (k ¼ 1 or k ¼ p  1) of inner knots. All codes are of similar architecture, were developed as sole programs in MATLAB environment and were executed on a usual PC. The Bézier and B-splines basis functions (Bi;n and N i;p ) as well as their derivatives were created using the ‘‘spcol” function, which exists in the Spline Toolkit. The collocation schemes include:    

Nodal collocation. Gaussian points. Chebyshev roots (1st and 2nd kind). Images of Demko’s abscissae. b n abscissae.  Images of Greville’s- T  Images of Greville’s- b S nþ2 abscissae. For the last three cases, for detailed definitions and explanations the reader is referred to [30,p. 189,34–36]. Demko’s abscissae were determined using the ‘‘chbpnt” function. The eigenvalues were calculated using the standard ‘‘eig” function. In addition to the collocation method, for comparison purposes the Global Galerkin-Ritz Method (Section 5.2) was applied using either Lagrange polynomials or Bézier and B-splines. In the figures below, it is denoted as CPM (Coons Patch Macroelement). Since analytical expressions of eigenvalues do not exist in plane stress conditions, the ‘exact’ solution (x2exact ) was derived using a fine finite element mesh. The quality of the numerical solution ~ 2 is evaluated in terms of the absolute error, which was calcux lated as follows:

5.2. Global Galerkin-Ritz method The elements mij of the mass matrix are products of two basis functions /i times /j , where each of them is piecewise polynomial of (px -th in x and of py -th in y) degree. In the particular case of a rectangular domain, the integrant becomes a piecewise polynomial of 2px -th degree in x and thus it requires ðpx þ 1Þ-point Gauss quadrature in x-direction, and similarly ðpy þ 1Þ-point Gauss quadrature in y-direction, i.e. ðpx þ 1Þ  ðpy þ 1Þ Gauss points per integration cell. Therefore, the applied Galerkin-Ritz Method concerns: (i) Lagrangian type elements, for which a global ðpx þ 1Þ  ðpy þ 1Þ Gaussian quadrature is used for the estimation of mass and stiffness matrices. A Dirichlet problem, leads to mass and stiffness matrices of dimensions

61

er ¼

~ 2  x2exact k kx

x2exact

 100 %

ð50Þ

6.1. Example 1: Fixed square sheet of uniform thickness We consider plane stress free vibration of a square metal sheet of unit length (a = 1) and constant thickness, with a ratio of Young’s modulus over mass density equal to E=q ¼ 1000 and Poisson’s ratio  1 ¼ 0), i.e. m ¼ 0:3. The entire boundary is considered to be fixed (U under Dirichlet type boundary conditions. This example was introduced in order to avoid any ambiguities related with the way of imposing the boundary conditions. In the lack of an exact solution, a convergence analysis has led to a fine mesh of 200  200 uniform subdivisions (40,401 nodes and 40,000 square elements). The finite element solution obtained

62

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

by the latter computational mesh in ANSYS is considered as a reference (‘exact’ solution). We did not prefer a much finer mesh because the first eigenvalue was found slightly to increase. In all cases a uniform mesh of (nx  nx ) subdivisions of breakpoints has been used. In case of using Lagrange polynomials or Bézier basis functions, the aforementioned nx subdivisions per side are associated to global polynomials of nx -th degree. The performance of several sets of collocation points for the first five calculated eigenvalues is shown in Fig. 1, where the horizontal axis refers to the common number of subdivisions nx ð¼ ny Þ. It is noted that, due to double symmetry, the second eigenvalue (not shown) coincides with the first one. For reasons of additional clarity, we also present selective results for the lower ten eigenvalues in Table 1, where the proposed GCM is successfully compared with Galerkin-Ritz for the same global shape functions. It is worth-mentioning that, Lagrange and Bézier polynomials resulted in identical results provided the same collocation points were used. Remaining in the framework of global tensor product Lagrange polynomials, when they were constructed in conjunction with (non-uniformly distributed) nodal points arranged at the position of the global Gauss points in a [(p  1)  (p  1)] scheme and then nodal collocation was applied, the mass matrix degenerated to the identity matrix (as theoretically was anticipated) and also the obtained results were found to be identical to those obtained using

a uniform mesh with collocation points at the global Gauss points. Similar results had been previously found in 1D problems [24]. Furthermore, in case of B-splines, piecewise polynomials of several degrees, p, are considered. The discretization of breakpoints comprises 1, 2, 4 and 8 uniform subdivisions, using polynomial degree p = 3, 4, 6 and 7, and convergence rate of the first calculated eigenvalue is presented in Fig. 2. For a given p, the interpolation depends on the number of uniform subdivisions nn (defined by the 1 þ nx breakpoints) and the total number of unknowns is neq ¼ 2½ðp  1Þ þ kðnx  1Þ2 . In this case, the horizontal axis refers to the total number of unknowns, neq . To get a better idea, the convergence of the fifth calculated eigenvalue for several values of the polynomial degree, p (= 4, 5, 6 and 7), is shown in Fig. 3. For reasons of additional clarity, we also present selective results for the lower ten eigenvalues in Table 2, where the classical B-spline collocation is again successfully compared with the Galerkin-Ritz for the same piecewise basis functions. 6.2. Example 2: long cantilever beam We consider a long two-dimensional cantilever beam ABCD of length L = 10 and height h = 1. The beam with Young’s modulus E = 10,000, mass density q ¼ 1, and Poisson’s ratio m ¼ 0 is analyzed using the GCM and FEM models shown in Fig. 4.

Fig. 1. Example 1. Convergence of Global Collocation Method (GCM) for the first five modes, using tensor product Lagrange polynomials (second eigenvalue coincides with the first one). The curves that appear in every plot are as follows: FEM: Conventional finite elements; CPM: Coons-Patch-Macroelement (Galerkin-Ritz); GCM-Nodal: Global Collocation Method (PDE at nodes); GCM-Gaussian: Global Collocation Method (PDE at Gauss points); DEMKO: Global Collocation Method (PDE at Demko abscissae); GREVILLE-Tn: Global Collocation Method (PDE at Greville Tn abscissae); GREVILLE-Sn + 2: Global Collocation Method (PDE at Greville Sn+2 abscissae); CHEBYSHEV-1: Global Collocation Method (PDE at Chebyshev roots of 1st kind); CHEBYSHEV-2: Global Collocation Method (PDE at Chebyshev roots of 2nd kind).

63

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

Table 1 Example 1 (Fixed square sheet of uniform thickness). Differences between the errors of the calculated eigenvalues using the proposed Global Collocation Method (GCM) in conjunction with global tensor product Lagrange polynomials, versus the Galerkin-Ritz for the same shape functions. Mode

‘Exact’ eigenvalues

x2 ½s2  (200  200 elements)

1 2 3 4 5 6 7 8 9 10

138898.9923 138898.9923 197110.0286 295533.0076 377234.0101 381874.3929 381874.3929 494475.2845 556112.7789 556112.7789

Error (in%) of calculated eigenvalues Proposed Global Collocation Method (GCM) Degree of polynomial (p)

Galerkin-Ritz Degree of polynomial (p)

4

6

8

10

12

4

6

8

10

12

1.22 1.22 34.95 43.17 32.50 30.89 65.46 65.11 61.94 61.94

0.04 0.04 0.35 1.32 0.64 1.85 0.26 1.30 0.61 0.61

0.01 0.01 0.01 0.00 0.00 0.11 0.11 0.01 0.14 0.14

0.00 0.00 0.01 0.01 0.01 0.02 0.02 0.01 0.02 0.02

0.00 0.00 0.01 0.01 0.01 0.02 0.02 0.01 0.02 0.02

0.15 0.15 3.84 8.05 10.11 17.49 17.49 8.90 16.18 16.18

0.00 0.00 0.02 0.21 0.13 0.57 0.57 0.25 0.99 0.99

0.00 0.00 0.01 0.00 0.01 0.01 0.01 0.00 0.01 0.01

0.00 0.00 0.01 0.01 0.01 0.02 0.02 0.01 0.02 0.02

0.00 0.00 0.01 0.01 0.01 0.02 0.02 0.01 0.02 0.02

Fig. 2. Example 1. Convergence of B-spline Global Collocation Method (GCM) for the first eigenvalue, for several piecewise polynomials of degree p (= 3, 4, 6 and 7). Curve with additional names than those appearing in Fig. 1 are: CPM (k = 1): Coons-Patch-Macroelement (Galerkin-Ritz) with multiplicity k = 1; GCM-Bsplines (k = 1): Global Collocation Method based on B-splines (k = 1, p – 1, is the multiplicity).

As can be found in [37,38], the theoretical transverse eigenvalues of the one-dimensional model are given by:

x2i ¼


4
 ki EI  ; L qA

i ¼ 1; 2; . . .

ð51Þ

where A denotes the cross-sectional area of the beam, I its second moment of inertia, k1 = 1.87510407, etc, are known constants. Also, the longitudinal eigenvalues of the one-dimensional model are given by:

x2i ¼


 2
 ki E ;  L q

i ¼ 1; 2; . . .

where  ki ¼ ð2i  1Þp=2;

i ¼ 1; 2; 3; . . ..

ð52Þ

Generally, it is accepted that the free vibration of a cantilever beam is the combination of the transverse and the longitudinal vibrations in ascending order. According to [16], the first five eigenmodes are as follows: (i) transverse, (ii) transverse, (iii) longitudinal, (iv) transverse and (v) transverse. Despite the fact that the cantilever at hand is long (aspect ratio 10:1), Eqs. (51) and (52) are not applicable ‘‘as is”. To get more reliable values as ‘‘exact” solutions (for reference of computed values), conventional 4-node bilinear finite elements were used. Following [16], the length, L, was progressively divided into nn = 4, 6, 8, 10, 20, 40, 60, 80, 100, 120, 200, 240 and 300 uniform segments while the height, h, into ng ¼ 4 segments. Asking for higher accuracy than previously in [16], now the eigenvalues were obtained for a fine FE mesh of nx  ny = 5000  500 subdivisions (2,505,501 nodes and 2,500,000 square elements), which are considered to compose

64

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

Fig. 3. Example 1. Convergence of fifth calculated eigenvalue using B-spline Global Collocation Method (GCM), for several piecewise polynomials of degree p (= 4, 5, 6 and 7).

Table 2 Example 1 (Fixed square sheet of uniform thickness). Differences between the errors of the calculated eigenvalues using the proposed Global Collocation Method (GCM) in conjunction with tensor product B-splines within four uniform spans per direction and several piecewise polynomial degrees p, versus the Galerkin-Ritz for the same shape functions. Mode

‘Exact’ eigenvalues

x2 ½s2  (200  200 elements)

1 2 3 4 5 6 7 8 9 10

138898.9923 138898.9923 197110.0286 295533.0076 377234.0101 381874.3929 381874.3929 494475.2845 556112.7789 556112.7789

Error (in%) of calculated eigenvalues B-spline collocation Degree of piecewise polynomial (p)

Galerkin-Ritz Degree of piecewise polynomial (p)

3

4

5

6

7

3

4

5

6

7

0.02 0.02 0.60 0.35 1.57 2.89 2.89 0.67 2.21 2.21

0.04 0.04 0.56 0.23 1.16 1.86 1.86 0.41 0.91 0.91

0.02 0.02 0.29 0.20 0.61 0.59 0.59 0.29 0.10 0.10

0.02 0.02 0.05 0.09 0.07 0.00 0.00 0.07 0.13 0.13

0.01 0.01 0.00 0.01 0.01 0.02 0.02 0.01 0.06 0.06

0.01 0.01 0.08 0.31 0.31 1.07 1.07 0.38 1.27 1.27

0.00 0.00 0.00 0.03 0.01 0.38 0.38 0.05 0.54 0.54

0.00 0.00 0.01 0.00 0.01 0.01 0.01 0.00 0.04 0.04

0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01

0.00 0.00 0.01 0.01 0.01 0.02 0.02 0.01 0.02 0.02

the ‘exact’ solution. Major differences between 1-D theory and 2-D FE solution were found for the fourth and fifth mode. In the sequence, the performance of the proposed Global Collocation Method (GCM) based on Lagrange polynomials was compared with the Lagrangian type (Galerkin-Ritz) macroelement (denoted by CPM) as well as with the FE solution, in both cases using the same number of internal and boundary nodes. Therefore, the beam length was progressively divided into nx = 4, 6, 8, 10 and 20 segments while the height constantly into ny = 4 ones, as shown in Fig. 4. The convergence behavior of the lowest calculated eigen-

values is shown in Fig. 5; the third eigenvalue, which is a longitudinal mode, was omitted because it is always calculated with an excellent accuracy. It can be there noticed that the GCM performs well but it is less accurate than the Galerkin-Ritz method, for the same number of unknowns. Moreover, in order to support the use of considering only the vanishing normal tractions tn at the corners, instead of performing a least-squares scheme, a comparison between these two approaches is shown in Table 3, where one may observe minor differences. The location of points at which the Neumann b.c.’s were

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

65

Fig. 4. Example 2. Discretization of a long cantilever beam.

imposed are the (nx  1) or (ny  1) Gauss points for the horizontal or the vertical edge, respectively, plus the corner point (for each of these two edges). Concerning Bézier interpolation, the obtained eigenvalues were of the same quality with those obtained using Lagrange polynomials (for the same polynomial degree p). It was found that when imposing the Neumann b.c.’s as previously, i.e. at the level of the domain collocation points (Gaussian points, plus the corners), the Bézier based numerical solution is identical with that obtained using Lagrange polynomials (they share the same functional space). Concerning B-splines interpolation, it was found that the implementation of the least-squares when imposing traction-free b.c.’s does not always ensures acceptable results. Therefore, the abovementioned technique of vanishing the normal traction tn at the cor-

ners is imperative (A22 and A2i are given by Eq. (42)) and has been thoroughly adopted. Moreover, the so far numerical experimentation showed that it is preferable to impose the traction-free conditions at the level of control points of the corresponding boundary edge (in this case they lie on the true boundary). Typical numerical results for polynomial degrees p = 4 and p = 6 are shown in Figs. 6 and 7, respectively. 6.3. Example 3: plane arch The two previous examples refer to rectangular single macroelements, and hence there is no use of the mapping characteristics of the blending and basis functions. Also, they consider the same number of nodes along the opposite sides of the patch, the latter also having the same length. This third example differs in

66

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

Fig. 5. Example 2. Convergence of first calculated eigenvalues using tensor product Lagrange polynomials.

Table 3 Example 2 (Long cantilever beam). Differences between the errors of the calculated eigenvalues using the proposed Global Collocation Method (GCM) in conjunction with global tensor product Lagrange polynomials, and two different formulations: (i) least-squares and (ii) condensation of tractions in two normal ones. Mode

‘Exact’ eigenvalues

x2 ½s2  (200  200 elements)

1 2 3 4 5 6 7 8 9 10

1.0169 37.0818 246.7401 261.4949 878.8279 2077.1800 2220.6620 3993.0648 6168.5106 6716.1267

Error (in%) of calculated eigenvalues (i) Least-squares, three tractions per corner GLOBAL COLLOCATION METHOD (GCM) Degree of polynomial (p)

(ii) Normal tractions only at corners GLOBAL COLLOCATION METHOD (GCM) Degree of polynomial (p)

84

10  4

12  4

14  4

16  4

84

10  4

12  4

14  4

16  4

0.15 0.10 0.00 0.36 0.67 6.91 30.67 44.31 0.00 79.95

0.07 0.10 0.00 0.11 0.14 2.53 0.00 3.13 0.00 55.46

0.03 0.08 0.00 0.09 0.15 0.23 0.00 0.24 0.00 7.01

0.01 0.07 0.00 0.09 0.15 0.22 0.00 0.32 0.00 0.15

0.00 0.05 0.00 0.08 0.15 0.23 0.00 0.33 0.00 0.40

0.00 0.03 0.00 0.37 0.63 6.91 29.78 41.85 0.00 79.95

0.00 0.03 0.00 0.08 0.15 2.52 0.00 3.11 0.00 54.83

0.00 0.03 0.00 0.08 0.14 0.23 0.00 0.27 0.00 7.00

0.00 0.03 0.00 0.08 0.14 0.22 0.00 0.32 0.00 0.16

0.00 0.03 0.00 0.08 0.14 0.23 0.00 0.33 0.00 0.40

all these issues and refers to a curvilinear patch concerning the arch depicted in Fig. 8a. The geometrical dimensions and physical constants of the arch are r0 = 2.5, r1 = 5.0, E=q = 104 and m = 0.25. In the lack of an exact solution, a convergence analysis has led to the conclusion that a finite element solution with a mesh of 1,443,001 nodes in 1,440,000 quadrilateral elements (2400 and 600 uniform subdivisions in the circumferential and radial directions, respectively) can be considered as a reference (‘exact’ solution). Concerning the performance of a single Coons-Gordon macroelement, Table 4 shows a rapid convergence when GCM is

applied to a Lagrangian type macroelement. Moreover, when varying a number of subdivisions along the opposite sides AB and CD (transfinite elements: nAB –nCD ) the solution is not that far different than using the same number of subdivisons (Lagrangian type element: nAB ¼ nCD ). For purposes of a better visualization, in Fig. 8 we illustrate a medium-size transfinite macroelement for the particular case nAB ¼ 6, nBC ¼ 3, nCD ¼ 10, nDA ¼ 3. One may observe that 7  2 internal nodes (i.e. the mean average of subdivisions in the opposite edges) are uniformly arranged in the domain (Fig. 8b), and this is the rule adopted for the finer meshes that are presented on the right part of Table 4.

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

67

Fig. 6. Example 2. Convergence of first calculated eigenvalues using B-splines of piecewise fourth degree (p = 4).

Moreover, using the same location of boundary points for imposing Neumann b.c.’s, all eigenvalues calculated on the basis of Bézier formulation were found to coincide with those based on Lagrange polynomials. Finally, in Fig. 9 the B-splines collocation (GCM) is again successfully compared with the usual FEM and the global GalerkinRitz method (CPM). Since now the outer control points do not belong to the circular boundary, it is clarified that the position at which the Neumann boundary conditions were applied was determined on the basis of the control points that lie in the reference interval [0, 1] for a hypothetical straight segment. 7. Discussion The Global Collocation Method (GCM) has been developed since 2005 (see [22, p.6704] and references therein), opposed to the Galerkin formulation on which the so-called ‘‘Coons Patch Macroelement” (CPM) method was based until that time instance. Experience gained at NTUA with the CPM method (initially called ‘C’-elements, published in 1989 [7]), showed that sometimes the CPM method required considerable CPU-time, especially when natural cardinal B-splines were used, and due to this shortcoming the GCM was conceived by the author as an alternative procedure to reduce the computer effort. In the period 2008–2014, a number of further works related to GCM are [23–29]. Nevertheless, the GCM has not been fully presented in 2D elastodynamics so far. Any well-intentioned researcher with considerable knowledge of CAD and FEA history, will easily recognize that the abovementioned CPM constitutes the first ‘‘CAD-based macroelement”, in

the sense that Coons-Gordon interpolation (1964–1967, 1969– 1976) is chronologically the first mathematical formula in the CAD theory (see, for example, [3]). And since in CAD theory it is accepted that main successors of Gordon-Coons patches are the well-known Bézier-, B-splines- and NURBS-patches [10], the coming of the so-called ‘‘isogeometric method” (based on NURBS, chronologically the fifth main formula in the CAD evolution) was a matter of time. A part of the above statement is also shared by Cottrell et al. [20, p. 8] according to which, ‘‘. . .other geometric technologies that may play a role in the future of isogeometric analysis include Gordon patches (Gordon, 1969), Gregory patches (Gregory, 1983), S-patches (Loop and DeRose, 1989), and A-patches (Bajaj et al., 1995) . . .”. In other words, for almost every novel interpolation that may be proposed (e.g., in a future issue of a CAD journal) or even a formerly developed patch, a corresponding ‘‘CAD-based’ macroelement can be developed. In the rapid advent of the ‘‘isogeometric method” (IGA) [20] that was introduced in the form of Galerkin-Ritz by Hughes and associates in 2005 (based on NURBS interpolation, today adopted by all CAD-systems), most previous methods than NURBS interpolation (i.e., Gordon-Coons, Bézier, and B-splines) have come off in the margins and stopped to attract the attention of the researchers; this gap is partially covered by this paper. Before we comment on the results of this study, the state-of-the-art of relevant (global) collocation methods will be presented below. It is widely known that the collocation method is highly influenced by the choice of the collocation points, and this topic has been studied particularly in the framework of B-splines (DeBoor

68

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

Fig. 7. Example 2. Convergence of first calculated eigenvalues using B-splines of piecewise sixth degree (p = 6).

Fig. 8. Example 3. Arch-like structure (a) Definition and (b) Discretisation (nAB ¼ 6, nBC ¼ 3, nCD ¼ 10, nDA ¼ 3).

69

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

Table 4 Example 3 (Arch-like structure). Differences between the errors of the calculated eigenvalues using the proposed Global Collocation Method (GCM) in conjunction with either a single classical Lagrange macroelement or a single Transfinite macroelement. Mode ‘Exact’ eigenvalues x2 ½s2  (2400  600 elements)

Error (in %) of calculated eigenvalues

nAB ¼ nCD ¼ 8 nBC ¼ nDA ¼ 8

nAB ¼ nCD ¼ 10 nBC ¼ nDA ¼ 8

nAB ¼ nCD ¼ 12 nBC ¼ nDA ¼ 8

nAB ¼ nCD ¼ 14 nBC ¼ nDA ¼ 8

nAB ¼ nCD ¼ 16 nBC ¼ nDA ¼ 8

nAB ¼ 10; nBC ¼ 8 nCD ¼ 12; nDA ¼ 8

nAB ¼ 12; nBC ¼ 8 nCD ¼ 14; nDA ¼ 8

nAB ¼ 14; nBC ¼ 8 nCD ¼ 16; nDA ¼ 8

1 2 3 4 5 6 7 8 9 10

1.00 0.78 1.15 3.26 8.56 10.18 7.19 6.09 35.03 19.05

0.09 0.03 0.14 0.47 2.05 0.32 5.09 1.16 24.83 3.56

0.06 0.02 0.01 0.07 0.05 0.02 0.83 0.17 4.10 1.13

0.07 0.01 0.01 0.05 0.05 0.01 0.06 0.00 0.20 0.37

0.08 0.01 0.01 0.06 0.07 0.01 0.05 0.01 0.02 0.01

0.68 0.40 1.17 0.00 1.06 1.06 5.08 2.81 21.57 2.05

0.08 0.05 0.25 0.10 0.52 0.29 0.01 0.90 0.71 1.26

0.07 0.01 0.04 0.03 0.12 0.06 0.09 0.34 0.47 0.47

197.6691 441.1749 1061.6384 1271.7454 2511.4892 2668.1392 4268.2837 5171.7252 6417.8669 8340.5939

Lagrange type macroelement Subdivisions along the edges of the patch ABCD

Transfinite macroelement Subdivisions

Fig. 9. Example 3. Arch-like structure. Convergence of B-spline Global Collocation Method (GCM) for the first calculated eigenvalues, for piecewise polynomials of seventh (p = 7).

[30]). For example, Kwok et al. selected the collocation points at the B-spline maxima [39, p. 514], whereas Johnson used the Greville abscissae [35,36]. In the framework of IGA-collocation, Auricchio et al. [34] performed numerical experiments with several sites (Demko, which are the extrema of the Chebyshev splines, and Greville abscissae). Two years later, Auricchio et al. [40] tested collocation points in the interior of a patch and then addressed the more subtle case of collocation points that belong to the patch boundary. In simple words, they proposed the collocation of both the partial differential equation and the Neumann boundary condition along C2 and then they properly condensed these two equations in one (this promising idea is tested in Appendix D). Quite recently, Anitescu et al. [41] improved the accuracy and the convergence rate

of the existing IGA collocation by selecting the collocation point at the zeros of a polynomial defined on a reference interval and scaled to each knot-span; in this way better results were obtained than those using the previous Greville abscissae. As the author reported in a previous study of 2009 concerned with the performance of GCM using Gordon-Coons macroelements (see, [24]), under Neumann type boundary conditions the results depend also on the boundary points to which these are imposed. Similarly, in the framework of IGA-collocation, in 2015 De Lorenzis et al. [42] showed that the strong imposition of Neumann boundary conditions (‘pure’ collocation) may lead to a significant loss of accuracy in some situations. Therefore they proposed two remedies, i.e. a ‘‘hybrid collocation-Galerkin” approach and an ‘‘en-

70

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

hanced collocation” approach. In the same year, Anitescu et al. [41] followed the two aforementioned remedies as well. Interestingly, well before the treatment of Neumann boundary conditions was of interest in CAD-based macroelements (mostly structured), is also an issue of research in meshless methods. For example, in 2000, Zhang et al. [43, p. 337] proposed a ‘‘Hermite type collocation” by considering and extra term in the series expansion of the unknown variable u (i.e. the derivative of the basis function that is associated to the node on the Neumann boundary). Not only that, but besides the internal points, the PDE was fulfilled at the Neumann boundary as well (not the Dirichlet one). The above brief literature survey shows that the imposition of Neumann boundary conditions, in all collocation methods, is currently an open topic for further research. Let us now focus upon the results of Section 6. In general, in all tests performed for this study the Global Collocation Method (GCM) was found to perform very well. All figures show that the rate of its convergence is substantially higher than that of the conventional bilinear FEM for the same mesh density. As anticipated, GCM is slightly worse than the Galerkin-Ritz formulation (Coons-Gordon Patch Macroelement: CPM) using the same global shape functions. An interesting finding of general importance concerning GCM (based on Eq. (42)) is that in all types of boundary conditions, all those DOFs related to the boundary are eliminated, and therefore the eigenvalue problem uses so many DOFs as those associated to the internal nodes (or ‘‘internal” control points). In more details, in full consistency with the usual FEM, in Dirichlet type boundary conditions GCM demands deletion for all those columns related to the boundary DOFs (in both the mass and stiffness matrices). However, in case of partially Dirichlet and partially Neumann b.c.’s, the GCM demands deletion of the columns related o Dirichlet b.c.’s and a proper elimination of the Neumann b.c.’s. As a result, the dimensions of final matrices equal twice the number of internal nodes (or control points). Obviously, in the extreme case of the free-free problem, the GCM demands elimination of all boundary DOFs (not a mere deletion), a feature that in is contrast to the usual FEM experience. A consequence of the above finding is that, although the boundary-only (Coons) approximation was applicable for the numerical solution of the eigenvalue problem in conjunction with the Ritz-Galerkin formulation (see, for example [13]) it is not the case for the collocation method. In other words, in GCM the use of internal nodes (or variables) is absolutely necessary. Nevertheless, further research should be conducted towards the establishment of a proper set of collocation points for which the fictitious complex eigenvalues disappear (see Appendix D), at least for the case of Lagrange polynomials where it was observed. Concerning Lagrange polynomials, the global collocation method is applicable in conjunction with Lagrangian type and also transfinite elements. In the latter, the internal nodes should not necessarily be at the same (n or g) position with those along the boundary but can be uniformly distributed within the domain using a somehow different density than that determined by the boundary. In addition to the numerical results shown in Section 6, the GCM was successfully applied for the rectangle of Example 1, of which the boundary was discretised using a larger (or smaller) number of nodes than those used in the interior. In this particular case, the numerical results showed that the problem is fully controlled by the internal nodes. Following previous findings in the context of Galerkin-Ritz formulation (see, [18] and appendices therein), in this study it was found that the GCM, in conjunction with either Lagrange polynomials or Bézier (Bernstein) ones, generally leads to identical results. This was anticipated because Lagrange and Bernstein polynomials share the same functional set of monomials (power series xiyj) [18]. It is noted that a full coincidence between the two latter

functional sets using the GCM had been previously reported in 1D eigenvalue analysis [27]. Closing the discussion, it is useful to examine the technical differences between the abovementioned Lagrangian, Bézierian and B-splines macroelements, when imposing the Neumann type boundary conditions (ti ¼ rij nj ¼ 0; i; j ¼ x; y). As an illustrative example, let us consider a rectangular beam ABCD (with nx = 6 and ny = 5 subdivisions) and focus upon the point P along the bottom side AB, as shown in Fig. 10. In this particular case, both Lagrange and Bézier approximations will be of polynomial degrees px = 6 and py = 5 (towards x- and y-direction, respectively). For both sets of polynomials the strain ex ¼ @u=@x at P depends on all polynomials (Lagrange and Bézier) associated to the nodes 1 up to 7 along the side AB. In contrast, the strain ey ¼ @v=@y at P in influenced by all Lagrange polynomials along the vertical line passing through P (nodes 4, 11, 18, 25, 32 and 39), whereas it is influenced by only two Bernstein (Bézier) polynomials associated to the nearest control points (nodes 4 and 11). Finally, concerning the shear strain cxy ¼ @u=@y þ @v=@x, in the Bézier formulation the first term of the sum is influenced by the two nearest nodes (4, 11) whereas the second term is influenced by all nodes (1 up to 7) along AB. Summarizing, upon the imposition of a Neumann type boundary condition at the point P, Lagrange polynomials influence all nodes along the horizontal (outer-tangential, t) and vertical (inward-normal, n) lines passing through P, whereas Bernstein (Bézier) polynomials influence again all nodes along the tangential direction and only a surplus of two nodes in the inward normal direction (i.e. the nodes 17, 18, 19, 24, 25, 26 inside the shadowed area EFGH in Fig. 10 are not affected). Using the aforementioned remarks in Bézier (as well as in B-spline formulation, see below), the matrix inversion [(A22)1] may not be blindly performed as a black-box but in a more efficient way. While Bernstein (Bézier) polynomials seem to reduce the overall computational effort, the situation may be reversed in favour of the cardinal Lagrange polynomials by shifting the internal points at the position of the collocation points (say Gauss points) and then performing nodal collocation, thus obtaining a lumped mass (identity matrix: M = I) [28]. In the latter case, only the eigenvalues of the stiffness matrix K have to be found. Concerning B-splines where the polynomial degrees px and py are generally smaller than the number of spans determined by the breakpoints (nx and ny, respectively), all three strain components (ex , ey , cxy ) are influenced by only a few control points in the neighbourhood of P in both the tangential and the normal (inward) direction. From the computational point of view, the fact

Fig. 10. Formation of tractions at point P.

71

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

of the compact support and the capability to deal with any large number of control points along one direction makes the use of Bsplines more attractive than Lagrange and Bézier polynomials. Nevertheless, transfinite macroelements still keep their value until an efficient B-spline unconstructed procedure is proposed. Finally, it is noted that this study did not include any domain decomposition aspects (for 2D potential problems we refer to [24]), as relevant studies in elasticity have been recently reported by others [40]. 8. Conclusions The proposed global collocation method was based on tensor products of several basis functions such as Lagrange polynomials, Bézier polynomials and B-splines. The method uses the same basis functions that had been previously used in the Galerkin-Ritz formulation. Not only structured tensor products but also slightly unstructured (transfinite: Coons-Gordon) meshes can be used, at least in conjunction with Lagrange polynomials. In contrast to the Galerkin-Ritz formulation, in the proposed global collocation method the estimation of mass and stiffness matrices does not need any domain integral to be computed thus the computer effort is highly reduced. The proposed global collocation method is applicable using either single or multiple internal knots. In general, multiple knots imply low continuities and are applicable in conjunction with a local arrangement of the collocation points; therefore a large number of degrees of freedom are required. In contrast, single knots can be utilised using either images of Demko’s or Greville’s abscissae or using more collocation points (for example, such as those used in multiple knots) in conjunction with the application of a least-squares scheme. Obviously, the extension of the proposed approach to other types of partial differential operators as well as to three-dimensional problems is straightforward. Appendix A. Coons-Gordon interpolation Let us consider the unit square ABCD in conjunction with nn þ 1 n-lines vertical to the n-axis at the points: ½n ¼ ½0  n0 ; n1 ; . . . ; nnn  1, and ng þ 1 g-lines vertical to the g-axis at the points: ½g ¼ ½0  g0 ; g1 ; . . . ; gnn  1. In this way, the internal points are arranged in ðnn  1Þ columns and ðng  1Þ rows within the macroelement ABCD. The blending functions Ei ðnÞ for i ¼ 0; 1;    ; nn are chosen so that they define a cardinal set of basis functions, i.e. Ei ðnk Þ ¼ dik (dik is Kronecker’s delta), where nk refers to the elements of the aforementioned vector ½n. The blending functions Ej ðgÞ for j ¼ 0; 1;    ; ng are chosen in an analogous manner. Among other possibilities (e.g. cardinal natural B-splines [17]), the blending functions can be chosen to be Lagrangian polynomials. Based on these blending functions, the three projections are given by:

Pn ðwÞ ¼

nn X wðni ; gÞ  Ei ðnÞ i¼0

Pg ðwÞ ¼

ng X wðn; gj Þ  Ej ðgÞ j¼0

Png ðwÞ ¼

ng nn X X wðni ; gj Þ  Ei ðnÞ  Ej ðgÞ

sides AB, BC, CD and DA. Considering that each of the aforementioned sides include q1 , q2 , q3 and q4 nodes, respectively, and further using a local numbering for the trial functions, then the global shape functions /j ðn; gÞ are obtained. Three characteristic cases are as follows. For the corner node A: AB /A ðn; gÞ ¼ E0 ðnÞBDA q4 ðgÞ þ E0 ðgÞB1 ðnÞ  E0 ðnÞE0 ðgÞ

For a node intermediate to the side AB:

/i ðn; gÞ ¼ E0 ðgÞBAB i ðnÞ;

i ¼ 2; . . . ; ðq1  1Þ

For an internal node:

/j ðn; gÞ ¼ Ek ðnÞEl ðgÞ;

k ¼ 2; . . . ; ðnn  1Þ;

l ¼ 2; . . . ; ðng  1Þ

Appendix B. Basic properties and remarks on Coons-Gordon elements For the sake of briefness, the following definitions are made: Definition 1. Let us call as Group-1 all those one-dimensional cardinal basis functions /j , which by definition fulfil the condition /j ðxi Þ ¼ dij ([0, 1]-type). Typical cases are piecewise-linear, piecewise quadratic etc., cardinal natural cubic B-splines and Lagrange polynomials, among others.

Definition 2. Let us also call as Group-2 the many existing ‘‘nodeless” 1D interpolations such as Bézier (Bernstein) polynomials, full B-splines, NURBS, Legendre, Laguerre or Chebyshev polynomials, and trigonometrical functions (1D Fourier series). Obviously, in the current formulation only the Group-1 can be used, because each member of this class fulfils the requirements of blending functions (Ei ðnk Þ ¼ dik ). Moreover, the trial functions, Bj , may belong to either of Group1 or Group-2. If the trial functions Bj belong to Group-1, the formulation is easier, in the sense that the involved mass and stiffness matrices refer directly to nodal displacement values, as the traditional finite element method does. In this case, the four blending functions (piecewise linear, cardinal natural cubic B-splines, and Lagrange polynomials) could be combined with three alternative sets of trial functions (piecewise linear, cardinal natural cubic B-splines, and Lagrange polynomials, among others) thus leading to 4  3 = 12 alternative formulations. A relevant study is [17]. If however the trial functions belong to Group-2, then nodeless variables, a, are introduced and the computation is performed in terms of a. For m  n uniformly placed nodes, below two lemmas are reminded.

Lemma #1. When the blending functions are piecewise linear polynomials, then the resulting global shape functions coincide with those called ‘‘hat”-functions, which are well-known in FEM literature (bilinear elements). The proof is trivial.

i¼0 j¼0

where wðn; gÞ represents either of the vectors xðn; gÞ or uðn; gÞ. In the sequence, the boundary values of the displacement vecBC tor are interpolated through the trial functions BAB j ðnÞ, Bj ðgÞ, DA BCD j ðnÞ and Bj ðgÞ that are applied along the corresponding four

Lemma #2. When the blending functions are Lagrange polynomials, then the Coons-Gordon macroelement coincides with the wellknown Lagrangian type finite element. The proof is achieved by inspection considering that the blending functions Ej coincide with the trial functions, Bj .

72

C.G. Provatidis / Computers and Structures 182 (2017) 55–73

Table 5 (1D bar fixed at one end). Calculated eigenvalues extracted by applying elimination, or without elimination, of the degree of freedom at the free end, when Lagrange polynomials are used. Mode

Exact eigenvalues x2 ½s2 

1 2 3 4 5 6 7 8 9 10

2.46740110 22.20660990 61.68502751 120.90265391 199.85948912 298.55553313 416.99078595 555.16524756 713.07891798 890.73179720

Applying elimination of the free DOF [Eq. (D2)]

Without elimination of the free DOF [Eqs. (D4) & (D5)]

Calculated eigenvalue

Error (%)

Calculated eigenvalue

Error (%)

2.46740110 22.20660990 61.68503257 120.90382191 199.80516704 297.58782188 452.95004744 908.59061476 3420.83645043 –

0.00 0.00 0.00 0.00 0.03 0.32 8.62 63.66 379.73 –

2.46740109 22.20446026 61.88334349 139.28  27.40i 139.28 + 27.40i 193.82535792 95.74  381.27i 95.74 + 381.27i 444.98811642 3418.56102041

0.00 0.01 0.32 – – – – – – –

Appendix C. Global derivatives The global derivatives of the shape functions, which are requires for the computation of mass and stiffness matrices as well as for imposing given tractions, are given by the following relationships:


2 2
2 2 2 @ /j @n @ /j @n @ g @ /j @ g @ /j @ 2 n @/j ¼ þ2 þ þ 2 2 @x @x @x @n@ g @x @x @ g2 @x2 @n @n

u1 = 0), for example that of the free node, un. As a result, the produced mass and stiffness matrices become of size (n  2)  (n  2) and therefore up to the (n  2)-th eigenvalue can be calculated. The stiffness matrix is given by

kij ¼ L00j ðxi Þ þ

2

þ

@ 2 g @/j @x2 @ g

ðC:1Þ

@ 2 g @/j @y2 @ g

ðC:2Þ

and


 2
 2 @ 2 /j @n @n @ /j @n @ g @ g @n @ /j þ ¼ þ @x @y @n2 @x @y @x @y @n@ g @x@y
 2 @ g @ g @ /j @ 2 n @/j @ 2 g @/j þ þ þ 2 @x@y @n @x@y @ g @x @y @ g The terms involved in the above equations, i.e. as well as

@g @2 g @g @2 g , , , @x @x2 @y @x@y

and

@2 g , @y2

L0n ðLÞ

L00n ðxi Þ;

ODE

ðC:3Þ @n @ 2 n @n @ 2 n , , , @x @x2 @y @x@y

and

@2 n , @y2

are calculated as usual, starting from

Appendix D. Handling the Neumann b.c. in 1-d collocation problems

2




p2 E ; j ¼ 1; 2; . . . 2 q 4L

ðD1Þ

Let us also assume that the abovementioned bar is uniformly divided into (n  1) segments with n nodal points, and relevant Lagrange polynomials Lj(x), of degree p = (n  1), are used to interpolate the axial displacement u(x). Concerning the collocation points, here we present two particular procedures, which are both based on the sites determined by (n  2) Gauss points in the interior of the bar. The second procedure requires an additional collocation point at the free end. The first procedure was presented by Provatidis [27] and consists of eliminating one of the n nodal DOFs (excepting the fixed one,

ðD3Þ

Setting unit magnitudes (E = 1, A = 1, L = 1) to avoid any necessary scaling factors for ensuring the same units, for every internal point:

kij ¼ L00j ðxi Þ;

i ¼ 2; . . . ; n  1;

j ¼ 1; . . . ; n;

ðD4Þ

whereas for the free point:

j ¼ 1; . . . ; n:

ðD5Þ

In all cases, the mass matrix consists of elements given by:

mij ¼ Lj ðxi Þ; Let us assume that a bar of length L is fixed at its left end (x = 0) whereas the other end (x = L) is allowed to perform free axial vibrations, of which the exact eigenvalues are given by the formula:

x¼L

BC

knj ¼ L00j ðxn Þ þ L0j ðxn Þ;

ð2j  1Þ

j ¼ 1; . . . ; n  1;

Now we deal with the possibility not to erase the Neumann end. In the early paper of 2008, the case of applying fulfilling the ODE at the free end has been dealt but not in a very successful way (see, [28, p.919]). The second case was later proposed by Auricchio et al. [40] in the framework of 2D and 3D problems, and it is here presented in full detail for this particular case of the 1D bar. In more details, the equations of stress equilibrium (governing ODE) and the Neumann boundary condition (BC) at the free end are added in parts (condensation), so the following equation is obtained:

EA 0 €Þ þ u ¼ 0;  ðEAu00  qAu |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |ffl{zffl} L

the inverse of the Jacobian matrix (for details we refer to [29]).

x2j ¼

i ¼ 2; . . . ; n  1;

ðD2Þ


2 2
2 2 2 @ 2 /j @n @ /j @n @ g @ /j @ g @ /j @ 2 n @/j ¼ þ 2 þ þ @y @y @y @n@ g @y2 @y @ g2 @y2 @n @n2 þ

L0j ðLÞ

i ¼ 2; . . . ; n;

j ¼ 1; . . . ; n;

ðD6Þ

Based on the above considerations, for n = 11 (i.e., Lagrange polynomials of degree p = 10), the results of Table 5 are obtained. There one may observe that the technique of eliminating the Neumann DOF un leads to smaller errors, whereas the technique of fulfilling the ODE and, at the same time, the Neumann boundary condition fails after the fourth mode and leads to complex eigenvalues, at least in the case of Lagrange polynomials. Similar findings were obtained for smaller and larger values of n. References [1] Coons SA. Surfaces for computer aided design of space form. Project MAC. MIT; 1964 [revised for MAC-TR-41 (1967). Available by CFSTI, Sills Building, 5285 Port Royal Road, Springfield, VA, USA]. [2] Gordon WJ, Hall CA. Transfinite element methods blending function interpolation over arbitrary curved element domains. Numer Math 1973;21:109–29.

C.G. Provatidis / Computers and Structures 182 (2017) 55–73 [3] Cavendish JC, Gordon WJ, Hall CA. Ritz-Galerkin approximations in blending function spaces. Numer Math 1976;26:155–78. [4] El-Zafrany A, Cookson RA. Derivation of Lagrangian and Hermitian shape functions for quadrilateral elements. Int J Numer Meth Eng 1986;23 (10):1939–58. [5] El-Zafrany A, Cookson RA. Derivation of Lagrangian and Hermitian shape functions for triangular elements. Int J Numer Meth Eng 1986;23(2):275–85. [6] Zhaobei Z, Zhiqiang X. Coon’s surface method for formulation of finite element of plate and shells. Comput Struct 1987;27:79–88. [7] Kanarachos A, Deriziotis D. On the solution of Laplace and wave propagation problems using C-elements. Finite Elem Anal Des 1989;5:97–109. [8] Kanarachos A, Provatidis C, Deriziotis D, Foteas N. A new approach of the FEM analysis of two-dimensional elastic structures using global (Coon’s) interpolation functions. In: Wunderlich J, editor. Proceedings, European conference on computational mechanics. Munich (Germany). [9] Provatidis CG. A review on attempts towards CAD/CAE integration using macroelements. Comput Res 2013;1(3):61–84. [10] Farin G, Hoschek J, Kim M-S. Handbook of computer aided geometric design. Amsterdam: North-Holland/Elsevier; 2002. [11] Provatidis C, Kanarachos A. Performance of a macro-FEM approach using global interpolation (Coons’) functions in axisymmetric potential problems. Comput Struct 2001;79:1769–79. [12] Provatidis CG. Analysis of box-like structures using 3-D Coons’ interpolation. Commun Numer Methods Eng 2005;21(8):443–56. [13] Provatidis CG. Three-dimensional Coons macroelements: application to eigenvalue and scalar wave propagation problems. Int J Numer Meth Eng 2006;65:111–34. [14] Provatidis CG. Coons-patch macroelements in two-dimensional parabolic problems. Appl Math Model 2006;30(4):319–51. [15] Provatidis C. Solution of two-dimensional Poisson problems in quadrilateral domains using transfinite Coons interpolation. Commun Numer Methods Eng 2004;20(7):521–33. [16] Provatidis C. Free vibration analysis of two-dimensional structures using Coons-patch macroelements. Finite Elem Anal Des 2006;42(6):518–31. [17] Provatidis CG. Two-dimensional elastostatic analysis using Coons-Gordon interpolation. Meccanica 2012;47(4):951–67. [18] Provatidis CG. Bézier versus Lagrange polynomials-based finite element analysis of 2-D potential problems. Adv Eng Softw 2014;73:22–34. [19] Höllig K. Finite element methods with B-splines. Philadelphia: SIAM; 2003. [20] Cottrell JA, Hughes TJR, Bazilevs Y. Isogeometric analysis: towards integration of CAD and FEA. Chichester: John Wiley; 2009. [21] Zieniuk E, Szerszen K. Triangular Bézier surface patches in modeling shape of boundary geometry for potential problems in 3D. Eng Comp 2013;29 (4):517–27. [22] Provatidis CG. Transient elastodynamic analysis of two-dimensional structures using Coons-patch macroelements. Int J Solids Struct 2006;43:6688–706.

73

[23] Provatidis CG. Global collocation method for 2-D rectangular domains. J Mech Mater Struct 2008;3:185–94. [24] Provatidis CG. Integration-free Coons macroelements for the solution of 2-D Poisson problems. Int J Numer Meth Eng 2009;77:536–57. [25] Provatidis CG, Ioannou KS. Static analysis of two-dimensional elastic structures using global collocation. Arch Appl Mech 2010;80:389–400. [26] Provatidis CG. B-splines collocation eigenanalysis of 2-D acoustic problems. J Mech Mater Struct 2014;9(3):259–85. [27] Provatidis CG. Free vibration analysis of elastic rods using global collocation. Arch Appl Mech 2008;78:241–50. [28] Provatidis CG. Time- and frequency-domain analysis using lumped mass global collocation. Arch Appl Mech 2008;78:909–20. [29] Filippatos A. Derivation of eigenfrequencies in acoustic cavities and elastic structures using the Global Collocation Method, Diploma Thesis. National Technical University of Athens, Department of Mechanical Engineering, October 2010 (supervised by C.G. Provatidis). [30] DeBoor C. A practical guide to splines. revised ed. New York: Springer-Verlag; 2001. [31] Piegl L, Tiller W. The NURBS book. Berlin: Springer-Verlag; 1995. [32] Zienkiewicz OC. The finite element method. London: McGraw-Hill; 1977. [33] DeBoor C. Private communication. November to December 2007. [34] Auricchio F, Beirão da Veiga L, Hughes TJR, Reali A, Sangalli G. Isogeometric collocation methods. Math Mod Meth Appl Sci 2010;20:2075–107. [35] Johnson RW. A B-spline collocation method for solving the incompressible Navier-Stokes equations using an ad hoc method: the Boundary Residual method. Comput Fluids 2005;34:121–49. [36] Johnson RW. Higher order B-spline collocation at the Greville abscissae. Appl Numer Math 2005;52:63–75. [37] Timoshenko SP, Young OH, Weaver W. Vibration problems in engineering. 4th ed. New York: Wiley; 1974. [38] Blevins RD. Formulas for natural frequency and mode shape. New York: Van Nostrand; 1979. [39] Kwok WY, Moser RD, Jiménez J. A critical evaluation of the resolution properties of B-spline and compact finite difference methods. J Comput Phys 2001;174:510–51. [40] Auricchio F, Beirão da Veiga L, Hughes TJR, Reali A, Sangalli G. Isogeometric collocation for elastostatics and explicit dynamics. Comput Meth Appl Mech Eng 2012;249–252:2–14. [41] Anitescu C, Jia Y, Zhang Y, Rabczuk T. An isogeometric collocation method using superconvergent points. Comput Meth Appl Mech Eng 2015;284:1073–97. [42] De Lorenzis L, Evans JA, Hughes TJR, Reali A. Isogeometric collocation: Neumann boundary conditions and contact. Comput Meth Appl Mech Eng 2015;284:21–54. [43] Zhang X, Song KZ, Lu MW, Liu X. Meshless methods based on collocation with radial basis functions. Comput Mech 2000;26:333–43.