Manifold-based mass lumping addressing the inertia representation for rotational/torsional degrees of freedom in Kirchhoff plate vibration analysis

Abstract

The vibration analysis of Kirchhoff plates requires robust mass lumping schemes to guarantee numerical stability and accuracy. However, existing methods fail to generate symmetric and positive definite mass matrices when handling rotational degrees of freedom, leading to compromised performance in both time and frequency domains analyses. This study proposes a manifold-based mass lumping scheme that systematically resolves the inertia matrix formulas for rotational/torsional degrees of freedom. By reinterpreting the finite element mesh as a mathematical cover composed of overlapping patches, Hermitian interpolations for plate deflection are derived using partition of unity principles. The manifold-based mass matrix is constructed by integrating the virtual work of inertia forces over these patches, ensuring symmetry and positive definiteness. Numerical benchmarks demonstrate that the manifold-based mass lumping scheme performance can be comparable or better than the consistent mass scheme and other existing mass lumping schemes. This work establishes a unified theory for mass lumping in fourth order plate dynamics, proving that the widely used row-sum method is a special case of the manifold-based framework. The scheme resolves long-standing limitations in rotational/torsional inertia conservation and provides a foundation for extending rigorous mass lumping to 3D shell and nonlinear dynamic analyses.

Graphical abstract

Keywords

Kirchhoff plates / manifold-based mass lumping / symmetric positive definite matrices / rotational/torsional inertia conservation / partition of unity / finite element vibration analysis

1 Introduction

Static, vibration and buckling analyses of plate structures play a pivotal role, because such structures are ubiquitous in engineering applications. Significant effort has been focused on the finite element method (FEM) [1] and its extensions [2], the boundary element method [3,4], and the meshfree method [5,6]. While FEM remains dominant for thin plate bending governed by Kirchhoff theory, the treatment of rotational degrees of freedom (DOFs) introduces fundamental challenges in mass matrix formulation. This article concerns vibration analysis of thin plate bending problems using the FEM in the presence of rotational/torsional DOFs.

Since proper representation of inertia properties can significantly reduce computational cost and yield reliably accurate results, different ad hoc mass representation schemes have been proposed and applied in vibration analysis [7–10]. The major two categories of mass matrices widely used in dynamic analysis of FEM are the consistent mass matrix (CMM) and diagonally lumped mass matrix (LMM) [11]. The template mass lumping technique were then proposed using algebraic forms with free parameters, and mass matrices are constructed to satisfy certain conditions [9,10,12] that are stated more intuitively rather than mathematically or physically.

Archer made large strides by formulating CMM in structural mechanics [13,14]. The CMM, derived using the same shape functions as the stiffness matrix, ensures optimal convergence properties but is inherently non-diagonal. CMMs have the following drawbacks [10]. First, in both the explicit direct integration method for transient response analysis and the eigenvalue problem for the free vibration analysis, using CMMs results in lack of efficiency. Second, in scenarios such as contact-impact analysis and wave propagation dynamics, CMMs have been found to result in spurious oscillations in the solution even if an implicit time integration scheme is applied. As a result, a LMM is generally preferable.

For the LMM, the total mass of an element is directly apportioned to the DOFs each nodal regardless of the inter-element coupling. The merits of diagonally LMM generally lie in computational efficiency, with less storage space and processing time. For the analysis of Kirchhoff plate bending, two mass lumping schemes are commonly used, namely, the row-sum method [15] and the diagonal scaling procedure [16].

As pointed out in [10], obtaining an LMM is more difficult. Since there is no straightforward approach to obtain rotational/torsional masses, a LMM derived by approaches such as the row-sum method has zero diagonal entries corresponding to rotational/torsional DOFs. The LMM that has lost positive definiteness cannot be employed directly in applications: 1) negative or zero components in lumped mass matrices lead to a divergence of explicit time integration schemes and the computation of zero or infinite eigenfrequencies in modal analysis [17]; 2) the singular LMM with zero entries is impossible to be used with Newmark method, and instead an implicit-explicit procedure is generally implemented to tackle this issue [18]; and 3) the singular LMM cripples the subspace iteration procedure as eigenvalue solvers, and to treat this, a static condensation process must be operated in advance [19]. Voet et al. [20] proposed geometry-aware mass lumping to suppress high-frequency outliers, enhancing CFL stability in explicit dynamics. However, their method targets second order systems. Consequently, the full potential of explicit and implicit time integrators cannot be tapped and the whole implementation becomes more challenging. To preserve the positive definiteness of LMMs for plate elements with rotational DOFs, Archer and Whalen developed a rationally consistent diagonal mass matrix by enforcing equality of the rotational inertias between the CMM and LMM [21], however the definition of the rotational inertia is based on intuition.The diagonal scaling procedure [16] was developed to generate positive definiteness of lumping matrices. Unfortunately, as is pointed out by Hughes [15], there is no strict mathematical theory to prop up this procedure and is likely to lead to reduced convergence rate [17].

Recent advances in mass lumping for second order problems, including template-based methods [10] and isogeometric dual-basis approaches [22], have improved efficiency but lack mathematical rigor when extended to fourth order plate dynamics. For instance, recent theoretical advancements in mass lumping, notably the algebraic matrix framework by Voet et al. [23], decouple diagonalization from Galerkin discretization through non-diagonal mass matrices constructed from Lagrange bases. Geevers and Maier [24] demonstrated second order convergence in energy norms using localized orthogonal decomposition. But their framework remains confined to second order systems. The dual-basis mass lumping method proposed by Hiemstra et al. [25] achieves high-order accuracy in second order systems but remains untested for fourth order problems with rotational/torsional DOFs. Similarly, Li et al. [26] developed an interpolatory basis method using Greville nodes for Kirchhoff plates, achieving fourth order frequency convergence with quartic splines. However, their displacement-centric formulation neglects rotational inertia conservation in the mass matrix. Notably, even in advanced fourth-order thin shell applications, such as FGM plates [27], composite vibrations [28], and cracked shell dynamics [29], or multiphysics problems like inverse deformation [30], thermal morphing [31], and swelling-driven deformation [32], mass matrices still rely on conventional diagonalization. These methods fail to ensure both positive definiteness and rotational inertia conservation, particularly for rotational/torsional DOFs, limiting their robustness in explicit dynamics and high-order vibration analysis. Recently, a new manifold-based lumping technique for finite element analysis was proposed with rigorous and unified theoretical basis and was proven to be able to generate positive-definite diagonal lumped mass matrices for second order problems, such as higher-order serendipity elements [33,34]. For example, Li et al. [35] developed an improved meshless manifold method (iMNMM) for second order fracture dynamics, achieving stable explicit time integration through aligned background quadrature and manifold-based mass lumping. The spectral decomposition theorem-driven mass lumping [36] and iMNMMs [37] highlight the potential of partition of unity (PU) principles and Cholesky-based block diagonalization. However, for dynamic analysis of fourth order problems with Kirchhoff plate elements containing rotational/torsional DOFs, more steps need to be performed in lumping to generate corresponding positive-definite diagonal mass matrices.

In this study, we continue to apply the manifold based method to lump mass of plate elements, which are also used in any problems of fourth order, where the virtual work done by the inertia force is calculated by means of integrals on manifolds. Different from solid elements in the analysis of second order problems where only translational DOFs are involved in each node, in FEA for plate bending following Kirchhoff’s theory any node has one translational and several rotational/torsional DOFs, causing the mass matrix to be a block lumped matrix {\mathit{M}}_{\mathrm{L}} that is symmetric and positive definite. The order of each diagonal block in {\mathit{M}}_{\mathrm{L}} is small and equals to the number of DOFs of a node. By carrying out a Cholesky decomposition of each diagonal block in {\mathit{M}}_{\mathrm{L}}, consequently, {\mathit{M}}_{\mathrm{L}}^{-1} can be easily obtained with calculations of O(n_{\mathrm{D}\mathrm{O}\mathrm{F}}), still in the same order of magnitudes as the inverse of a real diagonal matrix. Here, n_{\mathrm{D}\mathrm{O}\mathrm{F}} is the number of total DOFs in the mesh. The key contributions are threefold: 1) a mathematical framework addressing fourth order plate dynamics; 2) a block-diagonal mass matrix compatible with explicit/implicit solvers; 3) comprehensive validation across free vibration, transient, and harmonic response benchmarks.

2 Fundamentals on Kirchhoff plate bending

In this section, The weak form of the thin plate bending problem within the Kirchhoff plate theory framework is first derived, with kinematic relationships, governing equations, and strong form following established continuum mechanics principles [38]. Subsequently, the finite element discretization procedure is systematically recalled. While the scheme to be proposed is fit for other theories of plate deformation as well, the selection of the Kirchhoff plate is because more analytical/reference solutions available can be used to verify the effectiveness of the proposed scheme.

2.1 Kirchhoff plates and governing equations

Based on the Kirchhoff thin plate bending theory [1], the relation between the lateral deflection w\left(x,y\right) of the middle surface \left(\mathrm{z}=0\right) and the rotations \left({\theta }_x,{\theta }_y\right) can be given by:

(1)\{\begin{array}{rl}& {\theta }_x=\frac{\mathrm{\partial }w}{\mathrm{\partial }x},\\ & {\theta }_y=\frac{\mathrm{\partial }w}{\mathrm{\partial }y}.\end{array}

Under the coordinate system illustrated in Fig.1, the displacement field in a thin plate can be expressed as:

(2)\{\begin{array}{rl}& u\left(x,y,\mathrm{z}\right)=-\mathrm{z}\frac{\mathrm{\partial }w}{\mathrm{\partial }x},\\ & v\left(x,y,\mathrm{z}\right)=-\mathrm{z}\frac{\mathrm{\partial }w}{\mathrm{\partial }y},\\ & w\left(x,y,\mathrm{z}\right)=w\left(x,y\right).\end{array}

It is obvious that the transversal deflection of the middle plane of the thin plate w can be regarded as the field variables of the bending problem of thin plates.

The corresponding bending and twisting curvatures are the generalized strains expressed as:

(3)\mathit{k}=\left\{\begin{array}{c}{k}_{xx}\\ k_{yy}\\ k_{xy}\end{array}\right\}=-\left\{\begin{array}{c}{\displaystyle \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}}\\ {\displaystyle \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}}\\ 2{\displaystyle \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }x\mathrm{\partial }y}}\end{array}\right\}=\mathit{L}w,

where \mathit{L} being the differential operator defined as \mathit{L}=-{\left\{{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}}{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}}2{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }x\mathrm{\partial }y}}\right\}}^{\mathrm{T}}.

According to the Kirchhoff theory, the bending and twisting moments, depicted in Fig.1 are the generalized stresses defined as:

(4)\{\begin{array}{rl}& M_x=-D_0\left(\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}+\upsilon \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}\right),\\ & M_y=-D_0\left(\upsilon \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}\right),\\ & M_{xy}=M_{yx}=-D_0\left(1-\nu \right)\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }xy}.\end{array}

D_0={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}} is the bending rigidity, where E and \nu are the Young’s modulus and Poisson’s ratio, and h is the thickness of the thin plate.

Accordingly, Hookie’s law for the thin plate formulation can be further expressed in the matrix form as:

(5)\mathit{S}=\mathit{D}\mathit{k},

where {{\mathit{S}=(M}_x,M_y,M_{xy})}^{\mathrm{T}}, \mathit{D}=D_0\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& {\displaystyle \frac{1-\nu }{2}}\end{array}\right], and D_0 is the bending rigidity denoted as

(6)D_0=\frac{E{h}^3}{12(1-{\nu }^2)},

where E being Young’s modulus and \nu being Poisson’s ratio, respectively.

The planar domain of the middle surface is denoted by \mathrm{\Omega }. The Frenet coordinates of the boundaries are marked in Fig.2, with its tangent and normal direction denoted by n and s.

The external unit normal vector n is \mathit{n}={\left(n_x,n_y\right)}^{\mathrm{T}}, and the unit tangent vector is \mathit{s}={\left(-n_y,n_x\right)}^{\mathrm{T}}. The first-order derivatives about the normal and tangent direction can thus be:

(7)\{\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }n}=n_x\frac{\mathrm{\partial }}{\mathrm{\partial }x}+n_y\frac{\mathrm{\partial }}{\mathrm{\partial }y},\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }s}={-n}_y\frac{\mathrm{\partial }}{\mathrm{\partial }x}+n_x\frac{\mathrm{\partial }}{\mathrm{\partial }y}.\end{array}

The generalized strains and the stresses about the normal and tangent direction are written as:

(8)\{\begin{array}{rl}& {\varphi }_{\mathrm{n}}={\varphi }_x{n_x}^2+2{{\varphi }_{xy}{n}_{x}n}_y+{\varphi }_y{n_y}^2,\\ & {\varphi }_{\mathrm{s}}={\varphi }_x{n_y}^2-2{{\varphi }_{xy}{n}_{x}n}_y+{\varphi }_y{n_x}^2,\\ & {\varphi }_{\mathrm{n}\mathrm{s}}=\left({\varphi }_y-{\varphi }_x\right)n_x{n}_y+{\varphi }_{xy}\left({n_x}^2-{n_y}^2\right),\end{array}

where the variant \varphi is defined to represent either curvatures \mathit{k} and moments \mathit{M}.

The lateral shear force can also be defined as: Q_{\mathrm{n}} = Q_x{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}+Q_y{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}}.

Generally speaking, the boundary conditions can be classified into three parts, namely \mathrm{\partial }\mathrm{\Omega }={\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2\cup {\mathrm{\Gamma }}_3.

For clamped edge boundary {\mathrm{\Gamma }}_1, w=\tilde{w},{\displaystyle \frac{\mathrm{\partial }w}{\mathrm{\partial }n}}={\tilde{\theta }}_{\mathrm{n}}. \tilde{w} and {\tilde{\theta }}_{\mathrm{n}} are functions of arc length.

For simply supported edge boundary {\mathrm{\Gamma }}_2, w=\tilde{w}, M_{\mathrm{n}}={\tilde{M}}_{\mathrm{n}}. {\tilde{M}}_{\mathrm{n}} is the function of arc length, too.

For free edge boundary {\mathrm{\Gamma }}_3, M_{\mathrm{n}}={\tilde{M}}_{\mathrm{n}},{\displaystyle \frac{\mathrm{\partial }{M}_{\mathrm{n}\mathrm{s}}}{\mathrm{\partial }s}}+Q_{\mathrm{n}}=\tilde{q}, where \tilde{q} is a linearly distributed load.

2.2 The weak form of Kirchhoff plate bending

The object in this study is a thin plate occupying the domain \mathrm{\Omega } in the xy plane. The traction \overline{p}\left(x,y;t\right)\mathit{n} exerted perpendicular to the xy plane is assumed to be time dependent, where \mathit{n} refers to the unit vector along z-axis.

The deflection w\left(x,y;t\right) of the middle plane is generally taken as the primal variable, and state the problem in the following way.

Given initial deflection and velocity of the middle plane

(9)w\left(x,y;0\right)=\overline{w}\left(x,y\right),

and

(10)\dot{w}\left(x,y;0\right)=\overline{v}\left(x,y\right),

find the function w\left(x,y;t\right), such that the virtual work

(11)\left(\delta w,\ddot{w}\right)+a\left(\delta w,w\right)=l\left(\delta w\right),

holds for w at any time t\in (0,T] and any \delta w\left(x,y\right) in the function space W.

The boundary conditions of plate in investigation are assumed to be simply supported along the boundary {\mathrm{\Gamma }}_{\mathrm{S}}, clamped along the boundary {\mathrm{\Gamma }}_{\mathrm{C}}, and free from loading along free boundary {\mathrm{\Gamma }}_{\mathrm{F}}, which make up the boundary \mathrm{\partial }\mathrm{\Omega }={\mathrm{\Gamma }}_{\mathrm{S}}\cup {\mathrm{\Gamma }}_{\mathrm{C}}\cup {\mathrm{\Gamma }}_{\mathrm{F}} of the plate, but anyone of the three portions can be empty. Thus, W can be defined as

(12)W=\left\{v|v=0,v_{\mathrm{n}}\equiv \frac{\mathrm{\partial }v}{\mathrm{\partial }n}=0,\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{C}};v=0\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{S}}\right\}.

We mention that for any t\in (0,T], w\left(x,y;t\right) subordinates to W, which is a function of spatial coordinates x and y with t as a parameter.

The virtual work done by inertia force \left(\delta w,\ddot{w}\right) in Eq. (11) can be represented as:

(13)\left(\delta w,\ddot{w}\right)={\int }_{\mathrm{\Omega }}^{}\rho \delta w\ddot{w}h\mathrm{d}\mathrm{\Omega },

to be introduced in Section 4, where h=h\left(x,y\right) is the thickness of the plate at \left(x,y\right). The virtual work done by generalized stress \mathit{S} can be represented as:

(14)a\left(\delta w,w\right)={\int }_{\mathrm{\Omega }}^{}{\left(\delta \mathit{\kappa }\right)}^{\mathrm{T}}\mathit{S}\mathrm{d}\mathrm{\Omega },

with \mathit{S} to be denoted shortly, and the virtual work done by \overline{p} can be expressed as:

(15)l\left(\delta w\right)={\int }_{\mathrm{\Omega }}^{}\overline{p}\delta w\mathrm{d}\mathrm{\Omega }.

2.3 Spatial discretization

The mass lumping scheme to be proposed in the following is applicable to any plate element type. However, to avoid unnecessary distractions, we choose the ACM (Adini-Clough-Melosh) element [39] and the BFS (Bogner-Fox-Schmit) element [40] as the representative of incompatible elements and compatible elements, both are rectangles with four nodes. Although the two element types are very classic, it will be seen subsequently that they are able to achieve rather high precision comparable to or even better than those elements developed much later.

As an aside, while the ACM element and the BFS element are rectangular, they can be applied to plates of any shape if the analysis is conducted in the framework of the numerical manifold method (NMM) [2]. In this study, however, we will only concern ourselves about the implementation in FEM and avoid wading in NMM deeply, which is a more powerful and intriguing field.

It is supposed that the plate is discretized into n_e elements and n nodes. In element e, either ACM or BFS element, the Hermitian form approximation to deflection w expressed by

(16)w\left(x,y;t\right)=\sum _{i=1}^4\sum _{j=0}^D{N}_{i,e}^j\left(x,y\right)d_{i,e}^j\left(t\right),

for any \left(x,y\right)\in {\mathrm{\Omega }}^e. Here {\mathrm{\Omega }}^e is the element domain e; the number of DOFs at a node is D_1=D+1, with D_1 = 3 for ACM element and D_1 = 4 for BFS element, respectively. N_{i,e}^j with j=0,\dots ,D is the jth order shape function associated to the ith node of element e; accordingly, d_{i,e}^j defines the jth DOF associated to the ith node of element e, with j=0,\dots ,D. It needs to be noted that only d_{i,e}^0 represents the translational DOF, d_{i,e}^1,…,d_{i,e}^D corresponds to rotational and/or torsional DOFs. Equation (16) can be rewritten in a compact matrix form:

(17)\mathit{w}={\mathit{N}}_e{\mathit{d}}_e.

With {\mathit{d}}_e being the 4D_1-dimensions DOFs for all four nodes of element e:

(18){\mathit{d}}_e^{\mathrm{T}}=\left(\begin{array}{ccc}{\mathit{d}}_{1,e}^{\mathrm{T}};& \cdots ;& {\mathit{d}}_{4,e}^{\mathrm{T}}\end{array}\right).

On the other hand, the D_1-dimensional vector {\mathit{d}}_{k,e} (k=1,…, 4) can be denoted as the DOFs of the kth node of element e, namely,

(19){\mathit{d}}_{k,e}^{\mathrm{T}}=\left(\begin{array}{ccc}{d}_{k,e}^0,& \cdots ,& d_{k,e}^{\mathrm{D}}\end{array}\right).

Likewise, {\mathit{N}}_e defines the shape functions of element e with 4D_1-dimensions:

(20){\mathit{N}}_e=\left(\begin{array}{ccc}{\mathit{N}}_{1,e}& \cdots & {\mathit{N}}_{4,e}\end{array}\right),

and its D_1-dimensional components {\mathit{N}}_{k,e}, k=1,...,4, corresponds to the D_1 shape functions of the kth (k=1,\dots ,4) node of element e as

(21){\mathit{N}}_{k,e}=\left(\begin{array}{ccc}{N}_{k,e}^0,& \cdots ,& N_{k,e}^{\mathrm{D}}\end{array}\right).

The D_1 shape functions can be referred in plate element literatures such as Refs. [39,40]. Here it needs to be emphasized that the zeroth order shape function N_{k,e}^0 having the PU properties as follows:

(22)N_{k,e}^0\left(x_j,y_j\right)={\delta }_{kj},k,j=1,\dots ,4,

and

(23)\sum _{k=1}^4{N}_{k,e}^0\left(x,y\right)=1,\mathrm{\forall }\left(x,y\right)\in {\mathrm{\Omega }}^e,

which will help to construct the manifold-based mass lumping technique for plate element. Specifically, Fig.3 demonstrates the image of N_{1,e}^0\left(x,y\right) at first node of a BFS element, it can be observed that the zeroth order shape function vanish on the edges opposite to kth which ensures the continuity of deflection across element edges. The ACM element is very similar to the BFS element.

Substitution of Eq. (16) into Eq. (11) leads to the semi-discrete system of ordinary differential equations (ODEs),

(24)\mathit{M}\ddot{\mathit{d}}+\mathit{K}\mathit{d}=\mathit{f},

where \mathit{d} is the DOFs vector of all the nodes, \ddot{\mathit{d}} is the corresponding acceleration vector. \mathit{M}, \mathit{K} and \mathit{f} are the global mass matrix, global stiffness matrix and global force vector, respectively, which can be obtained by assembling the corresponding element matrices:

(25)\mathit{M}=\sum _e{\mathit{M}}_{\mathrm{C}}^e,\mathit{K}=\sum _e{\mathit{K}}^e,\mathit{f}=\sum _e{\mathit{f}}^e,

with

(26){\mathit{M}}_{\mathrm{C}}^e={\int }_{{\mathrm{\Omega }}^e}^{}\rho h{\mathit{N}}_e^{\mathrm{T}}{\mathit{N}}_e\mathrm{d}\mathrm{\Omega },

(27){\mathit{K}}^e={\int }_{{\mathrm{\Omega }}^e}^{}{\mathit{B}}_e^{\mathrm{T}}\mathit{D}{\mathit{B}}_e\mathrm{d}\mathrm{\Omega },

and

(28){\mathit{f}}^e={\int }_{{\mathrm{\Omega }}^e}^{}{\overline{p}\mathit{N}}_e^{\mathrm{T}}\mathrm{d}\mathrm{\Omega }.

{\mathit{B}}_e=\mathit{L}{\mathit{N}}_e. The subscript C in the element mass matrix refers to the fact that the mass matrix obtained in this fashion is the CMM. Our approach to instead construct a LMM is presented in Section 4.

In the above deduction, no damping force is considered. In the presence of damping, the system of ODEs turns into

(29)\mathit{M}\ddot{\mathit{d}}+\mathit{C}\dot{\mathit{d}}+\mathit{K}\mathit{d}=\mathit{f},

where the damping matrix \mathit{C} may be determined, for example, by the Rayleigh damping assumption, can be defined as

(30)\mathit{C}={\alpha }_{\mathrm{R}}\mathit{M}+{\beta }_{\mathrm{R}}\mathit{K},

where {\alpha }_{\mathrm{R}} and {\beta }_{\mathrm{R}} are two given constants.

3 Three problem types for plate vibration

In this section, we summarize three problem types commonly encountered in the vibration analysis of plates, which are best suitable for the response analysis of transient forcing functions, the mode analysis of free vibration, and the harmonic response analysis, respectively. All the three problem types are illustrated by examples in Section 5.

3.1 The Newmark method for direct time integration

The initial-value problem associated with ODEs in Eq. (29) consists of finding a deflection vector \mathit{d}\left(t\right), satisfying ODEs in Eq. (29) and the initial data

(31)\mathit{d}\left(0\right)={\mathit{d}}_0,

(32)\dot{\mathit{d}}\left(0\right)={\mathit{v}}_0.

Perhaps the most widely used method for solving Eq. (29) together with Eqs. (31) and (32) is the Newmark method that is stated briefly as follows.

It is supposed that {\mathit{d}}_n, {\mathit{v}}_n, and {\mathit{a}}_n are the approximations of \mathit{d}\left(t_n\right), \dot{\mathit{d}}\left(t_n\right), and \ddot{\mathit{d}}\left(t_n\right), respectively. At t=t_{n+1}=t_n+\mathrm{\Delta }t, writing

(33){\tilde{\mathit{d}}}_{n+1}={\mathit{d}}_n+\mathrm{\Delta }t{\mathit{v}}_n+\frac{\mathrm{\Delta }{t}^2}{2}\left(1-2\beta \right){\mathit{a}}_n,

(34){\tilde{\mathit{v}}}_{n+1}={\mathit{v}}_n+\mathrm{\Delta }t\left(1-\gamma \right){\mathit{a}}_n,

then the approximations of \mathit{d}\left(t_{n+1}\right) and \dot{\mathit{d}}\left(t_{n+1}\right) are

(35){\mathit{d}}_{n+1}={\tilde{\mathit{d}}}_n+\beta \mathrm{\Delta }{t}^2{\mathit{a}}_{n+1},

(36){\mathit{v}}_{n+1}={\tilde{\mathit{v}}}_{n+1}+\gamma \mathrm{\Delta }t{\mathit{a}}_{n+1},

where {\mathit{a}}_{n+1} is the solution to system

(37)\overline{\mathit{M}}{\mathit{a}}_{n+1}={\overline{\mathit{f}}}_{n+1},

with

(38)\overline{\mathit{M}}=\mathit{M}+\gamma \mathrm{\Delta }t\mathit{C}+\beta \mathrm{\Delta }{t}^2\mathit{K},

(39){\overline{\mathit{f}}}_{n+1}={\mathit{f}}_{n+1}-\mathit{C}{\tilde{\mathit{v}}}_{n+1}-\mathit{K}{\tilde{\mathit{d}}}_{n+1}.

To start the above process, {\mathit{a}}_0 is calculated from

(40)\mathit{M}{\mathit{a}}_0={\mathit{f}}_0-\mathit{C}{\mathit{v}}_0-\mathit{K}{\mathit{d}}_0.

To this point, we can see that if \mathit{M} is a DLMM (Diagonally Lumped Mass Matrix), solving {\mathit{a}}_0 requires \mathit{M} to be positive. Therefore, those lumping schemes that produce positive semidefinite DLMMs would encounter difficulties in applying the Newmark method. Furthermore, if damping is ignored, as is done frequently and \beta = 0, then \overline{\mathit{M}} reduces to \mathit{M}. Consequently, if \mathit{M} is a positive LMM or any other matrix whose inverse {\mathit{M}}^{-1} is easy to calculate, solving Eq. (37) for {\mathit{a}}_{n+1} becomes a mere trifle.

Since each node of the finite element mesh for the plate has both translational and rotational DOFs, the LMM \mathit{M} to be proposed subsequently is not a genuine diagonal matrix but a diagonal block matrix, namely,

(41)\mathit{M}={\mathit{M}}_{\mathrm{L}}=\left[\begin{array}{ccc}{\mathit{M}}_{\mathrm{L}}^1& \cdots & 0\\ & \cdots & \\ 0& \cdots & {\mathit{M}}_{\mathrm{L}}^n\end{array}\right],

where each diagonal block {\mathit{M}}_{\mathrm{L}}^i with i=1,\dots ,n is a symmetric positive definite matrix of order D_1; recalling that D_1 is the number of DOFs of a node and n is the number of all the nodes in the mesh. Therefore, the Cholesky decomposition of {\mathit{M}}_{\mathrm{L}}^i,

(42){\mathit{M}}_{\mathrm{L}}^i={\mathit{L}}_i{\mathit{L}}_i^{\mathrm{T}},

can be obtained with very few calculations, where {\mathit{L}}_i is the lower triangular matrix of {\mathit{M}}_{\mathrm{L}}^i.

Certainly we can direct solve ODEs in Eq. (29) with matrix {\mathit{M}}_{\mathrm{L}} replacing \mathit{M}, that is also easy to implement. However, to exploit the diagonal block structure of {\mathit{M}}_{\mathrm{L}} in Eq. (41), we propose to solve the transformed ODEs,

(43)\ddot{\tilde{\mathit{d}}}+\tilde{\mathit{C}}\dot{\tilde{\mathit{d}}}+\tilde{\mathit{K}}\tilde{\mathit{d}}=\tilde{\mathit{f}},

instead, which has much higher efficiency. Here \tilde{\mathit{C}}=\left[{\tilde{\mathit{C}}}_{ij}\right], \tilde{\mathit{K}}=\left[{\tilde{\mathit{K}}}_{ij}\right] and \tilde{\mathit{f}}=\left[{\tilde{\mathit{f}}}_i\right], i,j=1,\dots ,n, are defined block by block as follows

(44){\tilde{\mathit{C}}}_{ij}={\mathit{L}}_i^{-1}{\mathit{C}}_{ij}{\mathit{L}}_j^{-\mathrm{T}},

(45){\tilde{\mathit{K}}}_{ij}={\mathit{L}}_i^{-1}{\mathit{K}}_{ij}{\mathit{L}}_j^{-\mathrm{T}},

and

(46){\tilde{\mathit{f}}}_i={\mathit{L}}_i^{-1}{\mathit{f}}_i.

Once \tilde{\mathit{d}}=\left({\tilde{\mathit{d}}}_i\right) is solved from ODEs in Eq. (43), where {\tilde{\mathit{d}}}_i is the DOF vector of node-i, a D_1-dimensional sub-vector of vector \tilde{\mathit{d}}, we immediately have the solution \mathit{d}=\left({\mathit{d}}_i\right) by

(47){\mathit{d}}_i={\mathit{L}}_i^{-\mathrm{T}}{\tilde{\mathit{d}}}_i.

3.2 The generalized eigenvalue problem arising from free vibration analysis

Once we have the stiffness matrix and the mass matrix, for the same structure we have the following generalized eigenvalue problem

(48)\mathit{K}\varphi ={\omega }^2\mathit{M}\varphi ,

which yields the n_{\mathrm{D}\mathrm{O}\mathrm{F}} eigensolutions \left\{{\omega }_k^2,{\varphi }^{\left(k\right)}\right\}, k=1,…,n_{\mathrm{D}\mathrm{O}\mathrm{F}} arranged in the ascending order of {\omega }_k. Here n_{\mathrm{D}\mathrm{O}\mathrm{F}} is the number of all the DOFs in the system; the value {\omega }_k is the kth order natural frequency, and the lowest value {\omega }_1 is the fundamental natural frequency; the eigenvector {\varphi }^{\left(k\right)} denotes the mode shape vector corresponding to {\omega }_k.

Each pair of \left\{{\omega }_k,{\varphi }^{\left(k\right)}\right\} determines a mode of vibration of the plate. The pairs of \left\{{\omega }_k,{\varphi }^{\left(k\right)}\right\} are important in the mode method and evaluation of dynamic properties. However, finding out all the n_{\mathrm{D}\mathrm{O}\mathrm{F}} modes is an exceedingly time-consuming process. Therefore, only those of lower frequency are found out. The subspace iteration method is perhaps the most widely used solution to find out the first lower order modes.

However, if \mathit{M} is a diagonal block matrix as in Eq. (41), the generalized eigenvalue problem in Eq. (48) reduces to the standard eigenvalue problem,

(49)\tilde{\mathit{K}}\tilde{\varphi }={\omega }^2\tilde{\varphi },

where \tilde{\mathit{K}} is defined in Eq. (45) and \varphi =\left({\varphi }_i\right) is related to \tilde{\varphi }=\left({\tilde{\varphi }}_i\right) as

(50){\varphi }_i={\mathit{L}}_i^{-\mathrm{T}}{\tilde{\varphi }}_i,i=1,\dots ,n.

There are more off-the-shelf algorithms available to as standard eigenvalue solvers in Eq. (49). Moreover, under the same computational efforts and memory with the generalized eigenvalue problem, more eigenpairs can be extracted by solving the standard eigenvalue problems. Those higher order frequency eigenmodes play an important role in the analysis of vibration caused by blast or shock loads [18].

3.3 Harmonic response analysis for linear structures undergoing stead state vibration

The direct time integration procedures, such as the Newmark method, is best for transient forcing functions. If the force amplitude repeats itself regularly many times, that is to say the external load is determined by a periodic forcing function, any initial start-up transient response is not significant and is usually ignored [41]. Instead, the steady-state response that the structure settles down to is of interest. The most efficient solution method appears to be those algorithms in the frequency domain that calculate the steady-state response directly.

Since any periodic forcing function can be decomposed into a series of sine and cosine components by means of a Fourier series, the periodic forcing function can be assumed a cosine or sine function. Once the response of each component in the Fourier series is found, superposition is used to give the total response. And these steps constitute the harmonic response analysis for a linear structure undergoing steady-state vibration.

When the structure is in the steady-state response, all loads and displacements vary sinusoidally at the same frequency although not necessarily in phase. The steady-state response also obeys the system of ODEs in Eq. (29) but the real load vector \mathit{f} with the common factor \mathrm{s}\mathrm{i}\mathrm{n}\omega t among all the components, namely,

(51)\mathit{f}={\mathit{f}}^0\mathrm{s}\mathrm{i}\mathrm{n}\omega t,

is extended to a complex vector defined as

(52)\mathit{f}={\mathit{f}}^0{\mathrm{e}}^{\mathrm{i}\omega t},

where {\mathit{f}}^0 is an n_{\mathrm{D}\mathrm{O}\mathrm{F}}-dimensional vector; \mathrm{i}=\sqrt{-1}; and \omega is the imposed circular frequency (radians/time)

To determine the harmonic response, we first solve the system of complex linear algebraic equations

(53)\mathit{Z}\mathit{d}=\mathit{f},

and then take the imaginary part of solution vector, \mathrm{I}\mathrm{m}\left(\mathit{d}\right), as the displacement response \mathit{d}\left(\omega \right) corresponding to the imposed circular frequency \omega in the sine load vector defined in Eq. (51). Here \mathit{Z}=\mathit{Z}\left(\omega \right) is the impedance matrix at the frequency \phi,

(54)\mathit{Z}\left(\omega \right)=\mathit{K}-{\omega }^2\mathit{M}+\mathrm{i}\omega \mathit{C}.

Similarly, if the mass lumping scheme to be proposed subsequently is applied, then \mathit{K}, \mathit{M}, \mathit{C} and \mathit{f} in Eqs. (53) and (54) are replaced by \tilde{\mathit{K}}, \mathit{I}, \tilde{\mathit{C}}, and \tilde{\mathit{f}}, respectively; and the intermediate solution \tilde{\mathit{d}}\left(\omega \right) is transformed in the displacement response \mathit{d}\left(\omega \right) according to Eq. (47).

4 The manifold based mass lumping procedure

To propose the mass lumping scheme that yields the LMMs that are symmetric positive definite, while calculating the virtual work done by inertia force we first regard the finite element mesh as a cover of the plate which is made up of a series of subdomains called patches. Then, we view the Hermitian interpolation to the plate deflection from the perspective of the PU to retrieve the PU functions and the local approximations subordinate to the patches. This discussion about recovering PU functions from Hermitian interpolation functions constitutes the first portion of this section. After that is the exposition of the proposed mass lumping scheme.

4.1 The cover formed by a finite element mesh

The finite element mesh we are using can also be regarded as a cover. The cover is made up of a series of patches. Each patch corresponds to a node and is the union of all the elements connected to the node. As a result, the patch can be given an index identical to that of the node, and occupies a domain represented by {\tilde{\mathrm{\Omega }}}_i if the node index is i. Accordingly the number of the patches in the cover equals that of the nodes in the mesh, denoted by n. All the n patches overlap partially and constitute the cover of the plate domain \mathrm{\Omega }, written as \left\{{\tilde{\mathrm{\Omega }}}_i\right\}, i.e.,

(55)\mathrm{\Omega }=\bigcup _{i=1}^n{\tilde{\mathrm{\Omega }}}_i.

For ease of presentation, it might as well be supposed that all the elements in the mesh are 4-node quadrilateral elements. Then, any element in the mesh is covered by the 4 patches associated with the 4 nodes of the element.

The following is an illustrative example. Shown in Fig.4(a) is a plate discretized with a mesh of 27 rectangular elements and 40 nodes. So there are 40 patches associated with the 40 nodes. Let us look at the point P in element 13, which is covered by the 4 patches indexed by 23{18-19-13-12}, 24{19-20-14-13}, 17{13-14} and 16{12-13-9}. Here, the numbers in front of the curly braces are the patch indices, the numbers in the curly braces are the indices of the elements constituting the patches. For example, 16{12-13-9} means that patch 16 is made up of elements 12, 13, and 9. Fig.4(b) displays separately the four patches.

Associated with the ith patch {\tilde{\mathrm{\Omega }}}_i is a weight function denoted by p_i\left(x,y\right). All the weights \left\{p_i\right\} subordinate to the cover \left\{{\tilde{\mathrm{\Omega }}}_i\right\} satisfy the following properties

(56)\{\begin{array}{ll}{p}_i\left(x,y\right)=0,& \mathrm{\forall }\left(x,y\right)\notin {\tilde{\mathrm{\Omega }}}_i,\\ {\displaystyle \sum _{i=1}^n{p}_i\left(x,y\right)=1,}& \mathrm{\forall }\left(x,y\right)\in \mathrm{\Omega }.\end{array}

Due to the first property p_i is called to have support {\tilde{\mathrm{\Omega }}}_i; and the second property renders us to refer to \left\{p_i\right\} as the PU functions subordinate to the cover \left\{{\tilde{\mathrm{\Omega }}}_i\right\}.

As introduced in Subsection 2.2, the zeroth order shape functions of all the elements satisfy the PU properties, which can be collected as a simple and direct way to formulate \left\{p_i\right\}, and together with all elements constituting patch {\tilde{\mathrm{\Omega }}}_i can be defined,

(57)p_i\left(x,y\right)=\{\begin{array}{ll}{N}_{j,e}^0\left(x,y\right),& \mathrm{i}\mathrm{f}e\left[j\right]=i\mathrm{a}\mathrm{n}\mathrm{d}\left(x,y\right)\in {\mathrm{\Omega }}^e,\\ 0,& \mathrm{i}\mathrm{f}\left(x,y\right)\notin {\tilde{\mathrm{\Omega }}}_i,\end{array}

where {\mathrm{\Omega }}^e is one of the element domains constituting patch {\tilde{\mathrm{\Omega }}}_i; N_{j,e}^0 is the zeroth order shape function of the jth node of element e; and e\left[j\right] denotes the global index of the jth node of element e. Taking patch 16 shown in Fig.4(b) as an instance, we have

(58)p_{16}\left(x,y\right)=\{\begin{array}{ll}{N}_{2,9}^0\left(x,y\right),& \mathrm{i}\mathrm{f}\left(x,y\right)\in {\mathrm{\Omega }}^9,\\ N_{3,12}^0\left(x,y\right),& \mathrm{i}\mathrm{f}\left(x,y\right)\in {\mathrm{\Omega }}^{12},\\ N_{4,13}^0\left(x,y\right),& \mathrm{i}\mathrm{f}\left(x,y\right)\in {\mathrm{\Omega }}^{13},\\ 0,& \mathrm{i}\mathrm{f}\left(x,y\right)\notin {\tilde{\mathrm{\Omega }}}^{16},\end{array}

because the 3rd node of element 12 has the global index 16, and so on. Here {\tilde{\mathrm{\Omega }}}_{16}={\mathrm{\Omega }}^9\bigcup {\mathrm{\Omega }}^{12}\bigcup {\mathrm{\Omega }}^{13}.

If the weight functions p_i are defined in this way, we need to find out the proper approximations of the plate deflection over patches {\tilde{\mathrm{\Omega }}}_i. This needs to view from the PU perspective the Hermitian interpolation to the plate deflection w defined in Eq. (16), so as to retrieve the local approximations over patch {\tilde{\mathrm{\Omega }}}_i corresponding to the PU function p_i.

It is observed that the deflection w\left(x,y;t\right) at point \left(x,y\right)\in {\mathrm{\Omega }}^e can be approximated by the summation of contributions from the four patches {{\tilde{\mathrm{\Omega }}}_i} that covers element e. For patch {\tilde{\mathrm{\Omega }}}_i covering {\mathrm{\Omega }}^e, the global index of kth node for element e is assumed to be i. Thus, we have e\left[k\right]=i, then the subordinate function of patch {\tilde{\mathrm{\Omega }}}_i to w\left(x,y;t\right) in {\mathrm{\Omega }}^e can be written as

(59)\begin{array}{rl}{\mathit{N}}_{k,e}{\mathit{d}}_{k,e}\equiv & N_{k,e}^0{d}_{k,e}^0+\cdots +N_{k,e}^{\mathrm{D}}{d}_{k,e}^{\mathrm{D}}\\ \equiv & N_{k,e}^0\left(d_{k,e}^0+b_{k,e}^1{d}_{k,e}^1+\cdots +b_{k,e}^{\mathrm{D}}{d}_{k,e}^{\mathrm{D}}\right)\\ \equiv & N_{k,e}^0\left(d_i^0+d_i^1{b}_{k,e}^1+\cdots +d_i^{\mathrm{D}}{b}_{k,e}^{\mathrm{D}}\right),\end{array}

recalling that a node linked by all surrounding elements shares a unique DOFs denoted by d_i^j=d_{k,e}^j, j=0,…,D. Here b_{k,e}^j are defined from higher-order shape function of Hermitian interpolation:

(60)b_{k,e}^j\left(x,y\right)=\frac{N_{k,e}^j\left(x,y\right)}{N_{k,e}^0\left(x,y\right)},j=1,\cdots ,D.

Because the zeroth shape functions N_{k,e}^0\left(x,y\right), k=1,\cdots ,4, have the PU properties, we can conclude that the D_1 functions in the parentheses of Eq. (59): 1(one), b_{k,e}^1,…,b_{k,e}^{\mathrm{D}}, constitutes the local approximations basis over the subset {\mathrm{\Omega }}^e of patch {\tilde{\mathrm{\Omega }}}_i. Finally, each higher-order basis function b_i^1,\dots ,b_i^{\mathrm{D}} over the whole patch {\tilde{\mathrm{\Omega }}}_i can be calculated element-by-element according to

(61)b_i^j\left(x,y\right)=\{\begin{array}{ll}{b}_{k,e}^j\left(x,y\right),& \mathrm{i}\mathrm{f}(x,y)\in {\mathrm{\Omega }}^e,\\ 0,& \mathrm{i}\mathrm{f}\left(x,y\right)\notin {\tilde{\mathrm{\Omega }}}_i,\end{array}j=1,\cdots ,D.

That is to say, any local approximation over {\tilde{\mathrm{\Omega }}}_i can be written as

(62)w_i\left(x,y\right)=\sum _{j=0}^D{d_i^{j}b}_i^j\left(x,y\right)={\mathit{b}}_i{\mathit{d}}_i,

in which

(63){\mathit{d}}_i^{\mathrm{T}}=\left(\begin{array}{ccc}{d}_i^0& \cdots & d_i^{\mathrm{D}}\end{array}\right),

is the vector of D_1 DOFs attached to patch {\tilde{\mathrm{\Omega }}}_i or node-i; and

(64){\mathit{b}}_i=\left(\begin{array}{ccc}{b}_i^0& \cdots & b_i^{\mathrm{D}}\end{array}\right),

is the row vector of D_1 basis functions on patch {\tilde{\mathrm{\Omega }}}_i, defined element-by-element in Eq. (61), with b_i^0=1.

However, one major difference from implementation with other PU method [42] or the NMM [43], the local approximation is computed element-by-element constituting {\tilde{\mathrm{\Omega }}}_i, rather than over the entire patch.

4.2 The manifold-based mass lumping scheme

Let us return to the virtual work of inertia force and rewrite it here

(65)\left(\delta w,\ddot{w}\right)={\int }_{\mathrm{\Omega }}^{}\rho \delta w\ddot{w}h\mathrm{d}\mathrm{\Omega },

for convenience. Further, let the plate domain \mathrm{\Omega } be covered by \left\{{\tilde{\mathrm{\Omega }}}_i\right\} and \left\{p_i\right\} be the weight functions subordinate to \left\{{\tilde{\mathrm{\Omega }}}_i\right\}. By multiplying the unity \sum p_i with the integrand of \left(\delta w,\ddot{w}\right), it reduces to

(66)\begin{array}{rl}\left(\delta w,\ddot{w}\right)=& {\int }_{\mathrm{\Omega }}^{}\rho h\delta w\ddot{w}\sum _{i=1}^n{p}_i\mathrm{d}\mathrm{\Omega }=\sum _{i=1}^n{\int }_{\mathrm{\Omega }}^{}\rho h{p}_i\delta w\ddot{w}\mathrm{d}\mathrm{\Omega }\\ =& \sum _{i=1}^n{\int }_{{\tilde{\mathrm{\Omega }}}_i}^{}\rho h{p}_i\delta w\ddot{w}\mathrm{d}\mathrm{\Omega }=\sum _{i=1}^n\delta W_{\mathrm{I}}^i,\end{array}

because p_i=0 outside {\tilde{\mathrm{\Omega }}}_i; where

(67)\delta W_{\mathrm{I}}^i\equiv {\int }_{{\tilde{\mathrm{\Omega }}}_i}^{}\rho h{p}_i\delta w\ddot{w}\mathrm{d}\mathrm{\Omega }.

The derivation of Eq. (66) is based on the definition of integral of a scalar function over a manifold. The scalar function is \rho h\delta w\ddot{w}, and the manifold is the plate domain \mathrm{\Omega }, a trivial 2-dimensional manifold.

Prior to approximating \delta W_{\mathrm{I}}^i, it is observed that the integrand of \left(\delta w,\ddot{w}\right) includes merely the test function \delta w and the trial function w, with no spatial derivatives of them involved. This is different from a\left(\delta w,w\right) which includes the spatial second order derivatives of both \delta w and w (see Eq. (14)). Hence, the numerical computation of \left(\delta w,\ddot{w}\right) requires neither \delta w nor w to be as smooth as they are in the computation of a\left(\delta w,w\right). In other words, it is adequate to require \delta w,w\in L^2\left(\mathrm{\Omega }\right) in approximating \left(\delta w,\ddot{w}\right). Instead, the computation of a\left(\delta w,w\right) requires \delta w,w\in H^2\left(\mathrm{\Omega }\right), if the standard Garlerkin framework is followed. This case is somewhat similar to the analysis of soil consolidation problems, where what arises in the virtual work done by pore pressure p is p itself with no derivatives of p involved, and the approximation to p in FEA is required to be at least one order less smoother than displacement [15].

On the basis of the above discussions, if patch {\tilde{\mathrm{\Omega }}}_i is small enough, \delta w and \ddot{w} can simply be substituted by the local approximations \delta w_i and {\ddot{w}}_i over {\tilde{\mathrm{\Omega }}}_i, respectively, leading to

(68)\delta W_{\mathrm{I}}^i\left(t\right)\approx \delta {\tilde{W}}_{\mathrm{I}}^i\left(t\right)\equiv {\int }_{{\tilde{\mathrm{\Omega }}}_i}^{}\rho h{p}_i\delta w_i{\ddot{w}}_i\mathrm{d}\mathrm{\Omega }.

According to Eq. (62),

(69)\delta w_i(x,y)={\mathit{b}}_i(x,y)\delta ,

where \delta {\mathit{d}}_i is a D_1-dimensional vector of virtual DOFs, \delta d_i^0,…,\delta d_i^D; and

(70){\ddot{w}}_i\left(x,y;t\right)={\mathit{b}}_i\left(x,y\right){\ddot{\mathit{d}}}_i\left(t\right).

Substituting Eqs. (69) and (70) into Eq. (68), we have

(71)\delta {\tilde{W}}_{\mathrm{I}}^i\left(t\right)=\delta {\mathit{d}}_i^{\mathrm{T}}{\mathit{M}}_{\mathrm{L}}^i{\ddot{\mathit{d}}}_i\left(t\right),

with {\mathit{M}}_{\mathrm{L}}^i a D_1 × D_1 matrix that is symmetric

(72){\mathit{M}}_{\mathrm{L}}^i={\int }_{{\tilde{\mathrm{\Omega }}}_i}^{}\rho h{p}_i{\mathit{b}}_i^{\mathrm{T}}{\mathit{b}}_i\mathrm{d}\mathrm{\Omega },

and because basis functions b_i^j,j=0,…,D are linearly independent, the integral

(73){\int }_{{\tilde{\mathrm{\Omega }}}_i}^{}\rho h{\mathit{b}}_i^{\mathrm{T}}{\mathit{b}}_i\mathrm{d}\mathrm{\Omega }>0.

Hence, if

(74){\int }_{{\tilde{\mathrm{\Omega }}}_i}^{}{p}_i\mathrm{d}\mathrm{\Omega }>0,

based on the second mean value theorem for integrals, {\mathit{M}}_{\mathrm{L}}^i given by Eq. (72) is also positive definite. And this is satisfied by the zeroth shape functions of both ACM and BFS elements.

We can then yield the approximated virtual work done by inertia force by combining Eqs. (68) and (71) with Eq. (66), writen as

(75)\left(\delta w,\ddot{w}\right)=\sum _{i=1}^n\delta {\mathit{d}}_i^{\mathrm{T}}{\mathit{M}}_{\mathrm{L}}^i{\ddot{\mathit{d}}}_i=\delta {\mathit{d}}^{\mathrm{T}}{\mathit{M}}_{\mathrm{L}}\ddot{\mathit{d}},

where \delta \mathit{d} is a vector consisting of the n sub-vectors \delta {\mathit{d}}_1,…,\delta {\mathit{d}}_{\mathrm{n}}; {\mathit{M}}_{\mathrm{L}} is the block-diagonal LMM defined as

(76){\mathit{M}}_{\mathrm{L}}=\left[\begin{array}{ccc}{\mathit{M}}_{\mathrm{L}}^1& \cdots & 0\\ ⋮& \cdots & ⋮\\ 0& \cdots & {\mathit{M}}_{\mathrm{L}}^n\end{array}\right],

with {\mathit{M}}_{\mathrm{L}}^i defined in Eq. (72). However, each {\mathit{M}}_{\mathrm{L}}^i is a symmetric positive definite matrix with a dimension of 3 × 3 for the ACM element and 4 × 4 for the BFS element, hence its Cholesky decomposition

(77){\mathit{M}}_{\mathrm{L}}^i={\mathit{L}}_i{\mathit{L}}_i^{\mathrm{T}},

can be calculated using very few calculations.

However, all the deductions just serve for theoretical purpose, in which the global LMM {\mathit{M}}_{\mathrm{L}} accordingly will be obtained patch by patch. Considering that each element will simultaneously contribute to the four patches covering it, in practice, {\mathit{M}}_{\mathrm{L}}^e of element e can be computed first,

(78){\mathit{M}}_{\mathrm{L}}^e=\left[\begin{array}{ccc}{\mathit{M}}_{\mathrm{L}}^{1,e}& \cdots & 0\\ ⋮& \cdots & ⋮\\ 0& \cdots & {\mathit{M}}_{\mathrm{L}}^{4,e}\end{array}\right],

and then add {\mathit{M}}_{\mathrm{L}}^{k,e}, to the mass matrix {\mathit{M}}_{\mathrm{L}}^i of patch i=e\left[k\right], for k=1,…,4. In this way, the global LMM {\mathit{M}}_{\mathrm{L}} is formed by assembling all the element mass matrices in the similar element by element form as global stiffness matrix,

(79){\mathit{M}}_{\mathrm{L}}=\sum _e{\mathit{M}}_{\mathrm{L}}^e,

which enhance the computation efficiency since at each quadrature points, the Jacobian only entails to be calculated once rather than 4 times through a patch by patch assembling manner.

In Eq. (78),

(80){\mathit{M}}_{\mathrm{L}}^{k,e}={\int }_{{\mathrm{\Omega }}^e}^{}\rho h{N}_{k,e}^0{\mathit{b}}_{k,e}^{\mathrm{T}}{\mathit{b}}_{k,e}\mathrm{d}\mathrm{\Omega },

represents the contribution of the kth node of element e to the global mass matrix {\mathit{M}}_{\mathrm{L}}, or the mass matrix of this node, recalling the D_1-dimensional row-vector

(81){\mathit{b}}_{k,e}(x,y)=\left[1,b_{k,e}^1(x,y),\cdots ,b_{k,e}^{\mathrm{D}}(x,y)\right],

is the restriction of the basis vector of local approximations of patch e\left[k\right] onto element e, and b_{k,e}^j is defined by Eq. (60).

To this point, the row-sum method should be believed to be a rigorous method because it can be regarded as a particular case of this proposed mass lumping scheme by taking the basis vector over {\tilde{\mathrm{\Omega }}}_i as {\mathit{b}}_i = {1}. As a consequence, dropping the basis functions corresponding to rotational/torsional DOFs must cause rotational/torsional DOFs to zero mass in the resulting DLMM.

5 Numerical examples

In this study, the ACM element and BFS element are taken for example [33,34] respectively. They are almost the earliest developed plate elements and the representative of incompatible and compatible elements. Using the proposed mass lumping scheme, nevertheless, it will be seen shortly that they behave at least as well as or better than elements developed later.

In the following numerical examples, the same geometric and mechanical parameters are assigned as: length of the square plate a = 10 m and thickness h = 0.05 m, Young’s modulus E = 200 GPa, and Poisson’s ratio ν = 0.3; and the density ρ = 8000 kg/m³.

5.1 Free vibration of plates

First we investigate the accuracy and efficiency of the proposed manifold-based mass lumping scheme for extracting the natural frequencies of thin plates, under various boundary conditions.

Example 1: thin square plate with four boundaries completely free.

This benchmark is recommended by NAFEMS (National Agency for Finite Element Methods and Standards) and tested in Abaqus. The mesh is formed by a uniform partition of elements 10 × 10.

To evaluate the effectiveness and accuracy of this manifold-based mass lumping method, the non-dimensional natural frequency coefficients \tilde{\omega }=({\omega }^2\rho h{a}^4/D_0{)}^{1/4} are computed and compared, where a denotes the length size and D_0=E{t}^3/\left[12\left(1-{\upsilon }^2\right)\right] the flexural rigidity.

The results yielded by the ACM and BFS elements together with the three types of mass matrices, namely the CMM, the row-sum LMM and the proposed LMM, are listed in in Tab.1, in comparison with the analytical solutions in the first column.

According to the results listed in Tab.1, the natural frequencies yielded by the manifold-based mass lumping scheme using both ACM and BFS are in good agreement with the analytical solution [44], as good as the CMM in the lower modes. For the higher modes, more accurate results are obtained than the CMM scheme. In all the cases, the proposed scheme is more accurate than the row-sum method.

The first six mode shapes yielded by both ACM and BFS with the manifold-based mass lumping scheme agrees well with Abaqus results (for element type S4R5) [45].

Example 2: thin square plate with four boundaries simply supported.

The plate is simply supported along the four boundaries. The mesh is formed by a uniform partition of 20 × 20 elements. In addition to the exact solution, the shear-locking-free triangular element (DSG3) [46] is also used. The eight columns in Tab.2 list the first six frequencies scaled as in example 1.

According to Tab.2, using the three mass lumping schemes, both ACM and BFS elements can give rise to much more accurate results than those by DSG3 using a mesh of 22 × 22 triangular elements with the CMM.

The lower six mode shapes evaluated by ACM and BFS elements integrated with proposed mass matrices are depicted in Fig.5 and Fig.6, which accord with those mode shapes with the consistent mass matrices (not displayed for shortening length).

Example 3: thin square plate with four boundaries fully clamped.

In this example, the plate is fully clamped. The mesh is formed by a uniform partition of 20 × 20 elements.

The first six scaled natural frequencies assessed by the ACM and BFS elements with the three types of mass matrices are listed in Tab.3, with comparisons with the analytical solutions and the results by shear-locking-free triangular finite element (DSG3) formulation [46].

The same conclusions as the above example are true of this example.

The first six mode shapes evaluated by both ACM and BFS with the proposed mass matrices agree with those by the same elements with the consistent mass matrices.

Example 4: Thin square plate with simply supported-free boundaries.

In this example, the plate is fully clamped. The mesh is formed by a uniform partition of 20 × 20 elements.

The first six scaled natural frequencies assessed by the ACM and BFS elements with the three types of mass matrices are listed in Tab.4, with comparisons with the analytical solutions and the results by shear-locking-free triangular finite element (DSG3) formulation [46].

According to Tab.4, using the three mass lumping schemes, ACM element can give rise to much more accurate results than those by DSG3 using a mesh of 22 × 22 triangular elements with the CMM. However, a discrepancy is observed between the results from the Consistent and Row-sum methods for the BFS element and the analytical solution. This maybe due to the fact that by ensuring compatibility, the strict continuity requirements in compatible BFS element can sometimes lead to overly stiff responses, potentially causing locking phenomena or reduced numerical accuracy, especially when capturing higher-order eigenmodes. The ACM on the other hand, its intentional incompatibility can help avoid excessive stiffness.

Example 5: thin square plate with simply clamped-free.

In this example, the plate is fully clamped. The mesh is formed by a uniform partition of 20 × 20 elements.

The first six scaled natural frequencies assessed by the ACM and BFS elements with the three types of mass matrices are listed in Tab.5, with comparisons with the analytical solutions and the results by shear-locking-free triangular finite element (DSG3) formulation [46].

The same conclusions as the Example 4 are true of this example.

Moreover, it should be noted that both the ACM and BFS elements exhibit lower accuracy in predicting higher-order modes in this example.

Example 6: thin square plate with clamped-simply supported.

In this example, the plate is fully clamped. The mesh is formed by a uniform partition of 20 × 20 elements.

The first six scaled natural frequencies assessed by the ACM and BFS elements with the three types of mass matrices are listed in Tab.6, with comparisons with the analytical solutions.

The deviation in results for the Consistent and Row-sum methods of the BFS plate element is also observed in this case.

Example 7: cantilevered thin square plate.

This is also a standard NAFEMS benchmark investigated in Abaqus benchmarks manual [45]. In this example, the plate is clamped along the y-axis. The mesh is formed by an even partition of 20 × 20 elements. The first six scaled natural frequencies using the ACM and BFS elements with three types of mass matrices are listed in Tab.7, in comparison with the analytical solutions [44].

For this example, we have the same conclusions as the above examples.

5.2 Forced vibration of plates

The standard benchmarks concerning the forced vibration are investigated in this section to test the accuracy of the new mass lumping scheme. These tests are recommended by the NAFEMS. Harmonic and transient forcing functions are applied on the plate with the simply-supported boundary conditions. The results are compared with consistent and lumped mass formulations, and the Abaqus schemes and analytical solutions [44].

Example 8: harmonic response analysis of simply supported square plate.

The forcing function is described as follows. The uniformly distributed steady-state harmonic function: \overline{p}(x,y;t)={\overline{p}}_0\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega t\right), {\overline{p}}_0=100\mathrm{N}/{\mathrm{m}}^2 is exerted on the whole plate, with \omega =2\pi f. The frequency f varies from 0 to 4.16 Hz so as to capture the damped resonant frequency of the plate in considering the effect of damping, denoted by {\omega }_{\mathrm{R}}. Usually, the effect of damping makes {\omega }_{\mathrm{R}} slightly less than the fundamental natural frequency {\omega }_1 [41]. The damping factors in Eq. (30) for constructing the damping matrix \mathit{C} are {\alpha }_{\mathrm{R}} = 0.299 s^–1 and {\beta }_{\mathrm{R}}=1.339 × 10^–3 s. The central w and {\sigma }_x responses are calculated.

As for this benchmark, the results assessed by the three kinds of mass matrices using ACM and BFS elements are compared with the reference solution [47]. And the results evaluated by Abaqus using the S4R, S4R5, S8R, and S8R5 elements are also listed as comparisons in Tab.8.

According to Tab.8, elements developed much later than the ACM and BFS elements do not necessarily generate results better than the latter. For example, the peak stress calculated by S8R5 is 34.49, with a bigger error than the BFS element using the manifold-based LMMs. Once again using both ACM and BFS together with proposed mass lumping technique, it is verified from Tab.8 both ACM and BFS reproduce “{\omega }_{\mathrm{R}} is slightly less than {\omega }_1.” However, none of the later developed elements produces this result.

The deflection response curves in the frequency domain that are predicted by the ACM, BFS elements (using the proposed lumping scheme), S8R5 elements and the reference solution almost coincide. Therefore, only the deflection response curve by BFS elements in the frequency domain is demonstrated in Fig.7.

Example 9: vibration of thin plate under the transient load.

Assuming a suddenly applied pressure p_0=100\mathrm{N}/{\mathrm{m}}^2 is exerted on the simpl-supported square plate middle plane and kept. The damping factors are selected as {\alpha }_{\mathrm{R}} = 0.299 s^–1 and {\beta }_{\mathrm{R}}=1.339 × 10^–3 s.

In implementation, the evolution of w and {\sigma }_x at center of the plate are calculated and recorded. The Newmark method is employed by setting the time step \mathrm{\Delta }t= 0.002 s and corresponding parameters \gamma =\beta = 0.5 to ensure the implicit time scheme is unconditional stable.

Judging from the Tab.9, the results assessed by the ACM and BFS elements using the consistent and lumped mass matrices are both in excellent agreement with the reference solution within the 250 time intervals. The peak response occurs at the time t_{\mathrm{p}} = 0.2100 s.

Depicted in Fig.8 is the response curve of deflection predicted by the ACM elements using the manifold-based LMM. The results computed by ACM + Consistent mass, ACM + Proposed lumping scheme almost completely coincide with the reference response predicted by Abaqus [45].

Example 10: an undamped thin plate under forced vibration.

A sinuous time varying dynamic pressure F=F_0\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega t\right) is uniformly exerted to the middle plane of undamped, simply-supported thin plate, with F_0=100\mathrm{N}/{\mathrm{m}}^2 and \omega =4.

For evaluation of mass lumping techniques, the central deflection w, velocity \dot{w}, and acceleration \ddot{w} are specifically studied and calculated. The same parameters for Newmark method as last numerical examples are selected as time step \mathrm{\Delta }t= 0.002 s and \gamma =\beta =0.5 which is unconditional stable [48]. The analytic solution for this example is available as follows:

(82){\omega }_{ij}={\pi }^2\left(\frac{i^2}{a^2}+\frac{j^2}{b^2}\right)\sqrt{\frac{D}{\rho h}},

(83)\begin{array}{rl}w\left(x,y,t\right)=& \sum _{i=1,3,5,\dots }\sum _{j=1,3,5,\dots }\frac{16F_0}{{\pi }^2{D}_{0}ij\left({\pi }^4{\left({\displaystyle \frac{i^2}{a^2}}+{\displaystyle \frac{j^2}{b^2}}\right)}^2-\left({\displaystyle \frac{\rho h}{D}}\right){\mathrm{\Omega }}^2\right)}\\ & \cdot \mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{i}\pi x}{a}\mathrm{s}\mathrm{i}\mathrm{n}\frac{\mathrm{i}\pi y}{b}\left(\mathrm{s}\mathrm{i}\mathrm{n}\omega t-\frac{\omega }{{\omega }_{ij}}\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{ij}t\right)\right).\end{array}

According to Tab.10, the results yielded by ACM and BFS elements using the consistent and lumped mass matrices are both in excellent agreement with the analytical solution. BFS elements with both consistent mass and the proposed mass lumping formulations excel in the accuracy of peak deflection than the ACM elements. The evolution of central deflection, velocity and acceleration are demonstrated by selecting 1000-time intervals. The peak-time for BFS elements using manifold-based mass matrix matches with analytical occurring at 1.9920 s. In Tab.11, the velocity and acceleration at time t = 2 s is calculated. The numerical results demonstrate that the velocity and acceleration at t = 2 s gained by the proposed method agree well with those gained by the consistent mass formulations.

In Fig.9, Fig.10, and Fig.11, the evolution of central deflection, velocity and acceleration are demonstrated using the BFS elements with the manifold-based mass lumping scheme. The deflection and velocity results assessed by BFS + Consistent mass, BFS + Proposed lumping scheme are almost completely coinciding with the analytical results. The acceleration response curve also agrees well with the analytical acceleration response curve, albeit bigger error than the deflection and velocity results.

Now we report the central deflection errors in the period from 0 to 2 s in Fig.12. From it, the error given by BFS + Consistent mass is indeed the smallest among the four combinations. Nevertheless, BFS + Proposed lumping scheme behaves slightly inferior to BFS + Consistent, using computational costs and memory much less than the consistent mass formulation. For the central velocity errors and accelerations errors shown in Fig.13 and Fig.14, similar conclusions can be drawn. Considering lumped mass matrices are mainly used in scenarios where high speed and strong nonlinearity arise, the proposed lumping scheme should play its due role.

Finally, we studied the convergence of the proposed LMMs employing the Newmark method by comparing peak deflection and the occurrence time with different mesh density. Tab.12 lists the results corresponding to different mesh densities and evaluated by BFS + Consistent, BFS + Proposed and the analytic solution, which manifest that both the mass schemes converge very well with increasing mesh grids, and the proposed manifold-based LMM appears slightly better than the CMM.

6 Discussion and conclusions

In the present study, the challenge of accurately representing rotational/torsional DOFs in Kirchhoff plate dynamic analysis has been addressed. Unlike existing methodologies, this paper introduces a rigorous, manifold-based mass lumping scheme that yields strictly symmetric positive-definite diagonal mass matrices, effectively eliminating rank deficiencies and associated spurious zero-frequency modes.

It is shown that the commonly used row-sum method turns to be rigorous in theory because it is a particular case of the manifold-based mass lumping scheme by specifying constants as the local approximations for all the patches [49].

Typical benchmark examples are analyzed, including free vibration, harmonic response analysis, forced vibration subjected to transient and periodic forcing functions. For all the three problem types, the ACM and BFS plate elements using the manifold-based mass lumping scheme behave as well as or better than elements using the consistent mass or row-sum mass lumping scheme, and much better than those elements developed later.

As a consequence, it can be stated that the manifold-based mass lumping scheme is able to supersede not only the most commonly applied row-sum mass lumping scheme but also the consistent mass scheme.

References

[1] PetytM. Introduction to Finite Element Vibration Analysis. Cambridge: Cambridge University Press, 2010

[2] Zheng H , Liu Z , Ge X . Numerical manifold space of Hermitian form and application to Kirchhoff’s thin plate problems. International Journal for Numerical Methods in Engineering, 2013, 95(9): 721–739

[3] Guminiak M . Free vibrations analysis of thin plates by the boundary element method in non-singular approach. Scientific research of the Institute of Mathematics and Computer Science, 2007, 6: 75–90

[4] Useche J , Harnish C . A boundary element method formulation for modal analysis of doubly curved thick shallow shells. Applied Mathematical Modelling, 2016, 40(5–6): 3591–3600

[5] LiSLiuW K. Meshfree Particle Methods. Heidelberg: Springer Science and Business Media, 2007

[6] LiuG R. Meshfree Methods: Moving Beyond the Finite Element Method. Boca Raton: Taylor and Francis, 2009

[7] Hayri M N , Hakan E , Bekir A , Ömer C . A new eigenvalue problem solver for thermo-mechanical vibration of Timoshenko nanobeams by an innovative nonlocal finite element method. Mathematical Methods in the Applied Sciences, 2022, 45(5): 2592–2614

[8] Felippa C A . A historical outline of matrix structural analysis: A play in three acts. Computers and Structures, 2001, 79(14): 1313–1324

[9] FelippaC A. Recent advances in finite element templates. In: Computational Mechanics for the Twenty-First Century. Stirling: Civil-Comp Press, 2000, 71–98

[10] Felippa C A , Guo Q , Park K C . Mass matrix templates: General description and 1D examples. Archives of Computational Methods in Engineering, 2015, 22(1): 1–65

[11] ZienkiewiczO CTaylorR L. The Finite Element Method for Solid and Structural Mechanics. Oxford: Butterworth-Heinemann, 2005

[12] Felippa C A . Construction of customized mass-stiffness pairs using templates. Journal of Aerospace Engineering, 2006, 19(4): 241–258

[13] Archer J S . Consistent mass matrix for distributed mass systems. Journal of the Structural Division, 1963, 89(4): 161–178

[14] Archer J S . Consistent matrix formulations for structural analysis using finite-element techniques. AIAA Journal, 1965, 3(10): 1910–1918

[15] HughesT J R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. North Chelmsford: Courier Corporation, 2012

[16] Hinton E S , Rock T J , Zienkiewicz O C . A note on mass lumping and related processes in the finite element method. Earthquake Engineering and Structural Dynamics, 1976, 4(3): 245–249

[17] Duczek S , Gravenkamp H . Critical assessment of different mass lumping schemes for higher order serendipity finite elements. Computational Methods in Applied Mathematics, 2019, 350: 836–897

[18] Gresho P M , Lee R L , Sani R L . Advection-dominated flows, with emphasis on the consequences of mass lumping. Finite Elements in Fluids, 1978, 3: 335–350

[19] BatheK J. Finite Element Procedures. State of New Jersey: Printice Hall, 2006

[20] VoetYSandeEBuffaA. Mass lumping and outlier removal strategies for complex geometries in isogeometric analysis. Mathematics of Computation, 2025 (in press)

[21] Archer G C , Whalen T M . Development of rotationally consistent diagonal mass matrices for plate and beam elements. Computational Methods in Applied Mathematics, 2005, 194(6–8): 675–689

[22] Zou Z , Scott M A , Miao D , Bischoff M , Oesterle B , Dornisch W . An isogeometric Reissner–Mindlin shell element based on Bézier dual basis functions: Overcoming locking and improved coarse mesh accuracy. Computer Methods in Applied Mechanics and Engineering, 2020, 370: 113283

[23] Voet Y , Sande E , Buffa A . A mathematical theory for mass lumping and its generalization with applications to isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2023, 410: 116033

[24] Geevers S , Maier R . Fast mass lumped multiscale wave propagation modelling. IMA Journal of Numerical Analysis, 2023, 43(1): 44–72

[25] Hiemstra R R , Nguyen T H , Eisenträger S , Dornisch W , Schillinger D . Higher-order accurate mass lumping for explicit isogeometric methods based on approximate dual basis functions. Computational Mechanics, 2025, 76(1): 205–226

[26] Li X , Hou S , Wang D . An interpolatory basis lumped mass isogeometric formulation with rigorous assessment of frequency accuracy for Kirchhoff plates. Thin-walled Structures, 2024, 197: 111639

[27] Tan P , Nguyen-Thanh N , Rabczuk T , Zhou K . Static, dynamic and buckling analyses of 3D FGM plates and shells via an isogeometric-meshfree coupling approach. Composite Structures, 2018, 198: 35–50

[28] Huang J , Nguyen-Thanh N , Gao J , Fan Z , Zhou K . Static, free vibration, and buckling analyses of laminated composite plates via an isogeometric meshfree collocation approach. Composite Structures, 2022, 285: 115011

[29] Nguyen-Thanh N , Li W , Zhou K . Static and free-vibration analyses of cracks in thin-shell structures based on an isogeometric-meshfree coupling approach. Computational Mechanics, 2018, 62(6): 1287–1309

[30] Vu-Bac N , Duong T X , Lahmer T , Zhuang X , Sauer R A , Park H S , Rabczuk T . A NURBS-based inverse analysis for reconstruction of nonlinear deformations of thin shell structures. Computer Methods in Applied Mechanics and Engineering, 2018, 331: 427–455

[31] Vu-Bac N , Duong T X , Lahmer T , Areias P , Sauer R A , Park H S , Rabczuk T . A NURBS-based inverse analysis of thermal expansion induced morphing of thin shells. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 480–510

[32] Vu-Bac N , Rabczuk T , Park H S , Fu X , Zhuang X . A NURBS-based inverse analysis of swelling induced morphing of thin stimuli-responsive polymer gels. Computer Methods in Applied Mechanics and Engineering, 2022, 397: 115049

[33] Zheng H , Yang Y . On generation of lumped mass matrices in partition of unity based methods. International Journal for Numerical Methods in Engineering, 2017, 112(8): 1040–1069

[34] Yang Y , Zheng H , Sivaselvan M V . A rigorous and unified mass lumping scheme for higher-order elements. Computational Methods in Applied Mathematics, 2017, 319: 491–514

[35] Li W , Lin S , Wang Z , Guo H , Yu X . An explicit improved meshless numerical manifold method for dynamic crack propagation. Theoretical and Applied Fracture Mechanics, 2024, 130: 104293

[36] BäthgeFProvatidisC GJuhreDEisenträgerS. Novel Mass Lumping Approach Leveraging the Spectral Decomposition Theorem. 2025

[37] Lin S , Cao X , Zheng H , Li Y , Li W . An improved meshless numerical manifold method for simulating complex boundary seepage problems. Computers and Geotechnics, 2023, 155: 105211

[38] Ventsel E , Krauthammer T , Carrera E . Thin plates and shells: theory, analysis, and applications. Applied Mechanics Reviews, 2002, 55(4): B72–B73

[39] AdiniA. Analysis of shell strutures by the finite element method. Dissertation for the Doctoral Degree. Berkeley: University of California, 1961

[40] Schmit L A , Bogner F K , Fox R L . Finite deflection structural analysis using plate and shell discreteelements. AIAA Journal, 1968, 6(5): 781–791

[41] HitchingsD. A Finite Element Dynamics Primer. East Kilbride: National Agency for Finite Element Standards, 1992

[42] Babuška I , Melenk J M . The partition of unity method. International Journal for Numerical Methods in Engineering, 1997, 40(4): 727–758

[43] Zheng H , Liu F , Li C . The MLS-based numerical manifold method with applications to crack analysis. International Journal of Fracture, 2014, 190(1–2): 147–166

[44] AbassianFHawswellD JKnowlesN C. Free Vibration Benchmarks. East Kilbride: Department of Trade and Industry, National Engineering Laboratory, 1987

[45] DassaultSystemès. ABAQUS Version 6.14 User Documentation. Providence, RI: Dassault Systemes, 2014

[46] LiuG RNguyenT T. Smoothed Finite Element Methods. Boca Raton: CRC press, 2010

[47] MaguireJDawswellD JGouldL. Selected Benchmarks for Forced Vibration. East Kilbride: NAFEMS, 1989

[48] Newmark N M . A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, 1959, 85(3): 67–94

[49] Zhao S , Lin S , Dong M , Guo H , Zheng H . Mass lumping schemes fitted to MLS-based numerical manifold method in vibration of plates with cutouts using CPT and FSDT. Composite Structures, 2024, 330: 117815

RIGHTS & PERMISSIONS

Higher Education Press