Implementation and Evaluation of a Coupled Thermal-Structural Analysis Module for Laminated Composites in an Open-Source Finite Element Code by DOUGLAS BRIAN BEATON B.Eng., Dalhousie University, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2010 Douglas Brian Beaton, 2010 Abstract The aerospace industry invests heavily in the design and manufacture of composite materials. Aircraft components are produced by placing unprocessed composite materials in an autoclave and applying heat and pressure. The desired part geometry is achieved by forming raw composite materials around a tool, typically made of aluminum or other metal. Throughout the cure cycle, temperature changes cause the part and tool to expand at different rates. This differential expansion, combined with composite material properties that evolve over time, produces residual stresses in the part and leads to geometric instabilities (warpage) upon removal from the tool. Excessive warpage can render a part unusable. Errors of this nature can be quite costly, particularly in the aerospace industry where the tools created can be very large. A strong desire exists to predict the warpage and residual stresses imposed by the curing process and incorporate these stresses in the structural design of a component. To accomplish this goal for complex geometries, special additions to the finite element method are required. Commercial finite element programs provide some flexibility for users to implement custom elements and materials. Though, this flexibility has limits: some material models, such as non-local damage models, cannot be incorporated in proprietary software. This work selects an open-source finite element program and implements the ability to model curing processes of composite materials. The thermal and structural equations are solved in a coupled manner during each time step. This contrasts previous work by the UBC Composites Group, wherein the heat equation is solved over the entire model before the structural equations are considered. Numerous verification models are run to confirm the implementation, along with several example problems. Recommendations are made for further work to improve the process modeling and facilitate a link to subsequent structural models. Ultimately, the code produced represents the first step in seamlessly modeling composite structures during manufacturing processes through to in-service conditions. ii Table of Contents Abstract .............................................................................................................................. ii Table of Contents ............................................................................................................. iii List of Tables ..................................................................................................................... v List of Figures ................................................................................................................... vi Chapter 1 Introduction ................................................................................................. 1 1.1 Background and Motivation ................................................................................... 1 1.2 Scope ....................................................................................................................... 2 1.3 Software Survey ...................................................................................................... 2 Chapter 2 Overview of the OOFEM Package ............................................................ 7 2.1 Object-Oriented Programming................................................................................ 7 2.1.1 Notes on Nomenclature .................................................................................. 10 2.2 Development of OOFEM ...................................................................................... 10 2.3 Program Features .................................................................................................. 11 2.4 Class Hierarchy ..................................................................................................... 12 Chapter 3 Component Implementation .................................................................... 17 3.1 24-Node 3D Solid Element - Mechanical ............................................................. 17 3.2 24-Node 3D Layered Element – Mechanical........................................................ 19 3.2.1 Layup Class .................................................................................................... 20 3.2.2 Layer Scheme Class........................................................................................ 21 3.2.3 Non Uniform Layer Scheme Class ................................................................. 21 3.2.4 Layered Integration Rule Class ...................................................................... 22 3.3 24-Node 3D Solid Element – Thermal ................................................................. 22 3.3.1 3D – 24-Node Finite Element Interpolation Class ......................................... 23 3.4 24-Node 3D Layered Element – Thermal ............................................................. 24 3.5 CCA Material Class .............................................................................................. 25 3.5.1 Orthotropic Linear Heat Material Class ......................................................... 25 3.5.2 CCA Library Interface Class .......................................................................... 25 3.5.3 CCA Material Class – Thermal ...................................................................... 26 3.6 Coupled Thermal – Structural Analysis Module .................................................. 26 iii 3.6.1 Engineering Model – CTSM .......................................................................... 27 3.6.2 Elements – CTSM .......................................................................................... 28 3.6.3 Materials – CTSM .......................................................................................... 28 3.6.4 CCA Material Class – CTSM ......................................................................... 29 Chapter 4 Verification Models and Case Studies..................................................... 31 4.1 Structural Models .................................................................................................. 31 4.1.1 Cantilevered Beam ......................................................................................... 31 4.1.2 Eigenvalue Analysis ....................................................................................... 34 4.1.3 Circular Ring .................................................................................................. 35 4.1.4 Sandwich Beams............................................................................................. 37 4.2 Thermochemical Models ...................................................................................... 41 4.2.1 Logistic Cure Model ....................................................................................... 42 4.2.2 Single-Hold Cure Cycle Models..................................................................... 42 4.3 Coupled Thermal-Structural Models .................................................................... 48 4.3.1 Flat Composite Part on an Aluminum Tool.................................................... 48 Chapter 5 Conclusion and Future Directions........................................................... 66 5.1 Conclusions ........................................................................................................... 66 5.2 Future Directions .................................................................................................. 66 Bibliography .................................................................................................................... 68 Appendix A OOFEM Class Hierarchy ..................................................................... 69 Appendix B Nonstationary Transient Solving Algorithm ...................................... 80 Appendix C Logistic Cure Model .............................................................................. 83 Appendix D Surface Integration Using 3D Jacobian .............................................. 90 Appendix E CCA Models and Input Parameters .................................................... 94 iv List of Tables Table 1-1 Summary of open-source FEA software capabilities ......................................... 3 Table 4-1 Cantilever beam models - rectangular element normalized results for end-tip deflection........................................................................................................................... 32 Table 4-2 Cantilever beam models - parallelogram and trapezoidal element normalized results for end-tip deflection ............................................................................................. 34 Table 4-3 Circular ring models - normalized results for vertical deflection..................... 36 Table 4-4 Results of the 1D 3501 material ramp rate models .......................................... 45 Table 4-5 Results of the 1D 8552 material ramp rate models .......................................... 46 Table 5-1 Future work for seamless modeling in OOFEM .............................................. 67 Table E - 1 Constant material parameters in the CCA library, Material A ...................... 94 Table E - 2 CCA fibre specific heat capacity input parameters, Material A .................... 95 Table E - 3 CCA fibre thermal conductivity input parameters, Material A ..................... 95 Table E - 4 CCA composite lumped coefficient of thermal expansion input parameters, Material A ......................................................................................................................... 96 Table E - 5 CCA cure kinetics input parameters, Material A ........................................... 96 Table E - 6 CCA resin modulus development input parameters, Material A ................... 97 Table E - 7 Constant material parameters in the CCA library, Material B ....................... 98 Table E - 8 CCA resin specific heat input parameters, Material B................................... 98 Table E - 9 CCA resin thermal conductivity input parameters, Material B ..................... 98 Table E - 10 CCA fibre specific heat capacity input parameters, Material B .................. 98 Table E - 11 CCA fibre thermal conductivity input parameters, Material B .................... 99 Table E - 12 CCA composite lumped coefficient of thermal expansion input parameters, Material B ......................................................................................................................... 99 Table E - 13 CCA cure kinetics input parameters, Material B ....................................... 100 Table E - 14 CCA resin modulus development input parameters, Material B ............... 101 Table E - 15 CCA resin Poisson's ratio development input parameters, Material B ...... 101 v List of Figures Figure 2-1 Example of an object-oriented program code. .................................................. 8 Figure 2-2 Classes derived directly from FEMComponent .............................................. 13 Figure 2-3 Inheritance structure for the Element class ..................................................... 14 Figure 2-4 Inheritance structure for the Material class ..................................................... 15 Figure 2-5 Inheritance structure for the IntegrationRule class ......................................... 16 Figure 3-1 Node numbering for 24-node 3D solid element (taken from (Arafath, 2007) with permission)................................................................................................................ 18 Figure 3-2 In-plane and through-thickness integration points for 3D layered elements (reproduced from (Arafath, 2007) with permission)......................................................... 20 Figure 3-3 Addition of LayeredIntegrationRule to the class hierarchy ............................ 22 Figure 3-4 Addition of FEI3d24Node to the class hierarchy............................................ 24 Figure 4-1 Cantilevered beam geometry........................................................................... 31 Figure 4-2 Rectangular element model of the beam with 10 – 24-node elements (side view, not to scale) ............................................................................................................. 32 Figure 4-3 Parallelogram element model with 10 – 24-node elements (side view, not to scale) ................................................................................................................................. 33 Figure 4-4 Trapezoidal element model with 10 – 24-node elements (side view, not to scale) ................................................................................................................................. 33 Figure 4-5 Circular ring under vertical point loads .......................................................... 35 Figure 4-6 Circular ring model with 4 – 20-node elements in the tangential direction (not to scale) ............................................................................................................................. 35 Figure 4-7 Cantilevered sandwich beam - constant thickness .......................................... 37 Figure 4-8 Cantilevered sandwich beam – tapered ........................................................... 37 Figure 4-9 Constant thickness sandwich beam, solid element mesh ................................ 38 Figure 4-10 Constant thickness sandwich beam, layered element mesh .......................... 38 Figure 4-11 Deflection curves for constant thickness sandwich beam models ................ 39 Figure 4-12 Tapered sandwich beam, solid element mesh ............................................... 40 Figure 4-13 Tapered sandwich beam, layered element mesh ........................................... 40 Figure 4-14 Deflection curves for tapered sandwich beam models .................................. 41 vi Figure 4-15 1D thermal model with 3D solid elements and convective boundary conditions, side view......................................................................................................... 42 Figure 4-16 Temperature history at the mid-plane of Material A as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 1 °C/min .............................................. 44 Figure 4-17 Temperature history at the mid-plane of Material A as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 2 °C/min .............................................. 44 Figure 4-18 Temperature history at the mid-plane of Material A as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 3 °C/min .............................................. 45 Figure 4-19 Temperature history at the mid-plane of Material B as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 1 °C/min .............................................. 46 Figure 4-20 Temperature history at the mid-plane of Material B as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 2 °C/min .............................................. 47 Figure 4-21 Temperature history at the mid-plane of Material B as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 3 °C/min .............................................. 47 Figure 4-22 Flat part and tool models - geometry and structural boundary conditions .... 48 Figure 4-23 Flat part and tool models - meshing and structural boundary conditions ..... 49 Figure 4-24 Thin part model with specified temperature thermal boundary conditions – degree-of-cure results........................................................................................................ 51 Figure 4-25 Thin part model with specified temperature thermal boundary conditions ply E11 values .................................................................................................................... 52 Figure 4-26 Thin part model with specified temperature thermal boundary conditions ply G12 values.................................................................................................................... 52 Figure 4-27 Thin part model with specified temperature thermal boundary conditions axial stress through the part at x = 150 mm, t = 100 minutes ........................................... 53 Figure 4-28 Thin part model with specified temperature thermal boundary conditions axial stress through the part at x = 150 mm, t = 320 minutes ........................................... 54 Figure 4-29 Thin part models - temperature histories at x = 150 mm, y = 1 mm ............. 55 Figure 4-30 Thin part models – degree-of-cure history at x = 150 mm, y = 1 mm ......... 56 Figure 4-31 Thin part models - axial stress history at x = 150 mm, y = 1 mm ................. 57 Figure 4-32 Thin part models - axial stress through the part at x = 150 mm, t = 60 minutes ........................................................................................................................................... 57 vii Figure 4-33 Thick part model with specified temperature thermal boundary conditions axial stress through the part at x = 150 mm, t = 100 minutes ........................................... 59 Figure 4-34 Thick part model with specified temperature thermal boundary conditions axial stress through the part at x = 150 mm, t = 320 minutes ........................................... 60 Figure 4-35 Thick part models - temperature history at x = 150 mm, y = 7 mm .............. 61 Figure 4-36 Thick part models - temperature through the part at x = 150 mm, t = 80 minutes .............................................................................................................................. 62 Figure 4-37 Thick part models – degree-of-cure history at x = 150 mm, y = 7 mm......... 63 Figure 4-38 Thick part models - axial stress history at x = 150 mm, y = 7 mm ............... 63 Figure 4-39 Thick part models - axial stress through the part at x = 150 mm, t = 100 minutes .............................................................................................................................. 64 Figure 4-40 Thick part models - axial stress through the part at x = 150 mm, t = 320 minutes .............................................................................................................................. 65 Figure C - 1 Temperature, degree-of-cure, and normalized heat generation for the logistic cure model ......................................................................................................................... 85 Figure C - 2 Crank-Nicolson stencil ................................................................................. 86 Figure C - 3 Finite difference geometry stencil for time step n ........................................ 87 Figure C - 4 Temperature histories for logistic cure model, Excel™ finite difference model................................................................................................................................. 89 viii Chapter 1 1.1 Introduction Background and Motivation The UBC Composites Group is an interdisciplinary research group comprised of both Civil and Materials engineers. A main focus of the group’s research is numerical modeling of the physical behaviour of laminate composite structures. Two aspects of primary concern are process modeling and damage modeling. Over the years the group has developed advanced computational methods in both of these areas. These methods have been implemented either through additions to commercial finite element (FE) programs, such as LS-DYNA or ABAQUS, or via in-house computer programs, such as COMPRO (process modeling) and CODAM (damage mechanics constitutive models). The structure of commercial FE codes limits the number of additional modules that users can implement at any given time. Due to this limitation, a user could not run a model wherein many custom program components are used at once. Furthermore, since the source code of commercial packages is not available to licensees, would-be developers cannot see or alter the inner workings of a program. This lack of transparency means that some classes of numerical methods cannot be added to commercial codes by end users. For example, non-local damage models cannot be implemented in ABAQUS by anyone outside of SIMULIA (the company which maintains and develops ABAQUS). It has long been a goal of the Composites Group to seamlessly model composite structures throughout their entire life cycle. That is, beginning with process modeling, moving on to in-service performance, and concluding with failure or decommissioning. Such a model should take into account all induced warpage and residual stresses during manufacture and include any effects of damage during in-service conditions. To make this process seamless the use of multiple programs should be avoided, as should any computationally expensive procedures such as remeshing. Any platform capable of this seamless modeling would be a valuable tool in, for example, the aerospace industry: where material and testing costs are large, and difficulty or uncertainty in modeling the material behaviour necessitates high safety factors. 1 To realize this goal it is necessary to look beyond the realm of commercial FE software. Two options are present in the non-commercial realm: develop a new FE package including all desired capabilities, or add new functionality to an existing open-source code. Given the amount of work required to develop a new code it is much more efficient to utilize an existing package. 1.2 Scope This project establishes the first steps in seamlessly merging process modeling with structural analysis of composite materials. First a survey of available open-source FE packages is performed. Once a program is selected, necessary components are added, including solid-shell elements and an interface to a material processing library. Steps are then taken to verify this implementation. Once the implementation is verified a case study is performed to assess the effects of manufacturing parameters on a composite part. Some closing remarks and discussion of future directions conclude the thesis. 1.3 Software Survey A large number of open-source software packages are currently available for finite element modeling. These packages can vary in platform, programming language, capabilities, documentation, support and expandability, among other things. A thorough survey of available software was performed in selecting a package for this research. Potential programs were initially screened based on capability. Those found to have the required basic features were short-listed for further investigation. The initial search revealed twelve software packages which met the basic requirements. A summary of the capabilities of these packages is shown in Table 1-1 below. 2 Heat Tran sfer Fluid Mec hanic s Dyna mic A nalys is Linux Windows, Linux Fortran FENICS Linux C++ / Python NLFET Matlab Matlab Linux C++ Code_Aster FEAP OOFEM OpenFEM PSU TAHOE TOCHNOG Matlab Matlab Linux Windows, Linux C++ / Fortran Linux C++ Warp3D Windows, Linux Fortran Z88 Windows, Linux C 1 ion Nonli near Mate rials University of Tokyo EDF (Electricité de France) University of California, Berkeley University of Chicago (primary) Texas A&M University Czech Technical University in Prague INRIA / SDTools Sandia National Labs University of Twente University of Illinois at Urbana‐Champaign Deve lopm ent L ocat Nonl inear Ge Mesh visua liz atio n 3D El emen ts Windows, Linux Fract ure M echa nics 2D Ele ment s Writt en in ADVENTURE Platfo rm Python / Fortran Nam e 1D Ele ment s omet ry Table 1-1 Summary of open-source FEA software capabilities1 University of Bayreuth Empty cells indicate features that were uncertain or not specified. 3 Each program in the short list was downloaded and installation attempted on a Windows PC, Linux cluster, or a virtual machine running Linux. If installation was successful, a simple model was then run to explore the program’s user interface and verify basic functionality. Once these steps were taken the more advanced features of each program were also investigated. After considering all programs on the list, four finalists were chosen. Code_Aster: Aster is a very powerful open-source program for multiphysics modeling (Électricité de France, 2007). The program has an expanding community of users and many advanced features. Aster was developed by Électricité de France (EDF), the main electric utility in France. All source code for the program was written with comments and error messages in French. The majority of documentation for Aster is also written in French, though a portion is available in other languages. Although the difference in language may cause sporadic problems for English-speaking users, Aster continues to increase in popularity due to a wide range of capabilities. Unfortunately though, a source code written in French, and little English documentation, make Aster a poor candidate for expansion by non-French-speaking developers. FEniCS: The FEniCS project is a joint undertaking between several institutions, spearheaded by the University of Chicago. The vision of this project is to set a new standard in the automation of computational mathematical modeling, towards the goals of generality, efficiency, and simplicity, concerning mathematical methodology, implementation, and application (Logg, 2009). The project takes an innovative approach in that it automates the process of software generation. Typical finite element programs accept physical data as input and proceed to generate a mathematical solution to the problem. FEniCS accepts physical data and a governing differential equation as input and proceeds to generate an executable program to solve the problem. Thus each problem entered into FEniCS produces an executable file. This file, when run, will output a solution to the physical 4 problem. This approach works toward the goals of the FEniCS project and has potential for very fast run times. Although FEniCS may represent a new direction in the field of finite element modeling, the project is still in its infancy. This early stage of development, coupled with the sheer complexity of the program, make it a poor candidate for further expansion. A disappointing lack of documentation also makes FEniCS inaccessible to potential users and developers. TOCHNOG: TOCHNOG was originally developed at the University of Twente in Holland, and has since been adopted as a commercial code by FEAT. However, the open-source version remains available and thus TOCHNOG was included in the short-list of candidate packages. This program has a wide range of capabilities (Roddeman, 2001), though its primary focus is geotechnical modeling. Unfortunately, active development of TOCHNOG was halted in late 2001, corresponding to the program’s adoption as a commercial finite element code. This lack of active development, combined with insufficient programming documentation and a less favourable program structure, were sufficient reason to eliminate TOCHNOG as a candidate for extension. OOFEM: OOFEM is an open-source finite element code with object-oriented architecture (Patzák, 2008). It has been under development at the Czech Technical University in Prague since 1997. This package offers many analysis types and has an extensive material library. Although other packages are available with equal or greater modeling options, OOFEM was designed with extensibility in mind. Its object-oriented architecture makes the addition of new elements or materials relatively simple. OOFEM also has a sizeable user and developer community. Perhaps the most important feature of OOFEM is the level of support available for modellers and programmers. Several manuals are available covering modeling procedures and some program structure. The program structure is further outlined by extensive online documentation. Finally, a forum is available on the OOFEM 5 website where developers can post questions. Most often any questions are addressed in a timely manner by OOFEM’s original author. A large catalogue of features, combined with amicable program structure and ample support, made OOFEM the best choice for further development in this research. The code itself is explained in more detail in the following chapter. 6 Chapter 2 2.1 Overview of the OOFEM Package Object-Oriented Programming Object-oriented programming (OOP) differs considerably from procedural programming in several ways, the most obvious being code structure and data abstraction. When first learning a programming language it is natural for students to begin with procedural programming. Here programs are broken down into various functions (subroutines, etc.). There is a single flow control which calls the various functions as necessary. Any procedures that are repeated often can be separated into functions and there is some potential for code reuse. In OOP programs are better viewed as a system of objects with various interactions. These objects can all have various attributes, as well as corresponding behaviour. An example in OOFEM is an element. In OOFEM each element will be a separate object. Each element will have (among other things) an integer called ‘numberOfDofMans’, which gives the number of nodes the element has, as well as another integer called ‘material’, which denotes which material in the model is associated with that particular element. These are some of the attributes associated with the ‘Element’ object class. Object classes will also have functions to define their behaviour. One of the functions of the ‘StructuralElement’ class (a subclass of the Element class) is ‘computeStiffnessMatrix’ which is self-explanatory. These concepts are more easily illustrated by example. Figure 2-1 shows a sample of code written in C++. In C++, and in many languages, each type of object is called a class. Lines 7 - 25 in the code define a class called ‘Person’. This class has three attributes: an integer called ‘age’, a float-point number called ‘height’, and a string called ‘name’. It also has a number of functions including ‘giveAge’. In OOP, each time an object of a given class is made it is called an instance of that class. Calling the function ‘giveAge’ for a given instance of the class ‘Person’ would return the value of the integer ‘age’ for that particular instance. 7 Figure 2-1 Example of an object-oriented program code. Note that in the code in Figure 2-1 the class ‘Person’ may represent a real person, but it does not include all possible attributes of a person. Many attributes normally associated with people are missing from this class (eg. gender, eye colour, birthplace, etc.). It may be, however, that for the purposes of this program it is not necessary to know a person’s eye colour. By only taking into account the necessary attributes and functions of a person this code generates an abstraction of real persons. There are three main principles which define OOP. These are: encapsulation, inheritance, and polymorphism. Encapsulation: In OOP, the representation of an abstraction of a physical object as a self-contained set of internal variables and internal functions is termed encapsulation (Yevick, 2005). The fact that these sets of variables and functions are self-contained means that some information in an object is not accessible by other objects. Specifically, any variables or functions defined in the public section of a class are accessible to other objects. Any variables or functions defined in the private section of a class are only accessible by that specific object. Returning to the example in Figure 2-1, note that the only way for an 8 outside object to determine the integer ‘age’ for an object of the Person class, is to call the public member function ‘giveAge’. Similarly, there is no way for an outside object to access the string ‘name’, as there is no public member function which returns this string. Variables can be placed in the public section, making them accessible to any other object without the need to call a function. Functions can also be placed in the private section, meaning that these functions can only be called by other functions of the same object, not by outside objects. Inheritance: One of the most useful aspects of OOP is the code reuse enabled by inheritance. Through inheritance, new subclasses can be defined which acquire the variables and functions of an existing class2. Additional attributes can be defined for this new derived class to make it differ from the original base class. Functions in the base class can also be redefined in the derived class3. Finally, a convenient feature of inheritance is that a derived class can be substituted for a base class in any situation in a program. For example, a function that requires an object of the base class as an argument can be executed with an object of the derived class as an argument. Polymorphism: Polymorphism is the assignment of multiple context-specific behaviours to a single program entity (Yevick, 2005). This means that two different classes can have a function of the same name, but with different purposes. For example, in OOFEM the function ‘instanciateYourself’ appears in many different classes, and has slightly different behaviour in each. So the performance of a function of a given name will vary depending on the type of object calling it. The main benefit of using OOP for an open-source package such as OOFEM is that, if done properly, it will produce a code that is well structured and easily extensible. 2 Functions or variables that are public in the base class may not necessarily be public in the derived class. For a detailed explanation of public inheritance vs. private inheritance see (Yevick, 2005). 3 Depending on the situation, the original function from the base class my still be called from an object of the derived class. See (Yevick, 2005) for an explanation of virtual functions and virtual inheritance. 9 Software is also available which partially automates the process of documenting objectoriented programs, thus saving considerable time for the programmer. 2.1.1 Notes on Nomenclature Specific object classes are frequently discussed throughout this document. To maintain flow and readability, some liberties have been taken with regard to class nomenclature. The first time a class name appears in the text, it will be surrounded by single quotes, (ie. ‘CCAMaterial’). After the first appearance, quotes will be omitted when referring to a class by name. Secondly, terms such as ‘CCAMaterial class’ refer to the abstract class ‘CCAMaterial’, while terms such as ‘CCAMaterial object’ and ‘instance of the CCAMaterial class’ refer to the actual manifestations of the abstract class produced at runtime. To clarify, within OOFEM there will be only one ‘Brick24’ structural element class. But while running a model using these elements, there will be many Brick24 objects created. 2.2 Development of OOFEM According to the OOFEM website (Patzák, 2008), the program has been in constant development since 1997. The main architect of OOFEM is Dr. Bořek Patzák, an Associate Professor at the Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics. The project grew out of Dr. Patzák’s Ph.D. thesis and has since been adopted by various researchers and modellers for different purposes. Development of OOFEM was partially funded by: the Ministry of Education of the Czech Republic, Grant Agency of the Czech Republic, and CEC Brussels (Patzák, 2008). A number of unrelated developers have also made contributions in the form of code, documentation or by reporting bugs. At the time of this writing, the current release of OOFEM is version 1.8. This current version of OOFEM is available under the terms of the GNU General Public License version 2, see (Free Software Foundation, 1991) for a copy of this license. 10 2.3 Program Features OOFEM version 1.8 has many features desirable for finite element modeling. The following is a partial list of the features outlined on the OOFEM website (Patzák, 2008): Modular and Extensible FEM Kernel (OOFEMlib): Fully extensible – The kernel is extensible in any ‘direction’. New elements, materials, boundary conditions, numerical algorithms or analysis modules can all be added with modest effort. Parallel processing support – Many analyses can be run in parallel and good performance scalability can be obtained on various platforms. Efficient sparse solvers – Direct and iterative solvers are available. Interfaces to some third party linear solver libraries are also available. Adaptive analysis support – Multiple domain concept, error estimation and remeshing, primary unknown and internal variables mapping, and fast spatial localization algorithms are available. Structural Mechanics Module (sm): Many analysis procedures - linear static, linear dynamic (eigenvalue analysis, implicit and explicit direct integration), nonlinear static, nonlinear dynamic (explicit, parallel version). Large material library – includes damage models and nonlinear fracture mechanics Advanced modeling features – slave degrees-of-freedom, rigid arms, local coordinate systems Transport Module (tm): Analysis procedures – stationary and transient (linear and nonlinear) heat transfer and coupled heat and mass transfer problems Element library – axisymmetric, two and three dimensional elements 11 Staggered simultaneous solution – of heat transfer analysis and mechanical analysis Fluid Dynamics Module (fm): Analysis procedures – transient incompressible flow Element library – linear equal order velocity and pressure approximation triangle Postprocessing: Built-in X-windows postprocessor Export to VTK format supported 2.4 Class Hierarchy The overall class hierarchy of OOFEM is very large and will not be shown here in its entirety. Some key class structures of primary concern to this project will be outlined below. A full, textual version of OOFEM’s hierarchy is provided in Appendix A. The largest single class structure in OOFEM is the ‘FEMComponent’ class and its derivatives. The majority of classes which would concern a developer are included in this structure. Several important high-level classes such as ‘Element’, ‘Material’, and ‘IntegrationRule’ are directly derived from the FEMComponent base class. Figure 2-2 gives a graphical sample of this class structure. Note that each of the subclasses shown in Figure 2-2 will have one or many of its own subclasses, which in turn will have more subclasses, and so on. This diagram shows only the top two levels of the FEMComponent hierarchy. 12 Figure 2-2 Classes derived directly from FEMComponent Some of the subclasses shown in Figure 2-2 will be important to this project. Thus it is worth taking a closer look at these classes. Figures 2-3 and 2-4 depict the large element and material libraries provided with OOFEM release 1.8. Also, as each class shown requires both header and source files, the figures illustrate how much code is required in a program like OOFEM. Fortunately, OOFEM’s object-oriented structure allows developers to make additions to the program while ignoring much of the existing code. 13 Figure 2-3 Inheritance structure for the Element class 14 Figure 2-4 Inheritance structure for the Material class 15 Attention should also be drawn to the IntegrationRule class, which will be extended during this project. This particular class has only two derived classes, as shown in Figure 2-5. Figure 2-5 Inheritance structure for the IntegrationRule class 16 Chapter 3 Component Implementation This chapter describes the implementation of a number of components required for modeling laminated composite materials. OOFEM’s object-oriented structure makes adding new capabilities straightforward, with minimal change to existing code. For example, the general procedure for adding a new element to OOFEM is as follows: 1 - Create header and source files for the new element 2 - Add the new class name to the enum ‘ClassTypes’ in cltypes.h 3 - Add the ‘#include’ line for your new header file to usrdefsub.C 4 - Add the option for your new element to the function ‘CreateUsrDefElementOfType’ in usrdefsub.C 5 - Add your source file to the appropriate makefile and makefile.in (for future configuration) 6 - Make any other required additions or changes 7 - Clean existing object files (only necessary for major changes/additions) 8 - Recompile the program Steps 2-5 are quite simple and require only six new lines of code. Steps 7 and 8 rebuild the program with the new additions. Much of step 1 can be accomplished by copying and modifying the code for an existing element. The remaining work in step 1 and, if necessary, step 6 constitute the majority of work required in adding a new element. A comparison of the elements outlined in this section will reveal the impact that step 6 can have. 3.1 24-Node 3D Solid Element - Mechanical Analysis of large plate-like structures is commonly done using degenerated shell elements. Nodes within these 2D elements are coplanar and have translational and rotational degrees-of-freedom, thus making them compatible with 1D beam-column elements, but incompatible with 3D solid elements. Formulation of these shell elements is typically based on constant through-thickness properties, and ignores out-of-plane normal stresses. Laminated composite components subjected to a cure cycle will exhibit 17 considerable variation of properties in the through-thickness direction. Stacking orthotropic layers of prepreg at varying fibre angles is one source of this variation. Thick composite parts may also show non-uniform modulus development in the throughthickness direction. To capture this behaviour it is preferable to use 3D solid elements. The use of solid elements also provides many other advantages over shell elements, particularly in nonlinear and contact problems (Arafath, 2007). When modeling plate-like structures using 3D solid elements, it is desirable to use large element aspect ratios to reduce the number of elements required. 3D solid elements formulated to use high aspect ratios in modeling plate-like structures are termed solidshell elements. In formulating these elements, care must be taken to avoid shear-locking and any other phenomena which could reduce accuracy. To accomplish this, the elements implemented in OOFEM have 24 nodes as shown in Figure 3-1. This configuration permits higher aspect ratios than lower-order elements, thereby reducing the number of elements required to model plate-like structures and increasing computational efficiency (Arafath, 2007). Thickness direction 15 12 11 23 19 20 16 7 24 8 9 3 14 4 22 6 13 10 21 17 18 z 5 1 2 y x Figure 3-1 Node numbering for 24-node 3D solid element (taken from (Arafath, 2007) with permission) 18 The first addition to OOFEM was a 3D solid element with 24 nodes. Reduced-integration and full-integration are achieved using 8 or 27 integration points (IPs), respectfully. IPs are arranged in a 2x2x2 or 3x3x3 cube in natural coordinates as is typical for 3D brickshaped elements. The element was implemented in OOFEM’s mechanical (structural) module. OOFEM version 1.8 includes a 20-node solid element, implemented as the class ‘QSpace’. The QSpace class source code was copied and modified to create the new 24node element class, ‘Brick24’. The QSpace element included in OOFEM is similar to that shown in Figure 3-1 using only nodes 1 through 20. As this new element class was implemented by copying and modifying the QSpace class, it resides at the same level in the class hierarchy. 3.2 24-Node 3D Layered Element – Mechanical The Brick3D class outlined above models elements with uniform material properties throughout. To represent multiple plies with varying fibre angle in a single element, a layered integration scheme was developed in (Arafath, 2007). This approach divides a single element into different regions (layers) with constant material properties. Typically, each region would represent a single ply or group of plies with constant fibre angle. Each layer is given a specific number of IPs in the through-thickness direction while each element uses a constant number of IPs in the in-plane directions (see Figure 3-2). Such an approach captures changes in material properties due to different fibre angles as well as varying cure rates throughout the element. In the work presented, solid elements using this type of integration scheme are termed ‘layered’ elements. 19 7 4 3 3 tn Layer n 4 8 6 t2 Layer 2 1 2 t1 Layer 1 1 5 2 Figure 3-2 In-plane and through-thickness integration points for 3D layered elements (reproduced from (Arafath, 2007) with permission) The next addition to OOFEM was a 24-node 3D layered element designed for use with composite laminates. Again this implementation was part of the structural module. This process was more involved than adding the solid element. The first step was to determine how much information would be contained within the new element class and how much should be placed elsewhere in the OOFEM object hierarchy. As many layered elements in a model will have the same layup and layer scheme, it is not necessary to repeat these details in each element. In fact, doing so would create a great deal of redundant information in the model input file, and waste significant program memory. The decision was made to create two new classes to store this information: ‘Layup’, and ‘LayerScheme’. Each layered element would then have a number indicating its associated layup and layer scheme; analogous to how materials are associated with elements in OOFEM. This setup saves on memory allocated at run time and reduces the size of input files. 3.2.1 Layup Class The Layup class stores information about the layup of a laminated composite. Specifically, it contains an array of floating point numbers representing the fibre angle of each ply. Similar to the notation commonly used with laminates, shortcuts have been provided to avoid repeating entries when defining symmetric and anti-symmetric layups. 20 The Layup class was implemented as a base class and it does not inherit features from any other part of OOFEM. This is a very simple class and it is unlikely that any derived classes will need to be produced. Substantial modifications were made to the domain.C and domain.h to facilitate defining layup objects in the model input file. 3.2.2 Layer Scheme Class The LayerScheme class stores information about the thickness and number of throughthickness integration points in each layer. This class also provides services to the ‘LayeredIntegrationRule’ class for setting up integration points on a layered element. LayerScheme was implemented as a base class and does not inherit features from any other class in OOFEM. This base class only supports layer schemes with layers of uniform thickness. Layer schemes using nonuniform thicknesses are handled through a derived class (see Section 3.2.3 below). As with the Layup class, substantial modifications were made to domain.C and domain.h to facilitate defining layer scheme objects in the model input file. 3.2.3 Non Uniform Layer Scheme Class In many situations the thickness of certain layers in a composite laminate are not constant (tapered sandwich beams, ply drop offs, etc.).To model these structures using 3D layered elements it is necessary to vary the thickness of layers within the element. In OOFEM this is accomplished using ‘NULayerScheme’: a subclass of LayerScheme with nonuniform layer thicknesses. In this class layer thicknesses are defined at the four corners of the element and any thicknesses within are interpolated from these values. As this class is derived from LayerScheme, it provides the same services to LayeredIntegrationRule for setting up integration points. Thus to switch from constant layer thickness to variable layer thickness a modeller need only change the type of layer scheme used; the element declaration remains the same. This illustrates one of the benefits of inheritance in OOP. 21 3.2.4 Layered Integration Rule Class The LayeredIntegrationRule class provides services to layered elements that would be provided by ‘GaussIntegrationRule’ for solid elements. Its primary role is to set up the Gauss points for a given element, determining their coordinates and weighting coefficients. This procedure is not straightforward for layered elements, which have many more through-thickness integration points than typical solid elements. Determining the natural coordinates of these through-thickness integration points is particularly intricate. Consideration must be given to the position of a Gauss point within its own layer, as well as the position and relative size of the layer in the layup. As a substitute for GaussIntegrationRule the LayeredIntegrationRule class was implemented at the same level in the hierarchy. Figure 3-3 highlights the position of this new class. Figure 3-3 Addition of LayeredIntegrationRule to the class hierarchy 3.3 24-Node 3D Solid Element – Thermal Once the new structural elements were complete similar elements were added to the thermal4 module. The 24-node 3D solid element was implemented in the ‘TM_Brick24’ class. This new element shares the same nodes, integration scheme and shape functions as the corresponding structural element. Some adjustments were necessary to switch from 4 To avoid repetition, the terms ‘thermal’ and ‘transport’ are used interchangeably throughout this document. 22 the three displacement degrees-of-freedom in structural analyses to the single temperature degree-of-freedom in thermal analyses. Additional functionality was also required for computing convection flux vectors and other surface integrals. To match existing OOFEM code for computing shape functions and surface integrals, a new interpolation class was implemented for all 3D solid elements with the same 24-node formulation. 3.3.1 3D – 24-Node Finite Element Interpolation Class OOFEM interpolation classes provide services for geometric computations. Their main roles are evaluating shape functions, transforming global and local coordinates, computing surface integrals and determining the Jacobian. Each interpolation class is implemented for elements with a given geometry (node numbering) and set of shape functions. For example, the existing ‘FEI3dHexaLin’ class is appropriate for 8-node 3D brick elements with standard linear shape functions. To accompany the 24-node thermal elements, a new ‘FEI3d24Node’ interpolation class was implemented. Figure 3-4 below shows the positions of this new class in the hierarchy. This class was a straightforward extension of the interpolation class for brick elements, with a few exceptions. Considerable changes were necessary to go from surfaces with four nodes to surfaces with eight or nine nodes, as well as edges with two nodes to edges with three nodes. In some functions it was necessary to add an argument specifying which surface was under consideration, as some surfaces in the 24-node element have eight nodes while others have nine. 23 Figure 3-4 Addition of FEI3d24Node to the class hierarchy In contrast to OOFEM’s interface classes, where functionality is made available using multiple inheritance, interpolation classes will instead be an attribute of a given element class. Each element will instantiate a corresponding interpolation object, and call on this object to perform geometric calculations. 3.4 24-Node 3D Layered Element – Thermal Following the addition of a 24-node 3D solid element to the thermal module, a 3D layered version was also implemented. This layered element differs from previous formulations of solid-shell elements developed for processing modeling in the Composites Group. Where previous thermal elements used Fourier series in their shape functions with additional element degrees-of-freedom as terms in these Fourier series (Rasekh, 2007), this element uses the same shape functions as the 24-node structural elements and stores all degrees-of-freedom at the nodes. As such, this layered element also used the interpolation class outlined in Section 3.3.1. Matching the shape functions and integration scheme of this element to those of the layered structural element saved considerable programming effort when implementing the coupled thermal-structural layered element outlined in Section 3.6.2. 24 The 24-node 3D layered element in the thermal module was implemented as the ‘TM_Layered24’ class. With supporting code already in place from the structural implementation, copying and modifying the new solid thermal element into a layered element was a routine procedure. Significant portions of code were copied from the layered structural element, with only minor revisions necessary. The new thermal layered element has the same integration scheme as its structural counterpart, and employs associated Layup and LayerScheme objects in the same manner. 3.5 CCA Material Class The cure kinetics equations used in process modeling of composite materials are evaluated in OOFEM by an external library. This library was originally developed in (Arafath, 2007) as a module of the COMPRO Component Architecture (CCA), and is hereafter referred to as the CCA library. The version of the CCA library used here was compiled and provided by Convergent Manufacturing Technologies (Convergent Manufacturing Technologies Inc., 2009) and linked to OOFEM via an interface class. 3.5.1 Orthotropic Linear Heat Material Class As a prerequisite for the CCA material class, it was necessary to implement an orthotropic linear material for heat transfer analyses. This was a fairly straightforward extension of the isotropic heat transfer material. The only major changes were defining three values of material conductivity (in each principal direction) and adding the option of defining a local (material) coordinate system. These changes are implemented in the ‘OrthotropicHeatTransferMaterial’ class. 3.5.2 CCA Library Interface Class OOFEM interface classes contain code that is common to many other classes which do not inherit functionality directly from one another. The desired functions are instead defined in an interface class and passed on to subclasses via multiple inheritance. This 25 reduces the amount of redundant code in the program and facilitates faster revisions if a common function should need modification. Note also that OOFEM interface classes are only used to provide functionality to other object classes, and are never instantiated themselves. The ‘CCALibraryInterface’ class was implemented as a subclass of the more general ‘Interface’ class. This new object class provides the external function ‘compro3d’ which is used to call the CCA library. Any material class that uses the CCA library would therefore inherit the compro3d function from this interface class. 3.5.3 CCA Material Class – Thermal The CCA material class implemented in the Thermal module provides the functionality necessary to evaluate equations of cure kinetics for a composite material subjected to an autoclave cure cycle. This functionality is made available via a pre-compiled Fortran library (the CCA library). The OOFEM class ‘CCAMaterial’ provides the framework for calling this Fortran library. The Fortran library is called during every iteration to determine heat generation and thermal properties in the material. These properties are stored by instances of the ‘CCAMaterialStatus’ class at every Gauss point in an element. The CCAMaterialStatus class also tracks values of temperature and degree-of-cure for the current and previous time step, as these are the required input values for calling the CCA Fortran library. Using the properties stored in the material status classes, the CCAMaterial provides all the necessary functions of a typical transport material (ie. returning the material conductivity matrix, specific heat capacity, etc.). 3.6 Coupled Thermal – Structural Analysis Module To facilitate the modeling of coupled thermal-structural problems in OOFEM a new module was created. This module, named the Coupled Thermal-Structural Module (CTSM), was implemented on the same level as the structural, transport, and fluid mechanics modules (SM, TM and FM respectively). The decision to implement the 26 coupled thermal-structural solver in a separate module (rather than in either of the structural or transport modules) was based on the amount of functionality required. All material and elements classes in this new module require both thermal and structural functionality. To add this much functionality to either of the existing modules would result in complex and convoluted code. This would make further development by future programmers rather difficult. So the decision was made to implement coupled analysis capabilities in a separate module. With the desired functionality for the coupled analysis module already present in two other modules, an attempt was made to create the new CTSM classes using multiple inheritance from the structural and thermal modules. This approach was not successful. Since elements and materials in the SM and TM classes share common parent classes, ambiguities are created when attempting to use multiple inheritance. One optional workaround is to use virtual inheritance to eliminate these ambiguities. This approach was also attempted, but was found to produce scope-resolution issues in other areas of the code. Ultimately the decision was made to copy and paste appropriate pieces of code from the two original modules into the new coupled module, and make any necessary modifications. 3.6.1 Engineering Model – CTSM The first component required in this new module was an engineering model class. To facilitate modeling a composite part curing in an autoclave, the engineering model considered consists of a nonlinear transient thermal analysis coupled with a linear static structural analysis. The thermal portion of the analysis is similar to that in the TM module, while the structural portion varies considerably from the SM module. Although the structural analysis is linear within a given time step, the composite material will exhibit nonlinear behaviour over time. This nonlinear response is a result of mechanical material properties evolving as the material cures. Such behaviour means the resulting stresses and deformations at any time in the model will be path-dependent. Thus 27 when solving the structural problem at any given time step, incremental results must be obtained and added to results from the previous time step. This incremental tracking of results is similar to the methods used in nonlinear structural models in OOFEM. However, unlike the nonlinear analyses in the SM module, for the engineering model implemented here, there is no need to iterate when solving the structural problem. 3.6.2 Elements – CTSM At present, only 3D solid elements have been implemented in the CTSM module. Included are: 8-node brick element - ‘CTSMBrick8’ 24-node solid (brick) element -‘CTSM_Brick24’ 24-node layered element -‘CTSM_Layered24’ This module is not intended to function with 1D truss and beam elements or 2D plates and shells. Thus there is no requirement for an element to have an associated cross section. To simplify code readability and tracking program flow, the CrossSection class was not implemented in this module. Functions requiring the use of a CrossSection object in the SM module were revised to eliminate this requirement in the CTSM module. When combining the functionality from an element in the SM module with its counterpart in the TM module, minor modifications were required to ensure the correct version of similar functions is called at appropriate times (ie. returning the B-matrix for a thermal element during the thermal portion of an analysis, rather than the B-matrix for a structural element). 3.6.3 Materials – CTSM At present, there are three material models implemented in the CTSM module: Linear isotropic (all properties) - ‘CTSMIsotropicLinearElasticMaterial’ Linear orthotropic (all properties) - ‘CTSMOrthotropicLinearElasticMaterial’ 28 CCA material - ‘CTSMCCAMaterial’ The linear isotropic material is a straightforward combination of a linear structural material and linear transport material. The linear orthotropic material is also a straightforward extension of its structural and thermal counterparts. Similar to elements implemented in the CTSM module, the only substantial modifications required in the two linear materials were ensuring the appropriate thermal/structural version of a function was called during the appropriate part of the analysis. 3.6.4 CCA Material Class – CTSM The final class implemented in OOFEM was the CTSM module version of the CCA material. This class expanded on the version implemented in the TM module, adding support for the necessary structural material functions. The addition of structural material functions was relatively straightforward, with the typical consideration given to functions specific to one portion of the analysis (thermal vs. structural). One significant programming decision was required for this material. The CCA library returns only those material properties requested during each library call. In the TM module, only the relevant thermal properties were requested when calling the library. The CTSM module will require both thermal and structural material properties during the appropriate part of the analysis. Since material properties for a given time step are determined by the thermal analysis, structural properties could be calculated in one of two ways: 1) Return thermal and structural material properties during each call to the CCA library 2) Return only thermal properties during the thermal portion of the analysis. Make an additional call to the CCA library to calculate structural properties after the temperature results have converged. Option 1 makes each call to the library less efficient, which could slow things down considerably since the library may be called several times when solving the thermal problem. Option 2 allows each library call to be as efficient as possible, although each 29 time step will now require an additional call. Given the large number of properties required to solve the structural problem, it is likely that returning these properties during every iteration of the thermal problem would be less efficient than making one additional library call per time step. Thus the code implemented in OOFEM has been based on option 2. 30 Chapter 4 Verification Models and Case Studies This chapter outlines the verification and case study models investigated in this work. A number of structural, thermal and coupled models were run to test the code implemented in OOFEM. Input parameters and results from these models are provided in the following sections. 4.1 Structural Models 4.1.1 Cantilevered Beam The first structural models investigated were simple cantilever beams subject to either a uniformly distributed load or a point load at the beam’s end. Figure 4-1 shows the geometry and loads used. These models were designed to test the geometric limitations of the 24-node solid and layered elements while comparing their performance to the 20node solid element in OOFEM. Three types of models were considered, as outlined in the following sections. Each model was run using full and reduced integration. Figure 4-1 Cantilevered beam geometry Rectangular Elements The first set of models used rectangular elements and tested the maximum aspect ratio the elements can use while maintaining accurate results. Figure 4-2 shows an example of the rectangular element models. A single element was used in the width and thickness directions. The number of elements in the longitudinal direction was varied to test 31 different element aspect ratios. Models were run using full and reduced integration. Layered element models used a single layer with 5 through-thickness integration points. Table 4-1 gives the end-tip deflections from these models, normalized with respect to the theoretical result. The following parameters were used as input: Length, L = 1.0 m Width, W = 1.0 m Thickness, t = 10 mm Distributed Load, q = 1.0 kPa End Shear Load, V = 100 kN Modulus of Elasticity, E = 200 GPa Poisson’s Ratio, ν = 0.0 Theoretical (beam theory) deflection under distributed load, δUDL = 7.5 mm Theoretical (beam theory) deflection under end shear load, δV = 2.0 mm Figure 4-2 Rectangular element model of the beam with 10 – 24-node elements (side view, not to scale) Table 4-1 Cantilever beam models - rectangular element normalized results for end-tip deflection Element Type Integration 5 Scheme End Shear Load, V Distributed Load, q Number of Longitudinal Elements Number of Longitudinal Elements 1 2 4 8 10 1 2 4 8 10 20‐node Solid 3x3x3 0.750 0.938 0.985 0.996 0.998 n/s n/s n/s n/s n/s 24‐node Solid 3x3x3 9 in‐plane, 5 through‐ thickness 0.750 0.938 0.985 0.996 0.998 0.667 0.917 0.980 0.995 0.997 0.750 0.938 0.985 0.996 0.998 0.667 0.917 0.980 0.995 0.997 2x2x2 1.000 1.000 1.000 1.000 1.000 n/s n/s n/s n/s n/s 24‐node Layered 20‐node Solid 24‐node Solid 24‐node Layered 2x2x2 4 in‐plane, 5 through‐ thickness Unstable Results 1.000 1.000 1.000 1.000 Unstable Results 1.000 1.000 1.000 1.000 1.000 1.000 5 For 20 and 24-node solid elements, 3x3x3 yields 27 IPs (full integration), while 2x2x2 yields 8 IPs (reduced integration). For 24-node layered elements, the in-plane directions are fully integrated with 9 IPs and partially integrated with 4 IPs. 32 As shown in Table 4-1, models with four or more elements in the longitudinal direction returned results within 2 % of the theoretical value. This number of elements in the longitudinal direction corresponds to an aspect ratio of 25. Therefore, the 24-node elements, when used with rectangular geometry, will return accurate results for aspect ratios of 25 or less, regardless of the integration scheme used. At the time of this comparison the 20-node solid element in OOFEM did not support surface loads. Note that the 24-node solid elements returned unstable results when using a 2x2x2 integration scheme. This instability is addressed in Section 4.1.2. Parallel and Trapezoidal Elements The second two sets of models test the maximum angles an element side can have while still producing sound results. These tests follow those performed in (Arafath, 2007) to evaluate element performance. Element sides perpendicular to the longitudinal direction were sloped to create elements of either parallelogram or trapezoidal shape. Figures 4-3 and 4-4 show examples of these models. Figure 4-3 Parallelogram element model with 10 – 24-node elements (side view, not to scale) Figure 4-4 Trapezoidal element model with 10 – 24-node elements (side view, not to scale) Each parallelogram and trapezoidal element model used ten elements in the longitudinal direction. Only the end shear load was considered in these models. Table 4-2 gives the end-tip deflections of these models normalized with respect to the theoretical results. The following parameters were used as input: Length, L = 1.0 m 33 Width, W = 50 mm Thickness, t = 10 mm End Shear Load, V = 5 N Modulus of Elasticity, E = 200 GPa Poisson’s Ratio, ν = 0.0 Theoretical (beam theory) deflection under end shear load, δV = 2.0 mm Table 4-2 Cantilever beam models - parallelogram and trapezoidal element normalized results for end-tip deflection Element Type Integration Scheme 20‐node Solid 24‐node Solid Trapezoidal Models Element Side Angle, φ 0° 45° 75° 0.998 0.917 0.064 0.998 0.994 0.844 Parallel Models Element Side Angle, φ 0° 45° 75° 0.998 0.988 0.567 0.998 0.992 0.870 3x3x3 3x3x3 9 in‐plane, 24‐node Layered 0.998 0.994 0.844 0.998 0.992 0.869 5 through‐thickness 20‐node Solid 2x2x2 1.000 0.985 0.727 1.000 0.997 0.998 4 in‐plane, 24‐node Layered 1.000 0.999 0.891 1.000 1.000 0.985 5 through‐thickness With ten elements in the longitudinal direction, these elements have an aspect ratio of 5. As shown in Table 4-2, at this aspect ratio the 24-node solid and layered elements return values within 2 % of the theoretical result for side angles of 45° or less. Note that the 24node solid element was not run with a 2x2x2 integration scheme, for reasons outlined in the previous section. 4.1.2 Eigenvalue Analysis The cantilever beam models using 24-node solid elements with a 2x2x2 integration scheme produced unstable results. To confirm that this was a modeling problem and not a product of any coding in OOFEM an eigenvalue analysis was performed on the rectangular cantilever model in Section 4.1.1. This analysis computed the eigenvectors of the system’s global stiffness matrix, yielding the natural modes of deformation. The analysis also computed the eigenvalues corresponding to each natural mode. These 34 eigenvalues represent the amount of energy associated with an arbitrary amount of deformation in each corresponding mode shape. Several of the eigenvalues computed were zero or near-zero. A system with zero-energy modes such as these is unstable. Therefore the instability in models using a 2x2x2 integration scheme with the 24-node solid elements is a modeling problem and not an error in OOFEM’s coding. For this reason, all subsequent models using 24-node solid elements were only run with a 3x3x3 integration scheme. 4.1.3 Circular Ring A set of circular ring models were investigated to confirm the 24-node elements’ ability to model curved composite parts. These models consider a ring subjected to two opposing points loads. Applying symmetry conditions, only a quarter of the ring was modeled. Figures 4-5 and 4-6 show the geometry and the finite element model used. Figure 4-5 Circular ring under vertical point loads Figure 4-6 Circular ring model with 4 – 20-node elements in the tangential direction (not to scale) 35 The width of each circular ring model was adjusted to maintain element aspect ratios of 1 in the in-plane dimensions. Table 4-3 summarizes the results for the models considered. Again these values are normalized with respect to the theoretical result. The following parameters were used: Radius, r = 637 mm Thickness, t = 10 mm Vertical Load, P = 1 kN/m Modulus of Elasticity, E = 200 GPa Poisson’s Ratio, ν = 0.0 Theoretical vertical deflection, δ = 2.308 mm Table 4-3 Circular ring models - normalized results for vertical deflection Element Type 20‐node Brick 24‐node Brick 24‐node Layered 20‐node Brick 24‐node Layered Integration Scheme 3x3x3 3x3x3 9 in‐plane, 5 through‐thickness 2x2x2 4 in‐plane, 5 through‐thickness Number of Elements in the Longitudinal Direction 4 0.130 0.130 10 0.851 0.851 20 0.988 0.988 0.130 0.851 0.988 0.992 0.999 1.000 0.998 1.000 1.000 As shown in Table 4-3, the models with four elements produced only 13% of the theoretical deflection when fully integrated. Similar models with ten elements were within 15% of the theoretical value. All models using reduced-integration in-plane were within 1% of the predicted result. Models with four, ten and twenty elements in the tangential direction correspond to aspect ratios of 25, 10 and 5, respectively. Therefore, from the results shown, curved elements using reduced integration schemes can return accurate results for aspect ratios up to 25. Curved elements using full integration schemes should have an aspect ratio of 5 36 or less to maintain accuracy. These are in keeping with previous findings reported in the literature: (Dhondt, 2002), (Arafath, 2007), etc. 4.1.4 Sandwich Beams A set of cantilevered sandwich beam models was investigated to verify the behaviour of the 24-node layered elements. Two sandwich beams were considered: one of constant thickness, and one tapered beam with varying core thickness. Geometries of these two models are shown in Figures 4-7 and 4-8. Figure 4-7 Cantilevered sandwich beam - constant thickness Figure 4-8 Cantilevered sandwich beam – tapered Each beam was analyzed using two element arrangements. The first arrangement used three solid elements in the through-thickness direction: one for each different material. This arrangement was tested for 20 and 24-node solid elements. The second arrangement used a single 24-node layered element in the through-thickness direction, using layers to represent the different materials. Note that the layered elements implemented in OOFEM use a single material only. To account for this, an orthotropic material was created such 37 that its properties in the longitudinal direction matched those of the core material, when the fibre angle was 90°, and those of the skin material when the fibre angle was 0°. Constant Thickness Three models were examined for the cantilever sandwich beam of constant thickness. The first two models used solid elements with three elements in the through-thickness direction and 10 elements in the longitudinal direction. Figure 4-9 shows this element arrangement with the vertical scale exaggerated. The third model used a single layered element in the through-thickness direction. Figure 4-10 shows the element arrangement for this model, again with the vertical scale exaggerated. Figure 4-9 Constant thickness sandwich beam, solid element mesh Figure 4-10 Constant thickness sandwich beam, layered element mesh Figure 4-11 plots the deflection curves along the length of the beam for the three models considered. Each model used the following parameters: Length, L = 100 mm Skin thickness, tskin = 2 mm Core thickness, tcore = 5 mm Uniform load, q = 0.1 kN/m Skin material modulus of elasticity, Eskin = 170 GPa Core material modulus of elasticity, Ecore = 70 GPa 38 0 20 40 x (mm) 60 80 100 120 0.0 Single Layered Element 0.2 3-20 Node Element 0.4 3-24 Node Element Layer δ (μm) 0.6 0.8 1.0 1.2 1.4 1.6 Figure 4-11 Deflection curves for constant thickness sandwich beam models As shown in Figure 4-11 the deflection curves for all three models produce comparable results. The tip deflection from the model using 24-node solid elements is within 0.2 % of the model using 24-node layered elements. These results confirm the suitability of the layered elements with constant layer thickness for modeling the structural behaviour of laminated plates. Tapered Similar to the constant thickness beam, three models were also considered for the tapered sandwich beam. Again, two models were considered using 20- and 24-node solid elements while one model was considered using layered elements. Figures 4-12 and 4-13 show the element arrangements for the solid and layered element models, respectively. The vertical scale in these figures has been exaggerated for illustrative purposes. 39 Figure 4-12 Tapered sandwich beam, solid element mesh Figure 4-13 Tapered sandwich beam, layered element mesh Figure 4-14 plots the deflection curves along the length of the beam for the three models considered. Each model used the following parameters: Length, L = 100 mm Skin thickness, tskin = 2 mm Core thickness at fixed end, tcore_1 = 15 mm Core thickness at free end, tcore_2 = 5 mm Uniform load, q = 0.1 kN/m Skin material modulus of elasticity, Eskin = 170 GPa Core material modulus of elasticity, Ecore = 70 GPa 40 0 20 40 x (mm) 60 80 100 120 0.00 1-24 Node Layered Element 0.05 3-20 Node Solid Elements δ (μm) 0.10 3-24 Node Solid Elements 0.15 0.20 0.25 0.30 0.35 Figure 4-14 Deflection curves for tapered sandwich beam models As shown in Figure 4-14 the deflection curves for all three models produce comparable results. The tip deflection from the model using 24-node solid elements is within 3% of the model using 24-node layered elements. These results confirm the suitability of the layered elements with non-uniform layer thickness for modeling the structural behaviour of laminated plates. 4.2 Thermochemical Models This section outlines the thermochemical models considered. Several 1D models were constructed using a stack of 8-node 3D solid elements with boundary conditions applied to create 1D behaviour. These elements are formulated for thermal analyses; each node has a single temperature degree-of-freedom rather than displacement degrees-of-freedom. Figure 4-15 depicts a model with eight elements and convective boundary conditions. These 8-node elements were chosen to match the element types in proprietary software used for comparison purposes. For further details, see Section 4.2.2. 41 Figure 4-15 1D thermal model with 3D solid elements and convective boundary conditions, side view 4.2.1 Logistic Cure Model To verify the basic implementation of curing materials in OOFEM, as well as test OOFEM’s nonlinear transient solution algorithm, a 1D model with arbitrary cure kinetics was created. In this model, the degree-of-cure was defined purely as a function of time and did not depend on the cure (temperature) cycle. A finite difference model was also created to verify the results from OOFEM. With degree-of-cure defined over the course of the model, creation of a finite difference model was straightforward. Details and results from both the OOFEM and finite difference models are provided in Appendix C. 4.2.2 Single-Hold Cure Cycle Models Implementation of the CCA library interface was tested using a 1D model of a curing composite material with three different cure cycles. These models were run for each of the two materials provided with the CCA library, AS4-3501-6 and AS4-8552 system composites (Materials A and B, respectively). The 1D model is shown in Figure 4-15, with convective boundary conditions at both ends of a slab of material of thickness L. This geometry models a plate of infinite in-plane dimensions, with convective boundaries on both surfaces. The same temperature cure cycle is applied at each boundary. Each cure cycle consisted of a ramp-up phase, hold phase, and cool-down phase. The cure cycles are displayed graphically in Figures 4-16 to 4-21 along with results for the various models. Raven™ is a proprietary software package for 1D thermochemical analysis of composite materials. This program was used to compare temperature histories at the material centre 42 to verify the interface between OOFEM and the CCA library. The software was provided by Convergent, the same company that provided the CCA library. The linear elements used to model these 1D thermochemical problems were selected to match the behaviour of linear elements used in Raven™. Linear elements may be sufficient for representing the through-thickness dimension of a laminated composite part in a 1D thermochemical analysis, but they are not well suited for modeling the coupled thermal-structural problems in Section 4.3. Linear shape functions in the throughthickness direction require the use of elements with relatively low aspect ratios to avoid shear locking. Maintaining low aspect ratios while modeling shell-like structures would require the use of many elements, thereby increasing run-times. The 24-node solid-shell elements introduced in Section 3.1 have higher-order shape functions in the thickness direction. These elements can be used with much higher aspect ratios and are therefore better suited for the coupled thermal-structural models outlined in Section 4.3. For the problems considered here, however, linear elements are adequate. Material A: AS4-3501-6 CFRP The first material tested is described in the CCA interface as AS4-3501-6-2006. This represents a unidirectional composite material consisting of an AS4 carbon fibre and a 3501-6 type resin. Cure kinetics equations and input parameters for this material are provided in Appendix E. Figures 4-16 to 4-18 show temperature histories at mid-plane from the three ramp rate models (r = 1, 2, 3 °C/min) for Material A. The exotherms observed at the beginning of the hold phase in each model are of particular interest and are used for comparison between the OOFEM and Raven™ models. These results are outlined in Table 4-4 below. 43 200 OOFEM and Raven™ 180 160 Cure Cycle Temperature (oC) 140 120 100 80 60 40 20 0 0 50 100 150 200 250 300 350 400 Time (min) Figure 4-16 Temperature history at the mid-plane of Material A as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 1 °C/min 200 OOFEM and Raven™ 180 160 Cure Cycle Temperature (oC) 140 120 100 80 60 40 20 0 0 50 100 150 200 250 300 350 400 Time (min) Figure 4-17 Temperature history at the mid-plane of Material A as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 2 °C/min 44 200 OOFEM and Raven™ 180 Cure Cycle 160 Temperature (oC) 140 120 100 80 60 40 20 0 0 50 100 150 200 250 300 350 400 Time (min) Figure 4-18 Temperature history at the mid-plane of Material A as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 3 °C/min Table 4-4 Results of the 1D 3501 material ramp rate models OOFEM Raven™ OOFEM Raven™ 4.62 12.05 15.92 4.68 12.18 16.08 164.7 85.0 59.6 165.0 85.0 59.5 0.46 0.46 0.54 Ramp Rate r = 1 °C/min r = 2 °C/min r = 3 °C/min Time of Max Exotherm (min) Maximum Absolute Difference (°C) Max Exotherm (°C) As shown in Table 4-4, OOFEM and Raven™ produce similar results for the Material A ramp rate models. The time and magnitude of the exotherms are within 18 seconds and 0.16 °C, respectively. Temperature values for all three cases (ramp rates) are within 0.54 °C throughout the analyses. These results confirm the proper implementation of the CCA interface for thermal modeling. Material B: AS4-8552 CFRP The second material tested is described in the CCA interface as AS4-8552-2007. This represents a composite material with an AS4 carbon fibre along with an 8552 type resin. 45 Cure kinetics equations and input parameters for this material are again provided in Appendix E. Figures 4-19 to 4-21 show the temperature histories at mid-plane from the three ramp rate models for Material B. Again the exotherms observed at the beginning of the hold phase in each model are used for comparison purposes. These results are outlined in Table 4-5. Table 4-5 Results of the 1D 8552 material ramp rate models Time of Max Exotherm (min) Max Exotherm (°C) OOFEM Raven™ OOFEM Raven™ Maximum Absolute Difference (°C) 8.65 22.09 24.79 7.40 20.08 23.46 162.2 87.3 64.1 163.0 86.5 63.0 1.69 2.55 3.56 Ramp Rate r = 1 °C/min r = 2 °C/min r = 3 °C/min 220 OOFEM and Raven 200 180 160 Cure Cycle Temperature (oC) 140 120 100 80 60 40 20 0 0 50 100 150 200 250 300 350 400 Time (min) Figure 4-19 Temperature history at the mid-plane of Material B as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 1 °C/min 46 220 OOFEM and Raven 200 180 160 Cure Cycle Temperature (oC) 140 120 100 80 60 40 20 0 0 50 100 150 200 250 300 350 400 Time (min) Figure 4-20 Temperature history at the mid-plane of Material B as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 2 °C/min 220 OOFEM and Raven 200 180 160 Cure Cycle Temperature (oC) 140 120 100 80 60 40 20 0 0 50 100 150 200 250 300 350 400 Time (min) Figure 4-21 Temperature history at the mid-plane of Material B as predicted by Raven™ and OOFEM for a cure cycle with ramp rate r = 3 °C/min 47 As shown in Table 4-5, OOFEM and Raven™ results for the three ramp rates show slight differences. The time and magnitude of the exotherms differ by up to 66 seconds and 2.0 °C, respectively. Temperature values for all three models are within 3.56 °C throughout the analyses. The OOFEM implementation handles both materials A and B in the same manner. As such, the differences observed in the results for Material B are difficult to explain. At the time of writing, these differences remain unresolved. 4.3 Coupled Thermal-Structural Models 4.3.1 Flat Composite Part on an Aluminum Tool Four coupled thermal-structural models were considered in this work. Each model investigated a flat composite part on an aluminum tool in a 2D domain. The models represent each combination of either a thin part or thick part, and Dirichlet (specified temperature) or Von Neumann (convective) thermal boundary conditions. The part/tool assembly is geometrically symmetric about its centreline. As such, symmetric structural boundary conditions are used in all four models. Figure 4-22 shows the geometry and structural boundary conditions used in these models. Figure 4-22 Flat part and tool models - geometry and structural boundary conditions 48 Figure 4-23 shows the parametric mesh arrangement used for these models. Figure 4-23 Flat part and tool models - meshing and structural boundary conditions Note that all models in this section consist of a part made of Material B on an aluminum tool. Poisson’s ratios for both the tool and part were set to zero in these models. This condition was necessary to satisfy assumptions made in deriving the closed-form solution outlined below. The following parameters were used for models in this section: All Models Assembly half-length, L = 300 mm Part material properties: (see Appendix E) Tool material modulus of elasticity, Etool = 70 GPa Tool material shear modulus, Gtool = 70 GPa Tool material coefficient of thermal expansion, αtool = 23.6·10-6 °C-1 Tool material thermal conductivity, ktool = 180 W / (m°C) Tool material specific heat capacity, ctool = 1250 J/(kg°C) Tool material density, ρtool = 2710 kg/m3 Thin Part Models Part thickness, t2 = 1.6 mm Tool thickness, t1 = 5 mm Number of elements vertically through the part, npart = 4 Number of elements vertically through the tool, ntool = 2 49 Number of elements longitudinally, m = 13 Thick Part Models Part thickness, t2 = 12 mm Tool thickness, t1 = 25 mm Number of elements vertically through the part, npart = 6 Number of elements vertically through the tool, ntool = 3 Number of elements longitudinally, m = 13 The numbers of elements in both the longitudinal and thickness directions were varied to assess convergence of the results. Increasing the number of elements beyond those shown above did not have any significant effect on the results considered. Closed-Form Solution A closed-form solution derived by Arafath (Arafath, 2007), (Arafath, Vaziri, & Poursartip, 2008) was used for comparison to the models with specified temperature thermal boundary conditions. This solution calculated displacements and stresses for two fully-bonded beams subjected to a temperature change. The solution assumes a uniform temperature throughout the system. For this reason, thermal boundary conditions in the OOFEM model were applied such that all temperature degrees-of-freedom followed the applied cure cycle. Material properties for the part were determined using equations taken from(Johnston, 1997). In this framework, finite difference methods were used to determine degree-ofcure values throughout the analysis. These degree-of-cure results were input into the material property equations to determine the mechanical behaviour of the part over time. Stresses and strains in the system were updated incrementally over time. This combination of closed-form structural calculation with numerical time integration was implemented in GNU Octave, an open-source program designed to be compatible 50 with Matlab™ m-files. Results from this implementation are referred to as ‘Octave Results’ in the sections that follow. Thin Part – Specified Temperature Boundary Conditions The first coupled models investigated combined a thin composite part with specified temperature thermal boundary conditions. Figure 4-24 shows the applied cure cycle along with degree-of-cure results for both the OOFEM model and the closed-form solution implemented in Octave. Degree-of-cure results in OOFEM and Octave are within 0.45% of each other throughout the event. Figures 4-25 and 4-26 give the development of E11 and G12 for a single unidirectional ply over the course of the cure cycle. Values of E11 were within 0.02% throughout the event. Values of G12 differed by a maximum of 3.25% during the period of rapid modulus 200 1.0 180 0.9 160 0.8 140 0.7 120 0.6 100 0.5 Cure Cycle 0.4 80 Degreee of Cure Temperature (°C) development and returned to a numerically insignificant difference afterwards. 0.3 60 OOFEM and Octave 40 0.2 20 0.1 0 0.0 0 50 100 150 200 250 300 350 Time (min) Figure 4-24 Thin part model with specified temperature thermal boundary conditions – degree-ofcure results 51 122.5 200 180 122.0 140 121.5 120 Cure Cycle 100 121.0 80 60 OOFEM and Octave 40 Ply E11 (GPa) Temperature (°C) 160 120.5 20 120.0 0 0 50 100 150 200 250 300 350 Time (min) Figure 4-25 Thin part model with specified temperature thermal boundary conditions - ply E11 values 200 6.0 180 5.0 140 4.0 120 Cure Cycle 100 3.0 80 Ply G12 (GPa) Temperature (°C) 160 2.0 60 OOFEM and Octave 40 1.0 20 0.0 0 0 50 100 150 200 250 300 350 Time (min) Figure 4-26 Thin part model with specified temperature thermal boundary conditions - ply G12 values 52 Warpage of the composite part after removal from the tool is a function of the residual stresses present at the end of the cure cycle. The version of OOFEM used in this research (OOFEM version 1.8) was not capable of deleting elements or changing boundary conditions during an analysis. As such, the stress profile through the part was the output of primary interest in these models. Figure 4-27 plots the axial stress profile through the part at the mid-point of the length modeled. These are the axial stresses present at 100 minutes into the cure cycle, 20 minutes into the hold phase. The two curves show consistent profiles through the part thickness. The maximum difference in axial stress between the two models is 1.1% of the overall stress gradient observed through the thickness. Figure 4-28 plots the axial stress profile at the same location at the end of the cure cycle. These are the residual stresses present upon removal from the tool. Again, the two curves show consistent results throughout the part. The maximum difference in axial stress between the two models is 1.0% of the overall stress gradient observed through the thickness. 1.80 1.60 Octave Results Y‐Coordinate (mm) 1.40 OOFEM Results 1.20 1.00 0.80 0.60 0.40 0.20 0.00 150.0 200.0 250.0 300.0 350.0 Axial Stress (MPa) Figure 4-27 Thin part model with specified temperature thermal boundary conditions - axial stress through the part at x = 150 mm, t = 100 minutes 53 1.80 Octave Results Y‐Coordinate (mm) OOFEM Results 1.60 1.40 1.20 1.00 0.80 0.60 0.40 0.20 0.00 ‐150.0 ‐100.0 ‐50.0 0.0 50.0 Axial Stress (MPa) Figure 4-28 Thin part model with specified temperature thermal boundary conditions - axial stress through the part at x = 150 mm, t = 320 minutes Note that the curves in Figures 4-27 and 4-28 are similar in shape. As shown in Figures 4-24 to 4-26, degree-of-cure, E11 and G12 all experience major increases during the hold phase. In fact, for this assembly, all changes to the mechanical material properties occur during the hold phase. Therefore, the axial stress profile present at the beginning of the hold phase will have the same shape as that at the end of the cure cycle. The final residual stresses will be shifted toward the compressive side of the stress axis relative to those at the beginning of the hold phase. This is caused by differential contraction of the tool and part during the cool-down phase. Overall the closed-form solution results implemented in Octave match very well with those from OOFEM. These results confirm OOFEM’s implementation of the CCA library interface and its methodology of calculating composite material properties during a cure cycle. 54 Thin Part – Convective Boundary Conditions The second coupled models investigated used the same geometry as the first set, but with convective thermal boundary conditions applied. The top of the part was defined as a convective surface while all other surfaces were considered insulated. Figure 4-29 shows the temperature results for this model as well as the previous model with specified temperature BCs. These results are taken at the quarter-point of the assembly (x = 150 mm) at 0.6 mm from the top of the part (y = 1.0 mm). 200 180 Temperature (°C) 160 140 120 100 80 Convective Boundary Conditions 60 Specified Temperature Boundary Conditions 40 20 0 0 50 100 150 200 250 300 350 Time (min) Figure 4-29 Thin part models - temperature histories at x = 150 mm, y = 1 mm Temperatures in the model with convective BCs follow a very similar pattern to the model with specified temperature BCs (where all temperatures match the cure cycle). For the most part, temperatures in this second model match those from the first with a consistent lag in time. This lag in time is explained by the addition of the convective surface at the top of the part. Boundary conditions of this type mean that temperature changes in the assembly will lag behind ambient temperatures changes. A small exotherm is also observed at the beginning of the cure cycle hold phase. 55 With a thickness of 1.6 mm there is very little temperature gradient through the part. As such, this model behaves similarly to a model with a uniform temperature throughout (as was the case with specified temperature boundary conditions). Near-uniform temperatures occurring with a small time lag created conditions very similar to the first model. Therefore, it was not surprising that stress, displacement and degree-of-cure results from this model were almost identical to those from the first model with the same time-lag. Figure 4-30 shows the degree-of-cure evolution at the quarter-point location for this case 0.90 200 0.80 180 0.70 160 140 0.60 120 0.50 100 0.40 80 Convective Boundary Conditions Specified Temperature Boundary Conditions Cure Cycle 0.30 0.20 0.10 Temperature (°C) Degree‐of‐Cure as well as the previous one with specified temperature boundary conditions. 60 40 20 0 0.00 0 50 100 150 200 250 300 350 Time (min) Figure 4-30 Thin part models – degree-of-cure history at x = 150 mm, y = 1 mm As expected, the degree-of-cure results follow the previous case results very closely, with the most rapid cure rate occurring slightly later in the model. Figure 4-31 shows the axial stress development at x = 150 mm. Once again there is a time-lag between the axial stresses obtained from the two cases. During the exotherm, axial stresses obtained from the case with convective boundary conditions temporarily 56 exceed those obtained from the previous case with specified temperature boundary conditions. 250.0 200.0 Axial Stress (MPa) 150.0 100.0 50.0 0.0 0 50 100 150 200 ‐50.0 Convective Boundary Conditions ‐100.0 Specified Temperature Boundary Conditions ‐150.0 250 300 350 Time (min) Figure 4-31 Thin part models - axial stress history at x = 150 mm, y = 1 mm Figure 4-32 presents axial stresses through the part at t = 60 minutes. With temperatures through the part being nearly uniform, the stresses in both cases follow a similar profile. 1.80 Convective Boundary Conditions Specified Temperature Boundary Conditions 1.60 Y‐Coordinate (mm) 1.40 1.20 1.00 0.80 0.60 0.40 0.20 0.00 0.0 50.0 100.0 150.0 200.0 250.0 Axial Stress (MPa) Figure 4-32 Thin part models - axial stress through the part at x = 150 mm, t = 60 minutes 57 Overall the results from the thin part model with convective boundary conditions closely matched those with specified temperature boundary conditions, separated by an explainable time-lag. The few differences observed, including the small exotherm, were also expected effects of the boundary conditions applied. Thick Part – Specified Temperature Boundary Conditions The next coupled models combined a thick composite part with specified temperature thermal boundary conditions. OOFEM results again compared favourably to the Octave implementation of the closed-form solution. Given that temperature values are uniform throughout the assembly and follow the cure cycle exactly, the degree-of-cure results from this model are indistinguishable from those in Figure 4-24. Degree-of-cure results in OOFEM and Octave agree to within 1.14% throughout the model. With composite material properties a function of both temperature and degree-of-cure, the values of E11 and G12 in this model are also indistinguishable from those shown in Figures 4-25 and 4-26. Axial stresses through the part followed a different profile from that in the thin part model. Figure 4-33 shows these results from both Octave and OOFEM. Again the OOFEM and Octave results agree through the thickness. Note that the part used in this model is sufficiently thick that, as the cure cycle hold phase begins, the top fibre of the part has zero axial stress. This condition results from the composite material’s low shear modulus during the initial stages of the cure cycle. 58 14.0 Y‐Coordinate (mm) 12.0 10.0 Octave OOFEM 8.0 6.0 4.0 2.0 0.0 ‐50.0 0.0 50.0 100.0 150.0 200.0 250.0 300.0 350.0 400.0 450.0 Axial Stress (MPa) Figure 4-33 Thick part model with specified temperature thermal boundary conditions - axial stress through the part at x = 150 mm, t = 100 minutes Large increases in degree-of-cure during the cure cycle hold phase cause a rapid rise in the axial and shear stiffness of the part. As a result, the composite material achieves its final material properties before the beginning of the cool-down phase. The relative stress profile shown in Figure 4-33 is thereby ‘locked in’ and maintains its shape during the cool-down phase. Figure 4-34 shows the axial stress profile through the part at the end of the cure cycle. Note that the profile shown in Figure 4-34 is similar in shape to that in Figure 4-33 with a significant shift to the compressive side of the stress axis. This compression is imposed on the part by differential contraction of the tool during the cooldown phase. 59 14.0 12.0 Octave Y‐Coordinate (mm) 10.0 OOFEM 8.0 6.0 4.0 2.0 0.0 ‐300.0 ‐250.0 ‐200.0 ‐150.0 ‐100.0 ‐50.0 0.0 Axial Stress (MPa) 50.0 100.0 150.0 200.0 Figure 4-34 Thick part model with specified temperature thermal boundary conditions - axial stress through the part at x = 150 mm, t = 320 minutes The results shown above again verify OOFEM’s implementation of the CCA library interface and its methodology for modeling composite material cure kinetics. These models also illustrate the effects of increasing part thickness on the resulting throughthickness residual stress gradients. Thick Part – Convective Boundary Conditions The final coupled models investigated combined a thick part with convective boundary conditions. Figure 4-35 shows the temperature results for this model as well as the previous model with specified temperature boundary conditions. These results are taken at a point x = 150 mm, y = 7.0 mm (5.0 mm from the top surface of the part). 60 200 180 Temperature (°C) 160 140 120 100 80 Convective Boundary Conditions 60 Specified Temperature Boundary Conditions 40 20 0 0 50 100 150 200 250 300 350 Time (min) Figure 4-35 Thick part models - temperature history at x = 150 mm, y = 7 mm Temperatures in this model are a stark contrast to the model with specified temperature boundaries (where all temperatures match the cure cycle). The ambient heat applied by the cure cycle takes considerable time to penetrate the thick composite part. At a depth of 5.0 mm below the surface of the part the temperature peaks at 178.7 °C, having never reached the cure cycle hold temperature of 180.0 °C. With a thickness of 12.0 mm there is significant temperature gradient through the part. As such, degree-of-cure values and mechanical properties will also vary through the thickness. Figure 4-36 shows the temperature results through the thickness of the part as the cure cycle reaches the hold phase. 61 Y‐Coordinate (mm) 14.00 12.00 Convective Boundary Conditions 10.00 Specified Temperature Boundary Conditions 8.00 6.00 4.00 2.00 0.00 0.0 50.0 100.0 150.0 200.0 Temperature (°C) Figure 4-36 Thick part models - temperature through the part at x = 150 mm, t = 80 minutes As shown in Figure 4-36, the temperature at the top surface of the part lags behind the cure cycle temperature by 51.5 °C, while the temperature gradient through the part is 42.1 °C. Figure 4-37 shows the degree-of-cure results for this model as well as the previous model. These results were taken at the same point as those above. 62 200 0.90 180 0.80 160 0.70 140 0.60 120 0.50 100 0.40 80 0.30 60 Convective Boundary Conditions Specified Temperature Boundary Conditions Cure Cycle 0.20 0.10 Temperature (°C) Degree‐of‐Cure 1.00 40 20 0 0.00 0 50 100 150 200 250 300 350 Time (min) Figure 4-37 Thick part models – degree-of-cure history at x = 150 mm, y = 7 mm Here again the model with convective boundary conditions shows a stark contrast to the previous model. Degree-of-cure values lag well behind the model with specified temperature boundary conditions and, in fact, do not reach the same final value. 50.0 0.0 0 50 100 150 200 250 300 350 Axial Stress (MPa) ‐50.0 ‐100.0 ‐150.0 Convective Boundary Conditions Specified Temperature Boundary Conditions ‐200.0 ‐250.0 ‐300.0 Time (min) Figure 4-38 Thick part models - axial stress history at x = 150 mm, y = 7 mm 63 Figure 4-38 gives the axial stress at the same position in the part over the course of the cure cycle. Note that there was little to no axial stress development at this position in the assembly until late in the hold phase. Figures 4-39 and 4-40 reveal why this was the case. Figure 4-39 depicts the axial stresses through the part 20 minutes into the hold phase. These stresses are significant near the part/tool interface, but quickly die down as the ycoordinate increases. This behaviour is a result of the low shear modulus in the composite material at low values of degree-of-cure. 14.00 Convective Boundary Conditions Y‐Coordinate (mm) 12.00 Specified Temperature Boundary Conditions 10.00 8.00 6.00 4.00 2.00 0.00 ‐50.0 0.0 50.0 100.0 150.0 200.0 250.0 300.0 350.0 400.0 Axial Stress (MPa) Figure 4-39 Thick part models - axial stress through the part at x = 150 mm, t = 100 minutes Once the part has undergone a rapid increase in degree-of-cure, the axial stress profile is ‘locked in’ and the profile maintains its shape during the cool-down phase. Figure 4-40 gives the axial stress profile present at the end of the cure cycle (t = 320 minutes). As Figure 4-35 illustrates, temperatures in the model with convective boundary conditions did not cool to ambient conditions at the end of the model. Thus the stress profile given in Figure 4-40 does not represent the ultimate axial stress profile after the part has cooled to 20 °C. 64 Y‐Coordinate (mm) 14.00 12.00 Convective Boundary Conditions 10.00 Specified Temperature Boundary Conditions 8.00 6.00 4.00 2.00 0.00 ‐300.0 ‐200.0 ‐100.0 0.0 100.0 200.0 300.0 Axial Stress (MPa) Figure 4-40 Thick part models - axial stress through the part at x = 150 mm, t = 320 minutes 65 Chapter 5 5.1 Conclusion and Future Directions Conclusions The capability of creating coupled thermal-structural models of curing composite materials has been added to OOFEM. This implementation has been tested and verified to an acceptable level of accuracy for FE models. Given the time commitment and level of effort involved in adding these capabilities, the benefits of working with an open-source code are clear. With large portions of basic program structure already in place, selecting a relatively-mature open-source program for extension can save considerable time for academic researchers and developers. Further benefits arise if the program has an established support utility, such as OOFEM’s online user’s forum. Finally, researchers using open-source programs have the potential to gain from additions made to the program by other developers. This situation did occur over the course of this work. In one specific instance, an inquiry was made to the program’s maintainer, Dr. Patzak, asking if OOFEM was capable of reading a matrix directly from the input file. Dr. Patzak’s response was to program this capability himself before replying to the inquiry. Table 5-1 in Section 5.2 contains another example wherein contributions by other programmers will aid in the development of this project. This work has also demonstrated the suitability of analyzing composite process problems in a coupled manner. Previous work in this area solved the heat equations throughout the entire model first, thereby determining temperatures and material properties for all time steps in the model. Those results were then applied to solve the structural problem. In this work, both the thermal and structural equations are solved in each time step, before proceeding to the next step. 5.2 Future Directions The framework for coupled thermal-structural models implemented in this work represents a first step in seamlessly modeling composite structures from manufacture to 66 failure or decommissioning. A number of steps remain before the goal of seamless modeling can be realized. These steps are outlined in Table 5-1. Table 5-1 Future work for seamless modeling in OOFEM Task Comment This functionality has been included in OOFEM release 1.9. The code associated with this capability must be tested and adapted Add functionality to delete elements and change boundary such that it fits into the framework of the coupled analysis module. This capability will allow the direct prediction of warpage and conditions during a model. residual stresses imposed on the part after removal from the tool. Implement dynamic time‐ stepping. The ability to dynamically change the size of each time step will dramatically reduce run‐times. Small time steps are only required during periods of rapidly changing degree‐of‐cure. These periods represent a small portion of the overall model. Therefore, the ability the increase time step sizes throughout the bulk of the model will result in much shorter run times. Implement interface surface capabilities This functionality will enable the modeling of a part/tool interface that isn't perfectly bonded. Program ability to save the stress profile of a part at the end of a process model for use with future structural models. This functionality will enable the warpage and residual stress profile present at the end of a process model to be used for multiple subsequent structural models. With this framework established, process models for a given part will only need to be run once. The part can then be subjected to an arbitrary number of structural models. Implement materials with A necessary step to facilitate seamless modeling. damage mechanics capabilities. Implement ability to import a saved warped geometry and residual stress profile as the starting point for a structural model. Enables the results from a single process model to be used in multiple structural models. This emulates reality where a single finished composite part will be subjected to many load cycles over time. General code optimization. Implementation of the current code for composites modeling has focused on adding capability. The program would benefit from a project designed to optimize performance of the code. Add preprocessing and/or postprocessing facilities. OOFEM input and output is currently based on text files. The program would benefit from either dedicated facilities for preprocessing and postprocessing, or an interface to generic programs for preprocessing and postprocessing. A framework such as this would save users and developers from having to generate their own preprocessing and postprocessing procedures. 67 Bibliography Arafath, A. R. (2007). Efficient Numerical Techniques for Predicting Process-Induced Stress and Deformations in Composite Structures. Vancouver: The University of British Columbia. Arafath, A. R., Vaziri, R., & Poursartip, A. (2008). Closed-form solution for processinduced stresses and deformation of a composite part cured on a solid tool: Part 1 - Flat geometries. Composites: Part A , 39 (7), 1106-1117. Convergent Manufacturing Technologies Inc. (2009). Retrieved May 24, 2010, from Convergent: http://www.convergent.ca/ Dhondt, G. (2002). The Right Solid Element: a Challenge to Industry. Proceedings of the Fifth World Congress on Computational Mechanics (WCCM V). Vienna. Électricité de France. (2007, November). Code_Aster. Retrieved January 15, 2009, from http://www.code-aster.org/ Free Software Foundation. (1991, June). GNU General Public License v2.0. Retrieved January 15, 2009, from http://www.gnu.org/licenses/old-licenses/gpl-2.0.html Johnston, A. A. (1997). An Integrated Model of the Development of Process-Induced Deformation in Autoclave Processing of Composite Structures. Vancouver: The University of British Columbia. Logg, A. (2009, 01 12). The FEniCS Project. Retrieved January 15, 2009, from http://www.fenics.org/ Patzak, B. (2008, November 26). OOFEM - free object oriented finite element code. Retrieved January 15, 2009, from http://www.oofem.org/ Rasekh, A. (2007). Efficient Methods for Non-Linear Thermochemical Analysis of Composite Structures Undergoing Autoclave Processing. Vancouver: The University of British Columbia. Roddeman, D. (2001, August 19). TOCHNOG - a free explicit/implicit finite element program. Retrieved January 15, 2009, from http://tochnog.sourceforge.net/ Stewart, J. (2003). Calculus. Belmont, California, USA: Brooks/Cole-Thomson Learning. Yevick, D. (2005). A First Course in Computational Phsyics and Object-Oriented Programming with C++. Cambridge University Press. 68 Appendix A OOFEM Class Hierarchy Shown below is the textual class hierarchy for OOFEM, as generated by the program doxygen. This inheritance list is sorted roughly, but not completely, alphabetically: AList< T > ApplicationOrdering PetscOrdering_Base PetscNatural2GlobalOrdering PetscNatural2LocalOrdering c_fun CommunicationBuffer CommunicationPacket Communicator FETICommunicator ProblemCommunicator Communicator< PNlDEIDynamic > PNlDEIDynamicComunicator CommunicatorBuff ConnectivityTable ContextIOERR CreepData DataReader BufferedDataReader OOFEMTXTDataReader PlainTextDataReader Dictionary Dof HangingDof MasterDof RemoteMasterDof SharedMasterDof NullDof RigidArmSlaveDof SlaveDof Domain dynaList< T > dynaListIterator< T > EngngModel CBS EigenValueDynamic NonStationaryTransportProblem NLTransientTransportProblem StaggeredProblem StationaryFlow StationaryTransportProblem 69 StructuralEngngModel DEIDynamic DIIDynamic IncrementalLinearStatic LinearStability LinearStatic AdaptiveLinearStatic NonLinearStatic AdaptiveNonLinearStatic PLinearStatic NlDEIDynamic PNlDEIDynamic SUPG EngngModelContext ExportModule POIExportModule VTKExportModule ExportModuleManager FEInterpolation FEInterpolation1d FEI1dLin FEInterpolation2d FEI2dQuadLin FEI2dTrLin FEInterpolation3d FEI3dHexaLin FEI3dTrLin FEMComponent CrossSection EmptyCS HeatCrossSection StructuralCrossSection FiberedCrossSection LayeredCrossSection SimpleCrossSection DofManager ElementSide Node HangingNode RigidArmNode Element FluidDynamicElement FMElement CBSElement TR1_2D_CBS SUPGElement TR1_2D_SUPG 70 TR1_2D_SUPG2 TR1_2D_SUPG2_AXI TR1_2D_SUPG_AXI PPdeElement LTrElementPPDE StructuralElement Beam2d Beam3d Brick24 InterfaceElem1d InterfaceElem2dQuad LIBeam2d LIBeam3d LSpace LTRSpace NLStructuralElement Axisymm3d CCTPlate RerShell L4Axisymm LIBeam2dNL LIBeam3d2 LIBeam3dNL LIBeam3dNL2 PlaneStress2d TrPlaneStress2d TrPlaneStrRot Truss2d Truss3d Q4Axisymm QPlaneStress2d QSpace QTrPlaneStress2d Quad1PlaneStrain TrPlaneStrain Truss1d TransportElement Brick1_ht Quad1_ht QuadAxisym1_ht Tetrah1_ht Tr1_ht TrAxisym1_ht ErrorEstimator CombinedZZSIErrorEstimator HuertaErrorEstimator ScalarErrorIndicator 71 ZZErrorEstimator GaussPoint Microplane GeneralBoundaryCondition BoundaryCondition TractionPressureBC Load BodyLoad DeadWeight BoundaryLoad ConstantEdgeLoad ConstantSurfaceLoad LinearEdgeLoad NodalLoad PointLoad StructuralTemperatureLoad TF1 UserDefinedTemperatureField InitialCondition IntegrationRule GaussIntegrationRule LobattoIntegrationRule LoadTimeFunction ConstantFunction HeavisideLTF PeakFunction PiecewiseLinFunction PeriodicPiecewiseLinFunction UserDefinedLoadTimeFunction Material FluidDynamicMaterial BinghamFluidMaterial2 NewtonianFluidMaterial TwoFluidMaterial HeatMaterial IsotropicLinearHeatMaterial HydrationModel StructuralMaterial CebFipSlip90Material DeformationTheoryMaterial Concrete2 DruckerPragerPlasticitySM HellmichMaterial IsoInterfaceDamageMaterial IsoInterfaceDamageMaterial IsotropicDamageMaterial IsotropicDamageMaterial1 72 IDNLMaterial MazarsMaterial MazarsNLMaterial LinearElasticMaterial IsotropicLinearElasticMaterial AgingIsoLEMaterial OrthotropicLinearElasticMaterial MaxwellChainMaterial B3Material CebFip78Material DoublePowerLawMaterial MicroplaneMaterial MDM MicroplaneMaterial_Bazant M4Material MPlasticMaterial J2MPlasticMaterial RankinePlasticMaterial MPlasticMaterial2 J2Mat Masonry02 PerfectlyPlasticMaterial Steel1 PlasticMaterial J2plasticMaterial RCM2Material Concrete3 RCSDEMaterial RCSDNLMaterial RCSDMaterial TransportMaterial HeMoTKMaterial HydratingHeMoMaterial IsotropicHeatTransferMaterial HydratingIsoHeatMaterial MaterialInterface LEPlic MaterialStatus FluidDynamicMaterialStatus BinghamFluidMaterial2Status HydrationModelStatus StructuralMaterialStatus CebFipSlip90MaterialStatus Concrete2MaterialStatus DruckerPragerPlasticitySMStatus HellmichMaterialStatus IsoInterfaceDamageMaterialStatus 73 IsoInterfaceDamageMaterialStatus IsotropicDamageMaterialStatus IsotropicDamageMaterial1Status IDNLMaterialStatus MazarsMaterialStatus MazarsNLMaterialStatus M4MaterialStatus MaxwellChainMaterialStatus MDMStatus MPlasticMaterial2Status MPlasticMaterialStatus PerfectlyPlasticMaterialStatus PlasticMaterialStatus RCM2MaterialStatus RCSDEMaterialStatus RCSDNLMaterialStatus RCSDMaterialStatus TransportMaterialStatus HydratingTransportMaterialStatus NonlocalBarrier PolylineNonlocalBarrier SymmetryBarrier NumericalMethod GJacobi LineSearchNM SparseGeneralEigenValueSystemNM InverseIteration SubspaceIteration SparseLinearSystemNM DSSSolver FETISolver IMLSolver LDLTFactorization PetscSolver SpoolesSolver SparseNonLinearSystemNM CylindricalALM NRSolver NRSolver2 RemeshingCriteria CombinedZZSIRemeshingCriteria DirectErrorIndicatorRC HuertaRemeshingCriteria ZZRemeshingCriteria SpatialLocalizer DummySpatialLocalizer OctreeSpatialLocalizer 74 FETIBoundaryDofManager Field PrimaryField DofDistributedPrimaryField FieldManager FloatArray Column StressStrainBaseVector StrainVector StressVector Graph Graph::node GT_Exception InputRecord OOFEMTXTInputRecord IntArray Interface DirectErrorIndicatorRCInterface LTRSpace PlaneStress2d QTrPlaneStress2d Quad1PlaneStrain TrPlaneStrain TrPlaneStress2d Truss1d Truss3d EIPrimaryFieldInterface TR1_2D_CBS TR1_2D_SUPG TransportElement EIPrimaryUnknownMapperInterface LSpace LTRSpace PlaneStress2d QTrPlaneStress2d Quad1PlaneStrain TrPlaneStrain TrPlaneStress2d Truss1d FiberedCrossSectionInterface Beam3d LIBeam3d LIBeam3d2 HuertaErrorEstimatorInterface LSpace LTRSpace PlaneStress2d 75 Quad1PlaneStrain TrPlaneStrain TrPlaneStress2d Truss1d HuertaRemeshingCriteriaInterface LSpace LTRSpace PlaneStress2d Quad1PlaneStrain TrPlaneStrain TrPlaneStress2d Truss1d HydrationModelInterface HellmichMaterial HydratingHeMoMaterial HydratingIsoHeatMaterial HydrationModelStatusInterface HellmichMaterialStatus HydratingTransportMaterialStatus LayeredCrossSectionInterface Beam2d CCTPlate LIBeam2d LIBeam2dNL LEPlicElementInterface TR1_2D_SUPG MaterialModelMapperInterface IsotropicDamageMaterial1 MDM MMAShapeFunctProjectionInterface LTRSpace TrPlaneStrain TrPlaneStress2d Truss1d NodalAveragingRecoveryModelInterface Axisymm3d CCTPlate LSpace LTRSpace TR1_2D_CBS TR1_2D_SUPG TrPlaneStrain TrPlaneStress2d Truss1d NonlocalMaterialExtensionInterface StructuralNonlocalMaterialExtensionInterface IDNLMaterial 76 MazarsNLMaterial MDM RCSDNLMaterial NonlocalMaterialStatusExtensionInterface StructuralNonlocalMaterialStatusExtensionInterface IDNLMaterialStatus MazarsNLMaterialStatus MDMStatus RCSDNLMaterialStatus NonlocalMaterialStiffnessInterface IDNLMaterial SpatialLocalizerInterface Brick1_ht L4Axisymm LSpace LTRSpace PlaneStress2d QTrPlaneStress2d Quad1_ht Quad1PlaneStrain Tetrah1_ht TR1_2D_CBS TR1_2D_SUPG Tr1_ht TrPlaneStrain TrPlaneStress2d Truss1d SPRNodalRecoveryModelInterface Axisymm3d CCTPlate L4Axisymm LSpace LTRSpace PlaneStress2d QTrPlaneStress2d Quad1PlaneStrain TR1_2D_CBS TR1_2D_SUPG TrPlaneStrain TrPlaneStress2d ZZErrorEstimatorInterface LTRSpace TrPlaneStrain TrPlaneStress2d Truss1d ZZNodalRecoveryModelInterface Axisymm3d 77 CCTPlate L4Axisymm LSpace LTRSpace PlaneStress2d Quad1PlaneStrain TR1_2D_CBS TR1_2D_SUPG TrPlaneStrain TrPlaneStress2d Truss1d ZZRemeshingCriteriaInterface LTRSpace TrPlaneStrain TrPlaneStress2d Truss1d LEPlic::computeLEPLICVolumeFractionWrapper listItem< T > localIntegrationRecord Logger MaterialMappingAlgorithm MMAClosestIPTransfer MMALeastSquareProjection MMAShapeFunctProjection Matrix FloatMatrix SparseMtrx CompCol SymCompCol DSSMatrix DynCompCol DynCompRow PetscSparseMtrx Skyline SkylineUnsym mem_fun< T > MesherInterface FreemInterface T3DInterface Targe2Interface MetaStep NodalRecoveryModel NodalAveragingRecoveryModel SPRNodalRecoveryModel ZZNodalRecoveryModel NonisoData OOFEM_Terminate 78 oofemOctantRec OutputManager Pair Parser Parser::name PetscContext PlastData POIExportModule::POI_dataType Polygon Polygon::PolygonEdgeIterator Polygon::PolygonVertexIterator Preconditioner CompCol_ICPreconditioner CompCol_ILUPreconditioner CompRow_ILUPreconditioner DiagPreconditioner VoidPreconditioner PrimaryUnknownMapper EIPrimaryUnknownMapper ProcessCommunicator ProcessCommunicatorBuff Range RefinedElement RefinedMesh RowColumn SloanGraph SloanGraphNode SloanLevelStructure SloanNodalDegreeOrderingCrit StringReader TDictionary< Key, T > TDictionaryIterator< Key, T > TimeStep Tokenizer TPair< Key, T > Vertex 79 Appendix B Nonstationary Transient Solving Algorithm This appendix outlines the solving algorithm used for nonstationary transient transport problems in OOFEM. Details of this algorithm were provided via correspondence with Dr. Borek Patzák. The algorithm is used to solve Equation (B-1) below. , , , (B-1) where: = heat capacity = material density = thermal conductivity = temperature = heat generation term = spatial coordinate = time Approximate solutions to Equation (B-1) are produced by solving the discretized finite element Equation (B-2). (B-2) where: = capacitance matrix = conductivity matrix = temperature vector = internal source vector = time Time discretization is based on the generalized trapezoidal rule. For a given time solution vector the solution is approximated on time interval , ∆ with known by Equation (B-3). 80 (B-3) ∆ where: ∆ 0,1 The time derivative is expressed simply as Equation (B-4). (B-4) ∆ ∆ Substituting Equations (B-3) and (B-4) into (B-2) produces the following balance equation for time . ∆ (B-5) ∆ ∆ The system expressed in Equation (B-5) is generally nonlinear, so a Newton-Raphson iterative solving method is employed. To begin, Equation (B-5) is first expressed in residual form. ∆ ∆ ∆ ∆ (B-6) where: = residual vector at time Linearizing Equation (B-6) produces an expression for iterative change of the residual vector. ∆ ∆ ∆ ∆ ∆ ∆ (B-7) where: = iteration number Substituting Equation (B-6) into (B-7), evaluating and rearranging produces the iterative solution algorithm shown in Equation (B-8) below. 81 ∆ ∆ ∆ ∆ (B-8) ∆ ∆ ∆ ∆ The system of equations in (B-8) is applied iteratively until convergence is achieved. The system is considered converged once ∆ given time step OOFEM will set ∆ ∆ is within a predefined tolerance. For the first iteration in a ∆ . Alternatively the first iteration could take . 82 Appendix C Logistic Cure Model To verify the basic implementation of curing materials in OOFEM, as well as test OOFEM’s nonlinear transient solving algorithm, a 1D model with arbitrary cure kinetics was created. In this model, the degree-of-cure was defined purely as a function of time and does not depend on the cure cycle or temperature results. To simulate actual behaviour of materials in an autoclave, a logistic curve was chosen to model the degree-of-cure. Parameters were adjusted such that the slope of the curve (which is proportional to the heat generated) would be at maximum near the beginning of the cure cycle hold phase. Using this setup, the model produced an exotherm around the same time as models with material-specific cure kinetics. The model was created in OOFEM and compared with finite difference models from three different programs. The details of these models are outlined below. Model Setup This model solves the 1D heat equation with convective boundary conditions: , , , , , , (C-1) (C-2) , (C-3) where: = specific heat capacity = 978 (J / (kg oC)) = material density = 1500 (kg / m3) = thermal conductivity = 0.12 (W / (m oC)) = convection coefficient = 50.0 (W / (m2 oC)) = length of bar = 0.01 (m) = temperature 83 = heat generation term = spatial coordinate = time For models in this appendix the heat generation term is evaluated as per Equation (C-4): , , (C-4) where: = heat of reaction = 4.74·105 (J / kg) = degree-of-cure = heat generation The curing function used in these models and its first derivative are: (C-5) (C-6) where: = degree-of-cure (unitless) = initial degree-of-cure = 0.01 (unitless) = cure rate = 0.0405 (minutes-1) = time (minutes) As this model was only concerned with comparing results up to and after the exotherm occurred, the autoclave cure cycle used consists of ramp-up and hold phases, but omitted a cool-down phase. The cure cycle used was: T t T T r·t t t t t (C-7) 84 where: = autoclave temperature = initial temperature = 20.0 (oC) = hold temperature = 180.0 (oC) = ramp rate = 1.0 (oC/min) Figure C - 1 shows degree-of-cure and heat generation plotted with the applied cure cycle. Note that the heat generation values have been normalized with respect to their peak value. 200 1 Cure Cycle 180 160 Temperature (C) 0.9 Normalized Q 0.8 DOC 140 0.7 120 0.6 100 0.5 80 0.4 60 0.3 40 0.2 20 0.1 0 0 0 50 100 150 200 250 300 Time (min) Figure C - 1 Temperature, degree-of-cure, and normalized heat generation for the logistic cure model Finite Difference Model – Matlab™ This model solves Equation (C-1) using the finite difference method. Derivatives were approximated using central difference equations in both time and space. Examples of central difference approximations for derivatives are shown in Equations (C-8) to (C-9). (C-8) (C-9) 85 where: = function approximating solution of the differential equation = independent variable (space or time) = step size in the independent variable = discrete position in the independent variable The model used takes derivatives at time corresponding to the Crank-Nicolson method. When finite difference approximations are substituted into Equation (C-1) the following equation is produced: (C-10) where: = temperature at discrete time step n and discrete space step j = heat generated at time step n and space step j = time step size = space step size The relationships in Equation (C-10) can be visualized using the Crank-Nicolson stencil shown in Figure C - 2. With some simple manipulation, this equation can be rearranged to solve for . In the Matlab™ model, this formula is then used to calculate values at all interior nodes for any given time step. Since the repeated application of Equation (C-10) creates multiple circular references, the cells are reevaluated until their answers converge. All that remains is to define Figure C - 2 CrankNicolson stencil the boundary conditions. The boundary conditions in Equations (C-2) and (C-3) pertain to the first spatial derivative of temperature, evaluated at x = 0 and x = L respectively. Again using finite difference approximations to evaluate the boundary conditions at time for x = 0 and x = L the following equations are derived: 86 (C-11) (C-12) where: = temperature on the left-boundary at time step n = temperature on the right-boundary at time step n , = temperatures at fictitious nodes at time step n = autoclave temperature at time step n = thermal conductivity = convection coefficient In order to estimate the first spatial derivative on the boundaries, two fictitious nodes must be introduced. Temperatures at these nodes do not correspond to any position in space, but are necessary to model the boundary conditions. With the presence of fictitious nodes, as well as the ambient autoclave temperature, Figure C - 3 shows the space discretization used to model this 1D problem. Figure C - 3 Finite difference geometry stencil for time step n In this figure, square nodes represent the autoclave temperature at time step n, their values are determined using Equation (C-7). The triangular nodes are the fictitious nodes, their values are determined by solving and evaluating Equations (C-11) and (C-12) the appropriate boundary. Finally, the circular nodes are the interior nodes. These nodes correspond to positions in the 1D domain. Their values are determined by solving and evaluating Equation (C-10). Note that for all nodes other than the autoclave nodes, results at each time step are determined by applying the appropriate equations iteratively until convergence is achieved. Instead of iterating, the 87 temperature values at each time step could be determined by solving a set of linear simultaneous equations. However, for the small time steps used, an iterative approach was faster and easier to implement. Finite Difference Model – Excel™, Maple™ This model was also implemented in Excel™ using the same equations outlined above. Finally, the differential equations and boundary conditions were programmed into Maple and solved using Maple’s numeric PDE solver (which uses a finite difference formulation). Results from all three of the finite difference models agree throughout the cure cycle within a maximum difference of 0.013 °C. Finite Element Method – OOFEM To test its nonstationary transient transport solving algorithm, the logistic cure model was also implemented in OOFEM. This model used a similar implementation as the CCA interface described in section 3.5.3. Key differences between the CCA interface and logistic cure model are threefold: 1. All relevant functions for the logistic cure model are hard-coded into OOFEM. No external library is called. 2. Material properties are constant for the logistic cure model. 3. Heat generated over any given time step follows Equation (C-13) , , (C-13) Results The finite difference models were implemented with a space step of 1 mm and time step of 1 second. Time steps of 0.5 and 0.25 seconds were also investigated, with little noticeable effect on the results. The OOFEM models used 8-node 3D elements with boundary conditions applied such that temperature variations only occur in one dimension. This model used the same time steps as the finite difference models. Again the change in time step size had little effect on the results. 88 Results from OOFEM agreed with those from the finite difference models within a maximum difference of 0.02 °C. Figure C - 4 shows results from the Excel™ finite difference model. Results from the other models show no noticeable difference when plotted together. 250 Temperature (C) 200 150 Autoclave Temp 100 Centre Temp 50 0 0 50 100 150 200 250 300 Time (min) Figure C - 4 Temperature histories for logistic cure model, Excel™ finite difference model These results verify the appropriateness of OOFEM’s nonstationary transient solution algorithm, as well as the basic implementation of curing materials in OOFEM. 89 Appendix D Surface Integration Using 3D Jacobian As of OOFEM release 1.8, the only 3D solid elements available in the transport module are 8node brick elements. These elements only support surface integration if all nodes on the surface are coplanar. This requirement is a product of the integration method employed in OOFEM. When integrating over a surface, OOFEM uses a 2D integration scheme similar to that employed in 2D element classes. This 2D integration method sets up new 2D coordinates for the nodes and integration points, as well as a 2D Jacobian. The entire surface integration then takes place over a 2D domain. With the addition of 24-node solid elements into OOFEM’s transport module, it was desirable to integrate over non-planar surfaces (eg. convective boundary surface of a curved part). For obvious reasons, the 2D surface integration scheme in OOFEM is not suitable for a curved surface in 3D space. Thus, to facilitate these surface integrations, a scheme was devised in which the 3D Jacobian for the element volume integrals is utilized for surface integration. Evaluating convective boundary conditions for a given time step involves the surface integral in Equation (D-1) below. (D-1) where: = boundary convection matrix = shape function vector = convective boundary surface = convection coefficient Equation (D-2) below shows Equation (D-1) using Gauss integration and expressed in natural coordinates. In this case the top surface of the element ( 1) is considered. 90 , , , (D-2) where: , , = integration point location in natural coordinates = integration point weight coefficient , = number of integration points in each direction In order to calculate Equation (D-2) in the proposed integration scheme, it is necessary to evaluate the term using each element’s 3D Jacobian matrix. Equation (D-3) below shows the components of the 3D Jacobian matrix for transforming from Cartesian to natural coordinates. (D-3) where: , , = Cartesian coordinates , , = natural coordinates Calculation of the term begins by defining a parametric vector, ̂ , , that gives the Cartesian coordinates of any point on the surface, as a function of the element natural coordinates on that surface. Once this vector is defined, (Stewart, 2003) states that the relationship between the differential areas are as defined in Equation (D-4). ̂ ̂ (D-4) 91 Since the element natural coordinate axes don’t necessarily correspond to the model global Cartesian coordinate axes, each global coordinate within an element will be a function of all three of that element’s natural coordinates (ie. , , , , ). However, since each integral is evaluated on a specific element surface, the global coordinates on the surface will be a function of only two natural coordinates. In the case considered here , , , . Thus the following equations relate differential distances in the relevant coordinates. ̂ , , , , (D-5) Applying partial differentiation to the vector, and adopting a shorthand notation, yields Equation (D-6). ̂ (D-6) ̂ Substituting Equation (D-6) into Equation (D-4) and evaluating yields Equation (D-7). (D-7) Substituting elements of the 3D Jacobian matrix into Equation (D-7) produces Equation (D-8) below. , , , , , , , , , , , , Equation (D-8) applies to two surfaces in a hexahedral element, specifically surfaces 1 Equation (D-9) applies, while Equation (D-10) applies for surfaces (D-8) 1. For the 1. 92 , , , , , , , , , , , , , , , , , , , , , , , , (D-9) (D-10) 93 Appendix E CCA Models and Input Parameters This appendix outlines the calculations performed inside the CCA library to determine various physical properties of Materials A and B. The models outlined here were taken from (Johnston, 1997). The input parameters defined were taken from the CCA configuration files provided with the CCA library. The following symbols are consistent throughout this appendix: α = degree-of-cure (unitless) T = temperature (°K or °C) t = time (s) Material A: Constant Parameters Although all of the material parameters in the CCA library are setup such that they can vary with temperature and degree-of-cure, some of the parameters associated with Material A are set to be constant. These constant parameters are outlined in Table E - 1. Table E - 1 Constant material parameters in the CCA library, Material A Parameter Value Units Resin Density Resin Specific Heat Capacity Resin Thermal Conductivity 1265 kg/m3 1260 J/(kg °K) 0.167 W/(m °K) Resin Poisson’s Ratio 0.37 none Fibre Density 1790 kg/m3 Fibre E11 2.10E+11 Pa Fibre E33 1.72E+10 Pa Fibre G13 2.76E+10 Pa Fibre ν13 0.2 none Fibre ν23 Composite Fibre Volume Fraction Composite Resin Volume Fraction 0.25 none 0.558 none 0.442 none 94 Default Linear Parameter Development The formula shown in Equation (E-1)Error! Reference source not found. is the generic ‘default’ model used to determine a number of different material parameters inside the CCA library. This model is linear with respect to both temperature and degree-of-cure. (E-1) where: = the parameter to be determined = nominal parameter value , , , = constants used to calibrate the model Several properties of Material A in the CCA library follow this default linear model. These parameters are: fibre specific heat capacity, fibre thermal conductivity and composite lumped coefficient of thermal expansion. The model constants associated with these properties are outlined in Tables E - 2 to E - 4. Table E - 2 CCA fibre specific heat capacity input parameters, Material A Parameter Value 750 0 Units J/(kg °K) °C 2.05 n/a n/a J/(kg °K2) none none Table E - 3 CCA fibre thermal conductivity input parameters, Material A Parameter Value 7.69 2.4 0 Units W/(m °K) W/(m °K) °C 1.56E‐02 W/(m °K2) 5.07E‐02 W/(m °K2) n/a none n/a none 95 Table E - 4 CCA composite lumped coefficient of thermal expansion input parameters, Material A Parameter Value 0 Units 1/°K 2.00E‐09 0 2.25E‐05 1/°K2 1/°K 1/°K 2.50E‐08 ‐3.50E‐06 20 0 1/°K2 1/°K °C none Cure Kinetics, CCA Model 1 Equation (E-2) gives the first time derivative of degree-of-cure in the cure kinetics model for Material A. 1 1 (E-2) ∆ Table E - 5 gives values for the constants present in Equation (E-2). Table E - 5 CCA cure kinetics input parameters, Material A Parameter Value Units 3.502E+07 1/s 8.070E+04 J/mol ‐3.357E+07 1/s 7.780E+04 J/mol 3.267E+03 1/s 5.660E+04 J/mol 0.47 none 0.30 8.314472 none J/(°K mol) ∆ ∆ ∆ 96 Resin Modulus Development, CCA Model 1 Equation (E-3) gives the formulae used in calculating the resin’s elastic modulus for Material A. 1 1 1 (E-3) Table E - 6 gives values for the constants present in Equation (E-3). Table E - 6 CCA resin modulus development input parameters, Material A Parameter Value 3.50E+06 3.50E+09 0 0.35 1.00 Units Pa Pa none none none 0 °C 0 1/°K Material B: Constant Parameters Similar to Material A, some of the parameters associated with Material B are set to be constant. These constant parameters are outlined in Table E - 7. 97 Table E - 7 Constant material parameters in the CCA library, Material B Parameter Value Units Resin Density 1300 kg/m3 Fibre Density 1790 kg/m3 Fibre E11 2.10E+11 Pa Fibre E33 1.72E+10 Pa Fibre G13 2.76E+10 Pa Fibre ν13 0.2 none Fibre ν23 Composite Fibre Volume Fraction Composite Resin Volume Fraction 0.25 none 0.573 none 0.427 none Default Linear Parameter Development Several properties of Material B in the CCA library follow the default linear model outlined in Equation (E-1). These parameters are: resin specific heat capacity, resin thermal conductivity, fibre specific heat capacity, fibre thermal conductivity and composite lumped coefficient of thermal expansion. The model constants associated with these properties are outlined in Tables E - 8 E - 8to E - 12.E - 12 Table E - 8 CCA resin specific heat input parameters, Material B Parameter Value 1005 20 Units J/(kg °K) C 3.74 0 0 J/(kg °K2) none none Table E - 9 CCA resin thermal conductivity input parameters, Material B Parameter Value 0.148 0 Units W/(m °K) °C 3.43E‐04 0 6.07E‐02 W/(m °K2) none none Table E - 10 CCA fibre specific heat capacity input parameters, Material B 98 Parameter Value 750 0 Units J/(kg °K) °C 2.05 n/a n/a J/(kg °K2) none none Table E - 11 CCA fibre thermal conductivity input parameters, Material B Parameter Value 7.69 2.4 0 Units W/(m °K) W/(m °K) C 1.56E‐02 W/(m °K2) 5.07E‐02 W/(m °K2) n/a none n/a none Table E - 12 CCA composite lumped coefficient of thermal expansion input parameters, Material B Parameter Value 6.00E‐07 Units 1/°K 0 0 2.86E‐05 1/°K2 1/°K 1/°K 1.00E‐07 0 20 0 1/°K2 1/°K °C none Cure Kinetics, CCA Model 6 Equation (E-4) gives the first time derivative of degree-of-cure in the cure kinetics model for Material B. (E-4) ∆ 99 Table E - 13 gives values for the constants present in Equation (E-4). Table E - 13 CCA cure kinetics input parameters, Material B Parameter ∆ C Value Units 1.53E+05 1/s 6.65E+04 J/mol 0.8129 none 2.736 43.09 none none 5.48E‐03 1/°K ‐1.684 8.314472 none J/(°K mol) Resin Modulus Development, CCA Model 2 Equation (E-5) gives the formulae used in calculating the resin’s elastic modulus for Material B. (E-5) 1 · · Table E - 14 gives values for the constants present in Equation (E-5). 100 Table E - 14 CCA resin modulus development input parameters, Material B Parameter Value 4.67E+06 4.67E+09 ‐45.7 0 Units Pa Pa °C 1/°K ‐12 °K 268 220 20 0 °K °K °C 1/°K Resin Poisson’s Ratio Development, CCA Model 2 Equation (E-6) gives the formula used to calculate the resin’s Poisson’s ratio for Material B. 1 1 2 1 2 ∞ (E-6) ∞ Table E - 15 gives values for the constants present in Equation (E-6). Table E - 15 CCA resin Poisson's ratio development input parameters, Material B Parameter , Value Units 0.37 none As calculated in Equation (E-5) 101
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Implementation and evaluation of a coupled thermal-structural...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Implementation and evaluation of a coupled thermal-structural analysis module for laminated composites… Beaton, Douglas Brian 2010
pdf
Page Metadata
Item Metadata
Title | Implementation and evaluation of a coupled thermal-structural analysis module for laminated composites in an open-source finite element code |
Creator |
Beaton, Douglas Brian |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | The aerospace industry invests heavily in the design and manufacture of composite materials. Aircraft components are produced by placing unprocessed composite materials in an autoclave and applying heat and pressure. The desired part geometry is achieved by forming raw composite materials around a tool, typically made of aluminum or other metal. Throughout the cure cycle, temperature changes cause the part and tool to expand at different rates. This differential expansion, combined with composite material properties that evolve over time, produces residual stresses in the part and leads to geometric instabilities (warpage) upon removal from the tool. Excessive warpage can render a part unusable. Errors of this nature can be quite costly, particularly in the aerospace industry where the tools created can be very large. A strong desire exists to predict the warpage and residual stresses imposed by the curing process and incorporate these stresses in the structural design of a component. To accomplish this goal for complex geometries, special additions to the finite element method are required. Commercial finite element programs provide some flexibility for users to implement custom elements and materials. Though, this flexibility has limits: some material models, such as non-local damage models, cannot be incorporated in proprietary software. This work selects an open-source finite element program and implements the ability to model curing processes of composite materials. The thermal and structural equations are solved in a coupled manner during each time step. This contrasts previous work by the UBC Composites Group, wherein the heat equation is solved over the entire model before the structural equations are considered. Numerous verification models are run to confirm the implementation, along with several example problems. Recommendations are made for further work to improve the process modeling and facilitate a link to subsequent structural models. Ultimately, the code produced represents the first step in seamlessly modeling composite structures during manufacturing processes through to in-service conditions. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-13 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
IsShownAt | 10.14288/1.0062873 |
URI | http://hdl.handle.net/2429/29125 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_fall_beaton_douglas.pdf [ 1.89MB ]
- Metadata
- JSON: 24-1.0062873.json
- JSON-LD: 24-1.0062873-ld.json
- RDF/XML (Pretty): 24-1.0062873-rdf.xml
- RDF/JSON: 24-1.0062873-rdf.json
- Turtle: 24-1.0062873-turtle.txt
- N-Triples: 24-1.0062873-rdf-ntriples.txt
- Original Record: 24-1.0062873-source.json
- Full Text
- 24-1.0062873-fulltext.txt
- Citation
- 24-1.0062873.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0062873/manifest