University of Southern California Dissertations and Theses

Mathematical model and numerical analysis of a shape memory polymer composite cantilever beam for space applications

(USC Thesis Other)

Mathematical model and numerical analysis of a shape memory polymer composite cantilever beam for space applications

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MATHEMATICAL MODEL AND NUMERICAL ANALYSIS OF A SHAPE MEMORY POLYMER COMPOSITE CANTILEVER BEAM FOR SPACE APPLICATIONS by Dean Bergman A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (AEROSPACE ENGINEERING) May 2014 Copyright 2014 Dean Bergman Table of Contents Nomenclature v List Of Figures viii List Of Tables x Abstract xi Chapter 1: Introduction 1 Chapter 2: Mathematical Background 9 2.1 Hamilton's Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 Variational Principle . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.2 Useful Manipulations . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Newton's Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 Single Variable Newton's Method . . . . . . . . . . . . . . . . . . . 12 2.2.2 Multi-variable Newton's Method . . . . . . . . . . . . . . . . . . . 12 2.3 Fourth-Order Runge-Kutta Method . . . . . . . . . . . . . . . . . . . . . 13 Chapter 3: Analytical Beam Model 15 3.1 General Energy Functionals . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.1 Strain Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1.2 Kinetic Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1.3 Virtual Work by Nonconservative Forces . . . . . . . . . . . . . . . 19 3.2 Linear Beam Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2.1 Timoshenko Linear Beam Theory . . . . . . . . . . . . . . . . . . . 21 3.2.2 Euler-Bernoulli Beam Theory . . . . . . . . . . . . . . . . . . . . . 22 3.2.3 Rayleigh Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3 Nonlinear Beam Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.1 Small Shear and Extension . . . . . . . . . . . . . . . . . . . . . . 24 3.3.2 Pure Bending as a Function of Vertical De ection (w) . . . . . . . 26 3.3.3 Pure Bending as a Function of Local Rotation () . . . . . . . . . 30 ii Chapter 4: Shape Memory Polymer Composite Beam Material Model 34 4.1 1D Stress-Strain Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 2D Composite Material . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2.1 Unidirectional Laminae . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2.2 Equivalent Modulus of Elasticity for an SMP Composite Laminate 42 4.2.3 Other Types of Laminae . . . . . . . . . . . . . . . . . . . . . . . . 44 4.2.4 Constitutive Equation for a Beam . . . . . . . . . . . . . . . . . . 45 Chapter 5: Shape Memory Polymer Composite Cantilever Beam 49 5.1 Quasi-static Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.1.1 Quasi-Static Equilibrium Described by the Local Vertical De ection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.1.2 Quasi-Static Equilibrium Described by the Local Rotation Coordinate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.2 Temperature Dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2.1 Piecewise Constant Temperature Dependence for a Pure SMP Beam 54 5.2.2 Full Temperature Dependant Model . . . . . . . . . . . . . . . . . 56 5.3 Shape Fixation and Shape Recovery . . . . . . . . . . . . . . . . . . . . . 57 5.3.1 Cooling for Shape Fixation . . . . . . . . . . . . . . . . . . . . . . 57 5.3.2 Heating for Shape Recovery . . . . . . . . . . . . . . . . . . . . . . 60 5.3.3 Vibration after Shape Recovery . . . . . . . . . . . . . . . . . . . . 61 Chapter 6: Discretization and Numerical Methods 63 6.1 Finite Dierence Discretization in Terms of Vertical De ection . . . . . . 64 6.2 Finite Dierence Discretization in Terms of Rotation . . . . . . . . . . . . 67 6.3 Comparison of Methods Used to Model the Quasi- Static Nonlinear Cantilever Beam . . . . . . . . . . . . . . . . . . . . . . . 68 6.3.1 Coordinate Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 68 6.3.2 Solution Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.3.3 Calculation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.3.4 Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Chapter 7: Finite Element Derivations 76 7.1 Continuum Finite Element Equations . . . . . . . . . . . . . . . . . . . . 76 7.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 7.1.2 Updated Lagrangian Method . . . . . . . . . . . . . . . . . . . . . 79 7.1.3 Total Lagrangian Method . . . . . . . . . . . . . . . . . . . . . . . 81 7.1.4 Numerical Quadrature . . . . . . . . . . . . . . . . . . . . . . . . . 83 7.2 Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 7.2.1 Continuum-Based (CB) Beam . . . . . . . . . . . . . . . . . . . . . 83 7.2.1.1 Element Derivation . . . . . . . . . . . . . . . . . . . . . 84 7.2.1.2 Updated Lagrangian Beam Formulation . . . . . . . . . . 86 7.2.1.3 Total Lagrangian Beam Formulation . . . . . . . . . . . . 89 7.2.1.4 3D CB Beam Equations . . . . . . . . . . . . . . . . . . . 90 7.2.2 Shell Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7.2.2.1 Element Derivation . . . . . . . . . . . . . . . . . . . . . 91 iii 7.2.2.2 Updated Lagrangian Shell Formulation . . . . . . . . . . 95 7.2.2.3 Total Lagrangian Shell Formulation . . . . . . . . . . . . 97 7.2.2.4 Reduced Integration . . . . . . . . . . . . . . . . . . . . . 98 Chapter 8: Constitutive Model for Finite Element Analysis 100 8.1 Shape Memory Polymer Constitutive Model . . . . . . . . . . . . . . . . . 100 8.2 Composite Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 8.2.1 Unidirectional Laminae . . . . . . . . . . . . . . . . . . . . . . . . 103 8.2.2 Multilayer Laminate . . . . . . . . . . . . . . . . . . . . . . . . . . 110 Chapter 9: Finite Element Solution Methods 112 9.1 Direct Solution from Force Vector . . . . . . . . . . . . . . . . . . . . . . 112 9.2 Linearization of FE Equations . . . . . . . . . . . . . . . . . . . . . . . . . 114 9.2.1 Element Stiness Matrices . . . . . . . . . . . . . . . . . . . . . . . 115 9.2.2 Material Tangent Stiness with Stored Strain . . . . . . . . . . . . 119 9.2.3 Element Stiness Matrix for Shells . . . . . . . . . . . . . . . . . . 120 9.2.4 Global Stiness Matrix with Multilayering . . . . . . . . . . . . . . 122 9.3 Dynamic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 9.3.1 Explicit Dynamic Solvers . . . . . . . . . . . . . . . . . . . . . . . 124 9.3.2 Implicit Dynamic Solvers . . . . . . . . . . . . . . . . . . . . . . . 125 9.3.3 Mass Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Chapter 10: Shape Fixation and Shape Recovery for FE Model 127 10.1 Cooling for Shape Fixation . . . . . . . . . . . . . . . . . . . . . . . . . . 127 10.2 Heating for Shape Recovery . . . . . . . . . . . . . . . . . . . . . . . . . . 129 10.3 Vibration after Shape Recovery . . . . . . . . . . . . . . . . . . . . . . . . 130 Chapter 11: Simulation Results 131 11.1 Loading at T h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 11.2 Cooling from T h to T l with load applied . . . . . . . . . . . . . . . . . . . 134 11.3 Load removal at T l . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 11.4 Shape Recovery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 Chapter 12: Conclusions 142 Chapter 13: Recommendations for Future Work 145 References 148 iv Nomenclature D Equivalent Modulus of Elasticity for a Beam Normal Stress Tensor " Normal Strain Tensor " s Stored Strain Tensor " T Coecient of Thermal Expansion Tensor C Material Tensor E e Modulus of Entropic Deformation Tensor E i Modulus of Internal Energetic Deformation Tensor S e Inverse of Modulus of Entropic Deformation Tensor S i Inverse of Modulus of Internal Energetic Deformation Tensor Vriation Fiber Orientation Angle Shear Strain ^ s A Number Related to Timoshenko Beam Theory ^ y Perpendicular Coordinate to the Mid-surface Curvature T Kinetic Energy U Potential Energy W nc Work Done by Non-conservative Forces f Force vector u Unknown displacement vector Poisson's Ratio Density Shear Stress Angle With Respect to the Horizontal A Cross-sectional Area b Beam Width C Total Modulus of Elasticity c f Frozen Fraction Constant D Equivalent Modulus of Elasticity E Total Modulus of Elasticity E h Modulus of Elasticity at High Temperature E l Modulus of Elasticity at Low Temperature e b Distance from Mid-surface to Bending Surface E e Modulus of Entropic Deformation E i Modulus of Internal Energetic Deformation v f General Body Load G Shear Modulus H Internal Horizontal Shear Force I Second Moment of Area k Boltzmann's Constant L Total Length M Moment M e General Eective Applied Moment m z Applied Moment along the Body N Cross-Link Density N Force n Frozen Fraction Constant P x Applied Horizontal Point Load P y Applied Vertical Point Load Q e General Eective Applied Point Load q x Applied Horizontal Distributed Load q y Applied Vertical Distributed Load s Arc Length Coordinate T Temperature t Time T g Glass Transition Temperature T h High Temperature u Horizontal Displacement u 0 Horizontal Displacement of the Mid-surface V Internal Vertical Shear Force V f Fiber Volume Fraction V m Matrix Volume Fraction w Vertical De ection w 0 Vertical Displacement of the Mid-surface x Global Horizontal Coordinate y Global Vertical Coordinate z Out of Plane Coordinate Coecient of Thermal Expansion f Frozen Fraction Normal Stress " Normal Strain " s Stored Strain " g Stored Strain at T g N Number of Layers in a Laminate Subscripts ^ y Perpendicular to Arc Length Direction f Fiber m Matrix s Arc Length Direction vi 1 Principle Direction Along Fiber Orientation 2 Principle Direction Perpendicular to Fiber Orientation Superscripts 0 Extension Curvature int Internal s stored T Thermal ext External vii List Of Figures 1.1 (a) JHU/APL Hybrid In atable Antenna (b) JPL SMP Composite Truss 6 1.2 SMP forming and recovery process . . . . . . . . . . . . . . . . . . . . . . 7 3.1 General beam de ection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Cantilever beam with general loading conditions . . . . . . . . . . . . . . 19 3.3 Beam element under pure bending . . . . . . . . . . . . . . . . . . . . . . 26 4.1 Young's Modulus as a function of temperature . . . . . . . . . . . . . . . 36 4.2 (a) Composite material loaded in the 1 principle direction, (b) Composite material loaded in the 2 principle direction . . . . . . . . . . . . . . . . . 37 4.3 Coordinates of a Lamina . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.4 Laminate coordinates and numbering . . . . . . . . . . . . . . . . . . . . 41 5.1 Large deformation of an SMP cantilever beam . . . . . . . . . . . . . . . . 49 5.2 (a) Free-body diagram for an increment of a cantilever beam, (b) Cantilever Beam with applied body forces and point loads . . . . . . . . . . . . . . . 50 6.1 A cantilever SMPC beam with general loading conditions . . . . . . . . . 64 6.2 Cantilever beam discretized by a nite dierence approximation in local coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 7.1 A 2 node CB Beam element . . . . . . . . . . . . . . . . . . . . . . . . . . 84 7.2 A beam element with corotational unit vectors ^ e x and ^ e y , and director vectors p () . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 viii 7.3 An 8 node shell element . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 7.4 (a) Rotation described by a vector, (b) Rotation described by angles . . . 92 8.1 (a) Composite material loaded in the 1 principle direction, (b) Composite material loaded in the 2 principle direction . . . . . . . . . . . . . . . . . 103 8.2 Coordinates of a Lamina . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 8.3 Layers in the laminate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 11.1 Large de ection cantilever beam with a point load applied to the end . . . 131 11.2 Plot of normalized load versus displacement . . . . . . . . . . . . . . . . . 133 11.3 Plot of increasing normalized de ection of the beam as load increases . . . 134 11.4 Plot of normalized de ection of a pure SMP beam along its length comparing the cases of when the load is applied and when it is removed . . . . . . . 135 11.5 Plot of normalized de ection of SMP composite beam along its length comparing the cases of when the load is applied and when it is removed with 1 layer, 0 layup angle and varying ber volume fraction . . . . . . 136 11.6 Plot of normalized de ection of SMP composite beam along its length comparing the cases of when the load is applied and when it is removed with 1 layer (2 for [45 45 ] ), 10% ber volume fraction and varying layup angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 11.7 Plot of normalized de ection of SMP composite beam along its length comparing the cases of when the load is applied and when it is removed with 8 layers, 60% ber volume fraction and [90 90 90 0 ] s layup . . . . . 139 11.8 Plot of normalized de ection of SMP composite beam along its length comparing the cases of when the load is applied and when it is removed with 8 layers, 50% ber volume fraction and [90 90 90 0 ] s layup . . . . . 140 ix List Of Tables 6.1 Advantages and disadvantages of 2nd order global coordinates . . . . . . . 69 6.2 Advantages and disadvantages of 4th order global coordinates . . . . . . . 69 6.3 Advantages and disadvantages of 4th order local coordinates . . . . . . . . 70 6.4 Advantages and disadvantages of curvilinear coordinates . . . . . . . . . . 70 6.5 Advantages and disadvantages of using elliptic integrals . . . . . . . . . . 71 6.6 Advantages and disadvantages of using state-space and Runge-Kutte . . . 72 6.7 Advantages and disadvantages of using nite dierence method . . . . . . 72 6.8 Advantages and disadvantages of using MATLAB functions . . . . . . . . 73 6.9 Advantages and disadvantages of using single variable Newton's method . 73 6.10 Advantages and disadvantages of using multi-variable Newton's method . 74 6.11 Advantages and disadvantages of using continuation methods . . . . . . . 74 7.1 Reduced Integration in Shells . . . . . . . . . . . . . . . . . . . . . . . . . 99 11.1 Values of the beam parameters used in the simulation . . . . . . . . . . . 132 11.2 Values used to calculate f and E e . . . . . . . . . . . . . . . . . . . . . . 132 x Abstract Shape Memory Polymer Composite (SMPC) structures are an ideal solution for large space structure deployment problems. This is due to their ability to be formed into a desired light-weight loading shape and then transform back to their original shape by means of an applied stimulus. There is a wide array of constitutive and qualitative work being done on SMPC's but little or no development of dynamic equations. This dissertation documents a macroscopic model for the shape xation and shape recovery processes of a SMPC cantilever beam. In particular the focus is on a quasi-static model that can be used to describe the shape xation process. Numerical results are obtained in this regard by use of nite dierence approximation with Newton's method as well as a Finite Element Formulation. These models combine a nonlinear geometric model with a temperature dependent constitutive law. Additionally, the dynamic equations of the SMPC cantilever are derived. The Finite Element models are derived for a beam and a shell in order to keep the solution as general as possible. The shell element allows this to be expanded to thin or moderately thick structures of any shape. This study is a precursor to on-going work in modeling SMPC space structures. xi Chapter 1 Introduction This dissertation describes the nonlinear modeling and analysis of a cantilever beam made of Shape Memory Polymer Composite (SMPC) material. The intent is to eventually develop modeling and design tools for deployable structures in space applications. The current investigation is an extension of the previous work done by Yang [61] which studied the dynamics of a one-dimensional (transverse motion) Shape Memory Polymer (SMP) cantilever beam using linear geometry with a simplied material model. The current paper addresses the analysis of the SMP beam by considering nonlinear geometry, adding in a composite material and using the full temperature dependent material model. The beam model will be described in local coordinates, where the term \local coordinates" refers to the material, Lagrangian or undeformed coordinate system, whereby the beam is measured along its arc length. To be more focused on the nonlinearity, a quasi-static model will be considered, which neglects the inertia term in the beam governing equation. Quasi-static models are common in nonlinear structural mechanics. Here \quasi-static" means the beam de ection is temperature dependent. 1 In order to model the SMPC beam, a beam model needs to be chosen. An in depth study of beam models is presented in the report. The analysis of beams has been studied for over 200 years. A beam model can be used when the cross-sectional height and width are small with respect to the length of the element. Beams are a very useful modeling tool as the equations are much simpler to derive and solve than higher order element, such as plates or shells. This makes them ideal for obtaining an initial understanding of more complicated structures. A large amount of work has been done on dynamic linear beam models, where the de ections and strains are small with respect to the length. Euler and Bernoulli presented a model that only accounts for the vertical de ection [55]. Rayleigh then extended this to correct for the rotational inertia through de ection [51]. Timoshenko came up with a model for thicker beams whereby the shearing eect could not be ignored. In order to do this, however, the beam had to be described by the de ection and rotation of the beam [54]. An excellent review on linear beam theory is presented by Han et al [24]. There has also been a wide range of work done on nonlinear beams, which ranges back to the work of Euler in 1744 [19]. Nonlinear beams are needed when the de ection is large compared to the cross-sectional area of the beam causing the linear models to no longer be valid [20]. The equations needed to model a nonlinear beam can be obtained by elastica theory [28] and are often very complicated, such as Pai's models whereby all three de ection variables [43] or Euler angles [42] are used to describe the beam. As such it is common to implement a nite element solver to analyse these [8]. Beam models, however, are often intended to simplify the structural problem and as such a FEM model would be 2 overkill. It would be more benecial to use an analytical beam model similar to that of the linear beams. It is common in structural mechanics to model structures as static or quasi-static due to the complexity of the problem. Various analytical models have been presented to describe the beam using two or even all three of the beam coordinates [57]. However, it is convenient to simplify the problem further. This is done by assuming that the beam is in pure bending so that the model can be described in one coordinate rather than three. Bisshopp and Drucker presented a method to solve for large de ection with a point load applied using elliptic integrals [10]. This method can be extended to solve beams subjected to other simple loads, such as a moment or body force. It cannot however be used for combined loading. There are two other major problems with the method. Firstly, the method uses the moment equation to solve the problem. This means that it cannot be directly extended to a dynamic case. The other problem is that elliptic integrals are dicult to solve and generally need to be solved numerically anyway [38]. Holden [25] uses a similar approach to solve cantilever beam and buckling problems but he uses a Runge-Kutte solver to solve the equations instead of elliptic integrals. Batista et al. [6] uses elastica theory to model a static beam for an inclined end point load and as such nds the equilibrium positions. A fourth order static beam equation can be derived which describes the beam in terms of the vertical de ection measured along the horizontal (x) direction [20]. This model is complicated by the fact that end-shortening needs to be accounted for due to the large de ection. It also becomes a problem if the beam bends back on itself whereby there would be two de ections for one x coordinate. A better model, albeit more complicated, was presented by AL-Saddar whereby he described the 3 beam by only the vertical de ection along the arc length (s) [2]. The shortcoming with this model, though, is that there is a singularity in the moment equation. In order to overcome some of the shortcomings of the methods presented above, a similar model to the AL-Sadder beam is derived in this paper. This model describes the beam by the rotation along its arc length rather than the vertical de ection. This presents a more robust model, as it has no singularity, and it is easier to implement. The few shortcomings inherent in the model will be highlighted. Numerical simulations for this method will also be presented using the nite dierence method with the use of Newton's Method solver. Some of the static models stated above will also be extended to include the dynamic behavior of the beam. While a beam model is useful for gaining an initial understanding of the structural behavior, it is necessary to develop a Finite Element (FE) model. This is so that more general structures can be modeled. There has been much work done on developing beam, plate and shell FE models. A thorough study has been done in order to select the most appropriate model for the SMPC structure. FE models can be split into two groups: structural and continuum based models. A structural model takes the analytical energy equations for beams and plates, after all assumptions have been made, and converts these to a FE equations. Continuum based models, on the other hand, use the FE model based on continuum mechanics and then makes assumption after the fact to obtain a beam, plate or shell. The most commonly used structural beam models are the Euler-Bernoulli and Tim- oshenko beams. These use the same energy equations from the beam theories but they are converted into FE equations [18]. However, for application to an SMPC, it may be 4 more benecial to use a continuum based (CB) beam element. This theory was developed by Belytschko, whereby he took shell theory along with the nite element equations in continuum mechanics and applied it to beams [8]. This means that fewer assumptions are made compared to the Euler-Bernoulli or Timoshenko theories, creating a more general beam element. Belytchko uses a corotational coordinate system in order to utilize the Updated Lagrangian FE method. It was found, however, that this is unnecessary if the Total Lagrangian method is used instead. This paper also extends the CB model to 3D in order to be able to model composites. In modeling composites structures, it is important to recover the 2D planar and shearing eects. To do this it is more useful to use plates or shells rather than beam element. This will also help with modeling more complicated shaped structures. Plates and shells are very similar, the main dierence being that shells allow for membrane stresses [33]. Similar to the beams, there are various structural models for these, however, in order to keep this study as general as possible the focus will only be on continuum based shell elements. CB shell elements where rst introduced by Ahmad, Irons and Zienkiewicz [1]. This model was then converted to nonlinear form by Parisch [44] and then by Hughs and Liu [26,27]. The Parisch model uses rotations to describe the cross-section while the Hughes model uses vectors. Both methods will be highlighted here and then the more appropriate one will be selected. Zienkiewicz also investigated the reduced integration needed due to the inherent membrane and shear locking in shell elements [63]. Triangular shell elements have been shown to be more robust than their rectangular counterparts [11] but this was deemed to be beyond the scope of this research. 5 The nal stage in developing the geometric model is to implement the composite material. Barbero [3,64] and Reddy [31,32,40] with their respective collaborators are the most noted for this work. Each presents detailed studies of linear and nonlinear composite shell analysis. A similar approach will be detailed herein to model the SMPC structure. (a) (b) Figure 1.1: (a) JHU/APL Hybrid In atable Antenna (b) JPL SMP Composite Truss SMP materials have been discovered to be useful for many applications, such as heat-shrinkable tubes for electronics or lms for packaging [15], self-deployable sun sails in spacecraft [13], intelligent medical devices [56], implants for minimally invasive surgery [30], among various others [7]. The utility of SMP's in these applications is due to the ability to change from an \initial state" to a \nal state" after a stimulus is applied. This characteristic makes them useful for deployment of structures on a spacecraft such as large re ectors. SMP structure can be stowed in very compact shapes and then deployed into large structures. This has already been proven by JHU/APL with the Hybrid In atable Antenna (Fig. 1.1(a)) [58] and JPL with the Advanced Precipitation Radar Antenna Singly Curved Parabolic Antenna Re ector (Fig.1.1(b)) [34]. Another use for SMP's in 6 space structures is to use them for hinges. The advantage of SMP hinges is that they are light weight and do not require heavy mechanical devices, such as motors, to work [29]. An SMP structure can be formed into its \initial state" by forming the structure into the desired shape through a process of working, heating and cooling, known as programming. From this state the aforementioned stimulus, such as heat, can be applied to bring the structure back into its \nal state" [7]. This process is shown in Fig. 1.2, whereT g is the glass transition temperature of the polymer. The T g is the temperature at which the polymer changes between high and low stiness [35]. This is used to control the change from the polymer's \initial state" to its \nal state". T g can be controlled in order to design the polymer for the desired application. Heating Heating Cooling Shape Memory Deformation Shape Fixation Shape Recovery Figure 1.2: SMP forming and recovery process The SMP itself has a lower strength and stiness than would be required for space applications. However, by including composites in the SMP, the rigidity of the structure can be greatly improved. [47] A wide range of work has been done in this regard through the use dierent types of composites, such as: fabric weave composites [62]; chopped strand composites [41]; carbon nanotube composites [39]; and nanocomposites [36]. Most of the work done has been to analyse the recovery characteristics of the SMPC's. This study 7 will incorporate unidirectional composites into the SMP material model and mention will be made of how this can be extended to other types of composites, such as a plain weave. There has been a large amount of research done on SMP's and their composites, all of which relates to the material itself and its behavior under experimental conditions. Little to no work has been done on the development of dynamic equations. For this reason, the scope of this research includes deriving a dynamic theoretical geometric model and developing numerical techniques to simulate the performance of an SMPC beam. The scope does not include developing a material model for the SMP itself, as there is enough research being done on this topic. Thus a previously derived model can be used. The SMP material model, however, will need to be extended to account for a composite material. The importance of this research is to develop governing equations in order for the material to be properly modeled. To do this, a new composite material model is derived from an existing SMP material model. This is then used in a robust nonlinear beam and Finite Element model to simulate the behavior under dierent loading conditions. It is expected that this numerical model will be able to be used with any SMP material model to simulate the behavior of the SMPC structure through shape xation and recovery. The intent is that researchers and engineers may enjoy the benet of using this model to simulate a wide variety of SMPC structures for space application. 8 Chapter 2 Mathematical Background 2.1 Hamilton's Principle In order to nd the equations of motion for continuous structures in mechanics, it is useful to use Hamilton's principle. Hamilton's principle states that the energy in the system should be minimised for all time, namely [23] Z t 2 t 1 [TU +W nc ]dt = 0 (2.1) where, T is the kinetic energy functional, U is the potential energy functional, W nc is the energy functional due to non-conservative forces, and denotes the variation of the functional. 9 2.1.1 Variational Principle The variational operator () acts very much in the same way as a dierential operator (d), except that the variational does not dierentiate with respect to position [48], i.e. dF u 0 ;u;s = @F @s ds + @F @u du + @F @u 0 du 0 (2.2a) F u 0 ;u;s = @F @u u + @F @u 0 u 0 (2.2b) Thus some useful relations for variationals can be found, with F 1 =F 1 (u), F 2 =F 2 (u), G =G (u;v), and u and v are the dependant variables: (F 1 F 2 ) =F 1 F 2 (2.3a) (F 1 F 2 ) =F 1 F 2 +F 1 F 2 (2.3b) F 1 F 2 = (F 1 F 2 F 1 F 2 )=F 2 2 (2.3c) (F 1 ) n =n (F 1 ) n1 F 1 (2.3d) G = @G @u u + @G @v v (2.3e) It is also important to note that the variational operator is interchangeable with the dierential or integral operators, i.e.: du dx = d dx (u) (2.4a) Z udx = Z udx (2.4b) 10 2.1.2 Useful Manipulations It is often necessary to integrate by parts when nding the equations of motion in order to reduce the order of the dependant variable, for instance: I (u;v) = Z u dv dx dx =uv Z du dx vdx (2.5) It may also happen that the variational operator is inside another integral. If this is the case, then integration by parts is done to increase the order of the dependant variable, i.e.: I (u;v) = Z v Z ! udx dx = Z ! vdx Z ! udx Z Z ! vdx udx = Z Z vdx Z ! vdx udx (2.6) where ! is a subspace of . 2.2 Newton's Method There are many ways to solve single and systems of equations. The most common way is to use Newton's method. This can be done for a single equation or a system of equations as highlighted below [14]. 11 2.2.1 Single Variable Newton's Method For a given equation described by a single variable x f (x) = 0 (2.7) one can use Newton's method to solve for x. This is done by: x i+1 =x i f (x i ) f 0 (x i ) (2.8) To start the method, take a reasonable initial guess x 0 and substitute it into the equation above. The solution from this can then be taken and substituted back in, and this can be done until the solution converges to some desired tolerance. A `reasonable guess' is one such that the solution does not diverge. 2.2.2 Multi-variable Newton's Method The idea from the single variable problem can be extended to a system of n equations with n unknowns such that F (x) = 0 (2.9) where F = f 1 (x) f 2 (x) ::: f n (x) T (2.10a) x = x 1 x 2 ::: x n T (2.10b) 12 This can be solved by x i+1 =x i F (x i )J 1 (x i ) (2.11) where J is the Jacobian matrix of F such that J = 2 6 6 6 6 6 6 6 6 6 6 4 @f 1 @x 1 @f 1 @x 2 ::: @f 1 @xn @f 2 @x 1 @f 2 @x 2 ::: @f 2 @xn . . . . . . . . . . . . @fn @x 1 @fn @x 2 ::: @fn @xn 3 7 7 7 7 7 7 7 7 7 7 5 (2.12) This can then be solved in the same way as the single variable method, whereby a reasonable initial guess is taken and the method is iterated until convergence is reached. 2.3 Fourth-Order Runge-Kutta Method Initial value dierential equations can be solved in many ways. The most popular numerical method being the classical fourth-order Runge-Kutta method. This is performed as follows [14]. For some function dy dx =f (x;y) (2.13) This can be solved incrementally, starting at known values x 0 and y 0 , by y i+1 =y i + 1 6 ( 1 + 2 2 + 2 3 + 4 )h (2.14) 13 where h is the incrementation length and 1 =f (x i ;y i ) (2.15a) 2 =f x i + 1 2 h;y i + 1 2 1 h (2.15b) 3 =f x i + 1 2 h;y i + 1 2 2 h (2.15c) 4 =f (x i +h;y i + 3 h) (2.15d) This can be performed over the range of x values. This method is ecient and easy to implement. 14 Chapter 3 Analytical Beam Model Before the SMP structure can be analysed, a geometric model is needed. For simplicity sake, a cantilever beam will be used. Various linear and nonlinear models are given below. 3.1 General Energy Functionals To start, the total horizontal and vertical displacements, u andw respectively, need to be found. This can be done by analyzing Fig. 3.1. y x s u 0 u w 0 w θ ŷ ŷ Figure 3.1: General beam de ection Here, is the angle of a perpendicular to the cross-sectional face of the beam with respect to the horizontal at any given point along its length; the subscript 0 refers to the 15 de ections at the neutral surface; s and ^ y are the local (Lagrangian) coordinates, with s along the length of the beam and ^ y the perpendicular measure from the neutral surface; coordinates x andy are the horizontal and vertical global (Euler) coordinates respectively. From this, it is seen that u =u 0 ^ y sin (3.1a) w =w 0 ^ y + ^ y cos (3.1b) In order to nd the strains in the beam, the derivatives of the de ections are needed. Local coordinates will be used, as this makes it easier track a point along the beam. @u @s = du 0 ds ^ y d ds cos , @u @^ y = sin (3.2a) @w @s = dw 0 ds ^ y d ds sin , @w @^ y = cos 1 (3.2b) The strain components, e , are of the form [37] e = 1 2 (u ; +u ; +u ; u ; ) (3.3) 16 where greek indices run from direction s to ^ y and repeated indices are summed. Thus expanding and substituting (3.1), e s = du 0 ds + 1 2 " du 0 ds 2 + dw 0 ds 2 # ^ y d ds 1 + du 0 ds cos + dw 0 ds sin + 1 2 ^ y d ds 2 (3.4a) e ^ y = 0 (3.4b) s^ y = 1 2 dw 0 ds cos 1 + du 0 ds sin (3.4c) s^ y denotes the shear strain. In most cases, the ^ y 2 term in Eq. (3.4a) is ignored, but it will be left in for now in order to keep the equations as general as possible. For the sake of focusing on the beam equations, a homogeneous isotropic linearly elastic material will be used, such that the stresses, , and shear stress, s^ y , are [37] s =Ee s (3.5a) ^ y =Ee ^ y = 0 (3.5b) s^ y =G s^ y (3.5c) where E is the modulus of elasticity, and G is the shear modulus. A simplication in notation will be used in Eq.'s (3.4), such that u =u 0 , w =w 0 , and ( ) 0 =d( )=ds. This can be done since the total de ections are no longer needed. Hamilton's principle [23] will be used to nd the governing equations, namely Z t 2 t 1 [TU +W nc ]dt = 0 (3.6) 17 However, before this can be done, the variation of the energy functionals are needed. 3.1.1 Strain Energy The strain energy density, , is needed. This is found from [57] = Z (es;e s^ y ) (0;0) de =Ee 2 s +Ge 2 s^ y (3.7) The strain energy of the beam is U = 1 2 ZZZ V dV = 1 2 Z L 0 EA u 0 + 1 2 u 02 +w 02 2 ds + 1 2 Z L 0 EI 02 1 +u 0 cos +w 0 sin 2 + 1 2 u 0 + 1 2 u 02 +w 02 ds + 1 2 ZZZ V E 1 2 ^ y 0 4 dV + 1 2 Z L 0 ^ sGA w 0 cos 1 +u 0 sin 2 ds (3.8) where ^ s is the Timoshenko shear coecient which is dependent on the cross-sectional area of the beam [9], A is the cross-sectional area, and I is the second moment of area. 3.1.2 Kinetic Energy The kinetic energy in the beam is based on the de ection and rotation of the middle surface. Thus, with the notation simplication that _ ( ) =d( )=dt, the kinetic energy is given by T = 1 2 Z L 0 A _ w 2 + _ u 2 ds + 1 2 Z L 0 I _ 2 ds (3.9) 18 where is the linear density of the beam. The term I relates to the moment of inertia and can be substituted for J (the moment of inertia) or k 2 m as needed, where k is the radius of gyration and m is the mass of the beam [9]. Taking the variation T = Z L 0 A ( _ w _ w + _ u _ u)ds + Z L 0 I _ _ ds (3.10) Since Hamilton's principle will be used, it is convenient to take the time integral of this. Thus by integrating by parts with respect to time and noting that the variations at the beginning and ending time are equal to zero, the following is obtained: Z t 2 t 1 Tdt = Z t 2 t 1 Z L 0 A ( ww + uu)ds + Z L 0 I ds dt (3.11) 3.1.3 Virtual Work by Nonconservative Forces An example of a cantilever beam with general loading conditions is shown in Fig. 3.2. The beam has body forces in the horizontal (q x ) and vertical direction (q y ), body moment (m z ), end point loads in the horizontal (P x ) and vertical (P y ) direction, and end moment (M e ) applied. s y x M e P x P y q x m z q y Figure 3.2: Cantilever beam with general loading conditions 19 In this case, the energy from nonconservative forces will be W nc = Z L 0 (q x u +q y w +m z )ds +P x u s=L +P y w s=L +M e s=L (3.12) More generally, the point loads and moment at s =L can also be at s = 0 depending on the boundary conditions. The variation of this case is W nc = Z L 0 (q x u +q y w +m z )ds +P x u s=L +P y w s=L +M e s=L (3.13) A point load along the length of the beam can also be added here but it is easier to think of it as a body force, where for example q y =P y (ss p ) with (s) being the dirac delta function and s p the point of application. 3.2 Linear Beam Models The proposed nonlinear beam model can be reduced to the various existing linear beam models. In linear beam theory it is assumed that is small (small angle assumption) and that all nonlinear terms are negligable. This means that dxds. Thus, the variation of the strain energy functional becomes: U = Z L 0 EAu 0 u 0 +EI 0 0 + ^ sGA w 0 w 0 ds (3.14) 20 Integrating by parts: U =EAu 0 u L 0 +EI 0 L 0 + ^ sGA w 0 w L 0 Z L 0 EAu 00 u +EI 00 + ^ sGA w 00 0 w ^ sGA w 0 ds (3.15) There are two main theories for linear beams, namely the Timoshenko beam and the Euler-Bernoulli beam. The Rayleigh beam is an extension of the Euler-Bernoulli beam. 3.2.1 Timoshenko Linear Beam Theory In this theory, the linear beam is taken and the assumption that u 0 is introduced, such that the variations of the energy equations for the cantilever above become: U =EI 0 s=L + ^ sGA w 0 w s=L Z L 0 EI 00 + ^ sGA w 00 0 w ^ sGA w 0 ds (3.16a) Z t 2 t 1 Tdt = Z t 2 t 1 Z L 0 A wwds + Z L 0 I ds dt (3.16b) W nc = Z L 0 (q y w +m z )ds +P y w s=L +M e s=L (3.16c) Using Hamilton's principle and grouping terms, it is seen that the equations of motion for the Timoshenko beam are [54] I EI 00 + ^ sGA w 0 =m z (3.17a) A w ^ sGA w 00 0 =q y (3.17b) 21 with natural boundary conditions EI 0 s=L =M e (3.18a) ^ sGA w 0 s=L =P y (3.18b) 3.2.2 Euler-Bernoulli Beam Theory This theory uses the assumptions of the Timoshenko beam but adds the fact that the beam is slender with =w 0 , such that all shearing terms fall away. With this, the variations of the energy equations for the cantilever above become: U =EIw 00 w 0 s=L Z L 0 EIw 000 w 0 ds (3.19a) Z t 2 t 1 Tdt = Z t 2 t 1 Z L 0 A wwdsdt (3.19b) W nc = Z L 0 q y wds +P y w s=L +M e w 0 s=L (3.19c) The strain energy needs to be integrated by parts again to reduce the order on the variation: U =EIw 00 w 0 s=L EIw 000 w s=L + Z L 0 EIw 0000 wds (3.20) Therefore with Hamilton's principle and grouping terms, it is seen that the equation of motion for the Euler-Bernouli beam is [55] A w +EIw 0000 =q y (3.21) 22 with natural boundary conditions EI 0 s=L =M e (3.22a) EIw 000 s=L =P y (3.22b) 3.2.3 Rayleigh Beam As mentioned earlier, this model is an extension of the Euler-Bernoulli beam, whereby the rotational inertia term is included. To do this, Eq. (3.19b) is replaced by [24]: Z t 2 t 1 Tdt = Z t 2 t 1 Z L 0 A wwds + Z L 0 I w 0 w 0 ds dt = Z t 2 t 1 Z L 0 A wwds Z L 0 I w 00 wds dtI w 0 w s=L (3.23) such that the equation of motion becomes [51] A wI w 00 +EIw 0000 =q y (3.24) with natural boundary conditions EI 0 s=L =M e (3.25a) EIw 000 I w 0 s=L =P y (3.25b) 23 3.3 Nonlinear Beam Models Since the strain energy equation is so complicated, simplications need to be done in order to obtain the governing equations of the beam. If this is not desired then either care needs to be taken to account for all the terms or a nite element code will need to be implemented. This, however, is beyond the scope of this chapter. There are various ways of simplifying this equation, of which two will be highlighted here. The one is to assume small shear and extension by neglecting shearing eects and assuming the bracketed terms after the curvature (third line of Eq. (3.8)) are approximately equal to one. The other is to assume pure bending (the implications of this will be highlighted in the corresponding section). Both methods neglect the ^ y 4 term as the beam is assumed to be thin compared to its length. 3.3.1 Small Shear and Extension With the assumptions given above, namely 1 +u 0 cos +w 0 sin 2 + 1 2 u 0 + 1 2 u 02 +w 02 1 (3.26a) 1 2 Z L 0 ^ sGA w 0 cos 1 +u 0 sin 2 ds 0 (3.26b) The strain energy equation becomes U = 1 2 Z L 0 EA u 0 + 1 2 u 02 +w 02 2 ds + 1 2 Z L 0 EI 02 ds (3.27) 24 Taking the variation of this: U = Z L 0 EA u 0 + 1 2 u 02 +w 02 u 0 +u 0 u 0 +w 0 w 0 ds + Z L 0 EI 0 0 ds (3.28) Now, let N =EA u 0 + 1 2 u 02 +w 02 (3.29) and integrate by parts: U =N 1 +u 0 u +w 0 w L 0 +EI 0 L 0 Z L 0 N 0 1 +u 0 +Nu 00 u + N 0 w 0 +Nw 00 w ds Z L 0 EI 00 ds (3.30) Hamilton's principle can now be used along with Eq's (3.11) and (3.13) to get the equations of motion: I EI 00 =m z (3.31a) A w N 0 w 0 +Nw 00 =q y (3.31b) A u N 0 1 +u 0 +Nu 00 =q x (3.31c) with natural boundary conditions EI 0 s=L =M e (3.32a) Nw 0 s=L =P y (3.32b) N 1 +u 0 s=L =P x (3.32c) 25 These equations are a more general form of the model presented by Washizu [57]. His equations are obtained if the assumptions are made that: the u 02 term in N could be neglected due to very small extension; and the curvature 0 w 00 . Washizu's model is also a static case. Ignoring his assumptions adds an extra equation and 2 boundary conditions, as seen in Eq's (3.31) and (3.32), but the dynamic forms are less complicated. This method also keeps dynamic terms out of the boundary conditions. 3.3.2 Pure Bending as a Function of Vertical De ection (w) It would be convenient to describe the beam with one de ection instead of three, as beams are generally used as a rst order estimation. This can be done by assuming that there is no extension or shearing during deformation. This means that u 0 and become a function ofw such that the vertical de ection can be found directly from the equations [2]. These can be found from Fig. 3.3. ds ds dw dx du θ Figure 3.3: Beam element under pure bending = arcsinw 0 ! 0 = w 00 (1w 02 ) 1=2 (3.33a) u 0 = 1w 02 1=2 1 (3.33b) 26 It should be noted that here u 0 is negative. This is due to the fact that the beam is subjected to pure bending and as such u will increase in a negative direction, as in Fig. 3.3, i.e. du =dsdx. Now, the strain energy simply becomes: U = 1 2 Z L 0 EI w 00 (1w 02 ) 1=2 ! 2 ds (3.34) By writing the bending moment as M =EI w 00 (1w 02 ) 1=2 (3.35) and then taking the variation ofU and integrating by parts twice: U = M (1w 02 ) 1=2 w 0 L 0 M 0 (1w 02 ) 1=2 w L 0 + Z L 0 d 2 M ds 2 + w 0 w 00 1w 02 dM ds 1 (1w 02 ) 1=2 wds (3.36) It can be shown that the variables u and relate to w by = arcsinw 0 ! = w 0 (1w 02 ) 1=2 (3.37a) u = Z s 0 1w 02 1=2 1 d ! u = Z s 0 w 0 w 0 (1w 02 ) 1=2 d (3.37b) where Eq's (3.33a) and (3.33b) have been used. is a dummy variable such that w w (;t). 27 The kinetic energy term then becomes: Z t 2 t 1 Tdt = Z t 2 t 1 " Z L 0 A ww u Z s 0 w 0 w 0 (1w 02 ) 1=2 d ! ds + Z L 0 I w 0 (1w 02 ) 1=2 ds # dt (3.38) Integrating this by parts with respect to s and grouping terms: Z t 2 t 1 Tdt = Z t 2 t 1 Z L 0 A ww U w 0 w 0 (1w 02 ) 1=2 ! dsdt I (1w 02 ) 1=2 w L 0 + Z t 2 t 1 Z L 0 I d ds + w 0 w 00 1w 02 ! 1 (1w 02 ) 1=2 wdsdt (3.39) where = w 0 (1w 02 ) 1=2 + _ w 0 w 0 (1w 02 ) 3=2 (3.40a) U = Z L 0 uds Z s 0 ud (3.40b) u = Z s 0 w 0 w 0 1w 02 + _ w 02 (1w 02 ) 3=2 d (3.40c) Integrating the horizontal acceleration term by parts: Z t 2 t 1 Tdt = Z t 2 t 1 Z L 0 A w + d U ds w 0 + U w 00 1w 02 ! 1 (1w 02 ) 1=2 ! wdsdt + Z t 2 t 1 Z L 0 I d ds + w 0 w 00 1w 02 ! 1 (1w 02 ) 1=2 wdsdt +A Uw 0 (1w 02 ) 1=2 w L 0 I (1w 02 ) 1=2 w L 0 (3.41) 28 Similarly, the variation of the energy due to nonconservative forces becomes W nc = Z L 0 q x Z s 0 w 0 w 0 (1w 02 ) 1=2 d +q y w +m z w 0 (1w 02 ) 1=2 ! ds P x Z L 0 w 0 w 0 (1w 02 ) 1=2 ds +P y w s=L +M e w 0 (1w 02 ) 1=2 s=L (3.42) Integrating by parts and grouping terms: W nc = Z L 0 H w 00 + w 00 w 0 1w 02 q x w 0 +q y 1w 02 1=2 m z w 00 w 0 1w 02 dm z ds 1 (1w 02 ) 1=2 wds P x w 0 (1w 02 ) 1=2 P y ! w s=L +M e w 0 (1w 02 ) 1=2 s=L (3.43) where H is the internal horizontal shear force given by H = Z L 0 q x ds Z s 0 q x d +P x (3.44) Using Hamilton's principle, the equation of motion becomes A w 1w 02 1=2 + d U ds w 0 + U w 00 1w 02 ! I d ds + w 0 w 00 1w 02 ! + d 2 M ds 2 + w 0 w 00 1w 02 dM ds f (w;s) = 0 (3.45) 29 with natural boundary conditions w (0) =w 0 (0) = 0 (3.46a) M s=L =M e (3.46b) dM ds +A Uw 0 I s=L Q e (w (L)) = 0 (3.46c) where the forcing terms are given by f (w;s) =Hw 00 +H w 02 w 00 1w 02 q x w 0 +q y 1w 02 1=2 m z w 0 w 00 1w 02 dm z ds (3.47a) Q e (w (L)) =P x w 0 (L)P y 1w 02 (L) 1=2 (3.47b) 3.3.3 Pure Bending as a Function of Local Rotation () A similar approach can be done by making u 0 and w 0 a function of . This produces much simpler equations. From Fig. 3.3, w 0 = sin (3.48a) u 0 = cos 1 (3.48b) Now, the strain energy simplies to: U = 1 2 Z L 0 EI 02 ds (3.49) 30 Taking the variation and integrating by parts: U =EI 0 L 0 Z L 0 EI 00 ds (3.50) w and u are related to by integrating Eq. (3.48) and then taking the variation: w = Z s 0 sind ! w = Z s 0 cosd (3.51a) u = Z s 0 (cos 1)d ! u = Z s 0 sind (3.51b) The variation of kinetic energy then becomes: Z t 2 t 1 Tdt = Z t 2 t 1 Z L 0 A w Z s 0 cosd u Z s 0 sind ds + Z L 0 I ds dt (3.52) Integrating this by parts with respect to s and grouping terms: Z t 2 t 1 Tdt = Z t 2 t 1 Z L 0 A Z L 0 wds Z s 0 wd cos Z L 0 uds Z s 0 ud sin dsdt Z t 2 t 1 Z L 0 I dsdt (3.53) where w = Z s 0 cos _ 2 sin d (3.54a) u = Z s 0 sin + _ 2 cos d (3.54b) 31 Similarly, the variation of the energy due to nonconservative forces becomes W nc = Z L 0 q x Z s 0 sind +q y Z s 0 cosd +m z ds P x Z s 0 sind s=L +P y Z s 0 cosd s=L +M e s=L (3.55) Integrating by parts and grouping terms: W nc = Z L 0 Z L 0 q x ds Z s 0 q x d +P x sin + Z L 0 q y ds Z s 0 q y d +P y cos +m z ds +M e s=L (3.56) Using Hamilton's principle, the equation of motion becomes A +I EI 00 + Z L 0 q x ds Z s 0 q x d +P x sin Z L 0 q y ds Z s 0 q y d +P y cos =m z (3.57) and natural boundary condition EI 0 s=L =M e (3.58) where, from Eq. (3.53), is the acceleration due to translation and is given by = Z L 0 wds Z s 0 wd cos Z L 0 uds Z s 0 ud sin (3.59) 32 This is a new model which combines the moment equation solution method by Bisshopp and Drucker [10] with the pure bending relating to vertical de ection method by AL- Saddar and AL-Rawi [2] given above. This gives a more robust dynamic model for a pure bending beam. 33 Chapter 4 Shape Memory Polymer Composite Beam Material Model 4.1 1D Stress-Strain Relation For a shape memory polymer strip, its normal stress and normal strain (") are related by [35]: = f E i + 1 f E e 1 "" s Z T T h dT (4.1) where f is the frozen fraction that is only dependent on temperature; " s re ects the strain storage during cooling; E e and E i are the moduli of the entropic deformation and internal energetic deformation, respectively, which are usually found experimentally; is the coecient of thermal expansion; and T h is a reference temperature which is higher than the glass transition temperature T g . By writing 1 E = f E i + 1 f E e (4.2) 34 where E is referred to as the Young's modulus of SMP, Eq. (4.1) is reduced to =E "" s Z T T h dT (4.3) The strain storage can be found from the equation: d" s dT = E E e "" s Z T T h dT d f dT (4.4) The frozen fraction, based on experimental results, can be written as f = 1 1 1 +c f (T h T ) n (4.5) The modulus of entropic deformation is given by E e = 3NkT (4.6) where N is the cross-link density and k is Boltzmann's constant. The modulus of internal energetic deformation E i is constant. Shown in Fig. 4.1 is a typical plot of Young's modulus E versus temperature T. As seen from the plot, E is almost constant when T is below and above the glass transition temperature [35]. Therefore, for the sake of simplicity in analysis, it can be assumed that the SMP takes a valueE h when temperatureT is larger thanT g , and a value andE l when T is smaller than T g , where E l >>E h according to Fig. 4.1. Under this assumption, the 35 stress-strain relation (4.3) of the polymer is piecewise linear with temperature and is of the form =E h " Z T T h dT for TT g (4.7) =E h " Z T T h dT +E l "" g Z T T h dT for TT g (4.8) Of course, in a comprehensive study, Eq. (4.3) can be directly used instead of Eq.'s (4.7) and (4.8). 0.8 0.85 0.9 0.95 1 1.05 10 0 10 1 10 2 10 3 E (MPa) T/Tg Figure 4.1: Young's Modulus as a function of temperature 36 4.2 2D Composite Material The composite material in consideration is made up of an SMP matrix and some unidirectional ber. Without loss of generality, only the unidirectional lamina method will be highlighted herein. This involves a laminate made up of unidirectional laminae, whereby the orientation angle of each lamina can be chosen. The proposed method can easily be extended to a fabric, whereby a laminate is made up of laminae of woven bers set in the SMP matrix. To start, a plate theory is used before it is eventually simplied for the case of a beam. 1 2 σ 1 σ 1 Matrix Fiber (a) 1 2 σ 2 σ 2 W (b) Figure 4.2: (a) Composite material loaded in the 1 principle direction, (b) Composite material loaded in the 2 principle direction 4.2.1 Unidirectional Laminae Figures 4.2(a) and 4.2(b) show a composite material loaded in the 1 and 2 principle directions, respectively. This setting is used to nd the modulus of elasticity in each 37 principle direction. The constitutive equations for the ber and matrix, denoted by subscript f and m respectively are given by f =E f " f (4.9a) m =E m " m " s Z T T h dT (4.9b) where E m and " m are the E and " in Eq. (4.3). Assume that the stored strain is a constant and is not a function of " m , as found from previous knowledge. Because the strain is constant in the 1 principle direction (" f =" m =" 2 ), and because the stress is constant in the 2 direction ( f = m = 2 ), it is found that the combined moduli in each direction are given by [4]: 1 =E 1 " 1 E s 1 " s + Z T T h dT (4.10a) 2 =E 2 " 2 E s 2 " s + Z T T h dT (4.10b) with E 1 =E f V f +E m V m (4.11a) E s 1 =E m V m (4.11b) E 2 = E f E m E f V m +E m V f (4.11c) E s 2 = E f E m V m E f V m +E m V f (4.11d) 38 V f and V m are the ber and matrix volume fractions respectively, with V f +V m = 1. Similarly, it is easy to show that 12 =G 12 12 G s 12 s + Z T T h dT (4.12) where G 12 = G f G m G f V m +G m V f (4.13a) G s 12 = G f G m V m G f V m +G m V f (4.13b) The previous equations are for the cases when laminas are laid at 0 or 90 . The constitutive equations for dierent layup orientations can be derived by coordinate transformation. Consider a lamina (a plate segment) in Fig. 4.3, where 1 and 2 are the principle axes and s andz are the local Cartesian coordinates. The constitutive equations in the coordinates (s and z) are obtained by multiplying the stress and strains in the principle direction by a rotation matrix, thus obtaining [4]: z 1 2 s γ Figure 4.3: Coordinates of a Lamina 39 8 > > > > > > < > > > > > > : s z sz 9 > > > > > > = > > > > > > ; = Q 8 > > > > > > < > > > > > > : " s " z sz 9 > > > > > > = > > > > > > ; Q s 0 B B B B B B @ 8 > > > > > > < > > > > > > : " s s " s z s sz 9 > > > > > > = > > > > > > ; + 8 > > > > > > < > > > > > > : R T T h dT R T T h dT 0 9 > > > > > > = > > > > > > ; 1 C C C C C C A (4.14) where Q = [T ( )] 1 [Q] [R] [T ( )] [R] 1 (4.15a) Q s = [T ( )] 1 [Q s ] [R] [T ( )] [R] 1 (4.15b) with the matrices given by: [Q] = 2 6 6 6 6 6 6 4 E 1 0 0 0 E 2 0 0 0 G 12 3 7 7 7 7 7 7 5 (4.16a) [Q s ] = 2 6 6 6 6 6 6 4 E s 1 0 0 0 E s 2 0 0 0 G s 12 3 7 7 7 7 7 7 5 (4.16b) [T ( )] = 2 6 6 6 6 6 6 4 cos 2 sin 2 2 sin cos sin 2 cos 2 2 sin cos sin cos sin cos cos 2 sin 2 3 7 7 7 7 7 7 5 (4.16c) [R] = 2 6 6 6 6 6 6 4 1 0 0 0 1 0 0 0 2 3 7 7 7 7 7 7 5 (4.16d) 40 and the angle between thes and 1 axes (see Fig. 4.3). The stored strain can be obtained from the matrix form [35] of Eq. (4.4) as follows: df" s g k dT = [E e ] 1 f m g k d f dT = [E e ] 1 Q 0 f"g k Q 0s f" s g k + Z T T h dT k d f dT (4.17) where m is the stress in the matrix, subscript k denotes the k'th layer (as in Fig. 4.4), all vectors are in the s, z and sz direction, as dened following Eq. (4.14), and Q 0 = [T ( )] 1 2 6 6 6 6 6 6 4 E m 0 0 0 E 2 0 0 0 G 12 3 7 7 7 7 7 7 5 [R] [T ( )] [R] 1 (4.18a) Q 0s = [T ( )] 1 2 6 6 6 6 6 6 4 E m 0 0 0 E s 2 0 0 0 G s 12 3 7 7 7 7 7 7 5 [R] [T ( )] [R] 1 (4.18b) 1 2 k-1 N k . . . ŷ N ŷ k ŷ 2 ŷ 1 ŷ s Figure 4.4: Laminate coordinates and numbering 41 4.2.2 Equivalent Modulus of Elasticity for an SMP Composite Laminate Figure 4.4 shows a laminate of N layers with the corresponding local coordinates and laminate numbers. With the assumption that the plate is taken as being linear with no shear in the transverse (^ y) direction, the strain-displacement relations of the laminate are " s =" 0 s ^ y s (4.19a) " z =" 0 z ^ y z (4.19b) sz = 0 sz ^ y sz (4.19c) where " and are the normal and shear strains respectively, is the curvature, ^ y is the transverse distance from the mid-surface (see Fig. 4.4), superscript 0 denotes the measure at the mid-surface, and subscript s, z and sz relate to their respective directions in local coordinates. Substituted Eq's (4.19a) to (4.19c) into Eq. (4.14) yields the stress-strain relation fg = Q " 0 ^ yfg Q s " s0 + ^ yf" s g + " sT + Z T T h dT (4.20) where the stored strain is a sum of a tensile part " s0 , a curvature part " s , and a temperature part " sT in order to nd the stress as a function of strain at the mid-surface 42 and curvature. The vectors are dened in the s,z andsz directions as in Eq.'s (4.14) and (4.19). From the plate theory the force and moment equations for a laminate are given by fNg = N X k=1 Z ^ y k ^ y k1 fgd^ y (4.21a) fMg = N X k=1 Z ^ y k ^ y k1 ^ yfgd^ y (4.21b) Now substite Eq. (4.20) into Eq. (4.21a) and (4.21b) to obtain [4] fNg = [A] " 0 + [B]fg N int (4.22a) fMg = [B] " 0 + [D]fg M int (4.22b) where [A] = N X k=1 Q k Z ^ y k ^ y k1 d^ y = N X k=1 Q k (^ y k ^ y k1 ) (4.23a) [B] = N X k=1 Q k Z ^ y k ^ y k1 ^ yd^ y = 1 2 N X k=1 Q k ^ y 2 k ^ y 2 k1 (4.23b) [D] = N X k=1 Q k Z ^ y k ^ y k1 ^ y 2 d^ y = 1 3 N X k=1 Q k ^ y 3 k ^ y 3 k1 (4.23c) N int = N X k=1 Q s k (^ y k ^ y k1 ) " s0 k + " sT k + Z T T h dT k (4.23d) + 1 2 ^ y 2 k ^ y 2 k1 f" s g k (4.23e) M int = N X k=1 Q s k 1 2 ^ y 2 k ^ y 2 k1 " s0 k + " sT k + Z T T h dT k (4.23f) + 1 3 ^ y 3 k ^ y 3 k1 f" s g k (4.23g) 43 The components of stored strain are found from Eq. (4.17): d " s0 k dT = [E e ] 1 Q 0 " 0 Q 0s " s0 k d f dT (4.24a) df" s g k dT = [E e ] 1 Q 0 fg Q 0s f" s g k d f dT (4.24b) d " sT k dT = [E e ] 1 Q 0s " sT k + Z T T h dT k d f dT (4.24c) Thus, the combination of Eq.'s (4.22a) and (4.22b) yields the following equation 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : N s N z N sz M s M z M sz 9 > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > ; = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 A 11 A 12 A 16 B 11 B 12 B 16 A 22 A 26 B 21 B 22 B 26 A 66 B 61 B 62 B 66 D 11 D 12 D 16 D 22 D 26 D 66 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : " 0 s " 0 z 0 sz s z sz 9 > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > ; 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : N int s N int z N int sz M int s M int z M int sz 9 > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > ; (4.25) 4.2.3 Other Types of Laminae The unidirectional lamina method presented in the previous sections is for detonation of basic concepts of SMPC's. In practice, this method is less used because unidirectional laminae are dicult to manufacture and do not performance as well as some other types of composites. Nonetheless, the unidirectional lamina method provides useful guidance in dealing with fabric or chopped strand composites. 44 A fabric composite is similar to the unidirectional lamina, but since it is made up of a weave, the lamina will be bidirectional. The tows in the weave are generally perpendicular to each other so there will be no rotation around the ^ y axis. There will, however, be rotation matrices associated with the weaving around the s and z axes, as the tows will form sine waves, or similar, in these directions due to the weaving process. This means that now the equivalent modulus of elasticity will be a function of s and z [4,17]. The chopped strand composite method can be devised in a similar way if it is assumed that there will be a certain percentage of randomly distributed bers in each direction [4]. This is a simple way of modeling it but more complex methods do exist. Some popular composite types for use with SMP's are nanocomposites and carbon nanotube composites. These have shown to exhibit the best recovery rates out of all the methods [47]. In order to model these, a new constitutive model is needed but similar beam equations, as seen in section 5.3, can be used. 4.2.4 Constitutive Equation for a Beam 45 The focus of this work is on the governing equations of SMPC beams. For this, simplications will now be made to the plate equations in the previous sections [4]. Firstly, invert Eq. (4.25) to get an equation for strain: 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : " 0 s " 0 z 0 sz s z sz 9 > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > ; = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 11 12 16 11 12 16 22 26 21 22 26 66 61 62 66 11 12 16 22 26 66 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 0 B B B B B B B B B B B B B B B B B B B @ 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : N s N z = 0 N sz M s M z = 0 M sz 9 > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > ; + 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : N int s N int z = 0 N int sz M int s M int z = 0 M int sz 9 > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > ; 1 C C C C C C C C C C C C C C C C C C C A (4.26) For beam theory, allow N z +N int z and M z +M int z to equal zero as the beam's width is small compared to its length. Another simplication is that the composite must have a layup such that 16 16 16 0 so that the normal and shearing eects are uncoupled. It follows that Eq. (4.26) becomes 8 > > > > > > > > > > < > > > > > > > > > > : " 0 s s 0 sz sz 9 > > > > > > > > > > = > > > > > > > > > > ; = 2 6 6 6 6 6 6 6 6 6 6 4 11 11 0 0 11 11 0 0 0 0 66 66 0 0 66 66 3 7 7 7 7 7 7 7 7 7 7 5 0 B B B B B B B B B B @ 8 > > > > > > > > > > < > > > > > > > > > > : N s M s N sz M sz 9 > > > > > > > > > > = > > > > > > > > > > ; + 8 > > > > > > > > > > < > > > > > > > > > > : N int s M int s N int sz M int sz 9 > > > > > > > > > > = > > > > > > > > > > ; 1 C C C C C C C C C C A (4.27) 46 which can be expressed as 8 > > > > > > > > > > < > > > > > > > > > > : N s M s N sz M sz 9 > > > > > > > > > > = > > > > > > > > > > ; = 2 6 6 6 6 6 6 6 6 6 6 4 A B 0 0 B D 0 0 0 0 F C 0 0 C H 3 7 7 7 7 7 7 7 7 7 7 5 8 > > > > > > > > > > < > > > > > > > > > > : " 0 s s 0 sz sz 9 > > > > > > > > > > = > > > > > > > > > > ; 8 > > > > > > > > > > < > > > > > > > > > > : N int s M int s N int sz M int sz 9 > > > > > > > > > > = > > > > > > > > > > ; (4.28) Equation (4.28) shall be used for general beam problems. The extension eects can be further uncoupled from the curvature eects if bending occurs from the neutral axis of bending. This axis is where the curvature is equal to zero, but " 0 s 6= 0. By the rst two lines of Eq. (4.28) N s =A" 0 s N int s M s =B" 0 s M int s 9 > > = > > ; ) M s +M int s = B A N s +N int s =e b N s +N int s (4.29) where e b is the distance from the mid-plane to the neutral axis of bending. For pure bending, N s +N int s = 0, which by the rst line of Eq. (4.28) yields: " 0 s = B A s (4.30) Substituting Eq. (4.30) into the second line of Eq. (4.28) gives M s = D B 2 A s M int s (4.31) 47 Thus, by ignoring shearing eects, the internal forces are related to the strain components by 8 > > < > > : N s M s 9 > > = > > ; = 2 6 6 4 A 0 0 D 3 7 7 5 8 > > < > > : " 0 s s 9 > > = > > ; 8 > > < > > : N int s M int s 9 > > = > > ; (4.32) with D =De 2 b A (4.33) Furthermore, for pure bending in a beam with large de ection, the following substitution can be made in beam theory: EI Db (4.34a) EI" s M int s (4.34b) with b being the width of the beam. 48 Chapter 5 Shape Memory Polymer Composite Cantilever Beam 5.1 Quasi-static Equations A cantilever SMP beam is in large deformation, as shown in Fig. 5.1, where w is the displacement in the y-direction and u is the displacement in the x-direction. Assume that the beam deformation takes place very slowly in a quasi-static process. Without loss of generality, assume that the beam is in pure bending. In sequel, the governing equation quasi-static beam is derived in two cases of local coordinates. s y x w u Figure 5.1: Large deformation of an SMP cantilever beam 49 5.1.1 Quasi-Static Equilibrium Described by the Local Vertical De ection H V M ds q x q y m z M+ ds H+ ds V+ ds dM ds dH ds dV ds (a) s y x M e P x P y q x m z q y (b) Figure 5.2: (a) Free-body diagram for an increment of a cantilever beam, (b) Cantilever Beam with applied body forces and point loads Figure 5.2(a) shows a segment of the beam, where H, V and M are internal forces, and q x , q y and m z are external loads and moment. The quasi-static equilibrium equation with generalized loading conditions (f, M e and Q e ) is given by [2] d 2 M ds 2 + w 0 w 00 1w 02 dM ds =f (w;s) (5.1) and the boundary conditions are w (0) =w 0 (0) = 0 (5.2) 50 M s=L =M e (5.3) dM ds s=L =Q e (w (L)) (5.4) where w is the local vertical de ection of the beam, s is a local coordinate (Fig. 5.1), w 0 =dw=ds, and M =EI w 00 (1w 02 ) 1=2 (5.5) f (w;s) =Hw 00 +H w 02 w 00 1w 02 q x w 0 +q y 1w 02 1=2 m z w 0 w 00 1w 02 dm z ds (5.6) Q e (w (L)) =P x w 0 (L)P y 1w 02 (L) 1=2 (5.7) Here end loads (P x and P y ) or moment (M e ) at the beam tip has been applied, as seen in Fig 5.2(b). H is found to be the integral of q x along the length of the beam. 5.1.2 Quasi-Static Equilibrium Described by the Local Rotation Coordinate Because the beam is in pure bending, its deformation can be conveniently described by the local rotation (s) of the centerline with respect to a horizontal global axis at any point along the beam length. With the local rotation, the resulting governing equations are much simpler than those with the local vertical de ection. To show this, application of the principle of minimum potential energy yields the equilibrium equation. 51 dM ds + Z L 0 q y ds Z s 0 q y ds +P y cos Z L 0 q x ds Z s 0 q x ds +P x sin +m z = 0 (5.8) and the boundary conditions (0) = 0; M e M s=L = 0 (5.9) Note that in the equilibrium equation (5.8) dM ds =EI 00 (5.10) where 00 =d 2 =ds 2 . In general, the rotation is related to the vertical de ection w and the x coordinate by dw ds = sin; dx ds = cos (5.11) It follows that w = Z s 0 sinds (5.12a) x = Z s 0 cosds (5.12b) 52 For a cantilever beam, the vertical de ection and x coordinate of the beam at each node can be approximated by w k = k X i=1 s sin i (5.13a) x k = k X i=1 s cos i (5.13b) The utility of the rotation in modeling and analysis of an SMPC beam has pros and cons, as shown below: Advantages: 1. The equilibrium equation (5.8) are much easier to solve than Eq. (5.1). 2. There are no singularities in numerical simulations. Also, large beam de ections, including that of a beam bent back as shown in Fig. 5.1, can be easily computed and plotted, without having to use a continuation method. 3. The usage of permits various types of external forces, including pointwise moments. Disadvantages: 1. The vertical de ection w cannot be obtained directly from the equilibrium equation. 2. For accurate results, s in the approximated expressions (5.13a) and (5.13b) must be small enough, resulting in more computational eort. 3. Body forces need to be integrated before analysis and computation can be performed. 53 5.2 Temperature Dependence With the material model and governing equations given in the previous sections, the behaviors of the SMPC cantilever beam under dierent temperature and loading conditions can be investigated. To this end, consider three temperature parameters T g , T l and T h , where T g is the glass transition temperature of the SMP material, T l is a typical low temperature (T l <T g ) at which the SMP material is in solid state, and T h is a typical high temperature (T g < T h ) at which the SMP material is in rubbery state. Because only the quasi-static states of the cantilever are considered, the beam behaviors will be investigated at these temperatures. 5.2.1 Piecewise Constant Temperature Dependence for a Pure SMP Beam Assume that the SMP beam is in pure bending such that " =y (5.14) where is the curvature of the beam, and it is given by = d ds = w 00 (1w 02 ) 1=2 (5.15) 54 With the above discussion, the quasi-static equilibrium of the cantilever beam at the high and low temperatures is governed by d 2 M h ds 2 + w 0 h w 00 h 1w 02 h dM h ds =f (w h ;s) for TT g (5.16) d 2 M l ds 2 + w 0 l w 00 l 1w 02 l dM l ds E l E h d 2 M g ds 2 w 0 l w 00 l 1w 02 l E l E h dM g ds =f (w l ;s) for TT g (5.17) with the boundary conditions: Fixed end (s = 0): w h = 0;w 0 h = 0 for TT g (5.18a) w l = 0;w 0 l = 0 for TT g (5.18b) Free end (s =L): M h =M e ; dM h ds =Q e (w h (L)) for TT g (5.19a) M l E l E h M g =M e ; dM l ds E l E h dM g ds =Q e (w l (L)) for TT g (5.19b) where the moment terms at the three temperatures are given by M h =E h I w 00 h 1w 02 h 1=2 (5.20a) M l = E h +E l I w 00 l 1w 02 l 1=2 (5.20b) M g =E h I w 00 g 1w 02 g 1=2 (5.20c) 55 At the glass transition temperature, the beam transverse displacement is described by d 2 M g ds 2 + w 0 g w 00 g 1w 02 g dM g ds =f (w g ;s) (5.21) with boundary conditions: Fixed end (s = 0): w g = 0;w 0 g = 0 (5.22) Free end (s =L): M g =M e ; dM g ds =Q e (w g (L)) (5.23) It is easy to shown from Eq's. (5.16), (5.17) and (5.21) that at T =T g w h =w l =w g (5.24) Equation (5.24) assures displacement continuity through the entire temperature range. 5.2.2 Full Temperature Dependant Model Assume that the SMPC beam is in pure bending such that " =y (5.25) where the beam curvature is given by = 0 = w 00 (1w 02 ) 1=2 (5.26) 56 The quasi-static equilibrium of the cantilever beam at the high and low temperatures is described by Eq's. (5.1) to (5.4) or (5.8) and (5.9), with M =EI (" s ); dM ds =EI 0 d" s ds (5.27) where the strain storage " s is found by solving Eq. (4.24b), with Q 0 = Q 0 11 and Q 0s = Q 0s 11 . Note that EI = Db and EI" s =M int s for the composite beam. If Young's modulus is assumed piecwise constant at T h and T l , as in Eq's. (4.7) and (4.8), it is necessary to check for continuity of the displacement at T g (w h = w g = w l ), see Eq. (5.24). On the other hand, the nonlinear temperature dependency of Young's modulus (as shown in Fig. 4.1) can be simply used in the constitutive equations given in Section 4.1, with guaranteed displacement continuity. 5.3 Shape Fixation and Shape Recovery In this section, the physical behaviors of the SMPC beam in shape xation and shape recovery are examined. 5.3.1 Cooling for Shape Fixation In shape xation, the beam is rst loaded incrementally at a high temperature T h ; the beam then remains under xed loading while the temperature is decreased to a low temperature T l , where T h >T g >T l . After the temperature reaches T l , all loads will be removed. The cooling of the SMPC beam (from T l to T h ) is assumed to be a slow process 57 such that the quasi-static model given in Section 5.1 is sucient to describe the beam de ection. The governing equations for this process are given as follows. (i) At T =T h , external loads are applied. The beam de ection is governed by D h b 00 + Z L 0 q y ds Z s 0 q y ds +P y cos Z L 0 q x ds Z s 0 q x ds +P x sin +m z = 0 (5.28) subject to boundary conditions (0) = 0 (5.29a) D h b 0 s=L =M e (5.29b) (ii) From T h to T l , the external loads are held constantly as temperature decreases. The beam de ection is governed by D (T )b 00 dM int s ds + Z L 0 q y ds Z s 0 q y ds +P y cos Z L 0 q x ds Z s 0 q x ds +P x sin +m z = 0 (5.30) subject to boundary conditions (0) = 0 (5.31a) D (T )b 0 M int s s=L =M e (5.31b) 58 Here, D(T ) indicates the nonlinear temperature dependency of the equivalent Young's modulus. (iii) At T =T l , all forces are removed and the beam de ection is governed by D l b ^ 00 dM int s ds = 0 (5.32) subject to boundary conditions ^ (0) = 0 (5.33a) D l b ^ 0 M int s s=L = 0 (5.33b) As a special case, consider a pure SMP beam. Through the use of " E =E h =E l , and with the assumption that all the strain from T h has been stored as " s , Eq.'s (5.32) to (5.33b) are reduced to the equilibrium equation (1 +" E )E h I ^ 00 =E h I 00 = Z L 0 q y ds Z s 0 q y ds +P y cos + Z L 0 q x ds Z s 0 q x ds +P x sinm z (5.34) 59 subject to boundary conditions ^ (0) = 0 (5.35a) (1 +" E )E h I ^ 0 s=L =E h I 0 s=L =M e (5.35b) Because E h =E l << 1, Eq.'s (5.34) to (5.35b) are almost the same as Eq.'s (5.28) to (5.29b). Thus, it can be concluded that the de ection of the pure SMP beam at T l (with loads removed) is only slightly dierent from that of the beam at T h with loads applied. 5.3.2 Heating for Shape Recovery In this process, without loss of generality, assume that the beam is heated from T l to T h in a linear fashion T =T l +t; t 0 (5.36) The heating process is intrinsically a dynamic one. This means that the quasi-static beam model needs to be modied with an added inertia term. Thus, at the beginning of shape recovery, the equation of motion of the beam is given by [9]: A +I D (T )b 00 + dM int s ds = 0 (5.37) 60 for 0 t t 1 , with t 1 = (T h T l )=, where = (s;t), with t being the temporal parameter. The boundary conditions of the beam are (0;t) = 0 (5.38a) D (T )b 0 M int s s=L = 0 (5.38b) and the initial conditions are (s; 0) = ^ l0 (s); _ (s; 0) = 0 (5.39) 5.3.3 Vibration after Shape Recovery Once the shape recovery is completed, there is a nal process in which the beam is cooled down fromT h toT l and no loads are applied during this time period. The equation of motion of the beam after the cooling process is completed can be written as A +I D l b 00 + dM int s ds Z L 0 q y ds Z s 0 q y ds +P y cos + Z L 0 q x ds Z s 0 q x ds +P x sinm z = 0 (5.40) which is subject to boundary conditions (0;t) = 0 (5.41a) D l b 0 M int s s=L =M e (5.41b) 61 and the initial conditions (s; 0) = 0; _ (s; 0) = 0 (5.42) 62 Chapter 6 Discretization and Numerical Methods According to Bisshopp and Drucker [10], it is possible to obtain analytical expressions for the de ection of the cantilever beam at its tip by elliptic integrals. These analytical results, however, are only valid for a few simple cases, and are not applicable to the SMPC beam in the current study. For general loadings, as shown in Fig. 6.1, discretization of the governing equations in Section 5.1 must be utilized. In this study, a nite dierence scheme is developed, which eventually results in the following nonlinear matrix algebraic equation A (u;q;T ) = 0 (6.1) where vector A is a nonlinear function of displacement vector u, load vector q, and temperature T . Solution of Eq. (6.1) gives the de ection and internal stress of the SMPC beam at dierent temperatures. This scheme is described as follows. 63 s y x M e P x P y q x m z q y Figure 6.1: A cantilever SMPC beam with general loading conditions 6.1 Finite Dierence Discretization in Terms of Vertical De ection In Fig. 6.2, the beam is divided into segments by m nodes in local coordinates, where each node undergoes transverse de ection w. The global coordinate system is described in x and y, while the local coordinate system is described by the undeformed arc length coordinate s only. By central nite dierence spatial derivatives of w with respect to s are approximated by s y w i-1 w i h i m L θ i Δ x Figure 6.2: Cantilever beam discretized by a nite dierence approximation in local coordinates 64 w 0 i = w i+1 w i1 2h (6.2) w 00 i = w i+1 2w i +w i1 h 2 (6.3) where h = s. Similarly, the derivatives of the moment of the beam are expressed as dM ds = M i+1 M i1 2h (6.4) d 2 M ds 2 = M i+1 2M i +M i1 h 2 (6.5) where by Eq. 5.5 M i =EI w 00 i 1w 02 i 1=2 (6.6) Without loss of generality, assume that EI is constant through the length of the beam. The proposed discretization, however, is applicable to non-uniform beams, such as a fabric weave composite beam. Substituting Eq.'s (6.2) to (6.6) into Eq. (5.1) yields a set of nite dierence equations of the form of Eq. (6.1). These equations contain \virtual" displacements w 1 , w 0 , w m+1 and w m+2 . These quantities can be replaced by \true" displacements (w 1 , w 2 , ..., w m ) though examination of the cantilever boundary conditions of the beam as follows. By Eq. (5.2), it is found that w 0 = 0 (6.7) and w 1 =w 1 (6.8) 65 Also, discretization of Eq. (5.3) yields a nonlinear equation for w m+1 EI 2 (w m+1 2w m +w m1 ) h 4h 2 (w m+1 w m1 ) 2 1=2 M e = 0 (6.9) which can be solved by applying Newton's method. Instead of discretizingdM=ds in terms of w in Eq. (5.4) and nding w m+2 , simply discretize it in terms of Eq. (6.4) and solve for M m+1 . Thus, M m+1 =M m1 + 2h P x w 0 m P y 1w 02 m 1=2 (6.10) The equilibrium equation is solved iteratively from i = 1 to m. The above discretization is with respect to the local coordinates. The relation between s and the global coordinate x can be obtained according to Fig. 6.2. This can be done as follows x i =x i1 + q h 2 (w i w i1 ) 2 (6.11) Equation (6.11) assumes thatx i is greater thanx i1 . When the beam has large de ection or is signicantly bent as in Fig. 5.1, x i may be smaller than x i1 . In that case, discretization of dx=ds gives x i =x i2 + q 4h 2 (w i w i1 ) 2 (6.12) 66 6.2 Finite Dierence Discretization in Terms of Rotation Consider the m nodes in Fig. 6.2, where each node has a local rotation angle . By central dierence, spatial derivatives of the rotation coordinate are approximated by 0 i = i+1 i1 2h (6.13) 00 i = i+1 2 i + i1 h 2 (6.14) with h = s. It follows that Eq. (5.10) is used to EI 00 i + Z L 0 q y ds Z s 0 q y ds +P y cos i Z L 0 q x ds Z s 0 q x ds +P x sin i +m z = 0 (6.15) which represents a set of nite dierence equations, as Eq. (6.1). Again, these equations contain \virtual" rotation angles 0 and m+1 , which can be replaced by the \true" rotation angles as follows. By Eq. (5.9), it is found that 0 = 0 (6.16) Also, discretization of Eq. (5.9) yields an equation for m+1 m+1 = 2hM e EI + m1 (6.17) The nite dierence equations are solved iteratively from i = 1 to m. 67 The above discretization is with respect to the local coordinates. The relation between s and the global coordinate x and de ection w can be obtained according to Fig. 6.2. This can be done as follows x i =x i1 +h cos i1 (6.18) w i =w i1 +h sin i1 (6.19) 6.3 Comparison of Methods Used to Model the Quasi- Static Nonlinear Cantilever Beam Before numerical results are presented, it is important to compare the dierent methods that could be used to describe, model and solve the nonlinear beam. Various methods were analysed and tested. A description of each is given below. 6.3.1 Coordinate Systems The coordinate system relates to how the beam is described and how the equations are derived. (i) 2nd Order Global Coordinates (Moment Equation) [46] The eective moment applied to the beam can be described using a global coordinate system 68 Table 6.1: Advantages and disadvantages of 2nd order global coordinates Advantages Disadvantages Equations are easy to derive Easy to describe forces in x and y di- rections at any point along the beam Can be treated as an initial value problem rather than a boundary value problem Distributed loads and loads applied along the length of the beam need to be recalculated during simulation Need to account for end-shortening of the beam Applied moment along the length of the beam yields unstable results Cannot convert into dynamic simula- tion Cannot have two de ection values for one x-value (ii) 4th Order Global Coordinates [61] The equilibrium beam equation can be derived in terms of vertical de ection in a global coordinate system Table 6.2: Advantages and disadvantages of 4th order global coordinates Advantages Disadvantages Equations are easy to derive Easy to describe forces in x and y directions Can easily convert to a dynamic model Distributed loads and loads applied along the length of the beam need to be recalculated during simulation Need to account for end-shortening of the beam Applied moment to the end of the beam yields incorrect results Cannot have two de ection values for one x-value 69 (iii) Local Coordinates Described by w [2] The equilibrium beam equation can be derived in terms of vertical de ection in a local coordinate system Table 6.3: Advantages and disadvantages of 4th order local coordinates Advantages Disadvantages Can easily convert to a dynamic model Works for all loading conditions Will allow the beam to bend back on itself Equations are dicult to derive and implement There is a singularity at the turning point of the beam (iv) Local Coordinates Described by The equilibrium beam equation can be derived in terms of the rotation at each point along the beam Table 6.4: Advantages and disadvantages of curvilinear coordinates Advantages Disadvantages The equilibrium equations are much easier to solve There are no singularities The beam and force vs w plot can bend back on itself without using any continuation methods Any force can be applied to the beam, including a moment applied at any point along the beam w cannot be solved directly from the equilibrium equations Equations (5.12a) and (5.12b) are not very accurate so a smaller s is needed to improve the accuracy of the solution Body forces need to be integrated be- fore any analysis can be done 70 6.3.2 Solution Methods The solution method dictates how the beam equations will be solved, be it analytically or numerically. (i) Elliptic Integrals [10] Bisshopp and Drucker proposed a method for calculating the de ection of the tip of a cantilever beam with large de ection. Through the derivation, a closed form solution is obtained using elliptic integrals Table 6.5: Advantages and disadvantages of using elliptic integrals Advantages Disadvantages Closed form solution Elliptic integral are dicult to solve. Need to solve them numerically Can only solve for the tip de ection of the beam Dicult or impossible to solve for combined loading (ii) State-space and Runge-Kutte The equilibrium equations can be transformed into state space and solved using a Runge-Kutte scheme 71 Table 6.6: Advantages and disadvantages of using state-space and Runge-Kutte Advantages Disadvantages Easy to implement basic structure Cannot apply loads at a specic nodes Very dicult and inaccurate to apply end-shortening (iii) Finite Dierence The beam can be discretized into a number of nodes in order to calculate the displacement at each point. The derivative terms in the equilibrium equations can then be discretized with a central dierence formulation Table 6.7: Advantages and disadvantages of using nite dierence method Advantages Disadvantages Easy to implement Describes all points along the beam Introduces numerical error 6.3.3 Calculation Methods The calculation methods are the numerical methods that will be used to solve the nonlinear matrix described by Eq. (6.1). (i) MATLAB'S FSOLVE and FZERO functions MATLAB has built in solvers to nd unknown variable or variables for a function equal to zero 72 Table 6.8: Advantages and disadvantages of using MATLAB functions Advantages Disadvantages Very easy to implement Can be slow (ii) Single Variable Netwon's Method The 2nd order system can be solved by starting at the rst node and calculating each node thereafter based on the previous nodes. This can easily be done with Newton's Method Table 6.9: Advantages and disadvantages of using single variable Newton's method Advantages Disadvantages Easy to implement for a 2nd order system Slow when de ections get large Can only be used for a second order system (iii) Multi-variable Netwon's Method In order to combat the slowness of the single variable Newton's method and to be able to solve the 4th order system, the multi-variable Newton's method can be used 73 Table 6.10: Advantages and disadvantages of using multi-variable Newton's method Advantages Disadvantages Fast and ecient Dicult to calculate Jacobian matrix (iv) Continuation Methods [45] These methods are very similar to Newton's Method but they allow for the beam to bend back on itself Table 6.11: Advantages and disadvantages of using continuation methods Advantages Disadvantages Allows beam to bend back on itself Dicult to calculate Jacobian matrix May need to invert non-symmetric Jacobian matrix Dicult to implement 6.3.4 Selection A model for each of the described methods was made and tested in order to properly decide which methods would be best. It was decided that the local coordinates described by beam would be used to describe the beam due to the advantage of having the beam able to be formed to any shape. A nite dierence method was chosen as it could easily be implemented once the required beam equations were derived. As for the solver, a multi-variable Newton's method was used due to the ease of implementation. Also, due 74 to the fact that curvilinear coordinates were used, there was no need to implement a continuation method. The fourth-order Runge-Kutta method was used to solve the stored strain equation. 75 Chapter 7 Finite Element Derivations The analytical beam model is useful as a rst order approximation, however a more general robust model is needed to describe the SMPC under any condition. Hence, a nite element (FE) model is introduced. The following chapters will develope FE models for the SMPC structure. 7.1 Continuum Finite Element Equations 7.1.1 Overview The principle of virtual work in a Lagrangian system states that [8] Z ji u i;j d Z b i u i d n SD X i=1 Z t i t i u i d = 0 (7.1) where ij are the Cauchy stress components, u i are the displacement components, b i are the body forces, t i are the surface tractions, is the domain with boundary , is the density, is the variational operator, and n SD is the number of space dimensions. Here 76 the notation ( ) ;j is for a spatial derivative in the deformed coordinate system. This will be highlighted more later. Now, in order to implement a nite element formulation, approximate the undeformed coordinates X i to X i =X iI N I (7.2) with N I being the shape functions at node I, and X iI are the undeformed coordinates at node I. Also approximate the deformed coordinates to x i =x iI N I (7.3) wherex iI are the deformed coordinates at node I. Thus the displacement eld is given by u i =x i X i =u iI N I (7.4) with the nodal displacements u iI given by u iI =x iI X iI (7.5) Taking the variation of the displacement eld u i =u iI N I (7.6) 77 where u iI are the virtual nodal displacements. Substituting this into the virtual work equation (Eq. (7.1)) u Ii Z N I;j ji d u Ii Z N I b i d n SD X i=1 u Ii Z t i N I t i d = 0 (7.7) Considering this is valid for any arbitrary virtual work, this can be simplied to give the weak form Z N I;j ji d Z N I b i d n SD X j=1 Z t j N I t i d = 0 8 (I;i) = 2 u i (7.8) or f int iI f ext iI = 0 (7.9) where f int iI and f ext iI are the internal and external nodal forces, given by f int iI = Z N I;j ji d (7.10a) f ext iI = Z N I b i d + n SD X j=1 Z t j N I t i d (7.10b) From here it is important to decide if the analysis will be done with respect to the deformed coordinate system (Updated Lagrangian method) or with respect to the undeformed coordinate system (Total Lagrangian method). 78 7.1.2 Updated Lagrangian Method The Updated Lagrangian method uses the deformed coordinate system to perform the nite element analysis. To do this the Caucy stress () and the Almansi strain (") will be used. The Almansi strain is given by " ij = 1 2 @u i @x j + @u j @x i @u k @x i @u k @x j (7.11) The equilibrium equation is given as above in Eq. (7.9) with the forces as f int iI = Z dN I dx j ji d (7.12a) f ext iI = Z N I b i d + n SD X j=1 Z t j N I t i d (7.12b) In nite elements, it is easier to analyse the structure if the shape functions are described by a parent element coordinate system. This means that a mapping must be done to convert the cartesian coordinates to parent coordinates. The current (deformed) element domain will be and parent element domain will be . As such the deformed coordinate will be given by x i =x iI N I () (7.13) where are the parent element coordinates. Similarly the undeformed coordinate will be X i =X iI N I () (7.14) 79 and consequently the displacement and derivative of displacement is u i =u iI N I () (7.15a) @u i @x j =u iI dN I dx j () (7.15b) From Eq.'s (7.11) and (7.12a) it is clear that the derivatives on N I are needed. This is done by the chain rule: dN I dx j = dN I d k d k dx j = dN I d k @x j @ k 1 (7.16) Here, by Eq. (7.13) dx j d k =x jI dN I d k (7.17) For the mapping, the element Jacobian determinant (J ) is needed. This is found from J = det dx d (7.18) Thus mapping the internal force (Eq. (7.12a)) f int iI = Z dN I dx j ji d = Z dN I dx j ji J d (7.19) 80 7.1.3 Total Lagrangian Method The Total Lagrangian method uses the undeformed coordinate system to perform the nite element analysis. To do this the 2nd Piola-Kirchho stress (S) and the Green strain (E) will be used. The Green strain is given by E ij = 1 2 @u i @X j + @u j @X i + @u k @X i @u k @X j (7.20) The equilibrium equation is given as above in Eq. (7.9) with the forces converted to the undeformed domain ( 0 ) f int iI = Z 0 dN I dX k @X k @x j ji Jd 0 (7.21a) f ext iI = Z 0 N I b i d 0 + n SD X j=1 Z 0 t j N I t 0 i d 0 (7.21b) where J is the determinant of the Jacobian given by J = det @x @X (7.22) but Cauchy stress can be dened as PK2 stress by ji =J 1 @x j @X k S km @x i @X m ) dX k dx j ji J =S km @x m @X i =S km im + @u i @X m (7.23) therefore the internal nodal force is given by f int iI = Z 0 dN I dX k S km im + @u i @X m d 0 (7.24) 81 Similarly to the Updated Lagrangian formulation, the cartesian coordinates must be mapped to parent coordinates. The initial (undeformed) element domain will be 0 and parent element domain will be . The coordinates are given by Eq.'s (7.13) and (7.14) but the displacement and its derivative will be u i =u iI N I () (7.25a) @u i @X j =u iI dN I dX j () (7.25b) From Eq.'s (7.20) and (7.21a) it is clear that the derivatives on N I are needed. This is done by the chain rule: dN I dX j = dN I d k d k dX j = dN I d k @X j @ k 1 (7.26) Here, by Eq. (7.14) dX j d k =X jI dN I d k (7.27) For the mapping, the undeformed element Jacobian determinant (J 0 ) is needed. This is found from J 0 = det dX d (7.28) Thus mapping the internal force (Eq. (7.24)) f int iI = Z dN I dX k S km im + @u i @X m J 0 d (7.29) 82 7.1.4 Numerical Quadrature In order to perform the integration over the parent element it is useful to use Gauss quadrature. Since the parent element is structured such that it is over the interval [-1,1] in all dimensions, the quadrature states (in 3 dimensions for example) Z f ()d = Z 1 1 Z 1 1 Z 1 1 f ()ddd = n Q 1 X Q 1 =1 n Q 2 X Q 2 =1 n Q 3 X Q 3 =1 f ( Q 1 ; Q 2 ; Q 3 )w Q 1 w Q 2 w Q 3 (7.30) where, and are the parent element coordinates,Q i are the quadrature points,w Q i are the quadrature weights, and n Q i are the number of quadrature points in each dimension. 7.2 Elements 7.2.1 Continuum-Based (CB) Beam When modeling beams using nite element analysis, there are a number of dierent models that can be used to describe the beam. These fall into two categories: structural beam elements and continuum-based beam elements. The structural elements use the energy forms of the classical beam theories (such as Euler-Bernoulli, Timoshenko or Rayleigh beam models) to formulate the beam equations. A more general structural beam element could be derived from the beam equations previously presented. Continuum-Based (CB) beam elements use the equations derived above but make some assumptions in order to t with beam theory. The advantage that the CB beam has over the structural elements is that no assumptions need to be made based on the amount of extension or shearing in the beam. These elements also allow for the use of isoparametric elements, similar to 83 those used in the continuum modeling above. The structural elements are easily derived from the energy formulation previously presented using the same method as above to form a nite element equilibrium equation. The presentation following will present the CB beam element based o the work by Belytschko. [8] 7.2.1.1 Element Derivation In order to model the beam, the following assumptions need to be made: 1. The cross-section remains straight 2. The transverse normal stress is negligible (i.e. plane stress in y) 3. The beam is thin in the width direction (i.e. plane stress in z) 4. The cross-section does not deform (This does not have to be the case [50] but it makes the formulation easier) Beam theory requires that the beam be described by the neutral axis. To do this, the concept of master and slave nodes will be introduced. Any number of nodes can be used in a single beam element but for simplicity sake a 2 node beam will be used (as seen in Fig. 7.1) 4 * ,(1 + ) 3 * ,(2 + ) 1 ^ ,(1 - ) 2 ^ ,(2 - ) 4 * ,(1 + ) 3 * ,(2 + ) 1 ^ ,(1 - ) 2 ^ ,(2 - ) 1 2 ξ η θ 1 θ 2 Parent Element Figure 7.1: A 2 node CB Beam element 84 The master nodes are shown by the black circles in Fig. 7.1. The position of each of these nodes is described by three components: the horizontal de ection u; the vertical de ectionw; and the angle of the cross-section with respect to the horizontal . Using the assumptions that the cross-section remains straight and is non deformable, these master nodes can be extrapolated to slave nodes. Each master node will form 2 slave nodes, one on the top and one on the bottom face of the element. The numbering convention can be either: the number of the master node with a superscript \+" for the top node and \" for the bottom node; or the general counterclockwise numbering convention with a superscript \" for the top node and \ ^ " for the bottom node. The isoparametric map for the slave nodes can therefore be written as: x = n N X I + =1 x I +N I + (;) + n N X I =1 x I N I (;) = 2n N X I =1 x I N I (;) (7.31) where N I (;) are the shape functions for an isoparametric element for continua, and n N is the number of master nodes. The slave nodes can be found from the master nodes with use of the assumption that the cross-section is inextensible x I + =x I + 1 2 h 0 I p I x I =x I 1 2 h 0 I p I (7.32) with h 0 I being the initial beam thickness at node I, and p I the director at node I. The director at each node is given by p I = 1 h 0 I (x I +x I ) =e x cos I +e y sin I (7.33) 85 Here I is the rotation of the cross-section with respect to the horizontal at each master node. The forces at the master nodes f mast I are found using a transformation from the slave nodes to the master nodes, given by: f mast I = 8 > > > > > > < > > > > > > : f xI f yI m I 9 > > > > > > = > > > > > > ; mast =T T I 8 > > < > > : f I f I + 9 > > = > > ; slave = 2 6 6 6 6 6 6 4 1 0 1 0 0 1 0 1 y I y I x I x I y I y I + x I +x I 3 7 7 7 7 7 7 5 8 > > > > > > > > > > < > > > > > > > > > > : f xI f yI f xI + f yI + 9 > > > > > > > > > > = > > > > > > > > > > ; (7.34) where f xI and f yI are the forces in the horizontal and vertical directions on master node I, and m I is the moment on master node I. In order to use the assumptions of plane stress in the transverse (y) and width (z) directions, it is important that the solution method follow deformation of the beam. This can be done in two ways: either the Updated Lagrangian method be used with a corotational formulation implemented [8]; or the Total Lagrangian method be used. 7.2.1.2 Updated Lagrangian Beam Formulation To use the Updated Lagrangian method, it is necessary to use the laminar components of stress and strain. To do this a corotational formulation must be implemented. Thus, corotational unit vectors ^ e x and ^ e y are introduced, as seen in Fig. 7.2. 86 ê y ê x p(-1) p(1) Figure 7.2: A beam element with corotational unit vectors ^ e x and ^ e y , and director vectors p () where ^ e x is tangent to the arc, and ^ e y is normal to the arc. These are given by: ^ e x = x ; kx ; k = x ; e x +y ; e y x 2 ; +y 2 ; 1=2 (7.35a) ^ e y = y ; e x +x ; e y x 2 ; +y 2 ; 1=2 (7.35b) where x ; = X I x I N I ; (;) (7.36a) y ; = X I y I N I ; (;) (7.36b) The constitutive model in corotational system will then be of the form ^ =F (^ ") (7.37) 87 where the corotational strain (^ ") is given by ^ " =R T lam "R lam where R lam = 2 6 6 4 e x ^ e x e x ^ e y e y ^ e x e y ^ e y 3 7 7 5 (7.38) However, based on the continuum equations formulated above, f xI f yI int = n Q 1 X Q 1 =1 n Q 2 X Q 2 =1 N I ;x N I ;y 2 6 6 4 xx xy xy yy 3 7 7 5 bJ w Q 1 w Q 2 (7.39) whereb is the width dimension. This needs to be converted to the corotational coordinates. Doing this, the equation for internal forces is obtained: f xI f yI int = n Q 1 X Q 1 =1 n Q 2 X Q 2 =1 N I ;^ x N I ;^ y 2 6 6 4 ^ xx ^ xy ^ xy ^ yy 3 7 7 5 2 6 6 4 R x^ x R y^ x R x^ y R y^ y 3 7 7 5 bJ w Q 1 w Q 2 (7.40) This equation, however, will cause shear locking. This cannot be overcome using the general selective-reduced integration, instead a beam specic reduced integration must be used. This is done by only evaluating each element with a single-stack of quadrature points at = 0 (for the 2-node CB beam element). Plane stress must also be implemented by setting ^ yy = 0. Together this gives f xI f yI int = h h 0 n Q X Q=1 0 B B @ N I ;^ x N I ;^ y 2 6 6 4 ^ xx ^ xy ^ xy 0 3 7 7 5 2 6 6 4 R x^ x R y^ x R x^ y R y^ y 3 7 7 5 bJ w Q 1 C C A (0; Q) (7.41) 88 The factor h=h 0 accounts for the change in thickness of the beam at the single-stack of quadrature points. The weight factor w Q is a weight for the single stack (i.e. accounting for single point in ). This reduced integration is one of the most important features of the CB beam. Some commercial packages (such as ABAQUS) do not incorporate such reduced integrations and as such cannot model beams using 2D or 3D continuum elements. These packages typically use structural elements to model beams. 7.2.1.3 Total Lagrangian Beam Formulation The fact that the Total Lagrangian method uses a xed reference frame means that a corotational formulation is unnecessary. As such, the formulation is much simpler. The constitutive equation will be of the form S =G (E) (7.42) Based o the previous Total Lagrangian formulation, the internal forces are given by: f xI f yI int = n Q 1 X Q 1 =1 n Q 2 X Q 2 =1 N I ;X N I ;Y 2 6 6 4 S xx S xy S xy S yy 3 7 7 5 0 B B B @ 2 6 6 4 1 0 0 1 3 7 7 5 2 6 6 4 u x;X u y;X u x;Y u y;Y 3 7 7 5 T 1 C C C A bJ 0 w Q 1 w Q 2 (7.43) 89 But, based on the need for reduced integration and to enforce plane stress, as described above, the formulation yields: f xI f yI int = n Q X Q=1 0 B B @ N I ;X N I ;Y 2 6 6 4 S xx S xy S xy 0 3 7 7 5 0 B B B @ 2 6 6 4 1 0 0 1 3 7 7 5 2 6 6 4 u x;X u y;X u x;Y u y;Y 3 7 7 5 T 1 C C C A bJ 0 w Q 1 C C C A (0; Q) (7.44) 7.2.1.4 3D CB Beam Equations Both methods are easily extended into 3D by: 1. mapping the slave nodes out to form 8-node isoparametric quad elements 2. the reduced integration will be done in the and direction at = 0 3. enforce plane stress by letting ^ yy = 0 and ^ zz = 0 (orS yy = 0 andS zz = 0 in Total Lagrangian) and not in the constitutive model. 7.2.2 Shell Elements For various applications, the beam element may not be good enough. This could occur when 2D shearing eects are needed (especially with composites) or if a plate like structure is needed. To accomplish this, a plate element or a shell element can be used. A plate element can be derived from plate theory. This however causes it to be limited due to the assumptions made. A more general element is the shell element. This is analogous to the continuum-based beam element for plates. 90 7.2.2.1 Element Derivation Similarly to the CB beam element, various assumptions need to be made, namely: 1. The cross-section remains straight 2. The transverse normal stress is negligible (i.e. plane stress in y) Shell theory requires that the structure be described by the neutral plane. To do this, the concept of master and slave nodes will be introduced. Any number of nodes can be used in a single shell element but for simplicity sake an 8 node shell will be used (as seen in Fig. 7.3) 2* 6* 1^ 5^ 2 * 6 * 1 ^ 5 ^ 7 8 ξ η 3* 7* 4^ 8^ Parent Element ζ 9* 10^ 14* 16^ 13^ 12^ 15* 11* 3 * 7 * 4 ^ 8 ^ 4 9* 10^ 6 1 3 11* 12^ 2 5 14* 15* 16^ 13^ Figure 7.3: An 8 node shell element The master nodes are shown by the black circles in Fig. 7.3. The position of each of these nodes is described by three components: the horizontal de ection u; the vertical de ection w; and 2 rotations of the cross-section vector. Using the assumptions that the cross-section remains straight, these master nodes can be extrapolated to slave nodes. Each master node will form 2 slave nodes, one on the top and one on the bottom face of the element. The numbering convention can be either: the number of the master node with a superscript \+" for the top node and \" for the bottom node; or the general 91 counterclockwise numbering convention with a superscript \" for the top node and \ ^ " for the bottom node. This is equivalent to the beam element and as such, the derivation is very similar. The major change comes in Eq. (7.32), as now this must be done in 3D. There are various ways to do this. Two common ways of doing this are: to assign a vector connecting the top and bottom nodes [1], as in Fig. 7.4(a); or to describe the line connecting the two nodes by 2 rotations [44], as in Fig. 7.4(b). n i n i + n i - x V 1 V 2 V 3 α β (a) x n i - y z n i + n i θ 1 θ 2 (b) Figure 7.4: (a) Rotation described by a vector, (b) Rotation described by angles For the rst method (Fig. 7.4(a)), the vector connecting the top and bottom node V 2 can be constructed from rotations and about 2 orthogonal base vectors V 1 and V 3 respectively. It is obvious that there are an innite number of orthogonal vectors to V 2 and as such a specic scheme is need. The scheme generally used is as follows: Choose V 3 such that it is orthogonal to the x axis and V 2 V 3 =iV 2 (7.45) 92 V 1 can then be chosen to be orthogonal to V 2 and V 3 V 1 =V 2 V 3 (7.46) From here, it can now be seen that for each node V 2I =h 0 I v 1I ; v 3I 8 > > < > > : 9 > > = > > ; (7.47) whereh 0 I is the thickness at nodeI, V 2I is the change inV 2I at nodeI, andv 1I andv 3I are the normalised vectors of V 1I and V 3I at node I. This then yields from Eq. (7.32): p I = 1 h 0 I V 2I = 1 h 0 I V 0 2I + V 2I (7.48) where V 0 2I is the initial vector at node I. The problem with this method is that it is only valid for small rotations. This would make it useful for Updated Lagrangian method but not for Total Lagrangian. There is also a problem when the vector lies along the x-axis. In this case, one would use the same method except with the y-axis. For these reason it would be more useful to use a method described by angles. For the second method (Fig. 7.4(b)), 2 angles need to be chosen in order to describe the rotation of the line connecting the top and bottom nodes. The angles, as per the 93 diagram, are chosen in such a way that there is not a singularity when the line is in the xz-plane. Thus, Eq. (7.33) becomes p I = 1 h 0 I (x I +x I ) =e x cos 1I sin 2I +e y sin 1I sin 2I +e z cos 2I (7.49) This method can be thought of as being equivalent to the method used for the beam elements. The forces at the master nodes f mast I are found using a transformation from the slave nodes to the master nodes, given by: f mast I = 8 > > > > > > > > > > > > > > < > > > > > > > > > > > > > > : f xI f yI f zI m zI m xI 9 > > > > > > > > > > > > > > = > > > > > > > > > > > > > > ; mast =T T I 8 > > < > > : f I f I + 9 > > = > > ; slave = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 h 0 I 2 s 1 s 2 h 0 I 2 c 1 s 2 0 h 0 I 2 s 1 s 2 h 0 I 2 c 1 s 2 0 h 0 I 2 c 1 c 2 h 0 I 2 s 1 c 2 h 0 I 2 s 2 h 0 I 2 c 1 c 2 h 0 I 2 s 1 c 2 h 0 I 2 s 2 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 8 > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > : f xI f yI f zI f xI + f yI + f zI + 9 > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > ; (7.50) 94 where f xI and f yI are the forces in the horizontal and vertical direction on master node I, m I is the moment on master node I, s 1 = sin 1 , s 2 = sin 2 , c 1 = cos 1 , and c 2 = cos 2 . Most methods used to date perform all calculations at the master nodes [1,44]. By doing this an extra equation is needed to calculate the internal moments in the system. Rather, by performing the calculations at the slave nodes, Eq. (7.50) can be used to describe all forces and moments at the master nodes. This is the same method used for the CB beam elements above [8]. 7.2.2.2 Updated Lagrangian Shell Formulation To use the Updated Lagrangian method, it is necessary to use the laminar components of stress and strain. To do this a corotational formulation must be implemented. Thus, corotational unit vectors ^ e x , ^ e y and ^ e z are introduced. Any number of base vectors can be chosen, the only constraint is that ^ e y must be normal to the plane. From here, either choose ^ e x tangent to the =constant lines or ^ e z perpendicular to ^ e y and i. The latter case will be used here for simplicity sake. The bases are then: ^ e y = x ; x ; kx ; x ; k (7.51a) ^ e z = i ^ e y ki ^ e y k (7.51b) ^ e x = ^ e y ^ e z k^ e y ^ e z k (7.51c) 95 where x ; = x ; y ; z ; T (7.52a) x ; = X I x I N I ; (;;) (7.52b) y ; = X I y I N I ; (;;) (7.52c) z ; = X I z I N I ; (;;) (7.52d) The constitutive model in corotational system will then be of the form ^ =F (^ ") (7.53) where the corotational strain (^ ") is given by ^ " =R T lam "R lam where R lam = ^ e x ^ e y ^ e z (7.54) However, based on the continuum equations formulated above, f xI f yI f zI int = n Q 1 X Q 1 =1 n Q 2 X Q 2 =1 n Q 3 X Q 3 =1 N I ;x N I ;y N I ;z 2 6 6 6 6 6 6 4 xx xy xz xy yy yz xz yz zz 3 7 7 7 7 7 7 5 J w Q 1 w Q 2 w Q 3 (7.55) 96 This needs to be converted to the corotational coordinates. Doing this, the equation for internal forces is obtained: f xI f yI f zI int = n Q 1 X Q 1 =1 n Q 2 X Q 2 =1 n Q 3 X Q 3 =1 N I ;^ x N I ;^ y N I ;^ z 2 6 6 6 6 6 6 4 ^ xx ^ xy ^ xz ^ xy 0 ^ yz ^ xz ^ yz ^ zz 3 7 7 7 7 7 7 5 2 6 6 6 6 6 6 4 R x^ x R y^ x R z^ x R x^ y R y^ y R z^ y R x^ z R y^ z R z^ z 3 7 7 7 7 7 7 5 J w Q 1 w Q 2 w Q 3 (7.56) Plane stress must also be implemented by setting ^ yy = 0. 7.2.2.3 Total Lagrangian Shell Formulation The fact that the Total Lagrangian method uses a xed reference frame means that a corotational formulation is unnecessary. As such, the formulation is much simpler. The constitutive equation will be of the form S =G (E) (7.57) 97 Based o the previous Total Lagrangian formulation, the internal forces are given by: f xI f yI f zI int = n Q 1 X Q 1 =1 n Q 2 X Q 2 =1 n Q 3 X Q 3 =1 N I ;X N I ;Y N I ;Z 2 6 6 6 6 6 6 4 S xx S xy S xz S xy 0 S yz S xz S yz S zz 3 7 7 7 7 7 7 5 0 B B B B B B B @ 2 6 6 6 6 6 6 4 1 0 0 0 1 0 0 0 1 3 7 7 7 7 7 7 5 2 6 6 6 6 6 6 4 u x;X u y;X u z;X u x;Y u y;Y u z;Y u x;Z u y;Z u z;Z 3 7 7 7 7 7 7 5 T 1 C C C C C C C A J 0 w Q 1 w Q 2 w Q 3 (7.58) 7.2.2.4 Reduced Integration Similar to the CB beam elements, general reduced integration techniques are insucient to describe thin and thick shells. This, however, is simply overcome by reducing the number of integration points used [63]. The table below shows how many Gauss points should be used for dierent elements: 98 Table 7.1: Reduced Integration in Shells Lamina Bilinear Biquadratic Bicubic shape (4 nodes) (8 or 9 nodes) (12 or 16 nodes) functions Normal Gaussian 2 x 2 3 x 3 4 x 4 rule Reduced Gaussian 1 x 1 2 x 2 3 x 3 rule 99 Chapter 8 Constitutive Model for Finite Element Analysis 8.1 Shape Memory Polymer Constitutive Model The material model for a shape memory polymer describing the stress tensor () versus strain tensor (") is given by [35]: = ( f S i + (1 f )S e ) 1 : ("" s " T ) (8.1) where f is the frozen fraction which is only dependent on temperature;" s re ects the strain storage during cooling; S e is the inverse of modulus of the entropic deformation tensor (E e ) and S i the inverse of modulus of internal energetic deformation tensor (E i ) which are usually found experimentally; The thermal strain" T can be expressed as " T = Z T T h ()dI (8.2) 100 is the coecient of thermal expansion; andT h is a reference temperature which is higher than T g . Write C = ( f S i + (1 f )S e ) 1 (8.3) such that =C : ("" s " T ) (8.4) The strain storage can be found from the equation: d" s dT =S e (T ) : (T ) d f dT =S e (T ) :C (T ) : ("" s " T ) d f dT (8.5) Since the stresses and strains are symmetric, Eq's (8.4) and (8.5) can be simplied to: fg = [C] (f"gf" s gf" T g) (8.6a) df" s g dT = [S e (T )] [C (T )] (f"gf" s gf" T g) d f dT (8.6b) where the material parameter matrices are given by [C] = ( f [S i ] + (1 f ) [S e ]) 1 (8.7a) [S i ] = 1 E i [I] (8.7b) [S e ] = 1 E e [I] (8.7c) 101 The frozen fraction, based on experimental results, can be given by f = 1 1 1 +c f (T h T ) n (8.8) The modulus of entropic deformation can be given by E e = 3NkT (8.9) where N is the cross-link density and k is Boltzmann's constant. The modulus of internal energetic deformation (E i ) is constant. 8.2 Composite Model The composite material will be made up of a SMP matrix and some unidirectional ber, of which the actual material can be chosen as desired. Only the unidirectional lamina method will be highlighted here, as this is the simplest. This involves a laminate made up of unidirectional laminae, whereby orientation angle of each lamina can be chosen. This can easily be extended to a fabric, whereby a laminate is made up of laminae of weaved bers set in the SMP matrix. Usually the formulation is done using plate theory, however since the continuum model above requires a height component, the derivation will be done in 3D. 102 8.2.1 Unidirectional Laminae Figure 8.1 shows a composite material loaded in the 1 and 2 principle directions. The material loaded in the 3 principle direction would be identical to Fig. 8.1(b) with 2 replaced by 3. This is used to nd the modulus of elasticity in each principle direction. 1 2 σ 1 σ 1 Matrix Fiber (a) 1 2 σ 2 σ 2 W (b) Figure 8.1: (a) Composite material loaded in the 1 principle direction, (b) Composite material loaded in the 2 principle direction The constitutive equations for the ber and matrix, denoted by subscript f and m respectively are given by: f =E f " f (8.10a) m =E m " m " s Z T T h dT (8.10b) It is assume that the stored strain is a constant and not a function of " m since it is found from previous knowledge. Using the fact that the strain is constant in the 1 principle 103 direction, i.e. " f = " m = " 2 , and the stress is constant in the 2 and 3 directions, i.e. f = m = 2 , it is found that the combined moduli in each direction are [4]: 1 =E 1 " 1 E s 1 " s + Z T T h dT (8.11a) 2 =E 2 " 2 E s 2 " s + Z T T h dT (8.11b) 3 =E 3 " 3 E s 3 " s + Z T T h dT (8.11c) with E 1 =E f V f +E m V m (8.12a) E s 1 =E m V m (8.12b) E 2 =E 3 = E f E m E f V m +E m V f (8.12c) E s 2 =E s 3 = E f E m V m E f V m +E m V f (8.12d) and V f +V m = 1. Similarly, 23 =G 23 23 G s 23 s (8.13a) 12 =G 13 13 G s 13 s (8.13b) 12 =G 12 12 G s 12 s (8.13c) 104 with G 13 =G 12 = G f G m G f V m +G m V f (8.14a) G s 13 =G s 12 = G f G m V m G f V m +G m V f (8.14b) G 23 =G m V f + 4 V m 4 (1V f ) +V f G m =G f (8.14c) G s 23 =G m V m V f + 4 V m 4 V m +V f G m =G f (8.14d) 4 = 3 4 m +Gm=Gf 4 (1 m ) (8.14e) These are now the constitutive equations for when the composite is layed at 0 or 90 . It would be useful to nd the constitutive equations for dierent layup orientations. In order to transform from principle coordinates to Cartesian coordinates, the lamina will be assessed as a plate. The coordinate systems can be seen in Fig. 8.2, where 1 and 2 are the principle axes and s and z are the local coordinates. z 1 2 s γ Figure 8.2: Coordinates of a Lamina The constitutive equations in cartesian coordinates are obtained by multiplying the stress and strains in the principle direction by a rotation matrix, thus obtaining [4]: fg = Q f"g Q s f" s g + " T (8.15) 105 where Q = [T ( )] 1 [Q] [R] [T ( )] [R] 1 (8.16a) Q s = [T ( )] 1 [Q s ] [R] [T ( )] [R] 1 (8.16b) with the vectors are given by: fg = s ^ y z ^ yz sz s^ y T (8.17a) f"g = " s " ^ y " z ^ yz sz s^ y T (8.17b) f" s g = " s s " s ^ y " s z s ^ yz s sz s s^ y T (8.17c) " T = Z T T h ()d 1 1 1 0 0 0 T (8.17d) 106 and the matrices by [Q] = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 E 1 1 12 21 0 E 1 21 1 12 21 0 0 0 0 E 3 0 0 0 0 E 2 12 1 12 21 0 E 2 1 12 21 0 0 0 0 0 0 KG 23 0 0 0 0 0 0 G 12 0 0 0 0 0 0 KG 13 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; 21 = E 2 E 1 12 (8.18a) [Q s ] = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 E s 1 1 12 21 0 E s 1 21 1 12 21 0 0 0 0 E s 3 0 0 0 0 E s 2 12 1 12 21 0 E s 2 1 12 21 0 0 0 0 0 0 KG s 23 0 0 0 0 0 0 G s 12 0 0 0 0 0 0 KG s 13 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; 21 = E s 2 E s 1 12 (8.18b) [T ( )] = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 cos 2 0 sin 2 0 2 sin cos 0 0 1 0 0 0 0 sin 2 0 cos 2 0 2 sin cos 0 0 0 0 cos 0 sin sin cos 0 sin cos 0 cos 2 sin 2 0 0 0 0 sin 0 cos 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (8.18c) 107 [R] = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 2 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (8.18d) The stored strain can be obtained from the matrix form [35] of Eq. (8.5) as follows: df" s g k dT = [E e ] 1 f m g k d f dT = [E e ] 1 Q 0 f"g k Q 0s f" s g k + Z T T h dT k d f dT (8.19) 108 where, using a similar process as for Q and with Eq.'s (8.10b) and (8.13c) Q 0 = [T ( )] 1 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 Em 1 12 21 0 Em 21 1 12 21 0 0 0 0 E 3 0 0 0 0 E 2 12 1 12 21 0 E 2 1 12 21 0 0 0 0 0 0 KG 23 0 0 0 0 0 0 G 12 0 0 0 0 0 0 KG 13 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 [R] [T ( )] [R] 1 (8.20a) Q 0s = [T ( )] 1 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 E s m 1 12 21 0 E s m 21 1 12 21 0 0 0 0 E s 3 0 0 0 0 E s 2 12 1 12 21 0 E s 2 1 12 21 0 0 0 0 0 0 KG s 23 0 0 0 0 0 0 G s 12 0 0 0 0 0 0 KG s 13 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 [R] [T ( )] [R] 1 (8.20b) here 21 = 12 E 2 =E m , 21 = 12 E s 2 =E s m and K is a shear correction coecient, which is set to 5/6. 109 8.2.2 Multilayer Laminate The calculations for internal forces must be done at each layer independently. To do this, the mid-plane of the laminate must be extrapolated to the mid-plane of each layer. This is done by assuming that the cross-section does not deform: x Ik =x I + h 0 I 2 1 1 N lay 2 N k 1 N lay cos 1 sin 2 (8.21a) y Ik =y I + h 0 I 2 1 1 N lay 2 N k 1 N lay sin 1 sin 2 (8.21b) z Ik =z I + h 0 I 2 1 1 N lay 2 N k 1 N lay cos 2 (8.21c) whereN lay is the number of layers in the laminate, andN k is thek'th layer with numbering beginning at the top surface, as seen in Fig 8.3. The amount of rotation is the same at all points along the height of the cross-section. For the beam element, let 2 = 90 . 1 2 k-1 N k . . . ŷ N ŷ k ŷ 2 ŷ 1 ŷ s Figure 8.3: Layers in the laminate 110 The forces at the mid-plane can then be calculated by f int I = X N k f int Ik (8.22a) m int zI = X N k m int zIk + X N k h 0 I 2 1 1 N lay 2 N k 1 N lay sin 1 sin 2 f int xIk + cos 1 sin 2 f int yIk (8.22b) m int xI = X N k m int xIk + X N k h 0 I 2 1 1 N lay 2 N k 1 N lay cos 1 cos 2 f int xIk + sin 1 cos 2 f int yIk sin 2 f int zIk (8.22c) where the subscript k denotes the k'th layer. This method also works for beams and shells with non-rectangular cross sections, such as an I-beam. 111 Chapter 9 Finite Element Solution Methods Since the Finite Element Method is inherently a numerical analysis, it is important to rst describe the solution method before showing the structure at the dierent phases of the SMPC. There are numerous ways in which the FE equations can be solved but only two will be outlined here. The rst is a more direct approach whereas the second requires some extra derivation. 9.1 Direct Solution from Force Vector This method takes the equations directly as derived in the previous sections and uses them to solve the equilibrium equation f int (u)f ext =0 (9.1) where f int is the internal force vector which is a function of the unknown dispalcement vector u; and f ext is the external force vector, found from the applied loads. 112 The internal force vector f int (u) is built up element by element by directly using the Updated or Total Lagrangian equation for a beam or shell (Eq's (7.41), (7.43), (7.56) or (7.58)) along with the stress-strain relation (8.6). The major problem with this method is the diculty to solve Eq. (9.1). The simplest way to do this is to use Newton's method. This, however, is not robust since it needs an accurate initial guess; and is very slow as the Jacobian needs to be calculated at every iteration. This can be done numerically with a nite dierence algorithm. A quasi-Newton method, such as BFGS [12, 21, 22, 49], can be used to improve the speed of the solver by updating the inverted Jacobian with A k+1 = I y k u T k y T k u k T A k I y k u T k y T k u k + u k u T k y T k u k (9.2) whereA is the inverse Jacobian,y k =f k+1 f k ,k is the iteration number, and u comes from the Newton's method displacement updating step u k+1 = u k + u k . While this method is much faster than regular Newton's method, it still requires a good initial guess for the displacement and for the Jacobian. A lot of work has been done to improve on these methods. The best methods to do this involve the use of optimization algorithms. These work by minimizing f = 1 2 f T f (9.3) instead of solving Eq. (9.1). This is done by nding the gradiantrf and then either using the same Newton or quasi-Newton methods described above (replacing the Jacobian 113 with the Hessian) [59]; a steepest descent method, such as conjugate gradient [16]; or a combination of these [53,60]. While optimization methods are powerful, they are dicult to implement and can also be computationally expensive if the error gradientrf and the Hessian are found numerically. Instead, a common method used in Finite Elements is to linearize the equations. The derivation follows. 9.2 Linearization of FE Equations Linearization is commonly used in nonliear Finite Element analysis. This is because the equations are easier to solve and the method is much faster. Another advantage of this method is that the rst iteration yields the linear solution if the initial displacement guess is set to zero. This means that a good initial guess is not needed if the force increments through the linear range. Further derivation is needed before this can be implemented. This is done by substituting the displacement updating step and the incremental stress into all of the equations. This will only be done for the Total Lagrangian shell but can similarly be implemented into the Updated Lagrangian beam or shell. Namely, on the continuum level, substitute u k+1 iJ =u k iJ + u k ij (9.4a) S k+1 km =S k km + S k km (9.4b) into Eq's (7.24), (8.15) and (8.19). Then on the shell element level, substitute u k+1 =u k + u k (9.5) 114 into Eq's (7.50) and (8.21). Where u = u x u y u z 1 2 (9.6a) u = u x u y u z 1 2 (9.6b) for each unconstrained node. Superscript k refers to the iteration number. 9.2.1 Element Stiness Matrices The internal force equation (7.24) for each element becomes, ignoring higher order terms f int(k+1) iI = Z 0 dN I dX k S k km im +u k iJ dN J dX m d 0 + Z 0 dN I dX k S k km im +u k iJ @N J @X m d 0 + Z 0 dN I dX k S k km dN J dX m d 0 u k iJ + Z 0 dN I dX k S k km dN J dX m d 0 u k iJ (9.7) Using the strain-displacement relation, the stress increment is S k km =C kmln dN J dX l pn +u k pK dN K dX n u k pJ (9.8) After substituting this into Eq. (9.7), assuming no body forces applied, the internal force is f int(k+1) iI = G k ipIJ +K k ipIJ u k pJ +f int(k) iI (9.9) 115 where f int(k) iI = Z 0 dN I dX k S k km im +u k iJ dN J dX m d 0 (9.10a) K k ipIJ = Z 0 dN I dX k C kmln dN J dX l pn +u k pK dN K dX n im +u k iL dN L dX m d 0 (9.10b) G k ipIJ = Z 0 dN I dX k S k km dN J dX m pi d 0 (9.10c) Equation (9.9) can be further simplied into matrix form to obtain ffg int(k+1) = G k + K k fug k +ffg int(k) (9.11) where K and G are the material and geometric tangent stiness matrices respectively. Here the incremental strain equation is used to nd K " kl = 1 2 dN J dX l kp + @u p @X k + dN J dX k pl + @u p @X l u pJ (9.12) From this the vector form is f"g = [A] [B]fug = [D]fug (9.13) 116 where [A] = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 +u 1;1 0 0 u 2;1 0 0 0 u 1;2 0 0 1 +u 2;2 0 0 0 u 1;3 0 0 u 2;3 0 u 1;3 u 1;2 0 u 2;3 u 2;2 + 1 u 1;3 0 u 1;1 + 1 u 2;3 0 u 2;1 u 1;2 u 1;1 + 1 0 u 2;2 + 1 u 2;1 0 ::: ::: u 3;1 0 0 0 u 3;2 0 0 0 1 +u 3;3 0 u 3;3 + 1 u 3;2 u 3;3 + 1 0 u 3;1 u 3;2 u 3;1 0 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (9.14a) [B] = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 N 1;1 0 0 N 2;1 0 0 ::: N 16;1 0 0 N 1;1 0 0 N 2;1 0 0 ::: N 16;1 0 0 N 1;1 0 0 N 2;1 0 0 ::: N 16;1 0 0 0 N 1;1 0 0 N 2;1 0 ::: 0 N 16;1 0 0 N 1;1 0 0 N 2;1 0 ::: 0 N 16;1 0 0 N 1;1 0 0 N 2;1 0 ::: 0 N 16;1 0 0 0 N 1;1 0 0 N 2;1 ::: 0 0 N 16;1 0 0 N 1;1 0 0 N 2;1 ::: 0 0 N 16;1 0 0 N 1;1 0 0 N 2;1 ::: 0 0 N 16;1 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (9.14b) 117 with ( ) i;j =d( ) i =dX j . Now the material tangent stiness matrix K is given by K = Z [D] T [C] [D]dxdydz = Z 1 1 Z 1 1 Z 1 1 [D] T [C] [D]J ddd (9.15) The geometric tangent stiness matrix G is given as G = Z T T [G] T dxdydz = Z 1 1 Z 1 1 Z 1 1 T T [G] T J ddd (9.16) where T is a transformation matrix to convert the unknown incremental displacement vector from ordered by direction to ordered by node, i.e. from u 11 u 12 ::: u 1(16) j ::: j u 31 ::: u 3(16) T (9.17a) to u 11 u 21 u 31 j ::: j u 1(16) u 2(16) u 3(16) T (9.17b) and the unsorted geometric stiness matrix is [G] = 2 6 6 6 6 6 6 4 h ^ G i 0 0 0 h ^ G i 0 0 0 h ^ G i 3 7 7 7 7 7 7 5 (9.18) with h ^ G i = dN dX [S] dN dX T (9.19) 118 dN dX is the matrix of derivatives of the shape functions and [S] is the stress matrix. 9.2.2 Material Tangent Stiness with Stored Strain The incremented stored strain (" s + " s ) is derived in the same way by substituting the incremented strain equation (9.13) into Eq. (8.19). It then becomes df" s g dT = [E e ] 1 Q 0 [D]fug Q 0s f" s g d f dT (9.20) now assume that the incremtal stored strain " s is a function of u such that " s = [D s ] u. Now d [D s ] dT fug = [E e ] 1 Q 0 [D] Q 0s [D s ] d f dT fug (9.21) Using Eq's (8.15), (9.15) and (9.21) the material stiness matrix with stored strain becomes K = Z 1 1 Z 1 1 Z 1 1 [D] T Q [D] Q s [D s ] J ddd (9.22) 119 9.2.3 Element Stiness Matrix for Shells Equations (9.11), (9.16) and (9.22) are internal forces and corresponding tangent stiness matrices for a continuum element. To obtain the stiness matrix for a shell, substitute Eq. (9.5) into (7.50) ffg mast = 8 > > > > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > > > > : . . . f xI f yI f zI m zI m xI . . . 9 > > > > > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > > > > > ; mast = T T + T T G + K Tfug +ffg int (9.23) where T and T are matrices made up of T I (given in Eq. (7.50)) and T I along the diagonals respectively, with T I = h 0 I 2 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0 0 0 (c 1 s 2 1 +s 1 c 2 2 ) (s 1 c 2 1 +c 1 s 2 2 ) 0 0 0 (s 1 s 2 1 c 1 c 2 2 ) (c 1 c 2 1 s 1 s 2 2 ) 0 0 0 0 c 2 2 0 0 0 (c 1 s 2 1 +s 1 c 2 2 ) (s 1 c 2 1 +c 1 s 2 2 ) 0 0 0 (s 1 s 2 1 c 1 c 2 2 ) (c 1 c 2 1 s 1 s 2 2 ) 0 0 0 0 c 2 2 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (9.24) 120 Note that K , G andffg int in Eq. (9.23) must be reorder by bottom then top node in the order that the shell nodes (subscript I) are numbered. To get Eq. (9.23) into a useable form rearrange and ignore higher order terms such that ffg mast = T T G + K T + T T f int fug +T T ffg int (9.25) where the subscript I refers to the master node number and T I f int = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 T 1 T 2 0 0 0 T 2 T 4 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (9.26) with T 1 = h 0 I 2 h f xI f + xI c 1 s 2 + f yI f + yI s 1 s 2 i (9.27a) T 2 = h 0 I 2 h f xI f + xI s 1 c 2 f yI f + yI c 1 c 2 i (9.27b) T 4 = h 0 I 2 h f xI f + xI c 1 s 2 + f yI f + yI s 1 s 2 + f zI f + zI c 2 i (9.27c) 121 9.2.4 Global Stiness Matrix with Multilayering Finally, in order to stitch the global stiness matrix together for a multilayered laminate, use Eq. (9.5) along with (8.22). ffg mid = T T f lay fug +T T ffg lay (9.28) where, from (8.22) T I = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 0 0 hs 1 s 2 hc 1 c 2 0 1 0 hc 1 s 2 hs 1 c 2 0 0 1 0 hs 2 0 0 0 1 0 0 0 0 0 1 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (9.29a) T I f lay = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 T A T B 0 0 0 T B T D 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (9.29b) 122 with h = h 0 I 2 1 1 N lay 2 N k 1 N lay (9.30a) T A = h h f lay xI c 1 s 2 f lay yI s 1 s 2 i (9.30b) T B = h h f lay xI s 1 c 2 +f lay yI c 1 c 2 i (9.30c) T D = h h f lay xI c 1 s 2 f lay yI s 1 s 2 f lay zI c 2 i (9.30d) Finally, the total global stiness matrix [K] and global internal force vectorffg int are built up of the midlayer matrices and vectors which are given by [K] mid =T T G lay + K lay T + T T f lay (9.31a) ffg int;mid =T T ffg int;lay (9.31b) such that the global equilibrium equation is [K (u)]fug =ffg ext ff (u)g int (9.32) where [K (u)] is the global stiness matrix;fug is a vector of unknown incremental displacements;ff (u)g int is the internal force vector of known diplacements from the previous iteration; andffg ext is the external force vector. 123 9.3 Dynamic Equations The linearization method is very helpful for modeling a static or quasi-static structure. However, for shape recovery, it is necessary to have dynamic equations. Equation (9.1) is easily extended to allow for dynamics. The equation of motion for the structure is then given by M u +f int (u) =f ext (9.33) where M is the global mass matrix and ( ) denotes a double derivative in time (i.e. acceleration). This equation can either be solved explicitly or implicitly. 9.3.1 Explicit Dynamic Solvers The explicit solvers are much easier to implement than their implicit counterparts, however the implicit solvers are more robust and accurate. The two most popular explicit solvers are the nite dierence method and the Runge-Kutte method. The nite dierence method with a central dierence discretization is given by [M] u i+1 2u i +u i1 t 2 +ff (u)g int =ffg ext (9.34) wherei represents the increment in time and t is the chosen time step. The Runge-Kutte method solves the system 8 > > < > > : _ y 1 _ y 2 9 > > = > > ; = 8 > > < > > : y 2 [M] 1 ffg ext ff (y 1 )g int 9 > > = > > ; (9.35) 124 as per the equations in Section 2.3. To start both of these methods, the initial displacements and velocities are needed. The solution is then stepped through time until a desired stopping condition is met [8]. 9.3.2 Implicit Dynamic Solvers If the explicit solver is not adequate enough to solve the dynamic equations, due to stability or stiness issues, then an implicit solver is needed. There are various implicit solvers but the most popular is the Newton's Method. This method uses the equations outlined in Section 2.2.2 along with the linearized model. Implicit methods are inherently iterative and as such are slower than explicit methods [8]. This section, however, was deemed to be beyond the scope of this research. 9.3.3 Mass Matrix Generally for isoparametic elements the elemental mass matrix is found by [5] h ^ M i = Z v dN dX T dN dX dV (9.36) where is the density. This can then be converted to a beam or shell using the same transformation matrices as above. However, once the global mass matrix is assembled, it will not be a diagonal matrix. 125 To obtain a diagonal matrix instead, the lumped mass approach can be used [8,52]. This will produce a mass matrix of the form [M i ] = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 m i 0 0 0 0 0 m i 0 0 0 0 0 m i 0 0 0 0 0 I i 0 0 0 0 0 I i 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (9.37) herei denotes the node number. The advantage with this method is that the mass matrix is independant of the shape of the structure. This leads to a faster computational time. For a composite, the density is = f V f + m (1V f ) (9.38) where subscript f and m denote the ber and matrix respectively. 126 Chapter 10 Shape Fixation and Shape Recovery for FE Model In this section, the physical behaviors of the SMPC structure in shape xation and shape recovery are examined for the Finite Element (FE) Model. The linearized shell element with the Total Lagrangian method will be used. 10.1 Cooling for Shape Fixation In shape xation, the structure is rst loaded incrementally at a high temperature T h ; the structure then remains under xed loading while the temperature is decreased to a low temperature T l , where T h >T g >T l . After the temperature reaches T l , all loads will be removed. The cooling of the SMPC structure (from T l to T h ) is assumed to be a slow process such that the quasi-static model given in Section 9.2 is sucient to describe the de ection. The internal force vector and the global stiness matrix are built up element by element and layer by layer. The equations at the single layer element level are described at the dierent stages of shape xation as follows. 127 (i) At T =T h , external loads are applied. The de ection is governed by [K (u)]fug =ffg ext ff (u)g int (10.1) with tangent stiness matrix and stress vector given by K = Z 1 1 Z 1 1 Z 1 1 [D] T h Q h i [D]J ddd (10.2a) fg = h Q h i f"g (10.2b) (ii) From T h to T l , the external loads are held constantly as temperature decreases. The de ection is governed by [K (u)]fug =ffg ext ff (u)g int (10.3) with tangent stiness matrix and stress vector given by K = Z 1 1 Z 1 1 Z 1 1 [D] T Q [D] Q s [D s ] J ddd (10.4a) fg = Q f"g Q s f" s g + " T (10.4b) (iii) At T =T l , all forces are removed and the de ection is governed by [K (u)]fug =ff (u)g int (10.5) 128 with tangent stiness matrix and stress vector given by K = Z 1 1 Z 1 1 Z 1 1 [D] T h Q l i [D] h Q ls i [D s ] J ddd (10.6a) fg = h Q l i f"g h Q ls i f" s g (10.6b) 10.2 Heating for Shape Recovery In this process, without loss of generality, assume that the structure is heated from T l to T h in a linear fashion T =T l +t; t 0 (10.7) The heating process is intrinsically a dynamic one. This means that the quasi-static model needs to be modied with an added inertia term. Thus, at the beginning of shape recovery, the equation of motion is given by [9]: [M]f ug +ff (u)g int = 0 (10.8) for 0 t t 1 , with t 1 = (T h T l )=, where u = u(s;t), with t being the temporal parameter. The initial conditions are fu (s; 0)g =fu l0 (s)g; f _ u (s; 0)g = 0 (10.9) 129 10.3 Vibration after Shape Recovery Once the shape recovery is completed, there is a nal process in which the structure is cooled down fromT h toT l and no loads are applied during this time period. The equation of motion after the cooling process is completed can be written as [M]f ug +ff (u)g int =ffg ext (10.10) which is subject to initial conditions fu (s; 0)g = 0; f _ u (s; 0)g = 0 (10.11) 130 Chapter 11 Simulation Results The proposed quasi-static model of SMPC beams and numerical solution method are illustrated on the cantilever beam under a point load at the tip; see Fig. 11.1. Simulations of the beam with other loadings (as shown in Fig. 6.1) have also been performed. Due to limited space, only the tip load case is presented in this paper. The values of the parameters used in the simulation are given in Table 11.1, where the values of Young's moduli are from Reference 35. The temperature range was taken as 1:05T g >T g > 0:8T g . s y x P y P y Δ w Figure 11.1: Large de ection cantilever beam with a point load applied to the end 131 Table 11.1: Values of the beam parameters used in the simulation Variable Value E l 750 MPa E h 8:8 MPa E f 130x10 3 MPa T g 334 K b 10 mm h 0.25 mm L 100 mm P 10(EI) h =L 2 N m 101 Load steps at T h 300 Table 11.2: Values used to calculate f and E e Variable Value k 1:3086x10 17 mm 2 kg/s 2 K c f 2:76x10 5 n 4 132 11.1 Loading at T h To start the process, the beam is loaded incrementally at a high temperature (T h ). The computed results are shown in Figs. 11.2 and 11.3. Figure 11.2 plots the normalized displacement versus the normalized load. The end-shortening and linear results are also included in the gure. The numerical results are in good agreement with the analytical results and those results obtained using the nite element software ABAQUS. Figure 11.3 depicts the de ection shapes of the beam for every tenth load step. (There are 300 loading steps; see Table 11.1.) Since it was chosen that the load applied is a function of EI=L 2 at high temperature, the normalized shape of the beam does not change regardless of the composite layup or ber volume fraction, as seen in Fig.'s 11.2 and 11.3. 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 9 10 w/L or Δ/L PL 2 /(E h I) w Δ linear Figure 11.2: Plot of normalized load versus displacement 133 0 0.2 0.4 0.6 0.8 1 −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 y/L x/L Figure 11.3: Plot of increasing normalized de ection of the beam as load increases 11.2 Cooling from T h to T l with load applied The beam is now cooled from T h toT l with the tip load held constant. This is seen by the dash-dotted curve in Fig.'s 11.4 to 11.8 below. 11.3 Load removal at T l The beam de ection at T l after the load is removed is plotted in Fig.'s 11.4 to 11.8; see the solid curve. For comparison, the de ection of the beam with the load applied at T h is also included in the gures (dashed curve). To see the eect of making the SMP a composite, a plot of the results of the pure SMP beam was included in Fig. 11.4. 134 0 0.2 0.4 0.6 0.8 1 −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded Figure 11.4: Plot of normalized de ection of a pure SMP beam along its length comparing the cases of when the load is applied and when it is removed It should be noted here that Fig.'s 11.2 and 11.3 are valid for all composite and pure SMP cases, while Fig. 11.4 only relates to the pure SMP case (i.e. V f = 0). The composite model in this work allows for the variance of the ber volume fraction (V f ) of each layer, the ber layup angle ( ) of each layer, and the number of layers in the laminate. Figure 11.5 shows the results while varying the ber volume fraction. The lamina seems to react closer to the behavior of the ber due to the fact that the Moduli of carbon is so much higher than that of the SMP. This does, however, increase eective modulus of the composite substantially. There is good agreement between the nite dierence model (results on left) and that Finte Element model (results on right), except for V f = 1%. This is because as the 2D eects become important in this range. Figure 11.6 shows the results for a varying layup angle. When the bers are laid in the axial direction ( = 0 ), the lamina takes on the behavior much closer to that of the ber, 135 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (a) V f =1% (FD) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded Unloaded Linear (b) V f =10% (FD) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded Unloaded Linear (c) V f =50% (FD) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (d) V f =1% (FEM) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (e) V f =10% (FEM) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (f) V f =50% (FEM) Figure 11.5: Plot of normalized de ection of SMP composite beam along its length comparing the cases of when the load is applied and when it is removed with 1 layer, 0 layup angle and varying ber volume fraction 136 such that when loads are applied the shape of the beam does not change at all through the temperate range, but when the loads are removed the beam almost goes back to its original shape. This result was validated by solving the linear equation of a beam with the assumption that the only force applied is a force due to the stored strain. Namely: D l bw 0000 = d 2 M int s dx 2 (11.1) When the bers are perpendicular the axis of the beam ( = 90 ), the lamina behaves almost exactly like the pure SMP beam. In this case, however, the eective modulus does not increase much. To maintain symmetry, a 2 laminae layup of [45 -45 ] is used instead of a single 45 lamina. As can be seen in Fig. 11.6(b), the beam does not maintain its shape at low temperature. This is an error caused by the fact that the beam model is used. The stored strain is not adequately described in one dimension when the bers are not at 0 or 90 . To overcome this, a plate model must be used. However, the FE model experiences shear locking at any layup other than at 0 or 90 (see Fig. 11.6(e)). It is unclear why this occurs and further investigation is needed. The models do however perform well at 0 and 90 and show good agreement. In order to gain an increased eective modulus while maintaining the SMP behavior, a multilayered composite can be used (laminate). An example of a laminate with 8 layers is given in Fig 11.7. This behavior can further be improved by including layers of pure SMP [62] (as in Fig. 11.8). 137 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded Unloaded Linear (a) =0 (FD) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (b) =45 -45 (FD) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (c) =90 (FD) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (d) =0 (FEM) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (e) =45 -45 (FEM) 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (f) =90 (FEM) Figure 11.6: Plot of normalized de ection of SMP composite beam along its length comparing the cases of when the load is applied and when it is removed with 1 layer (2 for [45 45 ] ), 10% ber volume fraction and varying layup angles 138 0 0.2 0.4 0.6 0.8 1 −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (a) FD 0 0.2 0.4 0.6 0.8 1 −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded (b) FEM Figure 11.7: Plot of normalized de ection of SMP composite beam along its length comparing the cases of when the load is applied and when it is removed with 8 layers, 60% ber volume fraction and [90 90 90 0 ] s layup 139 It is important to remember that the shape of the beam at high temperature remains the same in all of these cases due to dierent magnitudes of point load being applied to each. The eective modulus in each case will be dierent but will always be higher than that of the pure SMP. By this, the rigidity of the beam is increased but as seen here, the ability to maintain the desired shape may be sacriced. 0 0.2 0.4 0.6 0.8 1 −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 x/L y/L High Temp Low Temp Loaded Low Temp Unloaded Figure 11.8: Plot of normalized de ection of SMP composite beam along its length comparing the cases of when the load is applied and when it is removed with 8 layers, 50% ber volume fraction and [90 90 90 0 ] s layup 11.4 Shape Recovery No numerical simulations were performed for the heating cases (Section 5.3.2). It is expected from the theory that the beam will recover its shape completely as no histeritic eects are included [35]. In order to use the model proposed here to illustrate the recovery 140 behavior seen in the experiments, a dierent recovery material model is needed or some sort of hysteretic eect needs to be included in the model. 141 Chapter 12 Conclusions This dissertation presents a quasi-static model of a Shape Memory Polymer Composite (SMPC) cantilever beam, which has temperature-dependent material nonlinearity and geometric nonlinearity due to large de ection. An analytical beam as well as a Finite Element model was presented to describe the geometry. A material model by Liu et al [35] was extended to be used for SMPC beams by a rule-of-mixtures approach. The beam model obtained can be used to study shape xation and shape recovery processes of SMPC structures in aerospace applications. A unied method for systematic establishment of dynamic equations for nonlinear beams has been presented. In this method, local (Lagrangian) coordinates are used to describe the elastic deformation of a beam and Hamilton's principle is applied to obtain the equations of motion. This method leads to a full nonlinear beam model with three independent local coordinates (u, w, ), which is readily useable in nite element discretization. While this work is focused on nonlinear analysis, the full beam model is shown to be able to produce linear beam models, such as Timoshenko beam, Rayleigh beam, and Euler-Bernoulli beam. Although only nonlinear strain-displacement relations 142 are considered in this paper, the proposed method is applicable to beams with both geometric and material nonlinearities. With appropriate assumptions, the full model yields three beam models: (i) a nonlinear beam with small shear and extension; (ii) a pure bending nonlinear beam described by the vertical de ection w; and (iii) a pure bending nonlinear beam described by the rotation . These models are simple, easy to be analyzed, and convenient to use for numerical solutions. The -based beam model is a new nonlinear beam model that is robust in computation, and avoids singularities due to large deformation, as with some other beam models. Liu's constitutive model was extended to an SMP composite material. Plane stress assumptions where made to reduce the model to that of a beam to obtain an equivalent modulus. This was then implemented into the analytical beam model to describe the structure in shape xation and recovery. Numerical analysis was then performed. For this, nite dierence discretization and related solution algorithm are also derived based on the beam model. The numerical results obtained by the proposed method reveal the behaviors of SMPC structures in quasi-static states. Although a cantilever beam is considered in this work, the results obtained in the current investigation can be extended to SMPC beams with arbitrary boundary conditions. A full Finite Element review was performed for geometrically nonlinear beams and shells using both the Total and Updated Lagrangian methods. Related solution algorithms were also highlighted but the focus was left on the linearization of the geometrically nonlinear shell using the Total Lagrangian method. This was then combined with a 3D 143 version of the previously derived composite form of Liu's SMP constitutive equation. This could then be used to numerically analyse the beam through shape xation and recovery. It was observed from numerical simulations that the shape xation behavior of the beam is highly dependent on the ber volume fraction and the layup angle of each lamina in the composite. It is important that the bers increase the rigidity of the SMP but there will be some penalty imposed on the xation of the beam. The model is seen to perform well for a 0 and 90 layup but fails at all other ber angles due to the lack of 2-D eects in the beam and shear locking in the FE model. The recoverability behavior as seen in experiments was not imposed in this model but further suggestions were made as to how this could be done. 144 Chapter 13 Recommendations for Future Work This report is a precursor of an ongoing project on the development of modeling and analysis tools of SMP composite space structures. It is recommended that further study address the following issues: 1. Dynamic modeling of SMP beams Dynamic equations for the SMPC beam have been derived but no numerical results have been done. These numerical results should be found in order to gain an understanding of the recovery process. The results would be more valuable with the addition of a hysteretic model but this would only be possible with experimental results or a new recovery material model. 2. Finite deformation model of SMP composite beams The constitutive equation proposed by Liu [35] for small strains. While this yields fairly accurate results for most of the cases that the SMPC space structure may encounter, a nite deformation model should be implemented for generality sake. 145 This would allow for large strains in the model. A similar approach for the beam equations may be used but a new constitutive equation would be required. 3. Nanocomposite or fabric composite model This paper presents a unidirectional laminated composite model and suggests meth- ods for implementing more complicated composites. While the unidirectional lami- nate is easy to model, it is not useful for real world applications. For this reason the more complicated composites, such as nanocomposites or fabrics, should be modeled. To further improve the useability of the model, thin wall beams can be used. This is easily extended from the method shown here for composite beams. 4. Mechanical behaviors of SMP composite structures in dierent stages of shape xation and shape recovery It is not very useful to only have one beam for application purposes. Generally a structure made up of multiple beams would be needed. These beam structures need to be analysed through the shape xation and recovery process in order to gain a better understanding of their performance characteristics. 5. Solve the Locking Problem One of the major shortcomings of the models presented is the fact that the composite beams can only have a 0 or 90 layup. Any other symmetric layup would lock due to the large dierence between normal and shearing strains at high temperature. It is imperative that this be solved. 6. Include Transient Temperature Eects 146 It was assumed that the beam was heated uniformly for the models presented within. It would be interesting to investigate the eects of heating dierent parts of the beam individually and possibly using this to further control shape recovery. It would also be benecial to include other heating eects such as space temperature or heating during launch. 147 References [1] Sohrabuddin Ahmad, Bruce M. Irons, and O. C. Zienkiewicz. Analysis of thick and thin shell structures by curved nite elements. International Journal for Numerical Methods in Engineering, 2(3):419{451, 1970. [2] S. AL-Sadder and R.A.O. AL-Rawi. Finite dierence scheme for large-de ection analysis of non-prismatic cantilever beams subjected to dierent types of continuous and discontinuous loadings. Archive of Applied Mechanics, 75:459{473, 2006. [3] E.J. Barbero. Finite Element Analysis of Composite Materials. CRC Press, 2008. [4] E.J. Barbero. Introduction to Composite Material Design. CRC Press, Boca Raton, FL, 2nd edition, 2011. [5] K.J. Bathe. Finite Element Procedures in Engineering Analysis. Prentice Hall, 1996. [6] M. Batista and F. Kosel. Cantilever beam equilibrium congurations. International Journal of Solids and Structures, 42(16{17):4663 { 4672, 2005. [7] M. Behl and A. Lendlein. Shape-memory polymers. Materials Today, 10(4):20 { 28, 2007. [8] T. Belytschko, W.K. Liu, and B. Moran. Nonlinear Finite Elements for Continua and Structures. John Wiley and Sons, LTD, England, 2000. [9] H. Benaroya and M.L. Nagurka. Mechanical Vibrations: Analysis, Uncertainty, and Control. CRC Press, Boca Raton, FL, 3rd edition, 2010. [10] K.E. Bisshopp and D.C. Drucker. Large de ection of cantilever beams. Quarterly Applied Math, 3(3):272{275, 1945. [11] Kai-Uwe Bletzinger, Manfred Bischo, and Ekkehard Ramm. A unied approach for shear-locking-free triangular and rectangular shell nite elements. Computers & Structures, 75(3):321 { 334, 2000. [12] C.G. Broyden. The convergence of a class of double-rank minimization algorithms. Journal of the Institute of Mathematics and Its Applications, 6:76{90, 1970. [13] D. Campbell. Elastic memory composite material: An enabling technology for future furable space structures. In 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Austin, 2005. 148 [14] S.C. Chapra and R.P. Canale. Numerical Methods for Engineers. McGraw-Hill, New York, NY, 5th edition, 2006. [15] A. Charlesby. Atomic Radiation and Polymers. Pergamon Press, 1960. [16] Y.Y. Chen and S.Q. Du. Nonlinear conjugate gradient methods with wolfe type line search. Abstract and Applied Analysis, 2013, 2013. [17] Abel Cherouat and Houman Bourouchaki. Numerical tools for composite woven fabric preforming. Advances in Materials Science and Engineering, 2013. [18] R.D. Cook, D.S. Malkus, M.E. Plesha, and R.J. Witt. Concepts and Applications of Finite Element Analysis. John Wiley & Sons, Ltd, 4th edition, 2002. [19] L. Euler. Methodus Inveniendi Lineas Curvas. Lausann, Genev, apud Marcum- Michaelem Bousquet & socios, 1744. [20] D. Fertis. Nonlinear Mechanics. CRC Press, Boca Raton, FL, 2nd edition, 1998. [21] R. Flecther. A new approach to variable metric algorithms. Computer Journal, 13(3):317{322, 1970. [22] D. Goldfarb. A family of variable metric updates derived by variational means. Mathematics of Computation, 24(109):23{26, 1970. [23] W. R. Hamilton. On a general method in dynamics. Philosophical Transactions of the Royal Society of London, 124:pp. 247{308, 1834. [24] S.M. Han, H. Benaroya, and T. Wei. Dynamics of transversely vibrating beams using four engineering theories. Journal of Sound and Vibration, 225(5):935 { 988, 1999. [25] J.T. Holden. On the nite de ections of thin beams. International Journal of Solids and Structures, 8(8):1051 { 1055, 1972. [26] Thomas J.R. Hughes and Wing Kam Liu. Nonlinear nite element analysis of shells: Part i. three-dimensional shells. Computer Methods in Applied Mechanics and Engineering, 26(3):331 { 362, 1981. [27] Thomas J.R. Hughes and Wing Kam Liu. Nonlinear nite element analysis of shells-part ii. two-dimensional shells. Computer Methods in Applied Mechanics and Engineering, 27(2):167 { 181, 1981. [28] S. Krenk. Non-linear Modeling and Analysis of Solids and Structures. Cambridge University Press, Cambridge, UK, 2009. [29] X. Lan, Y. Liu, H. Lv, X. Wang, J. Leng, and S. Du. Fiber reinforced shape-memory polymer composite and its application in a deployable hinge. Smart Materials and Structures, 18(2):024002, 2009. [30] A. Lendlein and R. Langer. Biodegradable, elastic shape-memory polymers for potential biomedical applications. Science, 296:1673, 2002. 149 [31] Chung-Li Liao and J.N. Reddy. Analysis of anisotropic, stiened composite laminates using a continuum-based shell element. Computers & Structures, 34(6):805 { 815, 1990. [32] Chung-Li Liao and J.N. Reddy. Analysis of anisotropic, stiened composite laminates using a continuum-based shell element. Computers & Structures, 34(6):805 { 815, 1990. [33] A. Libai and J.G. Simmonds. The Nonlinear Theory of Elastic Shells. Cambridge University Press, 2nd edition, 1998. [34] J.K. Lin, G.H. Sapna, S.E. Scarborough, and B.C. Lopez. Advanced pre- cipitation radar antenna singly curved parabolic antenna re ector development. In AIAA-2003-1651, 44th AIAA/ASME/ ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference and Exhibit AIAA Gossamer Spacecraft Forum, Norfolk, VA, 2003. [35] Y. Liu, K. Gall, M.L. Dunn, A.R. Greenberg, and J. Diani. Thermomechanics of shape memory polymers: Uniaxial experiments and constitutive modeling. International Journal of Plasticity, 22(2):279 { 313, 2006. [36] Y. Liu, K. Gall, M.L. Dunn, and P. McCluskey. Thermomechanics of shape memory polymer nanocomposites. Mechanics of Materials, 36(10):929 { 940, 2004. [37] A.K. Mal and S.J. Singh. Deformation of Elasticity Solids. Prentice Hall, Englewood Clis, NJ, 1991. [38] K Mattiasson. Numerical results from large de ection beam and frame problems analysed by means of elliptic integrals. International Journal for Numerical Methods in Engineering, 17(1):145{153, 1981. [39] Q. Ni, C. Zhang, Y. Fu, G. Dai, and T. Kimura. Shape memory eect and mechanical properties of carbon nanotube/shape memory polymer nanocomposites. Composite Structures, 81(2):176 { 184, 2007. [40] O.O. Ochoa and J.N. Reddy. Finite Element Analysis of Composite Laminates. Kluwer Academic Publishers, 1992. [41] T. Ohki, Q. Ni, N. Ohsako, and M. Iwamoto. Mechanical and shape memory behavior of composites with shape memory polymer. Composites Part A: Applied Science and Manufacturing, 35(9):1065 { 1073, 2004. [42] P. Frank Pai. Geometrically exact beam theory without euler angles. International Journal of Solids and Structures, 48(21):3075 { 3090, 2011. [43] P.F. Pai and A.N. Palazotto. Large-deformation analysis of exible beams. International Journal of Solids and Structures, 33(9):1335 { 1353, 1996. [44] H. Parisch. Geometrical nonlinear analysis of shells. Computer Methods in Applied Mechanics and Engineering, 14(2):159 { 178, 1978. 150 [45] S.A Ragon, Z G urdal, and L.T Watson. A comparison of three algorithms for tracing nonlinear equilibrium paths of structural systems. International Journal of Solids and Structures, 39(3):689 { 698, 2002. [46] M.A. Rahman and A.K. Kowser. Nonlinear analysis of cantilever shape memory alloy beams of variable cross section. Smart Materials and Structures, 16:531{540, 2007. [47] D. Ratna and J. Karger-Kocsis. Recent advances in shape memory polymers and composites: a review. Journal of Materials Science, 43:254{269, 2008. [48] J.N. Reddy. Applied Functional Analysis and Variational Methods in Engineering. McGraw-Hill, New York, NY, 1986. [49] D.F Shanno. Conditioning of quasi-newton methods for function minimization. Mathematics of Computation, 24(111):647{656, 1970. [50] J.L. Curiel Sosa and A.J. Gil. Analysis of a continuum-based beam element in the framework of explicit-fem. Finite Elements in Analysis and Design, 45(8{9):583 { 591, 2009. [51] J.W. Strutt. Theory of Sound. Macmillan Publications Co., Inc., London, England, 1877. [52] Karan S. Surana. Lumped mass matrices with nonzero inertia for general shell and axisymmetric shell elements. International Journal for Numerical Methods in Engineering, 12(11):1635{1650, 1978. [53] S. Taheri and M. Mammadov. Solving systems of nonlinear equations using a globally convergent optimization algorithm. Global Journal of Technology and Optimization, 3:132{138, 2012. [54] S.P. Timoshenko. On the correction for shear of the dierential equation for transverse vibrations of bars of uniform cross-section. Philosophical Magazine, 744, 1921. [55] S.P. Timoshenko. History of Strength of Materials. Dover Publications Inc., New York, NY, 1953. [56] H. M. Wache, D. J. Tartakowska, A. Hentrich, and M. H. Wagner. Development of a polymer stent with shape memory eect as a drug delivery system. Journal of Materials Science: Materials in Medicine, 14:109{112, 2003. [57] K. Washizu. Variational Methods in Elasticity and Plasticity. Pergamon Press, Oxford, 2nd edition, 1974. [58] C. Willey, R. Bucolic, W. Skulney, Cadogan D.P., and J.K. Lin. The hybrid in- atable antenna. In AIAA 2001-1258, 42nd AIAA/ASME/ASCE/AHS/ASC SDM Conference, 2001. [59] M.A. Wolfe. Numerical methods for unconstrained optimization : an introduction. Van Nostrand Reinhold, New York, NY, 1978. 151 [60] Pu yan Nie. AnfSQPg approach with line search for a system of nonlinear equations. Mathematical and Computer Modelling, 43(3{4):368 { 373, 2006. [61] B. Yang, G.L. Davis, and H. Fang. Mathematical formulation for dynamic deployment of shape memory polymer composite space structure. In 11th AIAA Gossamer Spacecraft Forum, Orlando, 2010. [62] C. Zhang and Q. Ni. Bending behavior of shape memory polymer based laminates. Composite Structures, 78(2):153 { 161, 2007. [63] O. C. Zienkiewicz, R. L. Taylor, and J. M. Too. Reduced integration technique in general analysis of plates and shells. International Journal for Numerical Methods in Engineering, 3(2):275{290, 1971. [64] R. Zinno and E.J. Barbero. Total lagrangian formulation for laminated composite plates analysed by three-dimensional nite elements with two-dimensional kinematic constraints. Computers & Structures, 57(3):455 { 466, 1995. 152

Abstract (if available)

Abstract Shape Memory Polymer Composite (SMPC) structures are an ideal solution for large space structure deployment problems. This is due to their ability to be formed into a desired light‐weight loading shape and then transform back to their original shape by means of an applied stimulus. There is a wide array of constitutive and qualitative work being done on SMPC's but little or no development of dynamic equations. This dissertation documents a macroscopic model for the shape fixation and shape recovery processes of a SMPC cantilever beam. In particular the focus is on a quasi‐static model that can be used to describe the shape fixation process. Numerical results are obtained in this regard by use of finite difference approximation with Newton's method as well as a Finite Element Formulation. These models combine a nonlinear geometric model with a temperature dependent constitutive law. Additionally, the dynamic equations of the SMPC cantilever are derived. The Finite Element models are derived for a beam and a shell in order to keep the solution as general as possible. The shell element allows this to be expanded to thin or moderately thick structures of any shape. This study is a precursor to on-going work in modeling SMPC space structures.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Static and dynamic analysis of two-dimensional elastic continua by combined finite difference with distributed transfer function method

PDF

Optimal design, nonlinear analysis and shape control of deployable mesh reflectors

PDF

Dynamic analysis and control of one-dimensional distributed parameter systems

PDF

Modeling and dynamic analysis of coupled structure-moving subsystem problem

PDF

Homogenization procedures for the constitutive material modeling and analysis of aperiodic micro-structures

PDF

Multi-scale modeling of functionally graded materials (FGMs) using finite element methods

PDF

Terrain-following motion of an autonomous agent as means of motion planning in the unknown environment

PDF

Control of two-wheel mobile platform with application to power wheelchairs

PDF

Control of spacecraft with flexible structures using pulse-modulated thrusters

PDF

An approximation of frequency response of vibration system in mid-frequency area

PDF

An approach to experimentally based modeling and simulation of human motion

PDF

On the synthesis of controls for general nonlinear constrained mechanical systems

PDF

Dynamic modeling and simulations of rigid-flexible coupled systems using quaternion dynamics

PDF

On the synthesis of dynamics and control for complex multibody systems

PDF

Modeling and vibration analysis of wheelchair-users

PDF

Design, modeling and analysis of piezoelectric forceps actuator

PDF

The development of a mathematical model of a hybrid airship

PDF

Development of composite oriented strand board and structures

PDF

New approaches to satellite formation-keeping and the inverse problem of the calculus of variations

PDF

Modeling, analysis and experimental validation of flexible rotor systems with water-lubricated rubber bearings

Asset Metadata

Creator Bergman, Dean (author)

Core Title Mathematical model and numerical analysis of a shape memory polymer composite cantilever beam for space applications

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Publication Date 07/09/2014

Defense Date 05/12/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag beam,composite,FEA,FEM,finite element,nonlinear,OAI-PMH Harvest,shape memory polymer,shells,SMP,SMPC

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Yang, Bingen (Ben) (committee chair), Flashner, Henryk (committee member), Shiflett, Geoffrey R. (committee member), Wellford, L. Carter (committee member)

Creator Email dbergie@gmail.com,dbergman@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-434948

Unique identifier UC11286647

Identifier etd-BergmanDea-2643.pdf (filename),usctheses-c3-434948 (legacy record id)

Legacy Identifier etd-BergmanDea-2643.pdf

Dmrecord 434948

Document Type Dissertation

Format application/pdf (imt)

Rights Bergman, Dean

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA