Next Article in Journal
Using the Theory of Functional Connections to Solve Boundary Value Geodesic Problems
Next Article in Special Issue
Modeling of Perforated Piezoelectric Plates
Previous Article in Journal
The Binomial–Natural Discrete Lindley Distribution: Properties and Application to Count Data
Previous Article in Special Issue
An Efficient Orthogonal Polynomial Method for Auxetic Structure Analysis with Epistemic Uncertainties
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Use of High-Order Shape Functions in the SAFE Method and Their Performance in Wave Propagation Problems

1
Department of Mechanical Engineering, The Ferdowsi University of Mashhad, Mashhad 9177948974, Khorasan Razavi Province, Iran
2
Institute of Mechanics, Otto von Guericke University Magdeburg, 39106 Magdeburg, Sachsen-Anhalt, Germany
3
The Department of Maritime and Transport Technology, Delft University of Technology, P.O. Box 5, 2600 AA Delft, Zuid-Holland, The Netherlands
4
School of Civil and Environmental Engineering, The University of New South Wales, Sydney, NSW 2052, Australia
*
Author to whom correspondence should be addressed.
Math. Comput. Appl. 2022, 27(4), 63; https://doi.org/10.3390/mca27040063
Submission received: 2 June 2022 / Revised: 10 July 2022 / Accepted: 13 July 2022 / Published: 25 July 2022

Abstract

:
In this research, high-order shape functions commonly used in different finite element implementations are investigated with a special focus on their applicability in the semi-analytical finite element (SAFE) method being applied to wave propagation problems. Hierarchical shape functions (p-version of the finite element method), Lagrange polynomials defined over non-equidistant nodes (spectral element method), and non-uniform rational B-splines (isogeometric analysis) are implemented in an in-house SAFE code, along with different refinement strategies such as h-, p-, and k-refinement. Since the numerical analysis of wave propagation is computationally quite challenging, high-order shape functions and local mesh refinement techniques are required to increase the accuracy of the solution, while at the same time decreasing the computational costs. The obtained results reveal that employing a suitable high-order basis in combination with one of the mentioned mesh refinement techniques has a notable effect on the performance of the SAFE method. This point becomes especially beneficial when dealing with applications in the areas of structural health monitoring or material property identification, where a model problem has to be solved repeatedly.

1. Introduction

The semi-analytical finite element (SAFE) method is an efficient tool for the analysis of wave propagation problems and the computation of dispersion curves. As the name suggests, the equations of motions are presented analytically in the direction of wave propagation, while the cross-section of the structure under investigation is discretized in a finite element sense [1]. Therefore, SAFE is often used in structures with non-varying cross-section such as pipes [1,2] and plates [1,3]. As a result of its formulation, the dimension of the spatial discretization is reduced by one or two and, therefore, the number of degrees of freedom (DOF) is drastically reduced in comparison to established numerical methods such as the finite element method (FEM), the spectral element method (SEM), or isogeometric analysis (IGA). Consequently, significant savings in terms of computational costs can be achieved by applying SAFE to suitable classes of problems. Note that, to ensure accurate numerical results, researchers have proposed to use at least 6 to 10 nodes [4,5,6] per wavelength if a piecewise linear approximation of the displacement field is employed. This is a rather rough estimate and according to the findings discussed in Refs. [7,8] more than 20 nodes per wavelength are required for accurate results. For a detailed guide on how to set up the spatial discretization for wave propagation problems including the analysis of high-order methods, the reader is referred to [8].
On the other hand, in comparison with purely analytical approaches, semi-analytical methods such as SAFE exhibit more flexibility in modeling structures with a higher degree of complexity. The existing analytical methods for modeling wave propagation in the structures are mainly based on either the transfer matrix method or the global matrix method [9,10,11,12,13]. These approaches are essentially limited to structures with simple geometries and low frequency excitation regimes. In addition, they typically require a numerical search algorithm for finding the roots of the analytical formulations [11,14].
Still, it must be noted that semi-analytical methods do not only inherit the best aspects of both worlds, but also some of their limitations. One of the major disadvantages of these kinds of methods is seen in the limitations related to the geometry of the structure. Here, the cross-section of the waveguide is required to be constant along the wave propagation direction. In the following, a comprehensive and unbiased discussion on recent advances regarding the extension of SAFE to more complicated structures is attempted. Additionally, the general benefits and drawbacks of this method are outlined.
The SAFE method for wave propagation problems in waveguides of arbitrary cross-section was first presented by Aalami [15] and Lagasse [16] in 1973. This method was previously used in structural engineering [17] and is closely related to both the strip and thin element methods [14]. Since its development, the SAFE method has been extended by other researchers to account for more complex structures. Examples including, but not limited to are:
  • plates with viscoelastic damping [1,3];
  • cylinders [3], thin-walled steel pipes, aluminum profiles, orthotropic tubular profiles, and angle steel beams [18];
  • wedges with different shapes [19];
  • rails [20,21,22];
  • pre-twisted rods and curved waveguides (e.g., rings with rectangular cross-section) [19];
  • viscoelastic axisymmetric waveguides [1,3];
  • springs, multi-wires [23];
  • finned tubes [24];
  • immersed or embedded (e.g., in soil or concrete) waveguides with various cross-sections [25,26,27,28,29,30,31,32].
Besides research devoted to extending SAFE in terms of its geometry approximation capabilities, another important aspect is the investigation of different approaches to spatially discretize the cross-section of the structure. Depending on the complexity of the cross-section, the SAFE method can be divided into one-dimensional (1D SAFE) and two-dimensional (2D SAFE) approaches where corresponding finite elements are used for the spatial discretization. In structures with simple geometries such as infinite plates and axisymmetric waveguides, in order to reduce the computational effort, the 1D SAFE method may be used [3]. In this case, three-noded elements with quadratic shape functions and three displacement DOF at each node are the most popular choice [1,2,11,12,18,27,30,33,34,35,36,37,38,39], while two-noded elements [40] and three-noded elements with two displacement DOFs at each node [1] are also considered in the literature. In the case of the 2D SAFE method, the elements mostly used are linear three-noded triangular and bi-linear four-noded quadrilateral elements [1,11,19,21,22,23,24,29,32,41,42,43,44,45,46,47,48,49,50]. In some cases, quadratic six-noded triangular elements, bi-quadratic eight-noded quadrilateral (Serendipity) elements [18,25,26,28,49,51,52,53,54], and nine-noded Lagrange quadrilateral elements [26] are also proposed. Steshedi et al. [49] compared three- and six-noded triangular with four- and eight-noded quadrilateral elements for the analysis of rails. They considered the convergence of the group velocity at 40 kHz for different element topologies to find the best spatial discretization. They recommended the eight-noded quadrilateral element due to its accuracy and convergence properties in comparison to the others. However, no data concerning the convergence rate or the attained accuracy for the different elements were presented.
Despite the advantages of using high-order shape function, in the wide body of literature concerning SAFE, the vast majority of research is related to low-order shape functions and only very few high-order implementations have been discussed. The advantages of using high-order shape functions over low order ones are seen in the following properties: (i) exponential rates of convergence [55], (ii) robustness with respect to mesh distortion, and (iii) negligible locking [55,56]. Kalkowski et al. [57] used high-order spectral elements to model an axisymmetric embedded waveguide using 1D elements along with perfectly matched layers (PML) for modeling the surrounding (infinite) medium. PMLs are of great use for modelling waves in structures embedded in a restraining medium. PML maps the original geometry into a new coordinate system which both ‘stretches’ the original geometry and attenuates the waves (as the new coordinate is complex) [57]. In their implementation, Kalkowski et al. used Lagrange polynomials over Gauss–Lobatto–Legendre (GLL) points. Xiao et al. [58] employed 1D elements to model guided wave propagation in an infinite functionally graded magneto-electro-elastic plate by the Chebyshev spectral element method. Here, the shape functions are Lagrange polynomials as well, but defined over Gauss–Lobatto–Chebyshev (GLC) points. Treyssede et al. [59] discretized the geometry with 2D spectral elements, and computed high-frequency low-leakage modes in structural waveguides surrounded by infinite solid media. Recently, Seyfaddini et al. utilized non-uniform rational B-splines (NURBS) basis functions (well-known from IGA) in the SAFE method for modeling wave propagation in isotropic plates [60], functionally graded plates immersed in fluids [61], and poroelastic plates [62]. In their study, they have compared the obtained results to a different SAFE approach utilizing high-order shape functions based on Lagrangian polynomials and equidistant nodes. According to their study, the semi-analytical approach based on IGA provides an improved performance.
According to the studies available in the literature, the SAFE method is an efficient method for structural health monitoring (SHM) and material property identification purposes. In material property identification problems [36,49], first, dispersion curves are derived using experimental methods. Thereafter, using model updating approaches, the material parameters are adjusted in the numerical model using an optimization strategy such that the computed dispersion curves agree with the experimental ones. Consequently, a large number of simulations is required. Since both the efficiency and accuracy of the solution procedure for such problems is of utmost importance, it is necessary to thoroughly investigate different numerical methods which can be straightforwardly combined with SAFE. To the authors’ knowledge, no research has been conducted yet, where the advantages of using various high-order shape functions and a thorough comparison in terms of their efficiency and computational costs have been attempted. Therefore, the overarching goal of this contribution is to implement high-order shape functions known from p-FEM, SEM, and IGA into SAFE. Furthermore, the ramifications of using different refinement techniques, e.g., h-, p-, and k-refinement, are investigated in detail.
The remaining paper is organized as follows: In Section 2, the theoretical derivation of the SAFE method is presented. Then, high-order shape functions are briefly discussed in Section 3. The numerical results for isotropic and composite plates are reported in Section 4, while important conclusions are drawn in Section 5.

2. SAFE Method

The SAFE method is a combination of an FE-based spatial discretization of the cross-section of the waveguide and an analytical solution in the direction of wave propagation. To this end, the analytical term e i ( ω t κ x ) is introduced in the solution, where κ represents the wavenumber, x denotes the direction of wave propagation, ω is the circular frequency, and t stands for time. In this method, the cross-section of the waveguide, depending on the complexity of the structure, is discretized employing 1D or 2D finite elements in order to achieve a three-dimensional model of the structure. Consequently, in comparison to numerical methods that fully discretize the computational domain, the generated models exhibit a significant reduction in the number of DOFs and the numerical costs required for performing the analysis. In order to emphasise the differences and similarities between the SAFE method and FEM, typical elements, nodal distributions, and shape functions are compared with each other. However, it should be noted that the derived stiffness and mass matrices are different from each other which will be addressed later. Typical 1D, 2D, and 3D finite elements with their corresponding DOFs are depicted in Figure 1.
For the case of the 3D element (see Figure 1c) with three displacement DOFs at each node
N l = N l ( x , y , z ) , l = 1 , 2 , , β ,
the displacement field at each point in the element may be expressed as:
u ( e ) ( x , y , z ) = u x ( e ) ( x , y , z ) u y ( e ) ( x , y , z ) u z ( e ) ( x , y , z ) = l = 1 β N l ( x , y , z ) u l x l = 1 β N l ( x , y , z ) u l y l = 1 β N l ( x , y , z ) u l z = N ( x , y , z ) q ( e )
with
N = N 1 0 0 N 2 0 0 N β 0 0 0 N 1 0 0 N 2 0 0 N β 0 0 0 N 1 0 0 N 2 0 0 N β
and
q ( e ) = [ u 1 x u 1 y u 1 z u 2 x u 2 y u 2 z u 3 x u 3 y u 3 z u β x u β y u β z ] T .
Here, β denotes the number of nodes per element, and u i j refers to the displacement of node number i in the direction j and the superscript ( e ) generally refers to an elemental property. As can be easily seen, the shape functions are functions of the spatial coordinates x, y, and z and the DOFs considered for each node are also available in the same coordinate directions.
For the cases of 1D and 2D elements (see Figure 1a,b, respectively), the shape functions are only functions of the spatial coordinate z
N l = N l ( z ) , l = 1 , 2 , , β ,
or the spatial coordinates y and z
N l = N l ( y , z ) , l = 1 , 2 , , β .
Despite the analogy with standard FEM, the displacement DOFs in the SAFE method are defined in a slightly different way. In the SAFE method, independent of whether we employ the 1D or 2D variant, all three displacement DOFs are available as we want to model a three-dimensional structure. That is to say, the shape function matrix N and the elemental displacement vector q ( e ) are formulated as in Equations (3) and (4), respectively. Note, however, that, for the 1D version of SAFE, the one-dimensional shape functions N l ( z ) according to Equation (5) are also used to populate the shape functions matrix, while in the 2D SAFE case the two-dimensional shape functions provided in Equation (6) are implemented. Thus, the displacement of each point in an element in the SAFE method is given as follows:
  • 1D SAFE method:
    u ( e ) ( x , z , t ) = u x ( e ) u y ( e ) u z ( e ) = l = 1 β N l ( z ) u l x l = 1 β N l ( z ) u l y l = 1 β N l ( z ) u l z e i ( ω t k x ) = N ( z ) q ( e ) e i ( ω t k x )
    with
    N = N 1 0 0 N 2 0 0 N β 0 0 0 N 1 0 0 N 2 0 0 N β 0 0 0 N 1 0 0 N 2 0 0 N β
    and
    q ( e ) = [ u 1 x u 1 y u 1 z u 2 x u 2 y u 2 z u 3 x u 3 y u 3 z u β x u β y u β z ] T
  • 2D SAFE method:
    u ( e ) ( x , y , z , t ) = u x ( e ) u y ( e ) u z ( e ) = l = 1 β N l ( y , z ) u l x l = 1 β N l ( y , z ) u l y l = 1 β N l ( y , z ) u l z e i ( ω t k x ) = N ( y , z ) q ( e ) e i ( ω t k x )
    with
    N = N 1 0 0 N 2 0 0 N 3 0 0 N β 0 0 0 N 1 0 0 N 2 0 0 N 3 0 0 N β 0 0 0 N 1 0 0 N 2 0 0 N 3 0 0 N β
    and
    q ( e ) = [ u 1 x u 1 y u 1 z u 2 x u 2 y u 2 z u 3 x u 3 y u 3 z u β x u β y u β z ] T
This conceptual difference results in the need to adapt the conventional FE-shape function matrix for the use in the SAFE method. According to the aforementioned explanations about the similarities and differences between the SAFE method and FEM, the formulation process for the case of the 1D SAFE method is briefly explained in the following. It should also be kept in mind that a three-dimensional model is generated in the 1D SAFE method in contrast to 1D FEM. The linear strain-displacement relation is presented as [1]:
ε ( e ) = L x x + L y y + L z z u ( e ) = ( β 1 + i k β 2 ) q ( e ) e i ( ω t k x )
with the matrices
L x = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 , L y = 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 , L z = 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0
and derived quantities
β 1 = L y N y + L z N z , β 2 = L x N .
The stress-strain relation is given as:
σ ( e ) = C ( e ) ε ( e ) ,
where C is the material’s elasticity matrix. Assigning the kinetic and potential energy to the variables K and Φ , respectively, the variation of the Hamiltonian may be written as:
δ H = t 1 t 2 δ ( Φ K ) d t = 0 ,
with
Φ = 1 2 V ε T C ε d V ,
K = 1 2 V u ˙ T ρ u ˙ d V .
Here, u ˙ refers to the derivative of displacement with respect to time and ρ denotes to the mass density of the material. By replacing the strain and displacement relations in Hamilton’s principle and following the same process as in the FEM to derive the weak form, the stiffness matrices ( k 1 ( e ) , k 2 ( e ) , and k 3 ( e ) ) and the mass matrix ( m ( e ) ) are derived [1]:
k 1 ( e ) = Ω e β 1 T C ( e ) β 1 d Ω e ,
k 2 ( e ) = Ω e β 1 T C ( e ) β 2 β 2 T C ( e ) β 1 d Ω e ,
k 3 ( e ) = Ω e β 2 T C ( e ) β 2 d Ω e ,
m ( e ) = Ω e N T ρ ( e ) N d Ω e .
As is evident from the above equations, unlike FEM, we obtain three stiffness matrices in the SAFE method. The matrix k 1 ( e ) relates to deformations in the cross section of the structure, which is located in the y z plane. The matrix k 2 ( e ) matrix couples the in-plane ( y z plane) deformations to out-of-plane (x axis) deformations. Finally, the matrix k 3 ( e ) represents out-of-plane deformations. The stiffness and mass matrices of each element, for the total number of n e l elements, are assembled in order to obtain the global stiffness matrices ( K 1 , K 2 , and K 3 ) and the global mass matrix ( M ) of the entire system:
K 1 = e = 1 n e l k 1 ( e ) , K 2 = e = 1 n e l k 2 ( e ) , K 3 = e = 1 n e l k 3 ( e ) , M = e = 1 n e l m ( e ) .
the corresponding eigenvalue problem, which needs to be solved to obtain the solution is given as:
( K 1 + i κ K 2 + κ 2 K 3 ω 2 M ) Q = 0 ,
where Q refers to the global vector of displacements at each node. Employing a change of variable, we may rewrite the second-order eigenvalue problem, cf. Equation (25), into a first-order one:
A κ B Q ^ = 0 .
To this end, the following matrices have been introduced [1]:
A = 0 K 1 ω 2 M K 1 ω 2 M K ^ 2 ,
B = K 1 ω 2 M 0 0 K 3 ,
Q ^ = Q κ Q ,
K ^ 2 = i T T K 2 T ,
T = i 1 0 1 i 0 1 1 .

3. High-Order Shape Functions

In order to increase the accuracy of classical numerical methods, either h-refinement or p-refinement or a combination of both may be used. The basic idea of h-refinement is to successively sub-divide the computational domain into smaller and smaller sub-domains, also known as elements, while the order of the polynomial shape functions is kept constant. On the other hand, p-refinement refers to an approach, where the number of elements is kept constant, and meanwhile the order of the polynomial shape functions is increased. Besides p- and h-refinement, the particular choice of shape functions is also of great importance. Typical examples of high-order shape functions which can be considered are:
  • Standard shape functions based on Lagrange polynomials with equidistant nodal distribution:
    • The Lagrange interpolation polynomials are defined as [63,64]
      L p ( ξ ) = k = 1 , k p p + 1 ξ ξ k ξ p ξ k .
      L p ( ξ ) is a p-th degree polynomial given by the product of p linear factors [65]. Moreover, ξ k represents the given set of points over a specified range, while ξ p is a particular point in that set. An equidistant distribution of the points over the interval [ 1 , 1 ] , which is a typical for a local coordinate system of an element, is calculated by [66]:
      ξ k = 1 + 2 k 1 p , k = 1 , 2 , , p + 1 .
  • Spectral shape functions based on Lagrange polynomials with Gauss–Lobatto–Legendre (GLL) and Gauss–Lobatto–Chebyshev (GLC) nodal distributions:
    • Standard shape functions based on Lagrange polynomials defined using an equidistant nodal distribution lead to large oscillations and consequently to ill-conditioning [63,67]. This poor behavior of polynomial interpolation in equally spaced points is known as Runge’s phenomenon, which infers that increasing the order of shape functions does not necessarily incur lower errors [63,68]. Therefore, in order to circumvent this problem, non-equidistant nodal distributions such as GLL points [69] and GLC [70] have been proposed. In Figure 2a–c, the shape functions of a fifth order Lagrange polynomial based on equidistant, GLL, and GLC nodal distributions are depicted, respectively. The GLL and GLC nodal distributions along with high-order Lagrange polynomials are used in the SEM [63], which was first presented by Patera et al. in the field of fluid dynamics [71]. Note that SEM with GLL points has a special characteristic besides low numerical dispersion. The mass matrix for this particular spectral element is approximated through reduced integration by a diagonal matrix [67,72]. However, the SEM suffers from the fact that, for increasing the order of the shape functions, an additional node should be inserted and consequently all other internal nodes are repositioned and the coefficient matrices must be re-computed [73,74]. To derive the GLL points, we choose the roots of the first derivative of a p-th degree Legendre polynomial plus the endpoints located at ξ 1 = 1 and ξ p + 1 = 1 [63]
      ξ k = { 1 k = 0 ξ k 2 k p + 1 k = p + 1 ,
      where ξ k denotes the aforementioned roots.
  • Hierarchic shape functions based on the normalized integrals of the Legendre-polynomials:
    • Hierarchic shape functions are mainly used in the p-version of FEM. The principal difference between hierarchical shape functions and other sets of shape functions lies in the fact that in the hierarchic case the new DOFs are added to the existing ones without changing the existing DOFs. Thus, all lower order shape functions are contained in the higher order basis as depicted in Figure 2d. Unlike the previously discussed shape functions, where all the shape functions are of the same order, the order of the hierarchic shape functions increases incrementally starting from the first two linear shape functions. The one-dimensional, hierarchic shape functions consist of two nodal shape functions and various internal shape functions
      N 1 ( ξ ) = 1 2 ( 1 ξ ) ,
      N 2 ( ξ ) = 1 2 ( 1 + ξ ) ,
      N n ( ξ ) = 1 l e n 2 x = 1 x = ξ l e n 2 ( x ) d x , n = 3 , 4 , 5 , , p ,
      which are normalized integrals of Legendre polynomials. Here, l e n 2 ( x ) is a Legendre polynomial of order n 2 . Note that, in p-FEM, not all DOFs represent actual nodal displacements [72], instead they are actually only unknowns of the high-order ansatz space [73].
  • Non-Uniform Rational Basis Splines:
    • In standard isoparamteric approaches, the shape functions used on an elemental level to discretize the geometry under investigation are only approximate in nature. Hence, numerical errors in geometries description are cause, which is especially important for geometries including curved surfaces. Thus, Hughes et al. [75] introduced the socalled isogeometric analysis in which nonuniform rational B-splines (NURBS) are used as shape functions. NURBS are frequently used for describing the geometry in a computer aided design (CAD) and computer aided manufacturing (CAM) applications [75,76]. NURBS are a generalization of B-splines and are derived by projecting B-splines from R d + 1 into R d , where d represents the dimensionality of the problem. A NURBS entity can be defined using knot vectors and a series of control points. The term ‘knot vector’ can be deceiving as a knot vector is simply a list of non-decreasing entries called knots and is not a real vector. A knot vector (see the example below),
      Ξ = ξ 1 ξ 2 ξ n + p + 1 = 0 0 0 1 / 3 2 / 3 1 1 1 ,
      defines the parametric space of a univariate NURBS entity (curve). In Equation (38), p and n represent the degree of the NURBS and the number of control points, respectively. A knot vector can range between any two real numbers, but it is very common to choose 0 and 1 as the lower- and upper-bounds of the parametric space. The space between two adjacent knots is called a knot span. Non-zero knot spans can be compared to elements in the classical finite element method. Knot vectors control the discretization of the geometry and also the continuity across the elements, whereas the control points are responsible for all the changes regarding the geometry. As a rule, a degree p NURBS has C p m -continuity across the knot spans inside a patch, where m is the knot multiplicity in the knot vector. (Remark: NURBS patches are geometrically simple domains within which the element types and material models are presumed to be uniform.) The first and the last knots in the knot vectors of ‘open NURBS’ which are considered in this study have a multiplicity of p + 1 , making their ends C 1 -continuous. For instance, the knot vector in Equation (38) represents a quadratic ( p = 2 ) univariate NURBS with two elements, where the continuity on knots 1 / 3 , and 2 / 3 is C 2 1 = C 1 . A univariate NURBS is a rational generalization of a B-spline curve. For a degree p NURBS curve, the corresponding shape function is given by (see [77,78])
      R i ( ξ ) = N i , p ( ξ ) w i W ( ξ ) = N i , p ( ξ ) w i i ^ = 1 N i ^ , p ( ξ ) w i ^ ,
      where W ( ξ ) is the weighting function, w i represents the weights associated with the control points, and, finally, N i , p ( ξ ) denotes the standard B-spline basis functions in the parametric space. Shape functions of basis splines are defined by the Cox–de Boor recursion formula:
      For p = 0 ,
      N i , 0 ( ξ ) = 1 , if ξ i ξ < ξ i + 1 0 , otherwise .
      For p = 1 , 2 , 3 , ,
      N i , p ( ξ ) = ξ ξ i ξ i + p ξ i N i , p 1 ( ξ ) + ξ i + p + 1 ξ ξ i + p + 1 ξ i + 1 N i + 1 , p 1 ( ξ ) .
      Refinement of NURBS entities is usually achieved by modifying the underlying B-spline (in R d + 1 ). There are three methods to refine a B-spline/NURBS entity:
      -
      Knot insertion: Additional knots are added to the knot vector, increasing the number of control points and basis functions. Knot insertion is closely related to h-refinement in the classical FEM, as inserting additional knots increases the number of finite elements in the resulting discretization. Each knot must be inserted p times, maintaining the C 0 -continuity between the elements, to resemble the classical h-refinement technique completely,
      -
      Degree elevation: The polynomial degree of the NURBS entity is increased, while preserving the continuity of the curve. In the case where the continuity between the elements is C 0 , degree elevation is equivalent to the p-refinement technique in the classical FEM.
      -
      k-refinement: Combination of knot insertion and degree elevation techniques. To better understand the k-refinement technique, one ought to note that knot insertion and degree elevation are not commutative. In other words, the order in which the refinement techniques are applied plays a crucial role in achieving the desired properties of the final discretization. k-refinement makes it possible to increase the polynomial degree of the NURBS object, while increasing its continuity as well. This leads to more accurate results using less control points. To achieve the highest level of its benefits, i.e., C p 1 continuity of the basis between the elements of the refined discretization, one needs to start with one element in the coarsest discretization, elevate the basis to the desired degree and then insert the additional knots to introduce the elements. On the other hand, inserting knots first on a lower degree and then elevating the degree will lead to lower continuity across the elements ( C 0 , if one starts from a linear discretization) and more control points. This technique does not have any equivalents in the classical FEM.
      It is important to note that none of these refinement techniques changes the original geometry of the NURBS object. For a more detailed introduction to NURBS and k-refinement, the interested reader is referred to the seminal monographs by Piegl and Tiller [78] and Cottrell et al. [77].

4. Results and Discussion

In order to investigate the effect of using the aforementioned shape functions and refinement techniques in the SAFE method, two case studies, i.e., an isotropic and a composite plate (consisting of eight layers), are considered. The efficiency and accuracy of various combinations of the introduced approaches are assessed through considering models with the same number of DOFs and consequently, eigenvalue problems with the same dimensions. From the obtained results, it is clear that the method attaining the highest accuracy should be favoured in future applications as a considerable number of DOFs can potentially be saved. The following combinations of shape functions and refinement techniques are investigated in the remainder of this article:
  • h-refinement approach using:
    • second order standard shape functions based on Lagrange polynomials with equidistant nodal distribution (FEM);
    • second order NURBS shape functions with C 1 -continuity (NURBS-based IGA);
  • p-refinement approach using high-order shape functions ( p > 2 ):
    • standard shape functions based on Lagrange polynomials with equidistant nodal distribution;
    • spectral shape functions based on Lagrange polynomials with Gauss-Lobatto-Legendre nodal distribution (SEM);
    • spectral shape functions based on Lagrange polynomials with Gauss-Lobatto-Chebyshev nodal distribution (SEM);
    • hierarchic shape functions based on the normalized integrals of the Legendre-polynomials (p-FEM);
    • NURBS with C 0 -continuity (NURBS-based IGA)
  • k-refinement approach using:
    • NURBS with C p 1 -continuity (NURBS-based IGA).
In order to calculate the attainable error, the SAFE results are compared to the analytical results derived using the Dispersion Calculator software, which is a free software provided by the German Aerospace Center (DLR) for calculating dispersion curves of guided waves in isotropic plates and multilayered composites [79]. Furthermore, for calculating the error, the following relation is used:
E r % = k ref k SAFE k ref × 100 ,
where k ref is the wavenumber derived using the analytical method and k SAFE is the wavenumber derived using the SAFE method. In order to calculate the k SAFE , real values were assigned to the circular frequency ω and the first-order eigenvalue problem, cf. Equation (26), was solved. More details about different approaches to deal with the eigenvalue problem and calculating phase velocity ( C p ), group velocity ( C g ) and cut-off frequencies are presented in Appendix A. The most important issue in dealing with the derived eigenvalues is correlating the eigenvalues to the corresponding wave modes. In order to fulfill this task, either a method based on the orthogonality property of modes or different mode tracking methods based on Taylor or Padé series expansions can be used. The details concerning various mode identification approaches are further discussed in Appendix B. In this paper, we have considered a Padé expansion of order [1/2] for extrapolating the wavenumber curves along with a Padé expansion of order [0/1] for extrapolating the phase velocity curves as suggested by Gravenkamp et al. [80]. It is worth mentioning that in the mode identification approach suggested by Gravenkamp et al. [80], k = 1 / C g is replaced in the Padé expansion for extrapolating the group velocity curve (see Equations (A13)–(A15) in Appendix B for further details) and the second and third derivatives of C g are considered for estimating the wavenumber value. This replacement has led to using the group velocity values which are already considered separately. Hence, to take full advantage of the proposed approach, we have considered the actual group velocities instead of replacing the value with C g . After deriving the predicted wavenumber using the extrapolation, the error corresponding to the wavenumber will be
Δ k = k p k j k p ,
where k p is the guessed wavenumber and k j is the wavenumber calculated in each step. Furthermore, in order to increase the robustness of this method, the phase velocity curves may also be utilized besides the wavenumber curves. Thus, in each step, the phase velocity is guessed using the previous phase velocities corresponding to the same mode, and the guessed value is compared to the the derived phase velocities as follows:
Δ C g = C g p C g j C g p ,
where C g p is the guessed phase velocity and C g j is the phase velocity calculated in each step and an error function is constructed including both wavenumbers and phase velocities as presented in Equation (45). The minimum error value may lead us to the wavenumber representing the same wave mode
Δ = Δ k 2 + Δ C g 2 .

4.1. Isotropic Case Study

The first case model problem is an aluminum plate with 1 mm thickness. The material properties of the plate are listed in Table 1.
In order to model wave propagation in the plate only one-dimensional, elements with two DOFs at each node (FEM) or control point (IGA) are taken into account. The DOFs considered herein are along the wave propagation direction (x) and along the plates thickness (z) as depicted in Figure 3. Note that the third DOF along the width of plate (y) can be avoided in infinite isotropic plates, since shear and Lamb waves can be easily decoupled. (Remark: A two-dimensional model is investigated in this example, where one of the in-plane directions is assumed as infinite (plane strain assumption). Therefore, only two DOFs per node/control point are used.) Consequently, we are only considering the existence of Lamb waves. As the first step, in order to illustrate the accuracy of the simulations performed (see Figure 4), the derived dispersion curves using the SAFE method and Dispersion Calculator software are presented. The dispersion curves derived using SAFE are denoted by ‘+’ and ‘o’ markers, which refer to symmetric and anti-symmetric Lamb waves, respectively. Note that the dispersion curves derived using Dispersion Calculator (which employs an analytical method to find the solution) are presented with dotted lines. Five elements with NURBS shape functions of degree p = 24 and C p 1 continuity along the cross section of the plate are considered to obtain the results presented below (which corresponds to the numerical results with the highest level of accuracy). It should be noted that a lower number of DOFs may be sufficient to derive the dispersion curves. The motivation behind selecting this excessive number of DOFs is twofold: (i) the error has already converged to the minimum possible error and (ii) an excellent agreement between the analytical and numerical dispersion curves even at high frequencies is obtained.
Before further discussing the different approaches considered here, it should be noted that, along with increasing the DOFs in modelling and consequently reducing the errors, two phenomena occur simultaneously which make the execution of a convergence study more complicated. As the accuracy of the model increases, along with decreasing the overall error in different frequencies, the existing error in estimating cut-off frequencies decreases as well. Consequently, except for first two modes which have zero cut-off frequencies, the other curves related to higher-order modes are shifted to the left in the diagram. (Remark: By increasing the number of DOFs, the structural behaviour becomes softer and thus, the cut-off frequencies are approximated from above.) Correspondingly, in order to compare the convergence and accuracy of different approaches, instead of comparing the results at a specified frequency, which is considered to be 5MHz for the first two modes, the error at the local maximum value of the third mode (‘A’ in Figure 5) is considered to assess the convergence behaviour. As clearly observed in Figure 5, the local extrema of the relative error curves related to the A 1 and S 1 modes, labeled by ‘A’ and ‘B’, move to the left and decrease as the number of elements increases for the case of FEM with second-order shape functions. The error for the local extremum at point ‘A’, which is located on the error curve of the anti-symmetric mode A 1 , decreases from 31.2 % to 2.48 % and then to 0.1 % by increasing the number of elements from one to three or seven elements, respectively. At the same time, the frequency at which the local extremum occurs moves from 6010 kHz to 4630 kHz and then to 4590 kHz. The same behavior is observed for the local extremum at point ‘B’, which is located on the error curve of the symmetric mode S 1 . By increasing the number of elements from 3 to seven, the error at extremum decreases from 14.2 % to 0.5 % , while the frequency reduces from 7570 kHz to 6630 kHz. As may be seen in Figure 5a, the local extremum ‘B’ cannot be seen in this figure, due to the fact that it occurs at frequencies higher than 10 MHz. This phenomenon is seen for all other modes and different numerical methods considered here. Therefore, considering the error at a specified frequency for higher modes can cause complications, when analysing the error changes by increasing the number of DOFs.
Since the other higher-order modes have approximately the same convergence behavior, and in order to avoid further repetitive figures, only the first three modes are considered in the comparison. Moreover, in Figure 6, the convergence rates and attainable accuracy (error level) of the various investigated approaches are illustrated.
Judging from the behaviour for the first three modes, i.e., A 0 , S 0 , and A 1 , the following conclusions may be drawn:
  • The highest convergence rate is achieved by the approach which takes advantage of NURBS shape functions of degree p with C p 1 continuity between elements and the k-refinement technique. Here, a total of five elements (non-zero knot spans) were considered, and, in each step, the degree of the shape functions and their continuity is increased.
  • Among the high-order shape functions with C 0 continuity between the elements, the hierarchical shape functions (p-FEM method; purple dash-dotted curve with ‘+’ markers) show the best performance. The attainable error level is slightly lower compared to the other methods. Here, no divergence is observed when the order of the shape functions is repeatedly increased. In the following order, the other C 0 continuous approaches can be ranked: NURBS shape functions (yellow dash-dotted curve with asterix markers), spectral shape functions (with GLL or GLC points; green dashed curve with ‘o’ markers or cyan dashed curve with triangle markers), and standard shape functions (dark blue dash-dotted curve with ‘x’ markers), respectively. As mentioned previously, it may also be observed that, as the shape function order increases in the standard shape functions with equidistant nodal distribution, the error increases, which is related to Runge’s phenomenon. Additionally, it should be noted that, in the thickness of the plate, a constant number of five elements were considered during the p-refinement approach.
  • Finally, the h-refinement technique shows the poorest performance among the other presented approaches. The methods which include h-refinement technique second order shape functions with C 0 continuity (FEM, blue dashed line with ‘o’ marker) and C 1 continuity (IGA, red dashed line with rectangular marker) were considered and the number of elements increased in each step. As it can be observed, the approach which takes advantage of C 1 continuity shows better performance than the standard FEM with C 0 continuity between the elements, which is the expected behavior.
It should be noted that the comparison of the computational time revealed that, although using different shape functions obviously has an effect on the efficiency of the approach, the current case studies are not suitable for an comprehensive assessment due to the small number of DOFs. For small scale investigations, only the number of DOFs or control points play a key role in the run time of the analysis. Thus, at this point, only the normalized simulation time is plotted against the number of DOFs, as depicted in Figure 7. As can be observed and was expected, the computational time increases significantly with the number of DOFs. Since this does not give us much additional insight, we also list the computational time when the relative error becomes less than 6 × 10 4 % for the different methods. For the specified criterion, the number of DOFs or control points and the corresponding computational time are summarized in Table 2.
As can be inferred from the data given in Table 2, the computational costs related to using the NURBS-enhanced method along with a k-refinement shows the best performance in comparison to other methods considered in this paper.
Note that these findings are in line with the results published in Ref. [8] for fully discretized models. The other important point which is worth mentioning is that the error at cut-off frequencies is the largest. It can be also noted that, as we move further away from the cut-off frequencies, the error decreases, but, again at higher frequencies, the error increases. The former point has been mentioned in previous studies [81] and using a higher number of DOFs was suggested in order to overcome this problem. Our study also reveals that, in case of using an insufficient number of DOFs, cut-off frequencies will also be estimated incorrectly. Consequently, the propagating waves will appear at higher frequencies than expected, which leads to higher errors in these frequency regimes. Overall, we have to stress that the obtained numerical results lead to the conclusion that the higher inter-element continuity is the decisive factor in the excellent performance of the NURBS-based technique.

4.2. Composite Case Study

The second case study is an eight layer composite plate with the ply lay-up sequence [ ± 45 , 0 , 90 ] s and 1 mm thickness. The material properties of a single layer in its principal direction are listed in Table 3.
As before, one-dimensional elements are considered. However, since shear and Lamb waves cannot be decoupled in this case, three DOFs at each node (FEM) or control point (IGA) at the cross-section of the plate are considered. In order to illustrate the accuracy of the different combinations of numerical schemes, the derived dispersion curves using the SAFE method and the Dispersion Calculator software are presented in Figure 8. The dispersion curves derived using the SAFE method are depicted with ‘+’ and ‘o’ markers, which refer to symmetric and anti-symmetric modes, while the dispersion curves derived using the analytical method are presented with dotted lines. For the simulation, five elements per layer with NURBS shape functions of degree p = 6 and C p 1 continuity along the cross section of the plate are considered. Note that the obtained numerical results are in excellent agreement with the reference solutions. As mentioned previously, this excessive number of DOFs is only considered to ensure convergence of the error to the minimum possible value, which has been found through the convergence study discussed in the remainder of this section. It should be noted that, depending on the mode shape of interest and the frequency, much fewer DOFs may suffice to derive the dispersion curves, especially for lower modes and frequencies.
In Figure 9, the convergence rates of the different approaches are illustrated. The analysis follows the same path as described in the previous subsection. Consequently, in order to compare the convergence and accuracy of the investigated approaches, the error at the local maximum value of an exemplary mode is considered.
From Figure 9, we observe that the calculated errors after convergence are approximately of order 10 4 , which is higher than the calculated error for the case of the isotropic plate. (Remark: For this case, it is conjectured that accuracy of the output obtained for the Dispersion Calculator software is limited, influencing the results of the convergence analysis.) As before, the highest rate of convergence is obtained by using a combination of NURBS shape functions of degree p with C p 1 continuity between elements and the k-refinement technique. However, the accuracy of different shape functions is the same for this case. Hierarchical shape functions (p-FEM method) (purple dash-dotted line with ‘+’ markers), and NURBS of order n with C 0 continuity (purple dash-dotted line with diamond marker), which take advantage of a p-refinement technique, show more robustness in comparison to the spectral shape functions (with GLL and GLC points). As expected, standard shape functions with an equidistant nodal distribution reveal some signs of instability as the order of the shape functions increases (see dark blue dash-dotted curve with ‘x’ markers). Finally, an h-refinement technique with C 0 continuity (FEM) shows the poorest performance among the presented approaches.

5. Conclusions

For the sake of improving the efficiency of the SAFE method, the effect of using different high-order shape functions is investigated. In order to compare the efficiency of different types of shape functions, the relative error of the SAFE method with respect to the analytical results is calculated. According to the obtained numerical results, it can be stated that, among the approaches considered here, NURBS based IGA in conjunction with the k-refinement approach is the most effective method. By utilizing the IGA approach, lower relative errors in comparison to other approaches using the same number of DOFs can be achieved. Due to the geometrical simplicity of the investigated model problems, the increase in accuracy is attributed to the higher inter-element continuity. This is the only property that sets the IGA-based approach apart from the other proposed numerical methods.
Furthermore, the simulation results reveal that using a lower number of points and consequently DOFs at the cross-section of the waveguide leads to inaccurate estimation of the cut-off frequencies. Due to this estimation error, the dispersion curves shift to the right on the frequency axis. Consequently, using an insufficient number of DOFs not only affects the accuracy of the dispersion curves at higher frequencies, but also causes higher errors in the dispersion curves of different modes around their cut-off frequencies.
Finally, it is worthwhile mentioning that although the SAFE method considerably decreases the computational requirements of wave propagation simulations in comparison to numerical methods that fully discretize the computational domain, such as the FEM and IGA, the improved performance will be even more pronounced when this method is used for material property characterization purposes, where a repeated number of simulations is required, as recently proposed by some researchers [82]. However, still the mentioned issues related to the approximations of complex waveguide geometries remain. To deal with these problems, the more versatile, but for the discussed problems closely related (basically identical) scaled boundary finite element method can be adopted [83].

Author Contributions

Conceptualization, S.E., L.P., E.M.K. and F.D.; methodology, E.M.K., S.E. and F.D.; software, E.M.K., S.E. and R.M.; validation, E.M.K.; resources, D.J., M.M. and J.R.; writing—original draft preparation, E.M.K.; writing—review and editing, E.M.K., J.R., F.D., L.P., R.M., D.J., M.M. and S.E.; supervision, J.R., D.J., M.M. and L.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
SAFESemi-Analytical Finite Element Method
DOFDegree of Freedom
FEMFinite Element Method
SEMSpectral Element Method
IGAIsogeometric Analysis
1D SAFEOne Dimensional Semi-Analytical Finite Element Method
2D SAFETwo Dimensional Semi-Analytical Finite Element Method
SHMStructural Health Monitoring
FEAFinite Element Analysis
GLLGauss–Lobatto–Legendre
GLCGauss–Lobatto–Chebyshev
p-FEMp-Version of Finite Element Method
NURBSNonuniform Rational B-Splines
CADComputer Aided Design
CAMComputer Aided Manufacturing

Appendix A. Basic methods and definitions

Appendix A.1. Eigenvalue Problem Solution and Analysis

Due to the existence of two unknown parameters, the eigenvalue problem may be solved using two approaches:
  • Assign real values to the circular frequency ω and solve either the second-order, cf. Equation (25), or the first-order eigenvalue problem, cf. Equation (26). The eigenvalues calculated in this case fall in one of the following three categories:
    (a)
    Real eigenvalues ( k = ± κ ) corresponding to propagating waves;
    (b)
    Complex eigenvalues ( k = ± κ real ± κ imag i ) corresponding to evanescent waves; (Remark: Evanescent waves are a type of waves whose energy is spatially concentrated in the vicinity of the source, i.e., after travelling a certain distance in the waveguide their amplitude becomes negligible [1].)
    (c)
    Imaginary eigenvalues ( k = ± κ i ) corresponding to standing waves. (Remark: Standing waves or non-oscillating evanescent waves are stationary in the structure and do not transport energy.)
  • Assign values to the wavenumber κ and solve the second-order eigenvalue problem given by Equation (25). It should be mentioned that, if real values are assigned to the wavenumber, only propagating modes are being considered, and, likewise, if we assign imaginary values to the wavenumber, only standing waves may be studied.

Appendix A.2. Phase Velocity and Group Velocity

The dynamic response of a waveguide is comprised of a superposition of many harmonics which may be propagating waves, evanescent waves, and standing waves. All these are combined to the observed wave packet in motion. The propagating velocity of each harmonic or each wave mode is called phase velocity and is defined as follows:
c p = ω κ real ,
while the group velocity is considered as the velocity of the observed packet of all traveling waves and may be presented as follows [39,84]:
c g = d ω d κ real = Q L T ( K ^ 2 + 2 κ K 3 ) Q R 2 ω Q L T M Q R ,
where Q R and Q L refer to the right and left eigenvector of the eigenvalue problem, respectively, which are derived using Equation (26).

Appendix A.3. Cut-Off Frequency

Cut-off frequencies may be obtained by inserting κ = 0 into Equation (25) and solving the following eigenvalue problem:
( K 1 ω 2 M ) Q = 0 .
The cut-off frequency is the frequency at which standing waves become propagating waves.

Appendix B. Mode Identification Approaches

The mode identification process is one of the most challenging tasks when dealing with eigenvalue problems. If only the propagating waves (real-valued wavenumbers) should be considered, the orthogonality property of modes may be exploited to fulfill this task [85]:
{ Q ^ m T B m Q ^ m 0 Q ^ n T B m Q ^ m = 0 ,
where n and m refer to the n- th and m- th eigenvectors derived using Equation (26) at the same frequency. If we assume that for two consecutive steps the orthogonality relation still holds, the orthogonality relation for the following step may be expressed as:
{ Q ^ m T ( ω r ) B m ( ω r ) Q ^ m + 1 ( ω r + 1 ) 0 Q ^ n T ( ω r ) B m ( ω r ) Q ^ m + 1 ( ω r + 1 ) = 0 ,
where Q ^ m and Q ^ n refer to the eigenvectors of the m-th and n-th mode, derived using Equation (26) with the frequency ω r , while Q ^ m + 1 and Q ^ n + 1 refer to the eigenvectors of the m-th and n-th mode, derived using Equation (26) with the frequency ω r + 1 . Equation (A5) simply relies on the assumption that the orthogonality relation between the modes for two consecutive frequency steps, namely ω r and ω r + 1 , still remains valid.
The second method for distinguishing different dispersion curves which belong to the same mode is mode tracking. In this method, by using the points existing on a curve, the next point that is expected to be on the same curve is extrapolated. By comparing the guessed point with the wavenumbers derived at the desired step and calculating the existing error, it will be decided which one of these wavenumbers belong to the same mode. In order to execute this approach, Lagrange polynomials [9], Taylor series expansions, and Padé expansions [80] can be used. Mode tracking is an essential part of the calculation of dispersion diagrams by which curves belonging to different modes can be distinguished. The fundamental idea of this scheme is to extrapolate points on the curve based on already known ones. The extrapolated value is then compared to the numerically determined ones and thus the results can be assigned to the correct modes. In the following, three different mode tracking methods are briefly discussed:
1
Quadratic polynomial [9]:
Considering the fact that in general one may fit any p + 1 points by a polynomial of p-th degree, Lagrange polynomials constitute a viable choice [63]. The fitted polynomial takes the form
P p ( x ) = i = 0 p f ( x i ) L i ( x ) ,
where the cardinal functions L i ( x ) are expressed as
L i ( x ) = j = 0 , j i p x x j x i x j .
According to this definition, if we consider a second-order Lagrange polynomial for the wavenumber k at the desired frequency ω n , it is extrapolated by using three points as follows:
k ( ω n ) = k ( ω n 3 ) L 1 ( ω n ) + k ( ω n 2 ) L 2 ( ω n ) + k ( ω n 1 ) L 3 ( ω n ) .
By considering a constant frequency step, Equation (A8) may be further simplified:
k ( ω n ) = k ( ω n 3 ) 3 k ( ω n 2 ) + 3 k ( ω n 1 ) .
2
Taylor series expansion:
Besides Lagrange polynomials, we can also make use of a Taylor series expansion to derive the extrapolating polynomial. If we consider z = x x 0 , the p-th order Taylor expansion of f is defined as [86]
f ( x ) = i = 0 p c i z i + R p ( x ) ,
where R p ( x ) is the remainder, and the factor c i is defined as
c i = f ( i ) ( x 0 ) i ! .
If we consider the terms up to the third-order, the wavenumber may be extrapolated using the following expression:
k ( ω n ) = k ( ω n 1 ) + k ( ω n 1 ) Δ ω + k ( Δ ω ) 2 + k ( ω n 1 ) ( Δ ω ) 3 ,
where Δ ω = ω n ω n 1 . By using the backward difference method, the derivatives of the function k ( ω n 1 ) are derived as follows:
k ( ω n 1 ) = k ( ω n 1 ) k ( ω n 2 ) Δ ω ,
k ( ω n 1 ) = k ( ω n 1 ) 2 k ( ω n 2 ) + k ( ω n 3 ) ( Δ ω ) 2 ,
k ( ω n 1 ) = k ( ω n 1 ) 3 k ( ω n 2 ) + 3 k ( ω n 3 ) k ( ω n 4 ) ( Δ ω ) 3 .
3
Padé series expansion:
Another extrapolation method, which has been proposed by Gravenkamp et al. [80] for tracing dispersion curves, uses Padé expansions. A Padé approximation of f is a rational function with a numerator of degree L and denominator of degree M, determined such that its power expansion matches up to including the power L + M . Considering z = x x 0 , the Padé expansion of f may be presented as follows:
f ( z ) = P L ( z ) Q L ( z ) + R ( z ) = [ L / M ] f ( z ) + R ( z ) ,
where R ( z ) is the remainder of order L + M + 1 , and [ L / M ] f ( z ) is called the Padé approximation of order [ L , M ] of function f. The coefficients of the Padé expansion may be derived by replacing the left term in Equation (A16) by a Taylor series using Equation (A10) [86]:
f ( z ) = c 0 + c 1 z + c 2 z + = a 0 + a 1 z + + a L z L b 0 + b 1 z + + b M z M ,
where b 0 = 1 . By multiplying both sides by the denominator of the term on the right-hand-side of the equation and comparing the coefficients of the terms with the same power, a set of linear algebraic equations is derived, which may be solved in order to calculate the Padé expansion coefficients. Using this method, the following set of equations may be derived for a Padé expansion of order [1/2] using Equation (A17):
a 0 = c 0 ,
a 1 = c 1 + c 0 b 1 ,
a 2 = c 2 + c 1 b 1 + c 0 b 2 ,
a 3 = c 3 + c 2 b 1 + c 1 b 2 + c 0 b 3 .
For a Padé expansion of order [1/2], we obtain a 2 = a 3 = b 3 = 0 . By solving the set of algebraic equations, the coefficients of a Padé expansion of order [1/2] are
a 0 = c 0 ,
a 1 = c 1 + c 0 b 1 ,
b 1 = c 2 + c 0 b 2 c 1 ,
b 2 = c 2 2 c 3 c 1 c 1 2 c 2 c 0 ,
where c i ( i = 0 , 1 , 2 , 3 ) are given by Equations (A11) and (A13)–(A15). Consequently, by replacing these coefficients in Equation (A16), the Padé expansion of order [1/2] of the wavenumber may be expressed as:
k ( ω n ) = a 0 + a 1 Δ ω 1 + b 1 Δ ω + b 2 Δ ω 2 .

References

  1. Bartoli, I.; Marzani, A.; di Scalea, F.L.; Viola, E. Modeling wave propagation in damped waveguides of arbitrary cross-section. J. Sound Vib. 2006, 295, 685–707. [Google Scholar] [CrossRef]
  2. Marzani, A.; Viola, E.; Bartoli, I.; di Scalea, F.L.; Rizzo, P. A semi-analytical finite element formulation for modeling stress wave propagation in axisymmetric damped waveguides. J. Sound Vib. 2008, 318, 488–505. [Google Scholar] [CrossRef]
  3. Marzani, A. Time–transient response for ultrasonic guided waves propagating in damped cylinders. Int. J. Solids Struct. 2008, 45, 6347–6368. [Google Scholar] [CrossRef]
  4. Schmiechen, P. Travelling Wave Speed Coincidence. Ph.D. Thesis, University of London, London, UK, 1997. [Google Scholar]
  5. Marburg, S.; Nolte, B. Computational Acoustics of Noise Propagation in Fluids—Finite and Boundary Element Methods; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar] [CrossRef]
  6. Bathe, K.J. Finite Element Procedures; Prentice Hall: Hoboken, NJ, USA, 2006. [Google Scholar]
  7. Bartoli, I.; Lanza di Scalea, F.; Fateh, M.; Viola, E. Modeling guided wave propagation with application to the long-range defect detection in railroad tracks. Ndt E Int. 2005, 38, 325–334. [Google Scholar] [CrossRef]
  8. Willberg, C.; Duczek, S.; Perez, J.V.; Schmicker, D.; Gabbert, U. Comparison of different higher order finite element schemes for the simulation of Lamb waves. Comput. Methods Appl. Mech. Eng. 2012, 241–244, 246–261. [Google Scholar] [CrossRef]
  9. Lowe, M. Matrix techniques for modeling ultrasonic waves in multilayered media. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 1995, 42, 525–542. [Google Scholar] [CrossRef]
  10. Rokhlin, S.I.; Wang, L. Stable recursive algorithm for elastic wave propagation in layered anisotropic media: Stiffness matrix method. J. Acoust. Soc. Am. 2002, 112, 822–834. [Google Scholar] [CrossRef]
  11. Astaneh, A.V.; Guddati, M.N. Dispersion analysis of composite acousto-elastic waveguides. Compos. Part Eng. 2017, 130, 200–216. [Google Scholar] [CrossRef]
  12. Joseph, R.; Li, L.; Haider, M.F.; Giurgiutiu, V. Hybrid SAFE-GMM approach for predictive modeling of guided wave propagation in layered media. Eng. Struct. 2019, 193, 194–206. [Google Scholar] [CrossRef]
  13. Su, Z.; Ye, L.; Lu, Y. Guided Lamb waves for identification of damage in composite structures: A review. J. Sound Vib. 2006, 295, 753–780. [Google Scholar] [CrossRef]
  14. Willberg, C.; Duczek, S.; Vivar-Perez, J.M.; Ahmad, Z.A.B. Simulation Methods for Guided Wave-Based Structural Health Monitoring: A Review. Appl. Mech. Rev. 2015, 67. [Google Scholar] [CrossRef]
  15. Aalami, B. Waves in Prismatic Guides of Arbitrary Cross Section. J. Appl. Mech. 1973, 40, 1067–1072. [Google Scholar] [CrossRef]
  16. Lagasse, P.E. Higher-order finite-element analysis of topographic guides supporting elastic surface waves. J. Acoust. Soc. Am. 1973, 53, 1116–1122. [Google Scholar] [CrossRef]
  17. Zienkiewicz, O.C. The Finite Element Method in Engineering Science; McGraw-Hill: London, UK; New York, NY, USA, 1971. [Google Scholar]
  18. Mazzotti, M.; Bartoli, I.; Miniaci, M.; Marzani, A. Wave dispersion in thin-walled orthotropic waveguides using the first order shear deformation theory. Thin-Walled Struct. 2016, 103, 128–140. [Google Scholar] [CrossRef]
  19. Hladky-Hennion, A.C. Finite element analysis of the propagation of acoustic waves in waveguides. J. Sound Vib. 1996, 194, 119–136. [Google Scholar] [CrossRef]
  20. Damljanović, V.; Weaver, R.L. Forced response of a cylindrical waveguide with simulation of the wavenumber extraction problem. J. Acoust. Soc. Am. 2004, 115, 1582–1591. [Google Scholar] [CrossRef]
  21. Loveday, P.W. Simulation of piezoelectric excitation of guided waves using waveguide finite elements. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2008, 55, 2038–2045. [Google Scholar] [CrossRef]
  22. Coccia, S.; Bartoli, I.; Marzani, A.; di Scalea, F.L.; Salamone, S.; Fateh, M. Numerical and experimental study of guided waves for detection of defects in the rail head. Ndt Int. 2011, 44, 93–100. [Google Scholar] [CrossRef]
  23. Treyssède, F. Elastic waves in helical waveguides. Wave Motion 2008, 45, 457–470. [Google Scholar] [CrossRef]
  24. Cong, M.; Wu, X.; Liu, R. Dispersion analysis of guided waves in the finned tube using the semi-analytical finite element method. J. Sound Vib. 2017, 401, 114–126. [Google Scholar] [CrossRef]
  25. Mazzotti, M.; Bartoli, I.; Marzani, A.; Viola, E. A coupled SAFE-2.5D BEM approach for the dispersion analysis of damped leaky guided waves in embedded waveguides of arbitrary cross-section. Ultrasonics 2013, 53, 1227–1241. [Google Scholar] [CrossRef] [PubMed]
  26. Mazzotti, M.; Marzani, A.; Bartoli, I. Dispersion analysis of leaky guided waves in fluid-loaded waveguides of generic shape. Ultrasonics 2014, 54, 408–418. [Google Scholar] [CrossRef] [PubMed]
  27. Treyssède, F.; Nguyen, K.; Bonnet-BenDhia, A.S.; Hazard, C. Finite element computation of trapped and leaky elastic waves in open stratified waveguides. Wave Motion 2014, 51, 1093–1107. [Google Scholar] [CrossRef]
  28. Nguyen, K.; Treyssède, F.; Hazard, C. Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods. J. Sound Vib. 2015, 344, 158–178. [Google Scholar] [CrossRef]
  29. Zuo, P.; Yu, X.; Fan, Z. Numerical modeling of embedded solid waveguides using SAFE-PML approach using a commercially available finite element package. Ndt Int. 2017, 90, 11–23. [Google Scholar] [CrossRef]
  30. Duan, W.; Kirby, R.; Mudge, P.; Gan, T.H. A one-dimensional, numerical approach for computing the eigenmodes of elastic waves in buried pipelines. J. Sound Vib. 2016, 384, 177–193. [Google Scholar] [CrossRef]
  31. Matuszyk, P.J. Modeling of guided circumferential SH and Lamb-type waves in open waveguides with semi-analytical finite element and Perfectly Matched Layer method. J. Sound Vib. 2017, 386, 295–310. [Google Scholar] [CrossRef]
  32. Zuo, P.; Fan, Z. SAFE-PML approach for modal study of waveguides with arbitrary cross sections immersed in inviscid fluid. J. Sound Vib. 2017, 406, 181–196. [Google Scholar] [CrossRef]
  33. tian Deng, Q.; chun Yang, Z. Propagation of guided waves in bonded composite structures with tapered adhesive layer. Appl. Math. Model. 2011, 35, 5369–5381. [Google Scholar] [CrossRef]
  34. Mei, H.; Giurgiutiu, V. Predictive 1D and 2D guided-wave propagation in composite plates using the SAFE approach. In Proceedings of the Health Monitoring of Structural and Biological Systems XII, Denver, CO, USA, 5–8 March 2018; Kundu, T., Ed.; International Society for Optics and Photonics: Bellingham, WA, USA, 2018; Volume 10600, pp. 215–225. [Google Scholar] [CrossRef]
  35. Duan, W.; Gan, T.H. Investigation of guided wave properties of anisotropic composite laminates using a semi-analytical finite element method. Compos. Part Eng. 2019, 173, 106898. [Google Scholar] [CrossRef]
  36. Cui, R.; di Scalea, F.L. On the identification of the elastic properties of composites by ultrasonic guided waves and optimization algorithm. Compos. Struct. 2019, 223, 110969. [Google Scholar] [CrossRef]
  37. Inoue, D.; Hayashi, T. Transient analysis of leaky Lamb waves with a semi-analytical finite element method. Ultrasonics 2015, 62, 80–88. [Google Scholar] [CrossRef] [PubMed]
  38. Duan, W.; Kirby, R. Guided wave propagation in buried and immersed fluid-filled pipes: Application of the semi analytic finite element method. Comput. Struct. 2019, 212, 236–247. [Google Scholar] [CrossRef]
  39. Han, X.; Liu, G.R.; Xi, Z.C.; Lam, K.Y. Characteristics of waves in a functionally graded cylinder. Int. J. Numer. Methods Eng. 2002, 53, 653–676. [Google Scholar] [CrossRef]
  40. Hayashi, T.; Inoue, D. Calculation of leaky Lamb waves with a semi-analytical finite element method. Ultrasonics 2014, 54, 1460–1469. [Google Scholar] [CrossRef]
  41. Gavrić, L. Computation of propagative waves in free rail using a finite element technique. J. Sound Vib. 1995, 185, 531–543. [Google Scholar] [CrossRef]
  42. Hayashi, T.; Song, W.J.; Rose, J.L. Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example. Ultrasonics 2003, 41, 175–183. [Google Scholar] [CrossRef]
  43. Damljanović, V.; Weaver, R.L. Propagating and evanescent elastic waves in cylindrical waveguides of arbitrary cross section. J. Acoust. Soc. Am. 2004, 115, 1572–1581. [Google Scholar] [CrossRef]
  44. Hayashi, T.; Tamayama, C.; Murase, M. Wave structure analysis of guided waves in a bar with an arbitrary cross-section. Ultrasonics 2006, 44, 17–24. [Google Scholar] [CrossRef]
  45. Predoi, M.V.; Castaings, M.; Hosten, B.; Bacon, C. Wave propagation along transversely periodic structures. J. Acoust. Soc. Am. 2007, 121, 1935–1944. [Google Scholar] [CrossRef]
  46. Loveday, P.W. Analysis of Piezoelectric Ultrasonic Transducers Attached to Waveguides Using Waveguide Finite Elements. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2007, 54, 2045–2051. [Google Scholar] [CrossRef] [PubMed]
  47. Loveday, P.W. Semi-analytical finite element analysis of elastic waveguides subjected to axial loads. Ultrasonics 2009, 49, 298–300. [Google Scholar] [CrossRef] [PubMed]
  48. Mazzotti, M.; Marzani, A.; Bartoli, I.; Viola, E. Guided waves dispersion analysis for prestressed viscoelastic waveguides by means of the SAFE method. Int. J. Solids Struct. 2012, 49, 2359–2372. [Google Scholar] [CrossRef]
  49. Setshedi, I.I.; Loveday, P.W.; Long, C.S.; Wilke, D.N. Estimation of rail properties using semi-analytical finite element models and guided wave ultrasound measurements. Ultrasonics 2019, 96, 240–252. [Google Scholar] [CrossRef] [PubMed]
  50. Loveday, P.W.; Long, C.S.; Ramatlo, D.A. Mode repulsion of ultrasonic guided waves in rails. Ultrasonics 2018, 84, 341–349. [Google Scholar] [CrossRef]
  51. Onipede, O.; Dong, S.; Kosmatka, J. Natural vibrations and waves in pretwisted rods. Compos. Eng. 1994, 4, 487–502. [Google Scholar] [CrossRef]
  52. Mazúch, T. Wave dispersion modelling in anisotropic shells and rods by the finite element method. J. Sound Vib. 1996, 198, 429–438. [Google Scholar] [CrossRef]
  53. Volovoi, V.; Hodges, D.; Berdichevsky, V.; Sutyrin, V. Dynamic dispersion curves for non-homogeneous, anisotropic beams with cross-sections of arbitrary geometry. J. Sound Vib. 1998, 215, 1101–1120. [Google Scholar] [CrossRef]
  54. Treyssède, F.; Laguerre, L. Investigation of elastic modes propagating in multi-wire helical waveguides. J. Sound Vib. 2010, 329, 1702–1716. [Google Scholar] [CrossRef]
  55. Düster, A. High Order Finite Elements for Three-Dimensional, Thin-Walled Nonlinear Continua. Ph.D. Thesis, Technical University Munich, Munich, Germany, 2002. [Google Scholar]
  56. Rezaiee-Pajand, M.; Gharaei-Moghaddam, N.; Ramezani, M. Higher-order assumed strain plane element immune to mesh distortion. Eng. Comput. 2020. ahead-of-print. [Google Scholar] [CrossRef]
  57. Kalkowski, M.K.; Muggleton, J.M.; Rustighi, E. Axisymmetric semi-analytical finite elements for modelling waves in buried/submerged fluid-filled waveguides. Comput. Struct. 2018, 196, 327–340. [Google Scholar] [CrossRef]
  58. Xiao, D.; Han, Q.; Liu, Y.; Li, C. Guided wave propagation in an infinite functionally graded magneto-electro-elastic plate by the Chebyshev spectral element method. Compos. Struct. 2016, 153, 704–711. [Google Scholar] [CrossRef]
  59. Treyssède, F. Spectral element computation of high-frequency leaky modes in three-dimensional solid waveguides. J. Comput. Phys. 2016, 314, 341–354. [Google Scholar] [CrossRef]
  60. Seyfaddini, F.; Nguyen, V.H. NURBS-enriched semi-analytical finite element method (SAFE) for calculation of wave dispersion in heterogeneous waveguides. In Proceedings of the 24e Congrès Français de Mécanique, Brest, France, 26–30 August 2019. [Google Scholar]
  61. Seyfaddini, F.; Nguyen-Xuan, H.; Nguyen, V.H. A semi-analytical isogeometric analysis for wave dispersion in functionally graded plates immersed in fluids. Acta Mech. 2021, 232, 15–32. [Google Scholar] [CrossRef]
  62. Seyfaddini, F.; Nguyen-Xuan, H.; Nguyen, V.H. Semi-analytical IGA-based computation of wave dispersion in fluid-coupled anisotropic poroelastic plates. Int. J. Mech. Sci. 2021, 212, 106830. [Google Scholar] [CrossRef]
  63. Boyd, J.P. Chebyshev and Fourier Spectral Methods, 2nd ed.; Dover: Mineola, NY, USA, 2001. [Google Scholar]
  64. Atkinson, K.E. An Introduction to Numerical Analysis, 2nd ed.; John Wiley & Sons: Hoboken, NJ, USA, 1989; p. xvi+693. [Google Scholar]
  65. Rao, S. The Finite Element Method in Engineering, 5th ed.; Elsevier Inc.: Amsterdam, The Netherlands, 2010. [Google Scholar] [CrossRef]
  66. Zienkiewicz, O.; Taylor, R. Finite Element Method: Volume 1—The Basis, 5th ed.; Butterworth-Heinemann: Oxford, UK, 2000. [Google Scholar]
  67. Seriani, G.; Oliveira, S. Numerical Modeling of mechanical wave propagation. Riv. Del Nuovo C. 2020, 43, 459–514. [Google Scholar] [CrossRef]
  68. Trefethen, L.; Weideman, J. Two results on polynomial interpolation in equally spaced points. J. Approx. Theory 1991, 65, 247–260. [Google Scholar] [CrossRef]
  69. Zrahia, U.; Bar-Yoseph, P. Space-time spectral element method for solution of second-order hyperbolic equations. Comput. Methods Appl. Mech. Eng. 1994, 116, 135–146. [Google Scholar] [CrossRef]
  70. Seriani, G.; Priolo, E. Spectral element method for acoustic wave simulation in heterogeneous media. Finite Elem. Anal. Des. 1994, 16, 337–348. [Google Scholar] [CrossRef]
  71. Patera, A.T. A spectral element method for fluid dynamics: Laminar flow in a channel expansion. J. Comput. Phys. 1984, 54, 468–488. [Google Scholar] [CrossRef]
  72. Duczek, S.; Gravenkamp, H. Mass lumping techniques in the spectral element method: On the equivalence of the row-sum, nodal quadrature, and diagonal scaling methods. Comput. Methods Appl. Mech. Eng. 2019, 353, 516–569. [Google Scholar] [CrossRef]
  73. Vu, T.H.; Deeks, A.J. Use of higher-order shape functions in the scaled boundary finite element method. Int. J. Numer. Methods Eng. 2006, 65, 1714–1733. [Google Scholar] [CrossRef]
  74. Zienkiewicz, O.; Taylor, R.; Zhu, J. The Finite Element Method Set: Its Basis and Fundamentals, 6th ed.; Butterworth-Heinemann: Oxford, UK, 2005; pp. 103–137. [Google Scholar] [CrossRef]
  75. Hughes, T.; Cottrell, J.; Bazilevs, Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput. Methods Appl. Mech. Eng. 2005, 194, 4135–4195. [Google Scholar] [CrossRef]
  76. Agrawal, V.; Gautam, S.S. IGA: A Simplified Introduction and Implementation Details for Finite Element Users. J. Inst. Eng. Ser. 2019, 100, 561–585. [Google Scholar] [CrossRef]
  77. Cottrell, J.A.; Hughes, T.J.R.; Bazilevs, Y. NURBS as a Pre-Analysis Tool: Geometric Design and Mesh Generation. In Isogeometric Analysis; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2009; pp. 19–68. [Google Scholar] [CrossRef]
  78. Piegl, L.; Tiller, W. The NURBS Book, 2nd ed.; Springer: New York, NY, USA, 1997. [Google Scholar]
  79. German Aerospace Center (DLR) Institute of Structures and Design, Center for Lightweight Production Technology Augsburg. The Dispersion Calculator: A Free Software for Calculating Dispersion Curves of Guided Waves in Multilayered Composites 2018. Available online: https://www.dlr.de/zlp/en/desktopdefault.aspx/tabid-14332/24874_read-61142/ (accessed on 1 June 2022).
  80. Gravenkamp, H.; Birk, C.; Song, C. The computation of dispersion relations for axisymmetric waveguides using the Scaled Boundary Finite Element Method. Ultrasonics 2014, 54, 1373–1385. [Google Scholar] [CrossRef]
  81. Mei, H.; Giurgiutiu, V. Guided wave excitation and propagation in damped composite plates. Struct. Health Monit. 2019, 18, 690–714. [Google Scholar] [CrossRef]
  82. Kudela, P.; Radzienski, M.; Fiborek, P.; Wandowski, T. Elastic constants identification of fibre-reinforced composites by using guided wave dispersion curves and genetic algorithm for improved simulations. Compos. Struct. 2021, 272, 114178. [Google Scholar] [CrossRef]
  83. Gravenkamp, H.; Saputra, A.A.; Duczek, S. High-order shape functions in the scaled boundary finite element method revisited. Arch. Comput. Methods Eng. 2019, 28, 473–494. [Google Scholar] [CrossRef]
  84. Finnveden, S. Evaluation of modal density and group velocity by a finite element method. J. Sound Vib. 2004, 273, 51–75. [Google Scholar] [CrossRef]
  85. Loveday, P.W.; Long, C.S. Time domain simulation of piezoelectric excitation of guided waves in rails using waveguide finite elements. Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems 2007; Tomizuka, M., Yun, C.B., Giurgiutiu, V., Eds.; International Society for Optics and Photonics: Bellingham, WA, USA, 2007; Volume 6529, pp. 316–325. [Google Scholar] [CrossRef]
  86. Simon Širca, M.H. Computational Methods in Physics, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 2018. [Google Scholar] [CrossRef]
Figure 1. (a) 1D, (b) 2D, and (c) 3D elements with the corresponding displacement DOFs at each node (FEM).
Figure 1. (a) 1D, (b) 2D, and (c) 3D elements with the corresponding displacement DOFs at each node (FEM).
Mca 27 00063 g001
Figure 2. Shape functions, p = 5 , (a) Lagrange polynomials based on an equidistant nodal distribution, (b) Lagrange polynomials based on a GLL nodal distribution, (c) Lagrange polynomials based on a GLC nodal distribution, (d) p-FEM, and (e) IGA—the graphical representation of the basis functions related to the knot vector in Equation (38).
Figure 2. Shape functions, p = 5 , (a) Lagrange polynomials based on an equidistant nodal distribution, (b) Lagrange polynomials based on a GLL nodal distribution, (c) Lagrange polynomials based on a GLC nodal distribution, (d) p-FEM, and (e) IGA—the graphical representation of the basis functions related to the knot vector in Equation (38).
Mca 27 00063 g002aMca 27 00063 g002b
Figure 3. 1D and 2D element with three displacement DOFs at each node in the SAFE method.
Figure 3. 1D and 2D element with three displacement DOFs at each node in the SAFE method.
Mca 27 00063 g003
Figure 4. Dispersion curves of the first six modes of propagating Lamb waves in a 1 mm thick aluminum plate; ⋯ Dispersion Calculator, + SAFE method (symmetric modes), o SAFE method (anti-symmetric modes); (a) wavenumber, (b) phase velocity, (c) group velocity vs. frequency.
Figure 4. Dispersion curves of the first six modes of propagating Lamb waves in a 1 mm thick aluminum plate; ⋯ Dispersion Calculator, + SAFE method (symmetric modes), o SAFE method (anti-symmetric modes); (a) wavenumber, (b) phase velocity, (c) group velocity vs. frequency.
Mca 27 00063 g004
Figure 5. Calculated percentage error for the aluminum plate using FEM, (a) one element, (b) three elements, and (c) seven elements with 2nd order shape functions.
Figure 5. Calculated percentage error for the aluminum plate using FEM, (a) one element, (b) three elements, and (c) seven elements with 2nd order shape functions.
Mca 27 00063 g005
Figure 6. Convergence curves derived using the shape functions along with the specified refinement technique in aluminum plate problem, (a) 1st mode, (b) 2nd mode, (c) 3rd mode, (d) legend.
Figure 6. Convergence curves derived using the shape functions along with the specified refinement technique in aluminum plate problem, (a) 1st mode, (b) 2nd mode, (c) 3rd mode, (d) legend.
Mca 27 00063 g006
Figure 7. Normalized time vs. number of nodes or control points.
Figure 7. Normalized time vs. number of nodes or control points.
Mca 27 00063 g007
Figure 8. Dispersion curves of the first six modes of propagating waves in a 1 mm thick composite plate; …. Dispersion Calculator, + SAFE method (symmetric modes), o SAFE method (anti-symmetric modes); (a) wavenumber, (b) phase velocity, (c) group velocity vs. frequency.
Figure 8. Dispersion curves of the first six modes of propagating waves in a 1 mm thick composite plate; …. Dispersion Calculator, + SAFE method (symmetric modes), o SAFE method (anti-symmetric modes); (a) wavenumber, (b) phase velocity, (c) group velocity vs. frequency.
Mca 27 00063 g008
Figure 9. Convergence curves derived using the mentioned shape functions along with the specified refinement technique in the composite plate problem, (a) 1st mode, (b) 2nd mode, (c) 3rd mode, (d) legend.
Figure 9. Convergence curves derived using the mentioned shape functions along with the specified refinement technique in the composite plate problem, (a) 1st mode, (b) 2nd mode, (c) 3rd mode, (d) legend.
Mca 27 00063 g009
Table 1. Material properties of aluminum.
Table 1. Material properties of aluminum.
Young’s Modulus ( GPa )Poisson’s RatioDensity ( kg / m 3 )
E ν ρ
73.10.332780
Table 2. Different approaches and their corresponding normalized time and required DOFs when the relative error is less than 6 × 10 4 %.
Table 2. Different approaches and their corresponding normalized time and required DOFs when the relative error is less than 6 × 10 4 %.
MethodDOFsNormalized Time
FEM (h-refinement, C 0 , p = 2)611.0000
NURBS (h-refinement, C 1 , p = 2)330.1962
NURBS (p-refinement, C 0 )210.0609
NURBS (k-refinement, C p 1 )100.0120
SEM (p-refinement, GLL)210.0609
SEM (p-refinement, CGL)210.0609
p-FEM (p-refinement)210.0609
FEM (p-refinement, equidistant)210.0609
Table 3. Material properties of T800/924 laminate (elastic constants in GPa) [1].
Table 3. Material properties of T800/924 laminate (elastic constants in GPa) [1].
C11C12C13C22C23C33C44C55C66
164.7085.4535.45311.3004.73911.3003.2806.0006.000
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Mirzaee Kakhki, E.; Rezaeepazhand, J.; Duvigneau, F.; Pahlavan, L.; Makvandi, R.; Juhre, D.; Moavenian, M.; Eisenträger, S. On the Use of High-Order Shape Functions in the SAFE Method and Their Performance in Wave Propagation Problems. Math. Comput. Appl. 2022, 27, 63. https://doi.org/10.3390/mca27040063

AMA Style

Mirzaee Kakhki E, Rezaeepazhand J, Duvigneau F, Pahlavan L, Makvandi R, Juhre D, Moavenian M, Eisenträger S. On the Use of High-Order Shape Functions in the SAFE Method and Their Performance in Wave Propagation Problems. Mathematical and Computational Applications. 2022; 27(4):63. https://doi.org/10.3390/mca27040063

Chicago/Turabian Style

Mirzaee Kakhki, Elyas, Jalil Rezaeepazhand, Fabian Duvigneau, Lotfollah Pahlavan, Resam Makvandi, Daniel Juhre, Majid Moavenian, and Sascha Eisenträger. 2022. "On the Use of High-Order Shape Functions in the SAFE Method and Their Performance in Wave Propagation Problems" Mathematical and Computational Applications 27, no. 4: 63. https://doi.org/10.3390/mca27040063

APA Style

Mirzaee Kakhki, E., Rezaeepazhand, J., Duvigneau, F., Pahlavan, L., Makvandi, R., Juhre, D., Moavenian, M., & Eisenträger, S. (2022). On the Use of High-Order Shape Functions in the SAFE Method and Their Performance in Wave Propagation Problems. Mathematical and Computational Applications, 27(4), 63. https://doi.org/10.3390/mca27040063

Article Metrics

Back to TopTop