# A Novel Highly Accurate Finite-Element Family

Authors: Giovanni Bernardini; Fabio Cetta; Luigi Morino
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

A novelNth order finite element for interior acoustics and structural dynamics is presented, with N 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 (p-convergence) 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 beam- and shell-like structures as three-dimensional small-thickness elements.


---

## Body

## 1. 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 fatigue-life of the structures and hence on the maintenance costs [1–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 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 p-FEM [18–20], and C1-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–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–38], with increasing level of algorithmic sophistication. This work is to be understood as a first step towards the development of a highly accurate finite-element 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–41] and (2) the high-order Lagrange and Hermite interpolation schemes [42–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 (p-convergence) just by changing the number of nodes (as one would for h-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.


## 2. Preliminaries: Lagrange and Hermite Interpolation

The formulation presented here involves a judicious combination of high-order Lagrange and Hermite interpolation polynomials [42–44], which are briefly reviewed here. First, we address the two-point interpolation polynomials (needed in Sections 3–5), and then we consider the N-point interpolation, with N 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–44]. Consider a function u=u(η), defined over the interval [0,2N]. Let us divide the interval [0,2N] into N 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 (1)uξ≃u+L+ξ+u-L-ξ,withξ∈-1,1,where u±=u±1 denotes the values of uξ at the end points, whereas the Lagrange first-order interpolation polynomials L±ξ are given by (2)L+ξ=121+ξ,L-ξ=121-ξ.The resulting interpolated function over the whole interval [0,2N] is continuous (class C0) and piecewise linear. It will be referred to as theone-dimensional first-order Lagrange interpolation.On the other hand, within each subinterval, the standard Hermite interpolation is(3)uξ=u+M+ξ+u-M-ξ+u+′N+ξ+u-′N-ξ,withξ∈-1,1,where u± denotes the values of uξ at ξ=±1 and u±′ denotes the values u′ξ=du/dξ at ξ=±1, whereas the Hermite interpolation polynomials M±ξ and N±ξ are given by (4)M±ξ=142±3ξ∓ξ3,N±ξ=14∓1-ξ±ξ2+ξ3.The resulting interpolated function in [0,2N] is continuous with its first derivative (class C1) and is piecewise cubic. It will be referred to as theone-dimensional third-order Hermite interpolation.The
P-Node Lagrange and Hermite Interpolation Polynomials. Here, we review the P-node Lagrange and Hermite interpolation polynomials [42, 44].Let us consider first theP-node Lagrange interpolation; namely, (5)uξ=∑k=0NukLkξ,where uk=u(ξk), (6)N=P-1,is the total number of subintervals between the P nodes ξk (k=0,…,N), and (7)Lkξ=∏j=0k≠jNξ-ξjξk-ξjk=0,…,N.[Throughout the paper ξ0 and ξN coincide with the end points of the interval under consideration.]In particular, ifN=2, ξ∈[-1,1], and ξ±=±1, we have L±ξ=1/21±ξ, in agreement with (2). Note that the Lagrange interpolation polynomials Lkξ satisfy the standard interpolation condition (8)Lkξj=δkjj,k=0,…,N.It is well known that, forN>7, the Lagrange interpolation polynomials are unstable if ξk are uniformly spaced [44]. However, the instability disappears if ξk 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 high-accuracy schemes, ξk used in this paper coincide with the Gauss–Lobatto abscissas.Next, let us consider theP-node extension of the two-point 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 (9)uξ=∑k=0NukMkξ+∑k=0Nuk′Nkξ,where the polynomials Mkξ and Nkξ are given by [43] (10)Mkξ=ξ-akξk-ak∏j=0k≠jNξ-ξj2ξk-ξj2=ξ-akξk-akLkξ2,Nkξ=ξ-ξk∏j=0k≠jNξ-ξj2ξk-ξj2=ξ-ξkLkξ2,with ak obtained by imposing (11)1ξk-ak+∑j=1k≠jN2ξk-ξj=0.These polynomials satisfy the Hermite interpolation conditions (12)Mkξj=δkj,Mk′ξj=0,Nkξj=0,Nk′ξj=δkj.Indeed, it is apparent that the polynomials Mkξ and Nkξ vanish with their first derivatives at ξ=ξj, for all j≠k (see (8)). Thus, we only have to verify what happens when j=k. It is apparent that Mkξk=1. In addition, taking the logarithmic derivative of Mkξ, we have (13)Mk′ξMkξ=ddξln⁡Mkξ=1ξ-ak+∑j=0k≠jN2ξ-ξj,which yields, Mk′(ξk)=0 (see (11)). Moreover, it is apparent that Nk(ξk)=0. In addition, we have (14)Nk′ξ=Lkξ2+2ξ-ξkLkξLk′ξ,which yields Nk′(ξk)=1 (see (8)). In summary, Mk(ξ) and Nk(ξ) satisfy all the conditions in (12).The polynomials have degree2N+1=2P-1. We will refer to this as the interpolation of order 2N+1.


## 3. 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-calledHermite 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 finite-element 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 thebrick, which are described by the mapping (15)x=xξα=xξ1,ξ2,ξ3=xξ,η,ζ,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 unknownu(ξα); namely, (16)uξα=∑susPsξα+∑β=13usβQsβξα+∑β,γ∈IβγusβγRsβγξα+us123Ss123ξα.The symbol ∑s, with s=(s1,s2,s3), where sα 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 u(β) denotes the partial derivative with respect to ξβ: u(β)=∂u/∂ξβ. Similarly, we have u(βγ)=∂2u/∂ξβ∂ξγ and u(123)=∂3u/∂ξ1∂ξ2∂ξ3. In addition, note that the Hermite interpolation scheme requires only the nodal values of the function and its first partial derivatives in each direction. No second repeated derivative with respect to the same variable arises. Accordingly, the second-derivative summation spans only over the mixed derivatives, and hence Iβγ=1,2;2,3;3,1. Moreover, the term 123 is the only mixed third-order derivative. Finally, Ps(ξα), Qsβ(ξα), Rsβγ(ξα), and Ss123(ξα) are suitable products of the Hermite polynomials in (4). For instance, (17)Psξα=Ms1ξ1Ms2ξ2Ms3ξ3,Qs1ξα=Ns1ξ1Ms2ξ2Ms3ξ3,Rs12ξα=Ns1ξ1Ns2ξ2Ms3ξ3,Ss123ξα=Ns1ξ1Ns2ξ2Ns3ξ3.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 ∂3u/∂ξ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 thethird-order 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 (18)uξα=∑j=1Mzjψjξα,where the unknowns zj comprise the values of u and of all the partial derivatives mentioned above, evaluated at the M=(N+1)3 nodes. On the other hand, ψj(ξα) are the global shape functions [20], obtained by assembling the local shape functions in (17).The Base Vector Discontinuity (BVD) Problem. The third-order 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 (19)gαξα=∂x∂ξαα=1,2,3,which are tangent to the coordinate lines, are discontinuous (either in magnitude or direction). In this case, the partial derivatives of a function u(ξα) with respect to ξβ are given by (20)∂u∂ξβ=gβ·grad⁡uβ=1,2,3,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⁡u. 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 three-dimensional Hermite interpolation. Similar considerations hold for the third-order derivatives. The problem discussed above will be referred to as theBVD 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–38] with the introduction of the Coons patches [39–41]. The corresponding approach is addressed here, along with recent developments.


## 4. 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 asblocks, which, in analogy with (15), are described by the mapping (21)x=xηα=xη1,η2,η3,withηα∈aα,bα.[This allows for sufficient generality, since this subdivision may be easily obtained for any region of interest in practical applications.] Then, the intervals [aα,bα] are divided into Nα subintervals (not necessarily uniform, in analogy with (5)). By doing this, each domain (block) has been subdivided into subdomains calledbricks, 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 (single-brick element).The various formulations presented are assumed to beisoparametric, 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–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 theCoons patch and is illustrated here, along with its three-dimensional extension.The Coons patch technique is also known as thetransfinite 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.Letx0ξ,η, with ξ,η∈-1,1, describe a generic topologically quadrilateral surface patch (see Figure 1). Let (22)x=x1±η=x0±1,η,x=x2±ξ=x0ξ,±1be the equations that describe the four edges of the patch, and let (23)x+-=x01,-1=x1+-1=x2-1,x++=x01,1=x1+1=x2+1,x-+=x0-1,1=x1-1=x2+-1,x--=x0-1,-1=x1--1=x2--1denote the four corner points.Figure 1
Coons patch.As indicated above, in the Coons patch algorithm the functionsx1±(η) and x2±(ξ) are assumed to be prescribed. Then, the algorithm generates a surface xC(ξ,η) 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) (24)xCξ,η=L+ξx1+η+L-ξx1-η+L+ηx2+ξ+L-ηx2-ξ-L+ξL+ηx++-L+ξL-ηx+--L-ξL+ηx-+-L-ξL-ηx--.The four edges of this surface indeed coincide with the four generating lines. For instance, we have (use L+(1)=1 and L-(1)=0) (25)xC1,η=x1+η+L+ηx2+1+L-ηx2-1-L+ηx++-L-ηx+-=x1+η,because x++=x2+(1) and x+-=x2-1 (see (23)).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(26)xξ,η,ζ=L+ξx1,η,ζ+L+ηxξ,1,ζ+L+ζxξ,η,1+L-ξx-1,η,ζ+L-ηxξ,-1,ζ+L-ζxξ,η,-1-L+ξL+ηx1,1,ζ-L+ξL+ζx1,η,1-L+ηL+ζxξ,1,1-L+ξL-ηx1,-1,ζ-L+ξL-ζx1,η,-1-L+ηL-ζxξ,1,-1-L-ξL+ηx-1,1,ζ-L-ξL+ζx-1,η,1-L-ηL+ζxξ,-1,1-L-ξL-ηx-1,-1,ζ-L-ξL-ζx-1,η,-1-L-ηL-ζxξ,-1,-1+L+ξL+ηL+ζx1,1,1+L+ξL+ηL-ζx1,1,-1+L+ξL-ηL+ζx1,-1,1+L+ξL-ηL-ζx1,-1,-1+L-ξL+ηL+ζx-1,1,1+L-ξL+ηL-ζx-1,1,-1+L-ξL-ηL+ζx-1,-1,1+L-ξL-ηL-ζx-1,-1,-1,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, x1,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.]


## 5. 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 three-dimensional 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 (27)x1±η=x±1,η=M+ηx±++M-ηx±-+N+ηx±+2+N-ηx±-2,x2±ξ=xξ,±1=M+ξx+±+M-ξx-±+N+ξx+±1+N-ξx-±1,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 M±(ξ) and N±(ξ) are the standard one-dimensional third-order Hermite polynomials (see (4)).Combining (24) and (27) yields(28)xCξ,η=L+ξM+η+M+ξL+η-L+ξL+ηx+++L+ξM-η+M+ξL-η-L+ξL-ηx+-+L-ξM+η+M-ξL+η-L-ξL+ηx-++L-ξM-η+M-ξL-η-L-ξL-ηx--+N+ξL+ηx++1+N+ξL-ηx+-1+N-ξL+ηx-+1+N-ξL-ηx--1+L+ξN+ηx++2+L+ξN-ηx+-2+L-ξN+ηx-+2+L-ξN-ηx--2,where the functions L±(ξ), M±(ξ), and N±(ξ) are given by (2) and (3). Introducing the functions (29)Ps1s2ξ,η=Ls1ξMs2η+Ms1ξLs2η-Ls1ξLs2η,where s1 and s2 stand for either + or −, (28) may be rewritten as (30)xξ,η=∑s1,s2Ps1s2ξ,ηxs1s2+∑s1,s2Ns1ξLs2ηxs1s21+∑s1,s2Ls1ξNs2ηxs1s22.


### 5.1. The Third-Order Brick

Next, let us turn to the three-dimensional extension (i.e., the expression for a brick). Combining (26) and (30) and setting xs(1)=∂x/∂ξ|xs, xs(2)=∂x/∂η|xs, and xs(3)=∂x/∂ζ|xs, one obtains (31)xbrickξ,η,ζ=xfacesξ,η,ζ-xedgesξ,η,ζ+xverticesξ,η,ζ,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(32)xfacesξ,η,ζ=∑sPs1s2ξ,ηLs3ζxs+∑sNs1ξLs2ηLs3ζxs1+∑sLs1ξNs2ηLs3ζxs2+∑sPs2s3η,ζLs1ξxs+∑sNs2ηLs3ζLs1ξxs2+∑sLs2ηNs3ζLs1ξxs3+∑sPs3s1ζ,ξLs2ηxs+∑sNs3ζLs1ξLs2ηxs3+∑sLs3ζNs1ξLs2ηxs1and 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=(s1,s2,s3), with sα=±, 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 (33)xedgesξ,η,ζ=∑sMs1ξxs+Ns1ξxs1Ls2ηLs3ζ+∑sMs2ηxs+Ns2ηxs2Ls3ζLs1ξ+∑sMs3ζxs+Ns3ζxs3Ls1ξLs2η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)xverticesξ,η,ζ=∑sLs1ξLs2ηLs3ζxsand is the contribution from the trilinear interpolation through the eight vertices.Combining these expressions, one obtains(35)xbrickξ,η,ζ=∑sxsRsξ,η,ζ+∑s∑α=13xsαSsαξ,η,ζ,where (36)Rsξ,η,ζ=Ps1s2ξ,ηLs3ζ+Ps2s3η,ζLs1ξ+Ps3s1ζ,ξLs2η-Ms1ξLs2ηLs3ζ-Ls1ξMs2ηLs3ζ-Ls1ξLs2ηMs3ζ+Ls1ξLs3ηLs3ζ;that is (use (29)), (37)Rsξ,η,ζ=Ms1ξLs2ηLs3ζ+Ms2ηLs3ζLs1ξ+Ms3ζLs1ξLs2η-2Ls1ξLs2ηLs3ζ,whereas (38)Ss1ξ,η,ζ=Ns1ξLs2ηLs3ζ,Ss2ξ,η,ζ=Ns2ηLs3ζLs1ξ,Ss3ξ,η,ζ=Ns3ζLs1ξLs2η.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 u(ξα), in terms of its nodal values, along with those of grad⁡u (isoparametric formulation). Accordingly, the unknowns consist solely of the values of u and of the three Cartesian components of grad⁡u at the nodes (to which the values of the partial derivatives u(α)=∂u/∂ξα 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.


## 5.1. The Third-Order Brick

Next, let us turn to the three-dimensional extension (i.e., the expression for a brick). Combining (26) and (30) and setting xs(1)=∂x/∂ξ|xs, xs(2)=∂x/∂η|xs, and xs(3)=∂x/∂ζ|xs, one obtains (31)xbrickξ,η,ζ=xfacesξ,η,ζ-xedgesξ,η,ζ+xverticesξ,η,ζ,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(32)xfacesξ,η,ζ=∑sPs1s2ξ,ηLs3ζxs+∑sNs1ξLs2ηLs3ζxs1+∑sLs1ξNs2ηLs3ζxs2+∑sPs2s3η,ζLs1ξxs+∑sNs2ηLs3ζLs1ξxs2+∑sLs2ηNs3ζLs1ξxs3+∑sPs3s1ζ,ξLs2ηxs+∑sNs3ζLs1ξLs2ηxs3+∑sLs3ζNs1ξLs2ηxs1and 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=(s1,s2,s3), with sα=±, 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 (33)xedgesξ,η,ζ=∑sMs1ξxs+Ns1ξxs1Ls2ηLs3ζ+∑sMs2ηxs+Ns2ηxs2Ls3ζLs1ξ+∑sMs3ζxs+Ns3ζxs3Ls1ξLs2η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)xverticesξ,η,ζ=∑sLs1ξLs2ηLs3ζxsand is the contribution from the trilinear interpolation through the eight vertices.Combining these expressions, one obtains(35)xbrickξ,η,ζ=∑sxsRsξ,η,ζ+∑s∑α=13xsαSsαξ,η,ζ,where (36)Rsξ,η,ζ=Ps1s2ξ,ηLs3ζ+Ps2s3η,ζLs1ξ+Ps3s1ζ,ξLs2η-Ms1ξLs2ηLs3ζ-Ls1ξMs2ηLs3ζ-Ls1ξLs2ηMs3ζ+Ls1ξLs3ηLs3ζ;that is (use (29)), (37)Rsξ,η,ζ=Ms1ξLs2ηLs3ζ+Ms2ηLs3ζLs1ξ+Ms3ζLs1ξLs2η-2Ls1ξLs2ηLs3ζ,whereas (38)Ss1ξ,η,ζ=Ns1ξLs2ηLs3ζ,Ss2ξ,η,ζ=Ns2ηLs3ζLs1ξ,Ss3ξ,η,ζ=Ns3ζLs1ξLs2η.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 u(ξα), in terms of its nodal values, along with those of grad⁡u (isoparametric formulation). Accordingly, the unknowns consist solely of the values of u and of the three Cartesian components of grad⁡u at the nodes (to which the values of the partial derivatives u(α)=∂u/∂ξα 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.


## 6. 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 Section4, so as to generate additional nodes for the block. Then, we simply replace the two-node interpolation schemes, with the Lagrange and Hermite P-node ones (see (5) and (9)). This way, by increasing the number N, we do not decrease h, which remains equal to the size of the element (block, not brick; no h-convergence). What changes is the order of the scheme, because, by using N bricks in each direction, we obtain polynomials of order 2N+1 in each direction (p-convergence).Specifically, as discussed at the beginning of Section4, 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 [a,b] are divided into N subintervals (not necessarily uniform). Specifically, let us start from the third-order brick (see (35)). For each block, we replace the two-node Lagrange and Hermite interpolations, with the P-node ones (see (5) and (9)). Then, we can extend (35), to have (39)xblockξ,η,ζ=∑h,j,k=0NxhjkRhjkξ,η,ζ+∑h,j,k=0N∑α=13xhjkαShjkαξ,η,ζ,where (40)Rhjkξ,η,ζ=MhξLjηLkζ+MjηLkζLhξ+MkζLhξLjη-2LhξLjηLkζ,whereas (41)Ss1ξ,η,ζ=NhξLjηLkζ,Ss2ξ,η,ζ=NjηLkζLhξ,Ss3ξ,η,ζ=NkζLhξLjη.[Note that these functions satisfy interpolation conditions analogous to those in (12). For instance, we have Rhjk(ξl,ηm,ζn)=δhlδjmδkn, whereas their partial derivatives vanish at all the nodes.]The order of the element (namely, the order of the interpolating polynomial) is2N+1. Indeed, the proposed interpolation involves functions that are continuous with all its derivatives.The comments regarding degenerate bricks (end of Section5.1) apply here as well.


## 7. 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 (Section7.1) and structural dynamics (Section 7.2). The theoretical formulations (variational principles, etc.) used are also discussed.


### 7.1. 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,t), which is such that v(x,t)=grad⁡φ. Setting φ(x,t)=ϕ(x)eıωt, we obtain the Helmholtz equation; namely, (42)∇2ϕ+ω2cS2ϕ=0,where cS=γp0/ρ0 denotes the undisturbed speed of sound.The pressure is related toϕ by the linearized Bernoulli theorem in the frequency domain, which states that (43)ıωϕ+pρ0=p0ρ0.Next, consider the boundary conditions. Assume thatS=S1∪S2, where on S1 we have p=p0 (a condition often used for the opening of wind and brass musical instruments), whereas S2 is a rigid wall. On S1, using (43), the boundary condition is given by (44)ϕ=0x∈S2.On the other hand, on S2, we have a zero-normal component of the velocity; namely, (45)∂ϕ∂n=0x∈S1.The problem in (42) may be stated in variational form as (46)12ω2cS2∫Vϕ2dV-12∫Vgrad⁡ϕ2dV=stationaryϕx.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, ∂ϕ/∂n=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)) (47)ϕ^ξα=∑j=1Mzjψjξα,where ψj(ξα) (j=1,…,M) are the global shape functions, defined within each brick/block by the local interpolation functions discussed in the preceding sections, whereas zj (j=1,…,M) comprises all the corresponding unknowns. [To be specific, zj 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, asN tends to infinity [28, 29].Following the Rayleigh–Ritz method, we combine the approximation forϕ(x) (see (47)), with (46), to obtain (48)12zTKz-λ^12zTMz=stationaryz,where z=zj is the vector of the unknowns and λ^=ω2/cS2, whereas the mass and stiffness matrices are, respectively, given by M=[mmq] and K=[kmq], with (49)mmq=∫VψmxψqxdV,kmq=∫Vgrad⁡ψmx·grad⁡ψqxdV.Equation (48) implies (50)Kz=λ^Mz,which is the equation used to obtain the unknown zj. Once the values zj 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 least-square error, defined by (51)ϵk=∫Vϕk-ϕ^k2dV1/2,where ϕk(ξα) and ϕ^k(ξα) denote, respectively, the exact and the approximate kth eigenfunction. Note that (see (47)) (52)ϕ^kξα=∑j=1Mzjkχjξα,where zjk denotes the jth component of the kth eigenvector of (48). The vectors zk satisfies the orthonormality relation zjTMzk=δjk, which is fully equivalent to the condition ∫V‍ϕ^jϕ^kdV=δjk, as one may verify. Similarly, for convenience, the exact eigenfunctions ϕk(ξα) are approximated with the same type of interpolation, namely, still using the shape functions χj(ξα), so as to have (53)ϕkξα=∑j=1Myjkχjξα,where yjk denotes the jth component of the kth vectors yk defined as follows: the components yjk are the nodal values corresponding to zjk but are obtained from the exact kth eigenfunction ϕk(ξα) and are normalized by ykTMyk=1. Combining (51)–(53) and using (49) as well as zkTMzk=ykTMyk=1, one obtains (54)ϵk2=∑h,j=1Myhk-zhk∫VψhψjdVyjk-zjk=yk-zkTMyk-zk=21-zkTMyk.Thus, in the following, the definition of the error used in the assessment of the eigenfunctions is (55)ϵk=21-zkTMyk.On the other hand, the error on the eigenvalues, denoted byek, is defined by (56)ek=λk-λ^k,where λk and λ^k denote, respectively, the exact and the approximate kth 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 a, b, and c. The walls are assumed to be rigid (natural boundary condition). In this case, we have an exact solution; namely, (57)ϕmnq=cos⁡mπxacos⁡nπybcos⁡qπzcm,n,q=0,1,…,(58)λmnq=ωmnq2cS2=m2π2a2+n2π2b2+q2π2c2m,n,q=0,1,….Here, we consider first a cube of sidel=π (namely, a=b=c=π), so as to have λmnq=ωmnq2/cS2=m2+n2+q2.Consider first theh-convergence. For the third-order results, the finite elements (blocks) coincide with the bricks (single-brick blocks). For the fifth-order results, we have N/2 finite elements (blocks) in each direction, with each block being composed of two bricks in each direction, for a total of N3/8 blocks, with eight bricks per block. Figures 2 and 3 illustrate the h-convergence rate (namely, the rate of convergence obtained by reducing the size h 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 h-convergence, the error behaves like hn, where n denotes the order of the accuracy of the scheme, whereas h=l/N 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 alog-log scale, with 1/N 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.]Figure 2
Cube:λ^111 versus 1/N (h-convergence).Figure 3
Cube:λ^222 versus 1/N (h-convergence).Next, let us consider thep-convergence, namely, the convergence obtained by increasing the order of the scheme. To study this, we use a single finite element (block), with N bricks in each direction, for total number of N3 bricks. Again, the order of the polynomials used in the scheme is 2N+1. Figures 4 and 5 present the eigenvalues λ111 and λ222, as functions of the order p=2N+1. Similar results hold for higher eigenvalues.Figure 4
Cube:p-convergence of λ^111.Figure 5
Cube:p-convergence of λ^222.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 a=b=c. Accordingly, the next problem considered is that of a cuboid of sides a=π, b=0.99π, and c=1.01π.Figures6–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 N=1,2,3 elements (namely, third, fifth, and seventh orders, resp.). These correspond to (2N+1)3=8, 27, and 64 nodes and hence 32, 108, and 256 total unknowns (four unknowns per node) and modes, respectively.Figure 6
Cuboid modes (3rd order,N=1).Figure 7
Cuboid eigenvalues (3rd order,N=1).Figure 8
Cuboid modes (5th order,N=2).Figure 9
Cuboid eigenvalues (5th order,N=2).Figure 10
Cuboid modes (7th order,N=3).Figure 11
Cuboid eigenvalues (7th order,N=3).In each figure, the modes and the related eigenvalues are ordered according to the number of waves. The label[k,l,m] relates the ith mode to its number of half-waves in each direction. Accordingly, we have (1) only one eigenvalue when k=l=m; (2) groups of three approximately equal eigenvalues when only two of them are equal (they would be identical if we had a=b=c); and (3) groups of six approximately equal eigenvalues when k, l, and m differ.Figures6, 8, and 10 estimate the projections of the numerical eigenfunction vectors zj, for j=1,…,(N+1)3 into each of the exact eigenfunction vectors yi, for i=1,…,(N+1)3, through δji-zjMyi, a generalization of (55). The darker the mark, the poorer the agreement. On the other hand, Figures 7, 9, and 11 depict λj for j=1,…,(N+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%).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 R=1 and length l=1 (its main dimensions are reported in Figure 12). The geometry has been dealt with as a single element (block), with N bricks along each direction. The walls are again assumed to be rigid.Figure 12
Cylindrical cavity.An exact solution is available also in this case:(59)ϕmnq=JmμnmrAmcos⁡mθ+Bmsin⁡mθcos⁡qπzl,where Jm(x) is the Bessel function of order m of the first kind. The corresponding eigenvalues are (60)λmnq=μnm+q2π2l2,where μn(m) are the solutions of Jm′(μn(m))=0.The type of mesh used in this case is depicted in Figure12. [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 Figures2–5 hold for Figures 13–16 as well. In particular, Figures 13 and 14 pertain to the h-convergence of λ^1,1,0 and λ^2,1,0, that is, those eigenvalues with m=1,2, n=1 and q=0. Figures 15 and 16 pertain to the corresponding p-convergence.Figure 13
Cylinder:λ^110 versus 1/N (h-convergence).Figure 14
Cylinder:λ^210 versus 1/N (h-convergence).Figure 15
Cylinder:p-convergence of λ^110.Figure 16
Cylinder:p-convergence of λ^210.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 l=1 (Figure 17). For each block, use N bricks in each direction, for a total of N3 bricks for each of the three blocks.Figure 17
L-shaped prismatic cavity.We indicate withλi the ith eigenvalue of the discretized problem.In Figures18-19 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 the present formulation. The corresponding results have been compared with those obtained by using the ANSYS Fluid 220 element. From the analysis of the figures, we see that, for the low-order case, both formulations (ANSYS and the present one) show similar levels of accuracy. However, the present results include higher accuracy ones.Figure 18
L-shaped cavity:λ^6 versus 1/NDOF.Figure 19
L-shaped cavity:λ^6 versus p.


### 7.2. 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(61)12ω2∫Vρu2dV-12∫VσαβϵαβdV=stationaryux,where ρ denotes the density, σαβ=2G[ϵαβ+ϵγγδαβν/(1-2ν)], and ϵαβ=(uα/β+uβ/α)/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 σαβ=cαβγδϵγδ, where cαβγδ 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 u discussed in the preceding sections (see (18)), where now the interpolating functions are vector functions Ψξα, yields (62)Kz=ω2Mz,where, again, z={zn} is the vector of the unknown nodal values. The mass and stiffness matrices are, respectively, given by M=[mmn] and K=[kmn], with (63)mmn=∫VρΨmξα·ΨnξαdV,kmn=∑α,γ∫V2GPαγPγα+ν/1-2νPααPγγdV,where Pαγ(ξρ)=∑β‍Ψm,α·gβ+Ψn,β·gα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 l=1, and rectangular cross section (h1=0.010 and h2=0.015) are presented. In this case, the eigenfunctions for the bending in a single plane are solutions of the differential problem (64)ϕIV=λϕ=ω2=μ4ϕ,where ϕ denotes the (one-dimensional and one-directional) mode along direction 1. The free-beam boundary conditions are ϕ′′(0)=ϕ′′′(0)=0 and ϕ′′(l)=ϕ′′′(l)=0 and the solution is given by (65)ϕx=Ccos⁡μx+cosh⁡μxcosh⁡μl-cos⁡μl+sin⁡μx+sinh⁡μxsin⁡μl-sinh⁡μl,where μ is the solution of the characteristic equation (66)cos⁡μcosh⁡μ=1.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 N bricks along its length (a single brick is used in the other two directions). The h-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 20-node element with three degrees of freedom per node, i.e., the three components of the displacement. The latter is an 8-node element with six degrees of freedom per node, i.e., the three components of the displacement along with the three rotations.]Figure 20
Free beam:λ^3 versus 1/NDOF (h-convergence).Figure 21
Free beam:λ^3 versus 1/N (h-convergence).Figure 22
Free beam:λ^4 versus 1/NDOF (h-convergence).Figure 23
Free beam:λ^4 versus 1/N (h-convergence).Hinged Rectangular Plate. Consider the results relative to a rectangular hinged plate with a=1.00,    b=1.01, and h=0.01. The eigenfunctions for a plate in bending are(67)ϕmnx,y=sin⁡mπxasin⁡nπyb.The corresponding natural frequencies are (68)ωmn=Dρmπa2+nπb2,where D denotes the bending stiffness of the plate. [Of course, the present results should be compared not to the eigenfunctions ϕmn=ϕmn(x,y), but rather to the three-dimensional modes Φmn(x,y,z)=ϕmn(x,y)k.]The exact eigenvalues are given by(69)λmn=ωmn2.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 λ^mn the computed eigenvalue, Figures 24–27 show the h-convergence for the eigenvalues for m=n=1 and m=n=2, respectively, using, again, the third-order element with a single element along the thickness, for a total of N2 bricks. The computed eigenvalues λ^mn are compared to the exact ones λmn 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 N=4 third-order bricks in each in-plane direction (only one brick is used in the normal direction). The approximate eigenvalues λ^mn are compared with the exact ones λmn. Similarly, Figure 29 pertains to the calculated eigenvalues obtained using the third-order element with N=8. The considerations presented for acoustics, regarding the number of accurate frequencies, apply to the present case as well. About one-half of the first (2N+1)2 modes are correctly captured. [Note that the total number of degrees of freedom is 4(2N+1)2. However, the last 3(2N+1)2 modes are related to the rotation of the normal and the curvature of the fiber along the thickness.]Figure 24
Hinged rectangular plate:λ^11 versus 1/NDOF (h-convergence).Figure 25
Hinged rectangular plate:λ^11 versus 1/N (h-convergence).Figure 26
Hinged rectangular plate:λ^22 versus 1/NDOF (h-convergence).Figure 27
Hinged rectangular plate:λ^22 versus 1/N (h-convergence).Figure 28
Eigenvalues (hinged rectangular plate; third-order brick,N=4).Figure 29
Eigenvalues (hinged rectangular plate; third-order brick,N=8).


## 7.1. 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,t), which is such that v(x,t)=grad⁡φ. Setting φ(x,t)=ϕ(x)eıωt, we obtain the Helmholtz equation; namely, (42)∇2ϕ+ω2cS2ϕ=0,where cS=γp0/ρ0 denotes the undisturbed speed of sound.The pressure is related toϕ by the linearized Bernoulli theorem in the frequency domain, which states that (43)ıωϕ+pρ0=p0ρ0.Next, consider the boundary conditions. Assume thatS=S1∪S2, where on S1 we have p=p0 (a condition often used for the opening of wind and brass musical instruments), whereas S2 is a rigid wall. On S1, using (43), the boundary condition is given by (44)ϕ=0x∈S2.On the other hand, on S2, we have a zero-normal component of the velocity; namely, (45)∂ϕ∂n=0x∈S1.The problem in (42) may be stated in variational form as (46)12ω2cS2∫Vϕ2dV-12∫Vgrad⁡ϕ2dV=stationaryϕx.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, ∂ϕ/∂n=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)) (47)ϕ^ξα=∑j=1Mzjψjξα,where ψj(ξα) (j=1,…,M) are the global shape functions, defined within each brick/block by the local interpolation functions discussed in the preceding sections, whereas zj (j=1,…,M) comprises all the corresponding unknowns. [To be specific, zj 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, asN tends to infinity [28, 29].Following the Rayleigh–Ritz method, we combine the approximation forϕ(x) (see (47)), with (46), to obtain (48)12zTKz-λ^12zTMz=stationaryz,where z=zj is the vector of the unknowns and λ^=ω2/cS2, whereas the mass and stiffness matrices are, respectively, given by M=[mmq] and K=[kmq], with (49)mmq=∫VψmxψqxdV,kmq=∫Vgrad⁡ψmx·grad⁡ψqxdV.Equation (48) implies (50)Kz=λ^Mz,which is the equation used to obtain the unknown zj. Once the values zj 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 least-square error, defined by (51)ϵk=∫Vϕk-ϕ^k2dV1/2,where ϕk(ξα) and ϕ^k(ξα) denote, respectively, the exact and the approximate kth eigenfunction. Note that (see (47)) (52)ϕ^kξα=∑j=1Mzjkχjξα,where zjk denotes the jth component of the kth eigenvector of (48). The vectors zk satisfies the orthonormality relation zjTMzk=δjk, which is fully equivalent to the condition ∫V‍ϕ^jϕ^kdV=δjk, as one may verify. Similarly, for convenience, the exact eigenfunctions ϕk(ξα) are approximated with the same type of interpolation, namely, still using the shape functions χj(ξα), so as to have (53)ϕkξα=∑j=1Myjkχjξα,where yjk denotes the jth component of the kth vectors yk defined as follows: the components yjk are the nodal values corresponding to zjk but are obtained from the exact kth eigenfunction ϕk(ξα) and are normalized by ykTMyk=1. Combining (51)–(53) and using (49) as well as zkTMzk=ykTMyk=1, one obtains (54)ϵk2=∑h,j=1Myhk-zhk∫VψhψjdVyjk-zjk=yk-zkTMyk-zk=21-zkTMyk.Thus, in the following, the definition of the error used in the assessment of the eigenfunctions is (55)ϵk=21-zkTMyk.On the other hand, the error on the eigenvalues, denoted byek, is defined by (56)ek=λk-λ^k,where λk and λ^k denote, respectively, the exact and the approximate kth 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 a, b, and c. The walls are assumed to be rigid (natural boundary condition). In this case, we have an exact solution; namely, (57)ϕmnq=cos⁡mπxacos⁡nπybcos⁡qπzcm,n,q=0,1,…,(58)λmnq=ωmnq2cS2=m2π2a2+n2π2b2+q2π2c2m,n,q=0,1,….Here, we consider first a cube of sidel=π (namely, a=b=c=π), so as to have λmnq=ωmnq2/cS2=m2+n2+q2.Consider first theh-convergence. For the third-order results, the finite elements (blocks) coincide with the bricks (single-brick blocks). For the fifth-order results, we have N/2 finite elements (blocks) in each direction, with each block being composed of two bricks in each direction, for a total of N3/8 blocks, with eight bricks per block. Figures 2 and 3 illustrate the h-convergence rate (namely, the rate of convergence obtained by reducing the size h 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 h-convergence, the error behaves like hn, where n denotes the order of the accuracy of the scheme, whereas h=l/N 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 alog-log scale, with 1/N 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.]Figure 2
Cube:λ^111 versus 1/N (h-convergence).Figure 3
Cube:λ^222 versus 1/N (h-convergence).Next, let us consider thep-convergence, namely, the convergence obtained by increasing the order of the scheme. To study this, we use a single finite element (block), with N bricks in each direction, for total number of N3 bricks. Again, the order of the polynomials used in the scheme is 2N+1. Figures 4 and 5 present the eigenvalues λ111 and λ222, as functions of the order p=2N+1. Similar results hold for higher eigenvalues.Figure 4
Cube:p-convergence of λ^111.Figure 5
Cube:p-convergence of λ^222.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 a=b=c. Accordingly, the next problem considered is that of a cuboid of sides a=π, b=0.99π, and c=1.01π.Figures6–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 N=1,2,3 elements (namely, third, fifth, and seventh orders, resp.). These correspond to (2N+1)3=8, 27, and 64 nodes and hence 32, 108, and 256 total unknowns (four unknowns per node) and modes, respectively.Figure 6
Cuboid modes (3rd order,N=1).Figure 7
Cuboid eigenvalues (3rd order,N=1).Figure 8
Cuboid modes (5th order,N=2).Figure 9
Cuboid eigenvalues (5th order,N=2).Figure 10
Cuboid modes (7th order,N=3).Figure 11
Cuboid eigenvalues (7th order,N=3).In each figure, the modes and the related eigenvalues are ordered according to the number of waves. The label[k,l,m] relates the ith mode to its number of half-waves in each direction. Accordingly, we have (1) only one eigenvalue when k=l=m; (2) groups of three approximately equal eigenvalues when only two of them are equal (they would be identical if we had a=b=c); and (3) groups of six approximately equal eigenvalues when k, l, and m differ.Figures6, 8, and 10 estimate the projections of the numerical eigenfunction vectors zj, for j=1,…,(N+1)3 into each of the exact eigenfunction vectors yi, for i=1,…,(N+1)3, through δji-zjMyi, a generalization of (55). The darker the mark, the poorer the agreement. On the other hand, Figures 7, 9, and 11 depict λj for j=1,…,(N+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%).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 R=1 and length l=1 (its main dimensions are reported in Figure 12). The geometry has been dealt with as a single element (block), with N bricks along each direction. The walls are again assumed to be rigid.Figure 12
Cylindrical cavity.An exact solution is available also in this case:(59)ϕmnq=JmμnmrAmcos⁡mθ+Bmsin⁡mθcos⁡qπzl,where Jm(x) is the Bessel function of order m of the first kind. The corresponding eigenvalues are (60)λmnq=μnm+q2π2l2,where μn(m) are the solutions of Jm′(μn(m))=0.The type of mesh used in this case is depicted in Figure12. [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 Figures2–5 hold for Figures 13–16 as well. In particular, Figures 13 and 14 pertain to the h-convergence of λ^1,1,0 and λ^2,1,0, that is, those eigenvalues with m=1,2, n=1 and q=0. Figures 15 and 16 pertain to the corresponding p-convergence.Figure 13
Cylinder:λ^110 versus 1/N (h-convergence).Figure 14
Cylinder:λ^210 versus 1/N (h-convergence).Figure 15
Cylinder:p-convergence of λ^110.Figure 16
Cylinder:p-convergence of λ^210.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 l=1 (Figure 17). For each block, use N bricks in each direction, for a total of N3 bricks for each of the three blocks.Figure 17
L-shaped prismatic cavity.We indicate withλi the ith eigenvalue of the discretized problem.In Figures18-19 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 the present formulation. The corresponding results have been compared with those obtained by using the ANSYS Fluid 220 element. From the analysis of the figures, we see that, for the low-order case, both formulations (ANSYS and the present one) show similar levels of accuracy. However, the present results include higher accuracy ones.Figure 18
L-shaped cavity:λ^6 versus 1/NDOF.Figure 19
L-shaped cavity:λ^6 versus p.


## 7.2. 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(61)12ω2∫Vρu2dV-12∫VσαβϵαβdV=stationaryux,where ρ denotes the density, σαβ=2G[ϵαβ+ϵγγδαβν/(1-2ν)], and ϵαβ=(uα/β+uβ/α)/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 σαβ=cαβγδϵγδ, where cαβγδ 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 u discussed in the preceding sections (see (18)), where now the interpolating functions are vector functions Ψξα, yields (62)Kz=ω2Mz,where, again, z={zn} is the vector of the unknown nodal values. The mass and stiffness matrices are, respectively, given by M=[mmn] and K=[kmn], with (63)mmn=∫VρΨmξα·ΨnξαdV,kmn=∑α,γ∫V2GPαγPγα+ν/1-2νPααPγγdV,where Pαγ(ξρ)=∑β‍Ψm,α·gβ+Ψn,β·gα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 l=1, and rectangular cross section (h1=0.010 and h2=0.015) are presented. In this case, the eigenfunctions for the bending in a single plane are solutions of the differential problem (64)ϕIV=λϕ=ω2=μ4ϕ,where ϕ denotes the (one-dimensional and one-directional) mode along direction 1. The free-beam boundary conditions are ϕ′′(0)=ϕ′′′(0)=0 and ϕ′′(l)=ϕ′′′(l)=0 and the solution is given by (65)ϕx=Ccos⁡μx+cosh⁡μxcosh⁡μl-cos⁡μl+sin⁡μx+sinh⁡μxsin⁡μl-sinh⁡μl,where μ is the solution of the characteristic equation (66)cos⁡μcosh⁡μ=1.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 N bricks along its length (a single brick is used in the other two directions). The h-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 20-node element with three degrees of freedom per node, i.e., the three components of the displacement. The latter is an 8-node element with six degrees of freedom per node, i.e., the three components of the displacement along with the three rotations.]Figure 20
Free beam:λ^3 versus 1/NDOF (h-convergence).Figure 21
Free beam:λ^3 versus 1/N (h-convergence).Figure 22
Free beam:λ^4 versus 1/NDOF (h-convergence).Figure 23
Free beam:λ^4 versus 1/N (h-convergence).Hinged Rectangular Plate. Consider the results relative to a rectangular hinged plate with a=1.00,    b=1.01, and h=0.01. The eigenfunctions for a plate in bending are(67)ϕmnx,y=sin⁡mπxasin⁡nπyb.The corresponding natural frequencies are (68)ωmn=Dρmπa2+nπb2,where D denotes the bending stiffness of the plate. [Of course, the present results should be compared not to the eigenfunctions ϕmn=ϕmn(x,y), but rather to the three-dimensional modes Φmn(x,y,z)=ϕmn(x,y)k.]The exact eigenvalues are given by(69)λmn=ωmn2.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 λ^mn the computed eigenvalue, Figures 24–27 show the h-convergence for the eigenvalues for m=n=1 and m=n=2, respectively, using, again, the third-order element with a single element along the thickness, for a total of N2 bricks. The computed eigenvalues λ^mn are compared to the exact ones λmn 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 N=4 third-order bricks in each in-plane direction (only one brick is used in the normal direction). The approximate eigenvalues λ^mn are compared with the exact ones λmn. Similarly, Figure 29 pertains to the calculated eigenvalues obtained using the third-order element with N=8. The considerations presented for acoustics, regarding the number of accurate frequencies, apply to the present case as well. About one-half of the first (2N+1)2 modes are correctly captured. [Note that the total number of degrees of freedom is 4(2N+1)2. However, the last 3(2N+1)2 modes are related to the rotation of the normal and the curvature of the fiber along the thickness.]Figure 24
Hinged rectangular plate:λ^11 versus 1/NDOF (h-convergence).Figure 25
Hinged rectangular plate:λ^11 versus 1/N (h-convergence).Figure 26
Hinged rectangular plate:λ^22 versus 1/NDOF (h-convergence).Figure 27
Hinged rectangular plate:λ^22 versus 1/N (h-convergence).Figure 28
Eigenvalues (hinged rectangular plate; third-order brick,N=4).Figure 29
Eigenvalues (hinged rectangular plate; third-order brick,N=8).


## 8. 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 calledblocks. [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 typex=x(ξα) (see Appendix). The geometry of each block is described by the nodal points xk1k1k3=x(ξkαα), along with the corresponding base vectors gβ(ξkαα), where kα=0,…,Nα (α=1,2,3), so as to have Pα=Nα+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 P-node Lagrange interpolation polynomials (of degree N) and the P-node Hermite interpolation polynomials (of degree 2N+1). The overall order of the scheme (namely, the degree of the polynomials used) is 2N+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 ξkα. [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 out-of-plane component of the base vector is allowed to be discontinuous. Nonetheless, the formulation is not affected by the BVD problem (Section3).] 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 (p-convergence) simply by increasing the number of nodes, as one would do for the h-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].



---

*Source: 7293682-2017-09-05.xml*