Optimal reduction of numerical dispersion for wave propagation problems. Part 2: Application to 2-D isogeometric elements

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268 www.elsevier.com/locate/cma

Optimal reduction of numerical dispersion for wave propagation problems. Part 2: Application to 2-D isogeometric elements A. Idesman ∗, B. Dey Department of Mechanical Engineering, Texas Tech University, Lubbock, TX 79409-1021, USA Available online 19 April 2017

Abstract Based on the optimal coefficients of the stencil equation, a numerical technique for the reduction of the numerical dispersion error has been suggested. New isogeometric elements with the reduced numerical dispersion error for wave propagation problems in the 2-D case have been developed with the suggested approach. By the minimization of the order of the dispersion error of the stencil equation, the order of the dispersion error is improved from order 2 p (the conventional isogeometric elements) to order 4 p (the isogeometric elements with reduced dispersion) where p is the order of the polynomial approximations. Because all coefficients of the stencil equation are obtained from the minimization procedure, the obtained accuracy is maximum possible. The corresponding elemental mass and stiffness matrices of the isogeometric elements with reduced dispersion are calculated with help of the optimal coefficients of the stencil equation. The analysis of the dispersion error of the isogeometric elements with the lumped mass matrix has also shown that independent of the procedures for the calculation of the lumped mass matrix, the second order of the dispersion error cannot be improved with the conventional stiffness matrix. However, the dispersion error with the lumped mass matrix can be improved from the second order to order 2 p by the modification of the stiffness matrix. The numerical examples confirm the computational efficiency of the isogeometric elements with reduced dispersion. The numerical results obtained by the new and conventional isogeometric elements may include spurious oscillations due to the dispersion error. These oscillations can be quantified and filtered by the two-stage time-integration technique developed recently. The approach developed in the paper can be directly applied to other space-discretization techniques with similar stencil equations. c 2017 Elsevier B.V. All rights reserved. ⃝

Keywords: High-order elements with reduced dispersion; Isogeometric elements; Wave propagation; Numerical dispersion

1. Introduction Part 2 of the paper is the extension of the approach for the reduction of the numerical dispersion error for wave propagation problems to the 2-D case with the application to the isogeometric elements (see also Part 1 [1] for the quadratic and cubic isogeometric elements with reduced dispersion in the 1-D case). Wave propagation in an isotropic ∗ Corresponding author. Fax: +1 806 742 3540.

E-mail address: [email protected] (A. Idesman). http://dx.doi.org/10.1016/j.cma.2017.04.008 c 2017 Elsevier B.V. All rights reserved. 0045-7825/⃝

236

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

homogeneous medium is described by the following scalar wave equation in domain Ω : ∂ 2u − c2 ∇ 2 u = 0, (1) ∂t 2 with the boundary conditions n · ▽u = g1 on Γ t and u = g2 on Γ u , and the initial conditions u(x, t = 0) = g3 , v(x, t = 0) = g4 in Ω . Here, u is the field variable, v = u˙ is the velocity, c is the wave velocity, t is the time, Γ t and Γ u denote the natural and essential boundaries, gi (i = 1, 2, 3, 4) are the given functions, n is the outward unit normal on Γ t . The application of the continuous Galerkin approach and the space discretization (e.g., the finite elements, spectral elements, isogeometric elements; see [2,3,4] and others) to Eq. (1) leads to a system of ordinary differential equations in time ¨ + c2 K U = 0, MU

(2)

with M=

∑

Me ,

K=

e

∑

Ke ,

(3)

e

where U(t) is the vector of the field variable, the global mass M and stiffness K matrices have a banded structure and are obtained by the summation of the corresponding local (element Ω e ) matrices Me and K e : ∫ Me = N T NdΩ e , (4) e Ω ] [ ] ∫ [ ∂N T ∂N e K = dΩ e . (5) ∂x ∂x Ωe are the shape matrix and its derivative with respect to the physical coordinate x; see [2,3,5]. Due to Here, N and ∂N ∂x the space discretization, the exact solution to Eq. (2) contains the numerical dispersion error; e.g., see [6,7]. A brief review of the numerical techniques used for the analysis and the reduction of the numerical dispersion error is reported in Part 1 of the paper; see [1]. Here, it is necessary to mention that we have not seen in the literature the numerical techniques that reduce the order of the numerical dispersion error for high-order elements in the multidimensional case (e.g., the numerical approaches presented in [8,9,10] improve the order of the dispersion error for high-order finite and isogeometric elements from order 2 p to order 2 p + 2 in the 1-D case only and for some specific directions in the 2-D case). In Part 1 of the paper (see [1]) the new isogeometric elements with reduced dispersion have been developed in the 1-D case. They reduce the dispersion error from order 2 p to order 4 p. In Part 2 of the paper the same improvement in the order of the dispersion error has been developed in the 2-D case. In contrast to [8,9,10], the improved order of the dispersion error is much higher and for the first time this is valid in the general 2-D case (the improved order of the dispersion error is independent of the propagation direction of harmonic waves). Based on the dispersion analysis, the quadratic isogeometric elements with reduced dispersion with the non-diagonal and diagonal mass matrices (the diagonal mass matrix can be used with explicit time-integration methods) have been developed in Section 2. Similar results have been obtained for the cubic isogeometric elements in Section 3. The order of the dispersion error has been improved from the order 2 p to the order 4 p for the new isogeometric elements with the non-diagonal mass matrix and from the second order to the order 2 p for the new isogeometric elements with the diagonal mass matrix. The numerical examples in Section 4 show the computational efficiency of the proposed approach. For the derivation of many analytical expressions presented below the computational program “Mathematica” has been used. The approach developed in the paper, can be directly applied to the reduction of the dispersion error of other space-discretization methods that have a similar structure of the stencil equation (e.g., to the high-order finite difference method). 2. Dispersion analysis of quadratic isogeometric elements in the 2-D case Assuming time-harmonic solutions of Eq. (1) in the form u(x, t) = exp(i ωt)u(x),

(6)

the wave Eq. (1) reduces to the Helmholtz equation ∇ 2 u + k 2 u = 0,

(7)

237

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

where ω is the angular velocity, x is the position vector, k = ω/c is the wave number, i = following exact solutions:

√

u(x) = exp(±ikn · x),

−1. Eq. (7) allows the (8)

where n is the unit normal to the wave front. After the space discretization, Eq. (7) reduces to (K − k 2 M ) U = 0.

(9)

Uniform meshes are used for the dispersion analysis considered below. The calculation of the conventional mass and stiffness matrices of the isogeometric elements is based on the univariate B-spline basis functions; see [3,5]. In the 1-D case they are defined recursively starting with p = 0 { 1, if ξi ≤ ξ < ξi+1 Ni, p (ξ ) = (10) 0, otherwise. For p ≥ 1: Ni, p (ξ ) =

ξi+ p+1 − ξ ξ − ξi Ni, p−1 (ξ ) + Ni+1, p−1 (ξ ), ξi+ p − ξi ξi+ p+1 − ξi+1

(11)

where a knot vector {ξ1 = 0, . . . , ξi , . . . , ξn+ p+1 = 1} is a set of non-decreasing real numbers representing coordinates in the parametric space of the curve, p is the order of the B-spline, n is the number of the basis functions, i = 1, 2, . . . , n + p + 1. In the 2-D case the basis functions can be constructed with the help of the tensor product as follows: pq

Mi j (ξ, η) = Ni, p (ξ ) N j,q (η),

(12)

where Ni, p (ξ ) and N j,q (η) are the basis functions of order p and q given by Eq. (11). The 2-D basis functions of the same order in the ξ and η directions ( p = q) are used in the following sections. 2.1. Conventional quadratic isogeometric elements Let us first start with the dispersion analysis of the conventional quadratic isogeometric elements based on the continuous Galerkin formulation; e.g., see [3,5]. The mass and stiffness matrices of a typical interior element in Eqs. (4) and (5) are (see also [9]) ⎛ ⎞ 36 78 6 78 169 13 6 13 1 ⎜ 78 324 78 169 702 169 13 54 13 ⎟ ⎜ ⎟ ⎜ 6 78 36 13 169 78 1 13 6 ⎟ ⎜ ⎟ ⎜ 78 169 13 324 702 54 78 169 13 ⎟ 2 ⎜ ⎟ h ⎜169 702 169 702 2916 702 169 702 169⎟ , Me = (13) ⎜ ⎟ 14400 ⎜ 13 169 78 ⎟ 54 702 324 13 169 78 ⎜ ⎟ ⎜ 6 13 1 78 169 13 36 78 6 ⎟ ⎜ ⎟ ⎝ 13 54 13 169 702 169 78 324 78 ⎠ 1 13 6 13 169 78 6 78 36 ⎛

12 ⎜ 10 ⎜ ⎜ −2 ⎜ ⎜ 10 1 ⎜ e ⎜−13 K = 360 ⎜ ⎜ −7 ⎜ ⎜ −2 ⎜ ⎝ −7 −1

10 −2 10 −13 −7 −2 −7 60 10 −13 −14 −13 −7 −26 10 12 −7 −13 10 −1 −7 −13 −7 60 −14 −26 10 −13 −14 −13 −14 108 −14 −13 −14 −13 10 −26 −14 60 −7 −13 −7 −1 10 −13 −7 12 10 −26 −7 −13 −14 −13 10 60 −7 −2 −7 −13 10 −2 10

⎞ −1 −7 ⎟ ⎟ −2 ⎟ ⎟ −7 ⎟ ⎟ −13⎟ ⎟, 10 ⎟ ⎟ −2 ⎟ ⎟ 10 ⎠ 12

(14)

where the sequence of h−spaced control points in the x and y directions with x A = h A and y B = h B is used. Considering the discretized Helmholtz equation (9) on an infinite plane with the sequence of h− spaced control

238

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

points in the x and y directions, the stencil equation for the degree of freedom u A,B can be calculated with the help of Eqs. (3), (9), (13) and (14) for 3 × 3 = 9 neighboring elements. It includes 25 degrees of freedom u i, j (i = A − 2, A − 1, A, A + 1, A + 2 and j = B − 2, B − 1, B, B + 1, B + 2) which are close to the degree of freedom u A,B (see Fig. A.24 in the Appendix A) and has the following form: k 2 h 2 [4356u A,B + 676(u (A−1),(B+1) + u (A+1),(B+1) + u (A−1),(B−1) + u (A+1),(B−1) ) + 1716(u (A−1),B + u (A+1),B + u A,(B+1) + u A,(B−1) ) + (u (A−2),(B−2) + u (A−2),(B+2) + u (A+2),(B+2) + u (A+2),(B−2) ) + 26(u (A+2),(B−1) + u (A−2),(B+1) + u (A−1),(B−2) + u (A+2),(B+1) + u (A−1),(B+2) + u (A+1),(B+2) + u (A+1),(B−2) + u (A−2),(B−1) ) + 66(u (A−2),B + u (A+2),B + u A,(B+2) + u A,(B−2) )] + 40[396u A,B − 52(u (A−1),(B+1) + u (A+1),(B+1) + u (A−1),(B−1) + u (A+1),(B−1) ) + 12(u (A−1),B + u (A+1),B + u A,(B+1) + u A,(B−1) ) − (u (A−2),(B−2) + u (A−2),(B+2) + u (A+2),(B+2) + u (A+2),(B−2) ) − 14(u (A+2),(B−1) + u (A−2),(B+1) + u (A−1),(B−2) + u (A+2),(B+1) + u (A−1),(B+2) + u (A+1),(B+2) + u (A+1),(B−2) + u (A−2),(B−1) ) − 30(u (A−2),B + u (A+2),B + u A,(B+2) + u A,(B−2) )] = 0.

(15)

Eq. (15) allows the following solutions (similar to Eq. (8)): u A,B = exp(±ikh (n 1 h A + n 2 h B)),

(16)

where kh is the numerical wave number, n 1 and n 2 are the components of the unit vector (n 21 + n 22 = 1). Inserting Eq. (16) into Eq. (15) the following relation between the exact k and numerical kh wave numbers can be found: √ √ 2 5 a k = (17) √ , kh (kh h) b a = −14 cos((kh h)(n 1 − 2n 2 )) − 52 cos((kh h)(n 1 − n 2 )) − cos(2(kh h)(n 1 − n 2 )) − 14 cos((kh h)(2n 1 − n 2 )) − 52 cos((kh h)(n 1 + n 2 )) − cos(2(kh h)(n 1 + n 2 )) − 14 cos((kh h)(2n 1 + n 2 )) − 14 cos((kh h)(n 1 + 2n 2 )) + 12 cos(n 1 (kh h)) − 30 cos(2n 1 (kh h)) + 12 cos(n 2 (kh h)) − 30 cos(2n 2 (kh h)) + 198, b = 676 cos(n 1 (kh h)) cos(n 2 (kh h)) + 26 cos(n 1 (kh h)) cos(2n 2 (kh h)) + 26 cos(2n 1 (kh h)) cos(n 2 (kh h)) + cos(2n 1 (kh h)) cos(2n 2 (kh h)) + 858 cos(n 1 (kh h)) + 33 cos(2n 1 (kh h)) + 858 cos(n 2 (kh h)) + 33 cos(2n 2 (kh h)) + 1089. √ Expanding the right-hand side of Eq. (17) into a Taylor series at small h ≪ 1 and using n 1 = 1 − n 22 , we get ( 4 ) ( 8 ) ( ) 3n 2 − 3n 2 2 + 1 (kh h)4 2n 2 − 4n 2 6 + 6n 2 4 − 4n 2 2 + 1 (kh h)6 k =1+ + + O (kh h)8 , (18) kh 1440 6720 i.e., in the 2-D case the conventional quadratic isogeometric elements yield the fourth order of the dispersion error for (k/kh -1). For convenience, the components of the unit normal in Eq. (17) can be expressed as n 1 = sin θ and n 2 = cos θ where θ is the propagation angle of harmonic waves. The numerical dispersion error e = khk−k calculated with the help of Eq. (17) is plotted in Figs. 1(a) and 2(a)–(c) in the logarithmic scale. As can be seen, the dispersion error decreases with the decrease in the mesh size kh h. However, the anisotropy of the numerical dispersion error (its dependence on angle θ of the propagating harmonic waves) is significant at any mesh size; see Figs. 2(a)– (c). The maximum dispersion error corresponds to the propagation of harmonic waves along the Cartesian axes with θ = 0, π/2, π, 3π/2. The minimum dispersion error corresponds to the propagation of harmonic waves with θ = π/4, 3π/4, 5π/4, 7π/4; see Figs. 2(a)–(c). 2.2. Quadratic isogeometric elements with reduced dispersion For the derivation of the quadratic isogeometric elements with reduced dispersion let us assume that the stencil equation for the degree of freedom u A,B includes the same degrees of freedom as that for the conventional elements

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

239

Fig. 1. The numerical dispersion error in the logarithmic scale Log10 |e| (with e = khk−k ) as a function of the mesh size kh h and angle θ . (a) and (b) correspond to the conventional and new quadratic isogeometric elements with the non-diagonal mass matrix.

in Eq. (15) with unknown coefficients m j and k j ( j = 1, 2, 3, 4, 5, 6) before each degree of freedom. It can be written as follows: k 2 h 2 {m 1 u A,B + m 2 [(u (A−1),(B+1) + u (A+1),(B−1) ) + (u (A+1),(B+1) + u (A−1),(B−1) )] + m 3 [(u (A−1),B + u (A+1),B ) + (u A,(B+1) + u A,(B−1) )] + m 4 [(u (A−2),(B−2) + u (A+2),(B+2) ) + (u (A−2),(B+2) + u (A+2),(B−2) )] + m 5 [(u (A+2),(B−1) + u (A−2),(B+1) ) + (u (A−1),(B−2) + u (A+1),(B+2) ) + (u (A−1),(B+2) + u (A+1),(B−2) ) + (u (A−2),(B−1) + u (A+2),(B+1) )] + m 6 [(u (A−2),B + u (A+2),B ) + (u A,(B+2) + u A,(B−2) )]} − {k1 u A,B + k2 [(u (A−1),(B+1) + u (A+1),(B−1) ) + (u (A+1),(B+1) + u (A−1),(B−1) )] + k3 [(u (A−1),B + u (A+1),B ) + (u A,(B+1) + u A,(B−1) )]

240

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

a

d –1

–2 –3

–2

–4 –5

–3

–6

–4

–7 1

2

3

4

5

6

b

1

2

3

4

5

6

1

2

3

4

5

6

1

2

3

4

5

6

e –2

–1

–3 –2

–4 –5

–3

–6 –4

–7 1

2

3

4

5

6

c

f –2

–1

–3 –2

–4 –5

–3

–6 –4

–7 1

2

3

4

5

6

Fig. 2. The numerical dispersion error in the logarithmic scale Log10 |e| (with e = khk−k ) as a function of the angle θ. (a, b, c) and (d, e, f) correspond to the conventional and new quadratic isogeometric elements with the non-diagonal mass matrix at kh h = 2 (a, d), kh h = 1 (b, e) and kh h = 0.5 (c, f).

+ k4 [(u (A−2),(B−2) + u (A+2),(B+2) ) + (u (A−2),(B+2) + u (A+2),(B−2) )] + k5 [(u (A+2),(B−1) + u (A−2),(B+1) ) + (u (A−1),(B−2) + u (A+1),(B+2) ) + (u (A−1),(B+2) + u (A+1),(B−2) ) + (u (A−2),(B−1) + u (A+2),(B+1) )] + k6 [(u (A−2),B + u (A+2),B ) + (u A,(B+2) + u A,(B−2) )]} = 0,

(19)

where the unknown coefficients m j and k j ( j = 1, 2, 3, 4, 5, 6) can be expressed in terms of the elemental mass and stiffness matrices by Eqs. (3) and (9) (at this point, these elemental matrices are not defined). For example, comparing Eqs. (15) and (19), it can be found that m 1 = 4356, m 2 = 676, m 3 = 1716, m 4 = 1, m 5 = 26, m 6 = 66, k1 = 15 840, k2 = −2080, k3 = 480, k4 = −40, k5 = −560, k6 = −1200 for the conventional quadratic isogeometric elements. Similar to the conventional isogeometric elements, the symmetry of the coefficients m j and k j in the stencil equation (19) is assumed for the degrees of freedom u j,l ( j = A − 2, A − 1, A, A + 1, A + 2 and l = B − 2, B − 1, B, B + 1, B + 2) symmetrically located with respect to the degree of freedom u A,B ; see Fig. A.24 in Appendix A. Similar to the derivations in the 1-D case considered in Part 1 of the paper (see [1]), by the insertion of

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

241

Eq. (16) into Eq. (19) the relation between the exact k and numerical kh wave numbers can be found. However, in the multidimensional case the analytical derivations of the optimal coefficients m j and k j ( j = 1, 2, 3, 4, 5, 6) in Eq. (19) can be significantly simplified by the insertion of Eq. (16) with kh = k into Eq. (19) and the reduction of the residual of the obtained expression (if the dispersion error is zero then the residual is exactly equal to zero); e.g., see [11,12] for the reduction of the dispersion error of the linear finite elements. It is necessary to note that the direct reduction of the dispersion error or the reduction of the residual of the dispersion equation yield exactly the same results. After the insertion of Eq. (16) with kh = k into Eq. (19), Eq. (19) can be reduced to the following form: R = R˜ exp(±ik(n 1 h A + n 2 h B)),

(20)

R˜ = k 2 h 2 {m 1 + 2m 2 [cos[k(n 1 h − n 2 h)] + cos[k(n 1 h + n 2 h)]] + 2m 3 [cos[kn 1 h] + cos[kn 2 h]] + 2m 4 [cos[2k(n 1 h + n 2 h)] + cos[2k(n 1 h − n 2 h)]] + 2m 5 [cos[k(2n 1 h − n 2 h)] + cos[k(n 1 h + 2n 2 h)] + cos[k(n 1 h − 2n 2 h)] + cos[k(2n 1 h + n 2 h)]] + 2m 6 [cos[2kn 1 h] + cos[2kn 2 h]]} − {k1 + 2k2 [cos[k(n 1 h − n 2 h)] + cos[k(n 1 h + n 2 h)]] + 2k3 [cos[kn 1 h] + cos[kn 2 h]] + 2k4 [cos[2k(n 1 h + n 2 h)] + cos[2k(n 1 h − n 2 h)]] + 2k5 [cos[k(2n 1 h − n 2 h)] + cos[k(n 1 h + 2n 2 h)] + cos[k(n 1 h − 2n 2 h)] + cos[k(2n 1 h + n 2 h)]] + 2k6 [cos[2kn 1 h] + cos[2kn 2 h]]} = 0.

(21)

For the simplification of the expressions in each parenthesis in Eq. (19), the following trigonometric formula is used: exp(ia) + exp(ib) = 2 exp(i(a + b)/2) cos[(a − b)/2]. The residual R of the dispersion equation can be decreased by ˜ see Eq. (20). Expanding the cosine functions in the expression for R˜ (see Eq. (21)) into a Taylor the reduction of R; series at small h ≪ 1, combining the terms with the same order of h, and equating zero the first five terms with the lowest orders of h, the coefficients m j and k j ( j = 1, 2, 3, 4, 5, 6) can be found from a system of the corresponding algebraic equations. They are: m 1 = (1650a1 − 3252a2 − 259a3 )/275, m 3 = (−440a1 − 57a2 − 184a3 )/110,

m 2 = (1320a1 − 783a2 + 124a3 )/495, m 4 = (330a1 − 54a2 + 37a3 )/1980, m 6 = a1 ,

m 5 = (−6600a1 − 351a2 − 392a3 )/9900, k1 = −3(213a2 + 5a3 )/11, k4 = (59a2 − 17a3 )/220,

k2 = −16(13a2 − 14a3 )/55, k 5 = a2 ,

k3 = 2(441a2 − 128a3 )/55,

k 6 = a3 ,

(22) √ where a1 , a2 and a3 are three arbitrary coefficients. At the derivation of Eq. (22) the expression n 1 = 1 − n 22 for the components of the unit vector is used where −1 ≤ n 2 ≤ 1 is an arbitrary number. Expanding the right-hand side of the dispersion equation (similar to the derivation of Eq. (18)) into a Taylor series at small h ≪ 1 with the coefficients given by Eq. (22) we get ( ) k d(kh h)8 =1+ + O (kh h)10 , kh 19051200(3a2 + a3 )

(23)

d = (6n 2 8 (38500a1 − 123a2 + 1884a3 ) − 12n 2 6 (38500a1 − 123a2 + 1884a3 ) + 55n 2 4 (4200a1 + 93a2 + 136a3 ) + n 2 2 (3824a3 − 5853a2 ) − 158(3a2 + a3 )), i.e., in the 2-D case the 8th order of the dispersion error for (k/kh -1) can be obtained by the suggested approach. It is necessary to note that for the harmonic waves propagating along the Cartesian axes (n 2 = 0 for the x-axis and n 2 = 1 for the y-axis), Eq. (23) reduces to ( ) k 79(kh h)8 2633(kh h)10 =1− − + O (kh h)12 , kh 9525600 2640496320

(24)

242

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

i.e., it is exactly the same as that for the quadratic isogeometric elements with reduced dispersion in the 1-D case (see [1]). It is interesting to note that the order of the dispersion error in Eq. (23) cannot be improved with the use of the three arbitrary coefficients a j ( j = 1, 2, 3) because for the particular directions of harmonic waves with n 2 = 0 or n 2 = 1, the leading term of the dispersion error is independent of the coefficients a j ; see Eq. (24). The leading term of the dispersion error in Eq. (23) is the polynomial function of n 2 with some undefined coefficients. Additional constraints can be imposed for the coefficients of this polynomial function. For example, it can be imposed that all coefficients containing any power of the term n 2 are zero. In this case Eqs. (22) and (23) can be rewritten as follows: m 1 = a1 ,

m 2 = 79708a1 /785673,

m 3 = 79708a1 /785673,

m 4 = 2437a1 /3142692,

m 5 = 2768a1 /785673,

k1 = 127845a1 /29099,

k2 = −46880a1 /261891,

k4 = −11465a1 /1047564,

k5 = −19120a1 /261891,

m 6 = 457a1 /87297, k3 = −56800a1 /87297, k6 = −9755a1 /87297,

(25)

and ( ) k 79(kh h)8 + d(kh h)10 + O (kh h)12 , =1− kh 9525600 d=

(26)

3424n 2 12 − 10272n 2 10 + 8300n 2 8 + 520n 2 6 − 9080n 2 4 + 7108n 2 2 − 2633 , 2640496320

i.e., in this case the leading term of the dispersion error is independent of n 2 and this reduces the numerical anisotropy. a1 in Eq. (25) is an arbitrary coefficient that can be easily found from the condition of the preservation of the kinetic energy (similar to Eq. (33)). Similar to the conventional isogeometric elements in Section 2.1, the numerical dispersion error e = khk−k is calculated for the new isogeometric elements with the coefficients given by Eq. (25) and is plotted in Figs. 1(b) and 2(d)–(f) in the logarithmic scale. As can be seen, in contrast to the conventional elements, the dispersion error for the new elements is much smaller at the same mesh size kh h and decreases much faster with the decrease in the mesh size kh h; see Fig. 2. Moreover, for the new elements the anisotropy of the numerical dispersion error (its dependence on angle θ of the propagating harmonic waves) becomes smaller with the decrease in the mesh size kh h; see Figs. 2(d)–(f) (because the leading term in Eq. (26) is independent of angle θ ). Remark 1. Along with the time-harmonic solutions of Eq. (1) given by Eq. (6), a time-independent linear function u(x, t) = co + c1 x + c2 y also meets Eq. (1). Let us show that the vector U calculated with the help of this function, also meet the semi-discrete system (2). In this case it is necessary to check that K U = 0. First, it is easy to check by inspection that the constant solution u j,l = c0 ( j = A − 2, A − 1, A, A + 1, A + 2 and l = B − 2, B − 1, B, B + 1, B + 2) with the coefficients k j ( j = 1, 2, . . . , 6) given by Eq. (22) meets the stencil equation (19) with m j = 0 ( j = 1, 2, . . . , 6) (for time-independent solutions there is no contribution to the stencil equation from the mass matrix). Finally, by the calculation of u j,l = c0 +c1 x j +c2 yl ( j = A−2, A−1, A, A+1, A+2 and l = B −2, B −1, B, B +1, B +2) and their insertion into the stencil equation (19) with m j = 0 ( j = 1, 2, . . . , 6), the result will be reduced to the previous case with u j,l = c0 due to the structure of Eq. (19). Remark 2. The new approach improves the order of the local truncation error in space by four orders (similar to the improvement of the order of the numerical dispersion error for the new quadratic isogeometric elements); see the Appendix B. Let us find the elemental mass and stiffness matrices that yield the stencil equation (19) with the coefficients m j and k j ( j = 1, 2, . . . , 6) given by Eq. (22) (it is assumed that the stencil equation for the new elements is calculated similar to that for the conventional isogeometric elements according to Eqs. (3) and (9)). The new elemental mass and stiffness matrices Me and K e are assumed to have the same form as that for the conventional elemental matrices with

243

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

unknown coefficients d j and f j ⎛ d1 d4 d2 ⎜d4 d7 d4 ⎜ ⎜d2 d4 d1 ⎜ ⎜d4 d8 d5 ⎜ e 2⎜ M = h ⎜d6 d10 d6 ⎜d5 d8 d4 ⎜ ⎜d2 d5 d3 ⎜ ⎝d5 d9 d5 d3 d5 d2 ⎛

f1 ⎜ f4 ⎜ ⎜ f2 ⎜ ⎜ f4 ⎜ e K =⎜ ⎜ f6 ⎜ f5 ⎜ ⎜ f2 ⎜ ⎝ f5 f3

f4 f7 f4 f8 f 10 f8 f5 f9 f5

f2 f4 f1 f5 f6 f4 f3 f5 f2

for all degrees of freedom: d4 d8 d5 d7 d10 d9 d4 d8 d5

f4 f8 f5 f7 f 10 f9 f4 f8 f5

d6 d10 d6 d10 d11 d10 d6 d10 d6 f6 f 10 f6 f 10 f 11 f 10 f6 f 10 f6

d5 d8 d4 d9 d10 d7 d5 d8 d4 f5 f8 f4 f9 f 10 f7 f5 f8 f4

d2 d5 d3 d4 d6 d5 d1 d4 d2 f2 f5 f3 f4 f6 f5 f1 f4 f2

⎞ d3 d5 ⎟ ⎟ d2 ⎟ ⎟ d5 ⎟ ⎟ d6 ⎟ ⎟, d4 ⎟ ⎟ d2 ⎟ ⎟ d4 ⎠ d1

d5 d9 d5 d8 d10 d8 d4 d7 d4 f5 f9 f5 f8 f 10 f8 f4 f7 f4

(27)

⎞ f3 f5⎟ ⎟ f2⎟ ⎟ f5⎟ ⎟ f6⎟ ⎟, f4⎟ ⎟ f2⎟ ⎟ f4⎠ f1

(28)

where the symmetry of the coefficients of the mass and stiffness matrices for the degrees of freedom contributing to these matrices is used; see Appendix A. In this case the matrices Me and K e depend on the 22 unknown terms d j and f j ( j = 1, 2, . . . , 11). Considering 3 × 3 = 9 elements that contribute to the stencil equation for the degree of freedom u A,B according to Eqs. (3) and (9), the coefficients of the stencil equation (19) can be expressed in terms of the coefficients of the matrices Me and K e as follows: m 1 = 4d1 + d11 + 4d7 , m 4 = d3 ,

k1 = 4 f 1 + f 11 + 4 f 7 , k4 = f 3 ,

m 2 = 2d6 + 2d8 ,

m 5 = 2d5 ,

m 3 = 2d10 + 4d4 ,

m 6 = 2d2 + d9 ,

(29)

k2 = 2 f 6 + 2 f 8 ,

k5 = 2 f 5 ,

k3 = 2 f 10 + 4 f 4 ,

k6 = 2 f 2 + f 9 .

(30)

Solving simultaneously a system of linear algebraic equations (29), (30) and (22), the coefficients of the mass and stiffness matrices in Eqs. (27) and (28) for the 2-D quadratic isogeometric elements with reduced dispersion can be found. They are: d1 = 271/1890 + 3a1 /2 + 3a2 /4 − a3 − a7 /4, d2 = (a1 − a5 )/2,

d3 = (1 + 126a1 + 21a2 )/756,

d4 = 19/3024 − a1 − 3a2 /8 − a6 /2, d5 = 13/15120 − a1 /3 − a2 /72, d6 = 29/756 + 4a1 /3 + 7a2 /18 − a4 , d7 = a3 ,

d8 = a 4 ,

d9 = a5 ,

d10 = a6 ,

d11 = a7 ,

(31)

and f 1 = 355/504 − a10 /4 + 9a2 /2 − a8 , f 2 = −124/567 − 13a2 /6 + a8 /2 + a9 ,

f 3 = −59/4536 − a2 /6,

f 4 = −337/2268 + a10 /8 − 7a2 /6 − a9 /2,

f 5 = −55/2268 − a2 /6,

f 6 = 52/567 + 8a2 /3 − a9 , f 7 = a8 ,

f 8 = a9 ,

f 9 = 248/567 + 16a2 /3 − a8 − 2a9 ,

f 10 = −52/567 − a10 /4 − 8a2 /3 + a9 ,

f 11 = a10 ,

(32)

244

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

where a j (j=1,2, . . . ,10) are ten arbitrary coefficients that do not affect the stencil equation (19). For example, these coefficients can be used for the minimization of the error in the case of the combination of the different order elements (this will be studied in the future). The dispersion error is independent of the values of these coefficients. At the derivation of Eqs. (31) and (32) the following additional algebraic equations (see Eqs. (33) and (34)) have been also used. Let us assume that the velocity v(x, t) = u(x, ˙ t) = v0 is the same for the entire domain. In this case the kinetic energy for one square element can be calculated with the help of the mass matrix and by the analytical formula for the considered square element. Equating these two expressions we get V 0T Me V 0 = v02 h 2 ,

(33)

where V 0 = (v0 , v0 , v0 , v0 , v0 , v0 , v0 , v0 , v0 )T is the vector with nine equal components. Next, let us assume that the field function u(x, t) = u 0 is the same for the entire domain. In this case the internal forces are zero. These forces for one element can be calculated with the help of the stiffness matrix: K e U0 = 0,

(34)

where U0 = (u 0 , u 0 , u 0 , u 0 , u 0 , u 0 , u 0 , u 0 , u 0 )T is the vector with nine equal components. Remark. According to the procedure for the calculation of the optimal coefficients m j and k j (see Eq. (22)) of the stencil equation (19), the 8th order of the dispersion error in Eq. (23) is the maximum order for all quadratic isogeometric elements in the 2-D case. 2.3. Quadratic isogeometric elements with lumped mass matrix and reduced dispersion In Part 1 of the paper (see [1]), it has been shown in the 1-D case that for the conventional stiffness matrix and the lumped mass matrix, the dispersion error cannot exceed the second order of accuracy. For the harmonic waves propagating in the 2-D case along the Cartesian axes, the dispersion errors for the conventional stiffness matrix and the lumped mass matrix are the same as those in the 1-D case; i.e., the dispersion error cannot exceed the second order of accuracy in the 2-D case as well. However, similar to the 1-D case (see [1]) the accuracy of the results with the lumped mass matrix can be improved by the modification of the stiffness matrix. Let us consider the stencil equation (19) with the diagonal mass matrix (m 2 = m 3 = m 4 = m 5 = m 6 = 0) and the undefined coefficients k j ( j = 1, 2, 3, 4, 5, 6). For finding the optimal coefficients m 1 and k j ( j = 1, 2, 3, 4, 5, 6), the approach described in the previous Section 2.2 is used. Expanding the cosine functions in the expression for R˜ (see Eq. (21)) into a Taylor series at small h ≪ 1, combining the terms with the same order of h, and equating zero the first three terms with the lowest orders of h, the coefficients m 1 and k j ( j = 1, 2, 3, 4, 5, 6) can be found from a system of the corresponding algebraic equations. They are: m 1 = a1 , k1 = (72a2 + 60a3 + 5a1 )/2, k2 = 2(12a2 + 12a3 − a1 )/3, k3 = −2(9a2 + 8a3 ), k4 = −a2 − a3 /2 + a1 /24, k 5 = a2 , k 6 = a3 ,

(35) √

where a1 , a2 , a3 are three arbitrary coefficients. At the derivation of Eq. (35) the expression n 1 = 1 − n 22 for the components of the unit vector has been used where −1 ≤ n 2 ≤ 1 is an arbitrary number. Expanding the right-hand side of the dispersion equation (similar to the derivation of Eq. (18)) into a Taylor series at small h ≪ 1 with the coefficients given by Eq. (35) we get ( ) k (kh h)4 (a1 (12n 2 4 − 12n 2 2 − 1) − 90n 2 2 (n 2 2 − 1)(3a2 + 2a3 )) =1+ + O (kh h)6 , (36) kh 180a1 i.e., for the lumped mass matrix in the 2-D case, the 4th order of the dispersion error for (k/kh -1) can be obtained by the proposed approach. For the harmonic waves propagating along the Cartesian axes (n 2 = 0 for the x-axis and n 2 = 1 for the y-axis), Eq. (36) reduces to ( ) (kh h)4 k =1− + O (kh h)6 , (37) kh 180 i.e., it is exactly the same as that for the new quadratic elements with the lumped mass matrix in the 1-D case; see [1]. It is interesting to note that the order of the dispersion error in Eq. (36) cannot be improved with the use of the

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

245

three arbitrary coefficients a j ( j = 1, 2, 3) because for the particular directions of harmonic waves with n 2 = 0 or n 2 = 1, the leading term of the dispersion error is independent of the coefficients a j ; see Eq. (37). The leading term of the dispersion error in Eq. (36) is the polynomial function of n 2 with some undefined coefficients a j ( j = 1, 2, 3). Additional constraints can be imposed for the coefficients of this polynomial function. For example, it can be imposed that all coefficients containing any power of the term n 2 are zero. In this case Eqs. (35) and (36) can be rewritten as follows: m 1 = a1 , k1 = (60a2 + 41a1 )/10, k2 = 8a2 /3 − 14a1 /45, k4 = a2 /6 − a1 /360, k5 = −2a2 /3 + 2a1 /45, k 6 = a2 ,

k3 = −4a2 − 4a1 /5, (38)

and ( ) (kh h)4 k =1− + d(kh h)6 + O (kh h)8 , kh 180 ( )2 ( ) 840a2 n 2 2 − 1 n 2 4 + a1 −32n 2 8 + 64n 2 6 − 40n 2 4 + 8n 2 2 + 5 , (39) d= 10080a1 i.e., in this case the leading term of the dispersion error is independent of n 2 and the numerical anisotropy is reduced. a1 and a2 in Eqs. (38) and (39) are arbitrary coefficients. Similar to the conventional and new isogeometric elements with the non-diagonal mass matrix in Sections 2.1 and 2.2, the numerical dispersion error e = khk−k is plotted in Figs. 3 and 4 in the logarithmic scale for the conventional and new isogeometric elements with the diagonal mass matrix. For convenience, for the new elements in Figs. 3b and 4(d)–(f) the coefficient a2 in Eq. (38) is calculated from the condition that the term containing the second power of n 2 in coefficient d of Eq. (39) is zero (in this case, the dispersion error is independent of coefficients a1 and a2 ). The dispersion error for the conventional elements with the lumped mass matrix has only the second order of accuracy (it is calculated similar to that for the conventional elements with the non-diagonal mass matrix in Section 2.1). As can be seen, the new elements with the lumped mass matrix significantly reduce the numerical dispersion error which is very large for the conventional isogeometric elements with the lumped mass matrix. Let us find the elemental mass and stiffness matrices that yield the stencil equation (19) with the coefficients m 1 and k j ( j = 1, 2, 3, 4, 5, 6) given by Eq. (35). Similar to the previous Section, using Eqs. (27) and (28) for the elemental mass and stiffness matrices, using Eq. (29) with m 2 = m 3 = m 4 = m 5 = m 6 = 0, Eqs. (30) and (35), the coefficients of the mass and stiffness matrices in Eqs. (27) and (28) can be found for the new 2-D quadratic isogeometric elements with the lumped mass matrix and the reduced dispersion. They have the following form: d1 = (1 − 4a3 − a7 )/4, d5 = 0, d6 = −a4 , d8 = a4 , d9 = a5 ,

d2 = −a5 /2, d3 = 0, d7 = a3 , d10 = a6 , d11 = a7 ,

d4 = −a6 /2, (40)

and f1 f2 f4 f6 f9

= (5 + 72a1 − 2a10 + 60a2 − 8a8 )/8, = −4a1 + 1/2(−7a2 + a8 ) + a9 , f 3 = 1/24 − a1 − a2 /2, = (−60a1 + 3a10 − 4(1 + 12a2 + 3a9 ))/24, f 5 = a1 /2, = −1/3 + 4(a1 + a2 ) − a9 , f 7 = a8 , f 8 = a9 , = 8(a1 + a2 ) − a8 − 2a9 , f 10 = 1/3 − 4a1 − a10 /4 − 4a2 + a9 ,

f 11 = a10 ,

(41)

where a j (j=1,2, . . . ,10) are ten arbitrary coefficients. The dispersion error is independent of the values of these coefficients. At the derivation of Eqs. (40) and (41), Eqs. (33) and (34) have been also used with U0 = (u 0 , u 0 , u 0 , u 0 , u 0 , u 0 , u 0 , u 0 , u 0 )T and V 0 = (v0 , v0 , v0 , v0 , v0 , v0 , v0 , v0 , v0 )T where U0 and V 0 are the vectors with nine equal components. It is necessary to mention that the elemental mass matrix is not necessarily diagonal (however, the global mass matrix is diagonal independent of the values of coefficients a j (j=1, 2, . . . ,9)). Remark. According to the procedure for the calculation of the optimal coefficients m 1 and k j (see Eq. (35)) of the stencil equation(19), the 4th order of the dispersion error in Eq. (36) is the maximum order for all quadratic isogeometric elements with the lumped mass matrix.

246

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

Fig. 3. The numerical dispersion error in the logarithmic scale Log10 |e| (with e = khk−k ) as a function of the mesh size kh h and angle θ . (a) and (b) correspond to the conventional and new quadratic isogeometric elements with the diagonal mass matrix.

3. Dispersion analysis of cubic isogeometric elements in the 2-D case 3.1. Cubic isogeometric elements with reduced dispersion For the derivation of the cubic isogeometric elements with reduced dispersion, the stencil equation is assumed similar to that for the conventional cubic isogeometric elements calculated with the help of 4 × 4 = 16 neighboring elements (the derivations are similar to those in Section 2.2 for the quadratic isogeometric elements). It includes 49 degrees of freedom u i, j (i = A − 3, A − 2, A − 1, A, A + 1, A + 2, A + 3 and j = B − 3, B − 2, B − 1, B, B + 1, B + 2, B + 3) which are close to the degree of freedom u A,B and has the following form: k 2 h 2 {m 1 u A,B + m 2 [(u (A−1),(B+1) + u (A+1),(B−1) ) + (u (A+1),(B+1) + u (A−1),(B−1) )] + m 3 [(u (A−1),B + u (A+1),B ) + (u A,(B+1) + u A,(B−1) )] + m 4 [(u (A−2),(B−2) + u (A+2),(B+2) ) + (u (A−2),(B+2) + u (A+2),(B−2) )] + m 5 [(u (A+2),(B−1) + u (A−2),(B+1) ) + (u (A−1),(B−2) + u (A+1),(B+2) )

247

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

a

d 1

2

3

4

5

–1.0

6

–1.5

–0.5

–2.0 –2.5

–1.0

–3.0 –1.5

–3.5

–2.0

b

1

2

3

4

5

6

1

2

3

4

5

6

1

2

3

4

5

6

e 1

2

3

4

5

–1.0

6

–1.5

–0.5

–2.0 –2.5

–1.0

–3.0 –1.5

–3.5

–2.0

c

f 1 –0.5 –1.0

2

3

4

5

6

–1.0 –1.5 –2.0 –2.5 –3.0

–1.5 –2.0

–3.5

Fig. 4. The numerical dispersion error in the logarithmic scale Log10 |e| (with e = khk−k ) as a function of angle θ. (a, b, c) and (d, e, f) correspond to the conventional and new quadratic isogeometric elements with the diagonal mass matrix at kh h = 2 (a, d), kh h = 1 (b, e) and kh h = 0.5 (c, f).

+ (u (A−1),(B+2) + u (A+1),(B−2) ) + (u (A−2),(B−1) + u (A+2),(B+1) )] + m 6 [(u (A−2),B + u (A+2),B ) + (u A,(B+2) + u A,(B−2) )] + m 7 [(u (A−3),B + u (A+3),B ) + (u A,(B−3) + u A,(B+3) )] + m 8 [(u (A−3),(B−3) + u (A+3),(B+3) ) + (u (A+3),(B−3) + u (A−3),(B+3) )] + m 9 [(u (A−3),(B−1) + u (A+3),(B+1) ) + (u (A−3),(B+1) + u (A+3),(B−1) ) + (u (A−1),(B−3) + u (A+1),(B+3) ) + (u (A+1),(B−3) + u (A−1),(B+3) )] + m 10 [(u (A−3),(B−2) + u (A+3),(B+2) ) + (u (A−3),(B+2) + u (A+3),(B−2) ) + (u (A−2),(B−3) + u (A+2),(B+3) ) + (u (A+2),(B−3) + u (A−2),(B+3) )]} − {k1 u A,B + k2 [(u (A−1),(B+1) + u (A+1),(B−1) ) + (u (A+1),(B+1) + u (A−1),(B−1) )] + k3 [(u (A−1),B + u (A+1),B ) + (u A,(B+1) + u A,(B−1) )] + k4 [(u (A−2),(B−2) + u (A+2),(B+2) ) + (u (A−2),(B+2) + u (A+2),(B−2) )]

248

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

+ k5 [(u (A+2),(B−1) + u (A−2),(B+1) ) + (u (A−1),(B−2) + u (A+1),(B+2) ) + (u (A−1),(B+2) + u (A+1),(B−2) ) + (u (A−2),(B−1) + u (A+2),(B+1) )] + k6 [(u (A−2),B + u (A+2),B ) + (u A,(B+2) + u A,(B−2) )] + k7 [(u (A−3),B + u (A+3),B ) + (u A,(B−3) + u A,(B+3) )] + k8 [(u (A−3),(B−3) + u (A+3),(B+3) ) + (u (A+3),(B−3) + u (A−3),(B+3) )] + k9 [(u (A−3),(B−1) + u (A+3),(B+1) ) + (u (A−3),(B+1) + u (A+3),(B−1) ) + (u (A−1),(B−3) + u (A+1),(B+3) ) + (u (A+1),(B−3) + u (A−1),(B+3) )] + k10 [(u (A−3),(B−2) + u (A+3),(B+2) ) + (u (A−3),(B+2) + u (A+3),(B−2) )]} = 0,

(42)

where the unknown coefficients m j and k j ( j = 1, 2, . . . , 10) can be calculated similar to the conventional elements in terms of the elemental mass and stiffness matrices (at this point, these matrices are not defined). As can be seen, the stencil equation (42) for the cubic elements can be obtained from the stencil equation (19) for the quadratic elements by the addition of 8 terms with the coefficients m j and k j ( j = 7, 8, 9, 10). Similar to the stencil equation for the quadratic elements, the symmetry of the coefficients m j and k j in the stencil equation (42) is assumed for the degrees of freedom u j,l ( j = A − 3, A − 2, A − 1, A, A + 1, A + 2, A + 3 and l = B − 3, B − 2, B − 1, B, B + 1, B + 2, B + 3) symmetrically located with respect to the degree of freedom u A,B . As mentioned in Section 2, the optimal coefficients m j and k j ( j = 1, 2, . . . , 10) in Eq. (42) can be calculated by the insertion of Eq. (16) with kh = k into Eq. (42) and the minimization of the residual of the obtained expression (if the dispersion error is zero then the residual is exactly equal to zero). After the insertion of Eq. (16) with kh = k into Eq. (42), Eq. (42) can be reduced to the following form: R = R˜ exp(±ik(n 1 h A + n 2 h B)),

(43)

2 2

R˜ = k h {m 1 + 2m 2 [cos[k(n 1 h − n 2 h)] + cos[k(n 1 h + n 2 h)]] + 2m 3 [cos[kn 1 h] + cos[kn 2 h]] + 2m 4 [cos[2k(n 1 h + n 2 h)] + cos[2k(n 1 h − n 2 h)]] + 2m 5 [cos[k(2n 1 h − n 2 h)] + cos[k(n 1 h + 2n 2 h)] + cos[k(n 1 h − 2n 2 h)] + cos[k(2n 1 h + n 2 h)]] + 2m 6 [cos[2kn 1 h] + cos[2kn 2 h]] + 2m 7 [cos[3kn 1 h] + cos[3kn 2 h]] + 2m 8 [cos[3k(n 1 h + n 2 h)] + cos[3k(n 1 h − n 2 h)]] + 2m 9 [cos[k(3n 1 h − n 2 h)] + cos[k(n 1 h + 3n 2 h)] + cos[k(n 1 h − 3n 2 h)] + cos[k(3n 1 h + n 2 h)]] + 2m 10 [cos[k(3n 1 h − 2n 2 h)] + cos[k(2n 1 h + 3n 2 h)] + cos[k(2n 1 h − 3n 2 h)] + cos[k(3n 1 h + 2n 2 h)]]} − {k1 + 2k2 [cos[k(n 1 h − n 2 h)] + cos[k(n 1 h + n 2 h)]] + 2k3 [cos[kn 1 h] + cos[kn 2 h]] + 2k4 [cos[2k(n 1 h + n 2 h)] + cos[2k(n 1 h − n 2 h)]] + 2k5 [cos[k(2n 1 h − n 2 h)] + cos[k(n 1 h + 2n 2 h)] + cos[k(n 1 h − 2n 2 h)] + cos[k(2n 1 h + n 2 h)]] + 2k6 [cos[2kn 1 h] + cos[2kn 2 h]] + 2k7 [cos[3kn 1 h] + cos[3kn 2 h]] + 2k8 [cos[3k(n 1 h + n 2 h)] + cos[3k(n 1 h − n 2 h)]] + 2k9 [cos[k(3n 1 h − n 2 h)] + cos[k(n 1 h + 3n 2 h)] + cos[k(n 1 h − 3n 2 h)] + cos[k(3n 1 h + n 2 h)]] + 2k10 [cos[k(3n 1 h − 2n 2 h)] + cos[k(2n 1 h + 3n 2 h)] + cos[k(2n 1 h − 3n 2 h)] + cos[k(3n 1 h + 2n 2 h)]]} = 0.

(44)

For the simplification of the expressions in each parenthesis in Eq. (42), the following trigonometric formula is used: exp(ia) + exp(ib) = 2 exp(i(a + b)/2) cos[(a − b)/2]. The residual R of the dispersion equation can be decreased by ˜ see Eq. (43). Expanding the cosine functions in the expression for R˜ (see Eq. (44)) into a Taylor the reduction of R; series at small h ≪ 1, combining the terms with the same order of h, and equating zero the first seven terms with the lowest orders of h, the coefficients m j and k j ( j = 1, 2, . . . , 10) can be found from a system of the corresponding algebraic equations. They are: m 1 = (16(1491994905a1 + 3781295700a2 + 113273090a3 + 209126636a4 ))/12314025, m 2 = (211116315a1 + 544893600a2 + 14818570a3 + 28750228a4 )/456075, m 3 = (2(17413880a1 + 42267825a2 + 1687640a3 + 1459931a4 ))/50675,

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

249

m 4 = (4(751395a1 + 2562600a2 + 24610a3 + 296644a4 ))/456075, m 5 = (6516240a1 + 9449475a2 + 1151120a3 + 815798a4 )/456075, m 6 = (8(893205a1 + 2359700a2 + 61765a3 + 76606a4 ))/50675, m 7 = (2(29464680a1 + 53138325a2 + 1543240a3 + 1744921a4 ))/12314025, m 8 = (1790745a1 + 6529200a2 + 35210a3 + 418484a4 )/12314025, m 9 = −a1 ,

m 10 = −a2 ,

k1 = (16(65(23328516(2a1 + 5a2 ) + 3067151a3 ) + 521338456a4 ))/7388415, k2 = (6033480a1 + 15083700a2 − 1169485a3 + 3695636a4 )/30405, k3 = −((16(30(4711(2a1 + 5a2 ) − 904a3 ) + 156587a4 ))/10135), k4 = (4(90960a1 + 227400a2 − 22825a3 + 6812a4 ))/30405, k5 = (160(−59907(2a1 + 5a2 ) − 6944a3 ) − 2212793a4 )/30405, k6 = −(16(41070a1 + 102675a2 + 18625a3 + 1859a4 ))/10135, k7 = −(4(72454440(2a1 + 5a2 ) + 8792000a3 + 12253031a4 ))/7388415, k8 = (−2603640a1 − 6509100a2 − 52783a3 − 538828a4 )/1477683, k9 = −a3 ,

k10 = −a4 ,

(45) √

where a1 , a2 , a3 and a4 are four arbitrary coefficients and the expression n 1 = 1 − n 22 is used at the derivation of Eq. (45). Expanding the right-hand side of the dispersion equation (similar to the derivation of Eq. (18)) into a Taylor series at small h ≪ 1 with the coefficients given by Eq. (45) we get ( ) k d(kh h)12 =1− + O (kh h)14 , kh (1169103936000(233a4 + 170(a3 + 30a2 + 12a1 )))

(46)

d = (6497910(a3 + 30a2 + 12a1 ) + 10n 22 (−1 + n 22 )(43153536(5a2 + 2a1 ) + 2a3 (3596128 + (−1 + n 2 )n 22 (1 + n 2 )(20288599 + 25939748n 22 (−1 + n 22 ))) + 3n 22 (−1 + n 22 )(344921222a1 + 451001224a1 n 22 (−1 + n 22 ) + 5a2 (172460611 + 42888182n 22 (−1 + n 22 )))) + a4 (8905959 + 2(−1 + n 2 )n 22 (1 + n 2 )(−11326943 + (−1 + n 2 )n 22 (1 + n 2 )(486773071 + 105394142n 22 (−1 + n 22 ))))),

(47)

i.e., in the 2-D case the 12th order of the dispersion error for (k/kh -1) can be obtained by the proposed approach. It is necessary to note that for the harmonic waves propagating along the Cartesian axes (n 2 = 0 for the x-axis and n 2 = 1 for the y-axis), Eq. (46) reduces to ( ) 12741(kh h)12 k =1− + O (kh h)14 , (48) kh 389701312000 i.e., it is exactly the same as that for the cubic elements with reduced dispersion in the 1-D case; see [1]. It is interesting to note that the order of the dispersion error in Eq. (46) cannot be improved with the use of the four arbitrary coefficients a j ( j = 1, 2, 3, 4) because for the particular directions of harmonic waves with n 2 = 0 or n 2 = 1, the leading term of the dispersion error is independent of the coefficients a j ; see Eq. (48). The leading term of the dispersion error in Eq. (46) is the polynomial function of n 2 with some undefined coefficients a j ( j = 1, 2, 3, 4). Additional constraints can be imposed for the coefficients of this polynomial function. For example, it can be imposed that all coefficients containing any power of the term n 2 are zero. In this case Eqs. (45) and (46) can be rewritten as follows: m 1 = 596306519225a2 /84850119, m 3 = 8111366000a2 /3142597, m 5 = 503422545a2 /3142597,

m 2 = 71866926375a2 /50281552, m 4 = 511593927a2 /12570388, m 6 = 750552715a2 /3142597,

m 7 = 293353820a2 /84850119,

m 8 = 200141201a2 /1357601904,

m 9 = 127330055a2 /50281552,

m 10 = a2 ,

250

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

a

b –4 –6 –8 –10

c

1

2

3

4

5

6

1

2

3

4

5

6

d –4

–4

–6

–6

–8

–8

–10

–10 1

2

3

4

5

6

Fig. 5. The numerical dispersion error in the logarithmic scale Log10 |e| (with e = khk−k ) of the new cubic isogeometric elements with the non-diagonal mass matrix as a function of the mesh size kh h and angle θ (a) as well as at kh h = 2 (b), kh h = 1 (c) and kh h = 0.5 (d).

k1 = 6283962950825a2 /254550357,

k2 = −29266535025a2 /50281552,

k3 = −3203349100a2 /3142597,

k4 = −3052027125a2 /12570388,

k5 = −3594049200a2 /3142597,

k6 = −5514016025a2 /3142597,

k7 = −24494241100a2 /254550357, k9 = −3432850225a2 /50281552,

k8 = −10994897725a2 /4072805712, k10 = −78665300a2 /3142597,

(49)

and ) ( 12741(kh h)12 d(kh h)14 k =1− + + O (kh h)16 , kh 389701312000 203190264076800

(50)

d = −1344161 − (−1 + n 2 )n 22 (1 + n 2 )(3261580 + (−1 + n 2 )n 22 (1 + n 2 )(607031 + 12(−1 + n 2 )n 22 (1 + n 2 )(−99116 + 214825n 22 (−1 + n 22 )))),

(51)

i.e., in this case the leading term of the dispersion error is independent of n 2 and the numerical anisotropy is reduced. a2 in Eq. (49) is an arbitrary coefficient that can be easily found from the condition of the preservation of the kinetic energy (see Eq. (33)). Similar to the quadratic isogeometric elements in Section 2.2, the numerical dispersion error e = khk−k is calculated for the new cubic isogeometric elements with the coefficients given by Eq. (49) and is plotted in Fig. 5 in the logarithmic scale. As can be seen, compared with the new quadratic isogeometric elements in Section 2.2, the dispersion error for the new cubic elements is much smaller at the same mesh size kh h and decreases much faster with the decrease in the mesh size kh h; see Figs. 2 and 5. Remark. According to the procedure for the calculation of the optimal coefficients m j and k j (see Eq. (45)) of the stencil equation (42), the 12th order of the dispersion error in Eq. (46) is the maximum order for all cubic isogeometric elements in the 2-D case.

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

251

3.2. Cubic isogeometric elements with lumped mass matrix and reduced dispersion It has been shown in Part 1 of the paper (see [1]) that for the conventional stiffness matrix and the lumped mass matrix in the 1-D case, the dispersion error cannot exceed the second order of accuracy. For the harmonic waves propagating in the 2-D case along the Cartesian axes, the dispersion errors for the conventional stiffness matrix and the lumped mass matrix are the same as those in the 1-D case; i.e., the dispersion error for the cubic elements with the conventional stiffness matrix and the lumped mass matrix cannot exceed the second order of accuracy in the 2-D case as well. However, similar to the 1-D case (see [1]) the accuracy of the results with the lumped mass matrix can be improved by the modification of the stiffness matrix. Let us consider again the stencil equation (42) with the diagonal mass matrix (m 2 = m 3 = m 4 = m 5 = m 6 = m 7 = m 8 = m 9 = m 10 = 0) and the undefined coefficients k j ( j = 1, 2, . . . , 10). For finding the optimal coefficients m 1 and k j ( j = 1, 2, . . . , 10), the approach described in the previous Section 2.2 is used. Expanding the cosine functions in the expression for R˜ (see Eq. (44)) into a Taylor series at small h ≪ 1, combining the terms with the same order of h, and equating zero the first four terms with the lowest orders of h, the coefficients m 1 and k j ( j = 1, 2, . . . , 10) can be found from a system of the corresponding algebraic equations. They are: m 1 = a1 , k1 = (1/90)(11520a5 + 540a2 + 29520a3 + 720a4 + 409a1 ), k2 = (1/15)(320a5 + 40a2 + 945a3 − 30a4 − 6a1 ), k3 = −62a5 − 4a2 − 162a3 − 2a4 − 9a1 /10, k4 = −32a5 /3 + a2 /6 − 18a3 − 2a4 − a1 /40, k5 = 29a5 /3 − 2a2 /3 + 18a3 + 2a4 + a1 /10, k 6 = a2 , k7 = −2a5 − 2a3 − 2a4 − a1 /90, k 8 = a3 , k 9 = a4 , k10 = a5 , (52) √ where a j ( j = 1, 2, . . . , 5) are five arbitrary coefficients and the expression n 1 = 1 − n 22 is used at the derivation of Eq. (52). Expanding the right-hand side of the dispersion equation (similar to the derivation of Eq. (18)) into a Taylor series at small h ≪ 1 with the coefficients given by Eq. (52) we get ( ) k (kh h)6 d =1− + O (kh h)8 , kh 3360a1 ( ) ( ( ) ( ) ) d = 280 n 2 2 − 1 n 2 2 8a5 5n 2 4 − 5n 2 2 + 3 + n 2 2 − 1 n 2 2 (−(a2 − 24a4 )) + 6(9a3 + a4 ) ( ) ( )( )2 + 3a1 4n 2 2 n 2 2 − 1 1 − 2n 2 2 + 1 ,

(53)

i.e., in the 2-D case the 6th order of the dispersion error for (k/kh -1) can be obtained by the proposed approach. For the harmonic waves propagating along the Cartesian axes (n 2 = 0 for the x-axis and n 2 = 1 for the y-axis), Eq. (53) reduces to ( ) (kh h)6 k =1− + O (kh h)8 , (54) kh 1120 i.e., it is exactly the same as that for the new cubic elements with the lumped mass matrix and the reduced dispersion in the 1-D case; see [1]. Remark. According to the procedure for the calculation of the optimal coefficients m 1 and k j (see Eq. (52)) of the stencil equation (42), the 6th order of the dispersion error in Eq. (53) is the maximum order for all cubic isogeometric elements with the lumped mass matrix in the 2-D case. 4. Numerical examples The 1-D wave propagation problems solved in [1] by the conventional and new isogeometric elements showed a much higher accuracy of the new isogeometric elements with reduced dispersion in the 1-D case. Here, the advantages of the new isogeometric elements with reduced dispersion will be demonstrated for three wave propagation problems in the 2-D case. The first two problems have continuous solutions and will be solved without the filtering stage. The third problem is related to impact loading and will be solved by the two-stage time-integration technique with the filtering stage; see [13,14].

252

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

Fig. 6. A 2-D plate considered in Sections 4.1 and 4.2.

4.1. A standing wave in a 2-D plate Similar to [15], the problem of a standing wave is used for the demonstration of the order of convergence of the new numerical technique. The wave velocity is chosen to be c = 1. A standing wave in the 2-D case can be given by the following exact solution to the wave equation: √ (55) u(x, y, t) = sin(aπ x)sin(aπ y)cos( 2aπ t), where a = 5 is used. Let us consider a square plate with the dimensions 1 × 1; see Fig. 6. The initial conditions at time t = 0 and the boundary conditions in terms of displacements at x = 0, x = 1, y = 0 and y = 1 are selected according to the exact solution, Eq. (55); i.e., the initial displacements are u(x, y, t = 0) = sin(aπ x)sin(aπ y), the initial velocities are zero and all boundaries of the plate are fixed. The observation time is selected to be T = 1/a. The problem is solved by the conventional and new isogeometric elements with the non-diagonal and diagonal mass matrices on meshes with uniformly spaced control points. For the time integration the trapezoidal rule for the elements with the non-diagonal mass matrices and the central-difference method for the elements with the diagonal mass matrices are used with very small time increments at which the error in time is very small and can be neglected. This means that the difference between numerical and analytical solutions is only related to the space-discretization error. Fig. 7 shows the convergence of the error in the displacement e = |u exact (x = 0.5, y = 0.5, t = T ) − u num (x = 0.5, y = 0.5, t = T )| at point O at mesh refinement where u exact (x = 0.5, y = 0.5, t = T ) and u num (x = 0.5, y = 0.5, t = T ) are the exact and numerical displacements in the center of the plate at the observation time T . We should mention that the maximum error in the displacements occurs in the center of the plate. h in Fig. 7 is the distance between two consecutive control points along the x- and y- axes. At the same h, the meshes with the conventional and new quadratic isogeometric elements with the diagonal and non-diagonal mass matrices include the same number of degrees of freedom. The results in Fig. 7 are plotted in the logarithmic scale. Therefore, the slopes of the curves at small h shown in Fig. 7 correspond to the order of convergence (the order of accuracy) of the conventional and new isogeometric elements. The new isogeometric elements with the non-diagonal and diagonal mass matrices yield much more accurate results than those obtained by the corresponding conventional elements (compare curves 1 and 2 as well as curves 3 and 4). The results in Fig. 7 are in good agreement with the theoretical order of accuracy of the conventional and new isogeometric elements with the non-diagonal and diagonal mass matrices reported in Section 2 of the paper.

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

253

Fig. 7. The error in displacement e = |u exact (x = 0.5, y = 0.5, t = T ) − u num (x = 0.5, y = 0.5, t = T )| at point O as a function of the mesh size h in the logarithmic scale; see Fig. 6 for the location of point O. Curves 1 and 3 correspond to the conventional quadratic isogeometric elements with the non-diagonal and diagonal mass matrices. Curves 2 and 4 correspond to the new quadratic isogeometric elements with the non-diagonal and diagonal mass matrices. The slopes of the curves at small h show the order of convergence of the corresponding techniques. Symbols □, ⃝, × and ⋄ correspond to the results on the uniform meshes used in calculations.

4.2. Propagation of sinusoidal pulse in 2-D plate Let us consider a square plate with the dimensions 1 × 1; see Fig. 6. The following exact solution is used for the wave propagation in the 2-D plate: √ u(x, y, t) = cos( 2kπ t + kπ x)cos(kπ y) (56) with k = 6. The initial conditions at t = 0 and the boundary conditions at x = 0, x = 1, y = 0, y = 1 are defined according to the exact solution, Eq. (56). The wave velocity is chosen to be c = 1. The observation time is chosen to be T = 3. The problem is solved by the conventional and new quadratic isogeometric elements with the non-diagonal mass matrix on meshes with uniformly spaced control points with 17 × 17 = 259 and 33 × 33 = 1089 degrees of freedom (see Figs. 8 and 9) as well as with the lumped mass matrix on meshes with uniformly spaced control points with 21×21 = 441, 33×33 = 1089, 51×51 = 2601, 101×101 = 10 201 and 151×151 = 22 801 degrees of freedom (see Fig. 10). For the time integration the trapezoidal rule is used for the isogeometric elements with the non-diagonal mass matrix and the central-difference method is used for the isogeometric elements with the diagonal mass matrix. Very small time increments are used in calculations; i.e., the error in time is very small and can be neglected. This means that the difference between numerical and exact solutions is only related to the space-discretization error. Fig. 8 shows the displacement distribution along lines M N , P Q and AC for the conventional and new quadratic isogeometric elements with the non-diagonal mass matrix; see Fig. 6 for M N , P Q and AC. As can be seen, on the same mesh with 17 × 17 = 259 degrees of freedom the new quadratic isogeometric elements (curves 2) yield much more accurate results compared to those for the conventional elements (curves 1); compare also with the exact solutions (curves 3). At mesh refinement, the numerical results in Fig. 9 converge to the exact solution for both the new (d,e,f) and conventional (a,b,c) isogeometric elements. Moreover, for the scale used in Fig. 9, the numerical results for the new quadratic isogeometric elements on the mesh with 33 × 33 = 1089 degrees of freedom practically coincide with the exact solution; see curves 2 and 3 in Fig. 9(d)–(f). In the case of the lumped mass matrix, the numerical results for the conventional isogeometric elements on coarse meshes are very inaccurate and mesh refinement is needed for the increase in accuracy; see Fig. 10(a)–(c). The new elements with the lumped mass matrix can be used with coarse meshes. At mesh refinement the numerical results obtained by the new elements with the lumped mass matrix fastly converge to the exact solution; see Fig. 10(d)– (f) (in contrast to the slow convergence of the numerical results obtained by the conventional isogeometric elements with the lumped mass matrix; see Fig. 10(a)–(c)). These numerical results are consistent with the theoretical findings related to the accuracy of the new and conventional isogeometric elements presented in the previous sections of the paper.

254

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

Fig. 8. The displacement distribution along lines M N (a), P Q (b) and AC (c) at the observation time T = 3 as a function of the distance R from points M (a), P (b) and A (c) for the conventional (curves 1) and new (curves 2) quadratic isogeometric elements with the non-diagonal mass matrix; see Fig. 6 for M N , P Q and AC. Curves 3 correspond to the exact solutions. The mesh with uniformly spaced control points with 17 × 17 = 259 degrees of freedom is used.

4.3. 2-D elastic plate under a suddenly applied concentrated load Let us consider a square plate of length 1 × 1 instantly loaded in the center; see Fig. 11(a). The concentrated loading is implemented in terms of the initial conditions for function u(x, t = 0) as follows: the initial value of one control variable corresponding to the basis function in the center of the plate equals one U O,O (t = 0) = 1, all other initial values of the control variables U are zero; the initial velocities are zero for the entire domain. Zero boundary conditions u(x, t) = 0 are applied along the entire boundary AB, BC, C D and AD. The wave velocity is chosen to be c = 1. The observation time is chosen to be T = 0.4. During this time the waves travel from point O to the boundary but do not reach it. Therefore, zero boundary conditions do not affect the solution at time T = 0.4. By symmetry, the exact solution to this problem at any point depends on the radius with the center at point O and is independent of the angle between the radius and the x-axis. Therefore, this is a good benchmark problem for the

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

255

Fig. 9. The displacement distribution along lines M N (a, d), P Q (b, e) and AC (c, f) at the observation time T = 3 as a function of the distance R from points M (a, d), P (b, e) and A (c, f) for the conventional (a, b, c) and new (d, e, f) quadratic isogeometric elements with the non-diagonal mass matrix; see Fig. 6 for M N , P Q and AC. Curves 3 correspond to the exact solutions. The meshes with uniformly spaced control points with 17 × 17 = 259 (curves 1) and 33 × 33 = 1089 (curves 2) degrees of freedom are used.

study of the numerical dispersion error of the spatial discretization in different directions. The problem is solved by the conventional and new quadratic isogeometric elements on meshes with uniformly spaced control points with 61 × 61 = 3721, 101 × 101 = 10 201, and 151 × 151 = 22 801 degrees of freedom. The non-diagonal and diagonal mass matrices are used. For comparison, the same problem is solved with the conventional quadratic finite elements on a mesh with 151 × 151 = 22 801 degrees of freedom and the consistent mass matrix. Due to the initial conditions that correspond to the impact loading, all modes are excited and spurious oscillations appear in the numerical results; see Figs. 11(b), 12(a), 13(a), 14(a), 16(a)–21(a). Therefore, the two-stage time-integration procedure with the basic computations and the filtering stage (as described in the papers [13,14,16]) is applied for accurate and non-oscillatory solutions. For the time integration at the stage of basic computations, the trapezoidal rule is used for the isogeometric elements with the non-diagonal mass matrix and the central-difference method is used for the isogeometric elements with the diagonal mass matrix. Very small time increments are used at basic computations; i.e., the error in time is very small and can be neglected. This means that the difference between numerical and exact solutions is only related to the space-discretization error. For the filtering of spurious oscillations, the implicit TCG method with large numerical

256

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

Fig. 10. The displacement distribution along lines M N (a, d), P Q (b, e) and AC (c, f) at the observation time T = 3 as a function of the distance R from points M (a, d), P (b, e) and A (c, f) for the conventional (a, b, c) and new (d, e, f) quadratic isogeometric elements with the lumped mass matrix; see Fig. 6 for M N , P Q and AC. Curves 1, 2 and 3 in (a, b, c) correspond to the numerical results on the meshes with uniformly spaced control points with 51 × 51 = 2601 (curves 1), 101 × 101 = 10 201 (curves 2) and 151 × 151 = 22 801 (curves 3) degrees of freedom. Curves 1, 2 and 3 in (d, e, f) correspond to the numerical results on the meshes with uniformly spaced control points with 21 × 21 = 441 (curves 1), 33 × 33 = 1089 (curves 2) and 51 × 51 = 2601 (curves 3) degrees of freedom. Curves 4 correspond to the exact solutions.

dissipation developed in [13] is used at the filtering stage. For all problems, 10 uniform time increments (5 positive plus 5 negative time increments) are used at the filtering stage. This means that there is no real time integration at the filtering stage (the sum of 10 time increments used at the filtering stage is zero). As shown in [13], this procedure is equivalent to the multiplication of each ( ) velocity and displacement of the uncoupled system of the semi-discrete 2

2

5

(3+m) +Ω equations by a factor of (3+m) (where Ω = ω j ∆t and ω j are the eigen-frequencies of the semi-discrete 2 +(2+m)2 Ω 2 system, ∆t is the time increment as well as m = 15 is used) and does not require the modal decomposition and the calculation of eigen-frequencies. As can be seen, this factor is close to zero for large Ω and is close to unity for small Ω . The size ∆t of time increments at the filtering stage indirectly defines the amount of numerical dissipation and the range of the spurious oscillations to be filtered as well as the range of the actual frequencies left in the numerical results. The details related to the filtering stage can be found in the papers [13,14,16].

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

257

Fig. 11. 2-D plate under suddenly applied load in the center O (a). The distribution of the velocity v(x, y, 0.4) (b, c) at the observation time T = 0.4 after basic computations (b) and after the filtering stage (c). The conventional quadratic isogeometric elements with the non-diagonal mass matrix on the mesh with uniformly spaced control points with 151 × 151 = 22 801 degrees of freedom are used.

Remark. The numerical experiments in [16] showed that except the linear finite elements with the diagonal mass matrix, the other linear and high-order finite elements, spectral elements and isogeometric elements yield more accurate results at small time increments at the stage of basic computations; i.e., in this case the error in time does not compensate the space-discretization error. The numerical results for the velocity distribution at time T = 0.4 obtained by the conventional and new quadratic isogeometric elements as well as by the conventional quadratic finite elements with the non-diagonal mass matrices are shown in Figs. 11–16 where Figs. 11 and 13 show the isolines of the velocity for the 2-D square domain, and Figs. 12 and 14–16 show the distribution of the velocity along lines O E and OC (see Fig. 11(a)) as a function of distance R from the center O. As can be seen, after basic computations all numerical results include spurious highfrequency oscillations; see Figs. 11(b), 12(a), 13(a), 14(a) and 16(a). In contrast to the exact solution, the isolines in Fig. 11(b) and 13(a) are not circumferential and the distributions of the velocity along lines O E and OC are very different due to large spurious oscillations. It is difficult to compare the accuracy of the numerical results after the stage of basic computations due to the spurious oscillations. After the filtering stage, the comparison of accuracy is much easier because similar to the exact solution, the numerical results in the polar coordinate system with the center at point O are practically independent of the angle (i.e., the filtering stage also removes the numerical anisotropy due to the different numerical dispersion in different directions); see Figs. 11(c), 12(b), 13(b), 14(b) and 16(b) (curves 1 and 2 in 12(b), 14(b) and 16(b) are practically coincide). Because there is no analytical expression for the exact solution, the convergence of the numerical results at mesh refinement is shown in Fig. 15 for the isogeometric elements with reduced dispersion (higher amplitudes correspond to more accurate results; see curve 3 in Fig. 15). It can be also clearly seen that the isogeometric elements with reduced dispersion significantly improve the accuracy of the numerical results compared with the conventional quadratic isogeometric and finite elements; compare curves 1 and 2 in Figs. 12(b), 14(b) and 16(b) (the amplitudes of the curves in Fig. 14(b) are higher than those in Fig. 12(b) and Fig. 16(b). The numerical results in Figs. 12, 14 and 16 have been obtained on the meshes with the same number of degrees of freedom and the same bandwidth of the matrix of the final algebraic equations. Figs. 12(c) and 16(c) also show the results after the filtering stage for the conventional quadratic isogeometric and finite elements when the same range of frequencies as that for the quadratic isogeometric elements with reduced dispersion has been filtered (the same time increments at the filtering stage are used for the new and conventional elements). In this case the numerical results for the conventional isogeometric elements still have spurious oscillations and the results along

258

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

Fig. 12. The velocity distribution along lines OC (curves 1) and O E (curves 2) as a function of distance R from the center O after basic computations (a) and after the filtering stage (b, c); see Fig. 11(a) for OC and O E. The observation time T = 0.4 and the conventional quadratic isogeometric elements with the non-diagonal mass matrix on the mesh with uniformly spaced control points with 151 × 151 = 22801 degrees of freedom are used. (c) corresponds to the filtering of the same range of frequencies as that for the isogeometric elements with reduced dispersion; see the text for more explanations.

lines O E and OC are different in the range 0.35 ≤ R ≤ 0.45; see curves 1 and 2 in Figs. 12(c) and 16(c). This also shows that the isogeometric elements with reduced dispersion are more accurate than the conventional elements (the range of accurate frequencies for the conventional elements is smaller for the meshes with the same number of degrees of freedom). Despite a higher accuracy, the isogeometric elements with reduced dispersion do not require additional computational costs compared with the conventional isogeometric elements. Remark 1. It was shown in [3,5] that uniformly spaced control points yield more accurate results for structural dynamics problems than those obtained with a linear parametrization and the constant Jacobian determinant. However, even with uniformly spaced control points, the numerical results include spurious oscillations in basic computations (see Figs. 11(b), 12(a), 13(a), and 14(a)) and require the filtering stage.

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

259

Fig. 13. The distribution of the velocity v(x, y, 0.4) (a, b) at the observation time T = 0.4 after basic computations (a) and after the filtering stage (b). The new quadratic isogeometric elements with the non-diagonal mass matrix on the mesh with uniformly spaced control points with 151 × 151 = 22 801 degrees of freedom are used.

Fig. 14. The velocity distribution along lines OC (curves 1) and O E (curves 2) as a function of distance R from the center O after basic computations (a) and after the filtering stage (b); see Fig. 11(a) for OC and O E. The observation time T = 0.4 and the new quadratic isogeometric elements with the non-diagonal mass matrix on the mesh with uniformly spaced control points with 151 × 151 = 22 801 degrees of freedom are used.

Remark 2. It is interesting to note that at the derivation of the mass and stiffness matrices of the new isogeometric elements, the basis functions have not been used. Therefore, after the solution of the global semidiscrete system related to the new quadratic isogeometric elements with the non-diagonal mass matrix, we have also plotted the numerical results assuming that the global control variables yield the values of the field function u(x, t) at the nodes of a uniform mesh (i.e., similar to the conventional finite elements, at these nodes one basis function equals unity, all other basis functions equal zero). The results with this post-processing procedure are shown in Fig. 17. As can be seen, in this case the field function u(x, t) after the filtering stage is practically the same as that in Fig. 14(b). Similar to Figs. 11–15 with the non-diagonal mass matrix, Figs. 18–21 show the numerical results for the conventional and the new quadratic isogeometric elements with the lumped mass matrix. As can be seen from these figures, the isogeometric elements with reduced dispersion are much more accurate than the conventional

260

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

Fig. 15. The velocity distribution along O E at the observation time T = 0.4 after the filtering stage. Curves 1, 2 and 3 correspond to the meshes with uniformly spaced control points with 151 × 151 = 22 801, 101 × 101 = 10 201, and 61 × 61 = 3721 degrees of freedom, respectively. The quadratic isogeometric elements with reduced dispersion and the non-diagonal mass matrix are used.

5 4 3 2 1 0 -1 -2 -3 -4 -5 0

0.1

0.2

0.3

0.4

0.5

0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2 -0.25 -0.3 0

0.1

0.2

0.3

0.4

0.5

0

0.1

0.2

0.3

0.4

0.5

0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2 -0.25 -0.3

Fig. 16. The velocity distribution along lines OC (curves 1) and O E (curves 2) as a function of distance R from the center O after basic computations (a) and after the filtering stage (b, c); see Fig. 11(a) for OC and O E. The observation time T = 0.4 and the conventional quadratic finite elements with the consistent mass matrix on a uniform mesh with 151 × 151 = 22 801 degrees of freedom are used. (c) corresponds to the filtering of the same range of frequencies as that for the isogeometric elements with reduced dispersion; see the text for more explanations.

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

261

Fig. 17. The velocity distribution along lines OC (curves 1) and O E (curves 2) as a function of distance R from the center O after basic computations (a) and after the filtering stage (b, c); see Fig. 11(a) for OC and O E. The observation time T = 0.4 and the global system of semidiscrete equations with 151 × 151 = 22 801 degrees of freedom (the same as that for the new quadratic isogeometric elements with reduced dispersion) are used. Curves 1 and 2 in (c) correspond to the velocity distribution along line O E from Fig. 14(b) and from (b); see the text for more explanations.

quadratic isogeometric elements with the lumped mass matrix. It is necessary to note that the numerically estimated stability limits for the central difference method with the conventional and new isogeometric elements are comparable. However, for accurate numerical results at the loading considered here, the time increments in basic computations should be smaller than the stability limits. In order to show the performance of the new approach on non-uniform meshes, the problem is also solved on a circular domain of radius r = 0.5 with zero boundary conditions u(x, t) = 0 along the entire boundary and the same initial conditions as those for the square domain. According to the boundary and initial conditions, the exact solutions for the square and circular domains at time T = 0.4 are the same. An example of a non-uniform mesh for the circular domain is shown in Fig. 22. This mesh is created using a NURBS discretization with the corresponding

262

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

Fig. 18. The distribution of the velocity v(x, y, 0.4) (a, b) at the observation time T = 0.4 after basic computations (a) and after the filtering stage (b). The conventional quadratic isogeometric elements with the diagonal mass matrix on the mesh with uniformly spaced control points with 151 × 151 = 22 801 degrees of freedom are used.

Fig. 19. The velocity distribution along lines OC (curves 1) and O E (curves 2) as a function of distance R from the center O after basic computations (a) and after the filtering stage (b); see Fig. 11(a) for OC and O E. The observation time T = 0.4 and the conventional quadratic isogeometric elements with the diagonal mass matrix on the mesh with uniformly spaced control points with 151 × 151 = 22 801 degrees of freedom are used.

refinement parameter (see [17] for details). Fig. 23 shows the numerical results obtained by the conventional and new approaches with the non-diagonal mass matrices on a non-uniform mesh with 130 × 130 = 16 900 degrees of freedom at time T = 0.4 after the filtering stage (after the filtering stage the velocity distribution along lines OC and O E are practically the same for any approach and any mesh). For the comparison, Fig. 23 shows the numerical results obtained for the square domain on a uniform mesh with 130 × 130 = 16 900 degrees of freedom (similar to those in Figs. 12(b) and 14(b)). As can be seen, the numerical results for the new and conventional quadratic isogeometric elements are more accurate on a uniform mesh (see curves 1 and 2 as well as 3 and 4 in Fig. 23). On a non-uniform mesh the

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

263

Fig. 20. The distribution of the velocity v(x, y, 0.4) (a, b) at the observation time T = 0.4 after basic computations (a) and after the filtering stage (b). The new quadratic isogeometric elements with the diagonal mass matrix on the mesh with uniformly spaced control points with 151 × 151 = 22 801 degrees of freedom are used.

Fig. 21. The velocity distribution along lines OC (curves 1) and O E (curves 2) as a function of distance R from the center O after basic computations (a) and after the filtering stage (b); see Fig. 11(a) for OC and O E. The observation time T = 0.4 and the new quadratic isogeometric elements with the diagonal mass matrix on the mesh with uniformly spaced control points with 151 × 151 = 22 801 degrees of freedom are used.

new quadratic isogeometric elements yield more accurate results compared with those obtained by the conventional quadratic isogeometric elements (see curves 3 and 1 in Fig. 23). 5. Concluding remarks The development of the isogeometric elements with reduced dispersion is considered in the paper in the 2-D case. The new findings of the paper can be summarized as follows: • By the minimization of the order of the numerical dispersion error for the stencil equation with arbitrary coefficients, the dispersion error has been reduced from the order 2 p for the conventional isogeometric elements

264

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

Fig. 22. NURBS discretization of the circular domain: (a) control net; (b) physical mesh.

0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2 -0.25 -0.3 0

0.1

0.2

0.3

0.4

0.5

Fig. 23. The velocity distribution along line O E as a function of distance R from the center O after the filtering stage. The observation time T = 0.4. Curves 1, 2, 3, 4 correspond to the numerical results obtained by the conventional (curves 1 and 2) and new (curves 3 and 4) quadratic isogeometric elements on the uniform mesh for the square domain (curves 2 and 4) and on the nonuniform mesh for the circular domain (curves 1 and 3). All meshes include 130 × 130 = 16 900 degrees of freedom.

•

• •

• •

to the order 4 p for the new elements ( p is the order of the polynomial approximations). Because all coefficients of the stencil equation can be found from the analysis of the dispersion error, the order 4 p of the dispersion error is maximum possible for the considered form of the stencil equation. For the known results with the linear finite elements in the multidimensional case (e.g., see [11,12,18]) and the high-order finite and isogeometric elements in the 1-D case (e.g., see [8,9,10]), the order of the dispersion error has been improved from the order 2 p to the order 2 p + 2 (i.e., much smaller than 4 p). In contrast to [8,9,10], for the first time the improved order of the dispersion error is valid for the high-order elements in the general 2-D case (it is independent of the propagation direction). This will lead to a significant reduction in the computation time. The isogeometric elements with the lumped mass matrix have the second order of the dispersion error (e.g., see [3,5]) and are computationally inefficient. A special iterative procedure that requires additional computational costs was suggested in [19] in order to obtain the same order of the dispersion error for the lumped and consistent mass matrices. It has been shown in the paper that independent of the procedures for the calculation of the lumped mass matrix, the order of the dispersion error cannot be improved with the conventional stiffness matrix and no additional iterations. It has been also shown that the dispersion error for the lumped mass matrix can be improved from the second order to the order 2 p by the modification of the stiffness matrix (without additional computational costs). By the solution of the 2-D wave propagation problems, the computational efficiency of the new isogeometric elements with reduced dispersion used with the non-diagonal and diagonal mass matrices has been demonstrated. By the application of the two-stage time-integration technique developed in our recent papers [13,14,16], the accurate numerical results can be obtained for impact problems solved by the new isogeometric elements. This technique quantifies and filters the spurious oscillations. Due to the same structure of the stencil equations, the numerical approach suggested in the paper can be directly applied to the dynamics of Euler–Bernoulli beams and Kirchhoff plates. A new technique with reduced numerical dispersion has been derived on uniform meshes. Practically all publications related to the analytical study of the numerical dispersion of different space-discretization techniques for wave propagation are based on uniform meshes. Generally, it is not clear how to define the

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

265

Fig. A.24. The spatial locations of the degrees of freedom u i, j (i = A − 2, A − 1, A, A + 1, A + 2 and j = B − 2, B − 1, B, B + 1, B + 2) contributing to the stencil equation for the degree of freedom u A,B of the new and conventional quadratic isogeometric elements. The nine elements shown are used for the calculation of the stencil equation for the degree of freedom u A,B according to Eq. (3).

numerical wave velocity and the numerical dispersion error on non-uniform meshes because the harmonic waves (functions) do not meet the corresponding discrete system after the space discretization in this case. Therefore, the performance of the new and conventional isogeometric elements on non-uniform meshes is compared only numerically. Another alternative for complicated geometries that can be considered in the future is to use relatively large elements with reduced dispersion for the subdomains with uniform meshes and small elements for the subdomains with non-uniform meshes. • For the derivation of the mass and stiffness matrices of the new isogeometric elements, the specific form of the basis functions is not used (see also the recent results in [9,10]). This is partly similar to the virtual element method (e.g., see [20,21]) where the specific form of the basis functions is not used for the derivation of the final equations as well. The approach developed in the paper can be applied to other space-discretization techniques with a similar form of the stencil equations; e.g., to the high-order finite-difference method in order to reduce the numerical dispersion error. Acknowledgments The research has been supported in part by the Air Force Office of Scientific Research (contract FA9550-16-10177) and by Texas Tech University. Appendix A. The structure of the elemental mass and stiffness matrices for the quadratic isogeometric elements on uniform meshes Here, uniform meshes with the same spacing h in the x and y directions are considered. Fig. A.24 shows the spatial locations of 25 degrees of freedom (or more precisely, the basis functions related to the corresponding control variables) for 3 × 3 = 9 neighboring elements contributing to the formation of the stencil equation for the degree of freedom u A,B with the conventional and new quadratic isogeometric elements. For the new quadratic isogeometric elements the form of the stencil equation is the same as that for the conventional quadratic isogeometric elements in [3,5]. To simplify the representation of the elemental mass or stiffness matrices, let us use the following order of the degrees of freedom in the local vector of the control variables for any element:

266

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

U = {u A−1,B−1 , u A−1,B+1 , u A+1,B+1 , u A+1,B−1 , u A−1,B , u A,B+1 , u A+1,B , u A,B−1 , u A,B }T (here, the degrees of freedom for the central element in Fig. A.24 are used). In this case the elemental mass or stiffness matrix can be expressed in terms of 11 coefficients ai as shown below in Eq. (A.1): ⎞ ⎛ a1 a2 a3 a2 a4 a5 a5 a4 a6 ⎜a2 a1 a2 a3 a4 a4 a5 a5 a6 ⎟ ⎟ ⎜ ⎜a3 a2 a1 a2 a5 a4 a4 a5 a6 ⎟ ⎟ ⎜ ⎜a2 a3 a2 a1 a5 a5 a4 a4 a6 ⎟ ⎜ ⎟ ⎟ (A.1) Me (or K e ) = ⎜ ⎜a4 a4 a5 a5 a7 a8 a9 a8 a10 ⎟ . ⎜a5 a4 a4 a5 a8 a7 a8 a9 a10 ⎟ ⎟ ⎜ ⎜a5 a5 a4 a4 a9 a8 a7 a8 a10 ⎟ ⎟ ⎜ ⎝a4 a5 a5 a4 a8 a9 a8 a7 a10 ⎠ a6 a6 a6 a6 a10 a10 a10 a10 a11 Here, the symmetry of matrices Me and K e as well as the symmetry of the location of the first and second four degrees of freedom in vector U have been taken into account (therefore the first and second four rows in Eq. (A.1) include the same coefficients). For the implementation of the matrices Me and K e into computer codes, the following order of degrees of freedom in the local vector of the control variables is often used: U = {u A−1,B−1 , u A−1,B , u A−1,B+1 , u A,B−1 , u A,B , u A,B+1 , u A+1,B−1 , u A+1,B , u A+1,B+1 }T . In this case the mass or stiffness matrices can be obtained from Eq. (A.1) by the corresponding relocations of coefficients a j in Eq. (A.1); see Eqs. (27) and (28). Appendix B. Local truncation error in space for the new and conventional quadratic isogeometric elements Based on the semi-discrete system of Eqs. (2) in the time domain and Eq. (3), the stencil equation for the degree of freedom u A,B in the time domain for the new and conventional quadratic elements has the following form (the derivation is similar to that for Eq. (19)): h 2 {m 1 u¨ A,B + m 2 [(u¨ (A−1),(B+1) + u¨ (A+1),(B−1) ) + (u¨ (A+1),(B+1) + u¨ (A−1),(B−1) )] + m 3 [(u¨ (A−1),B + u¨ (A+1),B ) + (u¨ A,(B+1) + u¨ A,(B−1) )] + m 4 [(u¨ (A−2),(B−2) + u¨ (A+2),(B+2) ) + (u¨ (A−2),(B+2) + u¨ (A+2),(B−2) )] + m 5 [(u¨ (A+2),(B−1) + u¨ (A−2),(B+1) ) + (u¨ (A−1),(B−2) + u¨ (A+1),(B+2) ) + (u¨ (A−1),(B+2) + u¨ (A+1),(B−2) ) + (u¨ (A−2),(B−1) + u¨ (A+2),(B+1) )] + m 6 [(u¨ (A−2),B + u¨ (A+2),B ) + (u¨ A,(B+2) + u¨ A,(B−2) )]} + c2 {k1 u A,B + k2 [(u (A−1),(B+1) + u (A+1),(B−1) ) + (u (A+1),(B+1) + u (A−1),(B−1) )] + k3 [(u (A−1),B + u (A+1),B ) + (u A,(B+1) + u A,(B−1) )] + k4 [(u (A−2),(B−2) + u (A+2),(B+2) ) + (u (A−2),(B+2) + u (A+2),(B−2) )] + k5 [(u (A+2),(B−1) + u (A−2),(B+1) ) + (u (A−1),(B−2) + u (A+1),(B+2) ) + (u (A−1),(B+2) + u (A+1),(B−2) ) + (u (A−2),(B−1) + u (A+2),(B+1) )] + k6 [(u (A−2),B + u (A+2),B ) + (u A,(B+2) + u A,(B−2) )]} = 0,

(B.1)

where the coefficients m j and k j ( j = 1, 2, . . . , 6) are given by Eq. (22) for the new elements with reduced dispersion and are m 1 = 4356, m 2 = 676, m 3 = 1716, m 4 = 1, m 5 = 26, m 6 = 66, k1 = 15 840, k2 = −2080, k3 = 480, k4 = −40, k5 = −560, k6 = −1200 for the conventional quadratic isogeometric elements; see Section 2.2. Let us analyze the local truncation error in space for the stencil equation Eq. (B.1) with the coefficients given by Eq. (22) (similar to the analysis in the 1-D case in Part 1 of the paper; see [1]). The local truncation error is the residual of the stencil equation when the unknowns u j,l ( j = A − 2, A − 1, A, A + 1, A + 2 and l = B − 2, B − 1, B, B + 1, B + 2) correspond to the exact solution. Because at the derivation of the coefficients Eq. (22) of the stencil equation the basis functions have not been used then we assume that the control variables u j,l yield the values of the field function u(x, t) at the nodes of a uniform mesh (similar to the conventional finite elements); see also the Remark 2 and Fig. 17

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

267

in Section 4. Under this assumption, let us expand the exact solution u j,l and u¨ j,l ( j = A − 2, A − 1, A + 1, A + 2 and l = B − 2, B − 1, B + 1, B + 2) into a Taylor series at small h ≪ 1 as follows: ∂u A,B ∂u A,B ∂ 2 u A,B (±i h)2 ∂ 2 u A,B (±i h)(± j h) (±i h) + (± j h) + + 2 ∂x ∂y ∂x 2! ∂ x∂ y 2! ∂ 2 u A,B (± j h)2 ∂ 3 u A,B (±i h)3 ∂ 3 u A,B (±i h)2 (± j h) + + + ∂ y2 2! ∂x3 3! ∂ y∂ x 2 3! 3 2 3 3 ∂ u A,B (±i h)(± j h) ∂ u A,B (± j h) + 2 + + ··· , ∂y ∂x 3! ∂ y3 3! ∂ u¨ A,B ∂ 2 u¨ A,B (±i h)2 ∂ 2 u¨ A,B (±i h)(± j h) ∂ u¨ A,B (±i h) + (± j h) + + u¨ A±i,B± j = u¨ A,B + ∂x ∂y ∂x2 2! ∂ x∂ y 2! 2 2 3 3 ∂ u¨ A,B (± j h) ∂ u¨ A,B (±i h) ∂ 3 u¨ A,B (±i h)2 (± j h) + + + ∂ y2 2! ∂x3 3! ∂ y∂ x 2 3! 2 3 3 3 ∂ u¨ A,B (± j h) ∂ u¨ A,B (±i h)(± j h) + + ··· + 2 ∂y ∂x 3! ∂ y3 3! with i, j = 0, 1, 2. The exact solution u A,B to Eq. (1) meets the following equations: u A±i,B± j = u A,B +

∂ 2 u A,B − c2 ∇ 2 u A,B = 0, ∂t 2

(B.2)

(B.3)

(B.4)

∂ 2i ∂ 2 j ∂ 2 u A,B ∂ 2i ∂ 2 j ∇ 2 u A,B − c2 = 0, (B.5) 2i 2 j 2 ∂ x ∂ y ∂t ∂ x 2i ∂ y 2 j with i, j = 0, 1, 2, 3, 4, . . .. Inserting Eqs. (B.2) - (B.5) into the stencil equation Eq. (B.1) we will get the following local truncation error in space enew for the new approach: ) ( 10 ) [ ( 10 c2 h 10 ∂ u A,B ∂ u A,B ∂ 10 u A,B ∂ 10 u A,B enew = + + (−8223a + 3034a ) + −158(3a2 + a3 ) 2 3 1386000 ∂ x 10 ∂ y 10 ∂ y2∂ x 8 ∂ x 2∂ y8 ( 10 )] ∂ u A,B ∂ 10 u A,B + 7(−3291a2 + 3028a3 + 33000a1 ) + + O(h 12 ). (B.6) 4 6 ∂y ∂x ∂ x 4∂ y6 and econv for the conventional quadratic isogeometric elements: ) ( 6 ∂ 6 u A,B ∂ u A,B + h 6 + O(h 8 ), (B.7) econv = −20c2 ∂x6 ∂ y6 i.e., the new approach improves the order of the local truncation error in space by four orders (similar to the improvement of the order of the numerical dispersion error for the new quadratic isogeometric elements). References [1] A.V. Idesman, Optimal reduction of numerical dispersion for wave propagation problems. part 1: Application to 1-d isogeometric elements, Comput. Methods Appl. Mech. Engrg. 317 (2017) 970–992. [2] T.J.R Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice- Hall, Englewood Cliffs, NJ, 1987. [3] 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 (49–50) (2008) 4104–4124. [4] G. Seriani, S.P. Oliveira, Dispersion analysis of spectral element methods for elastic wave propagation, Wave Motion 45 (6) (2008) 729–744. [5] J.A. Cottrell, A. Reali, Y. Bazilevs, T.J.R. Hughes, Isogeometric analysis of structural vibrations, Comput. Methods Appl. Mech. Engrg. (ISSN: 00457825) 195 (41–43) (2006) 5257–5296. [6] R. Vichnevetsky, Group velocity and reflection phenomena in numerical approximations of hyperbolic equations, J. Franklin Inst. B (ISSN: 00160032) 315 (5–6) (1983) 307–330. [7] L.N. Trefethen, Group velocity in finite difference schemes, SIAM Rev. 24 (2) (1982) 113–136. [8] M. Ainsworth, H.A. Wajid, Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration, SIAM J. Numer. Anal. (ISSN: 00361429) 48 (1) (2010) 346–371. [9] D. Wang, W. Liu, H. Zhang, Novel higher order mass matrices for isogeometric structural vibration analysis, Comput. Methods Appl. Mech. Engrg. 260 (2013) 92–108. [10] D. Wang, W. Liu, H. Zhang, Superconvergent isogeometric free vibration analysis of euler-bernoulli beams and kirchhoff plates with new higher order mass matrices, Comput. Methods Appl. Mech. Engrg. (ISSN: 00457825) 286 (2015) 230–267.

268

A. Idesman, B. Dey / Comput. Methods Appl. Mech. Engrg. 321 (2017) 235–268

[11] A. Idesman, M. Schmidt, J.R. Foley, Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error, Comput. Mech. 47 (2011) 555–572. [12] A.V. Idesman, D. Pham, Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error, Comput. Methods Appl. Mech. Engrg. 271 (2014) 86–108. [13] A.V. Idesman, Accurate time integration of linear elastodynamics problems, Comput. Model. Eng. Sci. 71 (2) (2011) 111–148. [14] A.V. Idesman, H. Samajder, E. Aulisa, P. Seshaiyer, Benchmark problems for wave propagation in elastic materials, Comput. Mech. 43 (6) (2009) 797–814. [15] C.-S. Chou, C. Shu, Y. Xing, Optimal energy conserving local discontinuous galerkin methods for second-order wave equation in heterogeneous media, J. Comput. Phys. 272 (2014) 88–107. [16] A.V. Idesman, D. Pham, J. Foley, M. Schmidt, Accurate solutions of wave propagation problems under impact loading bythe standard, spectral and isogeometric high-order finite elements.comparative study of accuracy of different space-discretization techniques, Finite Elem. Anal. Des. 88 (2014) 67–89. [17] A.-V. Vuong, Ch. Heinrich, B. Simeon, ISOGAT: A 2d tutorial {MATLAB} code for isogeometric analysis, Comput. Aided Geom. Design (ISSN: 0167-8396) 27 (8) (2010) 644– 655. [18] B. Yue, M.N. Guddati, Dispersion-reducing finite elements for transient acoustics, J. Acoust. Soc. Am. 118 (4) (2005) 2132–2141. [19] 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. (ISSN: 00457825) 249–252 (2012) 2–14. [20] L. Beirao Da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L.D. Marini, A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. (ISSN: 02182025) 23 (1) (2013) 199–214. [21] L.B. Da Veiga, F. Brezzi, L.D. Marini, A. Russo, F. Brezzi, G. Manzini, The hitchhiker’s guide to the virtual element method, Math. Models Methods Appl. Sci. (ISSN: 02182025) 24 (8) (2014) 1541–1573.