A Idesman

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

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

A numerical technique with the optimal coefficients of the stencil equation has been suggested. Based on this approach, new high-order isogeometric elements with the reduced dispersion error have been developed for wave propagation problems in the 1-D case. By the modification of the mass and stiffness matrices, the order of the dispersion error is improved from order 2 p (the conventional elements) to order 4 p (the new elements) where p is the order of the polynomial approximations. It was shown that the new approach yields the maximum order of the dispersion error for the stencil equations related to the high-order isogeometric elements. The analysis of the dispersion error of the high-order 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 for 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 o f t he n ew h igh-order i sogeometric e lements w ith r educed d ispersion. We h ave a lso s howed t hat n umerical results obtained by the new and conventional high-order 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 recently developed in our papers. The approach developed in the paper can be directly applied to other space-discretization techniques with similar stencil equations.

In the paper the reduction of the numerical dispersion error for wave propagation problems with the application to isogeometric elements in the 1-D case (Part 1) and the 2-D case (Part 2; see [START_REF] Idesman | Optimal reduction of numerical dispersion for wave propagation problems. Part 2: Application to 2-D isogeometric elements[END_REF]) is considered. Wave propagation in an isotropic homogeneous medium is described by the following scalar wave equation in domain Ω :

∂ 2 u ∂t 2 -c 2 ∇ 2 u = 0, (1)
with the boundary conditions n • ▽u = g 1 on Γ t and u = g 2 on Γ u , and the initial conditions u(x, t = 0) = g 3 , v(x, t = 0) = g 4 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, g i (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 [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF][START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF][START_REF] Seriani | Dispersion analysis of spectral element methods for elastic wave propagation[END_REF] and others) to Eq. ( 1) leads to a system of ordinary differential equations in time

M Ü + c 2 K U = 0, (2)
with

M =  e M e , K =  e K e , (3)
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 M e and K e :

M e =  Ω e N T NdΩ e , (4)
K e =  Ω e  ∂N ∂x  T  ∂N ∂x  dΩ e . (5)
Here, N and ∂N ∂x are the shape matrix and its derivative with respect to the physical coordinate x; see [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF][START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF][START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF]. Due to the space discretization, the exact solution to Eq. ( 2) contains the numerical dispersion error. Usually the analysis of the numerical dispersion error and its improvement for many space-discretization techniques such as the finite elements, spectral elements, isogeometric elements and others start with the analysis and modifications of the elemental mass and stiffness matrices; see [START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF][START_REF] Seriani | Dispersion analysis of spectral element methods for elastic wave propagation[END_REF][START_REF] Ahmadian | Minimization of the discretization error in mass and stiffness formulations by an inverse method[END_REF][START_REF] Ainsworth | Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration[END_REF][START_REF] Cherukuri | Dispersion analysis of numerical approximations to plane wave motions in an isotropic elastic solid[END_REF][START_REF] Dauksher | Solution of elastostatic and elastodynamic problems with Chebyshev spectral finite elements[END_REF][START_REF] Gabriel | Dispersion of elastic waves in the contact-impact problem of a long cylinder[END_REF][START_REF] Guddati | Modified integration rules for reducing dispersion error in finite element method[END_REF][START_REF] Gyrya | M-adaptation method for acoustic wave equation on square meshes[END_REF][START_REF] Krenk | Dispersion-corrected explicit integration of the wave equation[END_REF][START_REF] Marfurt | Accuracy of finite difference and finite element modeling of the scalar and elastic wave equation[END_REF][START_REF] Mullen | Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation[END_REF][START_REF] Seriani | Optimal blended spectral-element operators for acoustic wave modeling[END_REF][START_REF] Yue | Dispersion-reducing finite elements for transient acoustics[END_REF][START_REF] Zyserman | Dispersion analysis of a nonconforming finite element method for the three-dimensional scalar and elastic wave equations[END_REF][START_REF] Li | Seismic scalar wave equation with variable coefficients modeling by a new convolutional differentiator[END_REF][START_REF] He | Dispersion error reduction for acoustic problems using the edge-based smoothed finite element method (ES-FEM)[END_REF][START_REF] Wang | Overview of the discontinuous enrichment method, the ultra-weak variational formulation, and the partition of unity method for acoustic scattering in the medium frequency regime and performance comparisons[END_REF][START_REF] Idesman | Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error[END_REF][START_REF] Idesman | Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error[END_REF][START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF][START_REF] Wang | Superconvergent isogeometric free vibration analysis of Euler-Bernoulli beams and Kirchhoff plates with new higher order mass matrices[END_REF]. For example, one simple and effective finite-element technique for acoustic and elastic wave propagation problems is based on the calculation of the mass matrix M in Eq. ( 2) as a weighted average of the consistent and lumped mass matrices; see [START_REF] Krenk | Dispersion-corrected explicit integration of the wave equation[END_REF][START_REF] Marfurt | Accuracy of finite difference and finite element modeling of the scalar and elastic wave equation[END_REF][START_REF] Mullen | Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation[END_REF][START_REF] Seriani | Optimal blended spectral-element operators for acoustic wave modeling[END_REF] and others. For the 1-D case and the linear finite elements, this approach reduces the error in the wave velocity for harmonic waves from the second order to the fourth order of accuracy. However, for harmonic wave propagation in the 2-D and 3-D cases, these results are not valid (nevertheless, in the multi-dimensional case, the averaged mass matrix yields more accurate results compared with the standard mass matrix).

Some useful techniques for the reduction of the numerical dispersion of the finite element solutions for the Helmholtz equation and for the time-harmonic elastic waves have been developed for a given wave number; e.g., see [START_REF] Harari | Reducing dispersion of linear triangular elements for the Helmholtz equation[END_REF][START_REF] Harari | Stabilized finite elements for time-harmonic elastic waves[END_REF]. In these approaches, the parameters that improve the numerical dispersion depend on the wave number. Therefore, for general wave propagation problems in the time domain that include many other wave numbers, these approaches may not be very efficient.

An interesting technique with implicit and explicit time-integration methods is suggested in [START_REF] Yue | Dispersion-reducing finite elements for transient acoustics[END_REF] for acoustic waves. It is based on the modified integration rule for the calculation of the mass and stiffness matrices for the linear finite elements. In contrast to the averaged mass matrix, the use of the modified integration rule increases the accuracy for the phase velocity from the second order to the fourth order in the general multi-dimensional case of acoustic waves. A similar improvement in the order of the numerical dispersion error for the linear elements in the 1-D case yields the selective mass scaling technique developed in [START_REF] Tkachuk | Direct and sparse construction of consistent inverse mass matrices: general variational formulation and application to selective mass scaling[END_REF]. Using the mimetic finite difference approach, similar results for acoustic waves and electric waves in the 2-D case have been obtained in [START_REF] Gyrya | M-adaptation method for acoustic wave equation on square meshes[END_REF][START_REF] Bokil | Dispersion reducing methods for edge discretizations of the electric vector wave equation[END_REF] for the linear finite elements with reduced dispersion and reduced numerical anisotropy.

The dispersion reduction technique for the high-order finite elements have been suggested in [START_REF] Ainsworth | Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration[END_REF] for acoustic waves. This technique is based on the calculation of the mass matrix M in Eq. (2) as a weighted average of the consistent and lumped mass matrices for the high-order finite elements. It was also shown in [START_REF] Ainsworth | Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration[END_REF] that the same results can be obtained with the modified integration rule for the mass matrix. With this technique, the dispersion error is improved from the order 2 p to the order 2 p + 2 ( p is the order of polynomial approximations). However, even for the quadratic elements, the modified stencil equation derived in [START_REF] Ainsworth | Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration[END_REF] in the 1-D case is based on the replacement of the mid-side degree of freedom within the 1-D element by the degrees of freedom related to the nodes located at the ends of the same element (it seems that such a replacement is impossible in the 2-D or 3-D cases due to the belonging of the same mid-side degree of freedom to different elements).

New isogeometric elements with a higher order of continuity across elements are suggested in [START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF][START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF][START_REF] Auricchio | Isogeometric collocation for elastostatics and explicit dynamics[END_REF][START_REF] Benson | A generalized finite element formulation for arbitrary basis functions: From isogeometric analysis to XFEM[END_REF][START_REF] Benson | A large deformation, rotation-free, isogeometric shell[END_REF] for dynamics problems. It has been shown in these papers that the isogeometric elements yield more accurate numerical results for wave propagation problems compared with the high-order finite elements. The modification of the nondiagonal mass matrix in [START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF][START_REF] Wang | Superconvergent isogeometric free vibration analysis of Euler-Bernoulli beams and Kirchhoff plates with new higher order mass matrices[END_REF] for the isogeometric elements allows the increase in the order of the dispersion error from order 2 p to order 2 p + 2 in the 1-D case and for one specific direction of harmonic waves in the 2-D case (these techniques do not improve the order of the dispersion error in the general 2-D case). It is also necessary to mention that in contrast to the high order finite and spectral elements, the stencil equation for the isogeometric elements on uniform meshes is the same for all internal degrees of freedom (for the high order finite and spectral elements there are several stencil equations with different structures depending on the location of nodes). This simplifies the analysis of the numerical dispersion error for the isogeometric elements. However, along with the advantages, there are some issues with the isogeometric elements. For example, the lumped mass matrix for the isogeometric elements decreases the order of the dispersion error to the second order of accuracy; e.g., see [START_REF] Auricchio | Isogeometric collocation for elastostatics and explicit dynamics[END_REF]. Therefore, a special iterative procedure for linear dynamics problems was developed in [START_REF] Auricchio | Isogeometric collocation for elastostatics and explicit dynamics[END_REF] in order to obtain the same order of the dispersion error for the lumped and consistent mass matrices.

Despite the improved accuracy of the above-mentioned techniques, all of them may yield spurious high-frequency oscillations at high-frequency and impact loadings. These spurious oscillations may lead to divergent results at mesh refinement; e.g., see [START_REF] Idesman | Accurate time integration of linear elastodynamics problems[END_REF][START_REF] Idesman | Accurate finite element simulation and experimental study of elastic wave propagation in a long cylinder under impact loading[END_REF]. The existing techniques for the suppression of the spurious oscillations are based on the introduction of numerical dissipation (or artificial viscosity) for all time increments at the integration of Eq. ( 2); e.g., see [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF][START_REF] Bathe | Finite Element Procedures[END_REF] and many others. However, the existing techniques do not quantify the range of the spurious frequencies and also may yield inaccurate results. Recently a special approach for the accurate solution of acoustic and elastic wave propagation problems has been developed in [START_REF] Idesman | Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error[END_REF][START_REF] Idesman | Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error[END_REF][START_REF] Idesman | Accurate time integration of linear elastodynamics problems[END_REF][START_REF] Idesman | Benchmark problems for wave propagation in elastic materials[END_REF][START_REF] Idesman | A new exact, closed-form a-priori global error estimator for second-and higher-order time integration methods for linear elastodynamics[END_REF][START_REF] Idesman | Accurate finite element modeling of acoustic waves[END_REF][START_REF] Idesman | 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[END_REF][START_REF] Idesman | Accurate finite element modeling of wave propagation in composite and functionally graded materials[END_REF]. This approach includes the two-stage time-integration techniques with the quantification and the filtering of spurious oscillations for different space-discretization methods, an a priori global error estimator in time, the linear finite elements with reduced dispersion. The numerical results in [START_REF] Idesman | Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error[END_REF][START_REF] Idesman | Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error[END_REF] showed that the reduction of the dispersion error for the linear finite elements leads to a huge reduction in the computation time.

The derivation of the analytical formulas for the numerical dispersion error is difficult and very often intractable due to very long and complicated non-linear expressions with a large number of different parameters. This is an issue for high-order space-discretization methods in the multidimensional case. For example, in the recent paper [START_REF] Dede | Isogeometric numerical dispersion analysis for two-dimensional elastic wave propagation[END_REF] the numerical dispersion analysis of the conventional isogeometric elements for elastic waves in the 2-D case was presented numerically (without analytical expressions). To simplify the analytical study and improvement of the numerical dispersion error for high-order space-discretization methods, this paper starts the analysis with the corresponding stencil equation with unknown coefficients. It is important to note that all these unknown coefficients can be found by the minimization of the numerical dispersion error. Therefore, the accuracy of the approach considered in the paper cannot be improved with other coefficients of the stencil equation and is optimal. With the optimal coefficients the order of the numerical dispersion error is increased from order 2 p to order 4 p. In Sections 1.1 and 1.3 we will consider the dispersion analysis of the conventional high-order isogeometric elements in the 1-D case. Based on the dispersion analysis, in Sections 1.2 and 1.4 we will develop new quadratic and cubic isogeometric elements with reduced dispersion in the 1-D case. We will also show that by the modification of the stiffness matrix, the dispersion error of the high-order isogeometric elements with the lumped mass matrix can be improved; see Sections 1.2.2 and 1.4.1. The numerical examples in Section 2 will show the computational efficiency of the new approach. For the derivation of many analytical expressions presented below we have used the computational program "Mathematica".

Dispersion analysis in the 1-D case

Inserting time-harmonic solutions

u(x, t) = exp(i ωt)u(x) (6)
into the wave equation (1) in the 1-D case (∇ 2 u = ∂ 2 u ∂ x 2 ) leads to its reduction to the Helmholtz equation

∂ 2 u ∂ x 2 + k 2 u = 0, ( 7
)
where ω is the angular velocity, k = ω/c is the wave number, i = √ -1. Eq. ( 7) allows the following exact solutions:

u(x) = exp(±ikx). (8)
Fig. 1. The spatial locations of the degrees of freedom u i (i = A -2, A -1, A, A + 1, A + 2) contributing to the stencil equation for the degree of freedom u A of the new and conventional quadratic isogeometric elements. The three elements shown are used for the calculation of the stencil equation for the degree of freedom u A .

After the space discretization, Eq. ( 7) reduces to

(K -k 2 M) U = 0. (9)
For the dispersion analysis presented below we use uniform meshes. For the calculation of the mass and stiffness matrices of the conventional high-order isogeometric elements we use univariate B-spline basis functions defined recursively starting with p = 0

N i, p (ξ ) =  1, if ξ i ≤ ξ < ξ i+1 . 0, otherwise. (10)
For p ≥ 1:

N i, p (ξ ) = ξ -ξ i ξ i+ p -ξ i N i, p-1 (ξ ) + ξ i+ p+1 -ξ ξ i+ p+1 -ξ i+1 N i+1, p-1 (ξ ), (11)
where a knot vector {ξ 1 = 0, . . . , ξ 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 basis functions.

Dispersion analysis for the conventional quadratic isogeometric elements

Let us first start with the dispersion analysis of the conventional quadratic isogeometric elements (the detailed dispersion analysis for the conventional quadratic isogeometric elements can be also found in [START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF][START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF]). The elemental mass M e and stiffness K e matrices of a typical quadratic interior element are

M e = h 120   6 13 1 13 54 13 1 13 6   , K e = 1 6h   2 -1 -1 -1 2 -1 -1 -1 2   , (12)
where h = 1/n el = 1/(Np), n el is the number of elements, N is the number of control points, and p = 2 is the order of the basis functions for the quadratic elements. Considering the discretized Helmholtz equation ( 9) on an infinite line with the sequence of h-spaced control points x A = h A, the stencil equation can be calculated with the help of Eqs. ( 3) and ( 12) for three neighboring elements and has the following form (see Fig. 1 and [START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF]):

k 2 h 2 (u A-2 + 26u A-1 + 66u A + 26u A+1 + u A+2 ) + 20[u A-2 + 2(u A-1 -3u A + u A+1 ) + u A+2 ] = 0. ( 13
)
This equation allows the following solutions (similar to Eq. ( 8)):

u A = exp(±ik h h A), (14)
where k h is the numerical wave number.

Remark. In the dispersion analysis we do not include the solutions of Eq. ( 13) related to the spurious evanescent waves (see [START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF]). The numerical examples in [START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF] show that these solutions do not affect the numerical results.

Inserting Eq. ( 14) into Eq. ( 13) we can find the following relation between the exact and numerical wave numbers k and k h : . Curve 1 corresponds to the conventional quadratic isogeometric elements with the non-diagonal mass matrix. Curves 2 and 3 correspond to the new quadratic isogeometric elements with the 8th and 6th order of the numerical dispersion error (the non-diagonal mass matrix). The slopes of the curves in (b) show the order of convergence at the small mesh size h.

k k h = 1 (k h h) 2 √ 5 √ -2 cos(k h h) -cos(2(k h h)) + 3 √ 26 cos(k h h) + cos(2(k h h)) + 33 . (15)
The numerical dispersion error e = k h -k k = k h k -1 calculated with the help of Eq. ( 15) is plotted in Fig. 2(a); see curve 1. The same results are also plotted in Fig. 2(b) in the logarithmic scale in order to show the order of convergence at mesh refinement. As can be seen, the dispersion error decreases with the decrease in the mesh size k h h. Expanding the right-hand side of Eq. ( 15) into a Taylor series at small h ≪ 1, we get

k k h = 1 + (k h h) 4 1440 + (k h h) 6 6720 + 23(k h h) 8 4147200 + O  (k h h) 10  , (16)
i.e., the conventional quadratic isogeometric elements yield the fourth order of the dispersion error for (k/k h -1); see also Fig. 2(b).

A new approach for the quadratic isogeometric elements with reduced dispersion

For the derivation of new quadratic isogeometric elements with reduced dispersion let us start with the stencil equation with unknown coefficients for the degree of freedom u A that can be written as follows:

k 2 h 2 [m 1 (u A-2 + u A+2 ) + m 2 (u A-1 + u A+1 ) + m 3 u A ] -[k 1 (u A-2 + u A+2 ) + k 2 (u A-1 + u A+1 ) + k 3 u A ] = 0, (17)
where unknown coefficients m j and k j ( j = 1, 2, 3) are calculated in terms of the elemental mass and stiffness matrices, respectively (at this point, these matrices are not defined). For example, comparing Eqs. ( 17) and ( 13), we can find that m 1 = 1, m 2 = 26, m 3 = 66, k 1 = -20, k 2 = -40, k 3 = 120 for the conventional quadratic isogeometric elements. Similar to the conventional isogeometric elements we assumed the symmetry for coefficients m j and k j in the stencil equation [START_REF] Yue | Dispersion-reducing finite elements for transient acoustics[END_REF] for the degrees of freedom u j ( j = A -1, A + 1, A -2, A + 2) symmetrically located with respect to the degree of freedom u A . The idea of the new approach is to find the optimal coefficients of the stencil equation ( 17) that reduce the order of the dispersion error. Inserting Eq. ( 14) into Eq. ( 17) we can find the following relation between the exact and numerical wave numbers k and k h :

k k h = 1 (k h h) √ 2k 1 cos(2(k h h)) + 2k 2 cos(k h h) + k 3 √ 2m 1 cos(2(k h h)) + 2m 2 cos(k h h) + m 3 . ( 18
)
Expanding the right-hand side of Eq. ( 18) into a Taylor series at small h ≪ 1 and analyzing the first terms of a Taylor series we can find that the following constraint k 3 = -2k 1 -2k 2 should be met in order to get at least the first order of the dispersion error. In this case Eq. ( 18) can be rewritten as:

k k h = 1 (k h h) √ 2 √ k 1 cos(2(k h h)) + k 2 cos(k h h) -k 1 -k 2 √ 2m 1 cos(2(k h h)) + 2m 2 cos(k h h) + m 3 . ( 19
)
Expanding again the right-hand side of Eq. ( 19) into a Taylor series at small h ≪ 1 and equating the first term to one and the following three terms to zero, we can find the coefficients m j and k j ( j = 1, 2, 3) from a system of four algebraic equations. They are

m 1 = 23a 1 2358 , m 2 = 344a 1 1179 , m 3 = a 1 , k 1 = - 155a 1 786 , k 2 = - 320a 1 393 , k 3 = 265a 1 131 , (20)
where a 1 is the arbitrary coefficient. The numerical dispersion error e = k h -k k calculated with the help of Eqs. [START_REF] Li | Seismic scalar wave equation with variable coefficients modeling by a new convolutional differentiator[END_REF] and [START_REF] He | Dispersion error reduction for acoustic problems using the edge-based smoothed finite element method (ES-FEM)[END_REF] is plotted in Fig. 2; see curve 2. Expanding the right-hand side of Eq. ( 19) into a Taylor series at small h ≪ 1 with the coefficients given by Eq. ( 20) we get

k k h = 1 - 79(k h h) 8 9525600 - 2633(k h h) 10 2640496320 + O  (k h h) 12  , (21)
i.e., the stencil equation ( 17) with the coefficients given by Eq. ( 20) yields the 8th order of the dispersion error; see also curve 2 in Fig. 2(b).

Remark. The new approach also 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 Appendix.

Let us find the elemental mass and stiffness matrices that yield the stencil equation ( 17) with the coefficients m j and k j ( j = 1, 2, 3) given by Eq. ( 20). Similar to the conventional mass and stiffness matrices, let us assume the following form of the new mass and stiffness matrices:

M e = h   m 11 m 12 m 13 m 12 m 22 m 12 m 13 m 12 m 11   , K e = 1 h   k 11 k 12 k 13 k 12 k 22 k 12 k 13 k 12 k 11   , (22)
where we take into account the symmetry of the coefficients of the mass and stiffness matrices for the degrees of freedom contributing to these matrices. In this case the matrices M e and K e depend on 8 unknown coefficients m 11 , m 12 , m 13 , m 22 and k 11 , k 12 , k 13 , k 22 . With the help of (3), the coefficients of the stencil equation ( 17) can be expressed in terms of the coefficients of the matrices M e and K e as follows:

m 1 = m 13 , m 2 = 2m 12 , m 3 = 2m 11 + m 22 , k 1 = k 13 , k 2 = 2k 12 , k 3 = 2k 11 + k 22 . (23)
Solving simultaneously Eqs. ( 20) and ( 23) we can find the mass and stiffness matrices for the quadratic isogeometric elements with reduced dispersion. They have the following form:

M e = h 120        a 2 6880a 1 393 460a 1 393 6880a 1 393 120a 1 -2a 2 6880a 1 393 460a 1 393 6880a 1 393 a 2        , K e = 1 6h        a 3 - 320a 1 131 - 155a 1 131 - 320a 1 131 1590a 1 131 -2a 3 - 320a 1 131 - 155a 1 131 - 320a 1 131 a 3        , (24)
where a j ( j = 1, 2, 3) are three arbitrary coefficients. The dispersion error is independent of the values of these coefficients a j . Some coefficients a j can be found from the following conditions. Let us assume that the velocity is the same for the entire domain. In this case we can calculate the kinetic energy for one element with the help of the mass matrix and analytically for the considered element. Equating these two expressions we get

V T 0 M e V 0 = v 2 0 h, (25)
where V 0 = (v 0 , v 0 , v 0 ) T is the vector with three equal components. Calculating the left-hand side of Eq. ( 25) with the help of Eq. ( 24) we get that a 1 = 131/210. Next, let us assume that the displacement is the same for the entire domain. In this case strains, stresses and forces are zero. The forces for one element can be calculated with the help of the stiffness matrix:

K e U 0 = 0, (26)
where U 0 = (u 0 , u 0 , u 0 ) T is the vector with three equal components. Calculating the left-hand side of Eq. ( 26) with the help of Eq. ( 24) we get that a 3 = 95/42. For a 1 = 131/210 and a 3 = 95/42, Eq. ( 24) can be rewritten as follows:

M e = h 120       
       . ( 27
)
Remark. All unknown coefficients of the stencil equation ( 17) have been found from the minimization of the order of the dispersion error; see Eq. ( 20). The presence of the arbitrary coefficient a 1 in Eq. ( 20) is explained by the form of the homogeneous discretized Helmholtz equation ( 9); i.e., the multiplication of the mass and stiffness matrices by any coefficient does not change this equation and the stencil equation [START_REF] Yue | Dispersion-reducing finite elements for transient acoustics[END_REF]. The unknown coefficient a 1 can be found from Eq. ( 25). This leads to the following constraint for the coefficients m i of the stencil equation ( 17):

2m 1 + 2m 2 + m 3 = 1. ( 28
)
In this case it follows from Eqs. ( 20) and ( 28) that a 1 = 131/210. Because all coefficients of the stencil equation can be found from the analysis of the dispersion error, the 8th order of the dispersion error in Eq. ( 21) is maximum possible for the considered form of the stencil equation (that is related to the support of basis functions by a specific number of elements) for all quadratic isogeometric elements.

1.2.1. The 6th order of accuracy with the modified mass matrix and the conventional stiffness matrix It is interesting to note that in our papers [START_REF] Idesman | Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error[END_REF][START_REF] Idesman | Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error[END_REF] related to the linear finite elements with reduced dispersion, in the 1-D case we have modified just the mass matrix for the reduction of the dispersion error. Below we will show that it is possible to reduce the order of the dispersion error for the isogeometric elements with the conventional stiffness matrix and the modified mass matrix (partly similar to [START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF]). In order to find the corresponding modifications of the mass matrix, let us consider the stencil equation [START_REF] Yue | Dispersion-reducing finite elements for transient acoustics[END_REF] with k 1 = -20, k 2 = -40, k 3 = 120; see Eq. [START_REF] Krenk | Dispersion-corrected explicit integration of the wave equation[END_REF]. Similar to the previous section, with the help of Eqs. ( 14) and ( 17) we find the dispersion equation, expand it into a Taylor series and equate the first term of a Taylor series to unity, and the following two terms to zero. From these three algebraic equations we can find that m 1 = 7/720, m 2 = 19/90, m 3 = 67/120. The numerical dispersion error e = k h -k k in this case is plotted in Fig. 2; see curve 3. A Taylor series of the dispersion equation with these coefficients is

k k h = 1 + 11(k h h) 6 120960 + (k h h) 8 345600 + O  (k h h) 10  , (29)
i.e., we can get the 6th order of the dispersion error in this case (see also Fig. 2(b)). Finally, using the first three equalities in Eq. ( 23), we can find the following form of the modified mass matrix:

M = h 120        67 -a 1 2 38 3 7 6 38 3 a 1 38 3 7 6 38 3 67 -a 1 2        (30)
with one arbitrary coefficient a 1 . The dispersion error is independent of the value of this coefficient.

Remark. We should note that at a 1 = 54 2 3 the mass matrix in Eq. ( 30) coincides with that in [START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF]. The coefficients m 1 = 7/720, m 2 = 19/90, m 3 = 67/120 of the stencil equation are independent of a 1 and also coincide with those in [START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF]. This means that curves 3 in Fig. 2 correspond to the dispersion error of the isogeometric elements developed in [START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF]. In contrast to [START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF], we use a completely different derivation and also show that the 6th order of the dispersion error in [START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF] for the quadratic isogeometric elements with the conventional stiffness matrix is maximum possible (because the above-mentioned coefficients m i are optimal).

Numerical dispersion of the quadratic isogeometric elements with the lumped mass matrix

The approach presented in the previous sections can be also used for the analysis and the reduction of the dispersion error for the lumped mass matrix. In this case we should consider the stencil equation ( 17) with m 1 = m 2 = 0. For the conventional stiffness matrix with k 1 = 20, k 2 = 40, k 3 = -60, the coefficient m 3 can be found from the analysis of the dispersion equation (as described in the previous sections) and is equal to m 3 = 1. The numerical dispersion error e = k h -k k for the lumped mass matrix and the conventional stiffness matrix is plotted in Fig. 3(a); see curve 1. The same results are also plotted in Fig. 3(b) in the logarithmic scale in order to show the order of convergence at mesh refinement. A Taylor series of the dispersion equation has the following form

k k h = 1 - (k h h) 2 8 + 43(k h h) 4 5760 + O  (k h h) 6  , (31)
i.e., the dispersion error has the second order of accuracy for the lumped mass matrix (see also Fig. 3(b)). Using the first three equalities in Eq. ( 23), we can find the following form of the lumped mass matrix:

M e = h 120      120 -a 1 2 0 0 0 a 1 0 0 0 120 -a 1 2      (32)
with one arbitrary coefficient a 1 . The dispersion error is independent of the value of this coefficient. It is important to note that according to the analysis of the stencil equation with the lumped mass matrix, there are no additional degrees of freedom that allow the reduction of the dispersion error with the conventional stiffness matrix. Independent of the techniques for the calculations of the lumped mass matrix, the best accuracy corresponds to the second order of the dispersion error. However, 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 ( 17) with m 1 = m 2 = 0 and undefined coefficients k j ( j = 1, 2, 3). Using the approach similar to that in Section 1.2 and equating the first term of a Taylor series of the dispersion equation to unity, and the following one term to zero, we can find that

m 3 = 3/4a 1 , k 1 = a 1 /16, k 2 = -a 1 , k 3 = -2k 1 -2k 2 = 15/
8a 1 with one arbitrary coefficient a 1 . The numerical dispersion error e = k h -k k for the lumped mass matrix and the modified stiffness matrix is plotted in Fig. 3; see curve 3. In this case, a Taylor series of the dispersion equation has the following form i.e., we can get the fourth order for the dispersion error (see also Fig. 3(b)). Using Eq. ( 23), we can find the following form of the lumped mass matrix and the modified stiffness matrix:

k k h = 1 - (k h h) 4 180 + (k h h) 6 2016 + O  (k h h) 8  , (33)
M e = h 120     1 2 (15a 1 -a 2 ) 0 0 0 a 2 0 0 0 1 2 (15a 1 -a 2 )     , K e = 1 6h        1 16 (15a 1 -8a 3 ) - a 1 2 a 1 16
-

a 1 2 a 3 - a 1 2 a 1 16 - a 1 2 1 16 (15a 1 -8a 3 )        (34)
with three arbitrary coefficients a j ( j = 1, 2, 3) that do not affect the dispersion error. Some coefficients can be found with the help of Eqs. ( 25) and ( 26) and they are a 1 = a 3 = 8. With these coefficients, Eq. ( 34) can be rewritten as:

M e = h 120     60 - a 2 2 0 0 0 a 2 0 0 0 60 - a 2 2     , K e = 1 6h      7 2 -4 1 2 -4 8 -4 1 2 -4 7 2      . (35)

Dispersion analysis for the conventional cubic isogeometric elements

The analysis of the dispersion error for the conventional cubic isogeometric elements with p = 3 can be done similar to Section 1.1. The mass and stiffness matrices of a typical interior cubic element are

M e = h 5040     20
These matrices lead to the following stencil equation calculated with the help of four elements; see Fig. 4:

k 2 h 2 (u A-3 + 120u A-2 + 1191u A-1 + 2416u A + 1191u A+1 + 120u A+2 + u A+3 ) + 42 (u A-3 + 24u A-2 + 15u A-1 -80u A + 15u A+1 + 24u A+2 + u A+3 ) . (37)
Inserting Eq. ( 14) into Eq. ( 37) we can find the following relation between the exact and numerical wave numbers k and k h :

k k h = 1 (k h h) √ 42 √ -15 cos(k h h) -24 cos(2(k h h)) -cos(3(k h h)) + 40 √ 1191 cos(k h h) + 120 cos(2(k h h)) + cos(3(k h h)) + 1208 . ( 38
)
The numerical dispersion error e = k h -k k = k h k -1 calculated with the help of Eq. ( 38) is plotted in Fig. 5(a); see curve 1. The same results are also plotted in Fig. 5(b) in the logarithmic scale in order to show the order of convergence at mesh refinement. As can be seen, the dispersion error decreases with the decrease in the mesh size k h h. Expanding the right-hand side of Eq. ( 38) into a Taylor series at small h ≪ 1, we get i.e., the conventional cubic isogeometric elements yield the 6th order of the dispersion error for (k/k h -1); see also Fig. 5(b).

k k h = 1 + (k h h) 6 60480 + (k h h) 8 120960 + (k h h) 10 1064448 + 16519(k h h) 12 475517952000 + O  (k h h) 14  , (39)

A new approach for the cubic isogeometric elements with reduced dispersion

The derivation of new cubic isogeometric elements with reduced dispersion is similar to that for the quadratic elements in Section 1.2. Let us start with the stencil equation with unknown coefficients for the degree of freedom u A that can be written as follows:

k 2 h 2 [m 2 (u A-2 + u A+2 ) + m 3 (u A-1 + u A+1 ) + m 4 u A + m 1 (u A-3 + u A+3 )] -[k 2 (u A-2 + u A+2 ) + k 3 (u A-1 + u A+1 ) + k 4 u A + k 1 (u A-3 + u A+3 )] = 0, (40)
where the unknown coefficients m j and k j ( j = 1, 2, 3, 4) are calculated in terms of the elemental mass and stiffness matrices, respectively (at this point, these matrices are not defined). For example, comparing Eqs. [START_REF] Idesman | Accurate finite element modeling of wave propagation in composite and functionally graded materials[END_REF] and [START_REF] Idesman | A new exact, closed-form a-priori global error estimator for second-and higher-order time integration methods for linear elastodynamics[END_REF], we can find that m 1 = 1, m 2 = 120, m 3 = 1191, m 4 = 2416, k 1 = -42, k 2 = -1008, k 3 = -630, k 4 = 3360 for the conventional cubic isogeometric elements. Similar to the conventional isogeometric elements we assumed the symmetry for the coefficients m j and k j in the stencil equation [START_REF] Idesman | Accurate finite element modeling of wave propagation in composite and functionally graded materials[END_REF] for the degrees of freedom u j ( j = A -1, A + 1, A -2, A + 2, A -3, A + 3) symmetrically located with respect to the degree of freedom u A . Inserting Eq. ( 14) into Eq. ( 40) we can find the following relation between the exact and numerical wave numbers k and k h :

k k h = 1 (k h h) √ 2k 2 cos(2(k h h)) + 2k 3 cos(k h h) + k 4 + 2k 1 cos(3(k h h)) √ 2m 2 cos(2(k h h)) + 2m 3 cos(k h h) + m 4 + 2m 1 cos(3(k h h)) . (41)
Expanding the right-hand side of Eq. ( 41) into a Taylor series at small h ≪ 1 and analyzing the first terms of a Taylor series we can find that the following constraint k 4 = -2k 1 -2k 2 -2k 3 should be met in order to get at least the first order of the dispersion error. Expanding the right-hand side of Eq. ( 41) with k 4 = -2k 1 -2k 2 -2k 3 into a Taylor series at small h ≪ 1 and equating the first term to unity and the following five terms to zero, we can find the coefficients m j and k j ( j = 1, 2, 3, 4) from a system of six algebraic equations. They are

m 1 = a 1 , m 2 = 36774a 1 /619, m 3 = 329913a 1 /619, m 4 = 725308a 1 /619, (42)
k 1 = -49483a 1 /1857, k 2 = -261954a 1 /619, k 3 = -263655a 1 /619, k 4 = 3252620a 1 /1857, (43)
where a 1 is the arbitrary coefficient. The numerical dispersion error e = k h -k k calculated with the help of Eqs. ( 41)-( 43) is plotted in Fig. 5; see curve 2. Expanding the right-hand side of Eq. ( 41) into a Taylor series at small h ≪ 1 with the coefficients given by Eqs. [START_REF] Chou | Optimal energy conserving local discontinuous Galerkin methods for second-order wave equation in heterogeneous media[END_REF] and [START_REF] Hulbert | Space-time finite element methods for second-order hyperbolic equations[END_REF] we get

k k h = 1 - 12741(k h h) 12 389701312000 + O  (k h h) 14  . (44)
Eq. ( 44) corresponds to the 12th order of the dispersion error; see also Fig. 5(b).

Remark. The new approach also improves the order of the local truncation error in space by six orders (similar to the improvement of the order of the numerical dispersion error for the new cubic isogeometric elements); see Appendix.

Let us find the elemental mass and stiffness matrices that yield the stencil equation ( 40) with the coefficients m j and k j ( j = 1, 2, 3, 4) given by Eqs. ( 42) and [START_REF] Hulbert | Space-time finite element methods for second-order hyperbolic equations[END_REF]. Similar to the conventional mass and stiffness matrices, let us assume the following form of the new mass and stiffness matrices:

M e = h    
    , (45)
where we take into account the symmetry of the coefficients of the mass and stiffness matrices for the degrees of freedom contributing to these matrices. In this case the matrices M e and K e depend on 12 unknown coefficients m 11 , m 12 , m 13 , m 14 , m 22 , m 23 and k 11 , k 12 , k 13 , k 14 , k 22 , k 23 . With the help of (3), the coefficients of the stencil equation ( 40) can be expressed in terms of the coefficients of the matrices M e and K e as follows:

m 1 = m 14 , m 2 = 2m 13 , m 3 = 2m 12 + m 23 , m 4 = 2(m 11 + m 22 ), (46)
k 1 = k 14 , k 2 = 2k 13 , k 3 = 2k 12 + k 23 , k 4 2(k 11 + k 22 ). (47)
Solving simultaneously Eqs. ( 42), ( 43), ( 46) and (47) we can find the mass and stiffness matrices for the cubic isogeometric elements with reduced dispersion. They have the following form:

M e = h            a 1 329913a 2 1238 - a 3 2
a 3 2 a 1            , (48)
K e = 1 h            1626310a 2 1857 -a 4 -619a 5 -
           , (49)
where a j ( j = 1, 2, 3, 4, 5) are five arbitrary coefficients. The dispersion error is independent of the values of these coefficients. Some coefficients a j can be found from Eqs. ( 25) and ( 26) with U 0 = (u 0 , u 0 , u 0 , u 0 ) T and V 0 = (v 0 , v 0 , v 0 , v 0 ) T where U 0 and V 0 are the vectors with four equal components. Using Eq. ( 25) and the mass matrix from Eq. ( 48) we can find that a 2 = 619/1459920. Using Eq. ( 26) and the stiffness matrix from Eq. ( 49) we can find that a 4 = (25029 -69520a 5 )/139040. For a 2 = 619/1459920 and a 4 = (25029 -69520a 5 )/139040, Eqs. ( 48) and (49) can be rewritten as follows:

M e = h 5040          
          . ( 51
)
Remark. Similar to the explanations in Section 1.2, the 12th order of the dispersion error in Eq. ( 44) is the maximum order for all cubic isogeometric elements.

Numerical dispersion of the cubic isogeometric elements with the lumped mass matrix

The analysis of the numerical dispersion of the cubic isogeometric elements with the lumped mass matrix is similar to that for the quadratic elements in Section 1.2.2. Let us consider the stencil equation ( 40) with m 1 = m 2 = m 3 = 0. For the conventional stiffness matrix with k 1 = -42, k 2 = -1008, k 3 = -630, k 4 = 3360, the coefficient m 4 can be found by the analysis of the dispersion equation (as described in the previous Sections) and is equal to m 4 = 1. In this case, the dispersion equation has the following form

k k h = 1 (k h h) √ -15 cos(k h h) -24 cos(2(k h h)) -cos(3(k h h)) + 40 2 √ 15 . (52)
The numerical dispersion error e = k h -k k given by Eq. ( 52) is plotted in Fig. 3; see curve 2. It is interesting to note that this numerical dispersion error for the cubic elements is greater than that for the quadratic elements with the lumped mass matrix; see curve 1 in Fig. 3. A Taylor series of the right-hand side of Eq. ( 52) yields

k k h = 1 - (k h h) 2 6 + (k h h) 4 80 + O  (k h h) 6  , (53)
i.e., the dispersion error has the second order of accuracy for the lumped mass matrix (see also Fig. 3(b)). Using Eq. ( 46), we can find the following form of the lumped mass matrix:

M e = h         a 1 - a 2 2 0 0 - a 2 2 1 2 -a 1 a 2 0 0 a 2 1 2 -a 1 - a 2 2 0 0 - a 2 2 a 1         (54)
with two arbitrary coefficients a 1 and a 2 . The dispersion error is independent of the values of these coefficients. Due to Eq. ( 3) for the calculation of the global mass matrix, the coefficient a 2 in Eq. ( 54) should be zero in order to get the diagonal global mass matrix (or the elemental mass matrix in Eq. ( 54) should be modified for the elements close to the boundary). It is important to note that according to the analysis of the stencil equation [START_REF] Idesman | Accurate finite element modeling of wave propagation in composite and functionally graded materials[END_REF] with the lumped mass matrix, there are no additional degrees of freedom that allow the reduction of the dispersion error with the conventional stiffness matrix. Independent of the techniques for the calculations of the lumped mass matrix, the best accuracy corresponds to the second order of the dispersion error. However, the accuracy of the numerical results with the lumped mass matrix can be improved by the modification of the stiffness matrix. Let us consider again the stencil equation ( 40) with m 1 = m 2 = m 3 = 0 and undefined coefficients k j ( j = 1, 2, 3, 4). Using the approach similar to that in Section 1.4 and equating the first term of a Taylor series of the dispersion equation to unity, and the following three terms to zero, we can find that m 4 = 90a 1 ,

k 1 = -a 1 , k 2 = 27a 1 /2, k 3 = -135a 1 , k 4 = -2k 1 -2k 2 -2k 3 = 245a 1 with one arbitrary coefficient a 1 . The numerical dispersion error e = k h -k
k for the lumped mass matrix and the modified stiffness matrix is plotted in Fig. 3; see curve 4. In this case, a Taylor series of the dispersion equation has the following form

k k h = 1 - (k h h) 6 1120 + (k h h) 8 7200 + O  (k h h) 10  , (55)
i.e., we can get the 6th order for the dispersion error (see also Fig. 3(b)). Using Eqs. ( 46) and (47), we can find the following form of the lumped mass matrix and the modified stiffness matrix:

M e = h         a 1 - a 2 2 0 0 - a 2 2 -a 1 -45a 3 a 2 0 0 a 2 -a 1 -45a 3 - a 2 2 0 0 - a 2 2 a 1         , ( 56
)
K e = 1 h          1 2 (-2a 4 -245a 3 ) 1 2 (135a 3 -a 5 ) - 27a 3 4 a 3 1 2 (135a 3 -a 5 ) a 4 a 5 - 27a 3 4 - 27a 3 4 a 5 a 4 1 2 (135a 3 -a 5 ) a 3 - 27a 3 4 1 2 (135a 3 -a 5 ) 1 2 (-2a 4 -245a 3 )          (57)
with five arbitrary coefficients a j ( j = 1, 2, 3, 4, 5) that do not affect the dispersion error. The coefficients a 3 = -1/90 and a 4 = 1/40(27 -20a 5 ) can be found with the help of Eqs. ( 25) and [START_REF] Harari | Reducing dispersion of linear triangular elements for the Helmholtz equation[END_REF]. With these coefficients, Eqs. ( 56) and ( 57) can be rewritten as:

M e = h 5040     5040a 1 -2520a 2 0 0 -2520a 2 2520(1 -2a 1 ) 5040a 2 0 0 5040a 2 2520(1 -2a 1 ) -2520a 2 0 0 -2520a 2 5040a 1     , (58)
K e = 1 120h      
      . ( 59
)

Numerical examples

In order to show the advantages of the new high-order isogeometric elements with reduced dispersion, below we solve three 1-D wave propagation problems by the conventional and new high-order isogeometric elements. 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 [START_REF] Idesman | Accurate time integration of linear elastodynamics problems[END_REF][START_REF] Idesman | Benchmark problems for wave propagation in elastic materials[END_REF].

A standing wave in 1-D elastic bar

Similar to [START_REF] Chou | Optimal energy conserving local discontinuous Galerkin methods for second-order wave equation in heterogeneous media[END_REF], the problem of a standing wave is used for the demonstration of the order of convergence of the new numerical technique. Let us consider an elastic bar of length L = 1. The wave velocity is chosen to be c = 1. A standing wave in the 1-D case can be described by the following exact solution to the wave equation: where a = 5 is used. The initial conditions at time t = 0 and the boundary conditions in terms of displacements at x = 0 and x = L are selected according to the exact solution, Eq. ( 60); i.e., the initial displacements are u(x, t = 0) = Sin  aπ x L  , the initial velocities are zero and the two ends of the bar are fixed. The observation time is selected to be T = L/a. According to Eq. ( 60), the displacements at this observation time change the sign and the velocities at all points become zero again. The problem is solved by the conventional and new isogeometric elements on meshes with uniformly spaced control points. For the time integration the trapezoidal rule is 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. 6 shows the convergence of the error in the velocity e v = |v exact (x = L/2, t = T ) -v num (x = L/2, t = T )| in the center of the bar at mesh refinement where v exact (x = L/2, t = T ) and v num (x = L/2, t = T ) are the exact and numerical velocities in the center of the bar at the observation time T . We should mention that the maximum error in the velocity occurs in the center of the bar. h in Fig. 6 is the distance between uniformly spaced control points. At the same h, the meshes with quadratic and cubic isogeometric elements include the same number of degrees of freedom. The results in Fig. 6 are plotted in the logarithmic scale. Therefore, the slopes of the curves at small h shown in Fig. 6 correspond to the order of convergence (the order of accuracy) of the conventional and new isogeometric elements. The new isogeometric elements yield much more accurate results than those obtained by the conventional elements. Moreover, at the same number of degrees of freedom (at the same h), the new quadratic elements are more accurate compared to the conventional quadratic and cubic elements (see curves 1, 2 and 3 in Fig. 6). The results in Fig. 6 are in good agreement with the theoretical order of accuracy of the conventional and new isogeometric elements reported in the previous sections of the paper.

u(x, t) = Sin  aπ x L  Cos  aπt L  , (60)

Propagation of sinusoidal pulse in 1-D elastic bar

Let us consider an elastic bar of length L = 4 with free ends. The wave velocity is chosen to be c = 1. Zero initial displacements u(0, t) = 0 and the following distribution of the initial velocity are applied: v(x, 0) = 0 at x < 1.6 and x > 2 as well as v(x, 0) = Sin(Ω (x -1.6) -π/2) + 1 at 1.6 ≤ x ≤ 2 with frequency Ω = 5π ; see Fig. 7. According to the exact D'Alembert's solution to this problems, the initial velocity sinusoidal pulse splits into two velocity sinusoidal pulses of the same frequency Ω and a half amplitude that propagate in the opposite directions with wave velocity c = 1. For example, the exact solutions at the observation times T = 1 and T = 17 are shown in Fig. 8; see curves 1 (it can be checked that the exact solutions at the observation times T = 1 and T = 17 coincide). Due to the small shift of the center of the initial velocity pulse (located at x = 1.8; see Fig. 7) from the center of the bar the exact solution is not symmetric with respect the center of the bar.

The problem is solved by the conventional and new isogeometric elements on a mesh with uniformly spaced control points with 51, 101 and 201 degrees of freedom. For the time integration we use the trapezoidal rule 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. The numerical results for the velocity distribution at the observation times T = 1 and T = 17 obtained by the conventional quadratic and cubic isogeometric elements and by the new quadratic isogeometric elements with 51 degrees of freedom are ,(c,d) and (e, f) correspond to the numerical results obtained by the conventional quadratic and cubic isogeometric elements and by the new quadratic isogeometric elements with reduced dispersion, respectively (51 degrees of freedom are used for all results). Curves 1 and 2 correspond to the exact and numerical results, respectively.

shown in Fig. 8. As can be seen, despite the same exact solutions at the observation times T = 1 and T = 17, the numerical solutions at these times are different. At the smaller observation time T = 1, the numerical results are close to the exact solutions; see Fig. 8(a,c,e). At the larger observation time T = 17 the numerical results are much less accurate for all elements due to the numerical dispersion error (the space-discretization error is an increasing function of time even at the exact time-integration). The numerical dispersion leads to the shift in space of the center of the sinusoidal pulse as well as to the appearance of small spurious oscillations before and after the main oscillation; compare curves 1 and 2 in Fig. 8(b,d,f). The results also clearly show that at the same numbers of degrees of freedom, the new quadratic isogeometric elements with reduced dispersion significantly improve the accuracy of the numerical results compared with the conventional quadratic and cubic isogeometric elements; see also Fig. 9. Despite a higher accuracy, at the same numbers of degrees of freedom the new quadratic isogeometric elements with reduced dispersion require no additional computational costs compared with the conventional quadratic isogeometric elements and require less computation costs compared with the conventional cubic isogeometric elements (due to a smaller bandwidth of the stencil equation and the corresponding matrix of the system of algebraic equations for the quadratic elements compared to those for the cubic elements). The convergence of the numerical results obtained by the new and conventional quadratic isogeometric elements at mesh refinement is shown in Fig. 10. As can be seen from Fig. 10, the new quadratic isogeometric elements converge to the exact solution much faster compared to the conventional quadratic isogeometric elements. Fig. 11 shows that at the larger observation time T = 17, the new cubic isogeometric elements with reduced dispersion are much more accurate compared to the conventional cubic isogeometric elements (see curves 2 Remark. 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. 12; see curves 1. As can be seen, in this case the field function u(x, t) is very close to that obtained by the new isogeometric elements with the conventional isogeometric basis functions; see curves 2 in Fig. 12.

1-D elastic bar under suddenly applied load

The 1-D impact problem considered below is a good benchmark wave propagation problem; e.g., see [START_REF] Idesman | Benchmark problems for wave propagation in elastic materials[END_REF]. All low and high frequencies are excited for this problem. Despite a simple analytical solution, this problem is challenging for any numerical approach due to the appearance of spurious oscillations in numerical solutions; e.g., see [START_REF] Idesman | Accurate time integration of linear elastodynamics problems[END_REF][START_REF] Hulbert | Space-time finite element methods for second-order hyperbolic equations[END_REF] and Figs. 14(a,b) and 15(a, b). To treat this issue, we use the two-stage time integration technique developed in our papers [START_REF] Idesman | Accurate time integration of linear elastodynamics problems[END_REF][START_REF] Idesman | Benchmark problems for wave propagation in elastic materials[END_REF][START_REF] Idesman | 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[END_REF]. The filtering stage of this technique includes the time continuous Galerkin (TCG) method with 10 time increments (5 positive plus 5 negative time increments) the size of which is calculated according to the special formulas; see [START_REF] Idesman | Accurate time integration of linear elastodynamics problems[END_REF][START_REF] Idesman | Benchmark problems for wave propagation in elastic materials[END_REF]. We should also mention that due to a simple analytical solution, it is easy to compare the accuracy of the numerical solutions of this problem when the spurious oscillations are removed from the results. The numerical solutions of this problem obtained by the conventional high-order finite elements, spectral elements and isogeometric elements have been compared in our paper [START_REF] Idesman | 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[END_REF]. The numerical results in [START_REF] Idesman | 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[END_REF] also show the computational efficiency of high-order elements compared with the linear elements at the increase in the observation time (the numerical dispersion significantly affects the error increase at large observation times even for the exact time integration). Let us consider an elastic bar of length L = 4 under a suddenly applied load (see Fig. 13). The wave velocity is chosen to be c = 1. The following boundary conditions are applied: u(0, t) = t (which corresponds to the velocity v(0, t) = v 0 = 1) and u(4, t) = 0 (which corresponds to the velocity v(4, t) = 0). Zero initial conditions are used; i.e., u(x, 0) = v(x, 0) = 0. The observation time is chosen to be T = 18. During this time the waves travel within the bar with two reflections from each end of the bar. The analytical solution to this problem for the velocity at T = 18 is: v a (x, 18) = 1 for 0 ≤ x < 2 and v a (x, 18) = 0 for 2 < x ≤ 4 (at x = 2, a jump in the velocity occurs).

The problem is solved by the conventional and new quadratic and cubic isogeometric elements on a mesh with uniformly spaced control points with 101 degrees of freedom. For the time integration at the stage of basic computations we use the trapezoidal rule 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. The details related to the filtering stage can be found in our papers [START_REF] Idesman | Accurate time integration of linear elastodynamics problems[END_REF][START_REF] Idesman | Benchmark problems for wave propagation in elastic materials[END_REF][START_REF] Idesman | 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[END_REF]. The numerical results for the velocity distribution at time T = 18 obtained by the conventional and new quadratic and cubic isogeometric elements are shown in Figs. 14 and15. As can be seen, after basic computations all results include spurious high-frequency oscillations; see Figs. 14(a,b) and 15(a, b). These oscillations are smaller for the isogeometric elements with the smaller dispersion error. However, 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 simple due to the simple analytical solution of the 1-D impact problems (curves 1 in Figs. 14 and15). We can clearly see that the new high-order isogeometric elements with reduced dispersion significantly improve the accuracy of the numerical results compared with the conventional high-order isogeometric elements; compare curves 2 and 3 in Fig. 14(c,d) and curves 2 and 3, 4 in Fig. 15(c,d). We should also mention that the new quadratic isogeometric elements with reduced dispersion yield more accurate results than those obtained by the conventional cubic isogeometric elements; see curves 2 and 3 in Fig. 15(c,d). Despite a higher accuracy, the new high-order isogeometric elements with reduced dispersion do not require additional computational costs compared with the conventional high-order isogeometric elements.

Remark. It was shown in [START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF][START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] that uniformly spaced control points yield more accurate results for structural dynamics problems than a linear parametrization with the constant Jacobian determinant. However, even with uniformly spaced control points, the numerical results include spurious oscillations in basic computations (see Figs. 14(a,b) and 15(a, b)) and require the filtering stage.

Concluding remarks

The Part 1 of paper considers the development of the new high-order isogeometric elements with reduced dispersion in the 1-D case. The new findings of the paper can be summarized as follows:

• By the analysis of the dispersion error of the stencil equation with arbitrary coefficients we have shown that these coefficients can be found by the minimization of the order of the dispersion error. In this case, the order of the dispersion error can be increased from the order 2 p for the conventional high-order isogeometric elements 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 (that is related to the support of basis functions by a specific number of elements) for the corresponding high-order isogeometric elements. We should mention that we have never seen in the literature such a significant improvement of the order of the dispersion error for the conventional finite elements, the spectral elements or the isogeometric elements. For example, for the linear and high-order finite and isogeometric elements (e.g., see [START_REF] Ainsworth | Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration[END_REF][START_REF] Yue | Dispersion-reducing finite elements for transient acoustics[END_REF][START_REF] Idesman | Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error[END_REF][START_REF] Idesman | Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error[END_REF][START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF]), the order of the dispersion error has been improved from the order 2 p to the order 2 p + 2 and this leads to a significant decrease in the computation time at a given accuracy. The increase in accuracy for the new isogeometric elements with the order 4 p of the dispersion error is much higher compared with the known approaches. This will lead to a much more significant reduction in the computation time at a given accuracy. • Many wave propagation problems require small time increments for the time integration and the use of explicit time-integration methods with the lumped mass matrix is very computationally efficient. However, for the highorder isogeometric elements, the lumped mass matrix yields the second order of the dispersion error (e.g., see [START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF][START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF]) and makes the application of the high-order isogeometric elements computationally inefficient. Therefore, a special iterative procedure for linear dynamics problems that requires additional computational costs was suggested in [START_REF] Auricchio | Isogeometric collocation for elastostatics and explicit dynamics[END_REF] in order to obtain the same order of the dispersion error for the lumped and consistent mass matrices. In the paper we have shown that independent of the procedures for the calculation of the lumped mass matrix, the order of the dispersion error cannot not be improved with the conventional stiffness matrix and no additional iterations. We also showed 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). • For the linear basis functions, the conventional isogeometric elements coincide with the linear finite elements. The analysis of the stencil equation for the linear elements shows that the reduction of the numerical dispersion error for the linear elements from order two to order four in the known numerical techniques (e.g., see [START_REF] Krenk | Dispersion-corrected explicit integration of the wave equation[END_REF][START_REF] Marfurt | Accuracy of finite difference and finite element modeling of the scalar and elastic wave equation[END_REF][START_REF] Yue | Dispersion-reducing finite elements for transient acoustics[END_REF][START_REF] Idesman | Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error[END_REF]) is maximum possible for all linear elements and cannot be improved (the dispersion analysis of the linear elements is similar to that for the quadratic and cubic elements presented in Section 1). • By the solution of wave propagation problems in Section 2 we have shown the computational efficiency of the new high-order isogeometric elements with reduced dispersion. For example, the new quadratic isogeometric elements yield much more accurate results than those obtained by the conventional cubic isogeometric elements; e.g., see Figs. 6, 9 and 15. • Despite a good accuracy, the conventional and new high-order isogeometric elements yield spurious oscillations under high-frequency and impact loadings. In the paper we have shown that by the application of the two-stage time-integration technique developed in our recent papers [START_REF] Idesman | Accurate time integration of linear elastodynamics problems[END_REF][START_REF] Idesman | Benchmark problems for wave propagation in elastic materials[END_REF][START_REF] Idesman | 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[END_REF] we can get accurate numerical results with the isogeometric elements as well; e.g., see Figs. 14 and15. This technique quantifies and filters spurious oscillations. • We should also mention that the approach developed in the paper can be directly applied to other spacediscretization techniques that have the same form of the stencil equations (e.g., see Eqs. ( 17) and ( 40)); e.g., to the high-order finite-difference method in order to reduce the dispersion error.

The development of the new high-order isogeometric elements with reduced dispersion for wave propagation problems in the 2-D case with the corresponding numerical examples is considered in Part 2 of the paper; see [START_REF] Idesman | Optimal reduction of numerical dispersion for wave propagation problems. Part 2: Application to 2-D isogeometric elements[END_REF]. is similar to that for Eq. ( 17)):

h 2 [m 1 ( ü A-2 + ü A+2 ) + m 2 ( ü A-1 + ü A+1 ) + m 3 ü A ] + c 2 [k 1 (u A-2 + u A+2 ) + k 2 (u A-1 + u A+1 ) + k 3 u A ] = 0, (A.1)
where the coefficients m j and k j ( j = 1, 2, 3) are given by Eq. ( 20) for the new elements with reduced dispersion and are m 1 = 1, m 2 = 26, m 3 = 66, k 1 = -20, k 2 = -40, k 3 = 120 for the conventional quadratic isogeometric elements; see Section 1.1. Similarly, the stencil equation for the degree of freedom u A in the time domain for the new and conventional cubic elements has the following form (the derivation is similar to that for Eq. ( 40)):

h 2 [m 2 ( ü A-2 + ü A+2 ) + m 3 ( ü A-1 + ü A+1 ) + m 4 ü A + m 1 ( ü A-3 + ü A+3 )] + c 2 [k 2 (u A-2 + u A+2 ) + k 3 (u A-1 + u A+1 ) + k 4 u A + k 1 (u A-3 + u A+3 )] = 0, (A.2)
where the coefficients m j and k j ( j = 1, 2, 3, 4) are given by Eqs. ( 42) and ( 43 Because at the derivation of the coefficients m i and k i for the new approach (see Eq. ( 20) and Eqs. ( 42) and ( 43)) the basis functions have not been used we assume that the control variables u i 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 and Fig. 12 in Section 2.2. Under this assumption, let us expand u i (i = A-3, A-2, A-1, A+1, A+2, A+3) into a Taylor series at small h ≪ 1 as follows:

u A±l = u A + ∂u A ∂ x (±lh) + ∂ 2 u A ∂ x 2 (±lh) 2 2! + ∂ 3 u A ∂ x 3 (±lh) 3 3! + ∂ 4 u A ∂ x 4
(±lh) 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). Similarly, inserting Eqs. (A.3)-(A.5) into the stencil equation (A.2) we will get the following local truncation error in space e new for the new cubic isogeometric elements: i.e., the new approach improves the order of the local truncation error in space by six orders (similar to the improvement of the order of the numerical dispersion error for the new cubic isogeometric elements).

e new = 38223c

Fig. 2 .

 2 Fig. 2. The numerical dispersion errors e = k h -k k (a) and ln(|e|) (b) of the quadratic isogeometric elements as a function of the mesh size k h h (a) and ln(k h h) (b). Curve 1 corresponds to the conventional quadratic isogeometric elements with the non-diagonal mass matrix. Curves 2 and 3 correspond to the new quadratic isogeometric elements with the 8th and 6th order of the numerical dispersion error (the non-diagonal mass matrix). The slopes of the curves in (b) show the order of convergence at the small mesh size h.

Fig. 3 .

 3 Fig. 3. The numerical dispersion errors e = k h -k k (a) and ln(|e|) (b) of the quadratic and cubic isogeometric elements with the lumped mass matrix as a function of the mesh size k h h (a) and ln(k h h) (b). Curves 1 and 2 correspond to the conventional quadratic and cubic isogeometric elements with the lumped mass matrix. Curves 3 and 4 correspond to the new quadratic and cubic isogeometric elements with the lumped mass matrix. The slopes of the curves in (b) show the order of convergence at the small mesh size h.

Fig. 4 .

 4 Fig.4. The spatial locations of the degrees of freedom u i (i = A -3, A -2, A -1, A, A + 1, A + 2, A + 3) contributing to the stencil equation for the degree of freedom u A of the new and conventional cubic isogeometric elements. The four elements shown are used for the calculation of the stencil equation for the degree of freedom u A .

Fig. 5 .

 5 Fig. 5. The numerical dispersion errors e = k h -k k (a) and ln(|e|) (b) of the cubic isogeometric elements as a function of the mesh size k h h (a) and ln(k h h) (b). Curves 1 and 2 correspond to the conventional and new cubic isogeometric elements with the non-diagonal mass matrix. The slopes of the curves in (b) show the order of convergence at the small mesh size h.

Fig. 6 .

 6 Fig. 6. The error in velocity e v = |v exact (x = L/2, t = T ) -v num (x = L/2, t = T )| as a function of the mesh size h in the logarithmic scale. Curves 1 and 2 correspond to the conventional quadratic and cubic isogeometric elements. Curves 3 and 4 correspond to the new quadratic and cubic isogeometric elements. The slopes of the curves at small h show the order of convergence of the corresponding techniques.

Fig. 7 .

 7 Fig. 7. The distribution of the initial velocity in a 1-D bar of length L = 4.

Fig. 8 .

 8 Fig. 8. The velocity distribution in the 1-D bar at the observation times T = 1 (a, c, e) and T = 17 (b, d, f). (a, b),(c, d) and (e, f) correspond to the numerical results obtained by the conventional quadratic and cubic isogeometric elements and by the new quadratic isogeometric elements with reduced dispersion, respectively (51 degrees of freedom are used for all results). Curves 1 and 2 correspond to the exact and numerical results, respectively.

Fig. 9 .

 9 Fig. 9. The comparison of the velocity distribution in the 1-D bar at the observation time T = 17. Curves 1, 2, 3 and 4 correspond to the exact solution and the numerical results obtained by the conventional quadratic and cubic isogeometric elements and by the new quadratic isogeometric elements with reduced dispersion, respectively (51 degrees of freedom are used for all results); see also Fig. 8.

  ) and more accurate compared to the new quadratic elements at the same number of degrees of freedom (see curves 2 in Figs. 10(b) and 11(b)).

Fig. 10 .

 10 Fig. 10. The velocity distribution in the 1-D bar at the observation time T = 17 obtained by the conventional quadratic isogeometric elements (a, c, e) and by the new quadratic isogeometric elements with reduced dispersion (b, d, f). (a, b), (c, d) and (e, f) correspond to the numerical results obtained on the meshes with 51, 101 and 201 degrees of freedom, respectively. Curves 1 and 2 correspond to the exact and numerical results, respectively.

Fig. 11 .

 11 Fig. 11. The velocity distribution in the 1-D bar at the observation time T = 17 obtained by the conventional cubic isogeometric elements (a) and by the new cubic isogeometric elements with reduced dispersion (b). 51 degrees of freedom are used for all results. Curves 1 and 2 correspond to the exact and numerical results, respectively.

Fig. 12 .

 12 Fig. 12. The velocity distribution in the 1-D bar at the observation time T = 17 obtained by the solution of the global system of semidiscrete equations with 51 (a) and 101 (b) degrees of freedom (the same as for the new quadratic isogeometric elements with reduced dispersion). Curves 1 correspond to the post-processing similar to that for the finite elements (see the text for more explanations). Curves 2 correspond to the new quadratic isogeometric elements with reduced dispersion.

Fig. 13 .

 13 Fig. 13. Impact of an elastic bar of length L = 4 against a rigid wall.

Fig. 14 .

 14 Fig. 14. The velocity distribution along the bar for the 1-D impact problem at time T = 18 after several reflections from the ends of the bar (a-d). Curve 1 corresponds to the analytical solution. Curves 2 and 3 correspond to the numerical results without (a, b) and with (c, d) the filtering and are obtained with the conventional quadratic isogeometric elements and the new quadratic isogeometric elements with reduced dispersion, respectively (101 degrees of freedom are used for all results). (b, d) zoom (a, c) in the vicinity of the wave front (1.8 < x < 2.2).

Fig. 15 .

 15 Fig. 15. The velocity distribution along the bar for the 1-D impact problem at time T = 18 after several reflections from the ends of the bar (a-d). Curve 1 corresponds to the analytical solution. Curves 2, 3 and 4 correspond to the numerical results without (a, b) and with (c, d) the filtering and are obtained with the conventional cubic isogeometric elements and the new quadratic and cubic isogeometric elements with reduced dispersion, respectively (101 degrees of freedom are used for all results). (b, d) zoom (a, c) in the vicinity of the wave front (1.8 < x < 2.2).

  ) for the new elements with reduced dispersion and are m 1 = 1, m 2 = 120, m 3 = 1191, m 4 = 2416, k 1 = -42, k 2 = -1008, k 3 = -630, k 4 = 3360 for the conventional cubic isogeometric elements; see Section 1.3. Let us analyze the local truncation error in space for the stencil equations (A.1) and (A.2) with the corresponding coefficients m i and k i . The local truncation error is the residual of the stencil equation when the unknowns u i (i = A -3, A -2, A -1, A, A + 1, A + 2, A + 3) correspond to the exact solution.

4

 4

6 ∂ 6 u A ∂ x 6 h 6 +

 66 l = 1, 2 for the quadratic isogeometric elements and l = 1, 2, 3 for the cubic isogeometric elements. The exact solution u A to Eq. (1) meets the following equations:∂ 2 u A ∂t 2 -c 2 ∂ 2 u A ∂ x 2 = 0, (A.4) ∂ 2 j+2 u A ∂ x 2 j ∂t 2 -c 2 ∂ 2 j+2 u A ∂ x 2 j+2 = 0, (A.5) with j = 1,2, 3, 4, . . .. Inserting Eqs. (A.3)-(A.5) into the stencil equation (A.1) we will get the following local truncation error in space e new for the new quadratic isogeometric elements: e new = 79c 2 a 1 2971080 ∂ 10 u A ∂ x 10 h 10 + O(h 12 ) (A.6) and e conv for the conventional quadratic isogeometric elements: e conv = -c 2 O(h 8 ), (A.7)

  m 11 m 12 m 13 m 14 m 12 m 22 m 23 m 13 m 13 m 23 m 22 m 12 m 14 m 13 m 12 m 11 12 k 13 k 14 k 12 k 22 k 23 k 13 k 13 k 23 k 22 k 12 k 14 k 13 k 12 k 11

					 k 11 k
	   ,	K e	=	1 h	  

  2 a 1 247847600∂ 14 u A ∂ x 14 h 14 + O(h 16 ) (A.8)and e conv for the conventional cubic isogeometric elements:

	e conv =	6 c 2	∂ 8 u

A ∂ x 8 h 8

+ O(h 10 ), (A.9)

Acknowledgments

The research has been supported in part by the Air Force Office of Scientific Research (contract FA9550-16-1-0177) and by Texas Tech University.

Appendix. Local truncation error in space for the new and conventional 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 in the time domain for the new and conventional quadratic elements has the following form (the derivation


Mohammad Irfan Alam

Mohammad Irfan

MULTIDISCIPLINARY DESIGN OPTIMIZATION OF STRATOSPHERIC AIRSHIP

Keywords: High Altitude Airship, Sensitivity Analyses, MDO, Genetic Algorithm

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

There is a global interest in design and development of stratospheric airships [START_REF] Epley | Stratospheric aircraft, blimps, balloons and long endurance vehicles[END_REF], which can serve as a long endurance platform for deployment of equipment for several commercial and strategic applications e.g., next generation wireless broadband telecommunications [START_REF] Tozer | High-altitude platforms for wireless communications[END_REF], digital broadcasting [START_REF] Grace | Broadband communications via high-altitude platforms[END_REF], coastal surveillance [START_REF] Colozza | High-altitude, longendurance airships for coastal surveillance[END_REF], remote sensing and GPS augmented navigation systems [START_REF] Tsujii | A navigation/positioning service based on pseudolites installed on stratospheric airships[END_REF]. These airships are designed to be able to maintain a quasi-stationary position at altitudes of around 20 km, where ambient winds are of low magnitude. Such airships function as low-altitude satellites, but offer much shorter transmission distances and ranges with high resolution, and lesser signal propagation errors.

They are much more economical compared to satellites, as they can be relocated or brought down and refurbished with latest equipment.

The power requirements for an airship are two-fold, to meet the needs of the on-board mounted equipment, as well as the propulsive power needed to either maintain its station in the presence of ambient winds, or to travel at a required speed when it is to be relocated. Since such airships are to be deployed for long durations (several weeks or months), the only practical mechanism to address their power needs is a Solar Regenerative Fuel Cell (SRFC) system. In such a system, adequate amount of solar arrays are mounted at the top of envelope of the airship. During day time, the solar arrays generate enough power to meet the needs, and any excess power is used to charge the onboard batteries, which meet the power requirements at night time, or during lean periods. The working principle of a typical SRFC is depicted in Fig. 1.

Fig. 1. Working Principle of SRFC

Several

researchers have proposed methodologies and approaches for conceptual design and sizing of stratospheric airships [START_REF] Lobbia | A modular sizing model for high-altitude/long-endurance airships[END_REF][START_REF] Yu | Configurations analysis for highaltitude/long-endurance airships[END_REF][START_REF] Chen | Analysis to effects on conceptual parameters of stratospheric airship with specified factors[END_REF][START_REF] Zhang | Shape optimization of airship based on constrained particle swarm optimization[END_REF]. The shape of the envelope is one of the most critical elements in the design of such systems, and envelope shape optimization is a key area of research in this field [START_REF] Zhang | Shape optimization of airship based on constrained particle swarm optimization[END_REF]. Concurrent subspace optimization techniques have also been applied to the conceptual design and sizing of airships [START_REF] Liang | Conceptual design optimization of high altitude airship in concurrent subspace optimization[END_REF]. A critical review of all these methodologies has been carried out by Alam and Pant [START_REF] Alam | A methodology for sizing and optimization of high altitude airship[END_REF].

There is a need to develop a simple yet accurate and more realistic multidisciplinary design methodology, keeping in mind the operating scenario and meteorological conditions. For example, high winds lead to steep rise in drag, which results in higher power requirements, i.e., larger solar arrays. At the same time high wind can also maintain the efficiency of solar cells by restricting temperature rise. Therefore, a multidisciplinary design and analyses approach is needed to locate a truly optimal design.

This paper describes such a methodology, which can be used for initial sizing of stratospheric airship, to meet the user-specified requirements of payload weight and power.

A composite objective function is devised which takes into account various factors that influence airship performance, including aerodynamics, structures, energy and weight. The envelope shape is parameterized in terms of some geometry related parameters; and the optimum shape that minimizes this composite objective function is obtained. Constraints are imposed on the volume of the airship, to ensure comparability of the design results. Optimal solutions are obtained using an evolutionary technique, viz., Genetic Algorithm

System modeling for analyses and optimization

In the past, drag was often considered in airship optimization. When the coefficient of drag was smaller, the airship shape was considered more optimal [START_REF] Liang | Conceptual design optimization of high altitude airship in concurrent subspace optimization[END_REF]. In fact, the practical studies gradually show that many factors restrict the airship shape and resistance is only one of the several important factors, which also include weight, structural strength etc. After screening and contrasting to determine the most important factors on airship, coefficient of drag (CDV), surface area of the airship envelope (Ae), minimum hoop's stress (σmin) and area of the solar array (Asa) are employed in the present study.

Gertler Series-58 Shape Generator

For design, analyses and optimization of airship, Gertler Series 58 Shape Generator ( [START_REF] Gertler | Resistance Experiments on a Systematic Series of Streamlined Bodies of Revolution For Application to the Design of High-Speed Submarines[END_REF], [START_REF] Landweber | Mathematical formulation of bodies of revolution[END_REF]) (Eq. 1) has been considered, which is parameterized in terms of some geometry related parameters; and the optimum shape that minimizes composite objective function is obtained.

𝑦 2 (𝑥) = 𝐷 2 [𝑎 1 ( 𝑥 𝐿 ) + 𝑎 2 ( 𝑥 𝐿 ) 2 + 𝑎 3 ( 𝑥 𝐿 ) 3 … + 𝑎 6 ( 𝑥 𝐿 ) 6 ] (1)
Where, a1, a2, a3....a6 are shape coefficients whose values are driven by geometrical parameters of airship viz. m (location of maximum diameters), ro (Nose radius), r1 (tail radius) and Cp (prismatic coefficients).

Model of surface area

The surface area of the envelope 'Ae' can be calculated using Eq. ( 2) as:

𝐴 𝑒 = 2𝜋 ∫ 𝑦√ 1 + ( 𝑑𝑦 𝑑𝑥 ) 2 𝑑𝑥 𝐿 0 (2)

Model of volumetric drag coefficient

In order to seek the drag of the airship envelope, volumetric drag coefficient (CDV) is calculated as per the formula quoted by Cheeseman in [START_REF] Cheeseman | Airship Technology[END_REF].

𝐶 𝐷𝑉 = 0.172( 𝐿 𝐷 ) 1 3 +0.252( 𝐷 𝐿 ) 1.2 +1.032( 𝐷 𝐿 ) 2.7 𝑅𝑒 1 6 (3)
Where, L is the envelope length, D is the maximum diameter of the envelope and Re is the Reynolds number.

Model of Minimum Hoop Stress

Several researchers have estimated the circumferential hoop stress (σhoop) by assuming the airship body to approximate a thin cylinder with hemi-spherical ends. In this study, the pressures are obtained by considering the Elastic Engineering theory. A generalized moment equation is derived based on the Elastic theory with an assumption of linear distribution of mass along the length of the airship.

In order to maintain positive internal pressure, the minimum inner pressure (Δp) consisting of static pressure (pstatic), Munk pressure (pdyn) and internal differential pressure (pdiff) are calculated using Eqns. ( 4) to [START_REF] Yu | Configurations analysis for highaltitude/long-endurance airships[END_REF] as:

𝛥𝑝 = 𝑝 𝑠𝑡𝑎𝑡𝑖𝑐 + 𝑝 𝑑𝑦𝑛 + 𝑝 𝑑𝑖𝑓𝑓 (4)
The static pressure is caused by static bending moments, which is kept 15% higher to avoid nose caving.

𝑝 𝑠𝑡𝑎𝑡𝑖𝑐 = 1.15 ( 1 2 𝜌 𝑎 𝑣 2 ) (5)
Munk pressure caused by dynamic bending moment for a given pressure coefficient (cp) can be expressed as:

𝑝 𝑑𝑦𝑛 = 𝑐 𝑝 ( 1 2 𝜌 𝑎 𝑣 2 ) (6)
Differential pressure due to lifting gas and outside air is given by:

𝑃 𝑑𝑖𝑓𝑓 = (𝜌 𝑎-𝜌 𝑔 )𝑔 𝐷 2 (7)
The hoops stress can be calculated using Eq. ( 8) as:

𝜎 ℎ𝑜𝑜𝑝 = ∆𝑝 𝐷 2 (8)

Model of Power requirement

Total power requirement (Ptotal) of an airship consists of payload systems (Ppay), control systems (Pctrl) and propulsion systems (Pthrust).

The payload power and the control systems power are assumed to be constant at 10 W and 11 kW, respectively.

𝑃 𝑡𝑜𝑡𝑎𝑙 = 𝑃 𝑃𝑎𝑦 + 𝑃 𝑐𝑡𝑟𝑙 + 𝑃 𝑡ℎ𝑟𝑢𝑠𝑡 (9)
The energy used for propulsion system can be calculated as,

𝑃 𝑡ℎ𝑟𝑢𝑠𝑡 = 𝐷 𝑡𝑜𝑡𝑎𝑙 𝑣 𝜂 𝑝𝑟𝑜𝑝 (10)
Drag 'D' acting on the airship can be obtained using Eq. ( 11) as:

𝐷 𝑡𝑜𝑡𝑎𝑙 = 1 2 𝑁𝜌 𝑎 𝑣 2 𝐶 𝐷𝑉 𝑉 2 3 (11)
Where,

𝑉 = 𝜋 ∫ 𝑦 2 𝑑𝑥 𝑙 0 (12)
The total energy Qtotal is estimated using Eq. ( 13) as:

Q 𝑡𝑜𝑡𝑎𝑙 = 𝑃 𝑡𝑜𝑡𝑎𝑙 𝑡 𝑑𝑎𝑦 + 𝑃 𝑡𝑜𝑡𝑎𝑙 𝑡 𝑑𝑎𝑦 𝜂 𝑐𝑜𝑛𝑣 (13)
Where, tday and tnight are the duration of day and night time respectively and ƞconv is the conversion efficiency for storing the energy.

Model of Solar Area requirement

In the present study an approach similar to Wang et al. [START_REF] Wang | A methodology for optimization design and analysis of stratosphere airship[END_REF] is adopted. In which, for simplicity the airship geometry is assumed to approximate a cylinder. The radius of the envelope at the mass center rgc is taken as the radius of the cylinder.

So, area of solar array (Asa) required can be estimated by using Eq. ( 14) as:

𝐴 𝑠𝑎 = ∬ 𝑦(𝑥) (1 + ( 𝑑𝑦 𝑑𝑦 ) 2 ) 1/2 𝑑𝑥𝑑𝜉 (14)
Where, dx is elementary length along length and dξ is angular width of elementary area of solar panel along circumferential direction.

To find the total length of solar panel XSC and included angle ξ, as a limit of Eq. ( 15), energy supplied (Qsup) can be equated to total energy required (Qtotal) as:

𝑄 𝑠𝑢𝑝 = 𝑄 𝑡𝑜𝑡𝑎𝑙 (15)
Where, Qsup can be calculated by Eq. ( 16) as:

𝑄 𝑠𝑢𝑝 = ∭ ƞ 𝑆𝐶 𝐼 𝑑 𝑐𝑜𝑠𝜃𝑟 𝑔𝑐 𝑑𝑥𝑑ξdt (16)
Where, ƞ 𝑆𝐶 and 𝐼 𝑑 are solar conversion efficiency and normal solar irradiance values respectively. And 𝜃 is angle between surface normal vector and 𝐼 𝑑 . Eq. ( 16) was integrated for the whole day time (tday) to get total power supplied.

Wind Model

NASA Horizontal Wind Model (HWM07) was used to plot the wind conditions at different geographic locations. Fig. 2 illustrates the variation of wind over the year at the height of 20 km.

Fig. 2. Variation of wind over the year

The variations of wind speed with altitude for baseline cases are shown in Fig. 3. The figure illustrates the justification for the operational height of stratospheric airship to be around 20 km. However, a band of 17 km to 25 km needs to be explored to identify the optimum design with respect to a specified altitude of operation, for any specific location.

Solar Irradiance Model

Solar Irradiance Model by Dai and Fang [START_REF] Dai | A simple model to predict solar radiation under clear sky conditions[END_REF] has been integrated into the methodology for power output calculation from solar cell. Model estimates the direct solar irradiance as:

𝐼 𝑑 = 𝐼 𝑛 𝑒 [-0.103𝑚 𝑎 0.571 -0.081(𝜔𝑚 𝑟 ) 0.213 -𝜏 0.91 𝑚 𝑟 0.87 ] (17

)
𝑚 𝑟 = [𝑠𝑖𝑛 𝛼 + 0.15(3.885 + 𝛼) -1.253 ] -1 (18)

𝑚 𝑎 = 𝑚 𝑟 𝑝 𝑎 1013 (19)
Where, α is the solar elevation angle and mr and ma are the relative and absolute air mass respectively.

𝜔 = 𝜔 𝑚 𝑒 (-0.44Δ𝐻) (20)

𝜏 = 𝜏 𝑚 𝑒 (-0.691Δ𝐻) (21)
Where, τm and ωm are the values at measured site, and ∆H is the difference of altitude from measured site in km.

Objective Function and Design Variables

Several objective functions can be selected for determining the optimum envelope shape, e.g., minimum volumetric drag coefficient, minimum hoop stress and a composite objective function incorporating more than one objective function.

To consider the influences of various factors on the optimization of the shape, a composite objective function involving CDV, Ae, σmin and Asa is devised as follows: Design variables related to sizing are shown in the Fig. 5.

𝐹 𝑐𝑜𝑚 = 𝑤
For baseline case all the shape parameters are kept constant. For MDO all nine parameters have been given a bound as a range. Tab. II enlists upper and lower bound of all the design variables.

Results of Sensitivity Analyses

To understand the effects of various parameters sensitivity analyses have been carried out. Energy produced by the solar array of the airship are used by payload and system as well as it also fights with the wind to resist the drift and maintain quasi stationary position. However, performance of solar array is greatly affected by the temperature rise during day time. Upto 30% loss in efficiency of solar array can be observed for a typical operating condition and deployment site [START_REF] Alam | Modeling Transient Heat Transfer in Stratospheric Airships[END_REF]. Fig. 6 highlights the effect of efficiency of solar array on the sizing of airship. An attempt to highlight the seasonal variation on the sizing of airship due to the variation in availability of solar irradiance over the year. For same requirements of payload and power, significant variation in sizing can be observed in Fig. 7. Smallest airship fulfilling the needs of power and payload can be design at summer time due to ample availability of solar energy. While Fig. 7 keeps the operational height constant, interesting results can be found by varying altitude (as shown in Fig. 8). At different season, different operational heights act as optimum due to the fact of variation in wind speed.

Results of Optimization

In this section, the results after performing optimization using Genetic Algorithm (GA) are presented. Firstly, a baseline case was run keeping all the design variables related to shape of airship as constant to obtain the reference parameters. Then, optimization was run to obtain the optimal results. The comparison of optimized shape and reference shape for Fcomp listed in Eq. ( 17), is shown in Fig 10. It has been observed that, GA shows a good convergence to find the optimal solution for given constraints. It is also observed that the numerical values of shape coefficients (a1, a2, a3,…,a6) are very sensitive. For the given constraints, many combinations of design values to imaginary solution. To take care of the same, an additional penalty functions were imposed, to overlook the values which results in imaginary output of the shape function. The methodologies explained have some limitations especially in solar area calculation and structure modelling. There is a need to develop an all-encompassing methodology consisting of different models, which would be able to carry out sizing and arrive at the more realistic optimal configuration fulfilling any given requirements with low cost, high payload capability, and least size and weight. 8 Contact Author Email Address mailto:irfan@iitb.ac.in, irfan@aero.iitb.ac.in

Copyright Statement

The author confirms that he, and/or his company or organization, hold copyright on all of the original material included in this paper. The author also confirms that he has obtained permission, from the copyright holder of any third party material included in this paper, to publish it as part of their paper. The author confirms that he gives permission, or has obtained permission from the copyright holder of this paper, for the publication and distribution of this paper as part of the ICAS proceedings or as individual off-prints from the proceedings.

Fig. 3 .

 3 Fig. 3. Variation of wind with altitude

Fig. 4

 4 Fig. 4 plots the variations in available solar radiation for baseline cases i.e., for two different cities viz. Mumbai and Beijing, in Winter (15-Jan) and Summer (15-Jul) solstice.

Fig. 4 .

 4 Fig. 4. Solar Irradiance availability

Fig. 5 .

 5 Fig. 5. Airship profile with design variables

Fig. 6 .Fig. 7 .

 67 Fig. 6. Effect of Solar efficiency on sizing

Fig. 8 .Fig. 9 .

 89 Fig. 8. Optimum deployment altitude at different seasons

TABLE II :

 II Design variables with bounds

	Design Variable	Lower bound	Upper bound
	m	0.4	0.6
	r0	0	1
	.r1	0	1
	cp	0.6	0.7
	L/D	2	6
	XSC (m)	20	100
	ξ (deg)	50	100
	ALT (km)	17	25

TABLE I :

 I Input Parameters

	Input Parameter	Value
	Design Altitude	20,000 m
	Maximum design speed, v	25 m/s
	Payload mass	1000 kg
	Power required by Payload	10 kW
	Power required by Control Systems	11 kW
	Ratio of CDV, total / CDV, env	1.8
	Propulsive efficiency	90 %
	Solar array efficiency	12 %

Table

  III compares the results output for key parameters and improvements after optimization. It can have noticed that in almost each key parameters, optimized shape has better results.

TABLE III :

 III Output Parameters

	Ref. Shape	Opt. Shape	% Imprt.
	Solar Area	7889	7306	7.39
	Hull Area	31223	32514	-4.13
	Drag Coeff.	0.0365	0.0356	2.47
	Hoop Stress	1733.6	1649.4	4.86
	Comp. Func	100.0	97.3976	2.60
	Fig. 10. Comparing optimum shape with reference shape
	5 Conclusions
	A new and simpler methodology for
	multidisciplinary design optimization of
	stratospheric airship has been proposed, with
	lesser number of assumptions compared to other
	such methodologies reported in literature. The
	efficacy of the methodology in arriving at the
	output parameters, and in carrying out sensitivity
	analyses of key parameters has been
	demonstrated. Methodology couples the various
	model viz. solar, wind, aerodynamic, structure
	etc. for multidisciplinary design and analyses.