An Improved Differential Quadrature Time Element Method

A Differential Quadrature Time Element Method (DQTEM) was proposed by the author and co-worker, its drawback is the need of larger storage capacity since the dimension of the coefficients matrix for solution is the product of both spatial degrees of freedom and temporal degrees of freedom. To solve this problem, an improved DQTEM is developed in this work, in which the differential quadrature method is used to discretize both spatial and time domains, sequentially, and the dimension of the coefficients matrix is greatly reduced without losing solution accuracy. Theoretical studies demonstrate the improved DQTEM features superiorities including higher-order accuracy, adequate stability and symplectic characteristics. The improvement of DQTEM is validated by extensive comparisons of the present DQTEM with the original DQTEM.


Introduction
The Differential Quadrature Method (DQM) [1] is an effective numerical technique to be used as a global finite difference scheme or generalized collocation method.The study of DQM has been mainly focused on how to select grid points, how to impose initial or boundary conditions, how to determine weighting coefficients, and convergence of solutions.The relevant and representative publications are reviewed below in turn since DQM is applied to the discretization of spatial and temporal domains in the present work.
The extremum points of Legendre polynomial in the domain [−1, 1] [2] and the so-called Gauss-Lobatto-Chebyshev points [3] were taken as unequally aligned node points, and the latter has been widely accepted by far.
The efficient ways to deal with boundary conditions are formulating weighting matrix [4,5] except δ method [3,6], and modifying the weighting coefficient points [5,7], adding boundary degrees of freedom [8], employing Differential Quadrature Element Method (DQEM) [9], etc.The initial conditions can be applied via some representative ways, including expressing initial conditions as Differential Quadrature (DQ) analog equations at sampling points [2, 10,11], modifying DQ rule [10,12], modifying trial functions [13,14], and modifying weighting coefficient matrix [15,16].Notably, the advantage of modified DQ rule method is the accuracy of state variables at the end of time interval, and its accuracy can be improved up to 2n order by using only n + 1 sampling points with unconditional stability and non-dissipation.However, the power of the displacement polynomials is the same with that of the velocity, which violates the differential relation between displacement and velocity [10,17].
Xing and Guo changed coefficients matrix and generalized load vector to simplify the imposition of initial conditions [18].
As for the weighting coefficients, Quan and Chang obtained the first and second order explicit weighting coefficients by differentiating Lagrange polynomials [19,20].Shu and Richards derived high-order weighting coefficients from recursion formulae [21].
In addition, Bert and Wang pointed out that DQM convergence speed is equal to that of Fourier series solution while faster than that resulting from Ritz method; moreover, DQM is much simpler than Fourier series, which is easily implemented by coding and applicable to static and dynamic structural analysis [22].
In the past three decades, DQM has developed to an accurate and efficient numerical method for plenty of applications.In solid mechanics, DQM is usually used for plate or shell [23,24] structure analysis; many examples refer to mechanical behavior of a plate or panel [25,26], vibration of composite laminated plate [27,28] considering lower-order shear deformation, nonlinear bending problems of an orthotropic plate [29], nonlinear vibrations of a plate and nonlinear dynamics of viscoelastic beams [27,30], etc.In fluid mechanics, typical applications include incompressible viscous fluid flow analysis together with hydrodynamics and heat transfer [3], circulation around a circular cylinder at Reynold's number problems for incompressible viscous fluid flow [3,31], boundary layer equations solutions [32,33] and coupled heat transfer problems [2, 34,35].
In 2009, Xing and Liu [36] proposed a Differential Quadrature Finite Element Method (DQFEM) which was motivated by the complexity of imposing boundary conditions in DQM of strong form and the unsymmetrical element matrices in DQEM.In their later study, the explicit forms were presented for the C 0 and C 1 one-dimensional to three-dimensional structural element matrices and load vectors [37].
The innovative idea of discretizing time coordinates through DQ rule were initially given by Bellman et al. [1,38] and then extended to spatial domain.In 2012, Xing and Guo [18] formulated the Differential Quadrature Time Element Method (DQTEM), which is denoted as DQTEM0 in this paper for brevity, in which the time coordinate is also discretized by DQ rule and the spatial coordinate is discretized by the standard Finite Element Method (FEM).The DQTEM0 is much more accurate and efficient than conventional low-order time integration methods.Since Kronecker product of matrices are adopted in DQTEM0, the coefficients matrix is consequently the product of spatial and time Degrees of Freedom (DOFs), which results in considerable increase of the dimension of stiffness coefficient matrix.In this context, the DQTEM0 is improved in the present work, in which both spatial and temporal domains are discretized via DQM; therefore, the dimension of the coefficients matrix is greatly reduced without sacrifice of accuracy.
Generally speaking, it is agreed that a competitive method should be equipped with higher-order accuracy and broad stable range as well as high computational efficiency.Bathe and Wilson [39] presented a perfect and underlying process for the accuracy and stability analysis of time integration methods.By referencing to predecessors' procedure, this paper makes comprehensive researches on the improved DQTEM which is different from the Differential Quadrature Time Finite Element Method (DQTFEM) [40] proposed recently by the authors.DQTFEM solved the weak form of ordinary differential equations of motion based on the general form of Hamiltonian Variational Principle, the spatial domain was discretized by the standard FEM, and the Gauss-Lobatto-Legendre integration points were used in the calculation of the time derivatives through DQ rule.However, the improved DQTEM in this paper can be regarded as a combination of DQFEM and DQTEM0, in which the ordinary differential equations of motion in strong form is solved directly, and the structural matrices are calculated by using DQFEM.

Basic Rules and Formulations
In this paper, the spatial coordinate is discretized by using DQM, which is different from that in DQTEM0.The formulae of DQ rule and Gauss-Lobatto quadrature rule are the same for both spatial domain [−1, 1] and temporal domain [−1, 1], then for a better understanding of the present work, the spatial coordinate is used in the following for the brief introduction of DQ and Gauss-Lobatto quadrature rules.

The DQ Rule
The basic idea of DQ rule is to analogously obtain the derivative of a function or variable by a weighted linear sum of the function or variable values at all discrete points in the domain [41,42], the r-th derivatives of f (x) at the i-th point can be expressed by where ij are the weighting coefficients and m is the number of sampling points.The weighting coefficients have explicit forms [19] as Higher-order weighting coefficients can be obtained by the sum of products of lower orders as [21]

Gauss-Lobatto Quadrature Rule
The Gauss-Labatto quadrature formula for function f (x) with precision degree (2m − 3) is where the quadrature points are the zeros of polynomial as the quadrature coefficients C * j are given by

The Formulae of the Improved DQTEM
The ordinary differential equation in structural dynamics has the form as where M×M are the mass, damping and stiffness matrices which are formulated explicitly by using spatial DQFEM, and will be indicated below.The displacement vector x and load vector f are defined as By using the DQ rule, the velocities at sampling points can be obtained .
where n is the number of DOFs of time element, a dot over the letter denotes the time derivative; the first subscript i of x or .
x stands for the spatial DOF number, the second subscripts j and k stand for the time point number; and .
x ik and x ij represent the values of .
x i at time point k and x i at time point j, respectively.Note that Gauss-Lobatto-Chebyshev points [3] are employed as time sampling points in present work.Let based on the DQ rule, the velocity vector .x i and acceleration vector ..
x i can be written in a compact form as . ..
where A t and B t are the first and second order weighting coefficient matrices with respect to time, respectively.Substituting Equations ( 10) and ( 11) into Equation ( 7) yields and where has the similar definition as x i .According to matrix theory, Equation (12)  can be transformed into another form where Z = csX), R = cs(F), "cs" indicates column expansion of a matrix, and where "⊗" is Kronecker product, and E is an n × n unit matrix.Introducing initial conditions into Equation ( 14), one can solve displacements X, and the velocities and accelerations can be achieved through Equations ( 10) and ( 11).
For a system with Rayleigh damping, its damping matrix can be represented by where α and β are the pre-determined constants [43,44].When C = 0, Equation ( 12) can be solved accurately and efficiently by the Bartels and Stewart method [45].

DQ Rod and Beam Elements
For the comprehensive understanding to obtain K and M matrices based on DQ rule and the minimum principle of total potential energy, the C 0 DQ finite rod element and the C 1 DQ finite beam element are given below, respectively.

C 0 DQ Rod Element
Assume the longitudinal displacement function of a uniform rod element is where l i stands for Lagrange polynomials, u i = u(x i ) is the displacement at the Guass-Lobatto quadrature point, and x i is the coordinate of the Guass-Lobatto quadrature point.According to DQ rule and Gauss-Lobatto quadrature rule, the total potential energy functional and kinetic energy coefficient functional are where E and ρ are the Young's modulus and mass density, respectively; S and L are the across section area and length of rod element, respectively; and f i is the distributed force.The mass matrix M, stiffness matrix K, generalized load vector F and other variables are given by where the matrix A s is the first order weighting coefficient matrix with respect to space, and C * s is the Gauss-Lobatto quadrature coefficient diagonal matrix.

C 1 DQ Euler Beam Element
Assume the deflection of a uniform Euler beam element is where w i = w(x i ) is the deflections at the Guass-Lobatto quadrature nodes, similar to u i .The total potential energy and kinetic energy coefficient functions are where B s = (A s ) 2 , I is the moment of inertia, and In order to satisfy C 1 inter-element continuity requirement of Euler beam element, the displacement vector is modified to Based on DQ rule, the relation between w and w is readily given by w = Tw (30) where m,2 A (1) Substitute Equation (30) into Equations ( 26) and ( 27), the mass matrix M, stiffness matrix K and generalized load vector F are obtained as In fact, the mass matrix, stiffness matrix and generalized load vector of other structures with arbitrary geometry, such as plane , Kirchhoff plate and shells, etc., can be obtained by the same way [37].Note that DQ rule cannot be used directly for irregular geometric domain; it should be transformed from the Cartesian x-y coordinate system to the natural ξ-η coordinate system, so the DQ rule should be reformulated, and Bert and Malik [46] have done this work first, and Xing and Liu [36] gave a very detailed description about the transformation.

Discussion
(1) Unlike previous high-order methods, the present method is easily implemented, and the shape functions are no long necessary for constructing spatial element matrices, which are explicit and can be obtained by simple algebraic operations of weighing coefficient matrices of the DQ and Guass-Lobatto quadrature rules.(2) For any C 0 and C 1 problems, Lagrange polynomials are used as trial function; in addition, the Guass-Lobatto quadrature points are selected as nodes and collocated not uniformly.(3) The mass matrix of DQ finite element is diagonal, which is different from conventional diagonal lumped mass matrix, but the element mass matrix used in DQTEM0 is full.(4) Compared with DQTEM0, the dimension of the coefficient matrix G of the present method is reduced and the present method is more efficient than DQTEM0, which are validated by numerical results in Section 4. (5) To avoid confusion, the comparison among DQTFEM, DQTEM0 and the improved DQTEM is presented in Table 1.In addition, these DQ time element methods seem similar to the standard spectral method in some sense, since they are all based on polynomial interpolation.The standard spectral method is global interpolation, in which the Chebyshev polynomials, the Legendre polynomials and the Fourier series or harmonic series are employed in general.While the DQ time element methods employ Lagrange polynomials, and if only one element is employed in time and/or spatial domain, they are global interpolation; if the domain is divided into several elements, the DQ time element methods are local interpolation.x i = A t x i , . .
DQ finite element with non-uniform nodes .

Accuracy and Stability
By virtue of modal decomposition, the linear undamped multi-DOFs system can be equivalently transformed to a set of single-DOF systems; therefore, in the analysis of accuracy and stability of a new method, single-DOF system is used generally.Consider the following single-DOF system ..
x + (ωh) where ω is the angular frequency, h is time step size, here an over dot represents the differential with respect to the non-dimensional time τ = t/h, t is the real time.In this case, displacement vector x is where x i (i = 1, 2, . . ., n) denote the displacements at time sampling points, and According to the method of imposing initial conditions in Reference [18], the relationship between the n-th sampling point and the initial displacement and velocity is x 0 (36) and the elements of the Jacobi matrix J are where S = G −1 , then stable interval can be obtained by using criterion ρ(J) ≤ 1, ρ(J) is the spectral radius, which is shown in Figure 1, and the first stable intervals are listed in Table 2.It should be stressed that, for different sampling point n, ρ(J) > 1 or ρ(J) = 1.Before discussing the accuracy of the improved DQTEM, the definition of sympletic method is recalled first.If the Jacobi matrix J of a method satisfies the following conditions , , , J = J J J J I J J J J J I J J J J J 0 J J J J J 0 J J J J J 0 J J J J J 0 (38) then the method is symplectic which preserves energy and momentum of conservative system.The present improved DQTEM is proved to be highly symplectic (see Table 3), from which one can see Equation ( 37) is satisfied excellently, especially J3-J6 are rigorously satisfied.This implies the improved DQTEM has high-order accuracy for amplitude and phase, which will be validated below by numerical comparisons.Note: max (ωh) indicates the upper limit of ωh for the first stable interval in Table 2.

Numerical Analysis
This section aims at demonstrating the superiority of the improved DQTEM over the DQTEM0, the free vibrations of rod and beam are first solved, then the impact problem of rod.The used material properties are listed in Table 4.All numerical results are obtained by MATLAB.Before discussing the accuracy of the improved DQTEM, the definition of sympletic method is recalled first.If the Jacobi matrix J of a method satisfies the following conditions then the method is symplectic which preserves energy and momentum of conservative system.The present improved DQTEM is proved to be highly symplectic (see Table 3), from which one can see Equation ( 37) is satisfied excellently, especially J 3 -J 6 are rigorously satisfied.This implies the improved DQTEM has high-order accuracy for amplitude and phase, which will be validated below by numerical comparisons.

Numerical Analysis
This section aims at demonstrating the superiority of the improved DQTEM over the DQTEM0, the free vibrations of rod and beam are first solved, then the impact problem of rod.The used material properties are listed in Table 4.All numerical results are obtained by MATLAB.

Rod Euler Beam
Young

Frequency Comparisons
The natural frequencies of a rod and a Euler beam are obtained by the improved DQTEM (or using the spatial DQ finite element) and DQTEM0 (or using linear rod element and cubic beam element); the results are compared with exact ones to show the advantage of the proposed method.
For evaluating the computational accuracy of natural frequency, a relative error is defined where the exact frequencies are obtained by solving characteristic differential equations according to boundary conditions.For evaluating the effectiveness of both methods in calculating frequencies, a percentage is defined as where Merr5 is the number of frequencies with RelErr1 less than 5% and M is the total number of frequency.

Uniform Rod
This subsection aims to show the improved DQTEM is more accurate than DQTEM0, if the spatial DOFs used by both method are the same.
5 gives the comparisons of natural frequencies via constant strain rod element and DQ element for a fixed-free rod, in which the dimensionless frequency parameter λ = ωl/c is used, where c = E/ρ is the longitudinal wave propagation velocity.In Table 5, NE is the number of spatial element, and m is the DOF of each spatial element.
For this problem, DQTEM0 adopts 10 and 100 constant strain elements, respectively.For comparison using the same number of spatial DOFs, the improved DQTEM employs a DQ element with 11 and 101 unequally spaced points (Gauss-Lobatto-Chebyshev points), correspondingly.It can be concluded from Table 5 that the Nerr5 of DQTEM0 with NE = 10, m = 2 is about 30%, while the Nerr5 of the corresponding improved DQTEM with NE = 1, m = 11 is almost 70%, and the RelErr1 of the improved DQTEM are much smaller than DQTEM0.A similar conclusion is also apparent when NE = 100, m = 2 is used in DQTEM0 and NE = 1, m = 101 is used in the improved DQTEM.If more spatial DOFs are used, this conclusion can also be observable.
To further show the superiority, the improved DQTEM is also compared with High-Order Lagrange Element Method (HOLEM) with equally spaced nodes in Table 6.It can be seen in Table 6 that the Nerr5s of the improved DQTEM and HOLEM are the same since they use the same order of interpolation polynomials.However, the improved DQTEM is more accurate and efficient because, in the improved DQTEM, the unequally spaced nodes are used, and the mass and stiffness matrices are only the products of the explicit weighting coefficient matrix and the coefficient matrix of Gauss-Lobatto integration, while those of HOLEM depend on shape functions.
From the comparisons in Tables 5 and 6, one can see that the improved DQTEM is superior to DQTEM0 if the same number of spatial DOFs is used for the solution of natural frequencies, which is further validated by the results of the following beam.A simply supported Euler beam is considered here, the cubic element involving deflection and slope at each node is used in DQTEM0, but all DOFs of an element in the improved DQTEM are deflections.Non-dimensional frequency parameter λ = βl is introduced in comparison, where Similar to the rod, DQTEM0 adopts 10 and 50 cubic elements, respectively, whereas the improved DQTEM employs one element with 22 and 102 Gauss-Lobatto-Chebyshev points; apparently, they have the same number of spatial DOFs.It follows in Table 7 that the Nerr5 of DQTEM0 with NE = 10, m = 4 is nearly 45%, while the Nerr5 of the improved DQTEM with NE = 1, m = 22 is about 70%, and the RelErr1 of the improved DQTEM are much smaller than DQTEM0.A similar conclusion is achieved when DQTEM0 uses NE = 50, m = 4 and the improved DQTEM uses NE = 1, m = 102.
Note that, in the cubic beam element, half of the spatial DOFs is deflection, and half is slope, but all DOFs in the improved DQTEM are deflections.This motivates us to compare the case with the same degrees of deflections (see Table 8), in which the dimension of DQTEM0 model is nearly twice of that of the improved DQTEM.Although the results of DQTEM0 in Table 8 are more accurate than those in Table 7, it is also obvious that the improved DQTEM is much more accurate than DQTEM0.In other words, the improved DQTEM is more accurate than DQTEM0 when the same spatial DOFs are used.That is, the dimensions of G can be reduced by using the improved DQTEM without sacrifice of accuracy.

Comparison of Dynamic Responses
For evaluating the performance of the improved DQTEM in dynamic analysis, this section investigates the transient problem of a fixed-free uniform rod subjected to an impact at free end by a rigid body of mass 1.5ρSl.Initial conditions are zero except for the initial velocity of the free end colliding with the body is v.The relative error used in this section is defined by RelErr2 where "Exact" denotes the exact displacements from the analytical mode superposition method at each time sampling point; the infinity norm of "Exact" denotes the maximum absolute value of a vector, and the elements of this vector are the exact values at all time nodes; and "Approximate" denotes the displacements of DQTEM0 or the improved DQTEM at each time sampling point.
Consider the dimensionless time domain [0, 4000].Table 9 compares the efficiency and the used spatial DOFs for achieving the same accuracy, in which Δ is the dimensionless time step size.For this transient problem, the improved DQTEM requires much fewer spatial DOFs than DQTEM0; for different RelErr2, the spatial DOFs used in the improved DQTEM are two-fifths as many as that of the DQTEM0, i.e., the coefficients matrix G of the improved DQTEM is only 4/25 of the DQTEM0, which leads to considerable saving of computational space and CPU calculation time.All these superiorities of the improved DQTEM are due to the more finely description of spatial domain.By using M = 16, n = 15, we give the comparisons of the dimensionless displacement and velocity at free end or impact point with exact solutions, as shown in Figure 2. The dimensionless time τ for longitudinal wave travelling from impact point to fixed end and then back to impact point is 2, which is shown clearly in the time history of velocity as the wave characteristics of linear elastic impact of rod.

Conclusions
In this paper, the improved Differential Quadrature Time Element Method (DQTEM) was proposed, in which the DQ rule is employed to discretize both spatial and time domains.
It has been shown that, to achieve the same accuracy, the spatial DOFs in the improved DQTEM are much fewer than the original DQTEM due to the use of high-order spatial DQ finite element in the improved DQTEM; the improved DQTEM is a highly accurate symplectic method, which is suitable for long-term simulation, and the improved DQTEM is accurate and efficient for the analysis of both free vibrations and forced vibrations.
These advantages of the improved DQTEM were validated by comprehensive numerical comparisons.

Conclusions
In this paper, the improved Differential Quadrature Time Element Method (DQTEM) was proposed, in which the DQ rule is employed to discretize both spatial and time domains.
It has been shown that, to achieve the same accuracy, the spatial DOFs in the improved DQTEM are much fewer than the original DQTEM due to the use of high-order spatial DQ finite element in the improved DQTEM; the improved DQTEM is a highly accurate symplectic method, which is suitable for long-term simulation, and the improved DQTEM is accurate and efficient for the analysis of both free vibrations and forced vibrations.
These advantages of the improved DQTEM were validated by comprehensive numerical comparisons.

Figure 1 .
Figure 1.Spectral radii for different n of the improved DQTEM (Differential Quadrature Time Finite Element Method).

Figure 1 .
Figure 1.Spectral radii for different n of the improved DQTEM (Differential Quadrature Time Finite Element Method).

Table 5 .
Frequency parameter λ of rod by linear rod element and DQ (Differential Quadrature) element.
.x i = A t x i ,

Table 2 .
The first stable intervals with ρ(J) = 1 for different n.
Note: max (ωh) indicates the upper limit of ωh for the first stable interval in Table2.

Table 6 .
Frequency comparisons of rod using high-order rod elements.

Table 7 .
Frequency comparisons of beam with the same spatial DOFs (Degrees of Freedom) for both methods.

Table 8 .
Frequency comparisons of beam with same deflection degrees for both methods.

Table 9 .
Efficiency comparison of the improved DQTEM and DQTEM0.