Isogeometric Analysis of Curved Beams and Thin Shells Comparing the analysis technique to classic finite element analysis Master’s thesis in Structural Engineering and Building Technology PURIA SAFARI HESARI SARA ALMSTEDT Department of Applied Mechanics CHALMERS UNIVERSITY OF TECHNOLOGY Göteborg, Sweden 2017 MASTER’S THESIS IN STRUCTURAL ENGINEERING AND BUILDING TECHNOLOGY Isogeometric Analysis of Curved Beams and Thin Shells Comparing the analysis technique to classic finite element analysis PURIA SAFARI HESARI SARA ALMSTEDT Department of Applied Mechanics Division of Material and Computational Mechanics CHALMERS UNIVERSITY OF TECHNOLOGY Göteborg, Sweden 2017 Isogeometric Analysis of Curved Beams and Thin Shells Comparing the analysis technique to classic finite element analysis PURIA SAFARI HESARI SARA ALMSTEDT c© PURIA SAFARI HESARI , SARA ALMSTEDT, 2017 Master’s thesis 2017:47 ISSN 1652-8557 Department of Applied Mechanics Division of Material and Computational Mechanics Chalmers University of Technology SE-412 96 Göteborg Sweden Telephone: +46 (0)31-772 1000 Cover: The sculpture "Crest" was installed at the Victoria and Albert Museum in London the summer of 2014. The structural and architectural design is the result of collaborative work between BuroHappold Engineering, Zaha Hadid Architects and LiteStructures. It is used as case study subject in this thesis project. Photograph courtesy of c©Ed Reeve. Chalmers Reproservice Göteborg, Sweden 2017 Isogeometric Analysis of Curved Beams and Thin Shells Comparing the analysis technique to classic finite element analysis Master’s thesis in Structural Engineering and Building Technology PURIA SAFARI HESARI SARA ALMSTEDT Department of Applied Mechanics Division of Material and Computational Mechanics Chalmers University of Technology Abstract New analysis tools must be developed in order to realize designs of increasing geometric complexity. The iterative work flow that is practiced between engineers and designers today is not fully sufficient. Reports from several industries state that 80% or more of analysis time is spent on converting Computer Aided Design (CAD) data to analysis suitable models. Besides being highly time-consuming, this procedure is also prone to generate errors. When going from CAD models to analysis models a change of discretization is performed which introduces geometrical approximations. Thus, a common platform for designers and engineers to work on is needed. A technique that stands out in its ability to truly integrate analysis with CAD is Isogeometric Analysis (IGA). It is a young approach in Finite Element Analysis (FEA) that utilizes the same mathematical basis for description as CAD technology, i.e. Non-Uniform Rational B-Splines (NURBS). In this thesis, the benefits and drawbacks of using IGA for integrated structural and architectural design are investigated. The question posed is whether IGA gives more accurate results and higher convergence rate than classic FEA on a per-degree-of-freedom basis. To this end, isogeometric analysis models based on Timoshenko and Euler-Bernoulli beam theory are implemented and the results are compared with results from classic FEA. Additionally, an isogeometric analysis model based on Kirchhoff-Love shell theory is implemented. For shells, the structural mechanics module of the finite element analysis software COMSOL Multiphysics R©, the software suite Abaqus TM and the programs ETABS R© and Robot TM are used for comparing IGA with classic FEA. The findings from the studies of beams suggest that IGA gives more accurate results. Additionally, the ability to utilize higher continuity over element boundaries in IGA seems to give higher convergence rate compared to classic FEA. The findings from the studies of thin shells also suggest that IGA provides higher accuracy and convergence rate than classic FEA. A mesh refinement method, called k-refinement, that is available in IGA but does not exist in classic FEA proves an efficient way of improving the performance of IGA further. The implication of the results is that IGA outperforms classic FEA in terms of both accuracy and convergence rate. Furthermore, advantages in terms of time savings due to more efficient work flows are also discussed. The basis for discussion is analysis of a free-form thin shell sculpture. Having to return to the CAD model each time changes and mesh refinement is necessary in classic FEA proved to be a tedious task. Thus, the ability to perform adjustments and refinement directly on the isogeometric analysis model is a substantial asset. Keywords: Isogeometric analysis, Finite element analysis, Curved beam elements, Thin shell analysis, Timo- shenko beam theory, Euler-Bernoulli beam theory, Kirchhoff-Love shell theory, Integrated design and analysis i Sammanfattning För att designa strukturer med ökande geometrisk komplexitet måste nya analysverktyg utvecklas. Det iterativa arbetsflödet som praktiseras mellan ingenjörer och designers idag är inte tillräckligt. Rapporter från oberoende industrier konstaterar att 80% eller mer av analystiden spenderas på att konvertera Computer Aided Design (CAD) geometri till modeller som är lämpliga för analys. Denna process är, utöver att vara mycket tidskrävande, även benägen att generera modelleringsfel. Vid övergången från designmodeller till analysmodeller introduceras geometriska approximationer till följd av det byte av diskretisering som krävs för kompatibilitet med analysmjukvaror. En gemensam plattform för designers och ingenjörer att arbeta från är därför nödvändig. En teknik som utmärker sig genom sin förmåga att verkligen integrera analys med CAD är Isogeometrisk Analys (IGA). Det är en ung metod inom Finita Elementmetoder (FEM) som använder sig av samma matematiska bas för att beskriva geometri och basfunktioner som används inom CAD-teknologi, så kallad Non-Uniform Rational B-Splines (NURBS). I detta examensarbete undersöks fördelar och nackdelar med att använda IGA i integrerat konstruktionstekniskt och arkitektoniskt designarbete. Frågeställningen är huruvida IGA ger mer precisa resultat och högre konvergenshastighet än klassisk FEM för samma antal frihetsgrader. I det syftet implementeras isogeometriska analysmodeller för balkar baserade på Timoshenko och Euler-Bernoulli teori. Resultaten jämförs med resultat från klassisk FEM. Därtill implementeras en isogeoemtrisk analysmodell för tunna skal baserad på Kirchhoff-Love teori. För studierna av skal används modulen Structural Mechanics i analysprogrammet COMSOL Multiphysics R© och programmen Abaqus TM , ETABS R© samt Robot TM för att ta fram FEM-resultat att jämföra med. Resultaten från studierna av balkar tyder på att IGA ger högre precision än vad som åstadkoms med klassisk FEM med lika många frihetsgrader. Därtill ger möjligheten att utnyttja högre kontinuitet över elementgränser inom IGA högre konvergenshastighet jämfört med klassisk FEM. Resultaten från studierna av tunna skal tyder också på att IGA ger bättre approximationer och högre konvergenshastighet än klassisk FEM. En metod för nätförfining som finns tillgänglig i IGA men inte i klassisk FEM, så kallad k-refinement, visar sig vara ett effektivt sätt att ytterligare förbättra resultaten med IGA. Slutsatsen är att IGA ger både bättre precision och högre konvergenshastighet än klassisk FEM för samma antal frihetsgrader. Därtill undersöks även tidsbesparingar som kan göras till följd av mer effektiva arbetsflöden över yrkesgränser. Att vara tvungen att återgå till CAD-modellen varje gång förändringar och förfiningar av modellen krävs i klassisk FEM visar sig vara en utdragen process. Möjligheten att direkt införa förändringar och förfiningar i analysmodellen i IGA visar sig därför vara en väsentlig fördel. ii Preface This master’s thesis of 30 credits was carried out as the final part of the master’s program "Structural Engineering and Building Technology" at Chalmers University of Technology. The work took place during the spring of 2017 at the Department of Applied Mechanics. The examiner of the thesis was Prof. Fredrik Larsson. The work was supervised by Prof. Larsson together with Senior Lecturer Dr. Mats Ander. Acknowledgements To our examiner and advisor Prof. Fredrik Larsson. Thank you for steering us in the right direction from day one and for pushing us to excel. Your vast knowledge, especially within the field of shell theory, has been a major asset. We did not always understand the purpose of the tasks you gave us, but it has all come together beautifully in the end thanks to your keen eye. All of this in combination with your endless sharing and unwavering patience with all of our questions has made this project possible. To our advisor and strongest supporter Dr. Mats Ander. Thank you for popping by our cubicle and encouraging us when things did not look so bright. Thank you for your dedicated proof-reading of the report and for your sharp inputs. Your detailed knowledge of beam theory in combination with your commitment to the project has been greatly appreciated. Another special thanks goes out to Engineering Licentiate Vedad Alic from Lund University. We would not have gotten this far without your help. Thank you for debugging code, sending examples and for implementing the theory for handling arbitrarily curved boundaries. Because of your genuine interest in the project it has been a delight to turn to you in times of need. Civil Engineer Rasti Bartek from Cundall is someone else who has our warmest thanks. Having you on board and getting input from an industry perspective from you has been truly valuable. Thank you for providing us with a case study subject, for running finite element analyses on the same and for pushing this project forward. iii iv Contents Abstract i Sammanfattning ii Preface iii Acknowledgements iii Contents v 1 Introduction 1 1.1 Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Basics of Isogeometric Analysis 6 2.1 NURBS as geometric descriptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Defining NURBS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2 Basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.3 Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1.4 Refining NURBS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.5 Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 The isogeometric finite element formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.1 Analysis spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.2 Finite elements with IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.3 Shortcomings of NURBS based IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3 Practical aspects of implementation 21 3.1 Code structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Connectivity matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 Mesh refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.5 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 v 4 Curved Beams 24 4.1 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2 Timoshenko beam theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.3 Euler-Bernoulli beam theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.4 Shear and membrane locking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.5 Section forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.6 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.6.1 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.6.2 Analysis model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5 Thin shells 39 5.1 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.2 Kirchhoff-Love shell theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.2.2 The shell geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.2.3 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.2.4 The Kirchhoff-Love hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.2.5 Small strain kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.2.6 Equilibrium configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.2.7 Elastic constitutive relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.2.8 Galerkin Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.3 Membrane locking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.4 Section forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.5 Boundary conditions in the IGA shell model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.6 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.6.1 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.6.2 Analysis model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 6 Case Study 64 6.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 6.2 Limitations and simplifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 6.3 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 vi 6.5 Workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 7 Conclusions 79 7.1 Accuracy and efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 7.2 Potential for integrated design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 7.3 Suggestions for further work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 References 81 Appendices i A Second derivative of shape functions w.r.t. physical space coordinates i B Method for eliminating shear and membrane locking in classic FEA iii C Change of coordinate system in the beam models v D Change of coordinate system in the shell model vi E Enforcing fixed boundary and symmetry conditions in IGA vii vii viii 1 Introduction 1.1 Context There is a young approach in Finite Element Analysis (FEA). It is called Isogeometric Analysis (IGA) and is a viable alternative to standard polynomial-based FEA. The need for IGA originally sprung from the difference in how designers and engineers work. The models used for design and analysis are not inherently compatible because the analysis programs used by engineers are based on different mathematical discretizations than design software applications. The hinder in work flow over professional boundaries that is created by this leads both to geometric inaccuracy of analysis models and time consuming meshing processes. Interestingly, this gap exists in spite of the fact that analysis is underpinned by geometry and not the other way around. The explanation is that Computer Aided Design (CAD) as we know it today was developed later than classic FEA. However, by changing the way analysis is done through using the same underlying mathematical description as design applications use, the problem is dealt with before it even arises. Thus, the possibility to truly integrate the domains of design and analysis opens up [1], [2]. As stated by Hughes [3], who was the first to introduce the approach back in October 2005: "IGA expands the island of computational techniques to include a method that is both usable and ideal for design applications". Today, isogeometric analysis can play a large role in early schematic engineering work, for example in form-finding. Further down the road it has the potential to completely replace commercial FEA for various applications. Many fields stand to benefit from the method, such as the aviation and automobile industry as well as the building industry. One area that is highly impacted by the new technology is analysis of shells. The new method has made it possible to resume the numerical study of thin shells based on classical Kirchhoff-Love shell theory. It was abandoned early on in finite element formulations due to the high continuity over element boundaries that the theory requires, something that is easily overcome in IGA [3]. In structural engineering, we are particularly interested in shell structures in the context of the built environment. These structures already play an interesting role in the integration of architectural and structural design. The engineer’s understanding of structures is essential because the shapes of shell structures are directly derived from flow of forces. But the architect’s understanding of context and culture is what makes the building an addition to society. Hence, a collaborative effort is essential in order to achieve successful designs. The Hamburg Alster-Schwimmhalle, see Fig. 1.1, is an example of successful collaboration. The building’s roof consists of a hypar shell and the architect required that the edge beams should be cantilevering beyond the facade. According to classical theory for membrane dominated shells this is not possible, something that can be seen when looking at the church of San José Obrero in Fig. 1.2. In order to meet the architectural intent the shell was designed by superimposing the classical saddle surface for the surface loads with the straight line generators’ surface for edge loads. The result is a perfect structure that developed only because of the asserted architect and the skilled engineer [4]. L’Oceanografic is another shell example from the built environment. It is the re-elaboration of a thin shell by Felix Candela. The shape is intricate and it is impressive in its light expression. Fig. 1.3 shows the shell during the construction phase and Fig. 1.4 the finished building. Again, extensive engineering work goes into realizing a shell like this with its remarkable slenderness of 6 cm thickness over a 40 m span [5]. With the increased freedom of design that CAD software enables greater geometrical complexity is pursued. It is evident that in order to enable realization of these increasingly complex designs analysis must keep up. The situation today where meshes are generated from CAD data is not sustainable. Not only does it result in a different geometric description that is only approximate of the original design, it is also highly time-consuming. In several industries it is estimated that 80% or more of the total time modeling and only 20% on the analysis itself. Thus, the need for a technique that takes the exact geometry and effortlessly meshes its curves, surfaces or solids is urgently needed [3]. 1 Figure 1.1: Cantilevering shell of Hamburg Alster-Schwimmhalle. Photograph from [6]. Figure 1.2: Church of San José Obrero by Félix Candela. Note the need for vertical support of the edge beams. Photograph from [7]. 2 Figure 1.3: Roof structure of L’Oceanografic in Valencia during construction. Photograph from [5]. Figure 1.4: L’Oceanografic after completion. Photograph from [5]. 3 For a new technique to succeed commercially the underlying discretisation must be either directly compatible with or easily converted from CAD technology. In isogeometric analysis this is facilitated through use of the same mathematical basis for description as is used in dominating design software, i.e. Non-Uniform Rational B-Splines (NURBS) [8]. Independent of how successful the adoption of IGA eventually becomes for everyday engineering tasks it is a highly current topic with great potential and ideal for efforts of bringing architectural and structural design together. In this report, the basics of NURBS and how they are used in the implementation of isogeometric analysis are explained. Additionally, differences in the implementation of IGA compared to FEA are addressed. It is assumed that the reader has prior basic knowledge of finite element analysis based on the Bubnov-Galerkin method, in this report referred to as classic FEA [9], [10]. It is also necessary to have some understanding of differential geometry of surfaces to fully understand equations and derivations, a quick review of Ch. 2.4. of Kiendl’s thesis [11] is sufficient. After implementation a couple of basic examples are given and then the new-found knowledge is applied to analysis of an existing project. The report should be read as a summary of the thesis project as well as a "how-to" guide for basic implementation of IGA for anyone without previous knowledge of the subject. The main claim investigated is that IGA is more accurate on a per-degree-of-freedom basis compared to classic FEA. To this end IGA is applied to problems within the field of structural engineering. Structural elements that are difficult to analyze with classic FEA are studied and the results are compared with the traditional approach. 1.2 Purpose The purpose of the project can be split into two parts. Part one consists of a thorough guide to implementing IGA. This is something that the authors found lacking while working with the project. Most of the available literature on the topic deals with complex problems which makes it harder to get a comprehensive basic understanding of IGA. To achieve this the simple examples in 1D that were investigated while studying IGA and implementing the necessary algorithms for the final case study are presented in-depth. Part two consists of using the implemented IGA shell model to analyze an existing project. Throughout the project the advantages and disadvantages of using IGA over a standard FE-analysis are studied. For the final case study IGA with Kirchhoff-Love shell theory is implemented. Analysis according to this theory can not be performed using classic FEA but has been shown to be both robust and perform well under large deformations with IGA. Thus, by choosing IGA with Kirchhoff-Love theory the possibility to analyze problems that could not have been dealt with in the same way otherwise opens up [12]. Based on current research the hypothesis is that IGA will prove more time efficient and yield more accurate results than classic FEA on a per-degree-of-freedom basis. The main reasons include elimination of the step of converting geometry from CAD to classic FEA meshes and avoiding the geometrical error that is otherwise introduced by classic FEA basis. The sub-tasks devised to build theoretical and implementation knowledge in part one are: 1. Task: Analyze a simple bar with varying distributed axial load using IGA and classic FEA. Purpose: Master basic NURBS theory and differences in notation compared to classic FEA, as well as basic implementation of IGA (results and implementation are omitted from report). 2. Task: Derive equations for an arbitrarily curved beam in 1D using both Timoshenko and Euler-Bernoulli beam theory and analyze using IGA and classic FEA. Purpose: A first step towards shells through curved elements, handling of locking phenomena in IGA and a hint to whether the hypothesis will be confirmed (see Ch. 4). 3. Task: IGA with Kirchhoff-Love shell theory of benchmark shell problems. Purpose: Study Kirchhoff-Love shell theory and master implementation in 2D as well as verifying the developed model (see Ch. 5). 4 The case study is performed for a thin shell. The geometry consists of a single patch NURBS surface and it is analyzed with Kirchhoff-Love shell theory. The chosen project is called "Crest" and is a sculpture designed by Zaha Hadid architects in collaboration with BuroHappold Engineering and LiteStructures. Deflections and stresses under the impact of its own self-weight are studied. 1.3 Method Review of articles and other literature on the topic is done throughout the thesis work both to build theoretical knowledge and to provide guidance for analysis model decisions. For the subtasks the software Matlab R© is used to perform complex operations and to practice algorithmic thinking in relation to IGA. The main analysis program is built from scratch with some subroutines by supporting supervisor Alic and shape function routines and inspiration from the open source Matlab R© code package "igafem" by Nguyen et al. [8]. The software Rhinoceros R© by Robert McNeel & Associates and the plug-in Grasshopper are used to extract NURBS information through a self-implemented Grasshopper R© component. 1.4 Limitations For isogeometric analysis, the use of any other basis than NURBS is omitted from the project. Research has shown that T-spline basis exhibits better properties than NURBS basis, for example by enabling "watertight" compatibility between multiple patches [13]. However, using T-splines adds an extra level of complexity to the project and the basis is yet to be fully implemented in design software. Additionally, the advantages of using T-splines are not relevant for this project as only single patch geometries are analyzed. Other numerical methods than Gauss point integration are also omitted from the analysis models. Gauss point integration is neither optimal nor the most computationally efficient method for IGA, but it is commonly used for FEA and provides a good foundation for comparing the methods [8]. Additionally, the analyses are limited to linear elastic problems inflicted by point, line and/or body loads. 5 2 Basics of Isogeometric Analysis 2.1 NURBS as geometric descriptions There are three different types of mathematical descriptions of curves and surfaces: the explicit, the implicit and the parametric description. NURBS curves and surfaces belong to the parametric descriptions category. It is the most relevant group of mathematical descriptions when it comes to free-form geometries. NURBS as geometric description is a powerful tool that has become standard in CAD modelling due to its versatility in the ability to represent smooth free form shapes as well as linear shapes, sharp edges, geometric objects such as spheres and more. Since NURBS constitute the basis of isogeometric analysis, the first step towards understanding the analysis method is getting used to NURBS. Any NURBS geometry is defined by four parts: 1. a knot vector, 2. a list of control points, 3. corresponding weights, 4. and a polynomial degree. All of these will be briefly reviewed in the next section. A more detailed explanation is given by Piegl and Tiller [14]. 2.1.1 Defining NURBS Knot vector A knot vector is a set of ascending values corresponding to coordinates in parameter space. The vector is denoted by: Ξ = {ξ1, ξ2, ..., ξn+p+1} where ξi ≤ ξi+1 and ξi ∈ R . (2.1) The subscript i denotes the ith knot, n is the number of basis functions and p denotes the polynomial degree. The knot vector defines the parametric space and its knot intervals divides the same into so called knot spans. The knot spans of non-zero length represent elements in the physical space. Note that the term is not equivalent with elements in classic FEA. If the knots are equally spaced the knot vector is said to be uniform, otherwise it is non-uniform. If the same value appears m times, i.e. ξi = ξi+1 = · · · = ξi+m−1, the knot is said to have multiplicity m. By increasing the multiplicity of a knot the continuity at its location can be reduced. Most commonly knot vectors with multiplicity m = p+ 1 of the first and last knots are used. These are called open knot vectors. They are standard in CAD literature because they allow the start and end point of a curve to be specified by the designer [8],[2],[11]. An example: The NURBS curve shown in Fig. 2.1 (1a) has the knot vector Ξ = {0, 0, 0, 0, 1, 1, 1, 1}. It is uniform because the unique knots {0, 1} are equally spaced. The end knots have multiplicity m = 4. If it has degree p = 3 it is also an open knot vector since the multiplicity of the end knots is m = p + 1 = 4. The curve has one element since there is only one non-zero knot span. The same curve can be represented with more elements as shown in Fig. 2.1 (1b). Then the knot vector becomes Ξ = {0, 0, 0, 0, 0.5, 1, 1, 1, 1}. The unique knots {0, 0.5, 1} define two non-zero spans and thereby two elements. This configuration represents the exact same geometry. Note that knots are different from nodes in classic FEA and should not be used interchangeably. 6 Control points and weights Control points together with their corresponding weights control the NURBS geometry. A linear interpolation between the control points is called the control polygon. When defining the control points and weights of a NURBS geometry the information is usually stored in matrices. For this project the coordinates and weights are structured in a matrix with c = [x y z w] where x, y, z and w are column vectors containing the space coordinates and the weights. The row index of c correspond to the control point index. Control points can be grabbed and pulled and the geometry is thereby changed, see (2a) and (2b) in Fig. 2.1. The new configuration of the curve in (2a) can also be achieved by decreasing the value of the weight instead of changing the location of the control point. If the weight had been increased the curve would have been pulled towards the control point and if it had been lowered the curve would have moved in the opposite direction. The impact of both actions is compact and contained on certain intervals. The specifics will be further discussed in the section about basis functions. However, the ability to locally influence NURBS geometries by changing the position of control points is visualized in (2b). Note that only the left element is affected by the change in the control polygon. When a control point is pulled things move in an monotone fashion without oscillations. This monotone behaviour is critical in design and very important for analysis, furthermore it is not present in classic FEA [8],[2],[11]. Polynomial degree There are different definitions of the polynomial degree. The algorithms in this project are based on the definition p = k − n− 1 where k is the number of knots and n the number of control points. This seems to be the most common definition in IGA. The terms order and degree of NURBS are often used interchangeably in literature but in its strict definition the order equals p+ 1 [15]. Changing the order without changing the knot vector, control points or weights changes the geometry as shown in Fig. 2.1 (3a) [8],[2],[11]. 7 Figure 2.1: (1a): A NURBS curve with one non-zero knot span. (1b): The same NURBS curve with two non-zero knot spans or "elements". Note that the geometry remains the same when the number of elements is increased. (2a): The same curve as in (1a) but with a control point pulled downwards, dark lines show the new configuration. (2b): The same curve as in (1b) but with a control point pulled upwards. Note the local impact on the shape of the curve, the element to the right is not affected. (3a): The same curve as in (1a) but with the order of the curve lowered, the new configuration is shown with dark lines. 8 2.1.2 Basis functions The key to eliminating geometrical errors, as is the main the superior attribute of IGA, is to use NURBS for meshing and as basis functions for the analysis. Non-Uniform Rational B-Splines are, as the name suggests, a generalization of B-splines. A good base for NURBS-discussion is therefore to first address B-splines. In contrast to NURBS, B-splines do not have weights attached to their control points. However, they are defined by the knot vector, polynomial degree and control points just as NURBS. Fig. 2.2 displays the basis functions for a B-spline and a NURBS curve respectively. The only difference between the curves is that the weight of the fourth control point is changed for the NURBS curve. Figure 2.2: Basis functions for a B-spline curve on the left and NURBS basis for the same curve but with increased weight of the fourth control point on the right. B-spline basis functions Given a knot vector, the zeroth order B-spline basis functions are given by: Ni,0 = { 1 if ξi ≤ ξ < ξi+1 . 0 otherwise . (2.2) For polynomial orders p ≥ 1 the basis functions are defined recursively by the Cox-de Boor recursion formula: Ni,p(ξ) = ξ − ξi ξi+p − ξi Ni,p−1(ξ) + ξi+p+1 − ξ ξi+p+1 − ξi+1 Ni+1,p−1(ξ) . (2.3) The first derivative of B-spline basis functions are computed by: d dξ Ni,p(ξ) = p ξi+p − ξi Ni,p−1(ξ) + p ξi+p+1 − ξi+1 Ni+1,p−1(ξ) . (2.3) Basis function properties There are several important properties of the B-spline basis to take note of. Often highlighted attributes include [8], [2], [1]: 1. they constitute partition of unity, sumiNi,p(ξ) = 1 for any p, 2. each basis function is non-negative over the entire domain, Ni,p(ξ) ≥ 0, 9 3. they are linearly independent, 4. the support of a B-spline basis function Ni,p is compact and contained in the interval [ξi, ξi+p+1], 5. they exhibit Cp−m continuity across knot ξ where m is the multiplicity of the knot. Continuity The fifth property, continuity across knots, is visualized in Fig. 2.3 (1a). The knot vector is given on the x-axis and the continuity over element boundaries are given by Cp−m. From the knot vector it can be seen that the degree is p = 4 and that the end knots have multiplicity m = 5, thus the end knots have C4−5 = C−1 continuity. In the same way it can be seen that there is C4−1 = C3 continuity between element 1 and 2. 1. Continuity over element boundaries is, as illustrated, controlled by the knot vector. The ability to easily control continuity is one of the key aspects that for example enable implementation of the Kirchhoff-Love shell theory. What continuity means in terms of the geometry is visualized in Fig. 2.3 (2a)-(2d). For example, the C0 continuity between the fourth and fifth element means that the end of the fourth element has the same position as the start of the fifth element. This is how one introduces kinks in NURBS curves. Note that all basis functions except for one disappear at this location, the one that remains takes the value 1. Where there is C1 continuity, i.e. between the third and the fourth element, the tangent of the curve is the same on both sides of the element boundary in addition to positional continuity. NURBS basis functions NURBS basis functions are constructed from the B-spline basis and NURBS weights through: Rpi = Ni,p(ξ)wi∑n i=1Ni,p(ξ)wi . (2.4) The weights play an important role in defining the basis. Note that the NURBS basis coincides with the B-spline basis if all the weights are equal to each other (following the partition of unity property of B-splines). In other words, B-splines are a special case of NURBS. As one might expect from this, NURBS basis bear much in common with the B-spline basis. The continuity and support follows directly from the knot vector as before. Additionally, the basis still constitutes partition of unity and remains non-negative over the entire domain. Derivatives of B-splines basis functions are efficiently represented in terms of lower order B-spline basis, as follows naturally from its recursive definition. Given that the NURBS basis is constructed from the B-spline basis it also follows logically that NURBS basis derivatives are constructed from B-spline derivatives. Hughes et al. [2] give a good overview of B-spline and NURBS basis derivatives and Piegl and Tiller [14] present an in-depth analysis. Here the equation for the first derivative of the NURBS basis R with regards to the parametric coordinate ξ is given: d dξ Rpi (ξ) = wi W (ξ)N ′ i,p(ξ)−W ′ (ξ)Ni,p(ξ) (W (ξ))2 , (2.5) where N ′ i,p(ξ) is defined as: N ′ i,p(ξ) ≡ d dξ Ni,p(ξ) . (2.6) Furthermore, W ′ (ξ) is given by the sum: W ′ (ξ) = ∑n i=1N ′ i,p(ξ)wi . (2.7) 1Note that continuity at end knots is a matter of definition 10 Figure 2.3: (1a): Continuity over element boundaries in IGA. The basis functions are plotted with the knot vector listed on the horizontal axis. (2a): Discontinuity, C−1. (2b): Positional continuity, C0. (2c): Tangential continuity, C1. (2d): Continuity in curvature, C2. 11 2.1.3 Curves Figure 2.4: A NURBS curve with control points and control polygon displayed on the left. The knots mapped to the physical space (onto the curve) so that the non-zero knot spans divide the curve into elements on the right. B-spline curves B-spline curves are defined as a linear combination of control points and B-spline basis functions. The coefficients of the basis functions are the earlier mentioned control points. Given n basis functions and corresponding control points Bi ∈ Rd with i = 1, 2, ..., n a piece-wise polynomial B-spline curve is given by [2]: C(ξ) = ∑n i=1Ni,pBi . (2.8) NURBS curves Given a B-spline curve and its control points, the control points of the NURBS curve are given by: (Bi)j = (Bw i )j/wi, j = 1, ..., d . (2.9) The control point (Bi)j is the jth component of the vector (Bi) and wi is the ith weight. With Eqs. (2.4) and (2.9) a NURBS curve can be expressed as [2]: C(ξ) = ∑n i=1Ri,pBi , (2.10) this is the NURBS equivalence to Eq. (2.8). Fig. 2.4 shows a NURBS curve with the open non-uniform knot vector Ξ = {0, 0, 0, 0, 1, 2, 2, 2, 3, 3, 3, 3}. The control points are denoted by white circles in the left picture and linear interpolation of the control points give the control polygon. The right picture shows the knots mapped to the physical space as green squares, the line segments represent the elements. 12 2.1.4 Refining NURBS In order to use NURBS as basis for analysis in isogeometric analysis the possibility to refine meshes must exist. The refinement space for B-splines, and by extension NURBS, is rich. Changes can be done to the basis by controlling element sizes, continuity of basis over element boundaries and degree of curves while leaving the geometry both geometrically and parametrically intact. The three ways of refinement in IGA are: h-refinement, p-refinement and k-refinement. h-refinement In IGA, h-refinement is performed by insertion of knots. A new knot can be chosen or an existing one repeated. However, the latter leads to a decrease in continuity of the basis. Note that if an internal knot appears more than p times the curve becomes discontinuous. Given a knot vector as presented in Eq. (2.1) a new knot is defined by ξ̄ ∈ [ξk, ξk+1]. Figure 2.5: NURBS curve before h-refinement (left) and after (right). This gives n+ 1 new basis functions which are calculated using Eqs. (2.2) and (2.3). The new n+ 1 control points are calculated from the original control points through: B̄i = αiBi + (1− αi)Bi−1 , (2.11) where αi =  1 if 1 ≤ i ≤ k − p . ξ̄ − ξi ξi+p − ξi if k − p+ 1 ≤ i ≤ k . 0 if k + 1 ≤ i ≤ n+ p+ 2 . (2.12) An example is shown in Fig. 2.5 where a new knot ξ̄ = 0.5 is inserted in the knot vector Ξ = {0, 0, 0, 1, 1, 1}. The new curve is identical to the original one but more elements and control points are added. This leads to more basis functions of the same order, an enrichment of the base. The strategy is analogous to h-refinement in classic FEA [1]. p-refinement When the polynomial order of basis functions is increased without changing the curve geometrically or parametrically it is called p-refinement. The process involves subdividing the curve by knot insertion, elevating the order of each segment, and removing unnecessary knots before combining the segments again. Piegl and Tiller give an in-depth explanation [14]. The strategy is analogous to p-refinement in classic FEA [1]. 13 k-refinement The k-refinement technique differs from the other two in that there is no analogue refinement method in classic FEA. The order of the curve is elevated from p to q and unique knot values are inserted. This leads to a basis with q − 1 continuity at the new knots, i.e. the highest continuity that can be obtained with NURBS over element boundaries. The operations must be performed in the order described because the process of order elevation and knot insertion is not commutative, see Fig. 2.6. Hughes, Cottrell and Bazilevs [1] state that k-refinement is important and potentially superior for high-precision analysis compared to p-refinement. Figure 2.6: The right column displays order elevation and then knot insertion, the refinement method is called k-refinement and continuity over element boundaries is maintained. It is not commutative and the result from first doing knot insertion and then order elevation is not the same, as seen on the left. 14 2.1.5 Surfaces Figure 2.7: A NURBS surface and its control net. Original configuration on the top and configuration after h-refinement on the bottom. There are four types of NURBS surfaces. The general cylinder, for example, is a NURBS surface consisting of a NURBS curve sweeped a distance d along a vector. The ruled surface is constructed through linear interpolation between two NURBS curves with points of equal parameter values. The surface of revolution is constructed from a NURBS curve which is revolved through an arbitrary angle about an arbitrary axis. Additionally, there are NURBS surfaces that are simply made up from two or more NURBS curves. Again, a more elaborate description is given by Piegl and Tiller [14]. B-spline surfaces B-spline surfaces are defined as a tensor product of two curves. Given that a surface has two dimensions, two knot vectors and their respective polynomial orders are required. What was previously the control polygon is now a so called control point net [1]. Given the knot vectors Ξ = {ξ1, ξ2, ..., ξn+p+1} and H = {η1, η2, ..., ηm+q+1} with polynomial orders p and q, 15 as well as the control net {Bi,j}, i = 1, 2, ..., n, j = 1, 2, ...,m with n×m control points, a B-spline surface is given by: S(ξ, η) = ∑n i=1 ∑m j=1Ni,p(ξ)Mj,q(η)Bi,j . (2.13) The basis function notation Ni,p(ξ) is the same as was used Ch. 2.1.2. The corresponding basis functions for the knot vector H are denoted Mj,q(η). The tensor product nature gives the B-spline surface various properties. For instance, the combined basis is point-wise non-negative and constitutes partition of unity. Furthermore, the local support of basis functions follow directly from the one-dimensional curves that form them. For a surface an element is defined as the non-zero knot spans [ξi, ξi+1]× [ηj , ηj+1] [1]. NURBS surfaces NURBS surfaces are defined in the same way as B-spline surfaces, but with control points as given in Eq. (2.9) and with rational basis functions [2]: Rp,qi,j (ξ, η) = Ni,p(ξ)Mj,q(η)wi,j∑n i=1 ∑m j=1Ni,p(ξ)Mj,q(η)wi,j . (2.14) Giving the NURBS equivalent to Eq. (2.13): S(ξ, η) = ∑n i=1 ∑m j=1R p,q i,j (ξ, η)Bi,j . (2.15) Just as for curves, NURBS surfaces share properties with B-spline surfaces. Fig. 2.7 depicts a NURBS surface with corresponding control net B. Each of the NURBS curves that the surface is constructed from can be refined with the refinement methods described. The bottom picture shows the same NURBS surface but after h-refinement. Note how the control net approaches the surface. 16 2.2 The isogeometric finite element formulation Figure 2.8: An overview of analysis spaces and mapping between them in IGA. So far NURBS in their natural setting have been treated, but only touching upon what they are used for. NURBS are the foundation for bringing geometry and analysis together in IGA. As explained by Hughes, Cottrell and Bazilevs [1], the NURBS basis is part of the isogeometric finite element formulation in the following ways: • to generate a mesh that, for a single NURBS patch, is defined by one or a product of multiple knot vectors, • define elements by subdivisions of the domain created by knot spans, • generate basis functions with local support on a small number of elements, • define the geometry through control points, • through the isoparametric concept, meaning that the same basis is used for geometry and analysis, • and to create refined meshes. 2.2.1 Analysis spaces In IGA, there are three relevant analysis spaces: the physical space, the parametric/index space and the parent space. In Fig. 2.8 the spaces are visualized for a 2D NURBS surface. On the left the surface and the element subdivision is shown as well as the control polygon. The whole patch is transformed at once through geometrical mapping to parameter space, this is in contrast to classic FEA where transformations are done elementwise. Transformation to parent space is done elementwise through affine mapping. Integration by Gaussian quadrature is performed on one element at a time in parent space. The element is then pulled back to the parametric domain where the majority of the analysis takes place, e.g. calculating the basis functions, their gradients and the Jacobian determinant of the pullback. The index space is used to identify nonzero basis functions. The whole patch is then pulled back to the physical space. After solving the 17 system of equations, post-processing is required to obtain results on the geometry itself. Seldom do the control point coefficients coincide with what is sought, especially since control points differ from classic FEA nodes in that they do not lie on the geometry itself [2]. 2.2.2 Finite elements with IGA The isogeometric analysis in this project consist of NURBS based Galerkin finite elements. However, many different numerical solution methods have been explored for IGA and some are better suited than others. Galerkin finite elements are not necessarily the best method, but it is common in classic FEA and by choosing it parallels can more easily be drawn between IGA and classic FEA. Comparing IGA with classic FEA While IGA and classic FEA share a lot, many terms and their meaning differ between the two methods. By being similar but yet different, terms used in IGA can be confusing for the reader educated in FEA. Some of the things that seem similar but lead to faulty assumptions if used interchangeably are presented by Hughes, Cottrell and Bazilevs [2]: IGA Classic FEA Exact geometry Approximate geometry Control points Nodal points Control variables Nodal variables Basis not interpolatory between Basis interpolatory between control points and variables nodal points and variables NURBS basis Polynomial basis High, easily controlled continuity C0 continuity hpk -refinement space hp-refinement space Pointwise positive basis Basis not necessarily positive The isoparametric concept One of the things shared by IGA and classic FEA is the isoparametric concept. The basic idea is to use a basis that exactly models the geometry for the solution space of the numerical method. It is the foundation for IGA and it is common in FEA as well, only with a basis that approximates the geometry instead of one that mirrors it. In other words, in classic FEA the solution fields are imposed on the geometry whereas in IGA the geometry is imposed on the fields [2]: IGA : Geometry → Fields FEA : Geometry ← Fields 18 Numerical solution methods Most of FEA is based on the Bubnov-Galerkin method. The steps are well-known to the reader who has studied finite element methods, starting with defining the weak form and then using Galerkin’s method to construct finite-dimensional approximations by multiplication with a weighting function. This leads to a coupled system of linear algebraic equations the solution space consists of all linear combinations of a given set of basis functions. Beside Galerkin finite elements, collocation, least-square and meshless methods are some examples of numerical methods that have been adopted by IGA. It has been shown that IGA lends advantages over classical FEA to some of these methods. However, since none of these methods are incorporated in this thesis project the theory is not covered either. One issue that needs to be addressed in relation to this is efficient quadrature. Since Gauss point integration is used to solve the linearized system it should be discussed how IGA behaves with different quadrature rules [2]. Efficient quadrature Gaussian quadrature rules are optimal for polynomials in one dimension. They are commonly used for other elements as well even though more efficient spacial rules have been derived for several elements. The study of efficient quadrature rules for NURBS based IGA was initiated by Hughes, Reali and Sangalli [16]. Because IGA includes variable order and variable continuity within a patch numerical procedures are required to determine efficient quadrature rules. The key is to take precise account of the continuity at element boundaries and construct rules that span more than one element. The conclusion by Hughes et al. is the so called "half-thumb-rule" which state that the number of quadrature points should roughly equal half the number of degrees of freedom. This rule differs from other rules that circulate, such as the "equal to degree" or "degree plus one", in that it is independent of the polynomial degree. However, as stated under limitations efficient quadrature has been omitted from this project. It is mentioned here because insufficient quadrature led to oscillations of the results when implementing the theories covered in the coming chapters. Hence, the reader who wishes to follow this thesis for their own implementations should be aware of this. A simple solution is to increase the number of Gauss points when needed. Proof of convergence It is important to be able to prove that NURBS based IGA results in convergent methods. Sufficient conditions for a basic convergence proof for a wide class of problems are satisfied by a basis that is: • C1 on element interiors, • C0 on element boundaries, • and complete. The ease with which continuity over element boundaries can be controlled was discussed previously. Furthermore, it can be shown that the isoparametric concept and the partition of unity property are enough to ensure the third condition, a complete basis, both of which are satisfied in IGA by the NURBS basis as discussed earlier. The proof of several other theorems regarding convergence of NURBS based IGA are discussed by Hughes, Cottrell and Bazilevs [2]. 19 2.2.3 Shortcomings of NURBS based IGA Isogeometric analysis is not the only effort towards unifying the processes of design and analysis. Other efforts to integrate CAD and analysis include for example the finite cell, structured XFEM and immersed boundary methods. One technique that has had some success for shell applications is the use of subdivision surfaces [8]. However, NURBS based analysis remains highly promising in the area of integrated design and analysis. Among the downsides of using NURBS basis is that there is not watertight compatibility between multiple patches and there are difficulties in dealing with trimmed surfaces. The implication is that patches do not fit together tightly geometrically. One attempt at solving this is the use of immersed T-splines for basis. The inclusion of T-splines in CAD applications is quickly advancing. For shells, the local refinement capability of T-splines together with the discretization of Kirchhoff-Love shells shows promise as a highly efficient and accurate formulation. T-splines can reparameterize trimmed surfaces to form analysis suitable untrimmed T- spline representations [13]. Whether designers prefer immersed T-splines over NURBS is not of high importance though because both problems will be solved in IGA. There is also a lot of work to be done in the area of solids. Hierarchical B-splines, immersed T-splines and some boundary element methods are examples of what is being explored [3]. Efficient implementation of IGA is another area that needs to be further developed. It is one of the topics for this years international conference on isogeometric analysis (2017) [17]. Other topics that will be addressed during the conference include industrial applications of IGA, local refinement and adaptivity for IGA and isogeometric spaces over unstructured meshes. 20 3 Practical aspects of implementation 3.1 Code structure Given the similarities between classic FEA and IGA already discussed, it should not come as a surprise that the code structure for both analyses methods are very similar. As can be seen in Fig. 3.1, the main code structure for IGA is created by adding a couple of steps to the classic finite element analysis structure. In the figure these steps are visualized by grey boxes. Figure 3.1: Code structure for IGA. 21 The added steps consist of a different way of reading input data, i.e. data in the NURBS format, slight tweaks to the connectivity arrays used for analysis, a new routine for evaluation of basis functions and a final step where output data is post-processed to obtain results on the geometry evaluated. For the first step, reading input data, it is important to decide on a structure for NURBS data that is consistent throughout the code. Typically, it is convenient to store control points and weights in one matrix while the degree is stored as a separate variable and knot vectors as separate column vectors. However, the format is free as long as it is consistent. To quickly be able to test different geometries it was decided for this project to export NURBS data from Rhinoceros R© through a self-written Grasshopper R© component. This is an effective way for easily importing arbitrary geometries to Matlab R© as the control point matrix and the knot vector quickly increase in size with complex geometries. After refinement to desired mesh properties, connectivity matrices are constructed. The ones used in this thesis project are explained in Ch. 3.2. By looping over elements and quadrature points element stiffness matrices and load vectors are built. Pseudo code for shape function algorithms are stated by Piegl and Tiller [14]. This is the basis for most IGA shape function routines found as free downloads online. The routines by Nguyen et al. [8] in their open-source Matlab R© package "igafem" work well for implementation. Note that there is an inconsistency in their routine as one approaches upper element boundaries which might lead to numerical errors. For a more dependable routine the Matlab R© package by Alic [18] is suggested. It contains shape function routines for B-splines. By reviewing Eq. 2.4. in this report and the examples given in the download the NURBS basis can easily be obtained from these routines. Finally, the system is assembled to global form and solved for sought solution field. The solution field is given for the control points. The results must be post-processed to obtain the sought values since the control points are not attached to the geometry itself in IGA. Theoretically, it works the same way for FE only then the nodal values are located on the approximated geometry. Thus, the step is often omitted from classic FEA. An example of post-processing is given in Ch. 3.5. 3.2 Connectivity matrices There are several different connectivity matrices used in the implementation of IGA. The important thing though is not what they are called or how the data is structured but to understand the connectivity of NURBS. The ones mentioned here are partly used by Hughes [2] and partly used by Nguyen et al. [8]. NURBS coordinates array (INC) INC is a matrix that relates a global function number to uni-variate knot indices. The column index is a global function number and the row index a parametric direction, e.g. ξ or η. Out comes the corresponding NURBS coordinate of the index space. For example, global function A and parametric direction ξ gives index space coordinate i in direction ξ by INC(ξ, A) = i. Element functions array (IEN) IEN relates a global basis function number to its local ordering. It is often referred to the as the "element nodes array" in classic FEA. In IGA it takes the element number e and a local basis function number a and the output is the global basis function number A, e.g. IEN(e, a) = A. The above mentioned are enough for the beam models implemented in this project. However, the information is restructured according to the needs in the shell model. Three different connectivity matrices are used, namely element, index and sctr. In the element matrix the row index is the element number and the output is the index of the control points affecting that element. By inputting the element number in index the index of ξ is obtained from column 1 and η from column, this is used to get element ranges in respective directions. For degree of freedom per element sctr is used, again with row index as element number. 22 3.3 Mesh refinement There are two notions of a mesh in isogeometric analysis. There is the control mesh and then there is the mesh embedded in the geometry. The control mesh can be seen as a scaffold that enables manipulation of the geometry and the design. The embedded mesh consists of images of knot spans mapped onto the geometry. The number of knot spans in the knot vector will mirror the number of elements. Since a new knot can be inserted in an arbitrary location of the knot vector local refinement is possible. Local refinement is not covered in this thesis but it is easy to test by simply adding a knot in one of the knot spans of the knot vector and thereby locally turning one element into two. It is the dual mesh description that is key to implementation. It is true for classic FEA as well but due to the interpolatory basis the meshes coincide. The result is that one works on the control polygon in FEA and not the exact geometry [3]. The refinement methods described in Ch. 2.1.4. are used to refine the analysis mesh. These methods affect both the control mesh and the embedded mesh but not the geometry. Knot insertion corresponds to h-refinement in classic FEA. Additionally, p-refinement is similar in both FEA and IGA. However, k -refinement does not have an equivalent method in classic FEA. Throughout this project h-refinement is used for analyses in 1D and k -refinement for analyses in 2D. 3.4 Boundary conditions Since the dual mesh and the degrees of freedom are located on the scaffold instead of on the geometry, applying boundary conditions might seem diffuse. However, the procedure is very similar to that of finite element methods with control points instead of nodes being restricted. The main thing to take note of is whether the formulations are in local or global coordinates and what these coordinate systems are. Change of coordinate systems is covered in detail in relation to the examples given in Ch. 4 and 5. Something else that complicates how boundary conditions are dealt with is when the formulations are rotation- free. This is also covered in-depth in Ch. 4 for curves and in Ch. 5 for surfaces. 3.5 Post-processing Obtaining the solution fields in terms of node coefficients is enough in classic FEA because the control mesh and the embedded mesh coincide. However, additional steps post the solution are required in IGA because the control points do not intersect with the geometry. The steps begin with choosing a set of coordinates on the geometry in physical space for which solution fields are sought. Then the corresponding values ξ in parameter space are obtained. With ξ known the knot span and connectivity at that location can be found. Thus, control point coefficients for active control points and basis functions at ξ can be computed. Finally, coefficients are multiplied with shape functions and the sought value is given by the sum. Change of coordinate system might have to be applied depending on the problem. Note that this is not novel in any sense, it can be done in classic FEA as well. Example Say one is interested in the displacement in tangential direction at the middle of a NURBS curve representing a loaded curved beam, see subsequent discussion in Ch. 4.6.1 Fig. 4.4. The knot span and connectivity of the corresponding parameter space value is ξ = 0.5 if it is a single element NURBS curve. If the curve has degree p = 2 it will have p+ 1 = 3 control points affecting the geometry at this point. The tangential displacement coefficients for these control points can be extracted as a vector aet where t stands for tangential displacements and e for element in the sense that the vector contains the coefficient that affects that element. Thereafter the shape functions at ξ = 0.5 are computed. The displacement on the physical curve at the midpoint is then given by amidt = N(ξ = 0.5)aet . 23 4 Curved Beams Isogeometric analysis is especially advantageous when the subject for analysis is arbitrarily shaped with strong curvature. In classic FEA geometries are approximated with straight elements, leading to greater geometrical errors originating from the meshing process the larger the curvature is. As has been mentioned previously, this geometric approximation is avoided by using IGA. With the ultimate objective to analyze a shell in mind, implementation of IGA is initiated with routines for curved members in 1D. More specifically, curved beams are chosen as subject for analysis. While practicing implementation of IGA, there is also an interest in investigating whether the size of the geometrical error can be measured. Furthermore, the claim that IGA provides improved accuracy on a per-degree-of-freedom basis compared to classic FEA is studied. Results provide a good indication to whether the hypothesis presented in the introduction is correct before moving on to analysis in 2D and shells. The beam analysis models implemented in this chapter are applicable to linear elastic problems. Analysis models based on both the Euler-Bernoulli and the Timoshenko beam theory are implemented. Firstly, the beam theories are introduced. The theories are based on curved beam elements instead of straight elements, which typically are used in classic FEA. By using curved elements the solutions are brought as close to the exact geometry as possible and potential geometrical errors can be measured more accurately. Secondly, shear and membrane locking is discussed before moving on to an example. The analytical solution to the beam problem chosen is easy to find in literature, it is for example described by Cazzani et al. [19]. Thus, it is easy to verify the model. 4.1 Notations The models are valid for arbitrarily curved beams in one dimension, see schematic sketch in Fig. 4.1. The following restrictions and notations apply: Figure 4.1: Local and global coordinate system definition for arbitrarily curved beams in one dimension. • the centroid line is a planar curve, • the curve is parametrized by arc length s ∈ [0, l], • one of the principle inertia axes and the shear centre of the cross-section lie in the same plane, • the global reference system is denoted by (0 ; x, y), • the local reference system is denoted by (0 ; t, r) where t stands for tangential and r for radial, • and the radius of the curve is denoted by R. 24 4.2 Timoshenko beam theory Figure 4.2: Illustration of assumptions of the Timoshenko beam theory. The undeformed state is displayed on the left and the deformed state on the right. Deformations are exaggerated. The Timoshenko beam theory is based on the following assumptions: • all deformations and strains are small, • Hooke’s law applies for relating stresses and strains, • and planar cross sections remain planar to the deformed axis. The third property is shown in Fig. 4.2. The figure also shows that the Timoshenko beam theory, in contrast to the Euler-Bernoulli beam theory, allows rotation between the cross section and the deformed axis. To include shear deformations leads to a less stiff beam model compared to a model based on Euler-Bernoulli theory. The equilibrium, kinematic compatibility and constitutive relations for a curved beam element based on Timoshenko beam theory are stated in this section. The equations used to derive the finite element formulation are described by Lim and Kitipornchai [20]. Firstly, the equilibrium equations are stated. Here the normal force is denoted N , shear force T and bending moment M . The variable s refers to a coordinate on the beam centroid line. Additionally, distributed forces in tangential and radial direction are denoted qt and qr, whereas distributed moment is denoted m. dN ds + T R + qt = 0 (4.1) dT ds − N R + qr = 0 (4.2) dM ds − T +m = 0 (4.3) Secondly, the kinematic compatibility equations are expressed. Displacements in tangential and radial direction are denoted u and w respectively, whereas rotations are denoted ϕ. The generalized axial ε, shear γ and curvature bending χ strains are then given by: 25 ε = du ds + w R , (4.4) γ = dw ds − u R + ϕ , (4.5) χ = dϕ ds . (4.6) Finally, the constitutive relations are expressed. Young’s modulus is denoted E, the cross-section area A, area moment of inertia I, the shear modulus G and the shear reduction factor κ. N = EAε (4.7) T = GAκγ (4.8) M = EIχ (4.9) The principle of minimum potential energy is used to derive the FE-formulation. Based on Eqs. (4.1)-(4.9) the expression is defined as: u,w, ϕ = arg min û,ŵ,ϕ̂ { 1 2 ∫ L 0 (EAε̂2 +EIχ̂2 +GAκγ̂2) ds− ∫ L 0 (qtû+ qrŵ+mϕ̂) ds ∣∣∣ û(0) = ŵ(0) = ϕ̂(0) = 0 } (4.10) It gives the strain displacement matrix B, stiffness matrix K and load vector f . The variable Ni denotes the shape functions whereas R and s are the radius and coordinate on the centroid line respectively. Additionally, the Jacobian J that maps between parameter and physical space is introduced. In accordance with previous denotations, the variable ξ is a coordinate in the parameter space. The variable q denotes distributed load. B =  dN1 ds N1 R 0 · · · dNn ds Nn R 0 −N1 R dN1 ds N1 · · · −Nn R dNn ds Nn 0 0 dN1 ds · · · 0 0 dNn ds  (4.11) K = ∫ ξ BT EA 0 0 0 GAκ 0 0 0 EI B J dξ (4.12) f = ∫ ξ NTq J dξ (4.13) 26 4.3 Euler-Bernoulli beam theory Figure 4.3: Illustration of assumptions of the Euler-Bernoulli beam theory. The undeformed state is displayed on the left and the deformed state on the right. The classic beam theory, Euler-Bernoulli theory, is based on the same assumptions as Timoshenko beam theory. However, shear strain is neglected and it is assumed that the cross-sections remain normal to the deformed axis during deformation. This affects the kinematic compatibility equations, there is no longer a shear strain part. Additionally, the generalized shear stress T is no longer a part of the constitutive relations. It can be computed from the equilibrium equations if sought. Thus, the same equilibrium equations that were defined for the Timoshenko beam applies: dN ds + T R + qt = 0 , (4.14) dT ds − N R + qr = 0 , (4.15) dM ds − T +m = 0 . (4.16) The kinematic compatibility equations that follow are expressed without the rotation degree of freedom ϕ. There is no generalized shear strain γ, only axial ε and curvature bending χ strains:: ε = du ds + w R (4.17) χ = −d 2w ds2 + 1 R du ds (4.18) Note that the curvature bending strain χ includes the second derivative of the radial displacement w with regards to the centroid line coordinate s. This leads to a finite element formulation that require the second derivative of the shape functions with regards to s, as will be shown further down. 27 The constitutive relations excluding the generalized shear stress are: N = EAε (4.19) M = EIχ (4.20) Again, the principle of minimum potential energy is used to derive the FE-formulation. Based on Eqs. (4.14)-(4.20) the expression is defined by: u,w = arg min û,ŵ { 1 2 ∫ L 0 (EAε̂2 + EIχ̂2) ds− ∫ L 0 (qtû+ qrŵ +m ∂ŵ ∂s ds ∣∣∣ û(0) = ŵ(0) = ∂ŵ ∂s (0) = 0 } (4.21) Just as for the Timoshenko beam, the strain displacement matrix B, stiffness matrix K and load vector f are obtained. The same variable denotations apply. Note the new component d2Ni/ds 2 that is required for implementation of the strain displacement matrix. It has been shown previously that the second derivative with regards to parameter coordinate ξ easily is obtained in isogeometric analysis from the recursive formulation of shape function routines, but the second derivative with regards to physical space coordinate s is more complicated. For a full derivation and further explanation turn to Appendix A. B =  dN1 ds N1 R · · · dNn ds Nn R 1 R dN1 ds −d 2N1 ds2 · · · 1 R dNn ds −d 2Nn ds2  (4.22) K = ∫ ξ BT [ EA 0 0 EI ] B J dξ (4.23) f = ∫ ξ NTq J dξ 4.24) 28 4.4 Shear and membrane locking Geometrical errors are not the only source of inaccuracy in finite element solutions. Another important aspect to investigate is locking phenomena. As described by Bouclier et al. [21] locking occurs in thin beams due to inability of beam finite elements to reflect "shearless bending" and "inextensible bending". The result is an inaccurate answer that reflects much too small deformations and strains. It is the increasing slenderness that give rise to increasing shear-to-bending and membrane-to-bending stiffness that, in turn, introduces errors. Spurious shear and membrane contributions absorb a greater part of the strain energy which results in the tendency of the solution to lock to the transverse shear and/or membrane. Both membrane and shear locking are present in the Timoshenko beam model, while only membrane locking is an issue in the Euler-Bernoulli beam model. Locking phenomena in isogeometric analysis can often be reduced by using a higher degree to describe the NURBS geometry, so called k -refinement. This approach is used in the shell model that is implemented in the next chapter. However, since a close comparison between IGA and classic FEA beam models is the basis of this chapter the implementation is restricted to h-refinement. Since low degrees in isogeometric analysis leads to locking problems in the same way that classic FEA experience locking, methods for relieving low order solutions from locking are implemented. Two different methods to alleviate locking are used, one for the classic finite element analysis and another for the isogeometric analysis. They both consist of projecting shear and bending strains to a lower-order basis, but in the method used with IGA continuity over element boundaries is maintained. Locking in classic finite element analysis To deal with locking in classic FEA a stable discretization of a mixed problem is used. The concept of static condensation is used to solve the system. Practically, this means that the number of degrees of freedom are reduced on element level. Note that if this method is to be applied to isogeometric analysis the reduction must be done on global level since one control point affect multiple elements. Each element has three nodes before the mid node is eliminated. The shape functions living on this node are assigned the properties degree p = 1 and continuity C0. Remaining shape functions are given properties degree p = 2 and continuity C0. See Appendix B for full description. Locking in isogeometric analysis The method chosen to alleviate the IGA beam model from locking is called the B̄ strain projection-type method. As described by Bouclier, Elguedj and Combescure [21] advantage is taken of the distinctive properties of the NURBS basis [21]. More specifically, the approach utilizes the ability to easily control continuity over element boundaries. The basic idea is to produce maximum continuity on a patch and then projecting the shear and bending strain onto a basis of lower dimension. The method is based on a strain projection type method originally used in classic FEA to remove volumetric locking and transverse shear locking in plates. It is equivalent with mixed methods where the associated mechanical variables are searched in the lower dimension space. The discrete form of the mixed problem in the B̄ strain projection-type method is given by: Kf B̄T m B̄T s B̄m Mm 0 B̄s 0 Ms   a αN αT  = f 0 0  condensation ⇐⇒ KB̄ a = f (4.25) where f is the load vector previously defined and the displacement vector a is obtained by solving the system a = K−1 B̄ f . 1 For the Euler-Bernoulli beam model the third block column and the third block row in the stiffness matrix in Eq. (4.25) are omitted. The matrices Mm and Ms are computed with the basis of lower order, here denoted by N̄. This basis is produced from NURBS data with order p− 1 and continuity Cp−1−m. 1Note that K−1 B̄ formally denotes the inverse. However, in practice a linear system of equations is solved and not the inverse. 29 Mm = Ms = ∫ L 0 N̄T N̄ ds (4.26) Just as before s denotes a coordinate on the centroid beam line in physical space and R the radius. The matrix B̄m is obtained from: [B̄m] = [ [B̄1 m] [B̄2 m] [0] ] , (4.27) B̄1 mij = ∫ L 0 N̄T i dNj ds ds , B̄2 mij = ∫ L 0 − N̄i Nj R ds , (4.28) and the matrix B̄s is produced by: [B̄s] = [ [B̄1 s] [B̄2 s] [B̄3 s] ] . (4.29) B̄1 sij = ∫ L 0 − N̄T i Nj R ds B̄2 sij = ∫ L 0 N̄T i dNj ds ds B̄3 sij = ∫ L 0 − N̄T i Nj ds (4.30) The matrix KB̄ is given by the sum KB̄ = K̄m + K̄s + K̄f where the bending related parts of Eqs. (4.11) and (4.22) for the Timoshenko and Euler-Bernoulli theory respectively gives K̄f . The other two stiffness matrices K̄m and K̄s are given by: K̄m = (EA) B̄T mM−1 m B̄m (4.17) K̄s = (GκA) B̄T s M−1 s B̄s (4.18) 4.5 Section forces To complete the beam model the section forces are computed. This section covers how it is done in isogeometric analysis. An easy and quick way to extract them in classic FEA is to use the Matlab R© package CALFEM R©. The method for extracting the forces is a result of the chosen locking handling method. The mechanical variable αN is used to obtain the normal force for the Timoshenko and Euler-Bernoulli beam models. The variable αT is used to obtain the shear force for the Timoshenko beam model. These are given by: αN = (EA) M−1 m B̄T ma , (4.19) αT = (GκA) M−1 s B̄T s a . (4.20) The mechanical variables are used in the same way that control point coefficients are used to obtain displacements of points in physical space. The shape functions of a given evaluation point in parameter space N(ξ) are multiplied with active control point coefficients αN or αT to get the normal and/or shear force at a physical coordinate s(ξ). The moment is computed from the kinematic variable a and constitutive relation M = EIχ. 30 4.6 Example 4.6.1 Problem setup The problem chosen for analysis is a circular beam. More specifically, it consists of a quarter-circle that is fixed at one end and loaded with a point load transversely to the beam at the other end. The problem and the initial unrefined NURBS geometry is shown in Fig. 4.4. Figure 4.4: Quarter-circle beam fixed at one end and loaded with a point load at the other end on the left. Unrefined NURBS geometry on the right. The most basic form of this NURBS geometry is given by the knot vector Ξ = {0, 0, 0, 1, 1, 1}, the degree p = 2 and three control points with corresponding weights. These are the x-coordinates x = [1 1 0], the y-coordinates y = [0 1 1] and the weights w = [1 1/ √ 2 1]. The problem is solved for displacements and section forces. Both beam models are solved using IGA with curved beam elements and classic FEA with straight beam elements. The Timoshenko beam model is also solved with classic FEA using curved beam elements. This enables a comparison of two highly similar set-ups where the first is solved with IGA and the second with classic FEA, without any other differences except for the locking handling method. It is important because a comparison between classic FEA and IGA is the objective of this chapter. Setting up the Euler-Bernoulli beam model is mainly a step towards handling implementation of rotation-free formulations in 2D, which is the next step towards next chapter’s Kirchhoff-Love shell formulations. Therefore, it is sufficient for our purposes to implement curved beam elements in classic FEA only for the Timoshenko beam model. The analytical solution to the problem is given by Cazzani et al. [20]. It is given in terms of an angle θ instead of the centroid line coordinate s as shown in Fig. 4.4. The force P [N ] is the point load applied at the end and R is the radius. The variables E,A, I,G, κ are the material parameters that were previously defined in the beam theory sections. The tangential displacement u is given by: u = −P ( c1(sin(θ) + θcos(θ)) + c2sin(θ) ) , (4.21) c1 = 1 2 ( R EA + R GAκ + R3 EI ) , (4.22) c2 = ( R GAκ + R3 EI ) . (4.23) 31 Additionally, the analytical solution for the generalized stresses N , T and M are obtained from equilibrium: N = −P cos(θ) , (4.24) T = P sin(θ) , (4.25) M = −PR cos(θ) . (4.26) All of these entities are used to verify the beam models. Different heights h of the beam are tested but in general the material parameters used are: Variable Value Unit Young’s modulus, E 80e9 [Pa] Height of beam, h 0.001 [m] Width of beam, b 0.2 [m] Poisson’s ratio, υ 0.2 [-] Radius, R 1 [m] Point load, P −1 [N] The cross-section is rectangular and the area A is given by A = bh while the area moment of inertia I is given by I = bh3/12. The shear modulus G follows analytically for isotropic materials: G = E 2(1 + υ) . (4.27) The shear reduction factor κ is given by [22]: κ = 10(1 + υ) 12 + 11υ . (4.28) Now the only thing left to define are the boundary conditions. For the Timoshenko beam model it is self- explanatory. However, the Euler-Bernoulli beam model is a rotation-free beam formulation and keeping the beam from rotating is a bit more complicated. A degree of freedom on the control point adjacent to the control point at the clamped end must be fixed for the isogeometric analysis. Since the NURBS curve lies in the xy-plane it is the degree of freedom in the x-direction that must be locked into place to prevent rotation. 32 4.6.2 Analysis model Besides using different methods for treating locking the IGA and classic FEA analyses are close to identical. The additional steps required in the isogeometric analysis were explained in Ch. 3. With these added the following code structure of the IGA model for 1D beam analysis is obtained: Algorithm 1: General structure for isogeometric analysis of beams Input : NURBS data in chosen format 1 Load initial NURBS data. 2 Perform h-refinement. 3 Generate connectivity matrices. 4 for i = 1 to number of elements e do 5 Get the knot span of e. 6 Get the connectivity of e, i.e. indices of active control points. 7 Get coordinates of control points influencing e. 8 for j = 1 to number of Gauss points gp do 9 Compute Kf , B̄s, B̄m, Ms and Mm. 10 end 11 Assemble element contributions to global form. 12 end 13 Compute KB̄ , note that this is done with the global matrices. 14 Apply boundary conditions. 15 Apply point load. 16 Solve the system for a = K−1 B̄ F. 17 for k = 1 to number of evaluation points ξ do 18 Compute the section forces as described in Ch. 4.5. 19 Compute displacement in the physical domain by: – applying change of coordinate system as described in Appendix C, – and calculating the displacement as described in Ch. 3.5. 20 end Output : Displacements u,w, (ϕ) and section forces 33 4.6.3 Results Figure 4.5: Undeformed (blue line) and deformed (yellow line) state. Maximum deflection at free end is 0.58 m, this aligns well with analytic solution for both beam models. Verification of IGA beam model The results show that deflections and section forces computed with the Timoshenko and the Euler-Bernoulli beam models align well with the analytical solution. Here, results for the IGA based models are shown. As can be seen in Fig. 4.5 the maximum deflection for given data is 0.58 m. The same value is given by Eq. (4.21). Thus, the displacement in the tangential direction u is verified. The section forces are also compared to the analytical solution, given by Eqs. (4.24)-(4.26). See Fig. 4.6 on the next page for the results. The plots include moment, normal force and shear force diagrams for the Timoshenko beam model. The shear force diagram for the Euler-Bernoulli beam model is omitted for two reasons. Primarily because one can argue whether it actually is of interest to see the shear force plot for a theory that inherently neglects shear strain. Secondly because computing the shear force based on the equilibrium equations, see Eq. (4.16), is a time-consuming and tedious task. It requires the third derivative of the shape functions with regards to the physical space coordinate. In the force plots on the next page, the light blue lines represents the analytical solution. The values given with isogeometric analysis are represented by the dark blue lines. The two lines almost completely overlap. Note that the forces are plotted over a straight line as if the curved beam was unrolled with the clamped end on the left and the free end on the right. 34 Figure 4.6: Moment, normal force and shear force diagram for respective beam model. Dark blue lines represent numerical solution and light blue analytical solution. The distribution over the length of the beam is shown with the clamped end on the left at "0" and the free end on the right at "L". Here L denotes the length of the beam. 35 Comparison between classic FEA and IGA beam models The results include the beam problem analyzed with classic FEM with straight and curved beam elements as well as IGA solved with curved beam elements. Three error plots are produced based on the results. Figure 4.7: Relative error for the Timoshenko beam model. Figure 4.8: Relative error for the Euler-Bernoulli beam model. 36 The first two plots, Fig. 4.8 and 4.7, display the relative error in strain energy on a per-degree-of-freedom basis. The number of degrees of freedom is denoted ndof . The analyses are compared to the analytical solution. The classic FE analysis with straight beam elements is performed with CALFEM R©, an educational finite element Matlab R© package created at Lund University [23]. In the following, we study the strain energy of the solution on the form: E = 1 2 ‖u‖2 , (4.29) ‖ · ‖ = √ a(·, ·) , (4.30) where a(·, ·) denotes the "energy" product pertinent to the left-hand side of the weak form. For a discrete approximation, we then obtain: Eh = 1 2 a(uh, uh) = 1 2 aTKa , (4.31) With h denoting the approximate solution and ref the analytical reference solution, the relative error can be expressed as: erel = ∣∣∣Eref − Eh Eref ∣∣∣ . (4.32) The results obtained with the Timoshenko beam model in Fig. 4.7 clearly show that IGA performs better on both points that are investigated, i.e. more accurate results and higher convergence rate. The relative error is displayed on the vertical axis. The horizontal axis shows the number of degrees of freedom that gives said error. Classic FEA with curved beam elements gives the best approximation for few degrees of freedom, but somewhere around ndof = 10 the isogeometric solution outperforms classic FEA in this regard. IGA also shows a higher convergence rate than classic FEA with both straight and curved beam elements. Note the gap between the FEA curves. Classic FEA with curved beam elements give a more accurate solution for all degrees of freedom compared to classic FEA with straight beam elements. The additional error that comes from using straight elements reflect the geometrical error introduced by approximating a curved geometry with a mesh consisting of straight segments. The results obtained with the Euler-Bernoulli beam model are seen in Fig. 4.8. Again it is clearly shown that IGA continuously gives a better approximation compared to classic FEA with straight beam elements. These results are what was expected and indicate possible fulfillment of the initial hypothesis. The third error plot, see Fig. 4.9, is used to investigate accuracy of the Euler-Bernoulli and Timoshenko beam models with IGA as the height of the beam changes. A height variation gives the slenderness ratio on the horizontal axis. Given on the vertical axis is the difference in converged results for said slenderness of the two beam models. On the left end of the curve, where the beam is thin, the solutions approach each other. For thick beams, towards the right end of the curve, the error in the converged result of the Euler-Bernoulli beam increases. This is expected and verifies the model. The difference between the models is defined as: u∆ = uT − uEB uT , (4.33) where T denotes the Timoshenko beam solution and EB the corresponding Euler-Bernoulli solution. The height of the beam is denoted h and varies from 10 mm to 90 mm. The radius R is set to 1 m. 37 Figure 4.9: Accuracy of end displacement for the Euler-Bernoulli beam model uEB and the Timoshenko beam model uT for varying thickness of beam. 38 5 Thin shells Shell structures are defined by their geometry. In general, they are curved surfaces with a thickness that is very small compared to the other dimensions. The variations in curvature and shape are limitless, making it a popular structural element choice for a wide variety of applications. In architecture shells are found historically in for example masonry vaults. The shell structures designed today are often more organic in their architectural expression. As described by Williams [24], the general strive in structural design of shells is for structures that carry load in pure compression or through a combination of compression and tension. Pure tension structures are often omitted as the term shell implies something relatively rigid. Pure tension structures include for example fabric structures and they are considered a separate category even though they fulfill the geometrical definition of a shell. Shells can also carry load in bending. In that case the thickness is used in order to give the shell bending stiffness. In this thesis the focus is on the study of thin shells. They are, naturally given their slenderness, little or not at all dependent on bending. This leads to a less redundant structure that is more susceptible to geometric imperfections. A heightened geometric accuracy requirement is one of the reasons why isogeometric analysis is especially suitable for analysis of thin shells. 5.1 Notations Notations and variables are described as they appear but throughout the chapter the following applies: • the usual convention for Latin and Greek indices, i.e. i = 1, 2, 3 and α = 1, 2. • the Einstein summation rules. • a comma for partial differentiation, for example f,i = ∂f ∂xi . • subscripts for covariant components. • superscripts for contravariant components. • the subscript 0 for quantities in the reference configuration, for example a point on the reference middle surface is denoted ϕ0. 5.2 Kirchhoff-Love shell theory 5.2.1 Preliminaries Classic Kirchhoff-Love shell theory is widely accepted as the theory for analysis of thin shells. However, the Kirchhoff formulations are not as widespread in classic finite element analysis as formulations based on Reissner-Mindlin theory. This has to do with continuity requirements over element boundaries. In IGA it is easy to control continuity through the knot vector and available refinement methods, a property that becomes highly useful now. The shell model implemented in this chapter is valid for linear elastic analysis of arbitrarily shaped thin shells. Shearing and stretching of the normal to the mid-surface of the shell are omitted and the shell director is restricted to remain normal to the mid-surface in both the undeformed and the deformed state. This is analogous to the Bernoulli assumption about plane cross-sections remaining planar and perpendicular to the line of deflection. 39 The derivations of the discretized model equations follow Kiendl [11] and restrictions to small strain kinematics according to Millán et al. [25] are applied. First, the kinematics of the shell body is formulated. Secondly, the Kirchhoff-Love hypothesis and assumptions are introduced. Through the equilibrium configuration and elastic constitutive relations the weak form is obtained and the linearized system of equations for the shell model is given by Galerkin discretization. To be able to follow the derivations a basic understanding of differential geometry of surfaces is necessary. Bathe [10] provides a thorough description of the geometrical preliminaries. A short introduction to the shell geometry is also given below. 5.2.2 The shell geometry Figure 5.1: Surface in three-dimensional space with curvilinear coordinates, local basis and reference orthonormal basis. Basics A point in three-dimensional space is identified with its coordinates (x1, x2, x3) and unit vectors of a reference orthonormal basis such as (e1, e2, e3): x = x1e1 + x2e2 + x3e3 , (5.1) where x constitutes a position vector. In this example the basis ei denotes the global Cartesian basis, but for describing free-form surfaces it is helpful to use curvilinear coordinates and local bases. Fig. 5.1 shows an example of a surface in three-dimensional space with these coordinates and bases. Two important bases are the covariant basis gi and the contravariant basis gi. Corresponding covariant and contravariant coordinates are denoted ξi and ξi. The position vector expressed in terms of these coordinates x = x(ξ1, ξ2, ξ3) is related to the covariant basis by: dx = dξi gi . (5.2) Thus, the covariant base vectors are defined as: gi = ∂x ∂ξi = x,i . (5.3) The covariant basis is related to the contravariant base vectors by: 40 gi · gj = δji = { 0 i 6= j . 1 i = j . (5.4) A surface is a two-dimensional geometry and thus defined parametrically by two curvilinear coordinates (ξ1, ξ2). For the case of surfaces the first two covariant base vectors gα can thereby be computed as described in Eq. (5.3). The third base vector g3 is defined as a normalized vector orthogonal to g1 and g2: g3 = g1 × g2 ‖g1 × g2‖ . (5.5) From Eq. (5.3) it can be reasoned that the contravariant base vectors gα lie in the tangent plane spanned by the covariant base vectors gα. The third contravariant base vector must thereby be equal to the third covariant base vector g3 = g3. The metric tensor g, also called the identity tensor, is an important quantity for the description of surfaces. It can be expressed in both covariant and contravariant bases: g = gαβ gα ⊗ gβ = gαβ g α ⊗ gβ . (5.6) The metric coefficients gαβ are given by the scalar product of the covariant base vectors: gαβ = gα · gβ . (5.7) This expression is also called the first fundamental form of surfaces and contains important properties of the surface such as length of and angles between base vectors. The contravariant metric coefficients are given by: [ gαβ ] = [ gαβ ]−1 . (5.8) Given these metric coefficient matrices the contravariant base vectors can be computed from the covariant base vectors and the other way around: gα = gαβgα , (5.9) gα = gαβg α . (5.10) The curvature properties of the surface are described by the second fundamental form of surfaces with the curvature tensor coefficients defined by: bαβ = gα · g3,β = −gα,β · g3 . (5.11) 41 Figure 5.2: Mappings between undeformed Ω0 state, deformed Ω state and parametric space A . The shell body The mathematical model of a shell is analyzed as a midsurface with a parameter defining the thickness. Following the notations in Fig. 5.2 the midsurface is denoted Ω0 in the undeformed state, also called the reference state. The thickness is denoted h. The mapping denoted ϕ0 is the mapping between the parametric space of the midsurface, defined as A ⊂ R2, and the undeformed configuration. Furthermore, the unit vector field normal to the midsurface is denoted t0. An area element of the middle surface can be computed as dΩ0 = j̄0 dξ 1 dξ2 where j̄0 = ‖ϕ0,1 ×ϕ0,2‖. The shell body, here denoted ζ, lies in a three-dimensional reference domain such that ζ ⊂ R3. If a configuration x0 is defined as the mapping from the parametric domain A × [−h/2, h/2] into R3 it can be expressed by: x0(ξ1, ξ2, ξ3) = ϕ0(ξ1, ξ2) + ξ3t0(ξ1, ξ2) . (5.12) The body can then be geometrically defined as: ζ = { x0 ∈ R3 ∣∣ x0(ξ1, ξ2, ξ3) = ϕ0(ξ1, ξ2) + ξ3t0(ξ1, ξ2), −h 2 ≤ ξ3 ≤ h 2 , (ξ1, ξ2) ∈ A } (5.13) 42 5.2.3 Kinematics Kinematics are what describe the deformations of the shell body. The deformed state, also called the actual configuration, is defined in the same way as the undeformed state but without the subscript (·)0. The midsurface is denoted Ω and the unit vector field normal to the midsurface t. The mapping between the parametric space and the deformed configuration is given by ϕ and an area element on the deformed midsurface is computed as dΩ = j̄dξ1dξ2 with j̄ = ‖ϕ,1 ×ϕ,2‖. Thus, each point on the body can be described in both the undeformed state by ϕ0 and the deformed state by ϕ. The deformation u is described by these configurations such that: u = ϕ−ϕ0 . (5.14) Furthermore, the mapping from the reference configuration dϕ0 to the actual state dϕ can be used to describe a deformation gradient F : dϕ = F · dϕ0 . (5.15) This deformation gradient includes rigid body movements. It is defined by the base vectors in the reference and the actual configuration, for example by F = gi ⊗ gi0 and F T = gi0 ⊗ gi, and can be used for mapping between undeformed and deformed base vectors. The Jacobian of the mapping is given by J = det(F ) = j j0 where j = g3 · (g1 × g2) . (5.16) Since the gradient includes rigid body motions it cannot be used directly as a measure of strains. Different strain tensors exist but for this project the Green-Lagrange strain tensor E is used. It is defined by the deformation gradient and the identity tensor I: E = 1 2 (F T · F − I) = Eij g i 0 ⊗ g j 0 . (5.17) By inserting the expressions for F and F T into Eq. (5.17) and recalling that the identity tensor is identical to the metric tensor the following is obtained: E = 1 2 (gij − g0ij) g i 0 ⊗ g j 0 . (5.18) The Green-Lagrange strain coefficients Eij are given by the metric coefficients in the deformed and undeformed state: Eij = 1 2 (gij − g0ij) . (5.19) The coefficients refer to the undeformed state contravariant basis gi0 ⊗ g j 0 and not an orthonormal basis. 5.2.4 The Kirchhoff-Love hypothesis With the kinematics of the surface in place Kirchhoff-Love theory is introduced. The theory is based on the "direct approach" meaning that the shell is regarded as a two-dimensional surface and proper three-dimensional behaviour is assumed. Henceforth the derivations are restricted to thin Kirchhoff-Love theory of shells. As mentioned earlier it is analogous to Euler-Bernoulli beam theory as cross-sections are assumed to remain planar 43 during deformation. The result is a linear strain distribution through the thickness. It is also assumed that the cross-sections remain normal to the middle surface during deformation. This means that the director is always normal to the midsurface, hence the definition of an independent director is superfluous. The shell can be completely represented by its middle surface. The implication mechanically is that shear strains are neglected, like in the Euler-Bernoulli beam model. The assumption is reasonable for thin structures and a good guideline is that the shell radius to thickness ration should be around R/h > 20 [11]. Given that both normal strains and transversal shear strains are neglected only in-plane strain coefficients are considered. Eq. (5.18) can thereby be reduced to: E = Eαβ g α 0 ⊗ g β 0 , (5.20) giving the strain coefficients Eαβ : Eαβ = 1 2 (gαβ − g0αβ) . (5.21) Since every point on the shell can be described by the middle surface and normal vector alone the base vectors ai are given as in Eqs. (5.3) and (5.5) by: aα = ϕ,α , (5.22) a3 = ϕ,1 ×ϕ,2 ‖ϕ,1 ×ϕ,2‖ = t . (5.23) Note that aα differs from gα in that aα = ϕ,α(ξ1, ξ2) whereas gα = x,α(ξ1, ξ2, ξ3). Corresponding to Eqs. (5.7) and (5.11) the metric and curvature coefficients defined on the midsurface are given by: aαβ(ξ1, ξ2) = ϕ,α ·ϕ,β , (5.24) καβ(ξ1, ξ2) = ϕ,α · t,β = −ϕ,αβ · t . (5.25) The position vector x is then given by x = ϕ(ξ1, ξ2) + ξ3t(ξ1, ξ2) which gives the base vectors gα = ϕ,α + ξ3t,α and the metric coefficients gαβ : gαβ = aαβ − 2ξ3 καβ + (ξ3)2 t,α · t,β . (5.26) However, for thin shells the quadratic term in the equation above is omitted leaving: Eαβ = 1 2 (aαβ − a0αβ) + ξ3 (καβ − κ0αβ) . (5.27) The strains are represented by the metric and curvature coefficients of the middle surface. It can be seen from the equations that the strains consist of a constant and a linear part. The constant part correspond to the membrane strains εαβ and the linear part to the effect of bending ραβ . Eq. (5.27) can thus be rewritten as: Eαβ = εαβ + ξ3ραβ . (5.28) 44 5.2.5 Small strain kinematics Small displacements are assumed for the shell deformation field. Since the displacement vector of the middle surface of the shell is defined as u = ϕ−ϕ0 the linear membrane and bending strain measures can be derived from the kinematic variables up to first order in u on the form: εαβ = 1 2 (ϕ0,α · u,β +ϕ0,β · u,α) , (5.29) ραβ = −t0 · u,αβ −ϕ0,αβ ·∆t , (5.30) where t0,α = −j