A Novel Highly Accurate Finite-Element Family

Copyright © 2017 Giovanni Bernardini et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. A novelNth order finite element for interior acoustics and structural dynamics is presented, withN arbitrarily large. The element is based upon a three-dimensional extension of the Coons patch technique, which combines high-order Lagrange and Hermite interpolation schemes. Numerical applications are presented, which include the evaluation of the natural frequencies and modes of vibration of (1) air inside a cavity (interior acoustics) and (2) finite-thickness beams and plates (structural dynamics). The numerical results presented are assessed through a comparison with analytical and numerical results.They show that the proposed methodology is highly accurate. The main advantages however are (1) its flexibility in obtaining different level of accuracy (pconvergence) simply by increasing the number of nodes, as one would do for h-convergence, (2) the applicability to arbitrarily complex configurations, and (3) the ability to treat beamand shell-like structures as three-dimensional small-thickness elements.


Introduction
Interior cabin noise is a challenging problem in most aircraft design and has received considerable attention in the last decade.This is particularly true for rotary wing aircraft because the propulsive system induces direct acoustic disturbances and fuselage vibrations that, in turn, may cause unacceptable ride discomfort inside the cabin area hosting passengers, as well as significant impact on the fatiguelife of the structures and hence on the maintenance costs [1][2][3].One of the main characteristics of aircraft cabin noise is the wide range of frequencies of interest, due to different noise sources: fuselage boundary-layer, airborne, and structure-borne noise are among the most important ones [4].Boundary-layer noise, generated by the shaking of the fuselage-wall due to external turbulence pressure fluctuations, is a random, broadband, high-frequency signal.Airborne noise is a mid/low-frequency tonal noise associated with periodic pressure fluctuations from the propulsive system impinging the fuselage structure that, in turn, excites interior acoustics, whereas structure-borne noise is related to the acoustic energy generated by periodic vibratory loads (rotor hub-loads, engine vibrations, gearbox, etc.) acting on the airframe.
From the above considerations, one may infer that interior noise prediction is a challenging problem that requires the development of efficient formulations able to model accurately the interactions between the fuselage structural dynamics and the cabin interior acoustics, even in the presence of wide frequency band signals.An approach that may be used to address such a problem consists of coupling Finite-Element Methods (FEM) for structural dynamics with Boundary Element Methods (BEM) for interior acoustics [5].However, in several practical applications both standard FEM and BEM methods become inaccurate and computationally intensive when the number of elements needed to model the problem becomes large [6].In fact, the accuracy of these methods strictly depends on both the number of nodes used in the model and the order of the shape functions used to interpolate the solution within each element: a rule of thumb to maintain accuracy in both methods is to have six linear or three parabolic elements per wavelength [6].Moreover, in several practical applications FEM are more efficient than BEM for the numerical solution of interior acoustic problems [7].
In the last few years, considerable work has been done to improve the prediction capabilities of FEM algorithms in analyzing structural dynamics/interior acoustics problems 2 International Journal of Aerospace Engineering characterized by midfrequency signals [8,9].Among these, we can mention the Galerkin least-squares finite-element method for the solution of the two-dimensional Helmholtz equation [10], the Galerkin residual-free bubbles method [11,12], the smoothed FEM using cubic spline polynomial functions in hexahedral elements [13], and the isogeometric FEM [14].During the years, very accurate FEM methods suited for structural dynamics and fluid dynamics applications have been developed; among these, it is worth mentioning the Spectral Finite-Element Method (SFEM), based on Lagrange polynomials on the Gauss-Lobatto-Legendre grid [15], the Fourier transform-based and Wavelet transform-based spectral FEM [16,17], the hierarchical -FEM [18][19][20], and  1 -FEM methods based upon isoparametric Hermite elements [21,22].Since its origin, SFEM have been successfully applied in several fields of the physics and applied sciences: acoustics, fluid dynamics, heat transfer, and structural dynamics.In mechanical/aeronautical engineering they are widely used in wave propagation, in both homogeneous and inhomogeneous structures [15,23], in structural health monitoring applications [24], as well as in structural dynamics of rotating beams [25].
In this paper, we concentrate on finite elements.This methodology is widely used in engineering and science applications [26,27] and is extensively analyzed from both practical and theoretical points of view [20,[28][29][30].Specifically, we propose a new efficient and accurate finite-element methodology, suited for the analysis of both structural dynamics and interior acoustics.The formulation proposed here is a user-friendly evolution of the methodology developed in the past by the authors and their collaborators.This was presented in [31][32][33][34][35][36][37][38], with increasing level of algorithmic sophistication.This work is to be understood as a first step towards the development of a highly accurate finiteelement formulation for the evaluation of the natural modes of vibration of the air inside a cavity (interior acoustics) and/or of an elastic structure (structural dynamics), for relatively high spatial frequencies (specifically higher than those that may be efficiently obtained with the methodologies presently available), so as to make the technique useful for structural acoustics applications, which involve the coupling of structure and air.Specifically, the finite-element methodology presented here is based upon a combination of two important techniques: (1) the three-dimensional extension of the Coons patch technique [39][40][41] and (2) the high-order Lagrange and Hermite interpolation schemes [42][43][44].This combination is very powerful and yields the distinguishing feature of combining high efficiency (with the possibility of capturing relatively high spatial frequencies, as required in aeroacoustics applications), with user-friendliness, giving a finite element that is very accurate and computationally efficient.More interestingly, it has the unique feature of being flexible, in the sense that one may increase the order of the scheme accuracy (-convergence) just by changing the number of nodes (as one would for ℎ-convergence).Another important feature is its applicability to arbitrarily complex configurations, using always the same type of element.In particular, the formulation covers both interior acoustics and structures and is able to model beams and plates, which are to be treated as three-dimensional small-thickness structures.The above characteristics are important in that the technique proposed here is envisioned within the context of fully automated multidisciplinary optimal design [45].In addition the following two features (not all simultaneously present in other finite elements currently available for aeroacoustics) appear important: (1) the element captures modes with relatively high spatial frequencies (as essential for aeroacoustics utilization) and does so efficiently, so as to give good results with relatively few elements, an important feature when repeated calculations occur, as in automated optimization, and (2) the element is user-friendly.

Preliminaries: Lagrange and Hermite Interpolation
The formulation presented here involves a judicious combination of high-order Lagrange and Hermite interpolation polynomials [42][43][44], which are briefly reviewed here.First, we address the two-point interpolation polynomials (needed in Sections 3-5), and then we consider the -point interpolation, with  arbitrarily large (used in Section 6).
The Two-Point Lagrange and Hermite Interpolation Polynomials.Here, we discuss the two-point Lagrange and Hermite interpolation polynomials [42][43][44].Consider a function  = (), defined over the interval [0, 2].Let us divide the interval [0, 2] into  subintervals, for which we introduce a local coordinate , so as to have that the end points of each subinterval are given by  = ±1.Within each subinterval, the Lagrange linear interpolation is given by where  ± = (±1) denotes the values of () at the end points, whereas the Lagrange first-order interpolation polynomials  ± () are given by The resulting interpolated function over the whole interval [0, 2] is continuous (class C 0 ) and piecewise linear.It will be referred to as the one-dimensional first-order Lagrange interpolation.
On the other hand, within each subinterval, the standard Hermite interpolation is where  ± denotes the values of () at  = ±1 and   ± denotes the values   () = / at  = ±1, whereas the Hermite interpolation polynomials  ± () and  ± () are given by International Journal of Aerospace Engineering 3 The resulting interpolated function in [0, 2] is continuous with its first derivative (class C 1 ) and is piecewise cubic.It will be referred to as the one-dimensional third-order Hermite interpolation.
The -Node Lagrange and Hermite Interpolation Polynomials.
Let us consider first the -node Lagrange interpolation; namely, where is the total number of subintervals between the  nodes   ( = 0, . . ., ), and [Throughout the paper  0 and   coincide with the end points of the interval under consideration.]In particular, if  = 2,  ∈ [−1, 1], and  ± = ±1, we have  ± () = (1/2)(1 ± ), in agreement with (2).Note that the Lagrange interpolation polynomials   () satisfy the standard interpolation condition   (  ) =   (,  = 0, . . ., ) . ( It is well known that, for  > 7, the Lagrange interpolation polynomials are unstable if   are uniformly spaced [44].However, the instability disappears if   coincide with the Gauss quadrature abscissas.Unfortunately, the Gauss quadrature abscissas do not include the end points of the interval.The Gauss-Lobatto quadrature scheme on the other hand includes these points and is not affected by the above instability issue.Accordingly, in the case of highaccuracy schemes,   used in this paper coincide with the Gauss-Lobatto abscissas. Next, let us consider the -node extension of the twopoint Hermite interpolation presented in (3) and (4).There exist two possible approaches.The first is to increase the order of the derivatives at the end points.The other consists in using as interpolation parameters the function and its first derivative not only at the end points but also at additional interior points.As shown in [36], the latter is the most convenient (and definitely the most user-friendly).Accordingly, here we concentrate on such an approach.
In analogy with the third-order formulation (see ( 3)), we have where the polynomials   () and   () are given by [43] with   obtained by imposing These polynomials satisfy the Hermite interpolation conditions Indeed, it is apparent that the polynomials   () and   () vanish with their first derivatives at  =   , for all  ̸ =  (see (8)).Thus, we only have to verify what happens when  = .It is apparent that   (  ) = 1.In addition, taking the logarithmic derivative of   (), we have which yields,    (  ) = 0 (see (11)).Moreover, it is apparent that   (  ) = 0.In addition, we have which yields    (  ) = 1 (see (8)).In summary,   () and   () satisfy all the conditions in (12).
The polynomials have degree 2 + 1 = 2 − 1.We will refer to this as the interpolation of order 2 + 1.

Motivation: Hermite Brick and BVD Problem
In order to set this work in the proper context, we present some details of the formulation for the so-called Hermite brick.This allows us to motivate the introduction of the family of high-order elements proposed here.
Third-Order Hermite Brick.In order to place the finiteelement family proposed here in the proper perspective, let us consider a drawback of the Hermite interpolation.To this end, assume that the problem is defined in a topologically hexahedral region.Let us subdivide the region into topologically hexahedral subregions, here referred to as the brick, which are described by the mapping where   ∈ [−1, 1] ( = 1, 2, 3) are curvilinear coordinates, whereas x = x(  ) is a suitably smooth function.[Sometimes it is more convenient to use the symbols   ( = 1, 2, 3); other times we prefer to use (, , ).Accordingly, we use ( 1 ,  2 ,  3 ) and (, , ) interchangeably.] In this section, we use the third-order Hermite interpolation in all directions (for both geometry and unknown function).
This yields the following interpolation for the unknown (  ); namely, The symbol ∑ s , with s = ( 1 ,  2 ,  3 ), where   stands for either + or −, is understood as a sum that spans over the eight values of s, which correspond to the eight vertices of the brick.Furthermore, the symbol  () denotes the partial derivative with respect to   :  () = /  .Similarly, we have  () =  2 /    and  (123) =  3 / (1)  s s s In summary, for the three-dimensional third-order Hermite interpolation scheme, the finite-element unknowns are the nodal values of (1) the unknown function, (2) its three first-order partial derivatives, (3) its three second-order mixed derivatives, and (4) its third-order mixed derivative  3 / 1  2  3 .This yields a total of eight unknowns per node, out of which only one is the nodal value of the function, the other seven being related to the various derivatives.The finite element described above will be referred to as the thirdorder Hermite brick.
The expressions in (17) provide the local finite-element shape functions.They may be assembled to yield a global finite-element interpolation over the whole block, of the type where the unknowns   comprise the values of  and of all the partial derivatives mentioned above, evaluated at the  = ( + 1) 3 nodes.On the other hand,   (  ) are the global shape functions [20], obtained by assembling the local shape functions in (17).
The Base Vector Discontinuity (BVD) Problem.The thirdorder Hermite brick is seldom used, because of a problem that may arise at the interface of two or more bricks.Specifically, the issue arises when the coordinate lines of two adjacent bricks present a discontinuity, namely, when the covariant base vectors which are tangent to the coordinate lines, are discontinuous (either in magnitude or direction).In this case, the partial derivatives of a function (  ) with respect to   are given by and are discontinuous.As far as the first-order derivatives are concerned, the problem is removed by assuming as unknowns the values of the Cartesian coordinates of grad .
The first-order partial derivatives may then be obtained using (20).The problem, however, remains for the second-order derivatives, because, in order to express them in terms of Cartesian components, one needs the complete Hessian matrix, namely, all the second-order derivatives, and not only the mixed ones, which are the only ones utilized in the threedimensional Hermite interpolation.Similar considerations hold for the third-order derivatives.The problem discussed above will be referred to as the BVD problem (Base Vector Discontinuity problem).[It may be noted that the Hermite scheme is not limited to the third order.Higher order schemes may be easily introduced.However, these extensions are also subject to the BVD problem, and hence of little interest here.] In order to remedy the problem, various unsatisfactory attempts have been considered by the authors and their collaborators and are summarized in [31,32].The problem was resolved in [33][34][35][36][37][38] with the introduction of the Coons patches [39][40][41].The corresponding approach is addressed here, along with recent developments.

The Coons Patch and Its Three-Dimensional Extension
Throughout the rest of this paper, the region of definition of the problem under consideration consists of the union of topologically hexahedral domains, referred to as blocks, which, in analogy with (15), are described by the mapping [This allows for sufficient generality, since this subdivision may be easily obtained for any region of interest in practical applications.]Then, the intervals [  ,   ] are divided into   subintervals (not necessarily uniform, in analogy with ( 5)).By doing this, each domain (block) has been subdivided into subdomains called bricks, which are described by (15).In the rest of the paper, the term "element" is identified with the term "block."However, a block may consist of a single brick (singlebrick element).
The various formulations presented are assumed to be isoparametric, in the sense that the same type of interpolation scheme is used for the geometry and the unknown.
The finite-element family proposed here, which, as we will see, is not affected by the BVD problem, is based upon an approach introduced by Coons [39][40][41], for approximating a quadrilateral surface in terms of its four edges, namely, generating a suitable quadrilateral surface of which we know only its four edges.The resulting surface is known as the Coons patch and is illustrated here, along with its threedimensional extension.
The Coons patch technique is also known as the transfinite interpolation, introduced in [46,47].
The Coons Patch.The Coons patch was introduced to describe geometries.Accordingly, the scheme is discussed in relationship with the geometry, which is also less cumbersome to describe.This yields no loss of generality, because, as mentioned above, we use an isoparametric formulation.
Let x 0 (, ), with ,  ∈ [−1, 1], describe a generic topologically quadrilateral surface patch (see Figure 1).Let be the equations that describe the four edges of the patch, and let denote the four corner points.
As indicated above, in the Coons patch algorithm the functions x ± 1 () and x ± 2 () are assumed to be prescribed.Then, the algorithm generates a surface x  (, ) that has these lines as edges.Specifically, the Coons patch is obtained as follows: (1) consider the sum of the two linear interpolations between the two sets of opposite boundary lines; (2) from this subtract a bilinear interpolation through the four corner points.This yields (again with ,  ∈ [−1, 1]) The four edges of this surface indeed coincide with the four generating lines.For instance, we have (use  + (1) = 1 and  − (1) = 0) because

Three-Dimensional Extension of Coons Interpolation.
Here, we finally consider the three-dimensional extension of the Coons patch.Specifically, in this subsubsection, we begin by showing how to obtain a brick, by starting from its six generating faces.
Inspired by what is done to obtain (24), the function x(, , ), which describes the brick, is obtained as (1) the sum of the three linear interpolations between opposite faces, minus (2) the sum of three bilinear interpolations through three sets of four "parallel" edges, plus (3) a trilinear interpolation through the eight vertices.Stated in mathematical terms, we have International Journal of Aerospace Engineering where x(±1, , ), x(, ±1, ), and x(, , ±1) describe the geometries of the six faces, which for the time being are assumed to be prescribed.Similarly, x(1, 1, ) denotes the edge for  =  = 1, and x(1, 1, 1) denotes the vertex for  =  =  = 1.[In analogy with (25), the brick thereby generated does indeed have the prescribed generating faces.]

The Third-Order Formulation
Summarizing the results of the last subsubsection, the interpolation in (26) provides a mapping which describes the brick geometry in terms of the geometry of the six faces.We can do better.By combining the algorithm in (26) with that in (24) (namely, choosing that the six generating faces of the brick to be provided as Coons patches) one is able to generate a brick simply in terms of its twelve edges.We can do even better.We can obtain a brick that is described in terms of the location of its eight nodes, along with the base vectors there.To this end, let us begin by noting that the use of ( 24) can produce an interpolation scheme at most of order three, independently of how accurate the description the edges is.For, no matter what we do along the edges, in the Coons patch (see (24)) the fourth-order term  2  2 would be missing.Thus, without any loss of accuracy, we might as well limit ourselves to a cubic interpolation along the edges, namely, to the Hermite interpolation presented in (3).Accordingly, let us consider what happens to the threedimensional extension of the Coons patch technique, when the twelve generating edge lines are obtained by using the Hermite interpolation technique.
Coons-Hermite Patch.Let us begin by examining what happens to the Coons patch interpolation (see (24)), if we assume that the four edge lines of the patch are described by a Hermite interpolation of the type given in (3), which may be written as where x ±± are the four corner points of the patch, x (1)  ±± = x/| x ±± and x (2)  ±± = x/| x ±± are the corresponding base vectors evaluated at the four corner points, whereas  ± () and  ± () are the standard one-dimensional third-order Hermite polynomials (see ( 4)).
Combining ( 26) and ( 30) and setting x (1)  s = x/| x s , x (2)  s = x/| x s , and x (3)  s = x/| x s , one obtains x brick (, , ) = x faces (, , ) − x edges (, , ) where the three terms are defined below.The first term in the above equation, obtained by combining the first two lines on the right side of ( 26) with (30), is given by x (1)   s (32) and is the contribution from the three linear interpolations between the opposite faces.[For the first line on the right side of the above equation, use (30).For the other two, rotate indices.Again, the symbol ∑ s , where s = ( 1 ,  2 ,  3 ), with   = ±, is understood as a sum that spans over the eight values of s, which correspond to the eight vertices of the brick.]On the other hand, the second term in (31), obtained combining the next four lines in (26) with (30) (again rotate indices), is given by and is the contribution from three bilinear interpolations through "parallel" edges.Finally, the last term in (31) (last four lines in ( 26)) is given by (34) and is the contribution from the trilinear interpolation through the eight vertices.
Combining these expressions, one obtains where that is (use ( 29)), whereas Equation ( 35) describes an interpolation of x = x(  ) within the block, in terms of the values of x and of its partial derivatives x () = x/  , at the eight corners of the brick.As in the case of the Hermite brick (see (18)), the local interpolation (see (35)) may be recast in terms of global interpolation functions.
Comments.The same expression may be used to interpolate, within each brick, not only the geometry, but also the unknown function (  ), in terms of its nodal values, along with those of grad  (isoparametric formulation).Accordingly, the unknowns consist solely of the values of  and of the three Cartesian components of grad  at the nodes (to which the values of the partial derivatives  () = /  are related through (20)).As a consequence the above scheme is not affected by the BVD problem.
It should be emphasized, in order to avoid being misleading, that one or more vertices of a brick may coincide (degenerate bricks).This way, we can treat bricks that have less than eight nodes (e.g., a wedge-like brick) as if they were hexahedral bricks.In this case, the magnitudes of the corresponding covariant base vectors tend to zero and those of the contravariant base vectors tend to infinity.[This does not cause any problem, because to evaluate the integrals we use the Gaussian quadrature scheme, which does not use the values of the integrand at the end points of the interval.]The number of equations and unknowns are reduced accordingly.

The High-Order Element Proposed
Here, we want to show that the proposed formulation is not limited to the third-order scheme presented above.On the contrary, it is possible to have a whole family of finite elements that, like the third-order element, are not affected by the BVD problem.
Indeed, note that the third-order element is obtained by a combination of the two-point Hermite interpolation (namely, the third-order one) and the two-point Lagrange interpolation (namely, the first-order one).In order to introduce the high-order elements of the proposed family, we subdivide the block into bricks, as discussed at the beginning of Section 4, so as to generate additional nodes for the block.Then, we simply replace the two-node interpolation schemes, with the Lagrange and Hermite -node ones (see ( 5) and ( 9)).This way, by increasing the number , we do not decrease ℎ, which remains equal to the size of the element (block, not brick; no ℎ-convergence).What changes is the order of the scheme, because, by using  bricks in each direction, we obtain polynomials of order 2 + 1 in each direction (convergence).
Specifically, as discussed at the beginning of Section 4, for the high-order finite-element formulation under consideration, we still have that the region of definition of the problem consists of the union of blocks.We assume each block to be described by the mapping in (21).Then, the intervals [, ] are divided into  subintervals (not necessarily uniform).Specifically, let us start from the third-order brick (see (35)).
[Note that these functions satisfy interpolation conditions analogous to those in (12).For instance, we have  ℎ (  ,   ,   ) =   ℎ       , whereas their partial derivatives vanish at all the nodes.] The order of the element (namely, the order of the interpolating polynomial) is 2 + 1.Indeed, the proposed interpolation involves functions that are continuous with all its derivatives.
The comments regarding degenerate bricks (end of Section 5.1) apply here as well.

Validation, Assessment, and Applications
In this section, we validate, assess, and apply the methodology.For the sake of simplicity, here the applications of the methodology are limited to the evaluation of the natural frequencies and modes of vibration, for interior acoustics (Section 7.1) and structural dynamics (Section 7.2).The theoretical formulations (variational principles, etc.) used are also discussed.

Interior Acoustics.
Here, we present how the various finite elements formulations introduced above may be utilized to address problems in linear aeroacoustics, specifically, for the evaluation of the natural modes of vibration of the air inside a given cavity and the corresponding natural frequencies.This phenomenon is governed by the linearized wave equation for the velocity potential (x, ), which is such that k(x, ) = grad .Setting (x, ) = (x)  , we obtain the Helmholtz equation; namely, where  S =  0 / 0 denotes the undisturbed speed of sound.
The pressure is related to  by the linearized Bernoulli theorem in the frequency domain, which states that Next, consider the boundary conditions.Assume that S = S 1 ∪ S 2 , where on S 1 we have  =  0 (a condition often used for the opening of wind and brass musical instruments), whereas S 2 is a rigid wall.On S 1 , using (43), the boundary condition is given by On the other hand, on S 2 , we have a zero-normal component of the velocity; namely, The problem in (42) may be stated in variational form as .
As for any variational formulation, the boundary conditions in general may be divided into two types [28,29].The first type is the so-called essential or geometrical boundary condition (e.g., Dirichlet boundary condition,  = 0, (44)), which has to be satisfied by all the shape functions.The second type is the natural boundary condition (e.g., Neumann boundary condition, / = 0, (45)), which needs not be imposed explicitly, in that it is a direct consequence of (46), if the other type is not imposed.
Discretization of the Problem.For the finite-element formulation proposed here, set (see (18)) where   (  ) ( = 1, . . ., ) are the global shape functions, defined within each brick/block by the local interpolation functions discussed in the preceding sections, whereas   ( = 1, . . ., ) comprises all the corresponding unknowns.
[To be specific,   include the nodal values of  and grad .]Let us examine how the boundary conditions are implemented in the discretized formulation.For the first type of boundary conditions, it is sufficient to remove the unknown that vanishes, along with the corresponding equation.For the second type, no action is required, since they are automatically satisfied, in the limit, as  tends to infinity [28,29].
Following the Rayleigh-Ritz method, we combine the approximation for (x) (see (47)), with (46), to obtain where z = {  } is the vector of the unknowns and λ =  2 / 2 S , whereas the mass and stiffness matrices are, respectively, given by M = [  ] and K = [  ], with Equation (48) implies which is the equation used to obtain the unknown   .
Once the values   are obtained, the approximate modes of vibration are given by (47).
Definitions of Errors.Before presenting the numerical results, let us introduce the expression used to estimate the error in the evaluation of the eigenfunctions.Consider the leastsquare error, defined by where   (  ) and φ (  ) denote, respectively, the exact and the approximate th eigenfunction.Note that (see (47)) where   denotes the th component of the th eigenvector of (48).The vectors z  satisfies the orthonormality relation z T  Mz  =   , which is fully equivalent to the condition ∫ V φ φ dV =   , as one may verify.Similarly, for convenience, the exact eigenfunctions   (  ) are approximated with the same type of interpolation, namely, still using the shape functions   (  ), so as to have where   denotes the th component of the th vectors y  defined as follows: the components   are the nodal values corresponding to   but are obtained from the exact th eigenfunction   (  ) and are normalized by y T  My  = 1.Combining (51)-(53) and using (49) as well as z T  Mz  = y T  My  = 1, one obtains Thus, in the following, the definition of the error used in the assessment of the eigenfunctions is  On the other hand, the error on the eigenvalues, denoted by   , is defined by where   and λ denote, respectively, the exact and the approximate th eigenvalue.
Cuboidal Cavity.The first problem investigated is that of the evaluation of the natural frequencies and modes of vibration of the air inside a cuboidal cavity, with edges , , and .The walls are assumed to be rigid (natural boundary condition).
In this case, we have an exact solution; namely, (, ,  = 0, 1, . ..) . (58) Here, we consider first a cube of side ℓ =  (namely,  =  =  = ), so as to have   =  2  / 2 S =  2 +  2 +  2 .Consider first the ℎ-convergence.For the third-order results, the finite elements (blocks) coincide with the bricks (single-brick blocks).For the fifth-order results, we have /2 finite elements (blocks) in each direction, with each block being composed of two bricks in each direction, for a total of  3 /8 blocks, with eight bricks per block.Figures 2 and 3 illustrate the ℎ-convergence rate (namely, the rate of convergence obtained by reducing the size ℎ of the elements), as obtained with third and fifth-order schemes, for  111 and  222 , respectively, which are equal, respectively, to 3 and 12.The results are compared with those obtained using the ANSYS Fluid 220 element [27], used within the option "absence of air-structure coupling" ("structure absent" configuration).The ANSYS Fluid 220  element is a higher order 3D 20-node solid element that exhibits quadratic pressure behavior.In ℎ-convergence, the error behaves like ℎ  , where  denotes the order of the accuracy of the scheme, whereas ℎ = ℓ/ is the size of the element, when the elements are uniform (or the average size of the element, when the elements are nonuniform).Accordingly, the results are conveniently presented in a loglog scale, with 1/ in abscissa, as representative of the element size.The convergence rate for each scheme, as well as those obtained with the ANSYS element, is determined by its angular coefficients.[Note that the rate of convergence of the eigenvalues is much higher than the order of the interpolation scheme.Indeed, it is expected because the rate of convergence for the eigenvalues is double that of the eigenfunctions, which in theory is equal to the order of the polynomials used in the scheme.] Next, let us consider the -convergence, namely, the convergence obtained by increasing the order of the scheme.
To study this, we use a single finite element (block), with  bricks in each direction, for total number of  3 bricks.Again, the order of the polynomials used in the scheme is 2 + 1. Figures 4 and 5 present the eigenvalues  111 and  222 , as functions of the order  = 2 + 1.Similar results hold for higher eigenvalues.
Note that, in the above test case we have, for instance,  100 =  010 =  001 (degenerate eigenvalue).This phenomenon makes it difficult to compare the numerical eigenfunctions to the analytical ones (see (57)), because any linear combination of  100 ,  010 , and  001 is still an eigenfunction of the problem.In order to circumvent the problem it is sufficient to avoid having  =  = .Accordingly, the next problem considered is that of a cuboid of sides  = ,  = 0.99, and  = 1.01.
Figures 6-11 present the results for the modes (Figures 6, 8, and 10) and the corresponding eigenvalues (Figures 7, 9, and 11), as obtained using the proposed formulation with  = 1, 2, 3 elements (namely, third, fifth, and seventh orders, resp.).These correspond to (2+1) 3 = 8, 27, and 64 nodes and hence 32, 108, and 256 total unknowns (four unknowns per node) and modes, respectively.In each figure, the modes and the related eigenvalues are ordered according to the number of waves.The label [, , ] relates the th mode to its number of half-waves in each direction.Accordingly, we have (1) only one eigenvalue when  =  = ; (2) groups of three approximately equal eigenvalues when only two of them are equal (they would be identical if we had  =  = ); and (3) groups of six approximately equal eigenvalues when , , and  differ.
Figures 6, 8, and 10 estimate the projections of the numerical eigenfunction vectors z  , for  = 1, . . ., ( + 1) 3 into each of the exact eigenfunction vectors y  , for  = 1, . . ., ( + 1) 3 , through   − z  My  , a generalization of (55).The darker the mark, the poorer the agreement.On the other hand, Figures 7, 9, and 11 depict   for  = 1, . . ., ( + 1) 3 .Note that the higher the order of the element, the higher the number of eigenvalues captured.As a rule of thumb, one can state that, with respect to the total amount of degrees of freedom, one-third of the total eigenvalues have been captured accurately for this problem (the corresponding eigenfunctions present errors well below 5%).Proposed methodology (order 3) Exact eigenvalue     Cylindrical Cavity.The second problem analyzed is that of the natural frequencies and modes of vibration of the air inside a cylindrical cavity, with radius  = 1 and length ℓ = 1 (its main dimensions are reported in Figure 12).The geometry has been dealt with as a single element (block), with  bricks along each direction.The walls are again assumed to be rigid.An exact solution is available also in this case: where   () is the Bessel function of order  of the first kind.The corresponding eigenvalues are where  ()  are the solutions of    ( ()  ) = 0.The type of mesh used in this case is depicted in Figure 12.
[It may be noted that the bricks along the cylinder axis  are degenerate (wedge-like bricks).The comments regarding degenerate bricks (end of Section 5) apply here as well.]The considerations presented in connection with Figures 2-5 hold for Figures 13-16 as well.In particular, Figures 13 and 14 pertain to the ℎ-convergence of λ1,1,0 and λ2,1,0 , that is, those eigenvalues with  = 1, 2,  = 1 and  = 0. L-Shaped Prismatic Cavity.The third problem investigated pertains the evaluation of the natural frequencies and modes of vibration of the air inside the L-shaped prismatic cavity.Rigid-wall boundary conditions are assumed in this case as well.Contrary to previous test cases, exact solution is not available for this particular problem, so the results from the proposed formulation are only compared with those obtained by ANSYS.Numerical results are obtained by subdividing the region into three cubic blocks, having sides ℓ = 1 (Figure 17).For each block, use  bricks in each direction, for a total of  3 bricks for each of the three blocks.
We indicate with   the th eigenvalue of the discretized problem.
Proposed methodology (order 3) Proposed methodology (order 5) ANSYS Fluid 220   In  the convergence analysis with respect to the number of the degrees of freedom used in the FEM solver is presented.Specifically, we present the sixth natural frequency, which corresponds to the best result obtained with one) show similar levels of accuracy.However, the present results include higher accuracy ones.

Structural Dynamics.
The second problem, namely, the evaluation of the natural modes of vibration for an isotropic elastic material with continuous material properties, may be stated in variational form as where  denotes the density,    = 2[   +       ]/(1 − 2])], and   = ( / +  / )/2, with ⋅ ⋅ ⋅ / denoting covariant differentiation in the undeformed configuration, since the problem under consideration is linear.[The most general linear stress-strain relationship is given by   =     , where   includes only twenty-one independent coefficients, because of the symmetry of the stress and strain tensors, as well as from energy considerations.This general expression, used to obtain the results, is not further addressed here, because all the numerical results are limited to isotropic materials.] The boundary conditions typically considered are either those for a free surface (natural boundary condition, which require no action) or those for a fixed surface, for which the values of the nodal displacements (and their tangential derivatives) vanish.Here, for the validation of the formulation we have included the case of a free beam and a plate hinged at the boundary.
For all the numerical results regarding plates, we used a single brick along the thickness.This seems adequate to remove the shear locking phenomenon [48].Similarly, a single brick is used on the (rectangular) cross section of the beam.
Discretization of Problem.Substituting the approximation for  discussed in the preceding sections (see (18)), where now the interpolating functions are vector functions Ψ(  ), yields where, again, z = {  } is the vector of the unknown nodal values.The mass and stiffness matrices are, respectively, given by M = [  ] and K = [  ], with where    (  ) = ∑  [Ψ , ⋅ g  + Ψ , ⋅ g  ]  (  being the contravariant metric tensor components, in the undeformed configuration).
In order to assess the methodology, in the following we present some test cases that are particularly demanding, because of the above-mentioned shear locking phenomenon [48].These include plates, where one dimension is much smaller than the others, and even better beams, where two of them are much smaller than the third.These structures are especially important in aerospace and in other fields where light structures are needed.
Beams and plates are typically treated as one-dimensional and two-dimensional structures, respectively.However, in view of our desire to have a single type of element to describe the whole structure, we treat them as three-dimensional structures.
Free Beam, with Rectangular Cross Section.Preliminary results for a free beam, with length ℓ = 1, and rectangular cross section (ℎ 1 = 0.010 and ℎ 2 = 0.015) are presented.In this case, the eigenfunctions for the bending in a single plane are solutions of the differential problem where  denotes the (one-dimensional and one-directional) mode along direction 1.The free-beam boundary conditions are   (0) =   (0) = 0 and   () =   () = 0 and the solution is given by where  is the solution of the characteristic equation Then, one obtains the eigenvalue  =  4 by using (64).Again,  and λ indicate, respectively, the exact and the computed eigenvalue.The geometry is described by a sequence of  bricks along its length (a single brick is used in the other two directions).The ℎ-convergence analysis for the third and fourth eigenvalues is presented in Figures 20-23.The present results are obtained with formulations of order up to 7 and are compared with those obtained with ANSYS elements Solid 186 and Solid-Shell 190. [The former is a 20node element with three degrees of freedom per node, i.e., the three components of the displacement.The latter is an 8node element with six degrees of freedom per node, i.e., the The corresponding natural frequencies are where  denotes the bending stiffness of the plate.[Of course, the present results should be compared not to the eigenfunctions   =   (, ), but rather to the threedimensional modes Φ  (, , ) =   (, )k.]The exact eigenvalues are given by For small thicknesses, these frequencies are much lower than those relative to the in-plane motion and therefore are a good reference for the numerical results.Denoting with λ the computed eigenvalue, Figures 24-27 show the ℎ-convergence for the eigenvalues for  =  = 1 and  =  = 2, respectively, using, again, the third-order element with a single element along the thickness, for a total of  2 bricks.The computed eigenvalues λ are compared to the exact ones   and with those obtained by using the 8-node ANSYS Solid-Shell six degrees of freedom per node brick [this is specific for thin structures applications].The present formulation shows better accuracy as well.Figure 28 shows the eigenvalues, as obtained using  = 4 third-order bricks in each in-plane direction (only one brick is used in the normal direction).The approximate eigenvalues λ are compared with the exact ones   .Similarly, Figure 29  well.About one-half of the first (2 + 1) 2 modes are correctly captured.[Note that the total number of degrees of freedom is 4(2 + 1) 2 .However, the last 3(2 + 1) 2 modes are related to the rotation of the normal and the curvature of the fiber along the thickness.]

Concluding Remarks
A novel family of finite elements has been presented.The domain of interest is assumed to be the union of topologically hexahedral elements, which are called blocks.[It should be emphasized that the subdivision into blocks is important.For, it is not always convenient to have a global body-matched coordinate system.For example, consider a square-section beam, which (in analogy with a Möbius strip) is twisted by 90 ∘ degrees around the axial direction, say the -coordinate.Matching the two-end sections yields an interchange of the  and  coordinates.]The blocks are assumed to be described by a smooth mapping of the type x = x(  ) (see Appendix).The geometry of each block is described by the nodal points x  1  1  3 = x(    ), along with the corresponding base vectors g  (    ), where   = 0, . . .,   ( = 1, 2, 3), so as to have   =   + 1 nodes along the direction   .This way the element (block) is divided into subelements (bricks).The unknowns are the nodal values of the function and the Cartesian components of its gradient.The same interpolation scheme is used for both geometry and unknown (isoparametric element).The proposed interpolation family, applied to each block, is based upon a novel three-dimensional extension of the Coons patch.This extension requires a judicious combination of the -node Lagrange interpolation polynomials (of degree ) and the -node Hermite interpolation polynomials (of degree 2 + 1).The overall order of the scheme (namely, the degree of the polynomials used) is 2 + 1.In the case of high-accuracy schemes, in order to avoid the well-known instability connected with Lagrange polynomials through uniformly spaced points, the abscissas of the Gauss-Lobatto quadrature scheme are used to obtain    .[However, the standard Gaussian quadrature scheme is used for the interpolation, so as to avoid singularities on the degenerate faces (see the comments regarding degenerate bricks, the end of Section 5).] The resulting scheme guarantees the appropriate continuity between blocks.[Specifically, at the interfaces between blocks, the geometry is continuous (the in-plane components of the base vectors are also continuous).However, the outof-plane component of the base vector is allowed to be discontinuous.Nonetheless, the formulation is not affected by the BVD problem (Section 3).]This feature makes it applicable to highly complex configurations.
Numerical applications have been presented, which include the evaluation of the natural frequencies and modes of vibration of (1) air inside a cavity (interior acoustics) and (2) beams and plates, treated as three-dimensional structures (structural dynamics).[Beams and plates were used because of their importance in aerospace and because they are the most difficult to model (shear locking problem [48]).Treating them as three-dimensional structures is important, because it is envisioned that the proposed tool will be used within a Multidisciplinary Design Optimization (MDO) program [45], where resizing would be difficult, using one-dimensional models of beams and two-dimensional ones of plates.] The results indicate that the proposed methodology is highly accurate and capable of capturing efficiently relatively high frequencies, thereby resulting appropriate in design with acoustic constraints.
For low-order schemes, our results have a level of accuracy comparable to that of the code ANSYS.However, an important advantage of the proposed formulation lies in its flexibility in obtaining different level of accuracy (convergence) simply by increasing the number of nodes, as one would do for the ℎ-convergence.This fact makes it a good candidate for addressing problems pertaining to midrange frequencies.frequencies.[Preliminary results for very large values of N are presented in [36].]Further analysis in this direction is needed.Indeed, the geometries chosen for the illustrative examples are relatively simple.An L-shaped cavity has been used to validate the multiblock feature.Nonetheless, the results obtained are quite encouraging, and a deeper analysis is warranted, so as to validate the methodology for more complex configurations, before utilizing it to study the coupling of interior acoustics and structural dynamics (fluid-solid interaction), by using the methodology in [36,49,50].Only then the methodology may be considered within an MDO context [45].
Therefore, the question arises: "Can the geometry generated using the formulation above be used for the finite-element formulation of order 2+1?" "Almost, but not quite!"Indeed, the base vector along a coordinate line, say a -line, may be discontinuous at each node, albeit ever so slightly, provided that the function x(, , ) used to obtain the nodes for the first-order geometry is smooth.In this case, base vectors g  that are continuous at the nodes are obtained by taking the average of the two nodal values of g  .The final result is a block x(  ) that is continuous, with base vectors continuous at the nodes.
Of course, the base vectors are allowed to have a discontinuous out-of-plane component at the interfaces of the blocks.
Figures 15 and 16 pertain to the corresponding -convergence.