1 Introduction

Numerical methods play an important role in solving time dependent problems. Generally, dynamic problems can be categorized into structural dynamic and wave propagation problems.The majority of these methods produce spurious oscillations when they are applied for the numerical modeling of transient wave propagation problems^[1]. The standard results of the numerical simulations of transient wave propagation problems will be affected by the dispersion or dissipation errors caused by both spatial and temporal discretizations in a combined way, and thus, the poorly constructed discretizations will often induce a significant loss in the accuracy of the wave propagation solution^[2].

A large amount of research has been focused on the dispersion analysis of finite element solutions to the wave or Helmholtz equation^[3-7]. The standard finite element method may not be effective since very fine meshes are required. Moreover, even when such fine meshes are adopted, the solution still shows obvious spurious oscillations. Although enriched finite element methods with appropriate time integration schemes have been presented to overcome the above shortcomings^[8-9], numerical results demonstrate that there is still significant need for improved solution procedures. Some other space-discretization methods, e.g., spectral element method^[10-14], meshless method^[15-21], and isogeometric element method^{[8, 22-24]}, are also proposed for the improved solutions of wave propagation problems. It is noteworthy that isogeometric elements yield more accurate numerical results for wave propagation problems than high-order finite elements^[25-27]. Dispersion and frequency errors for isogeometric analysis (IGA) have been reported to decrease when the order of splines increases^[28-29]. For conventional isogeometric elements, the consistent mass matrix obtains a higher-order of dispersion errors than the lumped mass matrix^[28-30]. Despite the improved accuracy, the aforementioned spatial discretized elements may yield spurious high-frequency oscillations when high-frequency and impact loadings are considered, and these spurious oscillations often lead to inaccurate and even divergent results at mesh refinement^[31-32]. The improved isogeometric elements^[33-34] for wave propagation problems have thus been developed to obtain higher-order dispersion errors for both lumped and consistent mass matrices. Here, the order of the dispersion errors for the consistent mass matrix can be improved from order 2p (the conventional elements) to order 4p (the improved elements), where p is the order of the polynomial approximations, and the dispersion error for the lumped mass matrix can be improved from the second-order to the 2pth-order.

Actually, in the solution of wave propagation problems, the errors from the spatial and temporal discretizations often appear together and affect each other^{[3, 35-37]}, and an effective method for reducing the dispersion error of wave propagation problems is always used to select the appropriate direct time integration schemes, e.g., the explicit and implicit methods^[38-40] for different space-discretization elements. However, it is not easy to obtain the desirable schemes due to the difficulty of damping out the spurious high-frequency modes with no substantial accuracy loss in important low-frequency modes, i.e., the time integration scheme with good numerical dissipation characteristics is preferred. For example, when the lumped mass matrix is used for the standard finite element method (FEM), the non-dissipative central difference scheme is not recommended since the solution accuracy can be severely ruined by the dispersion errors in the high frequency. By comparison, the dissipative Noh-Bathe scheme^[41] shows good performance in the analysis of wave propagations due to its desirable numerical dissipation characteristics and second-order accuracy.

As for the isogeometric element method, there is few research in the literature devoted to the introduction of optimal numerical dissipation in explicit and implicit time integration schemes. Therefore, in this paper, in order to extend the isogeometric element method in explicit dynamics, the dispersion error analysis of the improved dissipative quadratic isogeometric elements and the central difference time integration scheme is first conducted to obtain the optimal Courant-Friedrichs-Lewy (CFL) numbers (i.e., optimal time sizes) for the solution of the wave propagation problems. The paper is organized as follows. First, the IGA solution procedure of wave propagation problems is briefly described in Section 2, where both the space-discretization elements (i.e., the improved quadratic isogeometric elements) and the time-discretization scheme (i.e., the central difference scheme) are presented. In Section 3, the stability property of the central difference scheme in the proposed IGA solution procedures is studied to obtain the critical time-step size. Further, the dispersion properties of the proposed IGA solution procedures are investigated analytically and numerically in Section 4, where the quantified effect of different time-step sizes on the dispersion error is obtained. Subsequently, in Section 5, the one-dimensional (1D) benchmark wave propagation problems are used to show the validity of the proposed IGA solution procedures. Finally, some concluding remarks of this paper are given in Section 6.

2 IGA solution procedure for wave propagation problems 2.1 Proposed quadratic isogeometric elements

In this paper, the considered wave propagation in an isotropic homogeneous medium is described by the following scalar wave equation in the domain Ω as follows:

(1)

where t is the time, ▽² is the Laplace operator, u is the solution variable, and c₀ is the wave propagation velocity. Here, body forces are not considered.

The application of the continuous Galerkin approach and the space discretization, i.e., the isogeometric element method, to Eq. (1) leads to a system of ordinary differential equations in time as follows:

(2)

where

(3)

d is the vector containing all control points related to the IGA, and is the second derivative of d with respect to time t. M and K are the global mass matrix and the stiffness matrix, respectively, which are obtained by the usual summation of the isogeometric element (Ω^e) matrices (see Eq. (3)).

For conventional isogeometric elements, K^e and M^e are defined by^[33]

(4)

(5)

where N and are the vector of shape functions and its derivative vector in terms of x, i.e., a position vector from the origin of the coordinate system. Then, M in Eq. (3) is the consistent mass matrix, and K in Eq. (3) is the sparse banded matrix.

In this study, to obtain the second-order dispersion error for the lumped mass matrix M, we adopt K^e and M^e for 1D quadratic isogeometric elements as follows^[33]:

(6)

where Δx is the characteristic length of the quadratic isogeometric element. K^e can be obtained from Eq. (4). a₁ is an arbitrary coefficient that does not affect the dispersion error. To ascertain a₁, we introduce the lumped matrix definition as follows:

(7)

From Eqs. (6) and (7), we have a₁ = 80.

The fourth-order dispersion error can be achieved by only modifying the element stiffness matrix K^e as follows^[33]:

(8)

For brevity, in the following analysis, we denote the aforementioned two proposed elements of the second-order (see Eq. (6)) and fourth-order (see Eq. (8)) dispersion errors, respectively, as IGA-s and IGA-f.

2.2 Central difference scheme for the IGA solution procedure

In terms of the high-order dispersion errors of the adopted isogeometric elements, it suffices to use the second-order central difference scheme based on^[39]

(9)

where Δt is the time-step size. Substituting the latter relationships of Eq. (9) into Eq. (2) at t_n leads to a system of algebraic equations for d_n+1 at t_n + Δt, i.e., t_n+1, as follows:

(10)

where the effective quantities are

(11)

Using Eqs. (2), (9), and (11), we obtain a linear multistep form of the central difference scheme in the modal basis

(12)

or for all equations

(13)

where the superscripts t + Δt, t, and t − Δt denote the discretized time, x is a modal degree of freedom, X is the vector of all modal degrees of freedom, ω is the natural frequency of a generic mode of the isogeometric element model, and Λ is the corresponding diagonal matrix listing all ω_i².

Moreover, with the eigenvector Φ of the problem

(14)

and the definition

we can rewrite Eq. (13) as the isogeometric element degrees of freedom as follows:

(15)

2.3 Stability analysis for the proposed IGA procedure

In this section, the von-Neumann stability criteria are employed for the stability analysis of the proposed IGA procedures, where two isogeometric elements, i.e., the IGA-s and IGA-f quadratic elements, are considered.

Assume a general solution form for Eq. (2) as follows:

(16)

where is any term of the field variable vector d corresponding to the time t = t_n and the spatial position x = x_j = jΔx. The wave number k = ω/c, where c is the numerical wave velocity. A_k is the complex coefficient corresponding to the spatial wave number k. ξ represents the time dependence of the solution, and is called the amplification factor.

Considering the discretized equation (15) on an infinite line with the sequence of Δx-spaced control points x_j = jΔx, with the help of Eqs. (3) and (6), we can decompose Eq. (15) for the IGA-s quadratic element as follows:

(17)

Substituting Eq. (16) into Eq. (17) and using e^ikjΔx = cos(kjΔx) + i sin(kjΔx) yield

(18)

For stability, we require |ξ|≤1. Then, the stability condition can be easily obtained as follows:

(19)

As for the IGA-f quadratic element, with the same procedure as Eqs. (17)–(19), we have

(20)

(21)

and the stability condition

(22)

Actually, with Eq. (13), we can express the stability conditions for the IGA-s and IGA-f quadratic elements in a unified form as follows:

(23)

Then, the adopted Δt for the IGA-s and IGA-f quadratic elements are

(24)

3 1D demonstrative dispersion analysis

In this section, we will investigate the dispersion errors of the solutions to the stencil wave propagation equation (1) with the consideration of the IGA-s and IGA-f quadratic elements and the central difference scheme. Here, we take the general solution to Eq. (2) as follows^[42]:

(25)

where is the complex coefficient corresponding to the spatial wave number k.

For the IGA-s element, substituting Eq. (25) into Eq. (17) yields

(26)

where is the relative wave velocity error, which also denotes the dispersion error.

For the IGA-f element, substituting Eq. (25) into Eq. (20) yields

(27)

In Eqs. (26) and (27), directly solving the relative wave velocity error is not feasible. Here, we prescribe that {(kΔx)₀, (kΔx)₁, …, (kΔx)_I, (kΔx)_I+1, … } is an increasing sequence of real numbers, where (kΔx)₀ = 0, which discretizes the continuous variable kΔx. Then, conducting the Taylor series expansion about Eqs. (26) and (27), we have

(28)

where is the calculated relative wave speed error corresponding to (kΔx)_I, l ∈ ℤ, l ≥ 2, and = 0.

For given (kΔx)_I+1 and N_CFL, we obtain a convergent iteration equation based on Eq. (28) as follows:

(29)

where I ∈ ℕ, J ∈ ℕ, and

In the above equation, J is the iteration number determined by

(30)

in Eq. (29) is given as follows:

(31)

where = 10⁻¹⁰, and l = 5.

With Eqs. (29)–(31), we can obtain , i.e., or , corresponding to (kΔx)_I+1 (see Figs. 1 and 2, where the relative errors of the IGA-s and IGA-f quadratic elements with the central difference scheme for various N_CFL are plotted in terms of the wavelength and the spatial element size , respectively). In Fig. 1, for the IGA-s quadratic element, the relative error augments when the CFL number N_CFL decreases. In Fig. 2, for the IGA-f quadratic element, the curves of the relative errors move downward when N_CFL decreases. It should be mentioned that the maximum N_CFL is 1.63 for the IGA-s element and is 0.86 for the IGA-f element, which are, respectively, obtained according to Eqs. (19) and (22). In Fig. 1, the case N_CFL = 1.63 represents the optimal dispersion error of the IGA-s quadratic element. In Fig. 2, the case N_CFL = 0.4 represents the minimal dispersion error of the IGA-f quadratic element shows minimal dispersion error.

4 Numerical simulation 4.1 1D elastic bar under suddenly applied load

In this section, a stencil wave propagation problem (see Fig. 3) is considered to test the validity of the proposed solution procedures. For this problem, all low and high frequencies are activated. In the figure, the elastic bar of length L = 4 under a suddenly applied load is simulated, and the wave velocity is chosen to be c₀ = 1 m·s⁻¹. The initial and boundary conditions are all given in Fig. 3.

In Fig. 4, the numerical velocities from the IGA-s quadratic element and the central difference scheme are illustrated to demonstrate the validity of the proposed solution procedures. In the figure, various values of N_CFL are considered. The CFL number N_CFL adopted here is defined as N_CFL = (see Subsection 2.2). It can be noticed from Fig. 4 that the case N_CFL = 1.6 shows the best results among all considered CFL cases for different spatial element numbers N_s = 200 and N_s = 400, which is in good agreement with the analytical dispersion error results shown in Fig. 1. In Fig. 4, the numerical results at larger observation time t = 10.000 s are less accurate than those at t = 2.000 s due to more activated high-frequency spurious modes. By comparison, in Fig. 5, the numerical results from the IGA-f quadratic element show higher calculation accuracy than those from the IGA-s quadratic element (see Fig. 4). According to the curves in Figs. 4 and 5 and the dispersion error analysis in Section 3, N_CFL = 1.6 and N_CFL = 0.4 are, respectively, suggested as the optimal CFL cases for the IGA-s and IGA-f quadratic elements.

4.2 Propagation of sinusoidal pulse in 1D elastic bar

In this section, the elastic bar model of L = 4 (see Fig. 3) is considered with two free ends. For brevity, the wave velocity is chosen to be c₀ = 1 m·s⁻¹. The initial displacement is given by u(x, 0) = 0, and the initial velocity is

(32)

where the frequency Ω = 5π.

In Figs. 6(a) and 6(b), when both N_s = 201 and N_s = 401 are, respectively, considered for the simulation, the numerical velocity curves illustrate that the solution procedure of the IGA-s quadratic element and central difference scheme gives clear numerical oscillations. For all considered CFL cases as shown in Fig. 6, the numerical results from the suggested case N_CFL = 1.6 are of the highest accuracy or smallest dispersion error, which roughly matches well with the dispersion error tendency as displayed in Fig. 1. As for the IGA-f quadratic element, Fig. 7 shows that the solution procedure of the IGA-f quadratic element gives far more accurate numerical results than the IGA-s quadratic element (see Fig. 6). Particularly, when the adopted spatial isogeometric element number increases from N_s = 201 to N_s = 401, obvious improvement in the calculation accuracy can be observed from Figs. 7(a) and 7(b). Different from the stencil wave propagation problem in Subsection 4.1, the time-step size Δt here is defined by

with the used central difference time integration scheme.

In Fig. 7, the case N_CFL = 0.4 of the IGA-f element demonstrates the best calculation accuracy among all considered CFL cases, which is consistent with the dispersion error results illustrated in Fig. 2.

5 Concluding remarks

In this paper, to develop the IGA application in explicit dynamics, two improved isogeometric quadratic elements coupled with the central difference scheme are used to form the solution procedures for transient wave propagation problems. Different from the extant solution procedures of wave propagation analysis, new solution procedures are first obtained by combining high-order isogeometric elements and an explicit time integration scheme, where the lumped matrices of the improved isogeometric quadratic elements are used. The theoretical dispersion error analysis of the proposed solution procedures is conducted to investigate the numerical dissipation/dispersion characteristics of the proposed solution procedures. 1D wave propagation simulations demonstrate the validity of the proposed solution procedures. In general, the major conclusions can be summarized as follows:

(ⅰ) Two improved isogeometric quadratic elements and the central difference scheme are first used to formulate the explicit solution procedures of transient wave propagation problems, and the stability characteristics of two isogeometric elements for the central difference scheme are first studied and obtained.

(ⅱ) The dispersion error analysis of the proposed solution procedure is conducted to ascertain the optimal CFL number N_CFL or the time-step size Δt for the presented second-order and fourth-order elements.

(ⅲ) Numerical simulations and dispersion error analysis confirm the effectiveness of the proposed solution procedures, and demonstrate that the improved isogeometric quadratic element of the fourth-order dispersion error (IGA-f) possesses desirable numerical dispersion characteristics.