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

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 C-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.


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 C 1continuity 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 C 1 -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, ℎrefinement of knot insertion, -refinement of order elevation, and the unique -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][13][14][15][16][17].In the field of plate analysis employing NURBSbased 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 Section 2; 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.

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 where  is the order of the B-spline and  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 ( + 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]: If each B-spline basis function  , () is given a weight   , respectively, NURBS basis function    () then can be defined as A NURBS curve can be constructed by using NURBS basis functions    () and control points P  : Given two knot vectors Ξ = { 1 ,  2 , . . .,  ++1 }, H = { 1 ,  2 , . . .,  ++1 } and control points P , , a NURBS surface is defined as where  , , (, ) is the tensor product of NURBS basis functions in two parametric dimensions  and : 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: Equation (8a) shows the mapping from the parametric space to the physical space, where  is the parametric coordinate, x is the physical coordinate, Φ  is the basis function and P  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.u  denotes the value of the unknown field at control point or node P  , and u ℎ 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.

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 where  0 and  0 are, respectively, the axial and transverse displacements of a point on the beam reference plane;  denotes the rotation of the cross section about  axis, and ℎ is beam thickness.
The linear strains with respect to the displacement field are When applying TSDT theory to beams, three degrees of freedom, , , and  0 , are considered.Equation (11) can then be transformed into the following form: where L is the differential operator matrix.Static discrete equation is as follows: The global stiffness matrix K is defined as in which where Φ  are the NURBS basis functions.D is a matrix of material constants: is elastic modulus, and  is shear modulus.
The global force matrix F is defined as when there is only transverse uniform pressure applied all along the beam.

Isogeometric Formulation for
Plates.The displacement field of a plate based on TSDT is as follows [7]: where ( 0 , V 0 ,  0 ) are the displacements of a point on the plate reference plane; (  ,   ) are the rotations of the normal to the reference plane with respect to and -axes, respectively; ℎ is the plate thickness.When applying TSDT theory to plates, five degrees of freedom, ,   ,   ,  0 , and V 0 are considered.The linear strains with respect to the displacement field are where u = [      0 V 0 ] T and the differential oper- The discretization method of plate displacement field is similar to that of beams.Static discrete equation is as follows: The global stiffness matrix K is defined as in which D is a matrix of material constants: is elastic modulus, and ] is poisson ratio.The global force matrix F is defined as if there is only transverse uniform pressure applied all over the plate.

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 where n = (  ,   ) 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: Then, ( 26) can be changed into where the differential operator matrix L  is By imposing essential boundary conditions using penalty method, the final static discrete equation is as follows: Matrix K is obtained by using in which is a diagonal matrix of the penalty parameter.In this paper, the penalty parameter is set to 10 4 , where  is elastic modulus [22].

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 modulus  = 10 4 MPa and poisson ratio V = 0.3.Numerical results are obtained using full ( + 1) Gauss integration points in beam analysis or ( + 1) × ( + 1) Gauss integration points in plate analysis, where  and  are the orders of the NURBS basis functions in parametric directions  and , respectively.

Beams under a Uniform Pressure.
A simply supportedsimply supported (SS) and a clamped-free (CF) beam subjected to a uniform pressure  = −1 N is considered.The length of the beam is  = 1 m and thickness ℎ = 0.1 m.
Boundary conditions are, respectively, for a simply supported end and 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.
Figure 1 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.

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.

Square Plate under a Uniform
Pressure.We consider a square plate with length  = 1m and subjected to a uniform load  = −1 N.Both simply supported and clamped boundary conditions are analyzed.
For simply supported plate, the boundary conditions are Penalty method proposed in Section 3.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 Figure 2, 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 ℎrefinement process.It can also be noted that increasing the order of NURBS basis functions can be very effective to alleviate shear locking.

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 be  = 1m and the uniform load  = −1 N.Both simply supported and clamped boundary conditions for a wide range of aspect ratios from 10 to 10 4 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 by -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.
The deflection shapes of both simply supported and clamped plates with aspect ratio 100 are plotted in Figure 4. Center displacements obtained from different orders of NURBS basis functions are normalized by analytical thinplate solutions [28]; the results are plotted against aspect ratios in Figure 5 for simply supported and clamped boundary conditions, respectively.
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.

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 pressure  = −1 N are considered.The length of the beam is  = 1 m and thickness ℎ = 0.1 m.The stress distributions along the thickness of the beam at location  = 0.3 are plotted in Figure 6 for normal stress component   and shear stress component   , 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.A simply supported square plate subjected to a uniform load  = −1 N is then considered.The length of the plate is  = 1 m and thickness ℎ = 0.1 m.The stress distributions along the thickness of the plate at location  = (0.3, 0.3) are plotted in Figure 7 for normal stress component   and shear stress component   , 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.

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 both ℎ-refinement and -refinement.
In the present approach, we only find out the phenomenon of shear locking and alleviate it by ℎ-refinement and -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][30][31].

Figure 1 :Figure 2 :
Figure 1: Numerical results compared with exact solutions of (a) a SS beam and (b) a CF beam.

Figure 5 :
Figure 5: Normalized center displacements versus aspect ratios for (a) simply supported circular plate and (b) clamped circular plate.

Figure 6 :Figure 7 :
Figure 6: The distributions of (a) normal stress   and (b) shear stress   along thickness in a SS beam at location  = 0.3.

Table 1 :
Normalized tip displacements for extreme thickness cantilever beams under uniform pressure modeled with 10 NURBS elements.