0
Review Article

# Weak Form Quadrature Element Method and Its Applications in Science and Engineering: A State-of-the-Art ReviewOPEN ACCESS

[+] Author and Article Information
Xinwei Wang

State Key Laboratory of Mechanics and
Control of Mechanical Structures,
Nanjing University of
Aeronautics and Astronautics,
Nanjing 210016, China

Zhangxian Yuan

School of Aerospace Engineering,
Georgia Institute of Technology,
Atlanta, GA 30332-0150

Chunhua Jin

State Key Laboratory of Mechanics and
Control of Mechanical Structures,
Nanjing University of
Aeronautics and Astronautics,
Nanjing 210016, China;
School of Architecture Engineering,
Nantong University,
Nantong 224019, China

1Corresponding author.

Manuscript received January 5, 2017; final manuscript received April 24, 2017; published online May 16, 2017. Assoc. Editor: Martin Schanz.

Appl. Mech. Rev 69(3), 030801 (May 16, 2017) (19 pages) Paper No: AMR-17-1002; doi: 10.1115/1.4036634 History: Received January 05, 2017; Revised April 24, 2017

## Abstract

The weak form quadrature element method (QEM) combines the generality of the finite element method (FEM) with the accuracy of spectral techniques and thus has been projected by its proponents as a potential alternative to the conventional finite element method. The progression on the QEM and its applications is clear from past research, but this has been scattered over many papers. This paper presents a state-of-the-art review of the QEM employed to analyze a variety of problems in science and engineering, which should be of general interest to the community of the computational mechanics. The difference between the weak form quadrature element method (WQEM) and the time domain spectral element method (SEM) is clarified. The review is carried out with an emphasis to present static, buckling, free vibration, and dynamic analysis of structural members and structures by the QEM. A subroutine to compute abscissas and weights in Gauss–Lobatto–Legendre (GLL) quadrature is provided in the Appendix.

<>

## Introduction

“Along with the ever-growing advancement of faster computing machines, the research into the development of new methods for numerical solution of problems in engineering and physical sciences also is an ongoing parallel activity. Such research interests, of course, remain motivated by needs of modern technology” [1]. The QEM, proposed by Striz et al. in 1994 [2], is gradually emerging as a distinct numerical solution technique for the initial-and/or boundary-value problems. In principle, the QEM is established based on the differential quadrature (DQ) law, which has the capability of obtaining highly accurate solutions with a minimal computational effort. It has successfully solved numerous problems in structural mechanics, which are the static, buckling, free vibration, and dynamic response of structural members as well as structures.

Currently, there are two kinds of quadrature element method. One, called the strong form quadrature element method [2,3], the differential quadrature element method (DQEM) [4,5], or the strong formulation finite element method [6], is formulated based on the governing as well as boundary equations [26], and the other, called the weak form quadrature element method, is formulated based on the principle of minimum potential energy [5,715]. The weak form quadrature element method is more flexible than the strong form quadrature element method in dealing with complex geometry since it is essentially a higher order FEM. The method combines the generality of the finite element method with the accuracy of spectral techniques. This paper focuses on the weak form quadrature element method. The method is called the QEM for simplicity.

In the early version of the QEM [79], any kind of grid points can be used as the element nodes. During formulations, however, it requires computing the inverse of a full Vandermonde matrix, and thus the number of elemental nodes should be fixed and cannot be greater than 15. In 2004, Franciosi and Tomasiello proposed a modified quadrature element method. Nonuniformly distributed points are used as the element nodes [15]. The method also belongs to the weak form QEM since the stiffness matrix is derived based on the potential energy.

To remove the deficiency existing in the early version of the QEM and raise the computational efficiency, GLL points are used as the element nodes and GLL quadrature is used to obtain the element stiffness matrix and mass matrix [1012]; thus the number of element nodes is not limited and the derivatives at nodal points can be calculated explicitly by using the available formulations in differential quadrature method (DQM). This reduces the difficulty in calculations of the derivatives for high-order elements, the round-off error, and the programming effort. To use the DQ law, however, the element nodal points should be the same as the quadrature points; thus, only GLL points can be used as the element nodes. In literatures, the Lobatto-type QEM is often regarded as the time-domain (SEM) or spectral finite element method (SFEM) [16]. The SEM is popular for simulating the wave propagation and an established technique in the wave and fluid fields. However, differences between the QEM and SEM do exist [79], for example, both Lagrange and Hermite-type shape functions are used in the QEM [5,714].

To remove the limitation on the selection of element nodes, a general way is proposed [13,14]. Any type of grid points, such as the expanded Chebyshev (EC) points [17,18] and the very good approximations to the Lebesgue-optimal grid points [19], can be used as element nodes and Gauss quadrature can be used during formulations of the stiffness and mass matrices to raise the accuracy and computational efficiency. The expanded Chebyshev point allows the longest time step [17,18]; thus, it is computationally more efficient in dynamic analysis by using explicit time integrations. Besides, a short program is available [5,20] for computing the abscissas and corresponding weights in Gauss quadrature. The research work presented in Refs. [13,14] makes the formulation of a weak form quadrature element general enough. Since any basis function can be considered by the developed QEM now, the formulation of the QEM is more general than the existing SEM or SFEM.

Since the SEM is a special case of the QEM [16], the quadrature element method also possesses the high accuracy of spectral methods with the flexibility of finite element method. Research results show that the QEM [5,714] exhibits the capability of producing highly accurate solutions with minimal computational effort. Therefore, the method is seemingly a potential alternative to the conventional numerical methods such as the finite difference method (FDM) and the FEM.

The progression on the weak form quadrature element method and its applications is clear from past research, but this has been scattered over many papers. This paper presents a state-of-the-art review of this subject. In what follows, the basic mathematical concepts underlying the QEM are presented first. Three types of elements, i.e., the Lagrange type, the Hermite type, and the mixed Lagrange and Hermite type, are formulated. Since the QEM is possibly not well known to the computational mechanics community due to its rather recent origin, this section is written in detail. In Sec. 3, the time domain spectral element method is briefly reviewed, and remarks on the similarities and differences between the QEM and SEM are given. After that, numerical examples on the application of the QEM are given. Then, a review of the chronological development of the weak form quadrature element method is presented. The paper is ended with conclusions.

## Weak Form Quadrature Element Method

In this section, the basic mathematical concepts underlying the weak form QEM are presented. The formulas of bar, Euler–Bernoulli beam, and Kirchhoff plate elements are given for illustrations. The procedures can be extended to formulate other types of weak form quadrature elements easily. Due to its rather recent origin, the method is under developing and possibly not well known to the research community of computational mechanics. Therefore, this section is written in a pedagogical manner to familiarize the readers with the QEM.

###### Lagrange-Type Bar Element.

For homogeneous isotropic materials, the kinematic energy (T) and strain energy (V) of a uniform bar element with length L and cross-sectional area A are given byDisplay Formula

(1)${T=12∫0LρA(∂u∂t)2dx=12(LρA2)∫−11(∂u∂t)2dξ V=12∫0LEA(∂u∂x)2dx=12(LEA2)∫−11(2L∂u∂ξ)2dξ$

where u(x,t) is the axial displacement, $ρ$ is the mass density, and E is the elastic modulus, respectively.

Denote $ξj∈[−1,1]$, the N node points that include the two end nodes and $xj=L(1+ξj)/2$. An N-node (N = 5) bar element is schematically shown in Fig. 1.

The displacement of the bar element is assumed asDisplay Formula

(2)$u(x,t)=∑j=1Nlj(x)u(xj,t)=∑j=1Nlj(ξ)uj(t)$

where the shape function $lj(ξ)$ is Lagrange interpolation function, defined byDisplay Formula

(3)$lj(ξ)=(ξ−ξ1)(ξ−ξ2)...(ξ−ξj−1)(ξ−ξj+1)...(ξ−ξN)(ξj−ξ1)(ξj−ξ2)...(ξj−ξj−1)(ξj−ξj+1)...(ξj−ξN)=∏k=1k≠jNξ−ξkξj−ξk$

Similar to the finite element method, the essential boundary condition can be easily applied since u is the nodal degree-of-freedom (DOF). The bar element belongs to the element of Lagrange type. Substituting Eq. (2) into Eq. (1) and performing the integration numerically yield [5]Display Formula

(4)$T=12{u˙}T[m]{u˙},V=12{u}T[k]{u}$

where the over dot denotes the first-order derivative with respect to time t, and the displacement vector ${u(t)}T=⌊u1(t),u2(t),...,uN(t)⌋$.

Note that $∑j=1Nlj(ξ)=1$. Thus, the diagonal mass matrix can be obtained by the row summation technique. The nonzero diagonal terms are given byDisplay Formula

(5)$mii=LρA2∑k=1MHkli(ξ¯k) (i=1,2,...,N)$

where $ξ¯k$ and $Hk$ are abscissas and weights of the quadrature, and $M$ ($M≤N$) is the number of integration points. If GLL quadrature is used, M = N and $li(ξ¯k)=δki$; thus, the mass matrix is the same as the one without using the row summation technique. In other words, the under-integrated mass matrix by using GLL quadrature equals to the resultant mass matrix after performing the row summation on the fully integrated mass matrix.

The stiffness matrix can be explicitly given byDisplay Formula

(6)$kij=2EAL∑k=1NHkAikAkj (i,j=1,2,...,N)$

where $Aik$ and $Akj$ are the weighting coefficients of the first-order derivative with respect to $ξ$.

If the element nodes are also the abscissas ($ξ¯i=ξi$), such as the element with GLL nodes and integration by using GLL quadrature, the explicit formulas to compute the weighting coefficients are available in the DQM and given byDisplay Formula

(7)$Aij={∏k=1k≠i,jN(ξ¯i−ξk)/∏k=1k≠jN(ξj−ξk)(i≠j)∏k=1k≠iN1(ξ¯i−ξk)(i=j) (i,j=1,2,...,N)$

If the element nodes are not coincided with the abscissas, i.e., $ξ¯i≠ξi$, then Gauss quadrature should be used for accuracy and computational efficiency considerations. The explicit expression for the assumed displacement such as Eq. (2) is required and the explicit formulas in the DQM cannot be directly employed. The novel method presented in Ref. [13] should be used to express kij explicitly. Aij are computed differently byDisplay Formula

(8)$Aij=dlj(ξ)dξ|ξ=ξ¯i=∑k=1NA¯iklj(ξ¯k)=∑k=1NA¯iklkj (i,j=1,2,...,N)$

whereDisplay Formula

(9)$A¯ij={∏k=1k≠i,jN(ξ¯i−ξ¯k) /∏k=1k≠jN(ξ¯j−ξ¯k)(i≠j)∏k=1k≠iN1(ξ¯i−ξ¯k)(i=j) (i,j=1,2,...,N)$

It is easy to show by either Eq. (9) or (8) that $Aij=A¯ij$ when $ξ¯i=ξi$ (i = 1,2,…,N), since $lkj=δkj$.

Expressing the stiffness matrix explicitly with the DQ law is regarded as one of the key differences between the QEM and the SEM. The explicit formula makes the implementation of an N-node quadrature element with ease, although the stiffness and mass matrices of the QEM and SEM are exactly the same if the nodes and integration method are exactly the same. Besides, instead of Eq. (6), Alk can be also computed by using harmonic differential quadrature law [21] asDisplay Formula

(10)$Alk={p∏i=1i≠lkN sin[p(ξ¯l−ξi)] /∏i=1i≠kN sin[p(ξ¯k−ξi)](l≠k)p∏i=1i≠lN cos[p(ξ¯l−ξi)] sin[p(ξ¯l−ξi)](l=k)$

where $p≤0.5$ generally and could be 0.25, 0.5 [22], or $π/8$ [5]. This is implicitly to use the shape function defined asDisplay Formula

(11)$lj(ξ)=∏k=1k≠jN sin[p(ξ−ξk)] sin[p(ξj−ξk)]$

It is seen that formulations of a Lagrange-type weak form quadrature element is flexible and more convenient than the ones of the SEM.

The procedures to formulate two-dimensional (2D) and three-dimensional (3D) weak form quadrature elements are similar to the ones of the SEM. The approximated displacement field of a 2D or 3D quadrature element is assumed by using a tensor product of the displacement function of the one-dimensional (1D) element. If the shape of the weak form quadrature element is irregular, the technique of formulating the subparametric element in conventional FEM should be adopted. The order of the geometric functions is usually much lower than the order of the displacement functions.

###### Hermite-Type Euler–Bernoulli Beam Element.

For homogeneous isotropic materials, the kinematic energy (T) and strain energy (V) of a uniform Euler–Bernoulli beam element with length L, cross-sectional area A, and bending rigidity EI is given byDisplay Formula

(12)${T=12∫0LρA(∂w∂t)2dx=12(LρA2)∫−11(∂w∂t)2dξ V=12∫0LEI(∂2w∂x2)2dx=12(LEI2)∫−11(∂2w∂x2)2dξ$

where w(x,t) is the transverse displacement, $ρ$ is the mass density, and E is the elastic modulus, respectively.

Denote $ξj∈[−1,1]$ the N node points which include the two end nodes and $xj=L(1+ξj)/2$. An N-node (N = 5) beam element is schematically shown in Fig. 2.

The displacement of the beam element is assumed as Display Formula

(13)$w(x,t)=∑j=1Nφj(x)w(xj,t)+ψ1(x)wx(x1,t)+ψN(x)wx(xN,t)=∑j=1N+2hjδj(t)$

where subscript x represents the first-order derivative with respect to x, $δj(t)=w(xj,t)(j=1,2,...,N)$, $δN+1(t)=wx(x1,t)$, $δN+2(t)=wx(xN,t)$, and the shape function $hj(x)$ is the Hermite interpolation function, defined by [23]Display Formula

(14a)$hj(x)=φj(x)={1xj−xN−j+1lj(x)(x−xN−j+1)−[lj(1)(xj)+1xj−xN−j+1]ψj(x)(j=1,N)1(xj−x1)(xj−xN)lj(x)(x−x1)(x−xN)(j=2,3,…,N−1)$
Display Formula
(14b)${hN+1(x)=ψ1(x)=1(x1−xN)l1(x)(x−x1)(x−xN)hN+2(x)=ψN(x)=1(xN−x1)lN(x)(x−xN)(x−x1)$
Display Formula
(14c)$l1(1)(x1)=∑m=2N1x1−xm,lN(1)(xN)=∑m=1N−11xN−xm$

where $lj(x)=lj(ξ)$ is the Lagrange interpolation function defined by Eq. (3) and superscripts (1) represent the first-order derivative with respect to x.

Similar to the finite element method, the essential boundary conditions can be easily applied since w and wx are nodal degrees-of-freedom (DOFs). The weak form quadrature beam element belongs to the element of Hermite type and is different from the existing spectral elements.

Substituting Eq. (13) into Eq. (12) and performing the integration numerically yield [5]Display Formula

(15)$T=12{δ˙}T[m]{δ˙},V=12{δ}T[k]{δ}$
where the over dot denotes the first-order derivative with respect to time t.

The stiffness matrix can be explicitly given byDisplay Formula

(16)$kij=LEI2∑k=1NHkBikBkj (i,j=1,2,...,N)$
where $Bik$ and $Bkj$ are the weighting coefficients of the second-order derivative with respect to x.

Gauss quadrature should be used to obtain a fully integrated stiffness matrix. The explicit expression for the assumed displacement such Eq. (13) is required and the explicit formulas of the DQM cannot be directly used. The novel method presented in Ref. [14] should be used to express kij explicitly. Bij is computed byDisplay Formula

(17a)$Bi(N+1)=1x1−xN{B̃i1(xi−x1)(xi−xN)+2Ãi1[(xi−x1)+(xi−xN)]+2δi1}=ψ1(2)(xi) (i=1,2,...,N)$
Display Formula
(17b)$Bi(N+2)=1xN−x1{B̃iN(xi−x1)(xi−xN)+2ÃiN[(xi−x1)+(xi−xN)]+2δiN}=ψN(2)(xi) (i=1,2,...,N)$
Display Formula
(17c)$Bij=1xj−xN−j+1[B̃ij(xi−xN−j+1)+2Ãij]−[lj(1)(xj)+1xj−xN−j+1]ψj(2)(xi) (j=1,N;i=1,2,...,N)$
Display Formula
(17d)$Bij=1(xj−x1)(xj−xN){B̃ij(xi−x1)(xi−xN)+2Ãij[(xi−x1)+(xi−xN)]+2δij} (j=2,3,...,N−1;i=1,2,...,N)$
where superscripts (1) and (2) represent the first- and second-order derivatives with respect to x, and $lj(1)(xj)$ is computed by Eq. (14c).

The weighting coefficients $Ãij$ and $B̃ij$ are computed byDisplay Formula

(18)${Ãij=2L∑k=1NA¯iklj(ξ¯k)=2L∑k=1NA¯iklkj B̃ij=4L2∑k=1NB¯iklj(ξ¯k)=4L2∑k=1NB¯iklkj (i,j=1,2,...,N)$

where $A¯ik$ is computed by Eq. (9) and $B¯ik=∑j=1NA¯ijA¯jk (i,k=1,2,...,N)$.

If the element nodes are also the abscissas ($ξ¯i=ξi$), such as the element with GLL nodes and integration by the GLL quadrature, the explicit formulas to compute the weighting coefficients, similar to Eq. (17), are available in the DQM but the stiffness matrix is under-integrated, since GLL quadrature is only accurate to $x(2N−3)$.

The consistent mass matrix is a full matrix and given byDisplay Formula

(19)$mij=ρAL2∑k=1N+2Hkhi(ξ¯k)hj(ξ¯k) (i,j=1,2,...,N+2)$

If a diagonal mass matrix is need, a different displacement function from Eq. (13) can be assumed to compute the kinematic energy (T), namely,Display Formula

(20)$w(x,t)=∑j=1Nlj(x)w(xj,t)=∑j=1Nlj(ξ)wj(t)$

The diagonal terms of the mass matrix are given either by Eq. (5) or byDisplay Formula

(21)$mii=LρA2Hi (i=1,2,...,N)$

where $Hi$ is the weight of GLL quadrature and Eq. (21) is exactly the same as Eq. (5). Note that $mii=0$ for $i=(N+1)$ and $(N+2)$.

###### Mixed Lagrange–Hermite-Type Kirchhoff Plate Element.

An elastic isotropic thin rectangular plate element with uniform thickness is considered. For a 2D rectangular plate element, the displacement can be assumed as the tensor product of Eq. (13), namely,Display Formula

(22)$w(x,y,t)=∑i=1N+2∑j=1N+2hi(x)hj(y)δij(t)$

At each corner point, four degrees-of-freedom are used, i.e., $w,wx,wy,wxy$, where subscripts x, y, and xy represent the first-order derivatives with respect to x, y, and the mixed second-order derivative with respect to x and y. Similar to the conventional finite element method, the introduction of the mixed second-order derivative degree-of-freedom ($wxy$) is not convenient for applying the boundary condition and thus it should be removed in formulating the weak form quadrature rectangular element. The degree-of-freedom $wxy$ can be easily removed by using the method of the mixed Lagrange–Hermite interpolations.

Denote a, b the side lengths of the plate and h is the plate thickness, and $E,ν,ρ$ are the elasticity modulus, Poisson's ratio, and the mass density, respectively. The strain energy of the rectangular plate element can be written asDisplay Formula

(23)$U=ab8∫−11∫−11{κ}T[D]{κ}dξdη$

where $w(x,y,t)$ is the deflection, $[D]$ is a 3 × 3 stiffness matrix of the material defined by Display Formula

(24)$[D]=Eh312(1−ν2)[1ν0ν1000(1−ν)/2]$

and ${κ}$ is the curvature vector defined byDisplay Formula

(25)${κ}T= ⌊∂2w∂x2∂2w∂y22∂2w∂x∂y ⌋= ⌊wxxwyy2wxy ⌋$

The kinetic energy of the isotropic rectangular plate element is given byDisplay Formula

(26)$T=ρhab8∫−11∫−11(∂w(x,y,t)∂t)2dξdη$

where $t$ is the time.

An N × N-node (N = 7) thin rectangular plate element in bending is schematically shown in Fig. 3. According to the five criteria for selection of displacement functions [24], three slightly different displacement functions can be assumed for the thin rectangular plate element in bending, namely [25],

Display Formula

(27)$w(x,y,t)=∑i=1N∑j=1Nli(x)lj(y)w(xi,yj,t)=∑i=1N∑j=1Nli(x)lj(y)wij(t)$
Display Formula
(28)$w(x,y,t)=∑i=1N+2∑j=1Nhi(x)lj(y)w⌢ij(t)$
Display Formula
(29)$w(x,y,t)=∑i=1N∑j=1N+2li(x)hj(y)w̃ij(t)$

where N is the number of nodal points in x and y directions, $li(x)$ and $lj(y)$ are Lagrange interpolation functions, $hi(x)$ and $hj(y)$ are Hermite interpolation functions. In Eq. (28), $w⌢ij(t)$ contains the nodal deflection $wij(t)$ (i, j = 1, 2,…, N) as well as the first-order derivative with respect to $x$ at nodes on edges $x=∓a/2$, i.e., $(wx)i1$ and $(wx)iN$ (i = 1, 2,…, N). In Eq. (29), $w̃ij(t)$ contain the nodal deflection $wij(t)$ (i, j = 1, 2,…, N) and the first-order derivative with respect to $y$ at nodes on edges $y=∓b/2$, i.e., $(wy)1j$ and $(wy)Nj$ (j = 1, 2,…, N). The plate element in bending contains $N2+4N$ degrees-of-freedom (DOFs), one DOF at all inner nodes, two DOFs at all boundary nodes, and three DOFs at the four corner nodes, shown in Fig. 3. It is seen that the mixed second-order derivative $wxy$ is not a degree-of-freedom at the plate corner point by using any of the three assumed deflections listed in Eqs. (27)(29).

To satisfy the five criteria for selection of displacement functions [24], one choice of the displacement function is that all three forms are used in formulations of the stiffness and mass matrices [25]. Equation (27) is used to calculate $wxy$ and the kinematic energy, Eq. (28) is used to compute $wxx$, and Eq. (29) is used to compute $wyy$. After integration by using Gauss quadrature, the stiffness matrix can be expressed as [5]Display Formula

(30)$[k]=ab4∑i=1N∑j=1NHiHj[B(x¯i,y¯j)]T[D][B(x¯i,y¯j)]$

where $Hi$ and $Hj$ are weights of Gauss quadrature. The strain matrix at any integration point $(x¯i,y¯j)$, i.e., $[B(x¯i,y¯j)]$, is given byDisplay Formula

(31)$[B(x¯i,y¯j)]{w¯}=[∑l=1N+2∑k=1NBilxljkyw¯lk∑l=1N∑k=1N+2lilxBjkyw¯lk2∑l=1N∑k=1NAilxAjkyw¯lk] (i,j=1,2,...,N)$

where $(x¯i,y¯j)=(aξ¯i/2,bη¯j/2)$, $ξ¯i,η¯j∈[−1,1]$ are abscissas of Gauss quadrature.

In Eq. (31), $lilx=ll(x¯i)$, $ljky=lk(y¯j)$. The degrees-of-freedom, $w¯lk$, contain the deflection at all nodes and the first-order derivative with respect to x at nodes on edges $x=∓a/2$ and the first-order derivative with respect to y at nodes on edges $y=∓b/2$, i.e., $wlk$, $(wx)l1$, $(wx)lN$, $(wy)1k$ and $(wy)Nk$, $(l,k=1,2,...,N)$. $Ailx$ and $Ajky$ are the weighting coefficients of the first-order derivative with respect to x and y, which can be computed by $Ailx=2Ail/a$ and $Ajky=2Ajk/b$, where $Ail$ and $Ajk$ are calculated by Eq. (8). $Bilx$ and $Bjky$ are the weighting coefficients of the second-order derivative with respect to x and y, which can be calculated by Eq. (17).

Substituting Eq. (27) into Eq. (26) and then performing numerical integration by using GLL quadrature results in a diagonal mass matrix. The nonzero diagonal terms $mII$ in the mass matrix areDisplay Formula

(32)$mII=ρhab4HiHj (i,j=1,2,...,N;I=N×(i−1)+j)$

where $Hi(i=1,2,...,N)$ and $Hj(j=1,2,...,N)$ are the weights of GLL quadrature.

It is seen that the formulation of an $N×N$ -node weak form rectangular plate element in bending is simple with the help of the DQ law. The novel method [5,13,14] should be employed to use the DQ law when the nodes are not coincided with integration points. Similar to the finite element method, the essential boundary conditions can be easily applied since w, wx, and wy are nodal degrees-of-freedom (DOFs).

When the element nodes are GLL points and the stiffness matrix is obtained by using GLL quadrature, Eq. (22) can be used to compute both $wxx$ and $wyy$, while Eq. (27) is still used to calculate $wxy$ to eliminate the degree-of-freedom $wxy$ at corner points. Exactly, the same stiffness matrix is obtained as the ones by using Eq. (28) to compute $wxx$ and Eq. (29) to compute $wyy$. Since $hi(xl)=hj(yk)=0(i,j=N+1,N+2)$, $ljky=δkj$, $lilx=δli$, $hi(xl)=δli$, and $hj(yk)=δkj$, the strain matrix at any integration becomesDisplay Formula

(33)$[B(x¯i,y¯j)]{w¯}=[∑l=1N+2Bilxw¯lj∑k=1N+2Bjkyw¯ik2∑l=1N∑k=1NAilxAjkyw¯lk] (i,j=1,2,...,N)$

This implies that different displacement functions yield the same stiffness matrix. However, the stiffness matrix is under-integrated by using GLL quadrature. If a full integration is used, the stiffness matrix is obviously different. Besides, the nodal degrees-of-freedom at corner points are also not the same since $wxy$ cannot be eliminated by using Eq. (22) to compute both $wxx$ and $wyy$.

Note that the formulations are quite different from the ones in Refs. [11,12]. In Ref. [11], Eq. (27) is used to calculate the strain energy, and the derivative degrees are introduced by using the DQ law. In Ref. [12], Eq. (22) is used to compute the strain energy and the mixed second-order derivative is computed approximately to eliminate the degree-of-freedom $wxy$ at the corner point. More details will be given in Sec. 5.

For formulating a weak form quadrature quadrilateral plate element in bending, the subparametric technique should be used to derive the stiffness matrix. Figure 4 schematically shows an N × N-node curved plate element in bending (N = 7) in both computational and physical domains. Detailed procedures can be found in Refs. [5,11], and [12] and many other papers listed in the references.

###### Various Types of Nodes.

There are several types of elemental nodes available. The widely used elemental nodes in spectral element method are the GLL points and Chebyshev–Gauss–Lobatto (CGL) nodes. For the QEM, other types of nodes such as the EC points [17,18] or very good approximations to the Lebesgue-optimal grid points [19] have been used as elemental nodes, since they allow the longest time step during explicit time integration in dynamic analysis. To distinguish them, they are called optimal (EC) node-1 and optimal node-2 in this paper. Other nodes may be also used, such as equidistant nodes [79] and grid V nodes [5]. If equidistant nodes are used, the number of nodes should be small and cannot be arbitrary large, however.

• (a)GLL nodesThe GLL points $ξi(i=1,2,...,N)$ are the roots of the following equation:Display Formula
(34)$(1−ξ2)dPN−1(ξ)dξ=0$
where $PN−1(ξ)$ is the (N − 1)th order Legendre polynomial. Note that $ξi∈[−1,1]$.
The weights corresponding to $ξi$ in the GLL quadrature, Hi, are calculated byDisplay Formula
(35)$Hi=2N(N−1) (i=1,N), Hi=2N(N−1)[PN−1(ξi)]2 (i=2,3,...,N−1)$
The GLL points and corresponding weights for N < 22 can be found in Ref. [5]. A short program should be written to calculate the GLL points and corresponding weights if N > 21. Since GLL nodes and GLL quadrature are widely used in the QEM, a fortran subroutine is included in the Appendix. The subroutine is also converted to the matlab function.
• (b)CGL nodesThe CGL points, also called the Chebyshev extrema [19], can be conveniently computed byDisplay Formula
(36)$ξk=−cos ((k−1)π/(N−1)) (k=1,2,...N)$
• (c)Optimal (EC) node-1Optimal (EC) node-1, the extended Chebyshev nodes [18], can be conveniently calculated byDisplay Formula
(37)$ξk=− cos((2k−1)π/(2N)) cos(π/(2N)) (k=1,2,...N)$
• (d)Optimal node-2Optimal node-2, also the modification of the Chebyshev nodes [19], can be calculated byDisplay Formula
(38)${ξN=−ξ1=1.0 ξk=− cos((2k−1)π/(2N)) cos(π/(2N)(1+0.25/ log(N))) (k=2,3,...N−1)$
• (e)Grid V nodesGrid V is the most reliable grids used in the DQM [5]. Thus, Grid V may be one type of nodes for the weak form quadrature element. Grid V is defined asDisplay Formula
(39)${ξ1=−1.0;ξN=1.0; ξj+1=ξ¯j (j=1,2,...,N−2)$
where $ξ¯j$ is the abscissa of Gauss quadrature.
As was mentioned before, a short program is available [5,20] to compute abscissas and weights of Gauss quadrature.
• (f)Equidistant nodesEquidistant nodes are widely used in conventional finite element method. The nodes can be calculated byDisplay Formula
(40)$ξk=−1+2(k−1)(N−1) (k=1,2,...N)$

Figure 5 shows the Lagrange shape functions (Lagrange interpolation functions) based on the six types of nodes. Figure 6 shows the Hermite shape functions (Hermite interpolation functions) based on the six types of nodes. The maximum and minimum values are listed in Table 1.

According to the criterion of the cardinal-function-minimizing (CFM) optimality [17,18], the GLL nodes-based Lagrange interpolation function is the best one and the equidistant nodes-based Lagrange interpolation function is the worst one. Optimal node-1 and optimal node-2-based Lagrange and Hermite interpolation functions work equally well. It is interesting to see that the equidistant nodes-based Hermite interpolation function is the best one, and Grid V, which works very well in the DQM, seems not a good choice for element nodes when the shape function is Hermite interpolation function. This finding was never reported before and may explain the conclusion given by Striz et al. [7] that the node type does not seem to be a crucial factor for the solution accuracy for problems of thin plates in bending.

It is pointed out by Brutman [17] and Boyd [18] that, however, CFM optimality is not a very discriminating criterion for selection of the nodal points.

## Similarities and Differences Between QEM and SEM

In literatures, the QEM is often regarded as the SEM or spectral finite element method [15,2634]. Thus, the spectral element method is briefly reviewed in this section. The similarities and differences between the QEM and the SEM are given.

According to the definition of the QEM given in this paper, the method presented in Refs. [2634] belongs to the QEM, since the DQ law is used during numerical integration to obtain the element stiffness matrix. Thus, they will be reviewed in Sec. 5. Although a vast body of the SEM paper is available, only a few of them are reviewed to give remarks on the similarities and differences between the QEM and the SEM, since this paper is focused on the QEM.

The spectral element method, originated by Patera [35] in 1984, combines the generality of the FEM with the accuracy of the pseudospectral method (PSM). The SEM is originally applied in the field of fluid dynamics [35]. Priolo and coworkers [36,37] first introduced the SEM to the wave field in 1994. Since then, the method has been gained popularity for simulating wave propagation due to delivering low numerical dispersion with respect to conventional finite element methods [38]. Now, the SEM has been widely applied in fields of fluid mechanics, heat transfer, and acoustics [3941]. Its accuracy and efficiency have been demonstrated in many papers.

Currently, there are two types of spectral elements, called Lobatto elements [4249] and Chebyshev spectral elements [35,36,39,5052]. Lobatto elements use GLL nodes and Chebyshev spectral elements use CGL nodes. The stiffness and mass matrices of Lobatto elements are almost uniquely obtained by using GLL quadrature, and thus the mass matrix is in a diagonal form [40]. Hybrid formulations are also used in obtaining the stiffness matrix [47,48], namely, GLL quadrature is used in integration of the in-plane variables and Gauss quadrature is used in the thickness direction to raise the accuracy of numerical integration. The mass matrix of Chebyshev spectral elements is a full matrix and can be diagonalized by using the row summation technique. With a diagonal mass matrix, a crucial reduction of the complexity and the cost of the numerical time integration can be achieved in simulating wave propagations.

It is first noted by Witkowski et al. [49] that the mass matrix of Lobatto elements is under-integrated. As was mentioned before, the diagonal mass matrix is actually equivalent to the matrix obtained by applying the row summation technique on the fully integrated mass matrix. For two-dimensional (2D) or three-dimensional (3D) Lobatto elements, the stiffness is also under-integrated by using GLL quadrature. To the best of the authors' knowledge, the notions of full or reduced integrations for formulations of the stiffness matrix of Lobatto elements by using GLL quadrature are not as clear as they are for Gauss quadrature.

Based on the contents presented in Sec. 2 and the brief review of the SEM, the similarities and differences between the QEM and SEM (in structural mechanics area) can be summarized as follows.

###### Similarities to the SEM

• (a)Both the weak form quadrature elements and spectral elements are high-order finite elements and formulated based on the variational principle.
• (b)For Lobatto and Chebyshev-type elements, the stiffness and mass matrices are exactly the same if the same quadrature method is used.
• (c)For 2D and 3D problems, the approximation space is built using a tensor product of orthogonal functions by both methods. Usually, subparametric elements are formulated, since the order of the displacement functions is much higher than the one of the geometric functions.

###### Major Differences From the SEM

• (a) The stiffness matrix of the QEM is written explicitly by using the formulas of weighting coefficients in the DQM. Thus, the implementation of a weak form quadrature element with a variable nodal number is very easy. There is no certain procedure to compute the stiffness by the SEM since several ways can be used. Even the spectral element with GLL nodes and integration by GLL quadrature, the explicit formulas for stiffness matrix is not given. Instead, the diagonal mass matrix is emphasized exclusively [40,45]. Sometimes, the stiffness matrix is not even used; instead, the dynamic equilibrium equations are written in terms of the vectors of inertia forces and internal forces [49].
• (b) The shape function of spectral elements is exclusively Lagrange interpolation function, while the shape functions of weak form quadrature elements can be either Lagrange or Hermite interpolation functions or mixed Lagrange–Hermite interpolation functions. For the thin quadrature plate element with GLL nodes and integration by GLL quadrature, the assumed displacement functions may not be uniquely given. Different displacement functions can yield the same stiffness matrix, since the stiffness matrix is under-integrated. This has been demonstrated in Sec. 2 by Eq. (33).
• (c)Besides the GLL and CGL nodes commonly used in SEM, other types of nodes, such as the approximation to the Lebesgue-optimal grid points [18,19], can be used as elemental nodes for the QEM. Two sets of grid points are given by Eqs. (37) and (38). These nodes allow the longest time step during explicit time integrations [18,31].
• (d)Lagrange interpolants for Chebyshev spectral elements, $lj(ξ)$ (j = 0,1,…, n), are defined differently from Eq. (3) and given by [35,36,39,5052]Display Formula
(41)$lj(ξ)=2n∑k=0n1cjckTk(ξj)Tk(ξ) (−1≤ξ≤1)$
where $cj=1 (j≠0,n),cj=2 (j=0,n)$, $Tk(ξ)$ is Chebyshev polynomial of the first kind with order k, and n = N − 1 in which N is the total number o