Lumped mass acoustic and membrane eigenanalysis using the global collocation method

The paper proposes a direct way to build lumped masses for performing eigenvalue analysis using the global collocation method in conjunction with tensor product Lagrange polynomials. Although the computational mesh is structured, it has a non-uniform density, in such a way that the internal nodes are located at the position of Gaussian points or the images of the roots of Chebyshev polynomials of second kind. As a result, the mass matrix degenerates to the identity matrix. In this particular nodal collocation procedure, no complex eigenvalue appears. The theory is successfully applied to rectangular and circular acoustic cavities and membranes. Subjects: Acoustical Engineering, Computer Science, Design, Engineering & Technology, Mathematics & Statistics for Engineers, Mechanical Engineering, Technology


Introduction
Transfinite elements were inspired in early 1970s (Gordon & Hall, 1973) for the purposes of CAD/CAE integration, but only much later were applied for the numerical solution of static and dynamic engineering problems using natural cubic B-splines (Kanarachos & Deriziotis, 1989;Kanarachos,

PUBLIC INTEREST STATEMENT
Physical phenomena are usually of timedependent character. Therefore, adequate computer effort is required to simulate their progress. One of the most important information needed to be extracted by an analyst (or a designer) is the resonance frequency. This frequency must be avoided because it is harmful or extremely annoying. This paper proposes a fast method for the numerical extraction of the resonant frequencies, particularly those of acoustic cavities (such as rooms, theaters, music halls, and musical instruments) and those of vibrating elastic membranes. The method eventually bypasses the computation of the mass matrix, and uses only the stiffness matrix. In this way, the total computation (CPU) time is highly reduced. In terms of mathematics, the proposed procedure is a global collocation method based on tensor product Lagrange polynomials. More generally, the same idea can be applied to the numerical solution of wave propagation problems such as acoustics and elastodynamics. Provatidis, Deriziotis, & Foteas, 1999;Provatidis & Kanarachos, 2001), as well as piecewise linear approximation (Provatidis, 2004) and Lagrange polynomials (Dimitriou, 2004;Provatidis, 2005Provatidis, , 2006aProvatidis, , 2006bProvatidis, , 2006cProvatidis, , 2009aProvatidis, , 2012, among others. However, due to the fact that the involved matrices are fully populated, the overall computational effort is not substantially less than that of the conventional finite element method for the same mesh density. This shortcoming was a source of inspiration to propose the substitution of the Galerkin-Ritz formulation by a global collocation method, in such a manner that the global character of the approximation is preserved (Provatidis, 2006b, p. 6704). At about the same time, global radial basis collocation schemes were proposed (Kee, Liu, & Lu, 2007), among other papers of the same team. Later, similar ideas were proposed in the context of isogeometric analysis (Auricchio, Da Veiga, Hughes, Reali, & Sangalli, 2010). For a review paper on CAD/CAE integration, which starts with Coons-Gordon interpolation and continues with the chronologically subsequent interpolations (i.e. Bézier, B-splines, and NURBS), the reader is referred to a recent paper (Provatidis, 2013).
The first attempts of using the global collocation methods in elliptic problems under pure Dirichlet-type (Provatidis, 2008a), as well as partially Dirichlet-type and partially Neumann-type (Provatidis, 2009b) boundary conditions, were encouraging. In problems of first (Dirichlet) type, no difficulty appears (Provatidis, 2008a). In contrast, when attempting to solve problems of second (Neumann) type in a domain owing some corners (surrounded by Neumann b.c.'s), it becomes necessary to treat them in a way somehow similar with that of the boundary element methods (see, e.g. Brebbia & Dominguez, 1992). In more details, there are two normal vectors at any corner, a fact that does not allow for the direct elimination of "slave" degrees of freedom (along the Neumann part of the boundary: Γ 2 ) since the involved matrices become non-square.
The second major problem is that when applying the so-called "nodal collocation" technique, according to which the collocation points coincide with the internal nodes (thus, no actual computation of the mass matrix is required), some of the obtained numerical eigenvalues are found to be complex. So far, this numerical shortcoming has been resolved by taking the collocation points at the position of Gaussian points (roots of Legendre polynomials) or at the images of the roots of Chebyshev polynomials (particularly of second kind) (Provatidis, 2008b); these choices constitute what is usually called "orthogonal collocation." Within this context, this paper shows first that in case of a structured mesh, the calculated eigenvalues do not depend on the particular location of a given number of internal nodes, provided a constant position of the collocation points. Second, it shows that when positioning the internal nodes at the location of either Gaussian points or the images of the roots of Chebyshev polynomials, no mass matrix is required to be calculated (it merely coincides with the identity matrix). In other words, for the first time, this paper achieves to successfully apply the nodal collocation in a way that avoids the (unrealistic) complex eigenvalues. Not only that but as previously explained, these eigenvalues are identical with those that are obtained using a uniform mesh in conjunction with collocation points at the aforementioned roots of the two polynomial types. Third, the theory is sustained by four examples in simple rectangular and circular shapes, which concern either acoustic cavities or membranes.

The governing equations
The governing equations of motion are as follows.
(1) Acoustic wave: where u is the acoustic pressure (or velocity potential) and c is the wave speed.
(2) Wave in membranes: where T is the in-plane tensile force, w is the deflection, p is the pressure, and ρ is the mass density per unit area of the membrane. Equation 2 implies that the wave speed is given by c = √ (T∕ ).
Due to the perfect analogy between Equations 1 and 2, henceforth we will refer mostly to Equation 1.

Matrix formulation
We consider a two-dimensional acoustic cavity (or a membrane) which is discretized in n x × n y (structured) intervals, i.e. n = (n x + 1)×(n y + 1) nodes. From these nodes, n b = 2(n x + n y ) nodes belong to the boundary of the cavity (or the membrane) and n i = (n x − 1)×(n y − 1) to the interior of the domain. Among the aforementioned n b nodes, n 1 are under Dirichlet boundary conditions (given acoustic pressure or fixed flexure, such as u = 0), while the rest n 2 are under Neumann boundary conditions (zero normal velocity or traction-free condition, implying that ⃗ n ⋅ ∇u = 0, where ⃗ n is the unit normal vector to the boundary).
Obviously, if we ignore the possible double nodes, it holds that: As usual, the approximate solution of Equation 1 is given by in which every global shape function N j (x, y) is a tensor product of two Lagrange polynomials.
Temporal and spatial derivatives of Equation 3 lead to the following two expressions: Substituting Equation 4 into Equation 1, the latter becomes: Although Equation 5 could be, in principle, applied to all the above-mentioned n i internal points (nodal collocation), previous numerical experimentation on the solution of one-dimensional problems using a uniform mesh has revealed the appearance of fictitious complex eigenvalues (Provatidis, 2008b). Therefore, one conservative choice is to locate the collocation points at another set of n i points such as the Gauss points or the roots of Chebyshev polynomials (first or second kind). (1) In any case, the collocation of Equation 5 at the above-mentioned n i points leads to the following matrix system: where the elements of the matrices [M] and [K] are given, respectively, by:

Boundary conditions
The boundary conditions on the boundary, Γ 1 +Γ 2 = Γ, are of two types: (a) Essential conditions, such as u =ū on Γ 1 (Dirichlet type) The fact that the matrices involved in Equation 6 are non-square requires further processing as follows. The vector of the acoustic pressure, u, splits into three groups with the first (u 1 ) referring to n 1 nodes (owing boundary conditions of Dirichlet type), the second (u 2 ) to n 2 nodes (owing boundary conditions of Neumann type), while the third (u i ) to the internal degrees of freedom. If the terms of matrices [M] and [K] are rearranged so that the lower index "1" corresponds to the n 1 nodes, the index "2" corresponds to the n 2 nodes, while the index "i" corresponds to the n i nodes, then Equation 6 becomes: For the needs of eigenvalue problem, we can accept that: and therefore Equation 8 takes the form: Concerning Equation 9, it is noted that the acoustic pressure u equals to zero when the part Γ 1 along the boundary of the acoustic cavity is open to the outer space. Alternatively, since the governing equation of a vibrating membrane is a two-dimensional form of the wave equation [cf. Equation 2], the condition u 1 = 0 (i.e. the deflection w = 0) means that the part Γ 1 of the boundary is fixed.
Unlike the procedure in finite elements method, we are not allowed to solve the eigenvalue problem by using all the terms involved in Equation 10. This is because, in essence, we have not taken into account the boundary conditions of Neumann type, since the n 2 partial derivatives equal to zero and they import an equal number of linear dependencies between the (n 2 + n i ) degrees of freedom, so that it can eventually arise only n i degrees of freedom (Provatidis, 2009b).
The process of elimination is as follows. Equation 3 becomes: where, again, ⃗ n represents the unit vector perpendicular to the boundary (its discontinuity at corners justifies the use of double nodes). Applying the zero value of Equation 11 at all the n 2 boundary nodes of Neumann type (free-free boundary conditions), we derive the following equation: Since u 1 = 0, Equation 12 is simplified as follows: where each element of the matrix appearing on the left-hand side of Equation 13 is given by:

Basis change
For a given discretization using n x × n y intervals, it will be shown that no matter where exactly the position of the n = (n x + 1)×(n y + 1) nodes is, provided the mesh is structured, but not necessarily at equal distances.
Actually, along each line normal to the sides in the reference square, the number of nodes does not change when these are all "uniformly" shifted and are finally arranged along ( = const. or = const.) directions, as shown in Figure 1. As a result, the polynomial degree p does not change but always remains equal to the number of subdivisions (p x = n x and p y = n y ).
The above statements are systematically justified, starting with a one-dimensional analog, as follows.
Proof Based on the above transformations, i.e. Equations 19 and 20, obviously we can successively write: The latter completes the proof.
An application example for a 1D analog is given in Appendix A.
In 2D problems, which is the topic of this paper, since the shape functions are tensor products of Lagrange polynomials, i.e. of the form N J (x, y) = L i (x)L j (y), the extension of the above-mentioned theorems is straightforward.

The proposed lumped mass approach
Since the calculated eigenvalues depend only on the position of the collocation points (Section 3.1), and not on the particular position of the nodal points, it is much effective to put the nodes at the collocation points that have been previously used for the successful numerical solution using a uniform mesh.
In this way, the mass matrix degenerates to the form: where I is the identity matrix. Therefore, no computation effort is needed, because the eigenvalues of the dynamic problems equal to the eigenvalues of the stiffness matrix K only, which of course are afterward multiplied by the factor (c 2 ).

Numerical examples
The theory is sustained by four examples, the first two referring to a rectangular domain, whereas the rest two to a circular domain. For every shape, both Dirichlet-and Neumann-type boundary conditions are tested.

Example 1: fixed rectangular membrane
Let us consider a rectangular membrane of dimensions a × b = 2.5 × 1.1 m, which is fixed along its entire boundary Γ (Dirichlet boundary conditions: u = 0).
The exact eigenvalues are given by the formula (Courant & Hilbert, 1966): and the corresponding eigenmodes v(x, y), with c denoting the corresponding wave velocity (e.g. c = 1 m/s).

K ���
mn 2 = 2 c 2 (m∕a) 2 + n∕b 2 , m, n = 1, 2, … v(x, y) = sin (m x∕a ⋅ sin (n y∕b), m, n = 0, 1, 2, … The longest side of the rectangular was discretized into n x = 4, 6, 8, …, 22 segments, whereas the shortest into n y = n x /2 ones. Two specific sets of collocation points were tested. The first concerns Gaussian points and the second Chebyshev of second kind. Both the "uniform" (mapping like that in Figure 1(a) and the "distorted" (mapping like that in Figure 1(b) meshes lead to identical results, which are shown in Figures 2 and 3, using Gaussian points and roots of Chebyshev of second kind, respectively.
In the sequel, we increase the number of segments from n x = 22 to n x = 50, with a step of "2". In addition, Figure 4 provides comparisons of CPU time with and without using the technique of lumped mass. It is observed that the larger the polynomial degree p x , the faster the lumped mass formulation is. For n x = 50, the proposed lumped mass (non-uniform mesh) is 3.8 times faster than the usual orthogonal (uniform mesh) formulation.

Example 2: rectangular acoustic cavity with hard walls
Concerning the Neumann problem, which corresponds to an acoustic cavity with hard walls, the exact eigenvalues are now given by the formula (Courant & Hilbert, 1966): and the corresponding eigenmodes v(x, y), As previously, Figure 5 shows that the numerical solution rapidly converges toward to exact analytical solution (Equation 23). Again, identical results were obtained with and without the lumped mass.

Example 3: fixed circular membrane
Consider a circular membrane of radius a = 1 m and wave velocity c = 1 m/s, which is fixed along its outer circumference. The theoretical eigenvalues are given by the formula (Courant & Hilbert, 1966): and the corresponding eigenmodes u(x, y) are: where J � m (ka) is the first derivative of the Bessel function J m (ka) of the first kind and order m and k = ∕c the wavenumber.
The two meshes, one "uniform" and another "distorted" are shown in Figure 6. Again, the global collocation solution based on the uniform mesh in conjunction with collocation points at the Gaussian points (orthogonal collocation) is identical with the nodal collocation method applied to the "distorted" mesh. The results are shown in Table 1.

Mode
Exact eigenvalues 2 (s −2 ) Errors (in %) of the calculated eigenvalues n x = n y = 6 n x = n y = 8 n x = n y = 10 n x = n y = 12

Example 4: circular acoustic cavity with hard walls
In case of a perfectly hard wall, the theoretical eigenvalues are given by the formula (Courant & Hilbert, 1966, p. 304): and the corresponding eigenmodes u(x, y), written in coordinates (r, ), by: where J � m (ka) is the first derivative of the Bessel function J m (ka) of the first kind and order m and k = ∕c the wavenumber. The quality of calculated eigenvalues is shown in Table 2 in the form of errors (in %) with respect to the exact analytical solutions (Equation 25). Once again, one may observe an excellent convergence when increasing the number of subdivisions n x (=n y ) per side.

Discussion and conclusions
In all examples of this study, the proposed global collocation method-based on tensor product Lagrange polynomials-was found to perform extremely well and rapidly converged to the exact solution. Identical results were found with and without using the technique of lumped masses, but in the second case, the CPU time was reduced by a factor that increases with the polynomial degree (mesh size).
The imposition of purely Dirichlet boundary conditions was a straight-forward procedure, whereas the imposition of purely Neumann conditions required a careful handling of the corners.
In general, the images of roots of Legendre polynomials (Gaussian points) lead to slightly better results than those obtained by the roots of Chebyshev polynomials of second kind.
Taking into consideration that the domains used in this study were of simple shape (rectangular or circle), we have to mention that for the analysis of components or domains of complicated shape, it becomes necessary to split them into a certain number of large convex subdomains. Then, these subregions can be coupled using compatibility and equilibrium conditions. The extension of the global collocation method, in both its usual and the proposed lumped mass versions, to transient scalar wave as well as to transient and eigenvalue elastodynamics is a straight-forward procedure. It is anticipated that its extension to 3D problems will perform equally well.

Mode
Exact eigenvalues 2 (s −2 ) Errors (in %) of the calculated eigenvalues n x = n y = 6 n x = n y = 8 n x = n y = 10 n x = n y = 12