# NURBS-Based Isogeometric Analysis of Beams and Plates Using High Order Shear Deformation Theory

Authors: Xinkang Li; Jifa Zhang; Yao Zheng
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Isogeometric analysis (IGA) based on nonuniform rational B-splines (NURBS) is developed for static analysis of beams and plates using the third-order shear deformation theory (TSDT). TSDT requires C1-continuity of generalized displacements; quadratic, cubic,  and quartic NURBS basis functions are utilized to satisfy this requirement. Due to the noninterpolatory nature of NURBS basis functions, a penalty method is presented to enforce the essential boundary conditions. A series of numerical examples of thick and thin beams and plates with different boundary conditions are presented. Results are compared with analytical solutions and other numerical methods. First-order shear deformation theory (FSDT) is also employed and compared with TSDT in the stress analysis. The effects of aspect ratios and boundary conditions are discussed as well.


---

## Body

## 1. Introduction

Various kinds of beams and plates have become the most significant applications in modern industries, which makes their analysis very important. However, theoretical analyses of complicated shapes of beams and plates are usually very difficult. Therefore, numerical techniques such as finite element method (FEM) have been developed and achieved remarkable success in the analyses of beams and plates.Based on different assumptions for displacement fields, several theories for the analysis of beams and plates have been developed. The Bernoulli-Euler beam theory and the Kirchhoff plate theory [1], which both can be regarded as fourth-order problems [2], can be applied if beams and plates are thin because transverse shear stresses are not considered. The first-order shear deformation theory (FSDT) which accounts for transverse shear effects [3, 4] can be applied for both moderately thick and thin beams and plates. However, in FSDT, the transverse shear strain is assumed to be constant along the thickness direction, and thus a shear correction factor (SCF) has to be considered. The SCF depends on many factors, such as material coefficients, stacking scheme, plate geometry and boundary conditions; therefore, the evaluation of the SCF remains a subject of research [5].To overcome the limitations of FSDT, a more accurate high order shear deformation theory, especially for thick beams and plates, termed the third-order shear deformation theory (TSDT) is proposed and widely used [6, 7]. TSDT matches the actual situation better, so the shear correction factors do not need to be calculated, and more accurate transverse shear stresses can be yield. However, TSDT requires C1-continuity of generalized displacements, which is not an easy task to achieve in FEM. To overcome this difficulty, several methods including some mesh-free approaches have been used, for example,  the moving least-squares (MLS) method [8] and the radial point interpolation method (RPIM) [9]. In this paper, we demonstrate that such a C1-TSDT formulation can be conveniently achieved by using NURBS-based isogeometric analysis.Isogeometric analysis (IGA) was first proposed by Hughes and coworkers [10, 11]; the basic idea is that basis functions, usually NURBS basis functions, which accurately represent the geometry, are directly used as the interpolation functions of the unknown field variables. Thus, one of its advantages over FEM is that it can provide an accurate geometry during analysis process. Mesh refinement, that is, h-refinement of knot insertion, p-refinement of order elevation, and the unique k-refinement of increasing the continuity of basis functions can be simply accomplished by reindexing the parametric space without interaction with the CAD system. This offers a possibility to build high order continuity basis functions and closes the existing gap between Computer Aided Design (CAD) and Finite Element Analysis (FEA).Since the pioneer work of Hughes et al. described above, IGA has been widely applied to various engineering problems [12–17]. In the field of plate analysis employing NURBS-based IGA, Shojaee et al.  [18] conducted the analysis of thin composite laminates based on the Kirchhoff plate theory; Thai  et al. [19] performed the analysis of composite laminates using FSDT; Tran et al. [20] first studied the behavior of functionally graded material plates based on TSDT. However, numerical examples of these previous researches usually confined themselves into a small limit range of plate thickness, and the effects of TSDT on the analysis of beams have not been discussed yet.This paper focuses on the NURBS-based isogeometric approach on static analysis of beams and plates by using TSDT. Different from the previous researches, primary focuses are placed on exploring the advantages of the present method in the analysis of extremely thin and thick beams and plates in a wide range of aspect ratios and also in the inner stresses analysis. The remainder of the paper is organized as follows: first, B-splines, NURBS, and NURBS-based IGA are introduced in Section2; according to TSDT, isogeometric formulations for both isotropic beams and plates are derived in Section 3; Section 4 presents a series of numerical examples of beams and plates with various geometries under different boundary conditions, and numerical results are compared with analytical solutions and those available data obtained by other numerical methods; and finally, some concluding remarks are given in Section 5.


## 2. B-Splines, NURBS, and NURBS-Based IGA

### 2.1. B-Splines and NURBS

Consider a parametric spaceξ∈[0,1]; a knot vector Ξ is defined as
(1)Ξ={ξ1,ξ2,…,ξn+p+1},
where p is the order of the B-spline and n is the number of basis functions (control points) necessary to describe it. A knot vector is said to be open if its first and last knots are repeated (p+1) times [21]. Open knot vectors are standard in the CAD literature. Unless otherwise stated, open knot vectors are always used.With the given knot vectorΞ, the univariate B-spline basis functions can be constructed by the Cox-de Boor recursion formula [11]:
(2)Ni,0(ξ)={0ifξi≤ξ<ξi+1,forp=0,1otherwise,(3)Ni,p(ξ)=ξ-ξiξi+p-ξiNi,p-1(ξ)+ξi+p+1-ξξi+p+1-ξi+1Ni+1,p-1(ξ),forp≥1.If each B-spline basis functionNi,p(ξ) is given a weight ωi, respectively, NURBS basis function  Rip(ξ)  then can be defined as
(4)Rip(ξ)=Ni,p(ξ)ωiW(ξ)=Ni,p(ξ)ωi∑i^=1nNi^,p(ξ)ωi^.A NURBS curve can be constructed by using NURBS basis functionsRip(ξ) and control points Pi:
(5)C(ξ)=∑i=1nRip(ξ)Pi.Given two knot vectorsΞ={ξ1,ξ2,…,ξn+p+1}, H={η1,η2,…,ηm+q+1} and control points  Pi,j, a NURBS surface is defined as
(6)S(ξ,η)=∑i=1n∑j=1mRi,jp,q(ξ,η)Pi,j,
where Ri,jp,q(ξ,η) is the tensor product of NURBS basis functions in two parametric dimensions ξ and η:
(7)Ri,jp,q(ξ,η)=Ni,p(ξ)Mj,q(η)ωi,j∑i^=1n∑j^=1mNi^,p(ξ)Mj^,q(η)ωi^,j^.


### 2.2. NURBS-Based IGA

IGA utilizes the concept of isoparametric, which refers to the use of the same basis functions for both geometry and unknown field discretization:(8a)x=∑i=1nΦi(ξ)Pi,(8b)uh(x)=∑i=1nΦi(ξ)ui.Equation (8a) shows the mapping from the parametric space to the physical space, where ξ is the parametric coordinate, x is the physical coordinate, Φi is the basis function and Pi is control point in IGA or node in FEM. Equation (8b) indicates that in an isoparametric formulation, the unknown field is approximated by the same shape functions used in constructing the geometry. ui denotes the value of the unknown field at control point or node Pi,  and uh is the approximate field value.In FEM, basis functions chosen to approximate the unknown field are then used to approximate the known geometry. In most cases, these basis functions are Lagrange interpolation functions, which makes geometries only approximated. However, in IGA, NURBS basis functions are first chosen to exactly capture the geometry and then be used to approximate the unknown field. Thus, the geometry can keep exact at the very beginning and during the mesh refinement process.


## 2.1. B-Splines and NURBS

Consider a parametric spaceξ∈[0,1]; a knot vector Ξ is defined as
(1)Ξ={ξ1,ξ2,…,ξn+p+1},
where p is the order of the B-spline and n is the number of basis functions (control points) necessary to describe it. A knot vector is said to be open if its first and last knots are repeated (p+1) times [21]. Open knot vectors are standard in the CAD literature. Unless otherwise stated, open knot vectors are always used.With the given knot vectorΞ, the univariate B-spline basis functions can be constructed by the Cox-de Boor recursion formula [11]:
(2)Ni,0(ξ)={0ifξi≤ξ<ξi+1,forp=0,1otherwise,(3)Ni,p(ξ)=ξ-ξiξi+p-ξiNi,p-1(ξ)+ξi+p+1-ξξi+p+1-ξi+1Ni+1,p-1(ξ),forp≥1.If each B-spline basis functionNi,p(ξ) is given a weight ωi, respectively, NURBS basis function  Rip(ξ)  then can be defined as
(4)Rip(ξ)=Ni,p(ξ)ωiW(ξ)=Ni,p(ξ)ωi∑i^=1nNi^,p(ξ)ωi^.A NURBS curve can be constructed by using NURBS basis functionsRip(ξ) and control points Pi:
(5)C(ξ)=∑i=1nRip(ξ)Pi.Given two knot vectorsΞ={ξ1,ξ2,…,ξn+p+1}, H={η1,η2,…,ηm+q+1} and control points  Pi,j, a NURBS surface is defined as
(6)S(ξ,η)=∑i=1n∑j=1mRi,jp,q(ξ,η)Pi,j,
where Ri,jp,q(ξ,η) is the tensor product of NURBS basis functions in two parametric dimensions ξ and η:
(7)Ri,jp,q(ξ,η)=Ni,p(ξ)Mj,q(η)ωi,j∑i^=1n∑j^=1mNi^,p(ξ)Mj^,q(η)ωi^,j^.


## 2.2. NURBS-Based IGA

IGA utilizes the concept of isoparametric, which refers to the use of the same basis functions for both geometry and unknown field discretization:(8a)x=∑i=1nΦi(ξ)Pi,(8b)uh(x)=∑i=1nΦi(ξ)ui.Equation (8a) shows the mapping from the parametric space to the physical space, where ξ is the parametric coordinate, x is the physical coordinate, Φi is the basis function and Pi is control point in IGA or node in FEM. Equation (8b) indicates that in an isoparametric formulation, the unknown field is approximated by the same shape functions used in constructing the geometry. ui denotes the value of the unknown field at control point or node Pi,  and uh is the approximate field value.In FEM, basis functions chosen to approximate the unknown field are then used to approximate the known geometry. In most cases, these basis functions are Lagrange interpolation functions, which makes geometries only approximated. However, in IGA, NURBS basis functions are first chosen to exactly capture the geometry and then be used to approximate the unknown field. Thus, the geometry can keep exact at the very beginning and during the mesh refinement process.


## 3. Isogeometric Formulations for Beams and Plates Based on TSDT

### 3.1. Isogeometric Formulation for Beams

The displacement field of a beam based on TSDT [7] is
(9)u(x,z)=u0(x)+zθ(x)-4z33h2(θ+∂w∂x),w(x,z)=w0,
where u0 and w0 are, respectively, the axial and transverse displacements of a point on the beam reference plane; θ denotes the rotation of the cross section about y axis, and h is beam thickness.The linear strains with respect to the displacement field are(10){εxxεxz}={∂u∂x∂u∂z+∂w∂x}={∂u0∂x+z∂θ∂x-4z33h2(∂θ∂x+∂2w∂x2)(1-4z2h2)(θ+∂w∂x)}.When applying TSDT theory to beams, three degrees of freedom,w, θ, and u0, are considered. Equation (11) can then be transformed into the following form:
(11){εxxεxz}=[-4z33h2∂2∂x2(z-4z33h2)∂∂x∂∂x(1-4z2h2)∂∂x1-4z2h20]×{wθu0}=Lu,
where L is the differential operator matrix.Static discrete equation is as follows:(12)KU=F.The global stiffness matrixK is defined as
(13)KIJ=∫LBITDBJdL,
in which
(14)BI=LΦI,
where ΦI are the NURBS basis functions.D is a matrix of material constants:
(15)D=[E00G].Eis elastic modulus, and G is shear modulus.The global force matrixF is defined as
(16)F=∫LpΦIdL,
when there is only transverse uniform pressure applied all along the beam.


### 3.2. Isogeometric Formulation for Plates

The displacement field of a plate based on TSDT is as follows [7]:
(17)u=u0+zφx-4z33h2(φx+∂w∂x),v=v0+zφy-4z33h2(φy+∂w∂y),w=w0,
where (u0,v0,w0) are the displacements of a point on the plate reference plane; (φx,φy) are the rotations of the normal to the reference plane with respect to y- and x-axes, respectively; h is the plate thickness.When applying TSDT theory to plates, five degrees of freedom,w, φx, φy, u0, and v0 are considered. The linear strains with respect to the displacement field are
(18)[εxxεyyεxyεxzεyz]T=[∂u∂x∂v∂y∂u∂y+∂v∂x∂u∂z+∂w∂x∂v∂z+∂w∂y]T=Lu,
where  u=[wφxφyu0v0]T and the differential operator matrix L is
(19)L=[-4z33h2∂2∂x2(z-4z33h2)∂∂x0∂∂x0-4z33h2∂2∂y20(z-4z33h2)∂∂y0∂∂y-8z33h2∂2∂x∂y(z-4z33h2)∂∂y(z-4z33h2)∂∂x∂∂y∂∂x(1-4z2h2)∂∂x1-4z2h2000(1-4z2h2)∂∂y01-4z2h200].The discretization method of plate displacement field is similar to that of beams. Static discrete equation is as follows:(20)KU=F.The global stiffness matrixK is defined as
(21)KIJ=∫ΩBITDBJdΩ,
in which
(22)BI=LΦI.D is a matrix of material constants:
(23)D=E1-ν2×[1ν000ν100000(1-ν)200000(1-ν)200000(1-ν)2].E is elastic modulus, and  ν is poisson ratio.The global force matrixF is defined as
(24)F=∫ΩpΦIdΩ,
if there is only transverse uniform pressure applied all over the plate.


### 3.3. Impose Essential Boundary Conditions in IGA

Similar to the moving least square (MLS) or reproducing kernel (RK) mesh-free shape functions, NURBS basis functions generally do not satisfy the Kronecker-Delta property. Special treatments of essential boundary conditions have to be employed taking the ideas in MLS/RK-based mesh-free method. Techniques available include Lagrange multiplier method [18], penalty method [22, 23], and  Nitsche's method [24, 25]. It is worth noting that the MLS/RK shape functions associated with the interior nodes may also cover the boundary nodes with nonzero measures; however, by using the open knot vectors, NURBS basis functions associated with the interior control points vanish at the boundary, thus the NURBS control points can be clearly partitioned into boundary and interior groups, which makes the treatments of essential boundary conditions in isogeometric analysis much easier than in mesh-free method [25].In this paper, we use penalty method to impose boundary conditions.Take a clamped boundary of a plate for example. The boundary conditions are(25)w=0,u0=0,v0=0,φn=0,φs=0,∂w∂n=0,
where n=(nx,ny) is the unit vector normal to the essential boundary and s is the unit vector perpendicular to n. φn and φs can be decomposed into the following equations:
(26)φn=nxφx+nyφy,φs=nxφy-nyφx.Then, (26) can be changed into
(27)[wφnφsu0v0∂w∂n]T=Lbu=0,
where the differential operator matrix Lb is
(28)Lb=[100000nxny000-nynx000001000001∂∂n0000].By imposing essential boundary conditions using penalty method, the final static discrete equation is as follows:(29)(K+K~)U=F.MatrixK~ is obtained by using
(30)K~=∫ΓbΨITαΨJdΓ,
in which
(31)ΨI=LbΦI.α is a diagonal matrix of the penalty parameter. In this paper, the penalty parameter is set to 104E, where E is elastic modulus [22].


## 3.1. Isogeometric Formulation for Beams

The displacement field of a beam based on TSDT [7] is
(9)u(x,z)=u0(x)+zθ(x)-4z33h2(θ+∂w∂x),w(x,z)=w0,
where u0 and w0 are, respectively, the axial and transverse displacements of a point on the beam reference plane; θ denotes the rotation of the cross section about y axis, and h is beam thickness.The linear strains with respect to the displacement field are(10){εxxεxz}={∂u∂x∂u∂z+∂w∂x}={∂u0∂x+z∂θ∂x-4z33h2(∂θ∂x+∂2w∂x2)(1-4z2h2)(θ+∂w∂x)}.When applying TSDT theory to beams, three degrees of freedom,w, θ, and u0, are considered. Equation (11) can then be transformed into the following form:
(11){εxxεxz}=[-4z33h2∂2∂x2(z-4z33h2)∂∂x∂∂x(1-4z2h2)∂∂x1-4z2h20]×{wθu0}=Lu,
where L is the differential operator matrix.Static discrete equation is as follows:(12)KU=F.The global stiffness matrixK is defined as
(13)KIJ=∫LBITDBJdL,
in which
(14)BI=LΦI,
where ΦI are the NURBS basis functions.D is a matrix of material constants:
(15)D=[E00G].Eis elastic modulus, and G is shear modulus.The global force matrixF is defined as
(16)F=∫LpΦIdL,
when there is only transverse uniform pressure applied all along the beam.


## 3.2. Isogeometric Formulation for Plates

The displacement field of a plate based on TSDT is as follows [7]:
(17)u=u0+zφx-4z33h2(φx+∂w∂x),v=v0+zφy-4z33h2(φy+∂w∂y),w=w0,
where (u0,v0,w0) are the displacements of a point on the plate reference plane; (φx,φy) are the rotations of the normal to the reference plane with respect to y- and x-axes, respectively; h is the plate thickness.When applying TSDT theory to plates, five degrees of freedom,w, φx, φy, u0, and v0 are considered. The linear strains with respect to the displacement field are
(18)[εxxεyyεxyεxzεyz]T=[∂u∂x∂v∂y∂u∂y+∂v∂x∂u∂z+∂w∂x∂v∂z+∂w∂y]T=Lu,
where  u=[wφxφyu0v0]T and the differential operator matrix L is
(19)L=[-4z33h2∂2∂x2(z-4z33h2)∂∂x0∂∂x0-4z33h2∂2∂y20(z-4z33h2)∂∂y0∂∂y-8z33h2∂2∂x∂y(z-4z33h2)∂∂y(z-4z33h2)∂∂x∂∂y∂∂x(1-4z2h2)∂∂x1-4z2h2000(1-4z2h2)∂∂y01-4z2h200].The discretization method of plate displacement field is similar to that of beams. Static discrete equation is as follows:(20)KU=F.The global stiffness matrixK is defined as
(21)KIJ=∫ΩBITDBJdΩ,
in which
(22)BI=LΦI.D is a matrix of material constants:
(23)D=E1-ν2×[1ν000ν100000(1-ν)200000(1-ν)200000(1-ν)2].E is elastic modulus, and  ν is poisson ratio.The global force matrixF is defined as
(24)F=∫ΩpΦIdΩ,
if there is only transverse uniform pressure applied all over the plate.


## 3.3. Impose Essential Boundary Conditions in IGA

Similar to the moving least square (MLS) or reproducing kernel (RK) mesh-free shape functions, NURBS basis functions generally do not satisfy the Kronecker-Delta property. Special treatments of essential boundary conditions have to be employed taking the ideas in MLS/RK-based mesh-free method. Techniques available include Lagrange multiplier method [18], penalty method [22, 23], and  Nitsche's method [24, 25]. It is worth noting that the MLS/RK shape functions associated with the interior nodes may also cover the boundary nodes with nonzero measures; however, by using the open knot vectors, NURBS basis functions associated with the interior control points vanish at the boundary, thus the NURBS control points can be clearly partitioned into boundary and interior groups, which makes the treatments of essential boundary conditions in isogeometric analysis much easier than in mesh-free method [25].In this paper, we use penalty method to impose boundary conditions.Take a clamped boundary of a plate for example. The boundary conditions are(25)w=0,u0=0,v0=0,φn=0,φs=0,∂w∂n=0,
where n=(nx,ny) is the unit vector normal to the essential boundary and s is the unit vector perpendicular to n. φn and φs can be decomposed into the following equations:
(26)φn=nxφx+nyφy,φs=nxφy-nyφx.Then, (26) can be changed into
(27)[wφnφsu0v0∂w∂n]T=Lbu=0,
where the differential operator matrix Lb is
(28)Lb=[100000nxny000-nynx000001000001∂∂n0000].By imposing essential boundary conditions using penalty method, the final static discrete equation is as follows:(29)(K+K~)U=F.MatrixK~ is obtained by using
(30)K~=∫ΓbΨITαΨJdΓ,
in which
(31)ΨI=LbΦI.α is a diagonal matrix of the penalty parameter. In this paper, the penalty parameter is set to 104E, where E is elastic modulus [22].


## 4. Numerical Examples

In this section, static behaviors of isotropic beams and plates are studied based on the TSDT. Unless otherwise stated, the material used in the numerical examples is isotropic with elastic modulusE=104 MPa and poisson ratio v=0.3. Numerical results are obtained using full (p+1) Gauss integration points in beam analysis or (p+1)×(q+1) Gauss integration points in plate analysis, where p and q are the orders of the NURBS basis functions in parametric directions ξ and η, respectively.


### 4.1. Beams under a Uniform Pressure

A simply supported-simply supported (SS) and a clamped-free (CF) beam subjected to a uniform pressureP=-1 N is considered. The length of the beam is L=1 m and thickness h=0.1 m.Boundary conditions are, respectively,(32)w=0,u0=0
for a simply supported end and
(33)w=0,u0=0,θ=0,∂w∂n=0,
for a clamped end. Penalty method proposed in Section 3.3 is used to impose boundary conditions.Two types of numerical analyses are conducted, one is IGA based on TSDT with quadratic NURBS basis functions; the other is FEA based on FSDT with quadratic Lagrange interpolation functions. Both of these two analyses utilize ten uniform elements and no special treatment such as reduced integration method is used.Figure1 shows the comparison of numerical results with their corresponding exact solutions [26] for a SS beam and CF beam, respectively. It can be seen that with the same basis function order and number of elements, IGA based on TSDT can get more accurate results than FEA based on FSDT. It can also be noted that under the same basis function order and number of elements, FEA needs more nodes to represent the geometry than the control points needed in IGA.Numerical results compared with exact solutions of (a) a SS beam and (b) a CF beam.
(a)
(b)


### 4.2. Beams at Extreme Thickness Situations

A test is performed to confirm the capability of the present method for analysis of both extremely deep and extremely thin beams [27]. Two cantilever beams under a uniform pressure are studied, one very deep with an aspect ratio of 4.0, and the other extremely thin with an aspect ratio of 100000. Normalized tip displacements [1] of these two cantilever beams modeled with different orders of NURBS basis functions are presented in Table 1. For deep beams, even quadratic NURBS basis functions can get accurate results, as shear locking is not playing an important role when the beams are deep. However, for thin beams, NURBS elements modeled with low order basis functions suffer from shear locking,  and with the order increases, shear locking is alleviated.Table 1
Normalized tip displacements for extreme thickness cantilever beams under uniform pressure modeled with 10 NURBS elements.
Order of basis functions
Very thick beam(aspect ratio = 4.0)
Very thin beam(aspect ratio = 100000)
Quadratic
0.9954
0.6599
Cubic
0.9976
0.8899
Quartic
0.9981
0.9980


### 4.3. Square Plate under a Uniform Pressure

We consider a square plate with lengthL=1 m and subjected to a uniform load P=-1 N. Both simply supported and clamped boundary conditions are analyzed.For simply supported plate, the boundary conditions are(34)w=0,u0=0,v0=0,φs=0.Penalty method proposed in Section3.3 is utilized to impose essential boundary conditions.A series of moderately thin square plates with aspect ratio of 1000 are examined to observe the convergence characteristics of IGA employing various densities of physical meshes. Quadratic, cubic and quartic NURBS basis functions are considered.  In Figure2, center displacements obtained from different orders of basis functions are normalized with the corresponding analytical Kirchhoff plate solutions [28]; the results are plotted versus the number of elements per side of the plates for simply supported and clamped plates, respectively. Convergence curves indicate that shear locking phenomenon is more obvious in clamped plates than in simply supported plates, and it can be improved during h-refinement process. It can also be noted that increasing the order of NURBS basis functions can be very effective to alleviate shear locking.Normalized center displacements versus the number of elements per side of (a) a simply supported square plate and (b) a clamped square plate.
(a)
(b)


### 4.4. Circular Plate under a Uniform Pressure

To demonstrate the applicability of IGA to model common engineering shapes, a circular plate under a uniform transverse pressure is studied in this numerical example. The radius of the circular plate is assumed to beR=1 m and the uniform load P=-1 N. Both simply supported and clamped boundary conditions for a wide range of aspect ratios from 10 to 104 are analyzed [27]. A 10 × 10 uniform mesh is used in the following numerical examples.Quadratic, cubic, and quartic NURBS basis functions are considered. The cubic and quartic cases are obtained from the quadratic case byk-refinement, which is relatively easy in IGA. Control points and physical mesh of a 10 × 10 uniform meshed circular modeled by quadratic NURBS basis functions are shown in Figure 3. Geometry can keep precise at the very beginning and during the mesh refinement process.A10×10 uniform meshed circular plate: (a) control points and (b) physical mesh.
(a)
(b)The deflection shapes of both simply supported and clamped plates with aspect ratio 100 are plotted in Figure4. Center displacements obtained from different orders of NURBS basis functions are normalized by analytical thin-plate solutions [28]; the results are plotted against aspect ratios in Figure 5 for simply supported and clamped boundary conditions, respectively.Deflection shapes of (a) simply supported circular plate and (b) clamped circular plate.
(a)
(b)Normalized center displacements versus aspect ratios for (a) simply supported circular plate and (b) clamped circular plate.
(a)
(b)For both boundary conditions, at a small aspect ratio of 10, that is, the plate is relatively thick, numerical results are a little larger compared to analytical ones because analytical solutions are based on Kirchhoff assumption of thin plates; however, basis functions with different orders can produce almost the same results. When the aspect ratio is  larger than 100, NURBS elements modeled with low order basis functions suffer from shear locking, and cannot produce precise results. With the order increases, shear locking can be significantly alleviated. It is worth noting that in comparison to simply supported circular plates, the clamped circular plates are much easier to be shear locked for thin plates of same aspect ratio.


### 4.5. Stress Analysis of Beams and Plates

In order to investigate the suitability of the present method for the simulation of inner stresses, stress analysis for both beams and plates is conducted in this section.First a simply supported-simply supported (SS) beam subjected to a uniform pressureP=-1N are considered. The length of the beam is L=1 m and thickness h=0.1m. The stress distributions along the thickness of the beam at location l=0.3L are plotted in Figure 6 for normal stress component σx  and shear stress component τxz, respectively. The results of FSDT are also shown in the pictures along with the corresponding exact solutions [26]. As shown in Figure 6(a), both TSDT and FSDT can predict the normal stress precisely. Figure 6(b) gives the shear stress distributions, and it can be seen that the results of TSDT coincide with the exact solutions very well, and zero shear stresses are obtained at both the bottom and the top surfaces, while the results of FSDT remain constant along the thickness and are not zero at the plate surfaces, which contradict to the real situation.The distributions of (a) normal stressσx and (b) shear stress τxz along thickness in a SS beam at location  l=0.3L.
(a)
(b)A simply supported square plate subjected to a uniform loadP=-1N is then considered. The length of the plate is L=1 m and thickness h=0.1 m. The stress distributions along the thickness of the plate at location l=(0.3L,0.3L) are plotted in Figure 7 for normal stress component σx  and shear stress component τxz, respectively. The results of FSDT are also shown in the pictures along with the corresponding exact solutions [26]. Same conclusions can be obtained just like the case in beam analysis.The distributions of (a) normal stressσx and (b) shear stress τxz along thickness in a simply supported plate at location l=(0.3L,0.3L).
(a)
(b)


## 4.1. Beams under a Uniform Pressure

A simply supported-simply supported (SS) and a clamped-free (CF) beam subjected to a uniform pressureP=-1 N is considered. The length of the beam is L=1 m and thickness h=0.1 m.Boundary conditions are, respectively,(32)w=0,u0=0
for a simply supported end and
(33)w=0,u0=0,θ=0,∂w∂n=0,
for a clamped end. Penalty method proposed in Section 3.3 is used to impose boundary conditions.Two types of numerical analyses are conducted, one is IGA based on TSDT with quadratic NURBS basis functions; the other is FEA based on FSDT with quadratic Lagrange interpolation functions. Both of these two analyses utilize ten uniform elements and no special treatment such as reduced integration method is used.Figure1 shows the comparison of numerical results with their corresponding exact solutions [26] for a SS beam and CF beam, respectively. It can be seen that with the same basis function order and number of elements, IGA based on TSDT can get more accurate results than FEA based on FSDT. It can also be noted that under the same basis function order and number of elements, FEA needs more nodes to represent the geometry than the control points needed in IGA.Numerical results compared with exact solutions of (a) a SS beam and (b) a CF beam.
(a)
(b)


## 4.2. Beams at Extreme Thickness Situations

A test is performed to confirm the capability of the present method for analysis of both extremely deep and extremely thin beams [27]. Two cantilever beams under a uniform pressure are studied, one very deep with an aspect ratio of 4.0, and the other extremely thin with an aspect ratio of 100000. Normalized tip displacements [1] of these two cantilever beams modeled with different orders of NURBS basis functions are presented in Table 1. For deep beams, even quadratic NURBS basis functions can get accurate results, as shear locking is not playing an important role when the beams are deep. However, for thin beams, NURBS elements modeled with low order basis functions suffer from shear locking,  and with the order increases, shear locking is alleviated.Table 1
Normalized tip displacements for extreme thickness cantilever beams under uniform pressure modeled with 10 NURBS elements.
Order of basis functions
Very thick beam(aspect ratio = 4.0)
Very thin beam(aspect ratio = 100000)
Quadratic
0.9954
0.6599
Cubic
0.9976
0.8899
Quartic
0.9981
0.9980


## 4.3. Square Plate under a Uniform Pressure

We consider a square plate with lengthL=1 m and subjected to a uniform load P=-1 N. Both simply supported and clamped boundary conditions are analyzed.For simply supported plate, the boundary conditions are(34)w=0,u0=0,v0=0,φs=0.Penalty method proposed in Section3.3 is utilized to impose essential boundary conditions.A series of moderately thin square plates with aspect ratio of 1000 are examined to observe the convergence characteristics of IGA employing various densities of physical meshes. Quadratic, cubic and quartic NURBS basis functions are considered.  In Figure2, center displacements obtained from different orders of basis functions are normalized with the corresponding analytical Kirchhoff plate solutions [28]; the results are plotted versus the number of elements per side of the plates for simply supported and clamped plates, respectively. Convergence curves indicate that shear locking phenomenon is more obvious in clamped plates than in simply supported plates, and it can be improved during h-refinement process. It can also be noted that increasing the order of NURBS basis functions can be very effective to alleviate shear locking.Normalized center displacements versus the number of elements per side of (a) a simply supported square plate and (b) a clamped square plate.
(a)
(b)


## 4.4. Circular Plate under a Uniform Pressure

To demonstrate the applicability of IGA to model common engineering shapes, a circular plate under a uniform transverse pressure is studied in this numerical example. The radius of the circular plate is assumed to beR=1 m and the uniform load P=-1 N. Both simply supported and clamped boundary conditions for a wide range of aspect ratios from 10 to 104 are analyzed [27]. A 10 × 10 uniform mesh is used in the following numerical examples.Quadratic, cubic, and quartic NURBS basis functions are considered. The cubic and quartic cases are obtained from the quadratic case byk-refinement, which is relatively easy in IGA. Control points and physical mesh of a 10 × 10 uniform meshed circular modeled by quadratic NURBS basis functions are shown in Figure 3. Geometry can keep precise at the very beginning and during the mesh refinement process.A10×10 uniform meshed circular plate: (a) control points and (b) physical mesh.
(a)
(b)The deflection shapes of both simply supported and clamped plates with aspect ratio 100 are plotted in Figure4. Center displacements obtained from different orders of NURBS basis functions are normalized by analytical thin-plate solutions [28]; the results are plotted against aspect ratios in Figure 5 for simply supported and clamped boundary conditions, respectively.Deflection shapes of (a) simply supported circular plate and (b) clamped circular plate.
(a)
(b)Normalized center displacements versus aspect ratios for (a) simply supported circular plate and (b) clamped circular plate.
(a)
(b)For both boundary conditions, at a small aspect ratio of 10, that is, the plate is relatively thick, numerical results are a little larger compared to analytical ones because analytical solutions are based on Kirchhoff assumption of thin plates; however, basis functions with different orders can produce almost the same results. When the aspect ratio is  larger than 100, NURBS elements modeled with low order basis functions suffer from shear locking, and cannot produce precise results. With the order increases, shear locking can be significantly alleviated. It is worth noting that in comparison to simply supported circular plates, the clamped circular plates are much easier to be shear locked for thin plates of same aspect ratio.


## 4.5. Stress Analysis of Beams and Plates

In order to investigate the suitability of the present method for the simulation of inner stresses, stress analysis for both beams and plates is conducted in this section.First a simply supported-simply supported (SS) beam subjected to a uniform pressureP=-1N are considered. The length of the beam is L=1 m and thickness h=0.1m. The stress distributions along the thickness of the beam at location l=0.3L are plotted in Figure 6 for normal stress component σx  and shear stress component τxz, respectively. The results of FSDT are also shown in the pictures along with the corresponding exact solutions [26]. As shown in Figure 6(a), both TSDT and FSDT can predict the normal stress precisely. Figure 6(b) gives the shear stress distributions, and it can be seen that the results of TSDT coincide with the exact solutions very well, and zero shear stresses are obtained at both the bottom and the top surfaces, while the results of FSDT remain constant along the thickness and are not zero at the plate surfaces, which contradict to the real situation.The distributions of (a) normal stressσx and (b) shear stress τxz along thickness in a SS beam at location  l=0.3L.
(a)
(b)A simply supported square plate subjected to a uniform loadP=-1N is then considered. The length of the plate is L=1 m and thickness h=0.1 m. The stress distributions along the thickness of the plate at location l=(0.3L,0.3L) are plotted in Figure 7 for normal stress component σx  and shear stress component τxz, respectively. The results of FSDT are also shown in the pictures along with the corresponding exact solutions [26]. Same conclusions can be obtained just like the case in beam analysis.The distributions of (a) normal stressσx and (b) shear stress τxz along thickness in a simply supported plate at location l=(0.3L,0.3L).
(a)
(b)


## 5. Conclusion

In this paper, NURBS-based isogeometric analysis (IGA) is developed for the analysis of beams and plates based on the third-order shear deformation theory (TSDT). One of the advantages of the present approach is that it is easy to build high continuity basis functions, which is required in the analysis of beams and plates based on TSDT. Compared to theoretical solutions, the presented numerical results show that IGA based on TSDT can get very accurate results in both thick and thin beams and plates with different boundary conditions. However, when the beams and plates become very thin, numerical results suffer from shear locking and it can be alleviated by bothh-refinement and k-refinement.In the present approach, we only find out the phenomenon of shear locking and alleviate it byh-refinement and k-refinement. However, the computational cost will be increased due to an excessive increment of control points and the high order of NURBS basis functions. Further studies will focus on the locking free IGA numerical methods for beams, plates, and shells [29–31].



---

*Source: 159027-2013-10-07.xml*