Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 11, N o 2, 2013, pp. 193 - 202 ISOGEOMETRIC STRUCTURAL ANALYSIS BASED ON NURBS SHAPE FUNCTIONS  UDC 515.3 Predrag Milić 1 , Dragan Marinković 2 1 University of Niš, Faculty of Mechanical Engineering, Niš, Serbia 2 Berlin Institute of Technology, Department of Structural Analysis, Berlin, Germany Abstract Finite element method (FEM) is a tool that is mostly used for the structural analysis. The method is based upon approximation of the actual geometry and displacement field of the structure over the finite element domain. In practice, Lagrange polynomials are mostly used as shape functions (approximate functions). They provide only C 0 -continuity at the element boundaries. This paper elucidates a possible approach based on the so-called isogeometric analysis to improve this property. The isogeometric FEM analysis applies the same functions to describe the CAD model, the element geometry and the element displacements. The paper represents developed algorithm and the obtained results for a relatively simple test case. Key Words: Isogeometric Structural Analysis, NURBS, FEM 1. INTRODUCTION A reliable design of sophisticated mechanical systems requires numerous computational analyses and simulations. Hence, the engineering design cannot be considered separately from the analysis. They are strongly related because the process of design optimization is evaluated through the analyses and, vice versa, the analyses drive the design process by revealing design weakness. In the 'classical' finite element method (FEM), the geometric model suitable for structural analysis is usually not the very same CAD geometry model. Most often a simplified geometry is used for generating the mesh of finite elements. Furthermore, the resulting domain of FEM models is only an approximation of the already simplified CAD geometry. These ap- proximations give rise to larger or smaller errors which produce several consequences. For example, the stability of thin-walled structures is very sensitive to geometric imperfections and,  Received December 5, 2013 Corresponding author: Predrag Milić Faculty of Mechanical Engineering, A. Medvedeva 14, 18000 Niš, Serbia E-mail: pmilic@masfak.ni.ac.rs Acknowledgements. This work is supported by the project TR35049 financed by the Ministry of Education and Science of Republic of Serbia. 194 P. MILIĆ, D. MARINKOVIĆ therefore, without the exact geometry, an adaptively formed FE mesh may yield results of insuf- ficient accuracy. The high accuracy of the geometry description is also important in design op- timization. Furthermore, the quality of the FEM geometry descriptions has a significant impact in contact problems. Development of the classical FEM is mainly based on the application of isoparametric finite elements, which implies mapping of finite elements from the real geometry into the so-called mas- ter element, i.e. from the global coordinate system into the natural coordinate one. These types of elements use the same shape function, typically Lagrange polynomials, to approximate both the displacement field and the geometry. The quality of geometry description directly depends on the degree of the selected polynomial and element size, i.e. the quality of the finite element mesh. One of the main disadvantages of this approach is that the Lagrange polynomials as shape functions provide only C 0 continuity at the element boundaries, while many elements require at least C 1 con- tinuity, such as beams, plates or shells. In relation to the aforementioned characteristics of the classical FEM, a significant step forward would be a FEM formulation that allows for geometry description done in exactly the same manner as in the CAD software packages. Such a FE model would ideally represent the analyzed geometry, regardless of the mesh quality. This FEM concept is known as isogeometric analysis [1, 2, 5]. As basis for setting the CAD model, different technologies of computational geometry can be used. Currently in the engineering design, NURBS (non-uniform rational B- spline) is mostly used. The main advantage of NURBS is a convenience for representing free forms, accurate representation of conical sections (circles, cylinders, spheres, ellipsoids), etc. The continuity of p-order NURBS is C p+1 . As opposed to B-spline and NURBS surfaces, a T- spline represents a surface based on the points of a control polygon and not on the control mesh [3, 4]. T-spline is non-uniform B-spline surface with T-junctions (for more information see [3]). T-splines are also used to merge multiple patches into one tight part. Immersed boundary meth- ods, which do not require body-fitted meshes, may also be applied within the isogeometric ap- proach. This allows a structural analysis of rather complex engineering parts [6]. 2. NURBS GEOMETRY Definition of complex geometric shapes is usually performed by using splines. Splines application in geometric modeling is related to the work of French engineers Pierre Bezier and Paul de Faget in the automotive industry in the sixties (Fig. 1). Fig. 1 Example of drawing a spline with lattice ISOGEOMETRIC STRUCTURAL ANALYSIS BASED ON NURBS SHAPE FUNCTIONS 195 As the basic type of spline for geometric modeling, B-splines and splines that use the same basis as B-splines are applied. This primarily refers to NURBS and T-splines, the latter representing a generalization of NURBS and also a part of the so-called PB spline (point based spline). A p-order NURBS curve can be represented as a rational function using the Bernstein basic functions by means of the following expression [1, 2, 8]: bua wuN PwuN uC n i ipi n i iipi      0 , 0 , )( )( )( (1) where Pi are the control points that form a control polygon, wi are the weights, and {Ni,p(u)} are the B-spline basic functions of order p defined on non-uniform knot vector U: },...,,,...,,,...,{ 11 bbuuaa pmp  U (2) Usually the knot vector is normalized, i.e. a = 0 and b = 1. Repetition of elements a and b in the knot vector depends on the order of spline and it is p+1. In this way, the discontinuity is realized at the ends of the spline. Rearranging Eq. (1), the NURBS curve can be represented by the following expression, where Ri,p are the basic rational functions of the p-order NURBS:      n i ipii n i n j jpj ipi PuRP wuN wuN uC 0 , 0 0 , , )( )( )( )( (3) Surfaces and solids can be defined in a similar way. A NURBS surface of order p in the u-direction and order q in the v-direction can be expressed as: 1,0 )()( )()( ),( 0 0 ,,, 0 0 ,,        vu wvNuN PwvNuN vuS n k m l lkqlpk n i m j ijijqjpi . (4) By substituting the basic rational function: 1,0 )()( )()( ),( 0 0 ,, ,,     vu wvNuN wvNuN vuR n i m j ijqjpi ijqjpi ij , (5) the equation for the surface takes the form:     n i m j ijij PvuRvuS 0 0 ),(),( . (6) The B-spline basic functions are used to determine the NURBS basic functions and they are determined by means of the Cox-de Boor's recursive formula and using the basic func- 196 P. MILIĆ, D. MARINKOVIĆ tions of lower order. For any given knot span, at least p+1 B-splines and NURBS basic functions are not equal to zero. Also, the basic function has values different from zero at p+1 knot spans. In Fig. 2, a NURBS surface is shown in the so-called physical space, index space and parameter space together with the corresponding basic functions of B-spline. The area defined by index coordinates ξ and η is known as a patch. A less complex geometry can be described by a single patch [7], while more complex geometric forms need to use more patches. Repeating the elements knot vector creates an open knot vector and, hence, in mod- eling complex geometries, the C 0 -continuity is obtained on the element borders because only the continuity of the polygon control point displacements is realized. As the NURBS basis functions are constructed from the B-spline basis functions, the de- rivatives of rational functions will clearly depend on the derivatives of their non-rational counterparts as well [1]. In literature, there are several algorithms for computing deriva- tives of the B-spline basic functions. Fig. 2 NURBS surface, index space, and B-spline basic function 3. BASICS ОF ISOGEOMETRIC ANALYSIS The isogeometric analysis is based on the fact that the same shape functions are used to describe the CAD geometry and the displacement field of the element. It implies that there is a set of basic functions, R, which perform geometric mapping from index space ξ to physical space x, i.e. there is parameter ξ such that mapping x: ̂ is realized as follows:                                 e a e a e ae a a z y x R z y x en 1 )(ξ ξ ξ ξ ξx (7) where xa e , ya e , za e are the components of the a th element of the vector that defines element geometry Ωe (the set of control points). In a similar way, other connections between the parameter and the physical space can be formed, such as connections between control point displacements di and the displacements within the element, i.e. R h  ˆ:)(û : ISOGEOMETRIC STRUCTURAL ANALYSIS BASED ON NURBS SHAPE FUNCTIONS 197 i n i i h np R d)()(û 1    ξξ (8) For a planar problem, the previous equation takes the form:                    enn i ee ii enen enen R vRvR uRuR v u 1 )()( 11 11 ... ... Rddu (9) where u and v are the components of the element displacement obtained by means of basic functions Ri and the corresponding displacements at control points di, and d (e) are the displacements of all the element control points. Assuming small displacements (linear approach), strains ε are further given as:               xy y x ε , where x v y u y v x u xyyx             ,, (10) Respectively, using the previously presented concept: )( 1 )( 1 e n i e ii n i i i i i i i i i enen v x R u y R v y R u x R BddBε                                (11) where B is the strain-displacement matrix. The element stiffness matrix is computed by integration over the element volume using the strain-displacement matrix, B, and the Hooke's matrix of the material, D. Element stiffness matrix Kij that relates the i-th and j-th point of the control polygon can be expressed as:        1 1 1 1 )()( ||),(),( ddt e j T i e ij JDBBK (12) where t is the thickness of the structure. Computation of the integral is performed by numerical integration, typically by means of the Gaussian quadratures [9]. In order to compute the system stiffness matrix, partial derivatives of the shape functions are needed so that the strain-displacement matrix and the Jacobi matrix can be determined for an adequate set of the Gaussian integration points. The number of integration points depends on the order of basic functions. For each element, there are at least p+1 main functions different from zero in one direction. This means that, for a solid element with the second-order NURBS basic functions, displacements in the element are determined by means of 27 main functions. This number of basic functions requires significant computational time. As the influence of a basic function spreads not only over the domain of a single element, but over p+1 elements, this means that any effect in one element is reflected in the p neighbouring elements. The consequence of this property is that the system stiffness matrix typically has a wider matrix band compared to the stiffness matrix in the classical finite element method. 198 P. MILIĆ, D. MARINKOVIĆ 4. SOFTWARE SOLUTIONS A new software solution has been developed in order to implement the isogeometric FEM formulation. There are basically two approaches to structural analysis by means of the isogeometric FEM. The first approach implies that the entire structure is described by means of one patch and it comes down to computing the stiffness matrix of a single patch. The second approach means solving a complex structure. In this case, once the stiffness matrices of individual patches are determined, the stiffness matrix of the whole system is assembled. The input data for the developed software are: geometry data (degree of inter- polation functions for each direction, knot vectors for all three directions, mesh of control points with coordinates and weights), material properties, loads and constraints. Fig. 3 The algorithm of the isogeometric analysis using NURBS shape functions For the example discussed in this paper, the assumption of linearly elastic, homogene- ous, isotropic material has been adopted. Loads can be defined as concentrated forces, surface loads (pressure or force components on the surface) and volume loads. The soft- ware implementation of loads depends on the type of load. The basic steps and informa- tion flow in the algorithm are shown in Fig. 3. The resulting system of equations yields the displacements of the polygon control points. The computation of element displacements and the corresponding element stress-strain states requires the values of basic functions and their derivatives at certain points of the element. The stresses are determined based on the known strains and the Hooke's matrix. ISOGEOMETRIC STRUCTURAL ANALYSIS BASED ON NURBS SHAPE FUNCTIONS 199 5. NUMERICAL EXAMPLE The application of the developed FEM formulation and algorithm will be demonstrated on a relatively simple example of a linear static analysis. The results have been compared with the solution obtained using commercial FE software package ANSYS. Fig. 4a depicts the quarter model of a plate with a hole that has been used as a test example. The constraints in the model are set so as to represent the symmetry of the modeled structure. The structure has been modeled as a solid and the depicted load is defined as a surface load of intensity 10 kN/cm 2 . The CAD model of the above described problem (polygon of control points and elements) is given in Fig. 4b as a NURBS patch. Fig. 4 Plate with a hole – a) physical model, b) NURBS geometry The exact analytical solution of the infinite plate with a circular hole under constant in-plane tension exists [10] and it can be represented in polar coordinates (r, θ) as: )2sin(321 2 ),( )2cos(31 2 1 2 ),( )2cos(341 2 1 2 ),( 4 4 2 2 4 4 2 2 4 4 2 2 2 2                                                             r R r R r r R r R r r R r R r R r x r xx xx rr (13) where   rrr ,, are the stress components in the polar coordinate system for corre- sponding coordinates r and θ. The maximum value of normal stress in the direction of tension computed with previ- ous equations yields 30 kN/cm 2 , i.e. the stress concentration factor is K = 3. For a similar case, but with finite plate dimensions, there is no exact analytical solution, but a higher stress concentration factor is expected. As a reference model for evaluating the accuracy of the isogeometric solution, a rather fine FEM model created in Ansys software package has been used (Fig 5.). 200 P. MILIĆ, D. MARINKOVIĆ Fig. 5 FEM model created in Ansys with 44235 quadratic elements (Plate183) The solution convergence has been investigated for this problem by means of isogeometric analysis based on two different geometry definitions. The geometric models and, consequently, the FEM models, differ in the degree of NURBS basic functions and definitions of knot vec- tors. The initial model has been prepared manually, while further model refinements have been done by using originally developed codes for knot insertion and degree elevation. The first al- gorithm increases the number of entries in the knot vector and, hence, the number of finite ele- ments. The second algorithm is used to obtain the NURBS basic functions of higher order. The results of analyses are shown in Table 1 and Fig. 6. Table 1 The dependence of results on the interpolation function degree and mesh size Case Number of elements Order Number of basis functions Translation in x direction [m] 1 2 2 24 -2.54158E-07 2 8 48 -2.69193E-07 3 32 120 -2.74696E-07 4 2 3 48 -2.70553E-07 5 8 80 -2.74405E-07 6 32 168 -2.75427E-07 Fig. 6 The result convergence with quadratic and cubic NURBS shape functions ISOGEOMETRIC STRUCTURAL ANALYSIS BASED ON NURBS SHAPE FUNCTIONS 201 Although the number of elements in the isogeometric model was significantly smaller (128 elements) compared to 44235 elements of the FEM model in ANSYS, the obtained results for displacements (Fig. 7) and stresses (Fig. 8) differ negligibly. This demonstrates the capability of the NURBS shape functions to describe the stress field very accurately even with a small number of elements. The stress concentration factor computed by means of the reference model is K=3.57 and, as anticipated, it is larger than the value for the theoretical case with infinite dimensions. The classical FEM model (the ANSYS model) needed a high number of elements to recover very accurately the high stress gradients in the vicinity of the hole edge. Fig. 7 The displacement field obtained by a) isogeometric structural analysis, b) commercial software package ANSYS (with 267032 DOF) Fig. 8 The normal stress component xx (N/cm 2 ) obtained by a) isogeometric structural analysis, b) commercial software package ANSYS 6. CONCLUSIONS The isogeometric FEM analysis is a relatively new direction of FEM development in the field of structural analysis and a rapid pace of its development is expected in the years to come. The world's leading vendors of CAD software have already included NURBS as the main tool to define and display more complex geometries. The main problem resides in a seamless integration of CAD geometries and FEM models for structural analysis. The 202 P. MILIĆ, D. MARINKOVIĆ isogeometric analysis is a significant step forward in this direction. One of the major disadvantages of NURBS is their inability for local mesh refinement, but this can be overcome by applying a special type of NURBS denoted as T-spline. This will be the focus of the authors' further work. It is expected that the field of structural dynamics will particularly benefit from the advantages offered by the isogeometric FEM analysis. REFERENCES 1. T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering, Volume 194, pp. 4135–4195, 2005. 2. J.A. Cottrell, T.J.R. Hughes, A. Reali, Studies of refinement and continuity in isogeometric structural analysis, Computer Methods in Applied Mechanics and Engineering, Volume 196, Issues 41–44, pp. 4160–4183, 2007. 3. T. W. Sederberg, J. Zheng, A. Bakenov, A. Nasri, T-splines and T-NURCCs, ACM Transactions on Graphics (TOG), Proceedings of ACM SIGGRAPH 2003 TOG Homepage, Volume 22 Issue 3, July 2003, pp. 477-484 4. P. Wanga, J. Xua, J. Denga, F. Chena, Adaptive isogeometric analysis using rational PHT-splines, Computer-Aided Design, 43, pp. 1438-1448, 2011. 5. R. Sevilla, S. Fernández-Méndez,A. Huerta, 3D NURBS-enhanced finite element method (NEFEM), International Journal for Numerical Methods in Engineering, Volume 88, pp. 103–125, 2011. 6. D. Schillinger, L. Dedè, M. Scott, J. Evans, M. Borden, E. Rank, T. Hughes, An isogeometric design- through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Computer Methods in Applied Mechanics and Engineering, Volumes 249–252, pp. 116–150, 2012. 7. E. Cohen, T. Martin, R.M. Kirby, T. Lyche, R. F. Riesenfeld, Analysis-aware modeling: Understanding quality considerations in modeling for isogeometric analysis, Computer Methods in Applied Mechanics and Engineering, Volumes 199, pp. 334-356, 2010. 8. L. Piegl, W. Tiller, The NURBS Book (Monographs in Visual Communication), Springer, 1996, ISBN: 3540615458 9. T.J.R. Hughes, A. Reali, G. Sangalli, Efficient quadrature for NURBS-based isogeometric analysis, Computer Methods in Applied Mechanics and Engineering, Volume 199, pp. 301–313, 2010. 10. W.Pilkey, D. Pilkey, Peterson's Stress Concentration Factors, third edition Wiley, 2008, ISBN-10: 0470048247 IZOGEOMETRIJSKA STRUKTURNA ANALIZA ZASNOVANA NA NURBS FUNKCIJAMA OBLIKA Metoda konačnih elemenata je najčešće korišćen alat u strukturnoj analizi. Metoda se zasniva na aproksimaciji stvarne geometrije i polja pomeranja strukture u domenu konačnih elemenata. U praksi se kao funkcije oblika (aproksimacione funkcije) najčešće koriste Lagrange-ovi polinomi. Oni omogućuju samo C 0 kontinuitet na granici elementa. U ovom radu će biti predstavljeno rešenje za poboljšanje kontinualnosti na granici elementa zasnovano na tzv. izogeometrijskoj analizi. Izogeometrijski pristup u okviru metode konačnih elemenata koristi iste funkcije oblika za definisanje CAD modela, geometrije konačnog elementa i za interpolacione funkcije pomeranja u polju konačnog elementa. Rad predstavlja razvijeni algoritam i dobijene rezultate ze relativno jednostavan test primer. Ključne reči: Izogeometrijska analiza, NURBS, FEM