# Finite Element Analysis of Structures UsingC1-Continuous Cubic B-Splines or Equivalent Hermite Elements

Authors: Christopher Provatidis
Journal: Journal of Structures
 (2014)

---

## Abstract

We compare contemporary practices of global approximation using cubic B-splines in conjunction with double multiplicity of inner knots (C1-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.


---

## Body

## 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 (h-version). Alternatively, keeping the number and the position of the nodal points unaltered, the numerical solution improves using polynomials of higher degree (p-version) [1].As an extension of the abovep-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-calledisogeometric 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 p ensures Cp-λ-continuity of the variable (here: displacement components) [7, 9]. Thus considering cubic B-splines (p=3) in conjunction with double inner knots (λ=2), C1-continuity is ensured. Höllig [5, page 93] has solved plane stress problems using B-splines of degree n=2, 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 (C1) 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λ=2. Given that both aglobal interpolation (B-splines) and alocal one (Hermite elements) ensure C1-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
(1)LTσ+b=ρu¨inΩ,
where σ=σxσyτxyT is the stress tensor, b=bxbyT is the body force, LT=∂/∂x0∂/∂y0∂/∂y∂/∂x is the stress operator, ρ is the mass density, and u¨=u¨v¨T is the acceleration vector.Considering the Hookean elasticity matrix in 2D problems (plane stress or plane strain),(2)E=d11d120d21d22000d33,
the relationship between the stress, σ, and the strain ε=εxεyγxyT is
(3)σ=Eε.
In addition, the relationship between the strain and the displacement vector, u=uvT, is
(4)ε=Lu.
Therefore, the final governing equation (Naviér-Kirchhoff) is written as follows:
(5)Du+b=ρu¨,
where
(6)D=LTEL
is a quadratic operator.


### 2.2. Boundary Conditions

Let us consider the field of displacementsux,y=ux,yvx,yT 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 asu=u-t on Γ1 (Dirichlet type),(b)
natural conditions of the typet=t-t on Γ2 (Neumann type), with t denoting the traction, the components of which are related to the stress tensor as follows:
(7)tx=σxn-x+τxyn-y,ty=τyxn-x+σyn-y,where n-x and n-y are the components of the outward normal vector on the boundary Γ2. The total boundary is Γ=Γ1+Γ2.


## 2.1. Governing Equations

In the case of an elastic, isotropic, and homogeneous 2D structure, the governing equation in the domainΩ is given by
(1)LTσ+b=ρu¨inΩ,
where σ=σxσyτxyT is the stress tensor, b=bxbyT is the body force, LT=∂/∂x0∂/∂y0∂/∂y∂/∂x is the stress operator, ρ is the mass density, and u¨=u¨v¨T is the acceleration vector.Considering the Hookean elasticity matrix in 2D problems (plane stress or plane strain),(2)E=d11d120d21d22000d33,
the relationship between the stress, σ, and the strain ε=εxεyγxyT is
(3)σ=Eε.
In addition, the relationship between the strain and the displacement vector, u=uvT, is
(4)ε=Lu.
Therefore, the final governing equation (Naviér-Kirchhoff) is written as follows:
(5)Du+b=ρu¨,
where
(6)D=LTEL
is a quadratic operator.


## 2.2. Boundary Conditions

Let us consider the field of displacementsux,y=ux,yvx,yT 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 asu=u-t on Γ1 (Dirichlet type),(b)
natural conditions of the typet=t-t on Γ2 (Neumann type), with t denoting the traction, the components of which are related to the stress tensor as follows:
(7)tx=σxn-x+τxyn-y,ty=τyxn-x+σyn-y,where n-x and n-y are the components of the outward normal vector on the boundary Γ2. The total boundary is Γ=Γ1+Γ2.


## 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 (x0,u0),(x1,u1),…,(xn,un) with x0<x1<⋯<xn, which we wish to interpolate through a multiply-defined function ux. The points x0,x1,…,xn are called “breakpoints.” For a cubic polynomial (p=3), the desired properties are as follows.(i)
In each intervalxi-1≤x≤xi, with i=1,2,…,n, function ux is a cubic polynomial.(ii)
Functionux and the first and second derivatives are continuous at the above points.Introducing the truncated power as
(8)x-xi+m=0,x≤xix-xim,x>xi,
which is Cm-1-continuous, the original expression consists of a power series in the form [11, 12]:
(9)ux=a0+a1x+a2x2+a3x3+∑j=1n-1bjx-xj+3.
It is apparent that (9) includes (n+3) terms and ensures C2-continuity, because for simplicity we considered m=3.Alternatively, (9) can be modified so as to include additional truncated polynomials of second degree; that is,
(10)ux=a0+a1x+a2x2+a3x3+∑j=1n-1bjx-xj+2+∑j=1n-1cjx-xj+3.
Obviously, (10) includes 2(n+1) terms and ensures C1-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(n+1)breakpoints(x0,x1,…,xn-1,xn) and for a certain polynomial degree “p,” we define theknot vector “V”:
(11)V=v0,…,vm,
which highly depends on the chosenmultiplicityλ of internal knots (usually single or double), as follows.(i)
Multiplicityλ=1:
(12)Vλ=1=x0,…,x0︸p+1,x1,x2,…,xn-1,xn,…,xn︸p+1.(ii)
Multiplicityλ=2:
(13)Vλ=2=x0,…,x0︸p+1,x1,x1︸2,x2,x2︸2,…,xn-1,xn-1︸2,xn,…,xn︸p+1.Therefore, (12) and (13) lead to the unified relationship:
(14)m=2p+1+λn-1-1,λ=1,2.Based on the abovementioned computed knot vectorV, the vector of control points is denoted by
(15)P=P0,…,Pnc,
where the number of control points (nc+1) is related to the number of elements in the knot vector (m+1) as follows:
(16)m=nc+p+1.
The ith B-splines basis functions of degree p (order p+1), denoted by Ni,pu, are defined as
(17)Ni,0u=1ifui≤u<ui+10otherwise,Ni,pu=u-uiui+p-uiNi,p-1u+ui+p+1-uui+p+1-ui+1Ni+1,p-1u.Then, for every positionx∈0,L, with normalized coordinate ξ=x/L∈0,1, we can determine the values of nc+1 basis functions, Ni,px or Ni,pξ,  i=0,…,nc, as follows.(i)
The Cartesian coordinate is approximated in terms of control points as(18)xξ=∑i=0nNiξ·xi.(ii)
The variable is approximated as(19)uξ=∑i=0nNiξ·ai.It is worth mentioning that the coefficientsai in (19) are generally different than the nodal values ui associated with the breakpoints, except for the ends where a0≡u0 and an≡un.


### 3.2. Piecewise Hermite Approximation

We refer again to the sequencex0<x1<⋯<xn, which was mentioned in Section 3.1. Let Ωe=xi-1,xi be an arbitrary element in a finite element partition of the interval a,b≡x0,xn. A C1-cubic element is obtained by implementing the function and derivatives at the ends xi-1 and xi of each element.The actual elementΩe is transformed linearly to the master element Ω=0,1 by the map ξ=1+r/2, where r=1+2x-xi-1+xi/xi-xi-1/2. On Ω, interpolating u and u′≡du/dξ at the end nodes j=1,2, the cubic Hermite expansion has the form:
(20)uξ=∑j=12ujEj0ξ+∑j=12uj′Ej1ξ,
where the Hermite basis functions Ej0,Ej1 satisfy the interpolation properties at end nodes ξ1=0,  ξ2=1:
(21)Ej0ξk=δjk,Ej1ξk=0,dEj0dξξk=0,dEj1dξξk=δjk,
for local nodal indices j=1,2 and k=1,2. Using properties (21) we may construct the Hermite cubics directly as
(22)E00ξ=1-3ξ2+2ξ3,E10ξ=3ξ2-2ξ3,E01ξ=ξ-2ξ2+ξ3,E11ξ=-ξ2+ξ3.


### 3.3. Contemporary versus Older Definitions

For the particular case ofp=3 (cubic splines) under this study, it is trivial to prove that(i)
when the multiplicity equalsone (λ=1), the (n+3) basis functions involved in (19) areequivalent (not identical) to those basis functions in (9);(ii)
when the multiplicity equalstwo (λ=2), the 2(n+1) basis functions involved in (19) areequivalent (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 λ=2, which corresponds only to (10).In order to get a better idea, the one-dimensional unit domain was divided into four subdivisions (Ω=0,1, nx=4), and then (i) local shape functions using Hermite polynomials of third degree and (ii) global basis functions using cubic B-splines (in conjunction with λ=2) are compared in Figure 1. Although each set of four local Hermite shape functions was separately plotted (in the intervals [0,1/4], [1/4,2/4], [2/4,3/4], and [3/4,1]), 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) Local and (b) global shape functions.
(a)
(b)


## 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 (x0,u0),(x1,u1),…,(xn,un) with x0<x1<⋯<xn, which we wish to interpolate through a multiply-defined function ux. The points x0,x1,…,xn are called “breakpoints.” For a cubic polynomial (p=3), the desired properties are as follows.(i)
In each intervalxi-1≤x≤xi, with i=1,2,…,n, function ux is a cubic polynomial.(ii)
Functionux and the first and second derivatives are continuous at the above points.Introducing the truncated power as
(8)x-xi+m=0,x≤xix-xim,x>xi,
which is Cm-1-continuous, the original expression consists of a power series in the form [11, 12]:
(9)ux=a0+a1x+a2x2+a3x3+∑j=1n-1bjx-xj+3.
It is apparent that (9) includes (n+3) terms and ensures C2-continuity, because for simplicity we considered m=3.Alternatively, (9) can be modified so as to include additional truncated polynomials of second degree; that is,
(10)ux=a0+a1x+a2x2+a3x3+∑j=1n-1bjx-xj+2+∑j=1n-1cjx-xj+3.
Obviously, (10) includes 2(n+1) terms and ensures C1-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(n+1)breakpoints(x0,x1,…,xn-1,xn) and for a certain polynomial degree “p,” we define theknot vector “V”:
(11)V=v0,…,vm,
which highly depends on the chosenmultiplicityλ of internal knots (usually single or double), as follows.(i)
Multiplicityλ=1:
(12)Vλ=1=x0,…,x0︸p+1,x1,x2,…,xn-1,xn,…,xn︸p+1.(ii)
Multiplicityλ=2:
(13)Vλ=2=x0,…,x0︸p+1,x1,x1︸2,x2,x2︸2,…,xn-1,xn-1︸2,xn,…,xn︸p+1.Therefore, (12) and (13) lead to the unified relationship:
(14)m=2p+1+λn-1-1,λ=1,2.Based on the abovementioned computed knot vectorV, the vector of control points is denoted by
(15)P=P0,…,Pnc,
where the number of control points (nc+1) is related to the number of elements in the knot vector (m+1) as follows:
(16)m=nc+p+1.
The ith B-splines basis functions of degree p (order p+1), denoted by Ni,pu, are defined as
(17)Ni,0u=1ifui≤u<ui+10otherwise,Ni,pu=u-uiui+p-uiNi,p-1u+ui+p+1-uui+p+1-ui+1Ni+1,p-1u.Then, for every positionx∈0,L, with normalized coordinate ξ=x/L∈0,1, we can determine the values of nc+1 basis functions, Ni,px or Ni,pξ,  i=0,…,nc, as follows.(i)
The Cartesian coordinate is approximated in terms of control points as(18)xξ=∑i=0nNiξ·xi.(ii)
The variable is approximated as(19)uξ=∑i=0nNiξ·ai.It is worth mentioning that the coefficientsai in (19) are generally different than the nodal values ui associated with the breakpoints, except for the ends where a0≡u0 and an≡un.


## 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 (x0,u0),(x1,u1),…,(xn,un) with x0<x1<⋯<xn, which we wish to interpolate through a multiply-defined function ux. The points x0,x1,…,xn are called “breakpoints.” For a cubic polynomial (p=3), the desired properties are as follows.(i)
In each intervalxi-1≤x≤xi, with i=1,2,…,n, function ux is a cubic polynomial.(ii)
Functionux and the first and second derivatives are continuous at the above points.Introducing the truncated power as
(8)x-xi+m=0,x≤xix-xim,x>xi,
which is Cm-1-continuous, the original expression consists of a power series in the form [11, 12]:
(9)ux=a0+a1x+a2x2+a3x3+∑j=1n-1bjx-xj+3.
It is apparent that (9) includes (n+3) terms and ensures C2-continuity, because for simplicity we considered m=3.Alternatively, (9) can be modified so as to include additional truncated polynomials of second degree; that is,
(10)ux=a0+a1x+a2x2+a3x3+∑j=1n-1bjx-xj+2+∑j=1n-1cjx-xj+3.
Obviously, (10) includes 2(n+1) terms and ensures C1-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(n+1)breakpoints(x0,x1,…,xn-1,xn) and for a certain polynomial degree “p,” we define theknot vector “V”:
(11)V=v0,…,vm,
which highly depends on the chosenmultiplicityλ of internal knots (usually single or double), as follows.(i)
Multiplicityλ=1:
(12)Vλ=1=x0,…,x0︸p+1,x1,x2,…,xn-1,xn,…,xn︸p+1.(ii)
Multiplicityλ=2:
(13)Vλ=2=x0,…,x0︸p+1,x1,x1︸2,x2,x2︸2,…,xn-1,xn-1︸2,xn,…,xn︸p+1.Therefore, (12) and (13) lead to the unified relationship:
(14)m=2p+1+λn-1-1,λ=1,2.Based on the abovementioned computed knot vectorV, the vector of control points is denoted by
(15)P=P0,…,Pnc,
where the number of control points (nc+1) is related to the number of elements in the knot vector (m+1) as follows:
(16)m=nc+p+1.
The ith B-splines basis functions of degree p (order p+1), denoted by Ni,pu, are defined as
(17)Ni,0u=1ifui≤u<ui+10otherwise,Ni,pu=u-uiui+p-uiNi,p-1u+ui+p+1-uui+p+1-ui+1Ni+1,p-1u.Then, for every positionx∈0,L, with normalized coordinate ξ=x/L∈0,1, we can determine the values of nc+1 basis functions, Ni,px or Ni,pξ,  i=0,…,nc, as follows.(i)
The Cartesian coordinate is approximated in terms of control points as(18)xξ=∑i=0nNiξ·xi.(ii)
The variable is approximated as(19)uξ=∑i=0nNiξ·ai.It is worth mentioning that the coefficientsai in (19) are generally different than the nodal values ui associated with the breakpoints, except for the ends where a0≡u0 and an≡un.


## 3.2. Piecewise Hermite Approximation

We refer again to the sequencex0<x1<⋯<xn, which was mentioned in Section 3.1. Let Ωe=xi-1,xi be an arbitrary element in a finite element partition of the interval a,b≡x0,xn. A C1-cubic element is obtained by implementing the function and derivatives at the ends xi-1 and xi of each element.The actual elementΩe is transformed linearly to the master element Ω=0,1 by the map ξ=1+r/2, where r=1+2x-xi-1+xi/xi-xi-1/2. On Ω, interpolating u and u′≡du/dξ at the end nodes j=1,2, the cubic Hermite expansion has the form:
(20)uξ=∑j=12ujEj0ξ+∑j=12uj′Ej1ξ,
where the Hermite basis functions Ej0,Ej1 satisfy the interpolation properties at end nodes ξ1=0,  ξ2=1:
(21)Ej0ξk=δjk,Ej1ξk=0,dEj0dξξk=0,dEj1dξξk=δjk,
for local nodal indices j=1,2 and k=1,2. Using properties (21) we may construct the Hermite cubics directly as
(22)E00ξ=1-3ξ2+2ξ3,E10ξ=3ξ2-2ξ3,E01ξ=ξ-2ξ2+ξ3,E11ξ=-ξ2+ξ3.


## 3.3. Contemporary versus Older Definitions

For the particular case ofp=3 (cubic splines) under this study, it is trivial to prove that(i)
when the multiplicity equalsone (λ=1), the (n+3) basis functions involved in (19) areequivalent (not identical) to those basis functions in (9);(ii)
when the multiplicity equalstwo (λ=2), the 2(n+1) basis functions involved in (19) areequivalent (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 λ=2, which corresponds only to (10).In order to get a better idea, the one-dimensional unit domain was divided into four subdivisions (Ω=0,1, nx=4), and then (i) local shape functions using Hermite polynomials of third degree and (ii) global basis functions using cubic B-splines (in conjunction with λ=2) are compared in Figure 1. Although each set of four local Hermite shape functions was separately plotted (in the intervals [0,1/4], [1/4,2/4], [2/4,3/4], and [3/4,1]), 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) Local and (b) global shape functions.
(a)
(b)


## 4. Tensor-Product Shape Functions

Let us consider a rectangular domainΩ=(ABCD)=0,a×0,b in R2. The axis origin is chosen at corner A, whereas the Cartesian axes x and y lie on the sides AB and AD, respectively. Without loss of generality, the sides (AB, CD) and (BC, DA) are uniformly divided into nx and ny segments, respectively, thus leading to (nx+1) breakpoints along AB or CD, as well as (ny+1) breakpoints along BC or DA. In this paper, we use a unique polynomial degree (cubic: p=3) for both x- and y-directions. As a result, the univariate function u(x,0) along the side AB can be interpolated through a piecewise B-splines polynomial of third degree in x, whereas the function u(0,y) along the side DA can be interpolated through a piecewise B-splines polynomial of third degree in y.


### 4.1. Tensor-Product B-Splines

Given the uniform subdivisionsnx and ny of the intervals [0,a] and [0,b], respectively, the breakpoints along each of the four sides (AB, BC, CD, and DA) are determined. Moreover, given the multiplicity of internal knots, λ=2, as well as the polynomial degrees p=3, the control points in the x- and y-directions are also determined. If the patch is curvilinear, then x- and y-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 nx+1×ny+1 nodal points xij lying at the intersections of ith and jth lines perpendicular to the axes and passing through thebreakpointsxi,yj, in this—let us say—“modern” formulation we have to deal only with the tensor-product of qcontrol points. Since the multiplicity of internal knots is λ=2, the tensor-product consists of q=4nx+1ny+1 coefficients aij, for each displacement component.Therefore, the two-dimensional global shape functions are given by (the double subscript is here to emphasize the two directions)(23)λ=2:ϕijx,y=Nix·Njy,i=0,…,2nx+1∧j=0,…,2ny+1.Since the multiplicity is equal to two, then the univariate approximation isC1-continuous (whereas for 2-D u∈C1,1Ωst);  Ωst=[0,1]×[0,1] is thestandard reference square.In general, theq control points are divided into two categories, that is, nc,in in the interior of the domain Ω and nc,b (generally)near the boundary (q=nc,in+nc,b). 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 (P0 and Pn) 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 cubicsEj0ξ,Ej1ξ for j=1,2 on ξ∈0,1 defined in (22) and introduce Ej0η, Ej1η in the orthogonal coordinate η. The tensor-product on Ω consists of the 16 functions:
(24)Eklmnξ,η=Ekmξ×Elnη
on Ωst=0,1×0,1 with m=0,1, and n=0,1, for nodes k=1,2 and l=1,2. The interpolant is
(25)Uξ,η=∑k,l=12‍∑m,n=01‍uklmnEkmξElnη,
where ukl00=uξk,ηl, ukl10=∂u/∂ξ at ξk,ηl, and so on.Using the mapsξ→x and η→y, we obtain the corresponding tensor-product shape functions on the rectangle Ωe. Considering the chain rule, the derivatives in the ξ and η directions in (25) transform to partial derivatives in x and y. There are now four element shape functions associated with each node to interpolate the corner values of u, ux≡∂u/∂x, uy≡∂u/∂y, and uxy≡∂2u/∂x∂y. For each variable u, the element has 16 degrees of freedom. The form of the element interpolant (25) transforms to
(26)uhex,y=∑k,l=12‍∑m,n=01‍uklmnEkmξElnη,
where uklmn={u(xk,yl),∂u(xk,yl)/∂x,∂u(xk,yl)/∂y,∂2u(xk,yl)/∂x∂y} and we have suppressed the superscript e on the element shape functions for notational simplicity. Details are given in the Appendix.This tensor-product basis ensuresC1-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 x and evaluating uxx,y on side x=xk yields a cubic function uxy on this side that is uniquely determined by the function value ux and its derivative uxy at the endpoints of that side.


## 4.1. Tensor-Product B-Splines

Given the uniform subdivisionsnx and ny of the intervals [0,a] and [0,b], respectively, the breakpoints along each of the four sides (AB, BC, CD, and DA) are determined. Moreover, given the multiplicity of internal knots, λ=2, as well as the polynomial degrees p=3, the control points in the x- and y-directions are also determined. If the patch is curvilinear, then x- and y-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 nx+1×ny+1 nodal points xij lying at the intersections of ith and jth lines perpendicular to the axes and passing through thebreakpointsxi,yj, in this—let us say—“modern” formulation we have to deal only with the tensor-product of qcontrol points. Since the multiplicity of internal knots is λ=2, the tensor-product consists of q=4nx+1ny+1 coefficients aij, for each displacement component.Therefore, the two-dimensional global shape functions are given by (the double subscript is here to emphasize the two directions)(23)λ=2:ϕijx,y=Nix·Njy,i=0,…,2nx+1∧j=0,…,2ny+1.Since the multiplicity is equal to two, then the univariate approximation isC1-continuous (whereas for 2-D u∈C1,1Ωst);  Ωst=[0,1]×[0,1] is thestandard reference square.In general, theq control points are divided into two categories, that is, nc,in in the interior of the domain Ω and nc,b (generally)near the boundary (q=nc,in+nc,b). 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 (P0 and Pn) 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 cubicsEj0ξ,Ej1ξ for j=1,2 on ξ∈0,1 defined in (22) and introduce Ej0η, Ej1η in the orthogonal coordinate η. The tensor-product on Ω consists of the 16 functions:
(24)Eklmnξ,η=Ekmξ×Elnη
on Ωst=0,1×0,1 with m=0,1, and n=0,1, for nodes k=1,2 and l=1,2. The interpolant is
(25)Uξ,η=∑k,l=12‍∑m,n=01‍uklmnEkmξElnη,
where ukl00=uξk,ηl, ukl10=∂u/∂ξ at ξk,ηl, and so on.Using the mapsξ→x and η→y, we obtain the corresponding tensor-product shape functions on the rectangle Ωe. Considering the chain rule, the derivatives in the ξ and η directions in (25) transform to partial derivatives in x and y. There are now four element shape functions associated with each node to interpolate the corner values of u, ux≡∂u/∂x, uy≡∂u/∂y, and uxy≡∂2u/∂x∂y. For each variable u, the element has 16 degrees of freedom. The form of the element interpolant (25) transforms to
(26)uhex,y=∑k,l=12‍∑m,n=01‍uklmnEkmξElnη,
where uklmn={u(xk,yl),∂u(xk,yl)/∂x,∂u(xk,yl)/∂y,∂2u(xk,yl)/∂x∂y} and we have suppressed the superscript e on the element shape functions for notational simplicity. Details are given in the Appendix.This tensor-product basis ensuresC1-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 x and evaluating uxx,y on side x=xk yields a cubic function uxy on this side that is uniquely determined by the function value ux and its derivative uxy 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 globalbasis functions ϕix,y,  i=1,2,…,q:
(27)u~x,y;t=∑j=1qαjt·ϕjx,y.
In this paper, the candidate functions are cubic B-splines ((23) in conjunction with generalized coefficients αj) or Hermite polynomials ((25) in conjunction with kinematic degrees of freedom uklmn). It should be clarified that the variable u~ 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ϕj involved in (27), we can apply the well-known Galerkin-Ritz method.Without loss of generality, the boundary consists ofn~1 breakpoints (which correspond to n1 control points) under Dirichlet and n~2 ones (which correspond to n2 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(28)Ma¨+Ca˙+Ka=ft,
where M is the mass matrix, C is the damping matrix, and K is the stiffness matrix, with
(29)M=∫ΩNTρNdΩ,K=∫ΩLNTρLNdΩ.
After imposing the boundary conditions, the dimensions of each matrix in (28) become neq×neq, where neq is thenumber of equations (i.e., the number of unrestricted degrees of freedom: coefficients or kinematic DOFs). For a certain choice of nx×ny subdivisions, neq 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 ofnx×ny subdivisions (cells), thus leading to nx+1×ny+1 nodes.Concerning the determination of the proper Gaussian quadrature, we start from the observation that the elementsmij of the mass matrix are products of two basis functions, each of piecewise pth (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.


## 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 globalbasis functions ϕix,y,  i=1,2,…,q:
(27)u~x,y;t=∑j=1qαjt·ϕjx,y.
In this paper, the candidate functions are cubic B-splines ((23) in conjunction with generalized coefficients αj) or Hermite polynomials ((25) in conjunction with kinematic degrees of freedom uklmn). It should be clarified that the variable u~ 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ϕj involved in (27), we can apply the well-known Galerkin-Ritz method.Without loss of generality, the boundary consists ofn~1 breakpoints (which correspond to n1 control points) under Dirichlet and n~2 ones (which correspond to n2 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(28)Ma¨+Ca˙+Ka=ft,
where M is the mass matrix, C is the damping matrix, and K is the stiffness matrix, with
(29)M=∫ΩNTρNdΩ,K=∫ΩLNTρLNdΩ.
After imposing the boundary conditions, the dimensions of each matrix in (28) become neq×neq, where neq is thenumber of equations (i.e., the number of unrestricted degrees of freedom: coefficients or kinematic DOFs). For a certain choice of nx×ny subdivisions, neq 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 ofnx×ny subdivisions (cells), thus leading to nx+1×ny+1 nodes.Concerning the determination of the proper Gaussian quadrature, we start from the observation that the elementsmij of the mass matrix are products of two basis functions, each of piecewise pth (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 functionsNi,p 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,uexact, the quality of the numerical solution u~ is evaluated in terms of the relative error, which was calculated as follows:
(30)er=u~-uexactuexact×100%.Example 1 (cantilever beam).
Consider a cantilever beam subject to end load as shown in Figure2. The following parameters were used: E=3.0×107 kPa; ν=0.30, D=12 m, L=48 m, P=1000 kN, and plane stress conditions.
The solution is given by Timoshenko and Goodier [14] as
(31)ux=-Py6EI32L-xx+2+νy2-D24,uy=P6EI3νy2L-x+4+5νD2x4+3L-xx2,
where I is the moment of inertia and for a beam with rectangular cross-section and unit thickness it is given by
(32)I=D312.
The rectangular domain is uniformly discretized using of nx×ny subdivisions along the x- and y-directions, respectively, as shown in Figure 2 (for ny=1,2). 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.Example1: cantilever beam subject to a parabolically distributed shear force P.
(a)
(b)
(c)
(d)
(e)Figure 3
Example1: convergence of flexure at the middle of the vertical side BC (shown in Figure 2).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:E/ρ=104 and ν=0.30.
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), neq.Figure 4
Example2: convergence diagram of the first six calculated eigenvalues.Example 3 (clamped plate in plane stress).
Consider a rectangular plate of dimensionsL×D=48×12m of uniform thickness under plane stress conditions, which is clamped along one of its smallest sides. The following parameters were used: E/ρ=104 and ν=0.20. Keeping a constant ratio of subdivisions (nx:ny=4) the plate is discretized into nx×ny 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.Figure 5
Example3: convergence diagram of the first six calculated eigenvalues.


## 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 C1-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 nx×ny subdivisions along the x- and y-directions, respectively, and considering local approximation within each of the aforementioned nx×ny 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 C1-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 C1-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 olderC1-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.



---

*Source: 754561-2014-11-13.xml*