A Moment-Fitted Extended Spectral Cell Method for Structural Health Monitoring Applications

Abstract

The spectral cell method has been shown as an efficient tool for performing dynamic analyses over complex domains. Its good performance can be attributed to the combination of the spectral element method with mesh-independent geometrical descriptions and the adoption of customized mass lumping procedures for elements intersected by a boundary, which enable it to exploit highly efficient, explicit solvers. In this contribution, we introduce the use of partition-of-unity enrichment functions, so that additional domain features, such as cracks or material interfaces, can be seamlessly added to the modeling process. By virtue of the optimal lumping paradigm, explicit time integration algorithms can be readily applied to the non-enriched portion of a domain, which allows one to maintain fast computing simulations. However, the handling of enriched elements remains an open issue, particularly with respect to stability and accuracy concerns. In addressing this, we propose a novel mass lumping method for enriched spectral elements in the form of a customized moment-fitting procedure and study its accuracy and stability. While the moment-fitting equations are deployed in an effort to minimize the lumping error, stability issues are alleviated by deploying a leap-frog algorithm for the solution of the equations of motion. This approach is numerically benchmarked in the 2D and 3D modeling of damaged aluminium components and validated in comparison with experimental scanning laser Doppler vibrometer data of a composite panel under piezo-electric excitation.

1. Introduction

The dynamic analysis of damaged domains is essential in several engineering applications, such as the study of impacts [1,2], seismic and hydraulic engineering [3,4], and structural health monitoring with guided waves [5,6]. Despite the progress accomplished in the last few decades, these technologies still involve a substantial computational cost, which limits their reach [5]. Typically, the necessity of representing high-frequency eigenmodes of a structure calls for high-resolution models, i.e., a fine discretization in both time and space. The discontinuous and/or singular character of features such as cracks or voids exacerbates these issues by complicating the modeling process, impacting the performance of traditional finite element (FE) models and severely affecting the conditional stability of explicit solvers [7,8,9]. In tackling these challenges, the engineer must adopt strategies pertaining to three main aspects of a model, which, as expressed by the underlying physics of wave propagation, are inevitably interlinked. These are: (i) a solution space appropriate for the dynamic phenomena of interest (space discretization); (ii) an efficient and versatile method to represent localized features such as defects; and (iii) an efficient and stable solution to the equations of motion (time discretization/integration).

For time integration, explicit solvers [10,11,12] can often outperform unconditionally stable, implicit [13,14,15] ones, due to their higher efficiency in terms of both memory storage usage and the number of operations. They require, however, the enforcement of stability conditions [16], as well as a diagonal mass matrix [17]. Although stability conditions in the form of a critical time step are coherent with the need for a fine time discretization, severe restrictions leading to excessive computational costs can arise from the introduction of very small elements in the meshing of details, as well as from discontinuities related to damage or interfaces. Since these features are often localized, possible remedies involve the adoption of leap-frog [18,19] or implicit–explicit (IMEX) [20,21,22] solvers. These strategies allow one to restrict the use of smaller time steps and implicit algorithms, respectively, to selected regions of the mesh, thus maintaining good performance for the majority of the domain, where explicit solvers with larger time steps can still be used. Secondly, mass lumping is a key aspect of explicit algorithms, as analysts are often willing to dispense with the variationally consistent formulation of the mass matrix in favor of accelerated solutions. With the exception of optimal lumping [23,24,25], it should be kept in mind that these methods cannot guarantee the convergence properties associated with the variational formulation and must therefore be intensively scrutinized.

Concerning space discretization, efficient methods for problems with smooth solutions can be achieved via high-order shape functions, as in the p-version of the FEM (p-FEM) [26,27,28], isogeometric analysis (IGA) [29,30,31], or the time-domain spectral element method (SEM) [25,32,33,34]. In the context of elastodynamics, these methods demonstrate similar capabilities in terms of evaluating the weak form in a variationally consistent sense [35]; however, only the SEM in a Gauss–Lobatto–Legendre (GLL) nodal configuration [36,37,38] offers a diagonal mass matrix within this formulation: an instance of optimal lumping, which is achieved via the nodal quadrature method [17,23,24]. This enables the deployment of explicit solvers without compromising accuracy, thus resulting in substantial performance gains. For these reasons, the SEM in conjunction with an explicit solver has seen widespread adoption in dynamic analysis [25,34,37,39,40].

Despite the effectiveness of high-order methods, when modeling damage, holes, or material interfaces, one is faced with the limitations of meshing software and the unsuitability of polynomial shape functions in approximating phenomena of a singular and/or discontinuous nature. Extended or generalized FE methods (XFEMs/GFEMs) [41,42,43,44] have enabled to overcome these issues by deploying mesh-independent geometrical descriptions [45] of these features, which are then used to locally supplement the solution space with partition-of-unity (PU) enrichment functions [46,47,48] and to generate appropriate quadrature rules for the evaluation of the weak form (e.g., [49,50]). In the presence of voids, the spectral cell method (SCM) [51,52,53,54,55,56] combines this kind of integration with a mass lumping procedure and a GLL spectral element (SE) mesh in order to perform dynamic analysis explicitly. On the other hand, the adoption of an XFEM in this context has been mainly driven by the study of dynamic crack propagation, which at first involved the use of IMEX [57] or fully implicit [58] solvers. This was due to the absence of mass lumping strategies for enriched elements and the near-zero mass coefficients [59] of nodes in the proximity of discontinuities. To overcome this, Menouillard et al. developed customized mass lumping routines for Heaviside enriched elements in 2D [7] and 3D [8] and demonstrated their excellent stability properties. This formulation, which assumed the conservation of discrete kinetic energy for rigid-body modes, was then generalized to arbitrary enrichment functions [60] and combined with the Chebyshev version of the SEM [9]. Asareh et al. [22] explored the possibility of enforcing energy conservation using only the standard mass coefficients; however, this necessitated the adoption of a node-based IMEX solver [61,62] to handle the resulting zero-valued enriched masses. It should be noted that these methods generally lead to a loss of convergence, as they introduce a customized expression of the mass matrix. In order to preserve optimal convergence rates, Sanchez-Rivandeira and Duarte [63] and Geelen et al. [64] applied Schweitzer’s variationally consistent lumping method [65], which could deliver block-diagonal mass matrices for non-negative shape functions.

In this contribution, we aimed to provide an efficient high-order method for explicit dynamic analysis in the presence of damage. Our main strategy consisted of extending the SCM with local PU enrichment functions in an attempt to combine the qualities of both approaches. Our discretization was similar to that of [9], with the difference that a GLL grid was employed in order to achieve the optimal lumping of standard SEs. Variationally consistent lumping via [65] of enriched SEs cannot be performed, since SE shape functions can assume negative values. Therefore, we developed a novel mass lumping method for enriched elements, which can be summarized as a nodal quadrature with customized weights. While integration at the element nodes delivered a (block-)diagonal mass matrix by construction, the weights were tailored to the enrichment functions and cut configuration of an element via moment-fitting equations. This setup is an established and versatile framework to generate quadrature rules. It can be used to integrate general polygons [66] and handle singular and discontinuous integrands [67,68,69,70]. To decrease the cost of the procedure, node locations could be chosen a priori [67,69,71,72]; however, one can also optimize them for accuracy [66,68,73]. In fictitious domain methods, moment fitting enables a drastic reduction in the number of integration points, with important performance implications [74,75], particularly in nonlinear applications [72,76]. In these methods, the positiveness of the resulting weights is a prerequisite for the stability of the models. For this reason, recent contributions have focused on non-negative moment fitting [77,78]. In the context of this work, the block-diagonal mass matrices of enriched elements ought to similarly be positive definite. Additionally, the Gauss points have to match the element nodes in order to perform nodal quadrature (see also [56]). To allow for these constraints, a residual error was allowed in the moment-fitting equations, and weights were computed via the solution of a quadratic programming problem in order to guarantee the properties of the mass matrix while minimizing the integration error.

The remainder of this work is organized as follows: In Section 2, the weak form of the elastodynamics problem is reviewed. In Section 3, the novel method is presented and studied, with a brief review of the SEM offered in Section 3.1, followed by a presentation of the PU enrichment functions used in this work in Section 3.2. Building on the nodal quadrature method (Section 3.4), the proposed mass lumping technique for enriched SEs is presented in Section 3.6, and its accuracy and stability are studied. In Section 4, this method is deployed in three numerical examples to highlight its capabilities. In Section 4.1 a cracked panel is used to benchmark the accuracy of the method in comparison with the consistent mass matrix (CMM) formulation. In Section 4.2 a similar comparison is offered through the study of a riveted aluminum plate, involving multiple complex damage configurations. In Section 4.3, the effectiveness of the numerical modeling is assessed by simulating composite delaminations of a glass-fiber-reinforced polymer (GFRP) plate, for which reference experimental data were obtained via a scanning laser Doppler vibrometer (SLDV). Concluding remarks are presented in Section 5.

2. Problem Statement and Main Strategy

In this contribution, we consider a linear elastic domain $\Omega$ in $d=\{1,2,3\}$ dimensions in the presence of complex geometrical features and/or damage (Figure 1a). These features are either voids or cracks, and their boundaries are considered stress-free; however, they are allowed to intersect. The governing equations can be expressed in weak form as

$${\int }_{\Omega }\rho \ddot{\mathit{u}}\left(t\right)·\mathit{v}\mathrm{d}\Omega +{\int }_{\Omega }\mathit{\sigma }\left(\mathit{u}\left(t\right)\right):\mathit{ϵ}\left(\mathit{v}\right)\mathrm{d}\Omega ={\int }_{{\Gamma }_{\mathrm{s}}}{\mathbf{p}}_{\mathrm{s}}\left(t\right)·\mathit{v}\mathrm{d}{\Gamma }_{\mathrm{s}}+{\int }_{{\Gamma }_{\mathrm{cs}}}{\mathbf{p}}_{\mathrm{c}}\left(t\right)·\mathit{v}\mathrm{d}{\Gamma }_{\mathrm{cs}},$$

(1)

where $\dot{\square }$ denotes time differentiation; $\mathit{u}\left(t\right)$ is the displacement solution at time t,

$${\mathcal{U}}_t=\left\{\mathit{u}\left(t\right)|\mathit{u}\left(t\right)\in {\left(H^1\left(\Omega \right)\right)}^d,\mathit{u}\left(t\right)=\overline{\mathit{u}}\mathrm{on}{\Gamma }_{\mathrm{u}}\right\},$$

(2)

conforming to Dirichlet boundary conditions $\mathit{u}\left(t\right)=\overline{\mathit{u}}$ defined on ${\Gamma }_{\mathrm{u}}$; v is the test or weighting function [17],

$${\mathcal{V}}^0=\left\{\mathit{v}|\mathit{v}\in {\left(H^1\left(\Omega \right)\right)}^d,\mathit{v}=0\mathrm{on}{\Gamma }_{\mathrm{u}}\right\};$$

(3)

$\rho$ is the material density; $\mathit{\sigma }$ is the Cauchy stress tensor; and $\mathit{ϵ}$ is the linear strain. In order to facilitate the modeling process, mesh-independent geometrical descriptions of the domain details and damage are introduced so that, in the ideal case, a structured SE mesh can be used to discretize the domain. In the finite cell method (FCM) [79,80,81] and the SCM, this might also be considered as complementing the physical domain $\Omega$ with a void domain ${\Omega }_{\mathrm{v}}$, which can become relevant in addressing the poor conditioning of the system of equations.

3. Moment Fitting for Enriched Spectral Elements

3.1. The Spectral Element Method

While low-order standard Lagrangian FE schemes have seen widespread adoption due to their robustness and simplicity, Gibbs (or Runge) phenomena (see, e.g., [82] Chapter 4.2 and [83] Chapter 3) must be overcome as the polynomial order is increased in pursuit of more efficient discretizations. In this context, the SEM might be considered as a version of the FEM with more sophisticated choices of interpolating polynomials, which are produced by embedding additional properties in the determination of their supporting nodes. For instance, a Chebyshev nodal distribution [32,84,85] can be derived by minimizing the interpolation error. With two nodes fixed at the boundaries of the domain, internal Gauss–Lobatto–Legendre (GLL) nodes can be derived as optimal quadrature points in one dimension [17,23]. It can also be shown that the node location corresponds with the maximum of the respective shape function, with a value of one [83,86]. We refer the interested reader to [33,83] for these derivations and their applications in the SEM. For reasons that are highlighted in Section 3.4, the GLL version of the SEM has been extensively applied in structural dynamics [25,34,87], including previous versions of the SCM [54,56,88], and is also adopted herein.

Consider a generic one-dimensional GLL SE of order $p>1$ with $n=p+1$ nodes. While the two vertices delimit the element boundaries in the reference system $\xi \in [-1,1]$, the remaining nodes are chosen as the roots of the Lobatto polynomial $L_{p-1}$ of order $p-1$. Thus, the element nodes can be defined as the solution of

$$\left(1-{\xi }^2\right)L_{p-1}\left(\xi \right)=0,$$

(4)

which is known as the completed Lobatto polynomial [89]. Then, shape functions can be constructed as Lagrangian interpolations supported at these nodes:

$$N_i\left(\xi \right)=\prod _{j=1,j\ne i}^{p+1}\frac{\xi -{\xi }_j^p}{{\xi }_i^p-{\xi }_j^p},$$

(5)

where i is the node number. The local coordinates and shape functions for quadrilateral and hexahedral elements can be generated by applying Equations (4) and (5) in each dimension of the reference system and taking the sparse products of the respective results. For a generic element in d dimensions, the unknown displacement field $\mathbf{u}(\mathit{\xi },t)$ at a given time t is interpolated from its nodal values ${\mathbf{u}}_i\left(t\right)$, which are collected in the vector ${\mathbf{u}}_{\mathrm{e}}\left(t\right)$:

$$\mathbf{u}(\mathit{\xi },t)=\sum _{i=1}^n{N}_i\left(\mathit{\xi }\right){\mathbf{u}}_i\left(t\right)=\mathbf{N}\left(\mathit{\xi }\right){\mathbf{u}}_{\mathrm{e}}\left(t\right),$$

(6)

with

$$\mathbf{N}\left(\mathit{\xi }\right)=\left[\begin{array}{cccc}{N}_1{\mathbb{I}}_d& N_2{\mathbb{I}}_d& \dots & N_n{\mathbb{I}}_d\end{array}\right]$$

(7)

and ${\mathbb{I}}_d$ being the identity matrix.

3.2. The Extended Finite Element Method

In the extended or generalized FEM (XFEM/GFEM) [41,90], the polynomial solution space of the underlying mesh is supplemented by additional enrichment functions, thus enabling one to model the singular and/or discontinuous character ofdamage [41,42] and material interfaces [45,91]. Herein, only damage in the form of cracks and delaminations is considered; thus, enrichment is limited to discontinuous and singular functions. Within the XFEM framework, cracks are typically represented implicitly by means of two level-set functions:

• The normal level set $\varphi$, representing the signed distance from the crack surface:

  $$\varphi \left(\mathbf{x}\right)=\underset{\overline{\mathbf{x}}\in {\Gamma }_c}{min}∥\mathbf{x}-\overline{\mathbf{x}}∥\mathrm{sign}\left({\mathbf{n}}^{+}·\left(\mathbf{x}-\overline{\mathbf{x}}\right)\right),$$

  (8)

  where ${\mathbf{n}}^{+}$ is the outward normal to the crack surface and $\mathrm{sign}\left(\square \right)$ denotes the sign function.
• The tangential level set $\psi$, representing a signed distance function satisfying the conditions

  $$\nabla \varphi ·\nabla \psi =0,$$

  (9)

  $$\left.\begin{array}{c}\varphi \left(\mathbf{x}\right)=0\\ \psi \left(\mathbf{x}\right)=0\end{array}\right\}\forall \mathbf{x}\in {\Gamma }_f,$$

  (10)

  where ${\Gamma }_f$ represents the crack front/tips.

Using these functions, a polar coordinate system, with its origin at the crack tip/front, can be defined:

$$r=\sqrt{{\varphi }^2+{\psi }^2},\theta =arctan\left(\frac{\varphi }{\psi }\right).$$

(11)

It is common practice for level-set functions to only be evaluated at nodal points and FE shape functions to be used to interpolate these values in element interiors. Since this provides a relationship between the level sets/polar coordinates and the local element coordinates $\mathit{\xi }$, in what follows, these terms will be used interchangeably. Then, enrichment functions can be defined using the level sets of polar coordinates. For cracks and delaminations, two types of enrichment are necessary: discontinuous functions to represent displacement jumps along crack/delamination faces; and asymptotic functions to represent the singularities at the crack tips/fronts. For the former, the Heaviside function is typically used:

$$H\left(\varphi \right)=\left\{\begin{array}{cc}-1\hfill & \mathrm{if}\varphi <0\hfill \\ 1\hfill & \mathrm{if}\varphi \ge 0.\hfill \end{array}\right.$$

(12)

For the latter, a set of four [41] or twelve asymptotic functions [92] are typically used for strong discontinuities in homogeneous materials or bimaterial interfaces, respectively. Although there is a vast body of research that applies these principles in the study of composite delamination (e.g., [93,94,95,96,97]), the present contribution is concerned with modeling small areas of damage at minimal cost, rather than studying the phenomena of fracture mechanics in detail. Therefore, the asymptotic behavior in the immediate vicinity of the crack/delamination is not of interest. It has been shown [98] that acceptable accuracy can also be obtained using only the first singular enrichment function or even with enrichment functions that do not include singularities [99]. For these reasons, a single enrichment function is used herein, which also allows us to simplify the mass lumping process:

$$F_{asy}\left(r,\theta \right)=\sqrt{r}sin\left(\frac{\theta }{2}\right).$$

(13)

To facilitate the assignment of different enrichment functions to nodes, the following nodal sets are defined:

${\mathcal{N}}^j$

is the set of nodes enriched with the discontinuous function (Equation (12)). This set includes all nodes belonging to elements that are split in two by cracks but not those that contain crack tips/fronts.

${\mathcal{N}}^t$

is the set of nodes enriched with the asymptotic function (Equation (13)). This set includes all nodes belonging to elements that contain crack tips/fronts.

The above definition implies a topological enrichment scheme, as opposed to geometrical enrichment, where all nodes within a certain distance from the crack tips/fronts would be enriched [100]. While geometrical enrichment can recover optimal convergence rates in static problems, topological enrichment is preferred herein, since the proposed mass lumping method is not variationally consistent, and thus we wish to minimize the number of enriched elements in order to benefit from optimal lumping over the largest possible portion of the domain. Moreover, since the nodal sets defined above are disjointed, only one enrichment function $F_i\left(\mathit{\xi }\right)$ is assigned to each node i:

$$F_i\left(\mathit{\xi }\right)=\left\{\begin{array}{cc}H\left(\mathit{\xi }\right)\hfill & \mathrm{if}i\in {\mathcal{N}}^j\hfill \\ F_{asy}\left(\mathit{\xi }\right)\hfill & \mathrm{if}i\in {\mathcal{N}}^t.\hfill \end{array}\right.$$

(14)

Then, the displacement approximation presented in Equation (6) is extended as follows:

$$\mathbf{u}(\mathit{\xi },t)=\sum _{i=1}^n{N}_i\left(\mathit{\xi }\right){\mathbf{u}}_i\left(t\right)+\sum _{j=1}^{n_{\mathrm{e}}}{N}_j\left(\mathit{\xi }\right)F_j\left(\mathit{\xi }\right){\mathbf{a}}_j\left(t\right),$$

(15)

where $n_{\mathrm{e}}$ is the number of enriched nodes within the element, i.e., nodes belonging to ${\mathcal{N}}^j$ or ${\mathcal{N}}^t$; $F_j$ is the (only) enrichment function applied to node j; and ${\mathbf{a}}_j\left(t\right)$ collects the unknown enrichment parameters relating to $F_j$, which now appear in ${\mathbf{u}}_{\mathrm{e}}\left(t\right)$ and in the solution vector ${\mathbf{u}}_{\mathrm{s}}\left(t\right)$. We should underline that ${\mathbf{a}}_j\left(t\right)$ is not established a priori, rather it results from the solution of the discretized equilibrium Equations (17) or (21), which will be introduced shortly. We refer the interested reader to the aforementioned references. According to Equation (15), the element shape function matrix (Equation (7)) is modified as follows:

$$\mathbf{N}\left(\mathit{\xi }\right)=\left[\begin{array}{ccccc}{\mathbb{I}}_d{N}_1& {\mathbb{I}}_d{N}_1{F}_1& {\mathbb{I}}_d{N}_2& \dots & {\mathbb{I}}_d{N}_n{F}_n\end{array}\right].$$

(16)

3.3. Discretized Equilibrium Equations

Upon the substitution of Equation (15) in the weak form (Equation (1)), the discretized equilibrium equations are obtained:

$$\mathbf{M}{\ddot{\mathbf{u}}}_{\mathrm{s}}\left(t\right)+\mathbf{K}{\mathbf{u}}_{\mathrm{s}}\left(t\right)={\mathbf{f}}_{\mathrm{s}}\left(t\right).$$

(17)

The mass matrix $\mathbf{M}$, the stiffness matrix $\mathbf{K}$, and the system’s force vector ${\mathbf{f}}_{\mathrm{s}}$ are assembled form the respective element contributions [17]:

$${\mathbf{M}}_{\mathrm{e}}={\int }_{{\Omega }_{\mathrm{e}}}\rho {\mathbf{N}}^{\mathrm{T}}\mathbf{N}\mathrm{d}{\Omega }_{\mathrm{e}}$$

(18)

$${\mathbf{K}}_{\mathrm{e}}={\int }_{{\Omega }_{\mathrm{e}}}{\mathbf{B}}^{\mathrm{T}}\mathbf{D}\mathbf{B}\mathrm{d}{\Omega }_{\mathrm{e}}$$

(19)

$${\mathbf{f}}_{\mathrm{e}}\left(t\right)={\int }_{{\Gamma }_{\mathrm{s}}}{\mathbf{N}}^{\mathrm{T}}{\mathbf{p}}_{\mathrm{s}}\left(t\right)\mathrm{d}{\Gamma }_{\mathrm{s}}+{\int }_{{\Gamma }_{\mathrm{cs}}}{\mathbf{N}}^{\mathrm{T}}{\mathbf{p}}_{\mathrm{c}}\left(t\right)\mathrm{d}{\Gamma }_{\mathrm{cs}}$$

(20)

where the expression in Equation (19) implies the application of Hooke’s constitutive law, with $\mathbf{B}$ collecting the nodal displacement derivatives and matrix $\mathbf{D}$ governing the stress–strain relations. Finally, the element loads result from integrating tractions ${\mathbf{p}}_{\mathrm{s}}$, ${\mathbf{p}}_{\mathrm{c}}$ on domain boundaries ${\Gamma }_{\mathrm{s}}$ and implicitly defined loading surfaces ${\Gamma }_{\mathrm{cs}}$, respectively.

For materials such as composites, a damping term $\mathbf{C}{\dot{\mathbf{u}}}_{\mathrm{s}}\left(t\right)$ might be added to the equations of motion:

$$\mathbf{M}{\ddot{\mathbf{u}}}_{\mathrm{s}}\left(t\right)+\mathbf{C}{\dot{\mathbf{u}}}_{\mathrm{s}}\left(t\right)+\mathbf{K}{\mathbf{u}}_{\mathrm{s}}\left(t\right)={\mathbf{f}}_{\mathrm{s}}\left(t\right).$$

(21)

As the physics described by the damping matrix $\mathbf{C}$ is complex, Rayleigh’s assumption [101,102] is often applied in order to simplfy the computations:

$$\mathbf{C}=\alpha \mathbf{M}+\beta \mathbf{K}$$

(22)

where $\alpha ,\beta \in {\mathbb{R}}^{+}$ are material-dependent constants. Since matrix $\mathbf{C}$ has the same sparsity structure as $\mathbf{K}$, a substantial increase in the memory requirements of the model and the number of operations necessary for time integration ensues. In explicit dynamics, this can be overcome by performing mass-proportional damping, i.e., setting $\beta =0$, thus producing a diagonal damping matrix.

3.4. Mass Lumping for Standard Spectral Elements

In addition to the construction of an effective Ansatz space, GLL nodes can also be used as evaluation points within the Lobatto, or Gauss–Lobatto, integration rule [103]. This is an important prerequisite of the nodal quadrature technique [23,24,25], which consists in applying an integration rule defined at the element nodes for the evaluation of Equation (18), thus enforcing the variational formulation whilst producing a diagonal mass matrix. In contrast with other mass lumping techniques, this approach has the great advantage of preserving the order of convergence of the SEM, provided that the quadrature is of sufficient accuracy. The effectiveness of the GLL configuration stems from the fact that, for an SE of degree p, the order of the corresponding rule is $2p-1$: two degrees less than the Gauss–Legendre (GL) rule (which is a ne plus ultra in one dimension [104]) but far better than what can be achieved with Newton–Cotes formulae [86,103] relying on equally spaced nodes. Still, this approach leads to a slight under-integration of the stiffness and mass matrices, and thus to some loss of convergence in the computation of distorted elements [105]. For this reason, we employ a GL rule for the evaluation of Equation (19), while nodal quadrature is applied only to the mass matrix (Equation (18)) for the purpose of lumping. It has been shown in several studies [25,86,106] that this is effective in preserving the convergence rates related to SE interpolation. To briefly review the method, let us denote

$${\int }_{{\Omega }_{\mathrm{e},\mathrm{GLL}}}f\left(\mathit{\xi }\right)\mathrm{d}{\Omega }_{\mathrm{e},\mathrm{GLL}}$$

(23)

as the numerical evaluation of the integral of a function f over the GLL-SE domain ${\Omega }_{\mathrm{e}}$ by means of a quadrature rule that employs its nodes. This technique exploits the fact that the shape functions are rendered orthogonal by the Kronecker delta property (i.e., $N_j\left({\mathit{\xi }}_k\right)={\delta }_{jk}$), as long as their support nodes coincide with the quadrature nodes in evaluating the following integral:

$${\int }_{{\Omega }_{\mathrm{e}}}{N}_j\left(\mathit{\xi }\right)N_k\left(\mathit{\xi }\right)\mathrm{d}{\Omega }_{\mathrm{e}}\approx {\int }_{{\Omega }_{\mathrm{e},\mathrm{GLL}}}{N}_j\left(\mathit{\xi }\right)N_k\left(\mathit{\xi }\right)\mathrm{d}{\Omega }_{\mathrm{e},\mathrm{GLL}}={\delta }_{jk}{w}_j\left|\begin{array}{c}{J}_{\mathrm{e}}\left({\mathit{\xi }}_j\right)\end{array}\right|,$$

(24)

where the approximation sign (≈) highlights the aforementioned slight under-integration of the product $N_j^2\left(\mathit{\xi }\right)$, ${\mathit{\xi }}_j$ are the reference nodal coordinates in the GLL configuration, $w_j$ are the corresponding integration weights, and $\left|\begin{array}{c}{J}_{\mathrm{e}}\end{array}\right|$ is the determinant of the element Jacobian. By combining the above with Equations (7) and (18), a diagonal mass matrix with the following coefficients is produced:

$$m_{jk}\approx \sum _{i=1}^n{w}_i\rho \left({\mathit{\xi }}_i\right)N_j\left({\mathit{\xi }}_i\right)N_k\left({\mathit{\xi }}_i\right)\left|\begin{array}{c}{J}_{\mathrm{e}}\left({\mathit{\xi }}_i\right)\end{array}\right|={\delta }_{jk}{w}_j\rho \left({\mathit{\xi }}_j\right)\left|\begin{array}{c}{J}_{\mathrm{e}}\left({\mathit{\xi }}_j\right)\end{array}\right|.$$

(25)

It is important to note that the application of this strategy is restricted to elements with a sufficiently smooth mass density and displacement field. When discontinuous integrands appear in Equation (18), customized quadrature rules must be employed (see Section 3.5), which lead to a fully populated mass matrix. When dealing with voids ($\rho \left(\mathit{\xi }\right)=0$), this can be overcome by additionally performing the diagonal scaling of this matrix [52,88] or using moment-fitting equations to adapt the weights $w_i$ to the element’s cut configuration [56]. In addition, when enrichment functions are present, one is faced with the task of establishing appropriate diagonal mass coefficients for the enrichment parameters (see Equations (15) and (16)). For this, a procedure based on the conservation of kinetic energy has often been used [7,8,60]. This approach has also been combined with the Chebyshev version of the SEM [9]. In this contribution, we follow the indications of Żak and Krawczuk [34] in choosing a GLL nodal distribution instead, which makes optimal lumping available for standard SEs.

3.5. Element Partitioning

The approach described in Figure 1b as well as in Section 3.2 effectively shifts the burden of modeling some discontinuities (voids and cracks in our case) from the meshing phase to the integration of the weak form. In other words, the integrals described by Equations (18)–(20) are now discontinuous and/or singular, and thus they must be evaluated by an appropriate integration rule. This subject has seen substantial research interest and progress in the last few decades (see, e.g., [49,69,72,107,108,109]). The solution used in this contribution is documented in [56] for the case of voids. For the present work, the same algorithm is applied recursively for each level set in order to produce crack-conforming element partitions.

3.6. Mass Lumping for Enriched Spectral Elements

In this section, a mass lumping method for enriched elements is derived by combining nodal quadrature with a moment-fitting procedure that tailors the integration weights to the character of the integrand. The starting point is to combine Equations (16), (18) and (24). In one dimension, the mass coefficients of an enriched node i are then given by:

$${\mathbf{M}}_i=\left[\begin{array}{cc}{m}_{i,11}& m_{i,12}\\ m_{i,12}& m_{i,22}\end{array}\right]\approx {\int }_{{\Omega }_{\mathrm{e},\mathrm{GLL}}}\rho \left(\mathit{\xi }\right)\left[\begin{array}{cc}{N_i}^2\left(\mathit{\xi }\right)& N_i\left(\mathit{\xi }\right)F_i\left(\mathit{\xi }\right)\\ N_i\left(\mathit{\xi }\right)F_i\left(\mathit{\xi }\right)& {N_i}^2\left(\mathit{\xi }\right)F_i^2\left(\mathit{\xi }\right)\end{array}\right]\mathrm{d}{\Omega }_{\mathrm{e},\mathrm{GLL}},$$

(26)

i.e., although not fully diagonal, the mass matrix of enriched elements assumes an advantageous block-diagonal structure. In the spirit of nodal quadrature, our goal is now to evaluate Equation (26) with the best possible accuracy. For this, the three mass coefficients $m_{i,jk}$, with $j,k\in \{1,2\}$, can be considered as distinct integrals of functions $f_{i,jk}$:

$$m_{i,jk}\approx {\int }_{{\Omega }_{\mathrm{e},\mathrm{GLL}}}{f}_{i,jk}\left(\mathit{\xi }\right)\mathrm{d}{\Omega }_{\mathrm{e},\mathrm{GLL}}=\sum _{l=1}^n{f}_{i,jk}\left({\mathit{\xi }}_l\right)w_{l,jk}=f_{i,jk}\left({\mathit{\xi }}_i\right)w_{i,jk},$$

(27)

where n is the number of coinciding element and quadrature nodes, and thus the Kronecker delta property can be applied to achieve the last expression on the right side. As these nodes (${\mathit{\xi }}_i$) are necessarily fixed, it only remains to establish three independent weight sets ${\mathbf{w}}_{jk}$ that are best-suited to integrate the respective function sets ${\mathcal{F}}_{jk}=\left\{f_{1,jk},f_{2,jk},\dots ,f_{n,jk}\right\}$. For this task, we first decompose these integrands into m basis functions ${\mathcal{G}}_{jk}=\left\{g_{1,jk},g_{2,jk},\dots ,g_{n,jk}\right\}$, such that, in the ideal case,

$$f_{i,jk}\in span\left\{g_{1,jk},g_{2,jk},\dots ,g_{m,jk}\right\}\forall i.$$

(28)

Then, by applying the above to Equation (27), the relation between a set of weights ${\mathbf{w}}_{jk}$ and basis functions ${\mathcal{G}}_{jk}$ can be expressed in terms of the moment-fitting equations

$$\left[\begin{array}{cccc}{g}_{1,jk}\left({\mathit{\xi }}_1\right)& g_{1,jk}\left({\mathit{\xi }}_2\right)& \dots & g_{1,jk}\left({\mathit{\xi }}_n\right)\\ g_{2,jk}\left({\mathit{\xi }}_1\right)& g_{2,jk}\left({\mathit{\xi }}_2\right)& & ⋮\\ ⋮& & \ddots & ⋮\\ g_{m,jk}\left({\mathit{\xi }}_1\right)& \dots & \dots & g_{m,jk}\left({\mathit{\xi }}_n\right)\end{array}\right]\left[\begin{array}{c}{w}_{1,jk}\\ w_{2,jk}\\ ⋮\\ w_{n,jk}\end{array}\right]=\left[\begin{array}{c}{\int }_{{\Omega }_{\mathrm{e}}}{g}_{1,jk}\left(\mathit{\xi }\right)\mathrm{d}{\Omega }_{\mathrm{e}}\\ {\int }_{{\Omega }_{\mathrm{e}}}{g}_{2,jk}\left(\mathit{\xi }\right)\mathrm{d}{\Omega }_{\mathrm{e}}\\ ⋮\\ {\int }_{{\Omega }_{\mathrm{e}}}{g}_{m,jk}\left(\mathit{\xi }\right)\mathrm{d}{\Omega }_{\mathrm{e}}\end{array}\right]\iff {\mathbf{A}}_{jk}{\mathbf{w}}_{jk}={\mathbf{b}}_{jk},$$

(29)

where the repetition of indices does not imply summation, and the choice of basis functions ${\mathcal{G}}_{jk}$ is reported in

Table 1. Definition of basis functions ${\mathcal{G}}_{jk}$ for three instances of the moment-fitting Equation (29), whose results are applied in the integration of the respective function sets ${\mathcal{F}}_{jk}$. $\rho$ is assumed as constant within ${\Omega }_{\mathrm{e}}$.

$\mathit{jk}$	11	12	22
${\mathcal{F}}_{jk}$	$\rho \times \left\{{N_1}^2,{N_2}^2,\dots ,{N_n}^2\right\}$	$\rho F\times \left\{N_1,N_2,\dots ,N_n\right\}$	$\rho F^2\times \left\{{N_1}^2,{N_2}^2,\dots ,{N_n}^2\right\}$
${\mathcal{G}}_{jk}$	$\left\{1,\xi ,{\xi }^2,\dots ,{\xi }^m\right\}$	$F\times \left\{1,\xi ,{\xi }^2,\dots ,{\xi }^m\right\}$	$F^2\times \left\{1,\xi ,{\xi }^2,\dots ,{\xi }^m\right\}$

. Three considerations are due at this point. Firstly, we should note that the number of basis functions m is limited by the number of nodes n, as the system would be overdetermined for $m>n$. This is in contrast to most other instances of this method, which exploit an under-determined system [71,77]. For this reason, only the basis monomials up to order m are considered, and the higher-order terms appearing in the diagonal integrands ${\mathcal{F}}_{11},{\mathcal{F}}_{22}$ will be under-integrated, as they are not represented in Equation (29). For standard SEs, however, this choice still leads to weights that match the GLL rule, and thus to only minor under-integration of the mass matrix [56]. Secondly, we should mention that, in this work, $\rho$ is considered constant within the physical portion of the element, and the effect of voids is accounted for by evaluating the vector of moments (${\mathbf{b}}_{jk}$) in a consistent sense, according to one of the procedures described in Section 3.5. If no voids are present, then ${\mathbf{w}}_{11}$ corresponds to the standard GLL weights.

A third important remark is that a positive definite mass matrix is required for explicit solvers to converge; therefore, negative weights must be avoided when $j=k$ in Equation (29). Similarly, certain values of the off-diagonal weights ${\mathbf{w}}_{12}$ might lead to complex eigenvalues for some of the blocks ${\mathbf{M}}_i$. These issues call for the application of box constraints on the weights, which cannot be applied directly to Equation (29). To overcome this, a residual is allowed in its evaluation:

$${\mathbf{r}}_{jk}={\mathbf{A}}_{jk}{\mathbf{w}}_{jk}-{\mathbf{b}}_{jk}.$$

(30)

Then, the computation of the weights can be expressed as an optimization problem, seeking to minimize the $L_2$ norm of this residual, leading to the following quadratic programming problem:

$$\underset{\mathbf{w}}{\mathrm{minimize}}\frac{1}{2}{\mathbf{w}}^{\mathrm{T}}\overline{\mathbf{A}}\mathbf{w}-{\mathbf{w}}^{\mathrm{T}}\overline{\mathbf{b}}$$

(31)

where $\overline{\mathbf{A}}={\mathbf{A}}^{\mathrm{T}}\mathbf{A}$ and $\overline{\mathbf{b}}={\mathbf{A}}^{\mathrm{T}}\mathbf{b}$. In Ref. [56], we discussed appropriate constraints for the set ${\mathbf{w}}_{11}$ that are necessary in the presence of voids. Constraints for the fitting problem in enriched elements are formulated herein. From the evaluation of Equation (26), we have

$${\mathbf{M}}_i\approx \rho \left[\begin{array}{cc}{w}_{i,11}& w_{i,12}{F}_i\\ w_{i,12}{F}_i& w_{i,22}{F}_i^2\end{array}\right],$$

(32)

where $F_i=F\left({\xi }_i\right)$. To ensure the positive definiteness of this block, its eigenvalues ${\lambda }_{1,2}$ must be positive [110]. They can be computed as

$${\lambda }_{1,2}=\frac{\rho }{2}\left[\begin{array}{c}{w}_{i,11}+w_{i,22}{F}_i^2\mp \sqrt{{(w_{i,22}{F}_i^2-w_{i,11})}^2+4w_{i,12}^2{F}_i^2}\end{array}\right].$$

(33)

Due to their impact on stability and the properties of the mass matrix, we wish to formulate the constraints for Equation (31) in terms of a smallest admissible eigenvalue ${\lambda }_{\min }\in {\mathbb{R}}^{+}$. From Equation (33), the following condition for ${\lambda }_{1,2}>{\lambda }_{\min }$ is derived in Appendix A:

$$w_{i,22}>\frac{{\lambda }_{\min }^2-{\lambda }_{\min }\rho w_{i,11}-{\rho }^2{w}_{i,12}^2{F}_i^2}{\rho F_i^2({\lambda }_{\min }-\rho w_{i,11})},$$

(34)

where $w_{i,11}$ is readily available, and $w_{i,12}$ can be computed with Equation (29). The natural upper bound for ${\lambda }_{\min }$ is ${\lambda }_{\min }<\rho w_{i,11}$, thus, appropriate values can be computed as follows:

$${\lambda }_{\min }={ϵ}_{\lambda }\rho \underset{0<i\le n}{min}\left\{w_{i,11}\right\},$$

(35)

where a tuning parameter ${ϵ}_{\lambda }\in \left[0,1\right)$, which will be established in Section 3.6.1 and Section 3.6.2, is multiplied by the mass density and the smallest weight of the set ${\mathbf{w}}_{11}$. In essence, Equations (30), (34) and (35) aim at a compromise between an advantageous mass matrix structure and an accurate evaluation of Equation (18) within the aforementioned restrictions. As long as the ordering of standard and enriched DOF related to a node is contiguous (see Appendix B), the mass matrix has a $2\times 2$ block structure in enriched regions of the domain and is diagonal elsewhere. In linear elastodynamics, Equations (18)–(20) can be assembled once and then reused; thus, the computational cost overhead introduced by the moment-fitting procedure is negligible.

These advantages, however, come with several caveats. This method is restricted to the use of one enrichment function per node, as the emergence of additional interaction coefficients in Equation (26) quickly complicates the formulation of constraints. Using different enrichment functions within an element is still possible, as long as these are applied on different nodes. Secondly, commonly used shifted enrichment functions [111,112] cannot be applied, since they would lead to zero mass coefficients. In general, it should be noted that the variational formulation is abandoned by introducing a residual in Equation (30). For this reason, the accuracy of the novel procedure will be studied at both the element level (in the upcoming Section 3.6.1) and the system level (in Section 4.1)). Finally, we expect the presence of discontinuities and the manipulation of the mass matrix to have an important impact on stability, which is studied in Section 3.6.2.

3.6.1. Accuracy

Figure 2 displays the reference system of a jump-enriched 1D element of unit density and order $p=5$, along with the standard ($N_3\left(\xi \right)$) and enriched ($N_3\left(\xi \right)H\left(\xi \right)$) interpolation functions related to its third node. In this section, we study the accuracy of the procedure in its intended application, namely the integration of the terms $m_{i,12}$ and $m_{i,22}$ (see Equation (27)). To this end, the integration error ${ϵ}_{\mathrm{I}}$ is computed with respect to a reference value $I_{\mathrm{ref}}$, which is evaluated with a boundary-conforming quadrature rule, according to Section 3.5:

$${ϵ}_{\mathrm{I}}=\left|\frac{f_{i,jk}\left({\mathit{\xi }}_i\right)w_{i,jk}-I_{\mathrm{ref}}}{I_{\mathrm{ref}}}\right|.$$

(36)

In Figure 3, ${ϵ}_{\mathrm{I}}$ is reported for different element orders p, integrands $f_{i,jk}={\left(N_{p,i}H\right)}^j$, crack locations $l_{\mathrm{c}}$, and a fixed value of ${ϵ}_{\lambda }=1-{10}^{-4}$. We should note that, by choosing H as an enrichment function, the integrands (as well as the moments in Equation (29)) are discontinuous only for the case with $j=1$. In general, it can be observed that the procedure is less successful whenever the discontinuity lies near an element node: an effect that is most pronounced in the interaction terms for $p=3$ (Figure 3a). With this exception, the accuracy achieved in the off-diagonal terms is satisfactory, as shown in Figure 3b,c. For diagonal terms, however, (Figure 3d–f) the errors are extremely elevated. In this context, we should underline that this strategy is first and foremost a mass lumping method, and secondarily an integration method. In other words, the quality of the weight set ${\mathbf{w}}_{22}$ is heavily restricted by the constraints expressed by Equation (34), which are necessary in order to maintain a positive definite mass matrix. We should also note that, with few exceptions [25], several mass lumping strategies provide little to no a priori guarantees in terms of accuracy, and the validity of the resulting model will be studied in Section 4.1, Section 4.2 and Section 4.3. In Figure 4, the convergence of the integration error is studied with integrands of the form $H{\left(\xi \right)}^j(1+{\xi }^{p_{\mathrm{I}}})$ and a fixed discontinuity at $l_c=1/3$. Although errors in the order of ${10}^{-10}$ are not as low as machine precision, Figure 4a suggests that an integration order of p can be achieved for the weight set ${\mathbf{w}}_{12}$. This is consistent with the performance of other moment-fitting procedures that do not optimize the Gauss point location [69,75] and thus perform similarly to Newton–Cotes quadrature rules [86,103]. In comparison, the error convergence for the weight set ${\mathbf{w}}_{22}$ (Figure 4b) is modest, and, given the simplicity of the integrand, can be again attributed to the presence of constraints, as the incompleteness of the fitting basis ${\mathcal{G}}_{22}$ (see Table 1) should not play a role in this case.

To gain further insight into this issue, we study the impact of the minimum mass eigenvalue ${\lambda }_{\min }$ on the accuracy of the weight set ${\mathbf{w}}_{22}$. Figure 5 displays the integration error ${ϵ}_{\mathrm{I}}$ [-] by varying ${ϵ}_{\lambda }$ for a selection of cases from Figure 3. For crack locations unfavorably near a node (cases (a) and (b)), this parameter has no significant impact, while in more favorable configurations (cases (c) and (d)), smaller values of ${ϵ}_{\lambda }$ allow for a better solution of Equation (31). Some local minima of ${ϵ}_{\mathrm{I}}$ [-] are sporadically detected for ${ϵ}_{\lambda }\in [0.1,1)$ (e.g., case (e)), which are likely related to the nonlinear character of Equation (34).

3.6.2. Stability

By intervening in the mass matrix, mass lumping methods have an impact on element eigenmodes and thus on the stability of explicit schemes. This concern is particularly relevant in the presence of discontinuities, which can lead to a large increase in eigenvalues, the extent of which we wish to quantify in this section. For central difference methods (CDMs), the Courant–Friedrichs–Lewy (CFL) condition is given conservatively by the smallest value of [113]:

$$\Delta t_{\mathrm{e}}=\frac{2}{{\omega }_{\max }}$$

(37)

across all elements. In Equation (37), $\Delta t_{\mathrm{e}}$ is the element critical time step, and ${\omega }_{\max }$ is the element spectral radius, i.e., its biggest eigenvalue obtained from [113]:

$$det\left({\mathbf{K}}_{\mathrm{e}}-{\omega }_{\mathrm{e}}^2{\mathit{M}}_{\mathrm{e}}\right)=0.$$

(38)

We study the stability of the procedure by considering a two-dimensional, unit-size SE traversed by a straight discontinuity parallel to the y axis and enriched by a modified Heaviside step function, similar to Figure 2. In Figure 6a,b, results are reported for different crack locations $l_{\mathrm{c}}$, comparing the CMM formulation with the proposed lumping method, where a fixed value of ${ϵ}_{\lambda }={10}^{-1}$ is chosen based on Figure 5. All results are scaled by $\Delta t_{\mathrm{e},0}$, i.e., the critical time step of a standard SE with $p=5$. While for the CMM the value of $\Delta t_{\mathrm{e}}$ remains in the same order of magnitude as $\Delta t_{\mathrm{e},0}$, the proposed procedure can lead to values that are from one to two orders of magnitude smaller. If the CDM were used for time integration, this would proportionally multiply the cost of a simulation, risking the negation of the benefits of mass lumping. However, if a leap-frog solver is employed, then this issue will affect only the enriched portion of the domain, which is usually rather small. It was shown in Ref. [56] that, with the CDM, a performance only slightly worse than that of the SEM could be achieved. It should also be noted that Figure 6a does not lead to a performance gain, since this formulation would require the use of an implicit solver. In Figure 6c the effect of ${ϵ}_{\lambda }$ is studied for a fixed configuration ($l_{\mathrm{c}}=0.3$). These results indicate that the constraint expressed in Equation (34) is active, as the spectral radius can be controlled by ${\lambda }_{\min }$. Finally, even if a leap-frog solver is used, setting ${ϵ}_{\lambda }<{10}^{-3}$ will lead to prohibitively expensive simulations for negligible accuracy gains (see Figure 5). For these reasons, we recommend using ${ϵ}_{\lambda }\in \left[{10}^{-2},{10}^{-1}\right]$, and a value of ${ϵ}_{\lambda }={10}^{-2}$ will be used in the following numerical examples.

4. Numerical Examples

In this section, three problems of increasing complexity are studied to assess the validity of the proposed method, which we label the extended spectral cell method (XSCM) for the remainder of this work. As mentioned previously and documented in Ref. [56], we deploy a second order leap-frog solver [18,19] whenever we use the (X)SCM. In the first 2D example (Section 4.1), a cracked panel in plane stress is considered, and the performance of the fitted mass lumping procedure is compared with a variationally consistent alternative. A true instance of the XSCM, involving both ’cut’ cells and enrichment functions, is then displayed in the 3D modeling of an aluminium plate with several details and complex damage configurations (Section 4.2), where a qualitative comparison with the SEM is offered. Finally, we propose an experimental validation of the method in Section 4.3, where our numerical model is compared with scanning laser Doppler vibrometer (SLDV) measurements of a glass-fiber-reinforced polymer (GFRP) panel under piezo-electric excitation. We refer the interested reader to [19] for a version of the leap-frog solver accounting for damped waves that is applied to the solution of this last example.

4.1. Example 1

Figure 7 displays an aluminum panel in plane stress with parameters $E=70$ GPa, $\nu =0.33$, and $\rho =2700\mathrm{kg}/{\mathrm{m}}^3$. Fixed boundary conditions are applied on its left end, while on the right-hand side the transient load

$${\mathbf{p}}_{\mathrm{s}}(\mathbf{x},t)=\mathbf{p}\left(\mathbf{x}\right)sin\left(\omega t\right){sin}^2\left(\frac{\omega t}{2n}\right),t\in \left[0;\frac{n}{f}\right]$$

(39)

is applied. In Equation (39), a spatially uniform load is given by $\mathbf{p}\left(\mathbf{x}\right)\equiv p{\mathbf{e}}_{\mathrm{x}}$, with ${\mathbf{e}}_{\mathrm{x}}={[1,0]}^T$ and $p={10}^6\mathrm{N}/\mathrm{m}$, while time modulation is produced by a Hann window of circular frequency $\omega =2\pi f$. The excitation window $n/f$ is determined by the carrier frequency $f=200$ kHz and the number of cycles $n=5$. This use of a smoothed signal is effective in producing narrow frequency band waves and thus limiting the dispersive effects [114]. The damage consists of an oblique, planar crack of length $\sqrt{2}a$, ($a=5$ cm), which is discretized according to Section 3.2. In this benchmark, the proposed method is compared to a CMM formulation paired with the implicit version of the Newmark algorithm [115]. For this, the displacement time history is recorded at sensors $s_1$ and $s_2$ positioned according to Figure 7, which are used to measure the accuracy of a simulation as follows:

$${ϵ}_{\mathrm{h},\mathrm{t}}=\sqrt{\frac{{\displaystyle \sum _{i=0}^{n_{\mathrm{t}}}}{\left[u_{\mathrm{h}}\left(i\Delta t\right)-u_{\mathrm{ref}}\left(i\Delta t\right)\right]}^2}{{\displaystyle \sum _{i=0}^{n_t}}{u}_{\mathrm{ref}}^2\left(i\Delta t\right)}},$$

(40)

where $n_{\mathrm{t}}$ is the number of time steps; $u_{\mathrm{h}}$ is the numerical solution at a sensor; and $u_{\mathrm{ref}}$ is a reference solution computed with the CMM formulation, discretized with $3.2\times {10}^6$ DOF, $p=6$, and $\Delta t={10}^{-10}$. The simulation time measures $0.1$ ms and, for all models to be tested, a fixed time step of $\Delta t={10}^{-9}$ s is used, which roughly represents the CFL condition of the finest meshes.

The convergence of ${ϵ}_{\mathrm{h},\mathrm{t}}$ at both sensors is reported for models of third to sixth orders with respect to the number of DOF in Figure 8a–d and the amount of time necessary for the solution of Equation (17) in Figure 8e–h. In the former set, the proposed method performed similarly to the CMM formulation; however, some loss of convergence emerged for the largest models. This likely coincided with the stage at which the integration errors discussed in Section 3.6 exceeded those of the SE interpolation and, thus, created an accuracy bottleneck. In the latter set (Figure 8e–h), the proposed approach increased the speed by roughly one order of magnitude for $p=3,4$ and more modestly for $p=5,6$. Although the CMM formulation paired with an implicit solver prevailed in terms of maximum accuracy, one should keep in mind that this required substantial additional memory storage for the factorization of the effective stiffness matrix. Finally, the same formulation was used for the reference solution, and thus the accuracy of the CMM could have been overstated for the largest models presented here.

4.2. Example 2

The aluminium plate illustrated in Figure 9a presents 34 cylindrical holes of radius $r_{\mathrm{h}}=5$ mm, whereas damage is represented by two planar cracks orthogonal to the xy plane. The bottom crack has a length of $21.2$ mm and an inclination of ${45}^{\circ }$, while the vertical, 10 mm long crack is connected to one of the holes. Piezo-electric excitation was modeled as a circular loading surface of diameter $2r_{\ell }=10$ mm. The loading followed a radial distribution in space with a maximum of p at the outside edge of the actuator and a value of zero at its center:

$$\mathbf{p}\left(\mathbf{x}\right)=\frac{p}{r_{\ell }}\mathbf{r}\left(\mathbf{x}\right),\parallel \mathbf{r}\left(\mathbf{x}\right)\parallel \le r_{\ell }$$

(41)

where $\mathbf{r}\left(\mathbf{x}\right)$ points from the center of the circle to $\mathbf{x}$, $\parallel \mathbf{r}\left(\mathbf{x}\right)\parallel$ is its Euclidean norm, and $p=2.5\times {10}^9$ Pa is the loading magnitude for $\parallel \mathbf{r}\left(\mathbf{x}\right)\parallel =r_{\ell }$. This assumption is similar to the pin-force model typically used in this context [116], where the loading is, however, applied only along the circumference of the actuator. Modulation in time was achieved via Equation (39) with $f=200$ kHz and $n=5$. Given a plate thickness of $t=1.5$ mm, with these parameters we expected to observe the first symmetric and anti-symmetric Lamb wave modes. It should also be noted that a high number of cycles n is effective in limiting wave dispersion, albeit at the cost of expanding the excitation window. This example highlights the benefits offered by the XSCM in the modeling of complex components and damage configurations. For the purpose of comparison with the SEM, a setup was chosen that still allowed for discretization with a conforming mesh. In both cases, a second-order mesh was produced with Gmsh [117] and then converted into a third-order GLL SE mesh via the shape functions of the original elements, leading to a sub-parametric representation of the domain. The differences between the meshes are shown in Figure 9b,c. With the SEM, a strucutred nonorthogonal mesh was employed in order to successfully discretize load surfaces, cracks, and holes with hexahedral elements. This required the use of a highly customized Gmsh script and would not be feasible in an automatized fashion should the location of a detail change. With the XSCM, an orthogonal mesh could be used instead, and all domain details, including loading surfaces, were defined at run-time. This enabled the optimal choice of the element size and automatized simulations for, e.g., damage detection applications [118,119,120,121,122]. In both cases, the plates were one element thick and had an element size of $h=2.6$ mm in the x and y directions, leading to 31 nodes per wavelength of the asymmetric mode. As reported in

Table 2. DOF and critical time steps $\Delta t_{\mathrm{c}}$ for the models presented in Figure 9 and Figure 10. A time step of $\Delta t=5\times {10}^{-9}$ was used to ensure good accuracy. In the XSCM, a local time refinement of $p_{\mathrm{t}}=34$ enabled us to comply with the CFL condition [56].

Discretization	SEM Pristine	SEM Damaged	XSCM Damaged
$\Delta t_{\mathrm{c}}\left({10}^{-8}\mathrm{s}\right)$	$1.5$	$1.5$	$0.015$
$n_{\mathrm{dofs}}\left({10}^6\right)$	$0.54$	$1.08$	$0.93$

, the resulting models had a size of $1.1\times {10}^6$ DOF for the SEM and $9.3\times {10}^5$ DOF for the XSCM. This difference, in spite of the additional enrichment parameters in the latter model, was due to the fact that, in the former model, the element size had to be frequently reduced to resolve the complex geometry of the modeled structure (see Figure 9b). An undamaged version of the plate, which was added for comparison (’SEM pristine’) enabled us to exploit the symmetry of the model along the x axis, thus halving its size. All solutions employed a time step of $\Delta t=5\times {10}^{-9}\left[\mathrm{s}\right]$. This is well below the critical time step of the conforming meshes to promote accurate time integration. The enriched model had a relatively small $\Delta t_{\mathrm{c}}$, which was expected given the results of Figure 6b. In order to accommodate for this, the leap-frog solver [18,19] was deployed with a local time refinement ratio of $p_{\mathrm{t}}=34$ [56].

Figure 10 collects the time histories of the vertical displacements at the sensors. The simulation with the SEM in the absence of cracks (dotted line) was included for context. For the damaged case, good agreement between the XSCM and the SEM could be observed, despite the differences in these models: enrichment functions and a mass lumping error were introduced in the proposed method, while explicitly meshed details and some element distortion were encountered in the SEM. In light of the complexity of this problem, and of these fundamental modeling differences, these results speak of the quality of both approaches.

4.3. Example 3

In this section, we study a glass-fiber-reinforced polymer (GFRP) panel with dimensions of $500\times 500$ mm and a thickness of about 2 mm. The specimen consisted of twelve layers of VV192T/202 prepregs designed by G. Angeloni with orientations of ${[0/90/0/90/0/90]}_s$. The glass fiber fabric had a twill weave with a density of $202\mathrm{g}/{\mathrm{m}}^2$. Delaminations were artificially introduced at the manufacturing stage using Teflon inserts between layers of prepregs. In this particular case, the inserts were located in the middle of the cross-section, i.e., between the sixth and seventh laminae. Four delaminations with different shapes and locations were introduced according to

Table 3. Delamination shapes and locations with origin in the center of the plate.

Shape	a  (mm)	b  (mm)	x  (mm)	y  (mm)
Ellipse	30	20	$-125$	1
Circle	20	20	$2.45$	$-125$
Square	20	20	126	$-1.5$
Rectangle	30	20	1	126

. Their positions differed slightly with respect to the specifications provided by the manufacturer—all delamination centers were designated to be $125\mathrm{mm}$ from the edge of the plate and aligned with the center of the plate. The actual delamination locations were identified through signal processing using wavenumber damage imaging (see Figure 11a).

A piezoelectric transducer (Noliac, NCE51) of diameter 10 mm and thickness $0.5$ mm located at the center of the plate was used for the excitation of guided waves. A Hann windowed signal with a central frequency of $50\mathrm{kHz}$ and 5 cycles was applied to the piezoelectric transducer (20 Vpp). An SLDV was used to acquire Lamb wave signals on a uniform spatial grid (497 × 497 points) on the surface of the plate. The measured area was $496\times 496$ mm, because a 2 mm margin along the plate edges was omitted. Out-of-plane particle vibration velocities were measured using one laser head (Polytec PSV-400) [123] placed perpendicularly to the surface of the plate.

The plate was modeled in 3D by a single layer of SEs of order $p=3$ and element size of $h=5$ mm. The time integration step measured $\Delta t=1.8\times {10}^{-7}$ s, corresponding to the CFL condition of the standard portion of the mesh. A local time refinement ratio of $p_{\mathrm{t}}=31$ was used in order to ensure the stability of the enriched regions of the model. Establishing appropriate parameters for the GFRP material was not a straight-forward task. After assuming an orthotropic material law, the elastic constants of the GFRP lamina were identified using the inverse method described in Ref. [124], and are reported in

Table 4. Model parameters for numerical simulation of composite panel. Unit is GPa unless specified otherwise.

Parameters	${\mathit{C}}_{11}$	${\mathit{C}}_{22}$	${\mathit{C}}_{33}$	${\mathit{C}}_{12}$	${\mathit{C}}_{13}$	${\mathit{C}}_{23}$	${\mathit{C}}_{44}$	${\mathit{C}}_{55}$	${\mathit{C}}_{66}$	$\mathit{\alpha }$ (1/s)	$\mathit{\rho }$ (kg/m3)
GFRP	$44.82$	$27.97$	$3.26$	$14.13$	$6.97$	$1.26$	$4.91$	$4.22$	$3.03$	8000	1750

. A layered integration technique based on laminae-conforming element partitions was used to correctly evaluate Equation (19), which was also discontinuous across laminae [125]. For simplicity, the square and rectangular delaminations were modeled by circular and elliptic level sets, respectively. The mass and stiffness of the piezo-electric actuator were also neglected, along with the piezo-electric effect (the mass of the sensor represented roughly $0.03\%$ of the panel mass). Instead, the same loading conditions as in the previous example were applied, with $p=7.2\times {10}^7\left[\mathrm{Pa}\right]$. This value, as well as the mass-proportional damping coefficient $\alpha$, was established empirically by matching the amplitudes of the measured and simulated wave velocities. Finally, this material introduced some additional uncertainty with respect to the plate thickness, which could have presented some significant variation from the specified value.

In Figure 12, the time histories of the vertical velocities are reported for the locations described in Figure 11c. The initial wave packet could be reproduced at most sensors, with $s_{2,3}$ showing the most pronounced differences in phase and amplitude. These sensors represent velocities purely in x and y directions, respectively, and show the limitations of the aforementioned assumptions. Reflections from the damaged area and/or panel boundaries could realistically be reproduced in most cases, with only $s_{3,6}$ showing a significant difference in time of flight. Overall, given all the aforementioned limitations, these results can be considered satisfactory.

5. Conclusions

When the SEM is complemented with PU enrichment functions in pursuit of fast simulations, high-oder accuracy, and relief from meshing issues, preserving and combining all these qualities hinges on the effectiveness of the mass lumping of enriched SEs. In this contribution, a novel method based on the moment-fitting framework and the nodal quadrature method was presented: while a block-diagonal mass matrix was produced by evaluating its integral at the element nodes, integration weights were delivered by a nonlinear moment-fitting procedure. This was intended to minimize the mass lumping error under the constraints of a positive definite mass matrix. Due to these constraints, the mass matrix was under-integrated, and the variational formulation was violated by a degree that was quantified in Figure 3, Figure 4 and Figure 5. This approach offered some room to balance accuracy and stability thanks to Equation (34); however, some further research could be devoted to improving its stability (Figure 6b,c).

The method yielded an accuracy that was comparable to that of the CMM formulation for several models, with convergence loss emerging only for relatively large models ($5\times {10}^5$ DOF in two dimensions) (Figure 8a–d). By alleviating stability issues via a leap-frog solver [18,19,56] and leveraging the fact that only a portion of the domain was cut or enriched, this approach resulted in efficient explicit dynamics simulations (Figure 8e–h). The XSCM applied in Section 4.2 stemmed from enhancing the SCM [51,52,53,54,55,56] with PU enrichment functions and enabled the automatic generation of complex damage configurations (Figure 9c) while delivering comparable results to the conforming SEM, thus saving meshing effort. These damage modeling capabilities were finally applied to composite delamination, where good agreement between the model and experimental data could be achieved.