Variationally consistent mass scaling for explicit time-integration schemes of lower- and higher-order finite element methods

In this paper, we propose a variationally consistent technique for decreasing the maximum eigenfrequencies of structural dynamics related finite element formulations. Our approach is based on adding a symmetric positive-definite term to the mass matrix that follows from the integral of the traction jump across element boundaries. The added term is weighted by a small factor, for which we derive a suitable, and simple, element-local parameter choice. For linear problems, we show that our mass-scaling method produces no adverse effects in terms of spatial accuracy and orders of convergence. We illustrate these properties in one, two and three spatial dimension, for quadrilateral elements and triangular elements, and for up to fourth order polynomials basis functions. To extend the method to non-linear problems, we introduce a linear approximation and show that a sizeable increase in critical time-step size can be achieved while only causing minor (even beneficial) influences on the dynamic response.


Introduction
Explicit time-integration techniques are suitable when the time scales of interest are inherently small and for partial differential equations that are not very stiff. Industry relevant example applications are those involving contact, fragmentation and penetration, such as explosions and crashes, or quasi-static analysis of highly nonlinear processes [1,2]. While each time step in an explicit time-integration scheme is relatively inexpensive to compute, stability requirements limit the permissible time-step size. This so-called "critical time step" is related to the maximum eigenfrequency of the mass-to-stiffness generalized eigenvalue problem. The specific relation depends on the time-integration technique that is employed. A typical example for undamped structural dynamics is: which holds for the central difference method as a special case of the Newmark method [3, p. 493]. Here, λ h max is the maximum eigenvalue of the discrete formulation and ω h max is the corresponding maximum eigenfrequency. Figure 1 shows the well-known ratios between the discrete and analytical eigenfrequencies for a one-dimensional linear wave equation discretized with standard C 0 -continuous finite elements of polynomial orders P of 1 to 4. As this graph implies, the numerical eigenfrequencies of a standard Galerkin approximation always overpredict the analytical ones [3]. Importantly, the low-frequency part of the spectra is accurate, and they become more accurate with increasing polynomial order [3,4]. The spectra corresponding to the higher-order elements exhibit jumps, leading to different "branches". The first branch is called the acoustic branch, and is the most accurate. The other branches are called the optical branches [5]. The maximum frequency, which limits the critical time step, is severely overpredicted in these optical branches and diverges with increasing polynomial order. Various techniques exist to reduce maximum eigenfrequencies. These mostly center around modifying the mass matrix. Classical examples are mass-lumping methods, which operate on the consistent mass matrix (in the sense of [6,7]) to produce a diagonal mass matrix [3,8]. While the primary advantage of a diagonal mass matrix is that its inversion is trivial, a secondary effect of mass lumping is a reduction of the maximal eigenvalues [4]. Various lumping strategies exist, such as row-sum lumping (which may lead to negative diagonal values for higher-order serendipity elements, which are catastrophic for numerical simulation) [8], diagonal scaling [9], lumping by nodal quadrature [10,11,12] and manifoldbased methods [13]. The approximations underlying these methods imply that the final formulations are not variationally consistent, in the sense of [10,Ch. 4]. As a consequence, the true analytical solution corresponding to the strong form no-longer satisfies the discrete problem statement, resulting in sub-optimal convergence.
Other approaches aim to add some form of additional mass to the mass matrix. Such "mass-scaling techniques" may or may not be variationally consistent depending on the origin of this fictitious mass. Examples are [14,15,16,17], where a weighting of some form of stiffness matrix is used as a mass scaling, [18], where the addition to the mass matrix follows variationally from a penalized Hamilton's principle, and [19], where nodal accelerations are modified to decrease the maximal eigenfrequency for solid shell elements.
In the current article, we explore a mass-scaling technique based on a symmetric penalization of the natural transmission conditions across element interfaces. Such techniques are effective at reducing outlier frequencies due to reduced continuity across patch interfaces and on domain boundaries in an isogeometric analysis framework [20,21]. Penalties on the lack of continuity are also a common theme in stabilized methods, such as ghost-penalty for immersed finite element methods [22], edge-stabilization for inf-sup stability of equal-order pairs in mixed formulations [23], and discontinuous Galerkin methods in general [24]. In each of those cases, though, the additional terms alter the stiffness matrix rather than the mass matrix. This serves to increase the smallest eigenvalues. To reduce the largest eigenvalues, we add a similar term to the mass matrix.
The remainder of this article is structured as follows. In Section 2, we propose our mass-scaling method for the linear wave equation. We also develop an estimation strategy for the involved penalty factor. In Section 3, we investigate the impact of our scaledmass formulation on discrete spectra for one-and two-dimensional domains. Afterwards, in Section 4, we perform dynamic computations and determine the impact on the convergence behavior. We then extend our formulation to a non-linear application in Section 5, and study the impact of necessary simplifications on the solution quality. In Section 6, we draw conclusions and make suggestions for further research.

Derivation of the mass-scaled semi-discrete formulation
We first develop our proposed finite element formulation with a variationally consistent scaled mass for the example of a linear wave equation. To get comparable results for different simulation cases, we propose a parameter estimation for the scaling factor that takes into account the material parameters, spatial dimension, mesh size and polynomial order.

Variational form of the wave equation
Consider a d-dimensional spatial domain Ω ⊂ R d that is partitioned into an arbitrary number of open subdomains Ω i , i = 1, · · · , N . In each subdomain, the linear wave equation reads: where u(t, x) is the dependent variable, ρ the positive mass per unit volume, and T the positive tensile pre-stress as load per unit surface area. The temporal domain is denoted T =]0, T [ and has end time T . At time t = 0 we require the initial conditions: Additionally, we require at all times and along the entire boundary of Ω, denoted ∂Ω, either Dirichlet (on ∂Ω D ) or Neumann (on ∂Ω N ) conditions: where n denotes the outward facing normal vector. We assume homogeneous Dirichlet conditions purely for ease of notation later on. Finally, the field u(t, x) is subjected to transmission conditions that couple the solution from subdomain to subdomain. We denote by Γ the collection of interfaces between the subdomains Ω i . At each point on Γ, we choose a unit normal vector n as either one of the outward facing unit normal vectors of the neighboring patches. Then, the appropriate transmission conditions across Γ are: where [[·]] denotes the jump operator: Equation (5a) represents the essential kinematic compatibility condition, and Eq. (5b) the natural interfacial equilibrium condition [25]. After multiplying Eq. (2) by a test function, integrating over Ω and performing integration by parts on each subdomain Ω i separately, we find that a weak solution u must satisfy the following statement almost everywhere in time: An appropriate weak formulation then follows after substitution of the natural conditions Eqs. (4b) and (5b), and by searching for u(t, x), which, for almost every t, lies in the Sobolev space H 1 0 (Ω), while its second time derivative is a member of the dual space H −1 (Ω). The space of test functions may also be widened to H 1 0 (Ω), resulting in [26, Ch. 6]: Find u ∈ H 1 T ; H 1 0 (Ω) s.t., for a.e. t, ∀ w ∈ H 1 0 (Ω) : where H 1 T ; H 1 0 (Ω) is the collection of H 1 0 (Ω)-valued functions that satisfy the smoothness requirements of H 1 in the temporal domain T , and where ·, · V ,V denotes a duality pairing between members of the space V and its dual.

Finite element formulation with variationally consistent mass scaling
To obtain a discrete formulation of the above problem, we use the method of lines. This requires a finite element approximation only in space. The usual finite element formulation of Eq. (8) follows from Galerkin's method as: where U h is a discrete subspace of H 1 0 (Ω) as any linear combination of basis functions defined on the mesh: Then, the trial and test functions can be represented as weighted sums of these basis functions, i.e.: and the vector of coefficients for the discrete initial conditions are based on an interpolation or projection of the fields from Eq. (3): We note that the space U h ensures that u h satisfies the essential conditions of Eqs. (4a) and (5a) by construction. We now propose to add variationally consistent terms that follow from the natural conditions of Eqs. (4b) and (5b) to improve the spectral properties of the formulation. We may take the second time derivative of each of these conditions to find: If we consider each element of the mesh a separate subdomain Ω i , then Γ becomes the collection of interior element facets. We then add the following two symmetrically weighted boundary integrals of the conditions in Eqs. (12a) and (12b) to the discrete formulation: with { {·} } the average operator: The additional multiplication with the physical parameters ρ and T in Eq. (13) ensures dimensional consistency. In this form, β needs to have units m 3 , for which we introduce a mesh-size dependency later on. The second term simplifies if ρ and T are continuous across element edges, but the current form is also variationally consistent if the mesh conforms to an internal material interface. If T is a tensor rather than a scalar (as in Section 5) then we define |T | as its largest eigenvalue. For now, we assume a constant scalar valued ρ and T . By adding the terms of Eq. (13) to the finite element formulation of Eq. (9), we obtain the following modified form: After substituting the representation of the test and trial functions per Eq. (10) into Eq. (15), we can rewrite it as a matrix system of equations: As this equation illustrates, the additional terms in the weak formulation only affect the mass matrix. The added matrix M Γ is symmetric positive semidefinite and serves to suppress the highest eigenvalues [21,27]. The variational consistency of the added terms ensures that the true solution also satisfies Eq. (15), which is an important property for ensuring optimal convergence [10,26].

Parameter estimation
The dynamic response of the solution advanced in time by Eq. (16) can be studied by examining the corresponding generalized eigenvalue problem. Of particular importance is the maximum eigenvalue, as this limits the critical time-step size for explicit time-integration algorithms. The generalized eigenvalue problems for the maximum eigenvalues with and without the mass-scaling term read: where λ * max is the target value that we aim to obtain by properly choosing β. We write this target as a factor of the true eigenvalue: λ * max = aλ max . Similarly, we write the eigenvalue of the system without mass scaling asλ max = bλ max . Premultiplying Eq. (17b) by the eigenvector from Eq. (17a) and making use of the symmetry of the matrices then gives: or: To obtain a closed and element-local estimation strategy for β we: 1. Choose an a to eliminate the unknown b from 1 a − 1 b . 2. Approximate λ max in terms of local mesh size h and polynomial order P .

Estimateξ
T K ξ ξ T M Γ ξ based on a simplified problem.

Choice of critical time step increase
First, we choose a, which represents the fraction between the maximum eigenvalue from Eq. (17b) and the true eigenvalue. We base our choice on a target improvement of the critical time step. The critical time step scales inversely to the square root of the eigenvalue. So, for increasing the critical time-step size by a factor ϑ, we require: We choose the following expression for ϑ: For a fixed c > 0, the factor ϑ is always larger than 1 and increases with b. This means that the critical time step will always increase due to the scaled mass, and it will increase more significantly when the discrete formulation without mass scaling severely overpredicts the largest eigenvalues. For a larger c, the time step is increased more, but this may come at the cost of decreased accuracy of the higher frequency response. The choice of Eq. (20) in combination with Eq. (21) also simplifies the expression for β: Based on Eq. (21) and the b-values observed in Fig. 1 (which range between 1.1 and 1.6), a value c = 1 should roughly corresponding to a 50% increase in critical time step, and c = 5 to a 150%-200% increase.

Element local maximum eigenvalue estimate
Next, we approximate λ max by considering the analytical values on simplified domains. On a one-dimensional domain of length L, a two-dimensional square domain with sides L and a three-dimensional cubic domain with sides L the true eigenvalues are: In 2D: λ m,n = T ρ In 3D: λ l,m,n = T ρ When these domains are discretized in the x, y and z directions with n equidistant elements of polynomial order P , then the maximum eigenvalue is the (P n) d -th eigenvalue. For that maximum eigenvalue, Eq. (23) may be written in terms of the element diameter h: In 2D: In 3D: and thus in general:

Estimate of stiffness to scaled mass ratio
Finally, we estimate the ratio betweenξ T Kξ andξ T M Γ ξ. The largest eigenvectors of Eq. (17) conflict most with the continuity requirement of Eq. (5b), and thereby correspond to the largestξ T M Γ ξ [28,21]. Based on this empirical observation, we propose the following For such a function, we can work out both the numerator and denominator as: such that we obtain the approximation: Collecting Eqs. (19), (22), (25) and (29) finally leads to the following estimation for β:

Effect on frequency spectra
Before we study the impact of the mass-scaling term on computations of transient response, we first focus on its effect on the discrete frequency spectra. We consider a onedimensional case and a two-dimensional case with a square domain, for both of which the true eigenvalues are known.

The one-dimensional case: a string
We begin by analyzing the discrete frequency spectrum for the one-dimensional case with Dirichlet boundary conditions on both ends. The eigenfrequencies are the square roots of the eigenvalues: The ordering of the numerical eigenfrequencies is determined based on an iterative scheme: starting with n = 0, we select the numerical eigenfrequency that corresponds to the numerical eigenfunction that yields the lowest L 2 -error compared to the n-th analytical eigenfunction. This numerical eigenfunction is removed from the selection set, and the process is repeated for the subsequent n.
We observe a significant improvement of the spectrum for all polynomial orders. For all choices of c, the eigenfrequencies in the high frequencies range, the last 25% of the spectrum, are suppressed. In many cases, the result is an underprediction of the analytical values, which is not a problem: the high frequencies do generally not contribute to the fidelity of dynamics simulations [3, Ch. 7 and 8], but do limit their critical time step. In the lowest 25% − 50% of the spectrum we observe that the mass scaling does not negatively affect the numerically computed eigenfrequencies, as can be seen clearly in the semi-log inset figures. This is crucial, since these eigenfunction-eigenfrequency pairs typically dominate the system response. For the intermediate range of frequencies, i.e. until roughly 75%, the normalized numerical eigenfrequencies are pushed closer towards the desired value of 1, improving the overall accuracy and smoothness of the discrete spectra.
Since Fig. 3 shows the ratio between the numerical and analytical eigenfrequencies, the results may give the false impression that an increase in c steadily decreases the maximum eigenfrequency. In Fig. 4 we show the non-normalized eigenfrequencies. For P = 1, we do observe that the maximum eigenfrequency is decreased at a steady rate for increasing c, but the same is not true for the cases of P = 2, 3 and 4. Past a threshold value for c the high frequencies are suppressed below the values of an intermediate frequency such that this become the new maximum. The intermediate frequencies only slowly decrease as c becomes larger. This is particularly clear for P = 4, where the eigenfrequencies in the before-last optical branch become dominant.
For higher-order C 0 finite elements, the lack of inter-element continuity is the cause of the different optical branches, and the final branch is due to the lack of first-order continuity (i.e., U h ⊂ C 0 \ C 1 ) [5]. This is the reason for the superior spectral properties of the higherorder continuous spline basis functions used in isogeometric analysis [29]. The mass-scaling term of Eq. (13) acts on this first-order inter-element continuity. It thus specifically targets the highest frequency branch. One potential approach for suppressing the frequencies from the lower optical branches is thus to add additional mass-scaling terms that involve jumps of higher-order derivatives. This concept is explored in an isogeometric analysis framework across patch-interfaces in [21], and on domain boundaries in [20]. However, in the current article we aim to maintain the variational consistency of the finite element formulation, which only permits us to add a penalty on the natural conditions of Eqs. (4b) and (5b), i.e., the lowest derivative.

A two-dimensional case: a square drum
Next, we consider a two-dimensional square domain with either Dirichlet or Neumann conditions around the entire circumference. The Neumann boundary case is particularly important, since the mass-scaling term of Eq. (13) includes a component specifically on the Neumann boundary. Of course, for pure Neumann conditions there is a single zero eigenvalue corresponding to the constant solution eigenfunction. In the following, we remove this eigenvalue from the spectrum.  such that the total number of degrees of freedom is kept at approximately 900. The overall behavior due to the mass scaling is very similar to that of the one-dimensional case, irrespective of the type of boundary condition and the type of elements used. We again observe that the lower end of the spectrum is unaffected, that accuracy of the frequencies in the medium range of the spectrum is significantly improved, and that the frequencies of the medium to higher end of the spectrum are decreased compared to the cases without mass scaling. The consistency between these results and the results for the one-dimensional case for various polynomial orders is remarkable, given the sensitivity to c and the predicted non-trivial dependency on P , d and h in Eq. (30).

Dynamic response: critical time step and convergence
The examples of Section 3 concern tensor-product domains and uniform meshes with regularly shaped elements. To understand if the improvement of the spectrum holds practical significance for more complex cases, we consider a square domain with a polygonal cut-out, and compute the dynamic response of a forced vibration. The improvement of the spectrum should permit a larger time step when explicit time integration is employed. To understand the generality of the impact on the critical time step, we consider two types of meshes with different cut-outs: a mesh with an octagonal cut-out and uniformly distributed elements and a mesh with a star-shaped cut-out with local refinement at the star points. Both meshes are illustrated in Fig. 7 for the refinement case with approximately 400 vertices.
We compute the critical time step with Eq. (1) for multiple refinement-levels and polynomial orders. The results are shown in Figs. 8 and 9 for the octagonal cut-out and the star-shaped cut-out, respectively. The left figures show the values of the critical time steps and the right figures show the factor by which they increase due to mass scaling. The black, blue, green and red lines indicate polynomial orders of 1, 2, 3 and 4, respectively, and the dotted line with the open markers concerns the case without mass scaling (i.e., c = 0). When mass scaling is added, we observe a uniform increase in critical time-step size for both domains, all refinements and all polynomial orders. With a larger value of c, the increase in time-step size is larger. The increase by 50% for c = 1 corresponds well to the prediction from Section 2.3.1. The obtained increase by 100% − 150% for c = 5 is lower than the predicted value and the prediction becomes less accurate for the higher polynomial orders. We anticipate that this overprediction is related to the reordering of the maximum eigenvalues. This was discussed in Section 3.1 and observed in Fig. 4. An increase of c (a) Hexagonal cut-out.
(b) Star-shaped cut-out with local refinement. suppresses the eigenvalues of the highest frequency eigenfunction most, which, at a certain c, causes the eigenvalue corresponding to a lower frequency eigenfunction to dominate. This reduces the effectiveness of the mass-scaling term. As a result, the approximation strategy of Section 2.3.1, which may be considered a first-order approximation around c = 0, becomes an overprediction.      Next, we compute the dynamic response for a system with the following true solution: from which we determine the initial condition, as well as the Neumann data g andg on the circumference of the cut-out. This true solution has a period of 2 seconds.
For the convergence study we focus on the case of the octagonal cut-out, as the refinement levels involve uniform refinement. We make use of the standard fourth-order explicit Runge-Kutta algorithm for the time integration. Figure 10a shows the convergence of the L 2 -error after T = 0.1 seconds and Fig. 10b shows the results precisely five periods later, at T = 10.1 seconds. The previously plotted critical time steps are used as time steps in the Runge-Kutta algorithm. This means that the data points for c = 1 and c = 5 required approximately 33% and 60% fewer time step computations than those without mass scaling. In Fig. 10a we observe almost identical errors and optimal rates of convergence for all cases, indicating that the proposed mass scaling does not negatively affect the solution quality. The convergence curves plotted in Fig. 10b are less clean. Nevertheless, the mass scaling never negatively affects the errors, despite the increase in time-step size. The exceptions are the last two datapoints for the case of P = 3 and c = 5, which could be due to the error in the time discretization dominating the complete solution error at this level of spatial resolution.

Extension towards non-linear solid mechanics
To illustrate the use of our method for a more complex scenario, we consider a beam oscillating at high amplitudes compared to its length. The non-negligible rotational deformation invalidates the linear strain assumption. As a consequence, a geometrically non-linear model needs to be considered. In the domain Ω of the reference configuration, the variational statement governing free vibration reads [25,30]: where u is the displacement vector, ρ 0 the density in the reference configuration and P (u) the first Piola-Kirchoff stress tensor. The first and second Piola-Kirchoff stress tensors relate to one-another via the deformation gradient: As a constitutive law we make use of the Saint Venant-Kirchhoff hyperelastic relation: with C the fourth order material response tensor, λ and µ the Lamé parameters and E(u) the Green-Lagrange strain tensor: All gradients in Eqs. (33), (34) and (36) operate in the reference configuration. Under the assumption that the solution is smooth enough we may perform integration by parts on the second term to reveal the implicitly assumed transmission conditions between two neighboring domains: so that it is the continuity of the first Piola-Kirchoff stress that represents the interfacial equilibrium condition when observed in the reference configuration [25, p. 204]. The theory outlined in Section 2 thus suggests to use a discrete formulation of Eq. (33) that also includes a mass-scaling term based on the second time derivatives of Eq. (37). However, due to its non-linearity, addition of a penalty of this traction jump to the mass matrix would require re-assembly and re-computation of the LU-decomposition of the stabilized mass matrix in each time step. The required additional computational expense to perform these operations is of same order of magnitude as the expense of a single implicit time-integration step, and this thus defeats the point of using an explicit time-integration method. Instead, we propose to add a linearized version of this penalty. Then, the extension of Eq. (15) to the current problem would produce: where we have assumed homogeneous Dirichlet and Neumann conditions. For |C| we take its maximum eigenvalue, that is max(µ, K) with K the bulk modulus as K = λ + 2 3 µ. For β, we make use of a straightforward extension of the expression from Eq. (30) from the scalar wave equation to a vector wave equation: The additional division by d 2 is due to each node carrying d-degrees of freedom, such that the index of the maximum eigenvalue as used in Section 2.3.2 now becomes (dP n) d .
The new terms that are added to the mass matrix involve the linear strain: This raises questions regarding the variational consistency of the formulation. The linearized strain measure yields erroneous (i.e., nonzero) strain values for rigid body rotation. However, this is not necessarily problematic: we add a penalty to the jump of this linearized stress, rather than the linearized stress itself. When the material parameters are continuous, and the true solution is smooth enough, Eq. (38) is still satisfied by the analytical solution. We investigate the impact of this model simplification by focusing on an example that involves non-negligible rotations. Consider the asymmetric beam illustrated in Fig. 11a. The material parameters are chosen as ρ = 2700kg/m 3 , E = 73GPa/m 2 and ν = 0.33, and correspond to an aluminum alloy. The beam is clamped at x = 0, which is the left-most face in the illustrations of Fig. 11. At t = 0, the beam is released from the following stationary initial displacement: for R = 0.15. This describes a circular shape in the y = 0 plane with a radius R with no stretching in the initial z = 0 plane. This initial displacement is shown in Fig. 11b. Figures 11c and 11d illustrate the coarse and fine meshes used in the following computations. First, we study the dependency of the critical time step on c, the remaining free parameter in the estimate of β in Eq. (39). To determine the eigenvalues for this nonlinear problem, we consider the stiffness matrix that arises from the linearization by a small perturbation around the initial displacement: (42b) Figure 12 then shows the ratio between the critical time step and the critical time step for the case without mass scaling (i.e., c = 0) as a function of c. These are computed for the coarse mesh. Overall, we observe significantly higher values for the increase in critical time step for this benchmark problem than for the benchmark problems of Sections 3 and 4. The ratio keeps increasing when linear basis functions are used, whereas it levels off for the quadratic basis functions. This corresponds to the observation in Section 3.1 that, for higher order polynomials, the suppression of the highest eigenvalues eventually leads to one of the lower eigenvalues dominating.
To study whether there is any quality loss when we make use of a c that produces a significant cost-reduction, we choose more extreme values for c than those in Section 4. For the computation with the linear basis functions we make use of c = 40 and c = 5000, which yield factors of increase of critical time step of approximately 6.5 and 22.2, respectively. For the computation with quadratic basis functions we use c = 10, as this already roughly produces the plateau value of a factor 4.0 increase in critical time step. Figure 13 shows the (a) Beam geometry.
(b) Initial displacement.   computed tip displacement as a function of time for the linear and quadratic cases computed on the coarse mesh of Fig. 11c. In both figures, the light-gray line is the reference solution computed with linear basis functions on the refined mesh that is depicted in Fig. 11d.
The results for both polynomial orders demonstrate that the impact of mass scaling with moderate choices of c is minor compared to the source of error originating from the insufficient mesh resolution. In both cases we even observe small improvements. In Fig. 13a, the case of c = 0 exhibits a severe overprediction of the dominant eigenfrequency. For c = 20, this overprediction is reduced slightly while requiring 6.5 times less time step computations. The case c = 5000 is chosen as an extreme as it shows that this behavior persists, and that the mass scaling can be tuned to yield the correct dominant frequency. Clearly, this comes at the cost of an apparent artificial damping, which reduces the accuracy of the high-and low-frequency response.
The quality of the coarse-mesh simulation is improved by the increase in polynomial order, as shown in Fig. 13b. For that simulation, the case of c = 5 does not noticeably affect the dominant frequency within these first 5.5 periods. However, it does appear to suppress the artificial high-frequency oscillations that are visible in the dotted line in the zoom-in, and which are not present in the result of the reference computation.
Due to the geometrical asymmetry, the beam response includes a component out of the y-plane. We compute the twist at the tip of the beam as a measure for this deformation. Figure 14 shows the frequency response thereof, as computed with a fast Fourier transform of the data from the first 0.01 seconds, windowed with a Blackman window [31]. Figure 14a shows the results for the computation with linear basis functions on the coarse mesh. The c = 5000 case has been removed from the plot to avoid clutter and since it produces a largely inaccurate spectrum due to the artificial damping. Figure 14b shows the results for the computation with quadratic basis functions on the coarse mesh, which again results in a considerable improvement in the reproduction of the reference spectrum. For both computations, we observe that the scaled mass only marginally affects the system response. The locations and magnitudes of the peaks are not affected, and the changes in magnitudes are minor. For the linear basis functions, the overpredicted magnitude of one of the dominant peaks is even reduced. For the quadratic basis functions, the two spectra only become distinguishable past the 6000Hz frequency, at which range there are no more significant peaks.

Conclusion and outlook
In this article, we have proposed a variationally consistent addition to the mass matrix to suppress the highest eigenvalues in structural dynamics related applications. Our proposed term follows from a symmetric penalty of the natural transmission condition across element interfaces. For typical solid mechanics applications, the natural transmission condition is the interfacial equilibrium condition, which says that the jump of the traction should be zero. We weigh the new term by a small penalty, for which we derived an expression that ensures consistent results irrespective of the mesh size, polynomial order and spatial dimension.
For the linear wave equation, we have shown that the mass scaling improves the spectra on simple domains considerably, for all considered polynomial orders (up to quartics). In all cases, the low frequencies are virtually unaffected, the intermediate frequencies become more accurate, and the high frequencies are reduced. Also on more complicated domains and on irregular meshes, the mass-scaling term permits increases of critical time steps by factors of 2.5 without negatively affecting solution accuracy or convergence behavior.
We have also proposed a simplification that makes the mass-scaling term suitable for explicit time-stepping methods in non-linear solid mechanics applications. In our example problem, the mass scaling permits a factor 6.5 and factor 4.0 increase in the critical timestep size for linear and quadratic basis functions, respectively, without any negative impact on the solution accuracy (arguably, we even observe improvements).
There are various potential directions for future research. The proposed mass-scaling term may be extended to target the higher-order derivatives across element interfaces for higher-order polynomial basis functions, i.e., an application of the formulation from [21] for C 0 -continuous basis functions. This would likely mitigate the plateauing behavior observed in Section 3.1 and Fig. 12. Secondly, one could investigate operator splitting techniques to move the mass-scaling term to the right-hand-side. This would guarantee the variational consistency of the formulation for a wide range of non-linear problems without the need for reassembly of the mass matrix in each time step. Additionally, operator splitting could provide a pathway for using our mass-scaling term in combination with a lumped mass matrix.