Novel higher order mass matrices for isogeometric structural vibration analysis

Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Novel higher order mass matrices for isogeometric structural vibration analysis Dongdong Wang ⇑, Wei Liu, Hanjie Zhang Department of Civil Engineering, Xiamen University, Xiamen, Fujian 361005, China

a r t i c l e

i n f o

Article history: Received 21 January 2013 Received in revised form 15 March 2013 Accepted 21 March 2013 Available online 3 April 2013 Keywords: Isogeometric analysis Structural vibration Reduced bandwidth mass matrix Higher order mass matrix Frequency error

a b s t r a c t A set of novel higher order mass matrices are presented for isogeometric analysis of structural vibrations using NURBS. The proposed method for the construction of higher order mass matrices contains two essential steps. Firstly based upon the standard consistent mass matrix a special reduced bandwidth mass matrix is designed. This reduced bandwidth mass matrix meets the requirement of mass conservation while simultaneously preserves the same order of frequency accuracy as the corresponding consistent mass matrix. Subsequently a mixed mass matrix is formulated through a linear combination of the reduced bandwidth mass matrix and the consistent mass matrix. The desired higher order mass matrix is then deduced from the mixed mass matrix by optimizing the linear combination parameter to achieve the most favorable order of accuracy. Both quadratic and cubic formulations are discussed in detail and it is shown that with regard to the vibration frequency, the proposed higher order mass matrices have 6th and 8th orders of accuracy in contrast to the 4th and 6th orders of accuracy associated with the quadratic and cubic consistent mass matrices. A generalization to two dimensional higher order mass matrix is further realized by the tensor product operation on the one dimensional reduced bandwidth and consistent mass matrices. A series of benchmark examples congruously demonstrate that the proposed higher order mass matrices are capable of achieving the theoretically derived optimal accuracy orders for structural vibration analysis. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction The landmark work on isogeometric analysis by Hughes et al. [1] has received significant attention during the past few years. This method employs the non-uniform rational B-splines (NURBS) as the basis functions and it is capable of seamlessly bridging the computer aided design and the finite element analysis. Exact geometry is preserved throughout the analysis regardless of model refinements. Moreover in isogeometric analysis the high order smoothing approximation can be constructed quite straightforwardly which is desirable of thin plate and shell or gradient elasticity analysis. An excellent summary of the main features and advantages of the isogeometric analysis over the conventional finite element method can be found in the monograph by Cottrell et al. [2]. After its inception, the isogeometric method has experienced very rapid advances and developments with a vast range of versatile applications such as the structural vibrations [3–9], flow simulations [10–14], fluid–structure interactions [15,16], shape optimization and design sensitivity analysis [17–19], smoothing contact formulations [20,21], higher order

⇑ Corresponding author. E-mail address: [email protected] (D. Wang). 0045-7825/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2013.03.011

differential equation problems [22–25], beam, plate and shell analysis [25–30], T-spline enhanced isogeometric analysis [31,32], PHT-based isogeometric analysis [33,34], NURBS-enhanced finite element method [35], NURBS-based meshfree method [36], isogeometric collocation analysis [37,38], isogeometric boundary element analysis [39,40], extended isogeometric analysis [41,42], etc. In this work we focus on the isogeometric structural vibration analysis. The structural vibration analysis was investigated in the isogeometric setting first by Cottrell et al. [3,4] and Reali [5]. They showed that the isogeometric analysis produces accurate frequency spectra than the typical higher order finite elements. The order of accuracy for isogeometric structural vibration analysis was derived in detail for the classical one dimensional rod problem and Euler beam problem. For the rod problem, it was proved that quadratic and cubic NURBS have 4th and 6th orders of convergence when consistent mass matrices are utilized, respectively. Subsequently Hughes et al. [6] further studied the dynamic performance of isogeometric analysis in comparison with the finite element method for structural vibrations and wave propagation problems. It was demonstrated that the NURBS-based isogeometric analysis has superior approximation properties and the entire spectrum converges with increasing order of approximation. Recently Shojaee et al. [7,8] performed isogeometric free vibration and buckling

93

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

analysis for thin plates and laminated composite plates with favorable solution accuracy. Meanwhile, superior performance was also observed by Thai et al. [9] in the static, free vibration, and buckling analysis of laminated composite Reissner–Mindlin plate using the isogeometric approach. It is noted that the aforementioned excellent dynamic performance is observed with respect to the consistent mass formulation, actually the lumped mass matrices obtained through the row-sum technique do not give satisfactory performance and the accuracy is limited to second order for both quadratic and cubic NURBS [3]. For the structural vibration analysis in the category of finite element method, the consistent mass matrix may not be optimal with regard to the vibration frequency and noticeable efforts were devoted to the development of higher order mass matrices [43–46]. A typical example is the higher order mass matrix of linear rod element obtained by averaging the consistent and lumped mass matrices [46]. Thus a natural question arises as that whether higher order mass matrices can be constructed for the isogeometric structural vibration analysis with high order basis functions, i.e., quadratic and cubic NURBS. In the present study we aim to address this issue and propose a rational method for the development of isogeometric higher order mass matrices. It is shown that a combination of two equal order mass matrices is critical to construct a higher order mass matrix by optimizing the combination parameter with respect to the frequency error. Thus a direct combination of the consistent and lumped mass matrices are not desirable since as mentioned previously, the order of accuracy for the consistent mass matrix elevates with increasing order of approximation, i.e., 4th order and 6th order for the quadratic and cubic NURBS, respectively, while for both cases, the lumped mass matrices only have 2nd order of accuracy [3]. A two-step method is proposed here to develop higher order mass matrices for isogeometric analysis. Firstly the consistent mass matrix is recast into a new mass matrix with smaller bandwidth and adjustable parameters. This reduced bandwidth mass matrix is designed to satisfy the mass conservation requirement, i.e., the sum of its all elements is equal to that of the consistent mass matrix. Following the standard procedure of frequency error analysis [3,5,46], the optimal values for the adjustable parameters in the reduced bandwidth mass matrix are obtained by maximizing the order of accuracy, which turns out to be equal to the accuracy order of the corresponding consistent mass matrix. Thereafter a mixed mass matrix is proposed by linearly combining the consistent and the reduced bandwidth mass matrices. Once again based the frequency error estimation, the combination parameter in the mixed mass matrix is optimized to achieve the highest order of accuracy. The mixed mass matrix with the optimized combination parameter is the desired higher order mass matrix that is shown to be two orders more accurate than its corresponding consistent mass matrix in the isogeometric structural vibration analysis. Both quadratic and cubic higher order mass matrices are discussed in detail for the classical rod model. Extension to two dimensional formulation via the tensor product operation is presented as well. The efficacy of the proposed method is verified through benchmark examples. The remainder of this paper is organized as follows. The isogeometric formulation for structural vibration analysis is summarized in Section 2, in which explicit mass and stiffness matrices are emphasized for subsequent development. Section 3 gives a detailed presentation of the development of one dimensional higher order mass matrix for rod problem. A generalization to the two dimensional higher order mass matrix for membrane vibration is elaborated in Section 4. In Section 5, the accuracy of the proposed higher order mass matrices are demonstrated in a series of vibration problems. Finally conclusions are drawn in Section 6.

2. Isogeometric formulation for structural vibration analysis 2.1. B-spline and NURBS basis functions A B-spline curve is a linear combination of the position vectors or control points and their related B-spline basis functions of a parameter n. For convenience of presentation, one often chooses n 2 [0, 1] and defines a knot vector kn to construct a set of n B-spline basis functions of order p, denoted by Nap. The knot vector kn takes the following form:

kn ¼ fn1 ¼ 0; . . . ; na ; . . . ; nnþpþ1 ¼ 1gT

ð1Þ

It is seen that kn has a length of (n + p + 1) and contains a monotonically increased series of real numbers with na 6 nb . A knot vector is said to be uniform if its knot values are evenly spaced, otherwise it is non-uniform. Moreover a knot vector kn is said to be open when the multiplicity of its first and last knot values is equal to (p + 1). With a given knot vector of kn, a set of n p-th order B-spline basis functions Nap are uniquely defined by the Cox-de Boor recursion formula [1,47]:

Na0 ðnÞ ¼



1;

na 6 n < naþ1

0;

n < na or n P naþ1

for p ¼ 0

ð2Þ

and

Nap ðnÞ ¼

naþpþ1  n n  na Naðp1Þ ðnÞ þ Nðaþ1Þðp1Þ ðnÞ for p P 1 naþp  na naþpþ1  naþ1 ð3Þ

The B-spline basis functions Nap’s defined by Eqs. (2) and (3) have a local compact support of [na, na+p+1], i.e., a typical Nap spans (p + 1) knot intervals or elements. The continuity of the B-spline basis functions are controlled by its order and the multiplicity m of the knot values in the knot vector kn. Generally Nap is Cpm continuous and thus it has an interpolatory characteristic if an interior knot is repeated p times. Meanwhile, the basis functions do interpolate the end values when an open knot vector is used. It is also noticed that the B-spline basis functions form a partition of unity as follows: n X Nap ðnÞ ¼ 1;

npþ1 6 n 6 nnþ1

ð4Þ

a¼1

Eq. (4) states that for an open knot vector, the partition of unity is met in the whole domain n 2 [0, 1] due to the (p + 1) multiplicity of the first and last knots. On the other hand the partition of unity condition is valid only for n 2 [np+1, nn+1] when a uniform knot vector is employed. Fig. 1 displays the recursive construction of cubic B-splines basis functions based on the following two types of knot vectors: (a) Uniform knot vector: uniform

kn

¼ fn1 ; n2 ; . . . ; n8 ; n9 gT  T 1 1 3 1 5 3 7 ¼ 0; ; ; ; ; ; ; ; 1 8 4 8 2 8 4 8

ð5Þ

(b) Open knot vector: open

kn

¼ fn1 ; n2 ; . . . ; n12 ; n13 gT  T 1 1 3 1 5 3 7 ¼ 0; 0; 0; 0; ; ; ; ; ; ; ; 1; 1; 1; 1 8 4 8 2 8 4 8

ð6Þ

It is observed from Fig. 1 that the B-spline basis functions generally have one maximum and are Cp1 continuous for non-repeated knots. Moreover the basis functions have a non-

94

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

Fig. 1. Recursive construction of B-spline basis functions: (a) uniform knot vector; (b) open knot vector.

interpolatory nature and for an open knot vector the resulting basis functions interpolate the first and last values. The basis functions using both uniform knot vector and open knot vector are identical and periodic for the internal uniform knots values, among them one can be obtained by translating or shifting another one. For convenience of subsequent development, following Eqs. (2) and (3) the explicit expressions of typical quadratic and cubic basis functions in a uniformly spaced knot vector are given by:

N a2 ðnÞ ¼

8 1 ðn  na Þ2 na 6 n < naþ1 > > 2h2n > > > > 1 < 2 ½ðn  na Þðnaþ2  nÞ þ ðnaþ3  nÞðn  naþ1 Þ naþ1 6 n < naþ2 2h n

> 1 > ðnaþ3  nÞ2 > > 2h2n > > : 0

naþ2 6 n < naþ3 otherwise ð7Þ

8 1 ðn  na Þ3 > > 6h3n > > " # > > > > 1 ðn  na Þ2 ðnaþ2  nÞ þ ðn  na Þðn  naþ1 Þðnaþ3  nÞ > > > 3 > 2 > > 6hn þðn  naþ1 Þ ðnaþ4  nÞ < " # 2 N a3 ðnÞ ¼ > 1 ðn  na Þðnaþ3  nÞ þ ðn  naþ1 Þðnaþ3  nÞðnaþ4  nÞ > > 6h3n > > þðnaþ4  nÞ2 ðn  naþ2 Þ > > > > 1 > > 3 ðnaþ4  nÞ3 > 6h > > : n 0

na 6 n < naþ1 naþ1 6 n < naþ2 naþ2 6 n < naþ3 naþ3 6 n < naþ4 otherwise ð8Þ

where hn is the length of the knot interval, i.e., na+i = na + ihn. By further introducing a weight wa into the B-spline basis function Nap, one dimensional non-uniform rational B-spline (NURBS) basis function Rpa ðnÞ is defined as:

Rpa ðnÞ

Nap ðnÞxa ¼ Pn b¼1 N bp ðnÞxb

ð9Þ

In two dimensional case, let the parametric coordinates be (n, g), by using the tensor product operation on two one dimensional Bspline basis functions Nap(n) and Nbq(g), a two dimensional NURBS basis function Rpq ab ðn; gÞ can be defined as:

Nap ðnÞNbq ðgÞwab Pn Pm Rpq ab ðn; gÞ ¼ c¼1 d¼1 N cp ðnÞN dq ðgÞwcd

ð10Þ

in which wab is the weight for geometry control. Nbq(g) is the q-th order basis function and m is the number of basis functions in the g dimension, respectively. According to this standard construction procedure, the NURBS basis functions and the B-spline basis functions would be identical when constant weights are used. The NURBS basis functions share the same continuity and locally compact support properties as those of the B-spline basis functions, i.e., Rpq ab ðn; gÞ has a non-negative value on the domain of [na, na+p+1]  [gb, gb+q+1],  is the dyadic symbol. Similarly, an element in the parametric space is defined by [na, na+1]  [gb, gb+1] and the corresponding physical element can be obtained through isogeometric mapping. Obviously the NURBS basis functions in Eqs. (9) and (10) constitute a partition of unity as well: n X Rpa ðnÞ ¼ 1;

n X m X Rpq ab ðn; gÞ ¼ 1

a¼1

a¼1 b¼1

ð11Þ

In order to analytically derive explicit expressions for novel high order mass matrices, we mainly focus on the periodic B-spline basis functions or the NURBS basis functions with unity weight. 2.2. Isogeometric structural vibration analysis 2.2.1. Model problems Here we consider the longitudinal vibration of one dimensional rod problem and the transverse vibration of two dimensional membrane problem as the model problems for subsequent development. (1) 1D rod problem. The equation of motion for an elastic rod is:

€ ðx; tÞ ¼ c2 u;xx ðx; tÞ; u

c¼

pffiffiffiffiffiffiffiffiffi E=q

ð12Þ

where u is the axial displacement, t denotes the time, q is the material density per unit length, E is the Young’s modulus and c is the wave speed. Assume a harmonic solution for u(x, t), i.e.,

uðx; tÞ ¼ u0 eiðkxxtÞ ;

pffiffiffiffiffiffiffi

i ¼ 1

ð13Þ

where u0 is the wave amplitude, k is the wave number, x is the circular frequency. Then substituting Eq. (13) into Eq. (12) yields the frequency of the continuum problem:

95

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

x ¼ kc

ð14Þ

(2) 2D membrane problem. The governing equation for the transverse vibration of an elastic membrane is:

€ ðx; tÞ ¼ c2 ½u;xx ðx; tÞ þ u;yy ðx; tÞ u

ð15Þ

with dab(t) being the unknown coefficient associated with the ab-th control point, R and d are the vectors consist of NURBS basis functions and control point coefficients. Interpreting the unknown field variable in Eq. (19) as displacement and standard arguments lead to the general discrete equations for structural vibration analysis:

In this case, u(x, t) is the time dependent transverse displacement of membrane structure. Similarly a harmonic solution of u(x, t) takes the following form:

€ þ Kd ¼ 0 Md

uðx; y; tÞ ¼ u0 eiðkx xþky yxtÞ

M ¼ A ½M e ;

ð16Þ

where kx and ky are the wave numbers in x and y directions, respectively. Plugging Eq. (16) into Eq. (15) gives the continuum frequency x for the membrane vibration:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

x ¼ c k2x þ k2y

ð17Þ

with nel

e¼1

xðnÞ ¼

n X m X Rab ðnÞxab ;

x ¼ fx; yg;

n ¼ fn; gg

ð18Þ

a¼1 b¼1

where xab denotes the control point for a given geometry, and n and m are the number of control points in n and g directions, respectively. Thus the total number of control points are mn. For brevity the subscripts of the NURBS basis functions are dropped. With Eq. (18), the problem domain X is decomposed into the knot spans of the NURBS basis functions, which are referred as the elements Xe, e = 1, 2, . . . , nel in isogeometric analysis. The essential step in the isogeometric analysis is to discretize the unknown field variable u(x, t) using the same basis functions as those of geometric representation in Eq. (18):

uh ðx; tÞ ¼

n X m X

Rab ðnÞdab ðtÞ ¼ Rd

nel

nel

K ¼ A ½K e ;

e

d ¼ A ½d 

e¼1

e¼1

Z Xe

ReT qRe dX;

Ke ¼

Z Xe

rReT ErRe dX

(1) Quadratic basis function:

2

3 6 13 1 qAh 6 7 Me ¼ 4 13 54 13 5 120 1 13 6

ð19Þ

x a (b+ 2)

x (a + 1)(b+ 2) x (a + 2)(b+ 2)

x (a - 2)(b+ 1) x (a - 1)(b+ 1)

x a (b+ 1)

x (a + 1)(b+ 1) x (a + 2)(b+ 1)

x (a - 2)b

x ab

x (a + 1)b

x (a - 2)(b- 1) x (a - 1)(b- 1)

x a (b- 1)

x (a + 1)(b- 1) x (a + 2)(b- 1)

x (a - 2)(b- 2) x (a - 1)(b- 2)

x a (b- 2)

x (a + 1)(b- 2) x (a + 2)(b- 2)

hy

hx

x (a - 1)b

ð22Þ

where Re is the element basis function vector and r is the gradient operator, r = d/dx for the 1D rod problem and r = {o/ox, o/oy}T for the 2D membrane problem, E is the material Young’s modulus. Next for convenience of the subsequent development the element mass and stiffness matrices are presented in detail. For an elastic rod with a cross section area A, by using the quadratic and cubic basis functions of Eqs. (7) and (8) that are also shown in Fig. 1, let h denote the element size, the corresponding element mass and stiffness matrices are given by:

a¼1 b¼1

x (a - 2)(b+ 2) x (a - 1)(b+ 2)

ð21Þ

where A is the assembly operator [46], de is the element vector containing the coefficients of the control points related a typical element e, Me and Ke are the element mass and stiffness matrices given by:.

Me ¼ 2.2.2. Isogeometric discrete formulation With the NURBS basis functions, a generic field point xðnÞ in a typical problem domain X is given by:

ð20Þ

x (a + 2)b

Fig. 2. Illustration of control points, control net and elements for 2D quadratic NURBS discretization.

ð23Þ

96

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

R (a - 1)(b + 1)

R a (b + 1)

(7)

R (a + 1)(b + 1)

(8)

R (a - 1)b

(9)

R (a + 1)b

R ab

(4)

(5)

R (a - 1)(b- 1)

(6)

R a (b- 1)

(1)

R (a + 1)(b- 1)

(2)

(3)

Fig. 3. 2D quadratic NURBS basis functions defined in an element.

2 Ke ¼

2

1 1

3

EA 6 7 4 1 2 1 5 6h 1 1 2

ð24Þ

Ke ¼

(2) Cubic basis function:

2

Me ¼

20

129

60

933

60 7 7 7 1188 129 5

1

60

129

933

ð25Þ

20 3

2

6 7 12 1 6 7 34 29 12 7 EA 7 6 Ke ¼ 7 6 120h 4 12 29 34 7 5 1

12

7

ð26Þ

6

In case of the 2D membrane problem, we consider a quasi-uniform discretization with rectangular elements as shown in Fig. 2, the element sizes in x and y directions be hx and hy, respectively. Then each element is associated with nine basis functions as shown in Fig. 3, the corresponding mass matrix is:

2

36

6 78 6 6 6 6 6 6 6 78 hx hy 6 e 6 169 M ¼ qs 14400 6 6 6 13 6 6 6 6 6 4 13 1

78

6

6

13

324 78 169 702 169 13

54

78

13

36

78 13

169

Es ½K es1 K es2 ; 720u

u¼

hy hx

ð28Þ

with the sub-matrices Kes1 and Kes2 being given by:

3

1

6 5040 4 60

qAh 6 6 129 1188

where s is the membrane thickness. The stiffness matrix takes the following form:

13

169

78

1

169 13 324 702

54

78 169

702 169 702 2916 702 169 702 169 78

54

702 324 13 169

13

1

78

169

54 13

13 169 702 169 78 324 6 13 169 78 6 78

13

36

78

1

2

12u2 þ 12 6 6 26u2  6 6 6 6 2 u2  6 6 6 6 6u2 þ 26 6 6 es1 2 K ¼6 6 13u  13 6 2 6 u  13 6 6 6 6u2 þ 2 6 6 6 13u2  1 4  u2  1

ð27Þ

2 u2  6

108u2 þ 12

26u2  6

26u2  6

12u2 þ 12

2

13u  13

u2  13

2

54u þ 26 13u2  13 13u2  13 6u2 þ 26 13u2  1

 u2  1

2

13u2  1

2

6u2 þ 2

54u þ 2 13u  1

6u2 þ 26

13u2  13

u2  13

13u2  13

3

7 13c2  13 54u2 þ 26 7 7 7 u2  13 13u2  13 7 7 7 12u2 þ 108 26u2  54 7 7 7 26u2  54 108u2 þ 108 7 7 7 2u2  54 26u2  54 7 7 7 6u2 þ 26 13u2  13 7 7 7 13u2  13 54u2 þ 26 7 5 ð29Þ

2

3

13 7 7 7 6 7 7 7 13 7 7 169 7 7 7 78 7 7 6 7 7 7 78 5 36

26u2  6

K es2

u2  13

6u2 þ 2

13u2  1

u2  1

3

7 6 6 13u2  13 13u2  1 54u2 þ 2 13u2  1 7 7 6 7 6 6 6u2 þ 26 u2  1 13u2  1 6u2 þ 2 7 7 6 7 6 6 2u2  54 6u2 þ 26 13u2  13 u2  13 7 7 6 7 6 ¼ 6 26u2  54 13u2  13 54u2 þ 26 13u2  13 7 7 6 7 6 u2  13 13u2  13 6u2 þ 26 7 6 12u2 þ 108 7 6 6 u2  13 12u2 þ 12 26u2  6 2u2  6 7 7 6 7 6 6 13u2  13 2 2 2 26u  6 108u þ 12 26u  6 7 5 4 6u2 þ 26

2u2  6

26u2  6

12u2 þ 12 ð30Þ

97

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

3. Construction of 1D higher order mass matrices

Finally the order of accuracy for the proposed algorithm can be deduced from Eq. (37) as:

We first present high order mass matrices for 1D rod problem, both quadratic cubic basis functions are considered herein. Without loss of generality, in the following presentation the 1D rod is assumed to have a unit cross section area A = 1.

xh 4b 18  b 1 ðkhÞ2 þ ðkhÞ4 240 2880 x

ð39Þ

Based on Eq. (39), clearly bopt = 4 leads to 4th order accuracy: h

x 7 1þ ðkhÞ4 1440 x

3.1. Quadratic basis function The consistent mass matrix for the quadratic basis function is shown in Eq. (23), to construct the higher order mass matrix, here we first propose the following new element mass matrix Mer:

2 M er ¼

7  b=2

qh 6 4 13 þ b=2

120

0

13 þ b=2 54  b 13 þ b=2

3

0 7 13 þ b=2 5 7  b=2

ð32Þ

h

da ðx; tÞ ¼ d0 eiðkxa x tÞ €a ðx; tÞ ¼ ðxh Þ2 d0 eiðkxa xh tÞ d

ð33Þ

where d0 is the amplitude of the harmonic wave, xh represents the discrete frequency. Then the solution da+j(x, t) for the neighboring control point xa+j = xa + jh is:

ð34Þ

aþj

Plugging Eqs. (33) and (34) into Eq. (32) gives: 2

ðxh hÞ ½ð26 þ bÞ cosðkhÞ þ 34  b  c2 ½3  cosð2khÞ  2 cosðkhÞ ¼ 0 20 ð35Þ where the identity of eikh + eikh = 2coskh is employed. Thus the discrete frequency xh can be obtained as:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xh 1 20½3  cosð2khÞ  2 cosðkhÞ ¼ ð26 þ bÞ cosðkhÞ þ 34  b x kh

ð36Þ

1þ

#g 2#

ð37Þ

20 ðkhÞ2

0

15

5

3

120

ð41Þ

Compared with the consistent mass matrix, this new mass matrix has the same order of accuracy but a relatively smaller halfbandwidth of 2. Next we show a novel higher order mass matrix can be rationally established by combining the consistent mass matrix Me of Eq. (23) and the reduced bandwidth mass matrix Mer of Eq. (41). Consider the following mixed combination of the consistent and reduced bandwidth mass matrices:

2

M

em

3

5þc 15  2c c qh 6 7 ¼ ð1  cÞM þ cM ¼ 15  2 c 50 þ 4 c 15  2c 5 4 120 15  2c 5þc c er

e

with c being a parameter in order to optimize the order of accuracy. With the mass matrix of Mem in Eq. (42), the discrete equation of motion at a generic control point xa becomes:

i h h € cda2 þ ð30  4cÞd€a1 þ ð60 þ 6cÞd€a þ ð30  4cÞd€aþ1 þ cd€aþ2 120 c2  ½da2 þ 2da1  6da þ 2daþ1 þ daþ2  6h ¼0 ð43Þ Further substituting Eq. (34) into (43) yields:

xh 1 20½3  cosð2khÞ  2 cosðkhÞ ¼ x kh c cosð2khÞ þ ð30  4cÞ cosðkhÞ þ 30 þ 3c

ð44Þ

which also can be expressed as the form of Eq. (37) with # and g given by:

(

43 # ¼ 60  15ðkhÞ2 þ 11 ðkhÞ4  336 ðkhÞ6 6

g ¼ 60  15ðkhÞ2 þ 2c4þ5 ðkhÞ4  2c24þ1 ðkhÞ6

ð45Þ

Therefore according to Eq. (37) the order of accuracy for the mixed mass matrix is:

xh 7  6c 29  28c 1þ ðkhÞ4  ðkhÞ6 x 1440 40320

ð46Þ

M eho ¼

1 ð7M e  M er Þ 6

ð47Þ

which has a 6th order of accuracy:

with

8 <

0

Consequently the optimal value of c is copt = 7/6 and it corresponds to the following new higher order quadratic mass matrix as indicated Eq. (42):

where the continuum frequency of Eq. (14) is used. Further invoking the Taylor’s expansion leads to [3]:

g

15

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

daþj ðx; tÞ ¼ eiðkjhÞ da ðx; tÞ €a ðx; tÞ € ðx; tÞ ¼ eiðkjhÞ d d

sffiffiffi #

5

ð42Þ

Let the discrete nodal solution take the following harmonic form:

xh  x

2

ð31Þ

i h h € € € ð26 þ bÞd a1 þ ð68  2bÞda þ ð26 þ bÞdaþ1 120 c2  ½da2 þ 2da1  6da þ 2daþ1 þ daþ2  6h ¼0

(

Up to this stage we construct a new reduced band mass matrix of Eq. (31) with bopt = 4, i.e.,

qh 6 7 M er ¼ 4 15 50 15 5

where b is a parameter to be adjusted for desired solution accuracy. Clearly the proposed mass matrix of Eq. (31) satisfies the requirement of mass conservation but with a reduced half-bandwidth by 1 compared with the consistent mass matrix given by Eq. (23). Next we show this reduced band mass matrix is capable of achieving the same order of accuracy as that of the consistent mass matrix. Based upon the stiffness of Eq. (24) and the proposed reduced bandwidth mass matrix of Eq. (31), the discrete equation of motion for a typical control point xa is:

(

ð40Þ

2

xh 11 1þ ðkhÞ6 120960 x

4

½3  cosð2khÞ  2 cosðkhÞ  60  15ðkhÞ þ 11 ðkhÞ  # 6

: ð26 þ bÞ cosðkhÞ þ 34  b  60  ð26þbÞ ðkhÞ2 þ 26þb ðkhÞ4  g 24 2 ð38Þ

ð48Þ

Meanwhile, letting c = 1 in Eq. (46) yields the order of accuracy for the quadratic consistent mass algorithm [3]:

98

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

xh 1 1þ ðkhÞ4 1440 x

ð49Þ

Thus compared with the consistent mass matrix of Eq. (23), the new higher order mass matrix of Eq. (47) improves the frequency accuracy by two orders. 3.2. Cubic basis function In case of the cubic basis function, based on the consistent mass matrix of Eq. (25), the following mass matrix with reduced halfbandwidth of 3 is proposed:

3 5 43 20 þ 604 h 129 þ 397 v 60 þ 12 b 0 6 43 v 1188 þ 297 h 933 þ 311 v 60 þ 12 b 7 qh 6 129 þ 397 7 604 397 M er ¼ 7 6 43 5040 4 60 þ 12 b 933 þ 311 v 1188 þ 297 h 129 þ 397 v5 397 604 2

0

60 þ 12 b

43 129 þ 397 v

5 20 þ 604 h

with a 6th order of accuracy:

xh 1 1 1 ðkhÞ6 þ ðkhÞ8 12096 36288 x

To achieve a higher order accurate algorithm, similar to the previous quadratic formulation we further propose the following mixed mass matrix Mem through a linear combination of the cubic consistent mass matrix of Eq. (25) and the reduced bandwidth mass matrix of Eq. (56): M em ¼ ð1  cÞM er þ cM e 2 3 cÞ cÞ 20 þ 25ð1 129  645ð1 63  3c c 151 397 6 7 6 7 cÞ cÞ cÞ 129  645ð1 1188 þ 1485ð1 933  4665ð1 63  3c 7 397 151 397 qh 6 6 7 ¼ 6 7 5040 6 63  3c cÞ 1485ð1cÞ 645ð1cÞ 7 933  4665ð1 1188 þ 129  6 397 151 397 7 4 5

c

63  3c

ð50Þ where b and v are the parameters to be optimized, h = (2  2b  2v). It is straightforward to show the corresponding stencil equation is:

"

€a2 þ ð1191 þ vÞd €a1 þ ð2418  2b  2vÞd €a ð120 þ bÞd h €aþ1 þ ð120 þ bÞd €aþ2 5040 þð1191 þ vÞd

#

c2  ðda3 þ 24da2 þ 15da1  80da þ 15daþ1 þ 24daþ2 þ daþ3 Þ ¼ 0 120h ð51Þ Subsequently substituting Eqs. (33) and (34) into Eq. (51) gives the following characteristic equation:

ðxh hÞ2 ½ð120 þ bÞ cosð2khÞ þ ð1191 þ vÞ cosðkhÞ þ ð1209  b  vÞ 42 2 ð52Þ þ c ½cosð3khÞ þ 24 cosð2khÞ þ 15 cosðkhÞ  40 ¼ 0 and thus the discrete frequency xh can be obtained as:

cÞ 129  645ð1 397

cÞ 20 þ 25ð1 151

ð58Þ

The stencil equation associated with the mixed mass matrix (58) and the standard stiffness matrix (25) reads:

" # cd€a3 þ ð126  6cÞd€a2 þ ð1176 þ 15cÞd€a1 þ ð2436  20cÞd€a h €aþ1 þ ð126  6cÞd €aþ2 þ cd €aþ3 5040 þð1176 þ 15cÞd 

c2 ðda3 þ 24da2 þ 15da1  80da þ 15daþ1 þ 24daþ2 þ daþ3 Þ ¼ 0 120h ð59Þ By virtue of Eqs. (33) and (34) we obtain:

" # ðxh hÞ2 c cosð3khÞ þ ð126  6cÞ cosð2khÞ 21 þð1176 þ 15cÞ cosðkhÞ þ ð1218  10cÞ

ð60Þ

2

þ c ½2 cosð3khÞ þ 48 cosð2khÞ þ 30 cosðkhÞ  80 ¼ 0 and

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

xh 1 ¼ x khsffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 

ð57Þ

42½40  cosð3khÞ  24 cosð2khÞ  15 cosðkhÞ ð120 þ bÞ cosð2khÞ þ ð1191 þ vÞ cosðkhÞ þ ð1209  b  vÞ ð53Þ

By standard arguments of Taylor’s expansion, Eq. (53) also reduces to the form of Eq. (37) in which # and g are:

8 > # ¼ 2520  840ðkhÞ2 þ 133ðkhÞ4  53 ðkhÞ6 þ 697 ðkhÞ8 > 4 720 < 1671    1 g ¼ 2520  2 þ 2b þ v2 ðkhÞ2 þ 1037 þ 23 b þ 24 v ðkhÞ4 8 > >  : 2957 4 v  v  2  240 þ 45 b þ 720 þ 315 b þ 40320 ðkhÞ6 þ 10637 ðkhÞ8 13440

xh 1 u u42½40  cosð3khÞ  24 cosð2khÞ  15 cosðkhÞ  ¼ u  x kh t c cosð3khÞ þ ð126  6cÞ cosð2khÞ þð1176 þ 15cÞ cosðkhÞ þ ð1218  10cÞ

ð61Þ

Thus through Taylor’s expansion Eq. (61) also reduces to Eq. (37) and the parameters # and g are now given by:

8 < # ¼ 2520  840ðkhÞ2 þ 133ðkhÞ4  53 ðkhÞ6 þ 697 ðkhÞ8 4 720     : g ¼ 2520  840ðkhÞ2 þ 133ðkhÞ4  2c þ 776 ðkhÞ6 þ 18 c þ 199 ðkhÞ8 240 ð62Þ

ð54Þ

Consequently the order of accuracy for the proposed mixed mass matrix is:

Therefore the order of accuracy for the proposed reduced bandwidth mass matrix is:

xh c  5=6 10  9c 1þ ðkhÞ6 þ ðkhÞ8 10080 x 362880

xh 9  4b  v 81  16b  v 1 ðkhÞ2 þ ðkhÞ4 x 10080 120960 669  64b  v 7121  256b  v 6

Therefore for the cubic basis function we can choose copt = 5/6 to construct a desired 8th order accurate cubic mass matrix:



3628800

ðkhÞ þ

203212800

ðkhÞ8

ð55Þ

in which choosing bopt = 6, vopt = 15 yields the desired reduced bandwidth mass matrix:

2

25 20 þ 151

129  645 397

63

0

3

7 645 1485 4665 63 qh 6 7 6 129  397 1188 þ 151 933  397 M er ¼ 7 6 4665 1485 645 5040 4 63 933  397 1188 þ 151 129  397 5 0

63

129  645 397

25 20 þ 151

ð56Þ

M eho ¼

1 ð5M e þ M er Þ 6

ð63Þ

ð64Þ

with the accuracy measured by

xh 1 1þ ðkhÞ8 145152 x

ð65Þ

At the same time, c = 1 in Eq. (63) also gives the order of accuracy for the consistent mass matrix formulation [3]:

xh 1 1þ ðkhÞ6 60480 x

ð66Þ

99

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

2

4. Generalization to 2D higher order mass matrix A direct extension of the previous methodology to multidimensional formulation is much more involved. Alternatively here we develop the higher order mass matrix via the tensor product approach by using the 1D high order mass matrices in each direction respectively. The motivation is the tensor product nature of the B-splines or NURBS basis functions. Consider a quasi-uniform discretization, i.e., rectangular mesh with a constant Jacobian for the isogeometric mapping as shown in Fig. 2, if all the weights of NURBS are unity and then NURBS are identical to the B-splines, under these conditions, the 2D mass matrix is equivalent to the tensor product of two 1D mass matrices in their respective directions, i.e.,

M eabcd ¼ ¼ ¼

Z ZX

e

ZX

e

qsRab ðn; gÞRcd ðn; gÞdX qsNa ðnÞN b ðgÞNc ðnÞNd ðgÞdX

¼ qsJ

Z Xen

Na ðnÞNc ðnÞdn

Z Xeg

e1g

ð67Þ

ða1Þðbþ2Þ

0

0 75 225 0 75 225 750 225

0 0

0 0

25

0

225

75

0

250 750

0

75 225

225 750 2500 750 225 750 75

0

750 250

0

225

0

75

225

25

75

0

225 750 225 75 250

0

0

225

14400

ð72Þ

0 75

0

75

3 0 0 7 7 7 0 7 7 7 0 7 7 225 7 7 7 75 7 7 0 7 7 7 75 5 25

½M ems1 M ems2 

ð70Þ

ða1Þðb1Þ

aðb2Þ

ðaþ1Þðb2Þ

aðb1Þ

ðaþ1Þðb1Þ

ðaþ1Þðbþ2Þ

ðaþ2Þðb2Þ

ðaþ2Þðb1Þ

€ € € þ 66cd ðaþ2Þb þ 26cdðaþ2Þðbþ1Þ þ cdðaþ2Þðbþ2Þ

ð73Þ

Gðdab Þ ¼ ðu2 þ 1Þdða2Þðb2Þ þ ð26u2 þ 2Þdða2Þðb1Þ þ ð66u2  6Þdða2Þb þ ð26u2 þ 2Þdða2Þðbþ1Þ þ ðu2 þ 1Þdða2Þðbþ2Þ þ ð2u2 þ 26Þdða1Þðb2Þ þ ð52u2 þ 52Þdða1Þðb1Þ þ ð132u  156Þdða1Þb þ ð52u2 þ 52Þdða1Þðbþ1Þ þ ð2u2 þ 26Þdða1Þðbþ2Þ þ ð6u2 þ 66Þdaðb2Þ þ ð156u2 þ 132Þdaðb1Þ þ ð396u2  396Þdab þ ð156u2 þ 132Þdaðbþ1Þ þ ð6u þ 66Þdaðbþ2Þ þ ð2u þ 26Þdðaþ1Þðb2Þ þ ð52u2 þ 52Þdðaþ1Þðb1Þ þ ð132u2  156Þdðaþ1Þb þ ð52u2 þ 52Þdðaþ1Þðbþ1Þ þ ð2u2 þ 26Þdðaþ1Þðbþ2Þ þ ðu2 þ 1Þdðaþ2Þðb2Þ þ ð26u2 þ 2Þdðaþ2Þðb1Þ þ ð66u2  6Þdðaþ2Þb þ ð26u2 þ 2Þuðaþ2Þðbþ1Þ þ ðu2 þ 1Þdðaþ2Þðbþ2Þ

ð69Þ

in which the sub-matrices Mems1 and Mems2 are given by: 3 2 25 þ 11c 75 þ 3c 6c 75 þ 3c 225  56c 6 75 þ 3c 250 þ 74c 75 þ 3c 225  56c 750  48c 7 7 6 7 6 7 6 6 c 75 þ 3 c 25 þ 11 c 13 c 225  56 c 7 6 7 6 13c 250 þ 74c 750  48c 7 6 75 þ 3c 225  56c 7 6 ems1 6 M ¼ 6 225  56c 750  48c 225  56c 750  48c 2500 þ 416c 7 7 7 6 225  56c 75 þ 3c 54c 750  48c 7 6 13c 7 6 6 6c 13c c 75 þ 3c 225  56c 7 7 6 7 6 54c 13c 225  56c 750  48c 5 4 13c c 13c 6c 13c 225  56c

ða1Þðb2Þ

€ € þ ð1800  84cÞd ðaþ1Þb þ ð900  224cÞdðaþ1Þðbþ1Þ € € € þ 26cd þ cd þ 26cd

Therefore a linear combination of the 2D consistent mass matrix of Eq. (27) and the 2D reduced bandwidth mass matrix of Eq. (68) gives the following 2D mixed mass matrix:

qshx hy

75 þ 3c

€ þ ð1800  84cÞd € € þ ð3600 þ 756cÞd ab aðbþ1Þ þ 66cdaðbþ2Þ € € þ 26cd þ ð900  224cÞd

ð68Þ

M em ¼ ð1  cÞM er þ cM ec ¼

250 þ 74c

€ € þ ð1800  84cÞd ða1Þb þ ð900  224cÞdða1Þðbþ1Þ € € € þ 26cd þ 66cd þ ð1800  84cÞd

where the material density q and the membrane thickness s are ase1g sumed to be constants. Me1n ac and M bd are the entries of the 1D element matrix discussed previously. The relationship of Eq. (67) is clearly demonstrated by the fact that the 2D mass matrix of Eq. (27) is the tensor product of two 1D mass matrices of Eq. (23) in n and g directions. In accordance with this philosophy, a 2D reduced bandwidth matrix can be built upon its 1D counterpart of Eq. (41) as follows:

0

225  56c 75 þ 3c

€ Þ  c2 Gðd Þ ¼ 0 Hðd ab ab

ða2Þðbþ2Þ

¼ qsJMe1n ac M bd

0

750  48c

€ Þ ¼ cd € € € € Hðd ab ða2Þðb2Þ þ 26cdða2Þðb1Þ þ 66cdða2Þb þ 26cdða2Þðbþ1Þ € € € þ cd þ 26cd þ ð900  224cÞd

e1g

0

13c 225  56c

6c

7 7 7 7 6c 7 7 13c 7 7 225  56c 7 7 7 75 þ 3c 7 7 7 6c 7 7 75 þ 3c 5 25 þ 11c 13c

Based on the mixed mass matrix of Eq. (69) and the stiffness matrix of Eq. (28) for the membrane structure with the discretization given in Fig. 2, the equation of motion associated with a typical control point xab = (xa, yb) takes the following form:

M bd

25 75 6 75 250 6 6 6 0 75 6 6 6 75 225 qshx hy 6 er 6 225 750 M ¼ 14400 6 6 6 0 225 6 6 0 0 6 6 0 4 0

c

54c

€ Þ and Gðd Þ are defined as: where Hðd ab ab

|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

2

13c

ð71Þ

20

Nb ðgÞNd ðgÞdg

M e1n ac

75 þ 3c

u2 h2x

qsNa ðnÞNb ðgÞNc ðnÞN d ðgÞJdndg e

Xn

M ems2

13c 6c 6 225  56c 13c 6 6 6 75 þ 3c c 6 6 54c 75 þ 3c 6 6 ¼6 6 750  48c 225  56c 6 13c 6 250 þ 74c 6 6 13c 25 þ 11c 6 6 75 þ 3c 4 225  56c

3

ð74Þ

For the 2D membrane vibration, similar to Eq. (16) the transverse deflectional coefficient dab related to the control point xab is assumed to have the following harmonic form:

(

dab ¼ dðxa ; yb ; tÞ ¼ d0 exp½iðkx xa þ ky yb Þ  ixh t € ¼ dðx € a ; y ; tÞ ¼ ðxh Þ2 d0 exp½iðkx xa þ ky y Þ  ixh t d ab b b

ð75Þ

where again xh represents the discrete frequency. Subsequently using of Eq. (75) yields the following expressions for the transverse deflectional coefficient d(a+j)(b+l) that is related the neighboring control point x(a+j)(b+l) = (xa + jhx, yb + lhy):

(

dðaþjÞðbþlÞ ¼ exp½iðkx jhx þ ky lhy Þdab € € d ðaþjÞðbþlÞ ¼ exp½iðkx jhx þ ky lhy Þdab

ð76Þ

100

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

Fig. 4. Comparison of the first three frequencies for the fixed–fixed and fixed-free rod problems using various mass matrices with periodic quadratic basis functions: (a)–(c)the fundamental, second, and third frequencies of the fixed–fixed rod; (d)–(f)-the fundamental, second, and third frequencies of the fixed-free rod.

 T ¼ ðu2  1Þcos 2ðkx hx þ ky hy Þ þ ð26u2  2Þcosð2kx hx þ ky hy Þ

Substituting Eqs. (75) and (76) into Eq. (72) gives:

ðxh uhx Þ2 S  c2 T ¼ 0 20

ð77Þ

þ ð66u2 þ 6Þcosð2kx hx Þ þ ð26u2  2Þcosð2kx hx  ky hy Þ þ ðu2  1Þcos½2ðkx hx  ky hy Þ þ ð2u2  26Þcosðkx hx þ 2ky hy Þ

with

þ ð52u2  52Þcosðkx hx þ ky hy Þ þ ð132u2 þ 156Þcosðkx hx Þ

S ¼ c cos½2ðkx hx þ ky hy Þ þ 26c cosð2kx hx þ ky hy Þ þ 66c cosð2kx hx Þ þ 26c cosð2kx hx  ky hy Þ þ c cos½2ðkx hx

þ ð52u2  52Þcosðkx hx  ky hy Þ þ ð2u2  26Þcosðkx hx  2ky hy Þ

 ky hy Þ þ 26c cosðkx hx þ 2ky hy Þ þ ð900  224cÞ

þ ð6u2  66Þcosð2ky hy Þ þ ð156u2  132Þcosðky hy Þ

 cosðkx hx þ ky hy Þ þ ð1800  84cÞ cosðkx hx Þ þ ð900

þ ð198u2 þ 198Þ

ð79Þ

 224cÞ cosðkx hx  ky hy Þ þ 26c cosðkx hx  2ky hy Þ þ 66c cosð2ky hy Þ þ ð1800  84cÞ cosðky hy Þ þ ð1800 þ 378cÞ

ð78Þ

Generally it is not straightforward to obtain an explicit error estimation for the discrete frequency xh from Eq. (77) as the previous 1D development due to the complexity of Eq. (77).

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

101

Fig. 5. Comparison of the first three frequencies for the fixed–fixed and fixed-free rod problems using various mass matrices with periodic cubic basis functions: (a)–(c)-the fundamental, second, and third frequencies of the fixed–fixed rod; (d)–(f)-the fundamental, second, and third frequencies of the fixed-free rod.

Consequently we further simplify the problem and consider a uniform discretization by setting with hx = hy = h, u = 1 and kx = ky = k, with these assumptions Eq. (77) becomes:

  ðxh hÞ2 c cosð4khÞ þ 52c cosð3khÞ þ ð900  92cÞ cosð2khÞ 20 þð3600  116cÞ cosðkhÞ þ 2700 þ 155c   2 cosð4khÞ  56 cosð3khÞ ¼0  c2 224 cosð2khÞ  8 cosðkhÞ þ 290

pffiffiffi in which since kx = ky = k, by using Eq. (17) we have x ¼ 2kc in this case. Once again after Taylor’s expansion Eq. (81) can be converted into the form of Eq. (37) where # and g are defined as:

(

# ¼ 14400  7200ðkhÞ2 þ 1760ðkhÞ4  1930 ðkhÞ6 g ¼ 14400 7 2 4  ð170 þ 100cÞðkhÞ6 7200ðkhÞ þ ð1500 þ 240cÞðkhÞ ð82Þ

ð80Þ Thus the discrete frequency xh is obtained as: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xh c u u40½ cosð4khÞ  28 cosð3khÞ  112 cosð2khÞ  4 cosðkhÞ þ 145   ¼ u c cosð4khÞ þ 52c cosð3khÞ þ ð900  92cÞ cosð2khÞ x xh t þð3600  116cÞ cosðkhÞ þ 2700 þ 155c ð81Þ

Following Eqs. (37) and (82), the order of accuracy for the 2D mixed mass matrix formulation is given by:

xh 13  12c 35c  37 1þ ðkhÞ4 þ ðkhÞ6 10080 x 1440

ð83Þ

Therefore selecting copt = 13/12 yields the following 2D higher order mixed mass matrix Meho for membrane structure:

102

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

Fig. 6. Comparison of the first three frequencies for the fixed–fixed and fixed-free rod problems using various mass matrices with open knot quadratic basis functions: (a)– (c)-the fundamental, second, and third frequencies of the fixed–fixed rod; (d)–(f)-the fundamental, second, and third frequencies of the fixed-free rod.

M eho ¼

1 ð13M ec  M er Þ 12

ð84Þ

with a 6th order of accuracy:

xh 11 1þ ðkhÞ6 120960 x

ð85Þ

On the other hand a 4th order of accuracy for the consistent mass matrix defined in Eq. (27) can be deduced from Eq. (83) by letting c = 1: h

x 1 1þ ðkhÞ4 1440 x

ð86Þ

In summary, similar to the 1D formulation the quadratic 2D higher order mass matrix of Eq. (84) also gains solution accuracy by two orders compared with its consistent counterpart of Eq.

(27). Of course it is noticed that at this time Meho in Eq. (84) is designed for the particular case of uniform discretization and equal wave numbers in the x and y directions, i.e., kx = ky = k. 5. Numerical results 5.1. 1D rod vibration To examine the proposed algorithm, the vibration of an elastic rod subjected to fixed–fixed, fixed-free boundary conditions is thoroughly studied. The analytical solutions for the fixed–fixed and fixed-free rod vibration problems are:

xi ¼



icp=L

fixed—fixed rod

ð2i  1Þcp=ð2LÞ fixed—free rod

ð87Þ

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

103

Fig. 7. Comparison of the first three frequencies for the fixed–fixed and fixed-free rod problems using various mass matrices with open knot cubic basis functions: (a)–(c)-the fundamental, second, and third frequencies of the fixed–fixed rod; (d)–(f)-the fundamental, second, and third frequencies of the fixed-free rod.

where L is the rod length. Without loss of generality, in this study the geometry and material properties for the elastic rod are chosen as: length L = 10, cross section area A = 1, material density q = 1, and Young’s modulus E = 1. Since the present higher order mass matrices are consistently derived based on a uniform and periodic NURBS discretization and do not take into account the boundary truncation effect induced by the open knot vector, we first consider the case that the rod is discretized by the periodic basis functions with uniform knot vector as shown in Fig. 1(a). It is noted that under this circumstance, the knots beyond the physical domain are needed to construct a periodic discretization on the physical domain of the rod problem. Appropriate boundary conditions can be effectively enforced via the transformation method proposed by Wang and Xuan [48]. Fig. 4(a)–(c) list the first three frequencies for the vibration of a

fixed–fixed elastic rod using quadratic basis functions, where three types of mass formulations are compared, i.e., the consistent mass matrix ‘‘CM’’, the reduced bandwidth mass matrix ‘‘RBM’’, and the higher order mass matrix ‘‘HOM’’. The results of the fixed-free vibration with quadratic basis functions are depicted in Fig. 4(d)– (f). These results agree very well with the orders of accuracy which are theoretically derived in this work, i.e., the orders of accuracy for the consistent, reduced bandwidth, and higher order mass matrices are 4, 4, and 6. In case of cubic basis functions, the convergence results for the first three frequencies are given in Fig. 5(a)–(c) for the fixed–fixed case and in Fig. 5(d)–(f) for the fixed-free case, respectively. These results again demonstrate that the proposed cubic higher order mass matrix has an expected 8th order of accuracy, which is an optimal combination of two 6th order accurate consistent and reduced bandwidth mass matrices.

104

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

Fig. 8. Periodic NURBS discretization of the square membrane problem: (a) control net; (b) physical mesh.

Fig. 9. Comparison of the first four frequencies for the square membrane problem using various mass matrices with periodic quadratic basis functions: (a) the fundamental frequency; (b) the second frequency; (c) the third frequency; (d) the fourth frequency.

In practice the open knot vector-based NURBS basis functions are often employed. In this case, the basis functions near the boundary are different from those in the interior domain as shown in Fig. 1(b), while the present formulation with the optimal combination parameters are designed for the periodic basis functions. Thus it is desirable to further investigate the performance of the proposed higher order mass formulation with regard to the open knot basis functions in order to comprehensively assess the effectiveness of the proposed higher order mass matrix formulation. In the computation, the proposed optimal parameters are directly

adopted to formulate the reduced bandwidth and the higher order mass matrix. The pervious fixed–fixed and fixed free rod vibration examples are re-examined. Fig. 6 plots the results of quadratic basis functions and the results of cubic formulations are shown in Fig. 7. It turns out that for the fixed–fixed rod problem, as shown in Figs. 6 and 7(a)–(c), the orders of accuracy for the proposed higher order mass matrices are in excellent accordance with the expected values, i.e., 6 for the periodic basis functions and 8 for the cubic basis functions. On the other hand, for the fixed-free rod problem, Figs. 6 and 7(d)–(f) show that the present higher

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

105

Fig. 10. Open knot NURBS discretization of the square membrane problem: (a) control net; (b) physical mesh.

Fig. 11. Comparison of the first four frequencies for the square membrane problem using various mass matrices with open knot quadratic basis functions: (a) the fundamental frequency; (b) the second frequency; (c) the third frequency; (d) the fourth frequency.

order mass matrices still produce the most accurate solutions with the highest convergence rates compared with those of other mass matrices, while the orders of accuracy are slightly below the theoretically derived optimal accuracy orders. The reason is that for the fixed–fixed rod problem, the first and last non-uniform basis functions associated with the two boundary points as shown in Fig. 1(b) are constrained in the computation, however for the fixed-free problem the non-uniform basis function associated with the free end participates into the rod vibration. As a whole, all the results uniformly demonstrate that the proposed higher order

mass matrices yield superior solution accuracy compared with those obtained by the consistent mass matrix.

5.2. Vibration of square membrane Consider a square membrane structure as shown in Fig. 8. The membrane is subjected to four side fixed boundary conditions. The analytical solution of the vibration frequency for this problem is [49]:

106

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

Fig. 12. NURBS discretization of the circular membrane problem: (a) control net; (b) physical mesh.

Fig. 13. Comparison of the first four frequencies for the circular membrane problem: (a) the fundamental frequency; (b) the second frequency; (c) the third frequency; (d) the fourth frequency.

qffiffiffiffiffiffiffiffiffiffiffiffiffi

xij ¼ pc=L i2 þ j2

ð88Þ

The geometric and material properties for this problem are: length L = 1, thickness s = 1, material density q = 1, Young’s modulus E = 1. The control net and mesh for a typical NURBS discretization with periodic quadratic basis functions are also plotted in Fig. 8. The first four numerical frequencies such as xh11 , xh12 , xh21 , and xh22 , are compared in Fig. 9(a)–(d) with their analytical counterparts given in Eq. (88). The results in Fig. 9(b) and (c) show that in case of i – j, say, xh12 , xh21 , although marginal improvement of

the convergence rates for the numerical frequencies is observed, the present higher mass matrix does produce much more accurate frequency solutions than it is consistent counterpart. As for the case of i = j in Fig. 9(a) and (d), for example, xh11 and xh22 , the vibration frequencies match very well with the theoretically derived 6th order of accuracy as elaborated in Eq. (85). Meanwhile, this problem is analyzed by using the open knot vector quadratic NURBS discretization as depicted in Fig. 10. The results are listed in Fig. 11(a)–(d), which are in harmonious agreement with those given by the periodic basis functions.

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108

5.3. Vibration of circular membrane As shown in Fig. 12, the last numerical example considered herein is the circular membrane vibration problem in order to demonstrate the robustness of the proposed method. The circular membrane is fixed at the surrounding boundary and its radius is r = 0.5. The material properties are the same as those of the previous square membrane problem. The analytical frequencies for this problem are the positive zeros of Bessel functions [49]. A typical NURBS discretization for the circular membrane [50] with open knot quadratic basis functions is also shown in Fig. 12. The proposed optimal parameters for the design of higher order mass matrix are directly adopted to construct the reduced bandwidth mass matrix and subsequently the higher order mass matrix based on the standard consistent mass matrix. The convergence results for the first four frequencies are listed in Fig. 13(a)–(d). The results clearly show that for this non-uniform discretization, the expected optimal order of accuracy may not be fully achieved for all the frequencies using the proposed algorithm, while the present method consistently performs superiorly compared with the standard consistent mass matrix. 6. Conclusions In the context of NURBS-based isogeometric method, novel quadratic and cubic higher order mass matrices were presented for structural vibration analysis. The higher order mass matrices were rationally constructed using a rational two-step method. The essential idea of the proposed method is to construct higher mass matrix through optimally combining the consistent mass matrix and another reduced bandwidth mass matrix with equal order of accuracy. This reduced bandwidth mass matrix is formulated by rephrasing the consistent mass matrix with smaller bandwidth and adjustable parameters. Then through the frequency error analysis, the adjustable parameters in the reduced bandwidth matrix can be decided by setting its order of accuracy the same as that of the consistent mass matrix. A mixed mass matrix follows through a linear combination of the consistent mass matrix and the reduced bandwidth matrix. The linear combination parameter is then obtained in order to get the optimal order of accuracy. It was shown that the quadratic and cubic higher order mass matrices have 6th and 8th orders of accuracy, which are two orders higher than those of their corresponding consistent mass formulations. A tensor product generalization of the one dimensional higher order mass matrices to two dimensional membrane vibration analysis was also described. Numerical results revealed the present higher order mass matrices are very effective and produce elevated order of accuracy for isogeometric structural vibration analysis. Acknowledgements The financial support of this work by the National Natural Science Foundation of China (11222221, 10972188) is gratefully acknowledged. References [1] 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 (2005) 4135–4195. [2] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons, Chichester, 2009. [3] J.A. Cottrell, A. Reali, Y. Bazilevs, T.J.R. Hughes, Isogeometric analysis of structural vibrations, Comput. Methods Appl. Mech. Engrg. 195 (2006) 5257– 5296. [4] J.A. Cottrell, T.J.R. Hughes, A. Reali, Studies of refinement and continuity in isogeometric structural analysis, Comput. Methods Appl. Mech. Engrg. 196 (2007) 4160–4183.

107

[5] A. Reali, An isogeometric analysis approach for the study of structural vibrations, J. Earthquake Engrg. 10 (2006) 1–30. [6] T.J.R. Hughes, A. Reali, G. Sangalli, Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: comparison of p-method finite elements with k-method NURBS, Comput. Methods Appl. Mech. Engrg. 197 (2008) 4104–4124. [7] S. Shojaee, E. Izadpanah, N. Valizadeh, J. Kiendl, Free vibration analysis of thin plates by using a NURBS-based isogeometric approach, Finite Elem. Anal. Des. 61 (2012) 23–34. [8] S. Shojaee, N. Valizadeh, E. Izadpanah, T. Bui, T.V. Vu, Free vibration and buckling analysis of laminated composite plates using the NURBS-based isogeometric finite element method, Compos. Struct. 94 (2012) 1677– 1693. [9] C.H. Thai, H. Nguyen-Xuan, N. Nguyen-Thanh, T.H. Le, T. Nguyen-Thoi, T. Rabczuk, Static, free vibration, and buckling analysis of laminated composite Reissner–Mindlin plates using NURBS-based isogeometric approach, Int. J. Numer. Methods Engrg. 91 (2012) 571–603. [10] Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali, G. Scovazzi, Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 197 (2007) 173–201. [11] Y. Zhang, Y. Bazilevs, S. Goswami, C. Bajaj, T.J.R. Hughes, Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2943–2959. [12] Y. Bazilevs, T.J.R. Hughes, NURBS-based isogeometric analysis for the computation of flows about rotating components, Comput. Mech. 43 (2008) 143–150. [13] Y. Bazilevs, I. Akkerman, Large eddy simulation of turbulent Taylor–Couette flow using isogeometric analysis and the residual-based variational multiscale method, J. Comput. Phys. 229 (2010) 3402–3414. [14] K. Chang, T.J.R. Hughes, V.M. Calo, Isogeometric variational multiscale largeeddy simulation of fully-developed turbulent flow over a wavy wall, Comput. Fluids 68 (2012) 94–104. [15] Y. Bazilevs, V.M. Calo, Y. Zhang, T.J.R. Hughes, Isogeometric fluid–structure interaction analysis with applications to arterial blood flow, Comput. Mech. 38 (2006) 310–322. [16] Y. Bazilevs, V.M. Calo, T.J.R. Hughes, Y. Zhang, Isogeometric fluid–structure interaction: theory, algorithms and computations, Comput. Mech. 43 (2008) 3–37. [17] W.A. Wall, M.A. Frenzel, C. Cyron, Full analytical sensitivities in NURBS based isogeometric shape optimization, Comput. Methods Appl. Mech. Engrg. 197 (2008) 2976–2988. [18] X. Qian, Isogeometric structural shape optimization, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2059–2071. [19] B. Koo, M. Yoon, S. Cho, Isogeometric shape design sensitivity analysis using transformed basis functions for Kronecker delta property, Comput. Methods Appl. Mech. Engrg. 253 (2013) 505–516. _ Temizera, P. Wriggers, T.J.R. Hughes, Contact treatment in isogeometric [20] I. analysis with NURBS, Comput. Methods Appl. Mech. Engrg. 200 (2011) 1100– 1112. [21] J. Lu, Isogeometric contact analysis: geometric basis and formulation for frictionless contact, Comput. Methods Appl. Mech. Engrg. 200 (2011) 726–741. [22] H. Gómez, V.M. Calo, Y. Bazilevs, T.J.R. Hughes, Isogeometric analysis of the Cahn–Hilliard phase-field model, Comput. Methods Appl. Mech. Engrg. 197 (2008) 4333–4352. [23] C.V. Verhoosel, M.A. Scott, T.J.R. Hughes, An isogeometric analysis approach to gradient damage models, Int. J. Numer. Methods Engrg. 86 (2011) 115–134. [24] P. Fischer, M. Klassen, J. Mergheim, P. Steinmann, R. Müller, Isogeometric analysis of 2D gradient elasticity, Comput. Mech. 47 (2011) 325–334. [25] J. Kiendl, K.U. Bletzinger, J. Linhard, R. Wuchner, Isogeometric shell analysis with Kirchhoff–Love elements, Comput. Methods Appl. Mech. Engrg. 198 (2009) 3902–3914. [26] D.J. Benson, Y. Bazilevs, M.C. Hsu, T.J.R. Hughes, Isogeometric shell analysis: the Reissner–Mindlin shell, Comput. Methods Appl. Mech. Engrg. 199 (2010) 276–289. [27] D.J. Benson, Y. Bazilevs, M.C. Hsu, T.J.R. Hughes, A large deformation, rotationfree isogeometric shell, Comput. Methods Appl. Mech. Engrg. 200 (2011) 1367–1378. [28] N. Nguyen-Thanh, J. Kiendl, H. Nguyen-Xuan, R. Wüchner, K.U. Bletzinger, Y. Bazilevs, T. Rabczuk, Rotation free isogeometric thin shell analysis using PHT-splines, Comput. Methods Appl. Mech. Engrg. 200 (2011) 3410– 3424. [29] R. Echter, B. Oesterle, M. Bischoff, A hierarchic family of isogeometric shell finite elements, Comput. Methods Appl. Mech. Engrg. 254 (2013) 170–180. [30] L. Beirao da Veiga, C. Lovadin, A. Reali, Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods, Comput. Methods Appl. Mech. Engrg. 241 (2012) 38–51. [31] 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 (2010) 229–263. [32] M.R. Döfel, B. Jüttler, B. Simeon, Adaptive isogeometric analysis by local hrefinement with T-Splines, Comput. Methods Appl. Mech. Engrg. 199 (2010) 264–275. [33] N. Nguyen-Thanh, H. Nguyen-Xuan, S.P.A. Bordas, T. Rabczuk, Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-

108

[34] [35] [36] [37]

[38]

[39]

[40]

[41]

D. Wang et al. / Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108 dimensional elastic solids, Comput. Methods Appl. Mech. Engrg. 200 (2011) 1892–1908. P. Wang, J. Xu, J. Deng, Adaptive isogeometric analysis using rational PHTsplines, Comput. Aided Des. 43 (2011) 1438–1448. R. Sevilla, S. Fernandez-Mendez, A. Huerta, NURBS-enhanced finite element method, Int. J. Numer. Methods Engrg. 76 (2008) 56–83. A. Shaw, D. Roy, NURBS-based parametric mesh-free methods, Comput. Methods Appl. Mech. Engrg. 197 (2008) 1541–1567. F. Auricchio, L. Beirao da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli, Isogeometric collocation methods, Math. Models Methods Appl. Sci. 20 (2010) 2075–2107. F. Auricchio, L. Beirao da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli, Isogeometric collocation for elastostatics and explicit dynamics, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 2–14. R. Simpson, S.P.A. Bordas, J. Trevelyan, T. Rabczuk, A two-dimensional isogeometric boundary element method for elastostatic analysis, Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100. M.A. Scott, R.N. Simpson, J.A. Evans, S. Lipton, S.P.A. Bordas, T.J.R. Hughes, T.W. Sederberg, Isogeometric boundary element analysis using unstructured Tsplines, Comput. Methods Appl. Mech. Engrg. 254 (2013) 197–221. D.J. Benson, Y. Bazilevs, E. De Luycker, M.C. Hsu, M.A. Scott, T.J.R. Hughes, T. Belytschko, A generalized finite element formulation for arbitrary basis

[42]

[43] [44] [45] [46] [47] [48]

[49] [50]

functions: from isogeometric analysis to XFEM, Int. J. Numer. Methods Engrg. 83 (2010) 765–785. S.S. Ghorashi1, N. Valizadeh, S. Mohammadi, Extended isogeometric analysis for simulation of stationary and propagating cracks, in: Int. J. Numer. Methods Engrg. 89 (2012) 1069–1101. C. Stavrinidis, J. Clinckemaillie, J. Dubois, New concepts for finite element mass matrix formulations, AIAA J. 27 (1989) 1249–1255. K. Kim, A review of mass matrices for eigenproblems, Comput. Struct. 46 (1993) 1041–1048. I. Fried, M. Chavez, Super accurate finite element eigenvalue computation, J. Sound Vib. 275 (2004) 415–422. T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, Mineola, New York, 2000. D.F. Rogers, An Introduction to NURBS with Historical Perspective, Academic Press, London, 2001. D. Wang, J. Xuan, An improved NURBS-based isogeometric analysis with enhanced treatment of essential boundary conditions, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2425–2436. N.H. Asmar, Partial Differential Equations with Fourier Series and Boundary Value Problems, second ed., Prentice Hall, New Jersey, 2004. A.V. Vuong, Ch. Heinrich, B. Simeon, ISOGAT: a 2D tutorial MATLAB code for isogeometric analysis, Comput. Aided Geom. Des. 27 (2010) 644–655.