Abstract
A novel th order finite element for interior acoustics and structural dynamics is presented, with 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 (-convergence) simply by increasing the number of nodes, as one would do for -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.
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 -FEM [18–20], and -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: () the three-dimensional extension of the Coons patch technique [39–41] and () 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 (-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: () 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 () 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 -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–44]. Consider a function , defined over the interval . Let us divide the interval into subintervals, for which we introduce a local coordinate , so as to have that the end points of each subinterval are given by . Within each subinterval, the Lagrange linear interpolation is given by where 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 is continuous (class ) 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 and denotes the values at , whereas the Hermite interpolation polynomials and are given by The resulting interpolated function in is continuous with its first derivative (class ) and is piecewise cubic. It will be referred to as the one-dimensional third-order Hermite interpolation.
The -Node Lagrange and Hermite Interpolation Polynomials. Here, we review the -node Lagrange and Hermite interpolation polynomials [42, 44].
Let us consider first the -node Lagrange interpolation; namely, where , is the total number of subintervals between the nodes (), and [Throughout the paper and coincide with the end points of the interval under consideration.]
In particular, if , , and , we have , in agreement with (2). Note that the Lagrange interpolation polynomials satisfy the standard interpolation condition
It is well known that, for , 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 high-accuracy schemes, used in this paper coincide with the Gauss–Lobatto abscissas.
Next, let us consider the -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 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 . In addition, taking the logarithmic derivative of , we have which yields, (see (11)). Moreover, it is apparent that . In addition, we have which yields (see (8)). In summary, and satisfy all the conditions in (12).
The polynomials have degree . We will refer to this as the interpolation of order .
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-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 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 the brick, which are described by the mapping where () are curvilinear coordinates, whereas is a suitably smooth function. [Sometimes it is more convenient to use the symbols (); other times we prefer to use . Accordingly, we use 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 , with , where stands for either + or −, is understood as a sum that spans over the eight values of , which correspond to the eight vertices of the brick. Furthermore, the symbol denotes the partial derivative with respect to : . Similarly, we have and . 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 . Moreover, the term 123 is the only mixed third-order derivative. Finally, , , , and are suitable products of the Hermite polynomials in (4). For instance,
In summary, for the three-dimensional third-order Hermite interpolation scheme, the finite-element unknowns are the nodal values of () the unknown function, () its three first-order partial derivatives, () its three second-order mixed derivatives, and () its third-order mixed derivative . 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 third-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 where the unknowns comprise the values of and of all the partial derivatives mentioned above, evaluated at the 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 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 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 . 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 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–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 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 (single-brick 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–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 three-dimensional 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 , with , 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 and are assumed to be prescribed. Then, the algorithm generates a surface that has these lines as edges. Specifically, the Coons patch is obtained as follows: () consider the sum of the two linear interpolations between the two sets of opposite boundary lines; () from this subtract a bilinear interpolation through the four corner points. This yields (again with ) The four edges of this surface indeed coincide with the four generating lines. For instance, we have (use and ) because and (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 , which describes the brick, is obtained as () the sum of the three linear interpolations between opposite faces, minus () the sum of three bilinear interpolations through three sets of four “parallel” edges, plus () a trilinear interpolation through the eight vertices. Stated in mathematical terms, we havewhere , , and describe the geometries of the six faces, which for the time being are assumed to be prescribed. Similarly, denotes the edge for , and denotes the vertex for . [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 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 where are the four corner points of the patch, and 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 (24) and (27) yieldswhere the functions , , and are given by (2) and (3). Introducing the functions where and stand for either + or −, (28) may be rewritten as
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 , , and , one obtains 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 byand 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 , where , with , is understood as a sum that spans over the eight values of , 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 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 within the block, in terms of the values of and of its partial derivatives , 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 (isoparametric formulation). Accordingly, the unknowns consist solely of the values of and of the three Cartesian components of 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.
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 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 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)). For each block, we replace the two-node Lagrange and Hermite interpolations, with the -node ones (see (5) and (9)). Then, we can extend (35), to have where whereas [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 . 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.
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 (Section 7.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 , which is such that . Setting , we obtain the Helmholtz equation; namely, where 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 , where on we have (a condition often used for the opening of wind and brass musical instruments), whereas is a rigid wall. On , using (43), the boundary condition is given by On the other hand, on , 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, , (44)), which has to be satisfied by all the shape functions. The second type is the natural boundary condition (e.g., Neumann boundary condition, , (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 () are the global shape functions, defined within each brick/block by the local interpolation functions discussed in the preceding sections, whereas () comprises all the corresponding unknowns. [To be specific, include the nodal values of and .]
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 (see (47)), with (46), to obtain where is the vector of the unknowns and , whereas the mass and stiffness matrices are, respectively, given by and , 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 least-square 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 satisfies the orthonormality relation , which is fully equivalent to the condition , 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 defined as follows: the components are the nodal values corresponding to but are obtained from the exact th eigenfunction and are normalized by . Combining (51)–(53) and using (49) as well as , 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,
Here, we consider first a cube of side (namely, ), so as to have .
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 finite elements (blocks) in each direction, with each block being composed of two bricks in each direction, for a total of 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 and , respectively, which are equal, respectively, to and . 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 log-log scale, with 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 bricks. Again, the order of the polynomials used in the scheme is . Figures 4 and 5 present the eigenvalues and , as functions of the order . Similar results hold for higher eigenvalues.
Note that, in the above test case we have, for instance, (degenerate eigenvalue). This phenomenon makes it difficult to compare the numerical eigenfunctions to the analytical ones (see (57)), because any linear combination of , , and 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 , , and .
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 elements (namely, third, fifth, and seventh orders, resp.). These correspond to , 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 () only one eigenvalue when ; () groups of three approximately equal eigenvalues when only two of them are equal (they would be identical if we had ); and () groups of six approximately equal eigenvalues when , , and differ.
Figures 6, 8, and 10 estimate the projections of the numerical eigenfunction vectors , for into each of the exact eigenfunction vectors , for , through , a generalization of (55). The darker the mark, the poorer the agreement. On the other hand, Figures 7, 9, and 11 depict for . 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 ).
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 and length (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 .
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 and , that is, those eigenvalues with , and . Figures 15 and 16 pertain to the corresponding -convergence.
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 (Figure 17). For each block, use bricks in each direction, for a total of bricks for each of the three blocks.
We indicate with the th eigenvalue of the discretized problem.
In Figures 18-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.
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 where denotes the density, , and , 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, is the vector of the unknown nodal values. The mass and stiffness matrices are, respectively, given by and , with where ( 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 , and rectangular cross section ( and ) 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 and and the solution is given by where is the solution of the characteristic equation
Then, one obtains the eigenvalue 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 -node element with three degrees of freedom per node, i.e., the three components of the displacement. The latter is an -node element with six degrees of freedom per node, i.e., the three components of the displacement along with the three rotations.]
Hinged Rectangular Plate. Consider the results relative to a rectangular hinged plate with , , and . The eigenfunctions for a plate in bending areThe 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 three-dimensional modes .]
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 and , respectively, using, again, the third-order element with a single element along the thickness, for a total of bricks. The computed eigenvalues are compared to the exact ones and with those obtained by using the -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 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 pertains to the calculated eigenvalues obtained using the third-order element with . The considerations presented for acoustics, regarding the number of accurate frequencies, apply to the present case as well. About one-half of the first modes are correctly captured. [Note that the total number of degrees of freedom is . However, the last modes are related to the rotation of the normal and the curvature of the fiber along the thickness.]
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 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 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 (see Appendix). The geometry of each block is described by the nodal points , along with the corresponding base vectors , where (), so as to have 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 ). The overall order of the scheme (namely, the degree of the polynomials used) is . 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 out-of-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 () air inside a cavity (interior acoustics) and () 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].
Appendix
High-Order Geometry Generation, from a First-Order One
As mentioned above, in the formulation presented in the main body of this work, we assume that the geometry consists of a collection of high-order topologically quadrilateral blocks, which requires a geometry preprocessor that provides the location of the nodes and the corresponding base vectors (see (19)). In this appendix we present a user-friendly algorithm that allows one to generate such a geometry, starting from a very simple one, which is available in any geometry preprocessor. Specifically, we assume that we know only the locations . Here, we show how to obtain the base vectors at the nodes of each brick and then those for a whole block.
The Third-Order Surface Generation. In order to arrive at the algorithm for a third-order brick, we begin by considering the face of a brick. For the sake of notational simplicity, the surface coordinates are denoted by and .
Consider the mapping , which denotes the geometry obtained by a bilinear interpolation of the four nodes (first-order geometry). The corresponding surface base vectors are given by and its normal by
In this subsection, we assume that we know the location of the four nodes of the surface, as well as the normal at the four corners of the patch. [How to obtain the normal at the four corners is the subject matter of the next subsection.] Specifically, we seek a surface described in the formwhere vanishes at the four corners, because we have chosen .
Such a surface is to be determined so as to have the prescribed normal at the four corners. To this end, we use the Coons–Hermite interpolation (see (28)), applied to the function and obtain (see Figure 30, for symbols definition)[The first four terms in (28) are not included here, because of (A.4).]
Consider the first term of the four terms on the right side of (A.5). This term vanishes () for , because , and () for , because . In other words, this term vanishes along three of the edges, whereas on the fourth edge (namely, for ) (A.5) simply provides the Hermite interpolation for the edge . Similar considerations hold for the other three contributions; each one of them affects the function only along one of the four edges.
Accordingly, the third-order surface is fully described by the two first-order partial derivatives of , evaluated at the four corners. How do we obtain these? By imposing that the normal is orthogonal to the base vectors. These are obtained by differentiating (A.3) with respect to and respectively, to obtain Recalling that (see (A.4)), the orthogonality condition yields, at the four corners which allows one to obtain and at the four corners.
Finally, the third-order surface is obtained by combining (A.3) with the expression for in (A.5), with and at the four corners obtained from (A.7). Then, the surface base vectors at the four corners are given by because (see (A.4)).
The fact that each edge is affected only by the derivative along the edge at its two-end points is noteworthy. This implies that the common edge of two adjacent surfaces is described by the same function. This guarantees continuity of the surface.
Third-Order Geometry. Let us now turn to a block. Let us consider a coordinate surfaces for the whole patch, say a surface. From (A.2), at each node we obtain four values for , namely, one value for each of the four first-order patches involved. The prescribed unit normal (as used in (A.7)) is obtained by taking the average of these four normals and normalizing it so as to have a unit vector.
Then, we can repeat the same formulation for each of the eight faces of the brick and obtain the three base vectors at each of its eight vertices. This completes the generation for the third-order geometry.
The Geometry of Order . Finally, note that once the three base vectors are available for each node of the patch, one may use the high-order element family proposed in Section 6. Therefore, the question arises: “Can the geometry generated using the formulation above be used for the finite-element formulation of order ?” “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 used to obtain the nodes for the first-order geometry is smooth. In this case, base vectors that are continuous at the nodes are obtained by taking the average of the two nodal values of . The final result is a block 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.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.