Abstract
We compare contemporary practices of global approximation using cubic B-splines in conjunction with double multiplicity of inner knots (-continuous) with older ideas of utilizing local Hermite interpolation of third degree. The study is conducted within the context of the Galerkin-Ritz formulation, which forms the background of the finite element structural analysis. Numerical results, concerning static and eigenvalue analysis of rectangular elastic structures in plane stress conditions, show that both interpolations lead to identical results, a finding that supports the view that they are mathematically equivalent.
1. Introduction
Structural analysis is usually performed using commercial codes that include finite elements of low (usually first or second) degree, where the accuracy of the calculations increases by mesh refinement (-version). Alternatively, keeping the number and the position of the nodal points unaltered, the numerical solution improves using polynomials of higher degree (-version) [1].
As an extension of the above -version “philosophy,” tensor-product Lagrange polynomials as well as CAD-based (Gordon-Coons) macroelements—based on several interpolations—have been used [2–4]. The aforementioned macroelements integrate the solid modelling (CAD: computer-aided-design) with the analysis (CAE: computer-aided-engineering). In more detail, these macroelements use the same global approximation for both the geometry and the displacement vector. In order to avoid the undesired numerical oscillations caused by Lagrange polynomials of high degree, the next generation of CAD-based macroelements replaced them with tensor-product B-splines [5]. Since 2005, the nonuniform-rational-B-splines (NURBS) interpolation has started to prevail [6].
A careful study of literature reveals that most of recent papers referring to the so-called isogeometric analysis (IGA) start with some essentials on the definition of B-splines and relevant recursive formulas due to de Boor [7]. It should be recalled that NURBS is an extension of B-splines (nonuniform and rational) modified on the basis of weighting coefficients, thus producing fully controlled sculptured surfaces [8, 9]. In a B-splines expansion, the multiplicity of the inner knots plays a significant role in the continuity of the variables. In general, the multiplicity of inner knots per breakpoint in combination with a piecewise polynomial of degree ensures -continuity of the variable (here: displacement components) [7, 9]. Thus considering cubic B-splines () in conjunction with double inner knots (), -continuity is ensured. Höllig [5, page 93] has solved plane stress problems using B-splines of degree , 3, 4, and 5, but his study is not a complete investigation on the influence of the multiplicity and corresponding continuity of variables involved.
On the other point of view, tensor-product Hermite elements of third degree have been proposed for the solution of fourth-order problems, such as plate-bending problems, using Galerkin-Ritz formulation. The need for smoother () global basis functions is also encountered in second-order problems when collocation finite element methods are utilized [10, page 66].
With these situations in mind, we next examine the relationship between particular Hermite elements of third degree and cubic B-splines elements with multiplicity . Given that both a global interpolation (B-splines) and a local one (Hermite elements) ensure -continuity, the question arises whether and how these are mutually equivalent.
Numerical examples concerning rectangular structures in plane stress conditions reveal that both approximations are equivalent to one another.
2. The General Elastodynamic Problem
2.1. Governing Equations
In the case of an elastic, isotropic, and homogeneous 2D structure, the governing equation in the domain is given by where is the stress tensor, is the body force, is the stress operator, is the mass density, and is the acceleration vector.
Considering the Hookean elasticity matrix in 2D problems (plane stress or plane strain), the relationship between the stress, , and the strain is In addition, the relationship between the strain and the displacement vector, , is Therefore, the final governing equation (Naviér-Kirchhoff) is written as follows: where is a quadratic operator.
2.2. Boundary Conditions
Let us consider the field of displacements on the domain , which is required to fulfil the governing equation (5). The boundary conditions that correspond to this problem are of two types as follows:(a)essential conditions, such as on (Dirichlet type),(b)natural conditions of the type on (Neumann type), with denoting the traction, the components of which are related to the stress tensor as follows: where and are the components of the outward normal vector on the boundary . The total boundary is .
3. Global and Local Interpolation
Below we present the two interpolations that will be compared to each other.
3.1. The B-Splines as a Global Functional Set
In the evolution of time, B-splines have appeared in two different forms.
3.1.1. Older Definitions
The concept of B-splines was published in 1946 and some years later by Schoenberg [11] and Schoenberg and Whitney [12]. It refers to the points with , which we wish to interpolate through a multiply-defined function . The points are called “breakpoints.” For a cubic polynomial , the desired properties are as follows.(i)In each interval , with , function is a cubic polynomial.(ii)Function and the first and second derivatives are continuous at the above points.Introducing the truncated power as which is -continuous, the original expression consists of a power series in the form [11, 12]: It is apparent that (9) includes terms and ensures -continuity, because for simplicity we considered .
Alternatively, (9) can be modified so as to include additional truncated polynomials of second degree; that is, Obviously, (10) includes terms and ensures -continuity.
3.1.2. Contemporary Procedures
The breakthrough was made in 1972, independently by de Boor [7] and Cox [13], who both achieved the B-spline function and its derivatives to be rapidly computed through recursive formulas. The latter procedures exist today in MATLAB (spline toolbox) under the name “spcol.”
In brief, for a nondecreasing set of breakpoints and for a certain polynomial degree “,” we define the knot vector “”: which highly depends on the chosen multiplicity of internal knots (usually single or double), as follows.(i)Multiplicity : (ii)Multiplicity :
Therefore, (12) and (13) lead to the unified relationship:
Based on the abovementioned computed knot vector , the vector of control points is denoted by where the number of control points () is related to the number of elements in the knot vector as follows: The th B-splines basis functions of degree (order ), denoted by , are defined as
Then, for every position , with normalized coordinate , we can determine the values of basis functions, or , , as follows.(i)The Cartesian coordinate is approximated in terms of control points as (ii)The variable is approximated as
It is worth mentioning that the coefficients in (19) are generally different than the nodal values associated with the breakpoints, except for the ends where and .
3.2. Piecewise Hermite Approximation
We refer again to the sequence , which was mentioned in Section 3.1. Let be an arbitrary element in a finite element partition of the interval . A -cubic element is obtained by implementing the function and derivatives at the ends and of each element.
The actual element is transformed linearly to the master element by the map , where . On , interpolating and at the end nodes , the cubic Hermite expansion has the form: where the Hermite basis functions satisfy the interpolation properties at end nodes , : for local nodal indices and . Using properties (21) we may construct the Hermite cubics directly as
3.3. Contemporary versus Older Definitions
For the particular case of (cubic splines) under this study, it is trivial to prove that(i)when the multiplicity equals one (), the basis functions involved in (19) are equivalent (not identical) to those basis functions in (9);(ii)when the multiplicity equals two (), the basis functions involved in (19) are equivalent (not identical) to those basis functions in (10);(iii)therefore, a lot of research conducted in 1960s using the older framework ((9) and (10)) may be repeated on the new unified framework of (19).In this paper, we are concerned with , which corresponds only to (10).
In order to get a better idea, the one-dimensional unit domain was divided into four subdivisions (, ), and then (i) local shape functions using Hermite polynomials of third degree and (ii) global basis functions using cubic B-splines (in conjunction with ) are compared in Figure 1. Although each set of four local Hermite shape functions was separately plotted (in the intervals , , , and ), Figure 1(a) shows that all of them (as a whole) present also a global character. Although both of the aforementioned sets have this global character, a first look at Figure 1 does not reveal any apparent relationship between the Hermite polynomials (shown in Figure 1(a)) and the cubic B-splines (shown in Figure 1(b)). In more detail, concerning the Hermite polynomials, those two shape functions that correspond to flexural DOFs vary between 0 and 1, whereas those two that correspond to slopes obtain both positive and negative values. In contrast, all B-splines basis functions are nonnegative and generally less than unity.
(a)
(b)
4. Tensor-Product Shape Functions
Let us consider a rectangular domain in . The axis origin is chosen at corner A, whereas the Cartesian axes and lie on the sides AB and AD, respectively. Without loss of generality, the sides (AB, CD) and (BC, DA) are uniformly divided into and segments, respectively, thus leading to breakpoints along AB or CD, as well as breakpoints along BC or DA. In this paper, we use a unique polynomial degree (cubic: ) for both - and -directions. As a result, the univariate function along the side AB can be interpolated through a piecewise B-splines polynomial of third degree in , whereas the function along the side DA can be interpolated through a piecewise B-splines polynomial of third degree in .
4.1. Tensor-Product B-Splines
Given the uniform subdivisions and of the intervals and , respectively, the breakpoints along each of the four sides (AB, BC, CD, and DA) are determined. Moreover, given the multiplicity of internal knots, , as well as the polynomial degrees , the control points in the - and -directions are also determined. If the patch is curvilinear, then - and -coordinates have to be replaced by the - and -normalized coordinates, respectively.
While in older B-splines formulation [11, 12] the degrees of freedom are associated with the nodal points lying at the intersections of th and th lines perpendicular to the axes and passing through the breakpoints , in this—let us say—“modern” formulation we have to deal only with the tensor-product of control points. Since the multiplicity of internal knots is , the tensor-product consists of coefficients , for each displacement component.
Therefore, the two-dimensional global shape functions are given by (the double subscript is here to emphasize the two directions)
Since the multiplicity is equal to two, then the univariate approximation is -continuous (whereas for 2-D ); is the standard reference square.
In general, the control points are divided into two categories, that is, in the interior of the domain and (generally) near the boundary (). In more detail, if a side of the quadrilateral ABCD (e.g., AB) is straight (as in this study) the corresponding control points lie on this side (AB). In contrast, if the side is curved, then the extreme control points ( and ) will belong to the boundary and even they coincide with the corners (e.g., A and B), whereas the rest will be either inside or outside the domain in accordance with the curvature of the curve AB.
4.2. Tensor-Product Hermite
Consider the Hermite cubics for on defined in (22) and introduce , in the orthogonal coordinate . The tensor-product on consists of the 16 functions: on with , and , for nodes and . The interpolant is where , at , and so on.
Using the maps and , we obtain the corresponding tensor-product shape functions on the rectangle . Considering the chain rule, the derivatives in the and directions in (25) transform to partial derivatives in and . There are now four element shape functions associated with each node to interpolate the corner values of , , , and . For each variable , the element has 16 degrees of freedom. The form of the element interpolant (25) transforms to where and we have suppressed the superscript on the element shape functions for notational simplicity. Details are given in the Appendix.
This tensor-product basis ensures -continuity, as is evident from the form of the global basis functions at an interface between two elements in the discretization. Equivalently, differentiating (26) with respect to and evaluating on side yields a cubic function on this side that is uniquely determined by the function value and its derivative at the endpoints of that side.
5. Galerkin-Ritz Procedure
5.1. General
For the given partial differential equation (5), we seek an approximate solution to (5) which is a linear combination of the bivariate global basis functions , : In this paper, the candidate functions are cubic B-splines ((23) in conjunction with generalized coefficients ) or Hermite polynomials ((25) in conjunction with kinematic degrees of freedom ). It should be clarified that the variable in (27) stands for either the horizontal or the vertical displacement components at any point of the elastic structure.
Based on the global shape functions involved in (27), we can apply the well-known Galerkin-Ritz method.
Without loss of generality, the boundary consists of breakpoints (which correspond to control points) under Dirichlet and ones (which correspond to control points) under Neumann boundary conditions. For the sake of brevity, below we limit the discussion in the two typical cases of boundary conditions, that is, Dirichlet-type and Neumann-type boundary conditions.
In the general elastodynamic problem the equations of motion are where is the mass matrix, is the damping matrix, and is the stiffness matrix, with After imposing the boundary conditions, the dimensions of each matrix in (28) become , where is the number of equations (i.e., the number of unrestricted degrees of freedom: coefficients or kinematic DOFs). For a certain choice of subdivisions, has the same value either the B-splines or the Hermite formulation is implemented.
5.2. Numerical Implementation
The numerical procedure is performed as follows. The domain is uniformly divided into a certain number of subdivisions (cells), thus leading to nodes.
Concerning the determination of the proper Gaussian quadrature, we start from the observation that the elements of the mass matrix are products of two basis functions, each of piecewise th (i.e., third) degree. In the particular case of a rectangular domain which is the topic of this paper, the integrant becomes of piecewise sixth degree and thus it requires four-point Gauss quadrature per direction, that is, sixteen Gauss points per integration cell.
6. Numerical Examples
Two MATLAB codes of very similar architecture were developed on a standard PC Pentium IV as follows:(i)a code based on rectangular 4-node (32-DOF) Hermite elements;(ii)a code based on contemporary cubic B-splines in conjunction with double inner knots; the basis functions and their derivatives were created using the “spcol” function, which exists in the Spline Toolkit.
The eigenvalues were calculated using the standard “eig” function.
The theory is now elucidated by three examples: one example deals with static analysis while the remaining two examples deal with the eigenvalue analysis of rectangular structures in plane stress conditions. With respect to the exact solution, , the quality of the numerical solution is evaluated in terms of the relative error, which was calculated as follows:
Example 1 (cantilever beam). Consider a cantilever beam subject to end load as shown in Figure 2. The following parameters were used: kPa; , m, m, kN, and plane stress conditions.
The solution is given by Timoshenko and Goodier [14] as
where is the moment of inertia and for a beam with rectangular cross-section and unit thickness it is given by
The rectangular domain is uniformly discretized using of subdivisions along the - and -directions, respectively, as shown in Figure 2 (for ). For these models, Figure 3 shows that the flexural displacement at the middle of the free side BC (see Figure 2) converges to the exact solution even for a very small number of subdivisions.
(a)
(b)
(c)
(d)
(e)
Example 2 (square plate in plane stress). Consider a square plate of uniform thickness under plane stress conditions, which is fixed along its entire boundary. The following parameters were used: and .
In the absence of an analytical solution, a parametric analysis using ANSYS (PLANE42 elements) shows that for a uniform mesh of 90 × 90 elements the solution has practically converged. The relative errors (cf. (30)) for the first six eigenvalues are shown in Figure 4. The horizontal axis corresponds to the previously mentioned number of equations (unrestricted DOF), .
Example 3 (clamped plate in plane stress). Consider a rectangular plate of dimensions m of uniform thickness under plane stress conditions, which is clamped along one of its smallest sides. The following parameters were used: and . Keeping a constant ratio of subdivisions () the plate is discretized into elements of square form. As an exact solution we have considered the numerical solution using 600 × 150 PLANE42 ANSYS elements. The results, shown in Figure 5, reveal monotonic convergence.
7. Discussion
Bicubic B-spline interpolation has been previously used for the finite element analysis of plates [5, 15] based on Galerkin-Ritz formulation, as well as in conjunction with collocation methods for potential problems [7, 10].
In the advent of the isogeometric analysis [6], the older piecewise interpolations have been substituted by knots and control points, but in some cases it may be a recast of older formulations. Within this framework, it was shown that B-splines interpolation in conjunction with double inner knots ensures -continuity. In other words, although it seems that cubic B-splines interpolation is a global approximation of both geometry and variable within the entire rectangular domain, it is equivalent of splitting the domain into subdivisions along the - and -directions, respectively, and considering local approximation within each of the aforementioned Hermite elements. A similar coincidence has been previously noticed in problems of one dimension, however, in conjunction with the collocation method [16].
According to Carey and Oden [10, page 68], unfortunately, there is one serious shortcoming regarding the abovementioned tensor-product Hermite elements—to retain -continuity, the discretization must be restricted to consist of rectangles in two dimensions and rectangular prisms in three dimensions. So far, two alternatives have been proposed: (i) the use of nonconforming and (ii) simplex elements. Alternatively, it is anticipated that the tensor-product B-splines formulation (in conjunction with double inner knots) will be applicable to problems with curved boundary, at the cost of determining the Jacobian and its inverse.
Remaining within the context of rectangular elements, higher order macroelements can be used. One possibility is to use similar global functions with those previously used for plate bending [17]. In addition, modified Hermite polynomials, a family of -rectangular elements of increasing degree, have been proposed (see, e.g., [10, page 67]). Based on the equivalency revealed in this paper, it is anticipated that similar results will be obtained when substituting the two aforementioned alternatives with proper breakpoints and corresponding multiplicities. For example, at inner points only the displacement components can be considered, whereas at extreme points their derivatives can be considered as well. While explicit expressions have been previously derived, equivalent B-splines with proper select of multiplicities are anticipated to be applicable.
8. Conclusions
The Galerkin-Ritz method using contemporary tensor-product cubic B-splines with double inner knots (global interpolation) was found to be numerically equivalent to the older -continuous Hermite elements of third degree (local interpolation) when applied to rectangular domains. This finding is valid for both static and eigenvalue problems. The superiority of the contemporary tensor-product cubic B-splines is probably that, in principle, they can work for nonrectangular structures as well. Moreover, ongoing research reveals that global B-splines interpolation can readily implement several modified Hermite alternatives such as higher degree or different number of function values and derivatives at internal and external points.
Appendix
Interpolation Using Cubic Hermite Polynomials
The cubic Hermite polynomials are where represents either - or -normalized coordinates.
Based on these polynomials we can construct their tensor products to derive the global functions. For a rectangular element ABCD having four degrees per node () in the -direction and another four DOFs () in the -direction, the tensor-product basis on master element consists of the 16 functions grouped in sets of 4 per node: Renaming the corners A, B, C, and D with the local numbers 1, 2, 3, and 4, respectively, and considering the abovementioned four degrees of freedom per node for the -displacement: and another four degrees of freedom per node for the -displacement: the vector of degrees of freedom is of the dimensions 32 × 1 as follows: with and so on.
Finally, for each element the interpolation of the components, and , will be
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.