Efficient Methods for Non-Linear Thermochemical Analysis of Composite Structures Undergoing Autoclave Processing by A L I R A S E K H B.Sc. (Civil Engineering), Sharif University of Technology, 1997 M.Sc. (Civil Engineering/Structural Engineering), Sharif University of Technology, 1999 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S FOR T H E D E G R E E OF D O C T O R OF P H I L O S O P H Y in T H E F A C U L T Y OF G R A D U A T E STUDIES (Civil Engineering) The University of British Columbia Apri l 2007 © A l i Rasekh, 2007 ABSTRACT Composite structures are increasingly being used in different industries. Their manufacturing imparts some challenges for the industry: most importantly prediction and control of the process to specification. The usual numerical solutions typically based on the use of the finite element analysis are not very suitable for large parts, especially when there is a need for quick estimation of the results for preliminary design and optimization. Therefore, there is a need both to have enhanced solutions that reduce the modeling effort for computer simulation of large and complex structures and also to simplify the solution and provide easy to use methodologies for quick estimations based on tables and charts. In the present work, a simple methodology is developed to estimate the temperature distribution in a thermoset polymer matrix composite slab placed on a tool and subjected to cycles of temperature ramp and hold leading to the curing of the composite and generation of heat due to the internal chemical reactions. Supplementary diagrams are also generated to set limits on the method. A modified finite element solution for heat transfer is also introduced that reduces the mesh generation and computational effort. This "higher order shell element" uses enhanced shape functions and efficient methods for spatial and temporal integrations in order to reduce the computational run times. The developed methods provide the design engineer with efficient analysis tools for predicting the temperature in a composite part. The simple diagrams and tables can be used for preliminary estimation of the temperature distribution in the part at each stage of the material development. The enhanced finite element methodology developed here can be used ii to reduce the amount of effort necessary in mesh generation and refinements necessary to achieve accurate solutions for thermochemical modelling of complex composite structures during processing. iii TABLE OF CONTENTS ABSTRACT II T A B L E OF CONTENTS IV LIST OF FIGURES VIII LIST OF T A B L E S XIII A C K N O W L E D G E M E N T S XIV DEDICATION X V 1 INTRODUCTION AND BACKGROUND 1 1.1 Autoclave Processing and Process Modeling 2 1.2 Research Objective and Thesis Outline 3 2 SIMPLE DIAGRAMS FOR T H E R M A L ANALYSIS 8 Nomenclature 8 Greek symbols 9 Subscripts 10 2.1 Introduction 12 2.2 One Layer of Material - Symmetric Boundary Conditions 14 2.2.1 Statement of the Problem 14 2.2.2 Applied Ramp Temperature 14 2.2.3 Applied Hold Temperature 21 2.2.4 Engineering Charts 21 2.3 Transition 22 2.4 One Layer of Material - Asymmetric Boundary Conditions 24 2.4.1 Statement of the Problem 24 2.4.2 Applied Ramp Temperature 25 2.4.3 Applied Hold Temperature 31 2.4.4 Approximate Solution 32 2.5 Thermal Analysis of a Composite-Tool Structure Using Simple Diagrams 35 2.5.1 General Formulation for Ramp Case .37 2.5.2 Applied Hold Temperature 44 2.5.3 Approximate Solution 45 2.6 Case studies 47 iv 2.6.1 Applied Ramp Temperature 48 2.6.2 Applied Hold Temperature 49 2.7 Temperature Hold and Heat Generation in the Composite 49 2.8 The Edge Effect 50 2.8.1 Case of Hold Boundary 51 2.8.2 Case of Heat Generation and/or Ramp Boundary 54 2.8.3 Asymmetric Boundaries and Two Layer Solids 56 2.9 Methods of Thermal Control 56 2.9.1 The Longitudinal Fin with a Rectangular Profile 57 2.10 Summary 59 2.10.1 Summary of the Approximate Solution Algorithms 59 2.10.2 Future Development 65 3 D E V E L O P M E N T OF A HIGHER ORDER S H E L L E L E M E N T FOR T H E R M O C H E M I C A L ANALYSIS 86 Nomenclature 86 Greek symbols 89 Subscripts 89 Superscripts 90 3.1 Introduction 91 3.2 Constitutive Equations 92 3.2.1 Heat Conduction and Diffusion in a Solid 92 3.2.2 Constitutive Equation for Exothermic Heat Generation 96 3.3 Parallelepiped Shell Element Formulation 97 3.3.1 Dirichlet Boundary Conditions 98 3.4 General 3D curved shell element formulation with stacking and tapering capabilities 103 3.4.1 Geometrical Shape Functions 103 3.4.2 Interpolation Functions 106 3.4.3 The Weak Form 108 3.5 Efficient Spatial Integration 111 3.6 Time integration 116 3.6.1 Time Integration of Heat Transfer Equation 117 3.6.2 Time Integration of Cure Kinetics Equation 119 3.7 Nonlinear Procedures 120 3.7.1 Newton-Raphson Iteration 121 3.7.2 Successive Substitution Method 124 3.8 Energy Norm 126 3.9 Errors and Adaptive Solution 127 3.10 Evaluation of Material Properties at Integration Points 129 3.11 Program Algorithm 131 3.11.1 Adaptivity Module 132 3.11.2 Boundary Contribution to the Stiffness and the Force Vector Module (Top and Bottom Surfaces) 133 3.11.3 Boundary Contribution to the Stiffness and the Force Vector Module (Edge Surfaces) 135 3.11.4 Cure Module 136 3.11.5 Database Module 138 3.11.6 Element Stiffness and Capacity Matrices Module 138 3.11.7 Element Temperature Distribution Module 140 3.11.8 Heat Vector Module 140 3.11.9 Integral Matrix Module 140 3.11.10 Jacobian Calculator Module 141 3.11.11 Place Finder Module 144 3.11.12 Smear Module 144 3.11.13 Superelement Module 144 3.11.14 Surface Shape Function Module 146 3.11.15 User Element Module 146 3.12 Summary and Future Work 147 4 CASE STUDIES 152 4.1 Introduction 152 4.2 Case I: Some Case Studies Using Simple Techniques 154 4.3 Case II: A Composite Slab Undergoing a Temperature Ramp with Asymmetric Boundary Conditions 155 4.4 Case III: A Two-Layer Material with Asymmetric Boundary Conditions Undergoing a Temperature Ramp 158 4.5 Case IV: Single Layer Behaviour during a Hold -Exotherm 161 4.6 Case V : Adaptive Solution using the M A T L A B Code 163 4.7 Case V I : Implementation of a Solution with Changing Number of Terms in ABAQUS164 4.8 Case VII: Subcycling Solution 165 4.9 Case VIII: Application of Higher Order Shell Elements in Modelling a Curved Composite Part on a Tool 166 4.10 Summary and Future Work 169 5 SUMMARY, CONCLUSIONS AND FUTURE W O R K 192 5.1 Summary 192 5.2 Conclusion 193 5.3 Future work 194 5.4 Contributions 195 REFERENCES 198 vi APPENDIX A. E X A C T T H R O U G H THICKNESS EIGEN FUNCTIONS FOR NON-DIRICHLET BOUNDARY CONDITIONS AND M U L T I L A Y E R PARTS 208 vii LIST OF FIGURES Figure 1-1 - Schematic of the autoclave temperature cycle (solid line) and the upper and lower limits set by the manufacturer on the thermal response of the composite (dashed lines) 6 Figure 1-2 - Component and tooling prepared in autoclave processing of thermoset composites (Johnston, [4]) 6 Figure 1-3 - Integrated sub-model approach for process modeling (Johnston et al., [1]). .7 Figure 2-1 - Typical autoclave process cycle 67 Figure 2-2 - Schematic representation of one dimensional transient heat conduction analysis with asymmetric boundaries 67 Figure 2-3 - Heisler chart for the temperature in the middle of a slab versus the Fourier number for applied constant temperature at the boundaries. Contours are for various Biot numbers 68 Figure 2-4 - The charts for the temperature in the middle of a slab versus the Fourier number for applied temperature ramp at the boundaries. Contours are for various Biot numbers 68 Figure 2-5 - Temperature distribution through the thickness of the slab in the steady state phase. Contours are for various Biot numbers 69 Figure 2-6 - The Fourier number at which a percentage of the steady condition state is acieved for the case of applied ramp temperature at the boundaries. Contours are for various Biot numbers 69 Figure 2-7 - Maximum error of strategy #1 in transient phase for the ramp heating of the autoclave; the error is presented as a function of BiB and hj I hB 70 Figure 2-8 - The global maximum errors extracted from curves in Figure 2-7 plotted as a function of BiB 70 Figure 2-9 - Maximum error of strategy #1 in transient phase for the case of temperature hold in the autoclave; the error is presented as a function of BiB and hj I hB. 71 Figure 2-10 - The global maximum of curves in Figure 2-9 as a function of BiB 71 Figure 2-11 - A as a function of hj I hB for 0.001 < BiB < 0.5 72 Figure 2-12 - A as a function of hj. lhB for 0.5 < BiB < 100 72 Figure 2-13 - Maximum values of A for all hj lhB as a function of BiB 73 Figure 2-14 - Maximum error of strategy #2 in transient phase for autoclave temperature ramp; the error is presented in terms of hy I hB for 0.001 < BiB < 0.5 73 viii Figure 2-15 - Maximum error of strategy #2 in transient phase for autoclave temperature ramp; the error is presented in terms of hj I hB for 0.5 < BiB < 100 74 Figure 2-16 - The global maximum of curves in Figure 2-15 and Figure 2-17 versus BiB. 74 Figure 2-17 - Maximum error of strategy #2 in transient phase for the case of temperature hold in the autoclave; the error is presented in terms of hj I hB for 0.001 < BiB < 0.5 75 Figure 2-18 - Maximum error of strategy #2 in transient phase for the case of temperature hold in the autoclave; the error is presented in terms of hj. I hB for 0 .5<5/ B < 100 75 Figure 2-19 - The global maximum of curves in Figure 2-17 and Figure 2-18 versus BiB. 76 Figure 2-20 - Typical two-layer composite studied in this paper ; 76 Figure 2-21 - Absolute values of the maximum errors incurred in estimating the temperature at the adiabatic line for the case of applied temperature ramp and for cases when the adiabatic line falls in the second layer (composite part) 77 Figure 2-22 - Absolute values of the maximum errors incurred in estimating the temperature at the adiabatic line for the case of applied temperature ramp and for cases when the adiabatic line falls in the first layer (tool) 77 Figure 2-23 - Maximum relative errors incurred in estimating the temperature at the adiabatic line for the case of applied temperature hold and for cases when the adiabatic line falls in the second layer (composite part) 78 Figure 2-24 - Maximum relative errors incurred in estimating the temperature at the adiabatic line for the case of applied temperature hold and for cases when the adiabatic line falls in the first layer (tool) 78 Figure 2-25 - Schematic of the model used to study the edge effect 79 Figure 2-26 - %ed as a function of Biot number for different error tolerances and for K/K2 =0.01 79 Figure 2-27 - ged as a function of Biot number for different error tolerances and for K/K2 =0.05 80 Figure 2-28 - %ed as a function of Biot number for different error tolerances and for K/K2 =0.1 80 Figure 2-29 - %ed as a function of Biot number for different error tolerances and for K/K2 =0.2 81 ix Figure 2-30 - %edge as a function of Biot number for different error tolerances and for K/K2=0.5 81 Figure 2-31 - £, r f as a function of Biot number for different error tolerances and for K/K2=\.0 82 Figure 2-32 - ged as a function of Biot number for different error tolerances and for K/K2=\0 82 Figure 2-33 - gedge as a function of Biot number for different error tolerances and for K/K2 =100 83 Figure 2-34 - <^edge_max as a function of K/K2 for different error tolerances. Markers show the maxima from the previous diagrams and trend lines are drawn through them 83 Figure 2-35 - ged as a function of K/K2 for different error tolerances for the case of Dirichlet boundaries (Bi-co) 84 Figure 2-36 - A typical longitudinal fin used to increase the efficient convective heat transfer coefficient 84 Figure 2-37 - Magnifying convective heat transfer coefficient by using fins with a rectangular cross section 85 Figure 3-1 - Schematic of the notions in heat transfer in a general solid as is presented in Section 3.2.1. The thinner line presents boundary condition of the first kind while the thicker line is the boundary condition of the second or third kind. 149 Figure 3-2 - A typical 4-noded heat transfer shell element. 149 Figure 3-3 - General 3D solid-like shell element in the physical coordinate system (top) and its mapped configuration in the parent coordinate system (bottom) 150 Figure 3-4 - Schematic of Gaussian integration points in surface (a) and through the thickness (b) of a shell element 150 Figure 3-5 - Schematic of typical layup (right side; full lines), integration points (<£, ,< 2^ , £ 3 ) and through thickness discretization points ( , < 2^, £3 > O defined in Section 3.5. This is drawn at integration point KGP with the thickness equal to DKGP. Also illustrated in this drawing is the tributary area of an integration point as is discussed in Section 3.10 151 Figure 4-1 - Magnified section of Figure 2-4 used to read 9m - Case II 178 Figure 4-2 - Magnified section of Figure 2-4 used to read 9m - Case III 178 Figure 4-3 - Rate of change of degree of cure as a function of degree of cure for different temperatures during isothermal curing of the composite - C A S E IV 179 x Figure 4-4 - Heat generated per unit volume as a function of time for different temperatures during isothermal curing of the composite - C A S E IV 179 Figure 4-5 - Error in exotherm prediction for single layer material using simple charts -C A S E IV 180 Figure 4-6 - Boundary conditions for testing the efficiency of the adaptive solution - C A S E V 180 Figure 4-7 - Comparison of the convergence of the series method with the F E M solution -C A S E V 181 Figure 4-8 - Time step sizes and number of terms used in the series solution - C A S E V.181 Figure 4-9 - Autoclave temperature, midplane temperature and degree of cure for Case VI . These are used as a guideline for choosing the patterns specified in Figure 4-10 : 182 Figure 4-10 - Different patterns of change for number of terms in case V I . These are mainly based on inspecting the solution for autoclave and midplane temperature and change in degree of cure in Figure 4-9 182 Figure 4-11 - Energy norm throughout the analysis for different patterns of number of terms as presented in Figure 4-10 - Case V I 183 Figure 4-12 - Mid-plane temperature in time for different patterns of number of terms as presented in Figure 4-10 - Case VI 183 Figure 4-13 - Final distribution of degree of cure through the thickness of the part for different patterns of number of terms as presented in Figure 4-10 - Case VI . 184 Figure 4-14 - Two element model of a composite layer modeled in Section 4.8 that is used to verify the subcycling solution - C A S E VII 184 Figure 4-15 - Autoclave temperature and mid-surface temperature for the case studied in Section 4.8 - Case VII 185 Figure 4-16 - A comparison of the errors in the midplane temperature for R U N 4 of Table 4-11 - C a s e VII 185 Figure 4-17 - A comparison of the errors in the energy norm for R U N 4 of Table 4-11 -Case VII 186 Figure 4-18 - A comparison of the errors in the midplane temperature for R U N 6 of Table 4-11 - C a s e VII 186 Figure 4-19 - A comparison of the errors in the energy norm for R U N 6 of Table 4-11 -Case VII 187 Figure 4-20 - 2D schematic of the model in Case VIII showing sections along which temperature distributions are sought 187 Figure 4-21 - A B A Q U S thermal elements are used to model both the tool and composite part-Case VIII 188 xi Figure 4-22 - Shell elements are used to model the composite and A B A Q U S thermal elements are used to model the tool - Case VIII 188 Figure 4-23 - Autoclave temperature, degree of cure, and temperature-time history both for cured composite (without exotherm) and curing composite(with exotherm) in the mid-plane of Section C (Figure 4-20) - Case VII 189 Figure 4-24 - Temperature distribution in the model at exotherm (t=173 min, temperature in degree C ) - C a s e VIII 189 Figure 4-25 - Temperature-time history on the top surface, mid surface and bottom surface of the composite from both the current shell element and A B A Q U S linear heat transfer elements. These pertain to Section C (as shown in Figure 4-20) - Case VIII 190 Figure 4-26 - Temperature-time history on the top surface, mid surface and bottom surface of the composite from both the current shell element and A B A Q U S linear heat transfer elements. These pertain to Section A (as shown in Figure 4-20) - Case VIII 190 Figure 4-27 - Through thickness temperature distributions for different sections of the composite ( Figure 4-20) from both the current shell element and A B A Q U S linear heat transfer elements at t=45 min - Case VIII 191 Figure 4-28 - Through thickness temperature distributions for different sections of the composite ( Figure 4-20) from both the current shell element and A B A Q U S linear heat transfer elements at t=173 min - Case VIII 191 Figure A - l - Typical temperature distribution in a three-layer composite 221 xii LIST OF TABLES Table 2-1 - Thermal properties for materials of interest 66 Table 4-1 - Temperature lag in aluminium tooling for different heat transfer coefficients, autoclave ramp rates, and tool thicknesses - C A S E 1 170 Table 4-2 - Temperature lag in invar tooling for different heat transfer coefficients, autoclave ramp rates, and tool thicknesses - C A S E 1 171 Table 4-3 - Temperature lag in composite tooling for different heat transfer coefficients, autoclave ramp rates, and tool thicknesses - C A S E 1 172 Table 4-4 - Temperature lag in steel tooling for different heat transfer coefficients, autoclave ramp rates, and tool thicknesses - C A S E 1 173 Table 4-5 - Composite cure kinetics equation and its variables 174 Table 4-6 - Convergence study for the current shell element - C A S E II 175 Table 4-7 - Convergence study for the linear element - C A S E II 175 Table 4-8 - Convergence study for the quadratic element - C A S E II 176 Table 4-9 - Error study for the minimum temperature in the plate from different methods of solution - C A S E II 176 Table 4-10 - Parameters used as input for C A S E III 177 Table 4-11 - Comparison of the efficiency of the subcycling method with the classical method of identical time-steping parameters - C A S E VII 177 xiii ACKNOWLEDGEMENTS I would like to acknowledge my sincere gratitude to several individuals who made this work possible. First, I would like to thank my supervisors Dr. Reza Vaziri and Dr. Anoush Poursartip. I really appreciate the support and technical guidance and help they gave me during the course of my PhD. I also enjoyed discussions with Dr Richardo Foschi and Dr. Michael Ward on the development of the models I used in this work. Many thanks to several past and present U B C Composites Group members and graduate students of U B C civi l engineering department for their friendship and help. Specially, I would like to thank Mr . Ahamed Arafath and Mr . Armin Bebamzadeh for the fruitful discussions and assistance they gave me on many occasions during the course of my work. I am also indebted to the special people in my life for their patience and persistent and unrelenting support they gave me during the course of my studies at U B C . Thanks to my dear daddy, mom, Saghi, Leila, Ehsan, Iman and Arman. Finally, I would like to acknowledge the Natural Sciences and Engineering Research Council of Canada and The University of British Columbia for their generous financial support. xiv DEDICATION I would/ IChe; to- cledAscate/ my fhe&Ly to-lovely mom/ cwid/ dx&d/ XV 1 INTRODUCTION AND BACKGROUND Fiber reinforced polymers are finding widespread applications as structural parts in many industries. Their popularity is due to many advantages that they have over the more traditional materials such as their light weight, high strength and stiffness and durability among others. Perhaps one of the main advantages of composite materials is the ability to manufacture very large and complex shape integral structures out of them. It is very important to predict the residual stresses and strains in the composite structural parts as they go through the manufacturing process. Such knowledge can be used in controlling and optimizing the manufacturing process to reduce the costs associated with assembly or scrap parts when they fail to meet the specified tolerances. Temperature profile through the thickness of the composite layer is an important factor in determining the acceptability of the process. If temperature at some point of the composite layer is outside the temperature range specified by the material manufacturer, then the process is not acceptable as the material properties specified by the manufacturer are recommended only for the material processed in the specified range of temperature. A schematic of the composite thermal response limits for a typical autoclave temperature cycle is presented in Figure 1-1. On the other hand, temperature distribution through the thickness of a composite part and its changes during processing in the autoclave is crucial in determining the residual stresses and deformations in the composite at the end of the cure process. The local temperature in the part is not only a function of autoclave heating (convection) but also a function of the internal heat generated in the resin material (a thermoset plastic) during 1 cure. The thermal simulation of both these phenomena is therefore a key part of the analysis used to estimate the residual stresses and deformations in a composite structure. Generally speaking, utilizing the computer simulations and numerical methods for analysis and optimization of large structures is not guaranteed to be a quick and easy procedure. In fact, it becomes quite expensive and time consuming i f it is not used intelligently. Computer simulations usually need huge amount of time to complete and it is often quite difficult (if not impossible) to scrutinize the accuracy of the results and determine the possible sources of errors. On the other hand, fast approximations of the solution are extremely useful in the initial stages of the design and therefore simple tools that provide approximate solutions are necessary both for the design stage as well as verification of more sophisticated analyses. In this chapter the autoclave processing of composite materials is briefly described and the objective of this research as well as the outline of the thesis is presented. 1.1 Autoclave Processing and Process Modeling Different methods can be used for manufacturing of thermoset fiber reinforced polymers. In the case of large structural components, autoclave processing is the preferred method. In order to manufacture these parts in an autoclave, thin layers of prepreg are cut and stacked on top of a tool to make a component. Prepreg is made of fibers impregnated in partially cured resin. The component could also contain honeycomb cores, pre-cured composites and adhesives. After assembly, the part is covered with bleeder and breather and then it is sealed inside a vacuum bag. This whole assembly is placed inside an 2 autoclave and undergoes a prescribed temperature and pressure cycle. This setup is schematically presented in Figure 1-2. When the above described composite component is cured inside an autoclave, a group of different physical and chemical phenomena occur within the material. These phenomena are triggered by the temperature rise and fall and cause final residual stresses in the component. Figure 1-3 presents these phenomena and the so-called integrated sub-model approach for their simulation [1]. 1.2 Research Objective and Thesis Outline The objective of this research is to develop and establish enhanced solutions for thermo-chemical analysis of fiber reinforced polymers in an autoclave. This is done in two fronts. One is developing simple diagrams that can be used for quick estimation of temperature distribution at different stages of cure. The second front is modifying the classical finite element solution based on general characteristics of the problem of interest so that it is less labour intensive to set up in terms of mesh generation and more efficient to run. The first goal is achieved through: • Generating simple 1-D thermal analysis diagrams for the case of a single layer of material with symmetric boundary conditions undergoing a temperature ramp. • Utilizing the above simple diagrams for thermal analysis of a single layer with asymmetric boundary conditions and estimating the error. • Introducing a method to use the same diagrams for two layers of material with asymmetric boundary conditions and estimating the error. 3 The second goal is to be accomplished by: • Introducing a higher order shell element with Fourier shape functions through the thickness. • Use of a subcycling solution for time integration. • Use of enhanced spatial integration for through thickness integrals. Meeting these objectives would provide a variety of numerical tools to the engineers both for simple and quick ways of analysis as well as efficient ways of carrying out more sophisticated and comprehensive analysis of heat transfer in composite structures during processing. The research presented in this dissertation is organized as follows: Chapter 2 - Simple Diagrams for Thermal Analysis: Presentation of the simple diagrams and analysis methods developed for thermal analysis of composites in an autoclave. Chapter 3 — Development of a higher order Shell Model for Thermochemical Analysis: This chapter describes the development of an enhanced thermochemical shell model for finite element analysis. Chapter 4 - Case Studies: In this chapter a few case studies are presented to apply and verify the solutions developed in Chapter 2 and Chapter 3. Chapter 5 - Summary, Conclusions and Future Work: A summary of the work presented in this thesis is presented in this chapter. The importance of the solutions presented in this work are then highlighted and recommendations for future improvements are given. 4 A P P E N D I X A - Exact Through Thickness Eigen Functions for Non-Dirichlet Boundary conditions and Multilayer Parts: Shape functions developed for other types of boundaries than the ones discussed in Chapter 2 and also for multilayer (more than 2 layers) parts are presented. The work described in this thesis is a continuation of several years of process model development in the U B C Composites Group. This work started with developing L A M C U R E (Smith, [2]) for heat transfer, resin cure, resin flow and fiber bed compaction. The work then continued with development of C O M P R O for thermochemical modeling, resin flow and residual stress analysis of composites (Hubert, [3]; Johnston, [4]). These efforts have continued through more current research in the group on viscoelastic constitutive model development (Zobeiry, [5]), this research on efficient solution techniques for thermochemical analysis, development of efficient solutions for stress analysis (Arafath, [6]), including contact surface and large strain deformation capabilities (Osooly, [7]). 5 Figure 1-1 - Schematic of the autoclave temperature cycle (solid line) and the upper and lower limits set by the manufacturer on the thermal response of the composite (dashed lines). Vacuum Bag Breather Bleeder Vacuum plug Sealant Figure 1-2 - Component and tooling prepared in autoclave processing of thermoset composites (Johnston, [4]). 6 Figure 1-3 - Integrated sub-model approach for process modeling (Johnston et ah, [I]). 7 2 SIMPLE DIAGRAMS FOR T H E R M A L ANALYSIS Nomenclature A, B, C, D, E, F parameters used in defining steady state distribution A cross sectional area of a fin Bi Biot number Bi2 Biot number in surface direction: Bi2 - hLjK2 Bie equivalent Biot number C circumferential length of a fin Et coefficients of the transient series solution Emax non-dimensional maximum value of the error of maximum lag (lead) EG maximum value of Emax for \<hT/hB<eo Fo Fourier number: aljL2 CP thermal capacity h convective heat transfer coefficient h equivalent convective heat transfer coefficient of a system of fins h equivalent convective heat transfer coefficient at the interface hss equivalent convective heat transfer coefficient at the interface at the steady state H exothermic heat generation per unit volume H adiabatic temperature rate K thermal conductivity in through thickness direction K2 thermal conductivity in surface direction L half thickness of slab Le effective half thickness of the slab nfm number of fins in unit length r temperature ramp in the fluid 8 R R number: (ti-r)l ja S thermal conductivity per unit length: K/L T slab temperature Terr if) m e difference between the exact and approximate temperature lag in the part Tlag (t) temperature lag in the part Terr Terr (t) at the moment when \Terr (t)/Tlag (t)\ maximizes ^hoia s t e P change in autoclave temperature T0 initial temperature of the system Tb temperature at the base of the fin Tb0 initial temperature at the base of the fin Tau,ocia,e autoclave temperature Tx fluid temperature Tss steady state temperature distribution in the slab TTr non-steady part of the transient solution u non-dimensional temperature distribution function through thickness u2 non-dimensional temperature distribution function in surface direction z through thickness coordinate Z, eigenfunctions Greek symbols a diffusivity in through thickness direction a2 diffusivity in surface direction PP parameter in estimating steady state Biot number p density of the material m non-dimensional length along a fin r\i non-dimensional eigenvalues for each layer S. non-dimensional parameter in eigenfunctions 9 dp non-dimensional location of the adiabatic line with respect to the midplane A width of a fin Xs eigenvalues cr the ratio of maximum norm of eigenfunctions in layer II to layer I E the Duhamel's auxiliary function E the Duhamel's auxiliary function - one dimensional £ non-dimensional coordinate through the thickness £ a location of adiabatic line in the first layer £ the coordinate in the second layer in 2 layers section % non-dimensional coordinate in surface direction in edge effect section £,a location of adiabatic line in the second layer < e^d the minimum £ from the edge to avoid edge effects £<fe*-max maximum value of 4edge for all Biot numbers q non-dimensional temperature parameter X non-dimensional distance from the adiabatic line 9 non-dimensional temperature lag (lead) 9m midplane, maximum or minimum 9 © reaction factor of the slab: 0.5 +1/2?/' I temperature at a point minus the ambient temperature u., £•,,/. parameters used in describing Ei (5P non-dimensional parameter A error involved in estimating the location of min (max) T by Sp A G maximum value of A for \<hj jhB < °o Y percentage of steady state achieved Subscripts B bottom surface of the slab T top surface of the slab 10 I, 1,2 layer number, 1: tool, 2:part (also 1,11) 1,2 heating ramp number e effective properties for the approximate methods P the value of the parameter at the adiabatic location n the value of the parameter at n percent of steady state 11 2.1 Introduction Engineers have long been interested in using simple temperature distribution diagrams for plates as solution to heat generation and distribution problems. Heisler [8] developed heat distribution diagrams for the case of one-dimensional unsteady heat conduction in a plate with symmetric convective boundaries and suddenly applied boundaries of constant value (temperature hold) in the fluid. Heisler diagrams are a result of a series of works by previous researchers, e.g. Gurney and Lurie [9] and Newman [10]. Heisler charts are used and referred to both in the literature and in the engineering textbooks like the classic book by Holman [11]. In development of Heisler charts, it is assumed that the material properties are constant. Suces and Hedge [12] have developed modification factors to compensate for the effect of change in material properties during the analysis so that Heisler charts are usable for a wider variety of problems. Another class of interesting problems are those involving constant internal heating, where there is a constant heat source per unit volume within the plate, which is independent of time and space. Fox [13] developed the governing equations for the solution and Young and Savino [14] prepared the corresponding charts. Other engineering charts that have been inspired from Heisler charts are the diagrams for the case of constant surface heat flux by Thibault [15] and the charts for transient temperature distribution for slabs subjected to thermal radiation by Zerkle and Sunderland [16] Bowley and Koenig [17] developed charts for some specific boundary conditions for composite cylinders and Fasiczka and Westermann [18] generated a number of diagrams 12 to be used in determining temperature distribution in slabs subjected to asymmetric convective heat transfer. Those diagrams are to be used alongside with convolution integrals to obtain the distribution for fluid temperature ramps. Compo [19] found these set of diagrams useful in predicting the transient conduction in materials with linear thermal properties. In this chapter, we first develop diagrams for the midplane temperature and distribution of temperature for the case of a temperature ramp applied to the symmetric convective boundaries of a slab with and without constant heat generation. Non-dimensional parameters wi l l be introduced in the simple equations developed for the steady state as well as the transient solutions. Secondly, we wi l l show that our diagrams and also the classical Heisler charts can be used in the case of asymmetric boundary conditions with a very small penalty in terms of absolute error. Next we investigate the application of the same diagrams for heat transfer and distribution in a two layer tool-composite hybrid. The errors would be more than the asymmetric case. It wi l l , however, still be useful for an approximate estimation. In the next section we look at the edge effect to see how much the two (and three) dimensional effects from the edges would affect the thermal distribution predicted by the simple diagrams which are based on one dimensional analysis. Finally, we present diagrams that are prepared to estimate the increase in the equivalent convective heat transfer coefficient on the tool side as a result of using fins. 13 These diagrams wi l l prove useful for engineers seeking immediate solutions to many practical heat transfer problems. They are ideal for what-if-exercises leading to preliminary design of materials and geometries for a given application. 2.2 One Layer of Material - Symmetric Boundary Conditions 2.2.1 Statement of the Problem Here we investigate the usefulness of engineering diagrams in one-dimensional heat transfer analysis of a slab of thickness 2L. The slab consists of a homogeneous material inserted in an autoclave and undergoes a heat cycle consisting of ramps and holds. The solution is obtained for a temperature ramp and hold with two identical convective heat transfer coefficients on both boundaries (symmetric case). 2.2.2 Applied Ramp Temperature A one dimensional heat transfer problem for a homogeneous one-layer solid with symmetric convective boundaries can be formulated as follows: dt U 2 M 2 - + H, t>0, \£\<l (2.1) 8T_ = ~{T-Tm) o n £ = +l,f>0 A. hi = +—(T-T„) onC = -\,t>0 (2.2) 8T T = oo T0+rt. (2.3) 14 Here H is the adiabatic temperature rate, defined as H = H, H is exothermic heat generated per unit volume, p is the density and Cp is the thermal capacity. T* = Tamodave i n t h e a b o v e equation is the autoclave temperature, r = tauloclme is the constant autoclave rate of heating (linear temperature ramp), and T0 is the initial temperature of the whole system. Also, a is the diffusivity, L is half-thickness and K is the thermal conductivity (note that a = ), and h is the convective heat transfer pCP coefficient, which is identical for both top and bottom surfaces. Finally, C, is the non-dimensional thickness coordinate with origin at the midplane of the slab. A n analytical solution to the problem above can be expressed in terms of an asymptotic solution and a transient one. The asymptotic solution is linear in time and is referred to as the steady state solution since it dominates the solution for large time. The transient solution is an inverse-exponential function of time and fades away after as time increases. 2.2.2-1 The Steady State Solution The steady state solution to Equations (2.1)-(2.3) has the following format: Tss=A{£)t + B{C) (2.4) Substituting this equation into Equations (2.1)-(2.3) leads to mathematical expressions to solve for A and B . The temporal and spatial derivatives of Tss can be written as: 15 57; d% K1 = A't + B' SS_ A" A"t + B" (2.5) 8Z ss dt where A' = dA/dC and A" = d2A/dC2 Substituting Equation (2.4) into Equation (2.1), we obtain: a (A"t + B") + H (2.6) As this equation should hold at all times (t), the following two identities must hold: A" = 0.0 B" = A-H 'a} Upon solving Equation (2.7) we obtain: A=CQ+D B = y6Cf+y2(D-H)£2 + E<; + F, (2.7) (2.8) where C, D, E , and F are integration constants. Now applying the boundary conditions (2.2) to the steady state solution, we get: A't + B'^-Bi(At + B-To0) on £ = +1 A't + B' = +Bi(At + B-Tx) on £ = -l (2.9) where Bi is the Biot number and is defined as: 16 K (2.10) Equations (2.9) should hold at all times and so we arrive at the following equations: A'= -Bi(A-r) on £ = +\ A' = +Bi(A-r) on £ = -1 (2.11) and B' = -Bi(B-T0) on ^ = +1 B' = +Bi(B-T0) on (2.12) Substituting the definition of A from (2.8) into Equation (2.11), one can see that C = 0 and Z) = r is the solution. Therefore A = r is the solution to Equation (2.11) and therefore we have: Substituting (2.13) reduces Equations (2.12) to the following: (2.13) +E -E ( 1 ^ + F fa \ 1 + = V Bi, U 2 / / 1 ^ + F fa \ 1+ — V Biy U 2 J {R@ + T0) (Re + T0), (2.14) where R,® are defined as (H-r)L2 a (2.15) 0 = 0.5 + — . Bi 17 The solution to Equations (2.14) is: E = 0 F = a (R® + T0) (2.16) Therefore, we have B = ~RC2+R® + T0 (2.17) So that the steady state temperature Tss becomes Tss=T0+rt + R\ U2+® T„+R 2 j (2.18) 2.2.2-2 Transient Solution In order to obtain the transient solution to Equations (2.1)-(2.3) we write the temperature as follows: T = Tss+TTr=T„+R-1 <Z2+® + T Tr (2.19) Upon substitution into Equations (2.1)-(2.3), we obtain dt f a \ d2T vL2j %-,\$\<Ut>0 (2.20) 8T Tr = -Bi-TTr on ^ = +1, f >0 = +Bi• TTr on £ = -\ , t>0 (2.21) 18 This problem is to be solved subject to some arbitrary initial condition. 7 ' 7 > = £ E / Z / ( ^ ) e x p ( - ^ 0 ) 1=1 (2.22) where, Z, = [ cos (40 i = 1,3,5,7,. I s i n ( ^ ) / = 2,4,6,8,. (2.23) and {4} is the set of eigenvalues obtained from the following equations: tan ( 4 ) = Bi r +— for i = 1,3,5,7,. A. - A for ,- = 2,4,6,8, I 5/ (2.24) Fo in Equation (2.24) is the Fourier number which is defined Fo = at (2.25) In order to get coefficients {E.} we should apply the initial condition for TTr, that TTr n=Ta-T„ =-R- —C+0 1 ~2 (2.26) or -R- ~C 2+e] = SE,z,(0 J 1=1 (2.27) Now using the orthogonality of {Z,} , we have E, = - i ? — V, +©£•, 2 ' (2.28) 19 where s. = ^ 1 J(z,)2*T - i - 1 J(z,)2<*r (2.29) Substituting all these into Equation (2.19) leads to: rp rjn -I 00 U, + ©£, 2 ' Z,(C)exp(-A, 2 Fo) (2.30) For the case of mid-plane temperature where C, = 0, the equation above is simplified to 0- = T~jn^rL=&-t(-U+&£) Z< (°) exp(-A 2 ^) (2.31) The right-hand-side of Equation (2.31) has only one parameter involved that is the Biot number and therefore this equation can easily be presented in an engineering chart. On the other hand, for large values of time when Fo —> oo, we would have ^ = © = 0.5 + Bi (2.32) 20 2.2.3 Applied Hold Temperature In this case the problem is formulated as in Equations (2.1) and (2.2) except that the autoclave temperature would be T„=TX (2.33) The solution is well established in heat transfer literature and can easily be calculated ^ = 2 > ' Z ' (C)exp(-A,. 2Fo) (2.34) - ' o ^1 '=1 The parameters and functions involved in the above equation are identical to those defined in Equations (2.23) to (2.28). 2.2.4 Engineering Charts The discussions in the previous section lead to preparation of some very useful engineering charts. The usefulness of Heisler charts has already been established in the field of heat transfer. Heisler chart (Figure 2-3) basically provides the midplane temperature in the case of the temperature hold problem; that is when the boundary conditions of the one-dimensional problem are constant specified temperature (hold). Our proposed chart (Figure 2-4) is based on Equation (2.31) and is a powerful tool in prediction of temperature lag in symmetric autoclave heating. It wi l l be shown in the next sections that this chart can also be used in asymmetric cases and two-layer structures with little error. 21 The heat distribution in the steady state is shown in Figure 2-5. In the case of symmetric boundaries, 0 < C, < 1. The rest of the diagram (i.e. £ > 1) wi l l be useful for asymmetric cases. 2.3 Transition Another important question that should be answered through the use of simplified diagrams is the "transition" time for a slab to react to a change in the autoclave temperature ramp and go from the steady state phase of the first ramp to the next. Specifically, the transition time for a change from hold to ramp or vice versa, is of interest to the autoclave process engineers. In order to calculate the transition time when the R number changes from i?, to R2, we track the solution that leads to Equation (2.26) with the following modifications to the initial condition and the steady state solution: r 0 = r . ^ + J U - k 2 + ® (2.35) ( 1 „, ^ And hence: 7 > l , = o J (2.36) In the midplane surface we have (^ = 0): (2.37) Therefore, the solution for E, is modified to: 22 E , = - ( * 2 - * . ) -O, +&£i (2.38) Now, it can be assumed that in the time close to the steady state phase of the second ramp, only the very first term in the transient solution is active and so the solution is simplified to: Tm-T^R2®-(R2-Rx) Tm ~TmSS =-{R2-R\) 1 V, + &£, 1 ' f 1 ^ V, +&£. v 2 1 Z,(0 .0)exp(-^ 2 Fo) (2.39) Z, (0.0)exp(-A? Fo) . Comparing with Equation (2.37) we obtain (Z , (0.0) = 1.0): 1 .0 -Y/100 : T -T 1m 1 mSS T -T mO mSS t=0 — u, + 0 f , exp( -^ Fo) _2 )_ 0 (2.40) Here, Y is the percentage of the steady state that is achieved. From Equation (2.40) we can obtain an approximation of the time (Fourier Number) at which the steady state has been attained: Foj --1.0 In f \ ( 1 . 0 - Y / 1 0 0 ) © -1.0 A2 l n ( l . 0 - Y / 1 0 0 ) + ln \ 0 1 V, +©£•. v 2 1 1 J ) (2.41) We can now consider the extreme cases. For large Biot numbers, from Equations (2.24) and (2.32) we have \= — , 0 = 0.5, and In Therefore f ^ 0 ^ " 2 u, + 0 f , 3.1546x10" 23 FoT =-0.4052847 x ln (1 .0 -Y /100 ) + 0.01278521. (2.42) On the other hand, for the cases with small Biot number we have A s yfWi, 0 s — and 4 Bi In 0 V 2 1 1 = 0.0. This gives that For = In ( l .0 - Y/100) . (2.43) Equation (2.41) is plotted in Figure 2-6 for the range of Biot numbers that are of interest in this study. 2.4 One Layer of Material - Asymmetric Boundary Conditions 2.4.1 Statement of the Problem This section investigates the usefulness of engineering diagrams in one-dimensional heat transfer analysis of a slab of thickness 2L. The slab consists of a homogeneous material under an autoclave heat cycle that consists of ramps and holds. The convective heat transfer coefficients might be different on the top and bottom surfaces and it is assumed throughout that hB < hr unless otherwise mentioned. The solution is obtained for different cases of boundary conditions; autoclave ramp and hold under different convective heat transfer coefficients on both sides (asymmetric case). 24 2.4.2 Applied Ramp Temperature 2.4.2-1 Steady State Solution A one dimensional heat transfer problem for a homogeneous one-layer solid with asymmetric convective boundaries on two sides can be formulated as follows: dT ( a\ d2T dt V1- J dt2 + H,\C\<\,t>0 (2.44) ^ = _ M ( r _ r ) 0 n 4 - = + i , / > o d£ K v ^L = +hsL(T-Ta>) on C = -ht>0 (2.45) (2.46) Here hB,hj. are the convective heat transfer coefficients for the top and bottom surfaces. This system reaches a steady state condition in which the temperature at each point of the solid changes linearly with time, that is: Tss=A{t)t + B{t). (2.47) In order to calculate the spatial functions A and B, Tss is substituted for T in Equation (2.44) to yield: A = (A"t + B") + H (2.48) thus A" = 0 A = C£ + D 25 (2.49) and, (A-H) B = /6C<;3+y2(D-H)C2 + E<; + F (2.50) where C, D , E , F are integration constants. Similarly, Tss is substituted for T in Equations (2.45) to obtain A't + B' = -BiT{At + B-Tm), on £ = +1 A't + B' = +BiB(At + B-T„), on ^ = -1 (2.51) which yields A' = -BiT(A-r), on£" = +l A' = +BiB(A-r), o n £ = - l . (2.52) Here BiB,BiT are Biot numbers corresponding to the bottom and top boundaries and are defined as B i B = ^ K hjL BiT = K (2.53) Now using Equation (2.49) leads to A = r (2.54) Another set of equations derived from Equation (2.51) is: B' = -BiT(B-T0), on £ = +1 B' = +BiB(B-T0), on £ = -1 (2.55) 26 Upon using Equations (2.50) and (2.54), we obtain +E f 1 ^ ( -E 1 + BiT 1 1 + V "*T J \ Bi + F = + F = I V 0.5 + B J R 0.5 + Bi 1_ BL + Tn T J \ (2.56) + T0 J J where R = (H -r) — , as defined previously in Equation (2.15). Solving the system of a Equations (2.56), we obtain E = F = R8„ T0 + R 0.5 + — 0, p J) (2.57) where /? , 8 are defined as: 1 1 2 ^ + — + — -BiT BiB BiT .BiB . i i 2 + — + — BiT Bia (2.58) 8P= B l \ B \ • (2-59) 2 + — + — Bij BiB Here 8 is actually the location of the adiabatic line where there is no heat flux. So the steady state solution is 27 f 1 ^ 0.5 + — V PPJ + R5PZ-V)R?.. (2.60) which can be written as T -T R 0.5 + -PP (2.61) The minimum value for the right hand side of Equation (2.61) occurs at C, = 5p. Shifting the origin of the coordinate system to this point leads to (£" = 5 + %): T -T 1SS 1<a _ R { 1 ^ 0.5 + — (2.62) which simplifies to T -T *ss - l ° o R 0.5 + ^ + ^ -PP 2 - V y2 (2.63) Now Bi. is defined as Bi 1 - ± + g - l PP 2 1 1 + + Bi ( BiB Bij .BiB ~ i i 2 + + -1 + — 2 Bi-,. Bir, 1 1 Bij Big n • r 2 + + — BL Bi V B J (2.64) which simplifies to 28 Therefore, Equation (2.63) can be rewritten as: R 00 (2.66) X in the equation above varies in the range \ -\-5p,+\-Sp J . As Equation (2.66) i is symmetric in %, and 5 changes in the range , then it suffices to have contours of Equation (2.66) for the range [0,+2]. Therefore Figure 2-5 shows the temperature distribution. The distribution in the range [0,+l + S ] would correspond to the length dT underneath the adiabatic line (i.e. the line corresponding to-— = 0) and the range Equations (2.61) and (2.66) can equally be used to describe the steady state distribution of temperature in the slab. 2.4.2-2 Transient Solution The transient solution can be determined from Equations (2.44)-(2.46). Using Equation (2.61), the temperature equation can be rewritten as follows: [ 0 , l - £ ] to the range above this line. (2.67) Substituting the above summation into Equations (2.44)-(2.46) yields 29 87^ dt 8% (2.68) 571 — = —BiT • TTr on t = +l,t>0 dT7 Z- = +BiB-TTr on £ = -\,t>0 (2.69) with initial condition T (<£", 0) = T0. This problem also has a classic solution based on Fourier series, given by ^ = l E , Z , ( C ) e x p ( - ^ F o ) , /=i (2.70) where Z , = cos (A, ( £ - * , ) ) . (2.71) Here {/I,} is the set of eigenvalues attained from the following equation: det -Ai sin {Aj) + Bir cos (A,) +Ai cos (A.) + BiT s'm(Aj) +Aj sin (A l) - BiB cos (A,) +A, cos ( A,-) + 5 / B sin ( l , ) = 0 (2.72) where A e v y 2 2 and each Ai is associated with a from the following equation: t a n ( ^ , ) = --At sin (A t) + Z?/;. cos (A,) +A, cos (/I,) + BIT sin ( l ( ) (2.73) In order to determine the coefficients {E,} we need to apply the initial condition for T7 Tr • This yields, 30 Tfr ,=0 - Tss (=0 - R • 0.5+—+SX-lAc2 (2.74) Upon using the orthogonality of {Z,.}, we obtain the coefficients E, as E, = -R r 1 ' — vt + V 2 0.5 + -PP (2,75) where jz , .z ,^ (2.76) Substituting all these into Equation (2.67) leads to: ryi rji 1 CO f i r „ . i ^ A 7? V 2 ' 0.5 + — v PPJ £ , + < y , Z , ( C ) e x p ( - ^Fo) (2.77) The right-hand-side of Equation (2.77) depends on pairs of Biot numbers (BiB,Bir). Therefore it would be difficult to present it as an engineering chart. In the next section a few simplified solutions are discussed that can reduce the number of parameters. 2.4.3 Applied Hold Temperature In this case the problem is formulated similar to Equations (2.44) and (2.45) except that the autoclave temperature is constant: T = T (2.78) 31 The solution in this case is: T-T. = £ * , z , ( 0 e x p ( - 4 2 F O ) . (2.79) The parameters and functions involved in this equation are taken from Equations (2.71)-2.4.4 Approximate Solution 2.4.4-1 Strategy #1 The solution obtained for the steady state can be used to give insight into an approximate solution for the transient problem. Comparing Equations (2.18) and (2.66) shows that using the equivalent Biot number given by Equation (2.65), the accurate steady state temperature can be calculated from the symmetric equations. This suggests that we use the same equivalent Biot number for the transient case as well which would mean using Figure 2-3 and Figure 2-4 with the equivalent Biot number given by Equation (2.65). The rest of the parameters involved would be calculated the same way as in the symmetric case. This method of solution is inexact in the transient phase. A n error study performed using a M A T L A B code was performed. Figure 2-7 shows the maximum error involved in each transient analysis. The y-axis shows the relative error defined as: (2.77). E, max (2.80) 32 where Tm is the minimum (maximum) temperature within the thickness (i.e. at the stationary location of the adiabatic line) and Tm is the approximate value of the same state variable using the equivalent Biot number in Figure 2-3. The error is shown for different BiB for a wide range of hj I hB. Figure 2-8 shows the maximums for each BiB or: E G = m a x ( E m a x ) (2.81) Figure 2-9 and Figure 2-10 show similar error for the case of hold temperature. E m a x is defined slightly different in this case, that is max (|r - f j ) E - Q ^ ° ° V I !Z (2.82) max It* T* The definition of E G is the same as that in Equation (2.81). 2.4.4-2 Strategy #2 As another strategy for an approximate solution, we may assume that the location of the adiabatic line (S ) as is specified by Equation (2.59) is constant throughout the transient dT analysis as well. In fact, this is not true and the adiabatic line, where — = 0, moves in dC, the transient phase. The assumption of non-moving adiabatic line leads to a good approximation of the solution in transient phase. This solution strategy would be accurate through the analysis 33 (both steady state and transient) not only in the case of symmetric boundaries but also in the extreme case of asymmetric boundaries, i.e. when BiB = 0.0. The mathematical validity of this strategy is related to the numerical observation that the minimum point of the steady state solution (5 ) which is the adiabatic line in that phase of analysis is located close to the minimum point (Sx) of the first trigonometric term; ^ = 0 0 5 ( 4 ( ^ - 4 ) ) (2-83) In order to establish the relation between these two parameters, we define: A = \Sp-S\ (2.84) The variation of A with the bottom Biot number and the ratio hr I hB is investigated in Figure 2-11 and Figure 2-12. Figure 2-13 shows the global maximum value of the parameter A A r = max (A) (2.85) which is a function of the bottom Biot number. Besides, similar trend is observed for the errors in temperature prediction as follows. In strategy #2, Sp is first calculated from Equation (2.59). The effective half-thickness, L , is then estimated as: Le=L(\-Sp) (2.86) 34 This effective half-thickness is used in the calculation of Fourier number, the R number and Biot number. With the advantage of these numbers, Figure 2-3 and Figure 2-4 can be used to obtain good estimate of the temperature lag. Figure 2-14 and Figure 2-15 show the error involved for different BiB and hj/hg ratios in the case of temperature ramp. Figure 2-16 shows the maximum error for any hjjhB ratio. Figure 2-17 and Figure 2-18 illustrate the error in the case of hold temperature and Figure 2-19 is the diagram for the global error in this case. The nomenclature for all these diagrams is the same as those in the previous section. A comparison of Figure 2-7 with Figure 2-14 and Figure 2-15 establishes that for the cases where h, and hB are of the same order, strategy #1 yields lower errors. Generally speaking, however, one would compare Figure 2-8 and Figure 2-16 to see that strategy #2 is more reliable as the error hits a global maximum of about 4% compared to 7% in strategy #1. Strategy #2 is also more physically sound and is recommended for engineering use. 2.5 Thermal Analysis of a Composite-Tool Structure Using Simple Diagrams The transient and steady state temperature distribution in a composite part resting on a tool is of practical importance in the process modeling of composite structures undergoing autoclave manufacturing procedure. Knowing the design temperature ramps and holds, the autoclave convective properties, the thermochemical properties of the resin and the thermal properties of the composite and the tool one can create an F E M model to obtain the temperature distribution in space and time. This approach is too complicated 35 and time-consuming and gives no insight into the solution of the problem. It is also a very expensive approach for design purposes, as each case study wil l need significant amount of computational effort. In this section we extend the use of simplified charts generated previously in this chapter for the case of one layer material to the cases of interest in the process modeling of a composite structure on a tool material. This wi l l prove to provide simplicity, insight and cost effective solution. Computation of thermal lags and leads as well as temperature distribution in a thermoset composite part during processing in an autoclave is of considerable interest to engineers in rapidly developing composites industry. The reason is that temperature distribution in the part during manufacturing directly affects the residual stress distribution in the part and hence it is an important factor in dimensional control of the part. On the other hand, the research on temperature distribution in a multi-layered slab is not restricted to composites industry. Many researchers have studied efficient analytical and numerical strategies to estimate the temperature distribution in such structures. Title [20] developed a one-dimensional orthogonal expansion for a two-layer composite and Padovan [21] developed a generalized Sturm-Liouville procedure for composite and anisotropic domains in transient heat conduction. These approaches as well as Green's function approach and the use of Laplace transform are discussed in the book by Ozisik [22]. Haji-Sheikh et al. [23] used the exact solution of heat conduction in composite materials for an inverse problem, where they seek the heating rate that causes a temperature distribution. More recently, de Monte [24] introduced a method that he calls "natural analytical approach" to make the eigenvalue solution more efficient. He then used this strategy in analytical solution of unsteady heat conduction of composites [25] and also in problems 36 of higher dimensions [26],[27]. Another paper by Haji-Sheikh and Beck [28] tries an efficient method of computing eigenvalues for heat conduction which was subsequently used in multi-dimensional multi-layer bodies [29], [30]. Sun [31] generalized the analytical solution of de Monte to three layers. L u et al. ([32]) have used novel analytical approach for conduction in composites slabs which was then generalized to multiple-dimensions ([33]). Here we examine the use of the same simple diagrams developed in previous sections for one layer material under symmetric and asymmetric boundary conditions for the case of two-layer material in one dimension. First, an analytical solution for both the steady state (asymptotic) and transient phase of the problem wi l l be presented. The transient solution is based on eigenvalue analysis and orthogonality of eigenfunctions. A novel technique is then introduced that computes the equivalent boundary conditions for an equivalent single-layer problem. This computation is approximate only and leads to a quick use of the transient charts developed previously. This method is then fully tested for the cases of interest in the manufacturing process of composite materials. Error diagrams are then prepared to show the limits of the validity of the solution. 2.5.1 General Formulation for Ramp Case The P D E for one-dimensional heat transfer analysis for two-layers of solid with general convective boundaries as shown in Figure 2-20 is as follows: 37 5T dt d2T • a dz d2T 2dz2 0 < z < 2L + H2 2Lx<z<2(Lx+L2) (2.87) K = +hB.lT_T\ o n 2 = 0 dz KXK ' ^—^(T-TJ on z = 2(Lx + L2) oz K2 (2.88) T„=T0+rt (2.89) The solution must also satisfy the following continuity equation at the interface of the two layers: dT dz 5z z = 2 i , + (2.90) In order to solve the above system of equations we define separate local coordinate for each layer, namely £ and £ , which are defined as follows: z-(2Lx+L2) (2.91) and so the P D E and boundary equations transform: 8TX _ M d2Tx dt ~ dC2 dT2 _ ( aj d% dt ~ ^ 2 ; d? • + Hn -\<C<+\ - 1 < £ < + 1 (2.92) 38 8TX 84 +BiB{Ti-Ta)onC—\ (2.93) --Bir(T2-Tx) on 4 = where we have Bi. = k ^ 1 (2.94) BiT K2 Here Tx, T2 are the values of the temperature field T in the regions 1 and 2. Therefore the continuity equations can be rewritten as: T\ = T \ Lx 8C('X ~ L2 84 'f=_1 2.5.1-1 The steady state (asymptotic) solution We now assume that the steady state solution has the following structure: Tlss=T„+Bl(t) = rt + T0+Bl(<4) T2SS=TK+B2{4) = rt + T0 + B2{4). Substituting the above into Equations (2.92) results in 39 B 1 2 - ^ +*2$ + fl (2.97) where tf2 = v a2 j (H2-r). (2.98) Using the boundary conditions (Equations (2.93)), we can write + ex - +BiB ( 1 ^ V 2 ( 1 ^ R2+e2=-Bir -~R2+e2+f2 v 2 ; (2.99) While the continuity Equations (2.95) lead to the following: ~Rl+el+fl=-^R2-e2+f2 Sl(-Rl+ei) = S2(R2+e2) where (2.100) ^2 (2.101) Equations (2.99) and (2.101) can be rewritten as 40 f 1 ^ 1+-1 + Bi, 1 Bj \ Bij e2+f2=R2 1 f \ 1 ^ — + — ^ 2 5/y J SB ST (2.102) el+fi+e2-f2=-(Rl-R2) (Sl)el-(S2)e2=SlRi+S2R2 The above equations can be solved for the cases of interest using a simple computer code, M A T H C A D , or a hand calculator. Obtaining the solution of the above simultaneous system of equations, we can use them in calculating the parameter of interest to us which is the heat flux at the interface between two materials: i K l ) 8T{ 8T2 K L\ ) <r=+i \L2 ) 3# (2.103) This parameter divided by the temperature difference between the interface and the autoclave gives us a measure of the equivalent convective heat transfer coefficient at the interface, that is: h = -T -T. oo i (2.104) Here qi can be positive or negative and so is h . If h > 0 then it can be used as an equivalent heat transfer coefficient for the top of the bottom layer, otherwise, it is the equivalent heat transfer coefficient at the bottom of the top layer (layer 2). In the case of the steady state (asymptotic) solution, one can calculate the equivalent convective heat 41 transfer coefficient from the solution to the Equations (2.102) and use of Equation (2.104), which yields BiBBiTSxS2(/?, (1 +1 / BiB)-R2(\ + \lBiT)) SXRX (1 + BiB )(1 + 2BiT ) + S2R2 ( l + BiT ) ( l + 2BiB ) ' (2.105) This parameter results in an accurate solution for the steady state case. It can also be used to obtain approximate h value for the transient case. 2.5.1-2 The transient solution The general solution to Equations (2.87)-(2.89) can be written in the following format: T = T +T 1i 1iSS T -M7V -*2 2SS 27> ' (2.106) Substituting into Equations (2.92) and (2.93) we obtain dt ( aAd2T ITr dT 2Tr dt \L2 J 2Tr - 1 < £ < + 1 - 1 < £ < + 1 . (2.107) dt • = +BiB.Tm@t = -\ dT 2 _ dt = -BiT.T2rr@4 = +l. (2.108) The continuity relations are that T = T Kl dTlTr I _ K2 dT2Tr Lx dt |f=1~ L2 dt K - 1 (2.109) 42 The solution has the following form r ™ = E % ( £ ) e x p ( - 4 2 0 ;=1 TTr2=fjE<Zi2({)exp(-Al2t) (2.110) Substituting into Equation (2.107) results in: z" + V «i J \ a2 J Z „ = 0 -1 Z 2 , = 0 - ! < £ < + ! This system has the following solution: z w =cos( j7„(£-<5i ( ) ) Z 2 l . = < r , . c o s f a 2 , ( # - £ 2 , ) ) where %i=A,2L22/a2 (2.111) (2.112) (2.113) Substituting Equation (2.112) into the boundary and continuity Equations (2.108), (2.109) we obtain det Vi,-sin fa,)- 5 / f l c o s f a „ ) 0 cos fa,) .sin fa,)' ? i ,cosfa , ) + 5 / f l s i n f a , ) 0 sinfa,,) cosfa,) 0 -7 2 , s infa 2 , ) + i? / r cosfa 2 , ) -cosfa 2 l . ) -7/2, sin fa,) 0 ^ , c o s f a , ) + J g / r s i n f a , ) s infa 2 / ) -77,, cosfa , ) (2.114) 0 43 The roots of (2.114) determine the set {A,}, while {£,,}, {S2j} are calculated from the following: >7i,sin(j7„)- BIB cos (77,,) tan(/7,A) = -tan(/72,^2,) = 7i,cos(/7„) + 5 / B s i n ( % ) -Vzi sin (tin ) + B I T c o s (VI,) (2.115) /72/cos(?72,) + 5/ 7 ,sin(77 2,) On the other hand, we have: ^ _ c o s f a Q . O - ^ ) ) cos (^ 2 , ( -1 .0 - J 2 , ) ) " (2.116) The coefficients {£.) are calculated from the following as a result of the orthogonality of the trigonometric functions: +i +i (pCP\Ly\{T0-Txss)cos(rju {£-Su))d£ + a, (pCP\ L2. \(T0 -T2SS)cos(rj2l (<f-S 2 i ))d£ E.= =1 - =L (pCP \ Lx.jcos2 (Vu (<T - Su )) dt + a 2 (PCP )2 L2.jcos2 (rj2l (<f - S2>)) d^ - i - i (2.117) 2.5.2 Applied Hold Temperature This is the case where the fluid temperature jumps to a constant temperature TA*T0. In this case the steady state temperature is obviously Txss = T2SS = TA throughout the composite and the transient temperature is calculated from Equations (2.106) and (2.110) with the parameters from Equations (2.114)-(2.117). 44 2.5.3 Approximate Solution In order to obtain a quick estimate of the extremum temperature in a bi-layer solid and also the temperature distribution, one can use the equivalent heat transfer coefficient from Equation (2.105) throughout the analysis. This wil l lead to a one-layer asymmetric solution, which can be further simplified by locating the steady state location of the adiabatic line as described previously for the asymmetric cases. The following is the solution algorithm for calculating the temperature lag in the part (layer 2 or top layer): 1) Evaluate BiB,BiT,RuR2,Sl,S2 according to Equations (2.94), (2.98), (2.101). 2) Calculate the interface heat transfer coefficient h from Equation (2.105). 3) If h > 0 then h is the equivalent heat transfer coefficient for the first layer, that where Bin is the equivalent Biot number for the upper surface of the first layer. The location of the adiabatic line is is: hri = h (2.118) Bi. ri Bi, B (2.119) 2 + + BI T l And so we have 45 max ( a / , , I ) Re=Rr(\ + \Q2 (2.120) Foe^/(l + \Ca\f. These equations are to be used to read the extremum temperature and the temperature distribution in the first layer according to Figure 2-3, Figure 2-4 and Figure 2-5. Note that in the case of temperature hold a fictitious heating rate (like r -1) could be used to get the location of the adiabatic line. 4) If h < 0 then -h is the equivalent heat transfer coefficient for the second layer, or: hB2 = -h BiB2 = -hL2 (2-121) K2 where BiB2 is the equivalent Biot number for the bottom surface of the second layer. The location of the adiabatic line in this layer become: J 1_ B t \ B I b \ • (2-122) 2 + + -Bi, BiB2 And so we have: Bie=max(BiB2,BiT)x(] + \£a\) Re=R2x(l + \{a\f (2.123) Foe=^/(l + \{a\y. L2 46 which are to be used to read the extremum temperature and the temperature distribution in the second layer according to Figure 2-3, Figure 2-4 and Figure 2-5. Again, in the case of temperature hold a fictitious heating rate should be used to get the location of the adiabatic line. 2.6 Case studies Extensive case studies were performed to verify the simplified solution for the cases of interest in process modelling of thermoset composites. In these cases, the first layer is taken to be the tool and the second layer is the composite part. The autoclave temperature consists of successive ramps and holds and the composite material generates exothermic heat at some stage during the cycle. The error in using the simplified method is guaranteed to be zero for the steady state phase. Therefore, all the errors occur in the transient phase of the analysis. A tool-part combination was studied undergoing different autoclave temperature ramps and also temperature hold. In addition, different tool materials and a reasonable range of tool and part thicknesses were studied. Besides, different values for autoclave convective heat transfer coefficients on the top and bottom surfaces of the tool-part set were considered in these studies. Table 2-1 lists the material properties for different tools of interest as well as the thermoset composite. Tool and part thicknesses in this study are in the range of 1mm-10cm and the autoclave convective heat transfer coefficient varies between 5 WI m2K and 200 W lm2K. 47 2.6.1 Applied Ramp Temperature For these cases, ramp rates of r = 1,2,3 C / min were studied. The autoclave temperature was allowed to vary between 20C and 25 OC and the case with temperature lags more than 100C were ignored. The ranges of thicknesses, material properties and convective heat transfer are as mentioned above. Figure 2-21 shows the absolute values of the errors (Terr) involved versus the autoclave temperature (Tauloc!aw) and the temperature lag (Tlag) for the cases where the (steady state) adiabatic line occurs in the second layer (h < 0). The temperature lag (Tlag ) considered here is the difference between the autoclave temperature and the temperature at the location of the adiabatic line. Figure 2-22 is a similar diagram for the cases where h > 0 and the (steady state) adiabatic line is in the tool. The error at any time during the transient analysis (Terr (t)) is the difference between the minimum temperature in the part calculated using the exact method (Equations (2.106)-(2.117)) and the approximate solution (Equations (2.118)-(2.123)). The error in every approximate transient analysis (T e , , ) in these diagrams is recorded at the time when the ratio \Terr {t)JTlag (?)| reaches a maximum, where Tlag (t) is the accurate temperature lag. At the same moment, Tlag and Tautocalve are also recorded for the ordinates. Both diagrams show errors equal or less than 7.5% for the ranges of parameters considered in the study. For the cases where temperature lag is less than 20C, both diagrams show errors of about 0.5C or less which proves them useful for practical cases. 48 2.6.2 Applied Hold Temperature Similar error diagrams were prepared for the case when an abrupt constant temperature ( A r W d ) is applied in the autoclave. Figure i2-23 and Figure 2-24 show the error for the cases where h < 0 and h > 0 respectively. Both diagrams show the lag temperature versus the error in temperature lag and are normalized with respect to AThold . These parameters are recorded in the diagrams for different cases of interest at the moment when \Terr (t)/Tlag (?)| reaches a maximum. Both diagrams show relative errors of less than 10% of AThold for the transient solution. 2.7 Temperature Hold and Heat Generation in the Composite Another case of interest in the process modeling of composite materials is when there is a constant heat generation in the composite and the temperature in the autoclave is constant in which case Equation (2.105) simplifies to ~ s s = BigS, (1 + 25/ , ) -h~L This equation shows that the composite could be modeled with a BiB2 in the K2 steady state phase. This approximation didn't prove to be reliable for the transient phase. Therefore, there is no general way using the simple charts for the case of heat generation in the composite part. 49 2.8 The Edge Effect The one-dimensional solution discussed above has some limitations for practical purposes. One important limit is the edge effect where the heat transfer through the edges of the slab affects the closed form solution presented in the diagrams of the previous In order to study the edge effect, a two dimensional case of autoclave hold temperature is considered here with constant convective heat transfer on the edge. The in-plane dimensions are assumed to be much bigger than the thickness and hence the solution presented here is for the case of infinite in-plane size as is illustrated in Figure 2-25. The part is homogeneous and single layer although conductivity and diffusivity could be different in the in-plane direction compared with the through thickness direction. In order to retain the notation that is used in this chapter so far, the through thickness conductivity and diffusivity remain K and a while K2 and a2 wi l l denote the in-plane direction conductivity and diffusivity. The other non-dimensional parameters defined These parameters are called the non-dimensional in-plane coordinate and the in-plane Biot number in this section. sections. here are i - T BL = • (2.125) 50 The theory in this section uses the heat transfer solution for semi-infinite regions and also the Duhamel integral for heat transfer solution from classic heat transfer literature. These concepts are well presented in the book by Ozisik ([22]). 2.8.1 Case of Hold Boundary The two dimensional equation for heat transfer in the solid of interest here is dT (a\d2T ( d 2 T dt 8T_ dT_ di 8T_ dt; d? - + KL2j d? = -Bi(T-Ta) on £ = +\ = +Bi{T-Tj) on £ = -1 = +Bi2{T-T„) on | = 0 = 0 on ^ = oo (2.126) The autoclave temperature Tx and the initial condition T0 are assumed to be constant. It is common in the literature to separate the variables for two dimensional problems. Here we apply this concept to a semi-infinite case and write the temperature lag as: T-T„ T0-Tm = u((Z,t)-u2 (l,t) (2.127) Now we use the following transformation of the variable: T = T„-(T„-T0)-u(£,t)-u2(l,t) (2.128) where u and u2 satisfy the following PDEs: 51 f a \ d2u dt2 du ~dt' du „ . / - 1 — = -Biu on C = +i dt du _. — = +Biu on c=—\ dt « ( ^ , 0 ) = 1 (2.129) and 9m, 9 m. 9 £ 2 *2 _ +Bi2 u2 on £ = 0 - - 0 on 4 ~ 00 (2.130) 9/ 9m. 9w. M 2 ( ^ , 0 ) = 1 Back-substituting (2.128) into (2.126) and simplifying the result using (2.129) and (2.130) we see that the solution to the latter equations combined using (2.128) would give the solution to the two dimensional Equation (2.126). We have already discussed the solution to (2.129) which is demonstrated by the Heisler chart. The equation for the midplane temperature is: ( (0 ,0= £ exp(-4 2 Fo) 2 sin (A,.) i=l,3,5,.. Xi + sin (/l,) cos (/I,) (2.131) On the other hand, the solution for a "semi-infinite region" is also well established in the literature (Ozisik ([22])) and for the case of convective boundary, as is the case in (2.130) would be 52 u2 (l,t) = erf \ K 2 J + exp Bit \ K 2 J J erfc ZNFo + Bi2 Fo \ K 2 J In the case of an isotropic material (K = K2), this simplifies further to (2.132) u2 (l,t) = erf tl4F~o + exp{^Bit)erfc + Bi2 Fo J) (2.133) So, in mathematical terms, the edge effect question would be whether or not u is a valid estimation o f the product uu2. The error involved is the difference between the estimation u and the product uu2 and thus to maintain a certain limit of error or tolerance, err, we must be a certain distance away from the edge which would be 4 Bi,—~,err v J •• minl^ max (u-uuA<err\ . (2.134) This equation is validated using a M A T L A B code and the results are depicted in Figure 2-26 to Figure 2-35. Figure 2-26 shows how g~edge varies as a function of Bi for err = 0.01,0.02,0.05 and K/K2 = 0.01. In Figure 2-27 to Figure 2-33 the value of K/K2 would be 0.05, 0.1, 0.2, 0.5, 1.0, 10, and 100. Figure 2-34 shows the maximum value of ted among all values of Biot number defined as: ec/ge-max = max {tedge} (2.135) 53 The values of K/K2 which are less than one are of special interest in practice as in-plane conductivity of composites is higher. The value of this parameter depends on the layup we are using. Figure 2-35 demonstrates ^ e d for the case of Dirichlet boundary (Bi = oo). 2.8.2 Case of Heat Generation and/or Ramp Boundary In the case of heat generation and ramp boundary conditions, we have the following set of equations to solve: • + d2T d? + H dT _(a\d2T 8t \J})dC — = -Bi(T-Tj) on £ = +l d<; v ' ?L = +Bi{T-Ta) on C = -l ^ - = +Bi2(T-Tx) on | = 0 dT_ = 0 on £ = 0 (2.136) where T=T0+rt , (2.137) and the initial condition is T0 and r is the rate of autoclave heating. Now defining T = Rg + Ht + T0 , (2.138) we can write 54 dg _ d2g + 1 d2g dFo dC ' f K_\ dl2 K K 2 , ^ - = -Bi(g + Fo) on C = +l ^ = +Bi(g + Fo) on ^ = -1 ^ = +Bi2(g + Fo) on £ = 0 5<T = 0 on £ = < g (Fo = 0) = 0 . (2.139) Now using the Duhamel's Theorem, we form the following auxiliary problem: dL 5 2 E 1 5 2 E + -dFo ' { K\ dig2 \ K 2 J = - 5 / ( S + l) on £ = +1 51 — = +5/(2 + 1) on £ = -1 = + f l / 2 ( l + l ) @ | = 0 (2.140) — = 0 @ # = oo 5 | -E ( F o = 0) = 0 Comparing Equation (2.140) with Equations (2.129)-(2.130), one can infer that Z = u u2 -1 . (2.141) or in the case of 5/ 2 = 0.0 E = « - l (2.142) 55 And hence the error is identical as both cases are normalized with respect to a maximum absolute value of one. Namely The solution of Equation (2.139) could be calculated using the auxiliary Equation (2.140) by This equation holds at all spatial points. So i f 2 could be approximated with a maximum tolerance of % at a point for all values of the Fourier number, so does g. Therefore, Figure 2-26 to Figure 2-35 set an upper limit for all types of heating. 2.8.3 Asymmetric Boundaries and Two Layer Solids For the cases with asymmetric boundaries and two layers, the equivalent Biot number presented in the previous sections should be used in combination with Figure 2-26 to Figure 2-35. K/K2, on the other hand, should be evaluated for the layer of interest. 2.9 Methods of Thermal Control In order to avoid complicated three-dimensional modelling for heat transfer analysis of real tool-part setups, one needs to estimate the contribution of spines, legs and fins attached to the convective heat transfer coefficient of tool. This area of engineering is well established in heat transfer literature under thermal control techniques. In this section, the longitudinal fin with a prismatic (rectangular) cross section is investigated. The solution given here is a modification of the typical solutions in the classical heat E - E = w - ww2 (2.143) (2.144) 56 transfer literature (look at [34], for instance). For other thermal control solutions one can refer to the mechanical and heat transfer handbooks like [35]. 2.9.1 The Longitudinal Fin with a Rectangular Profile In this section, we assume that the base of the longitudinal fin (Figure 2-36) is maintained at a temperature Tb =Tb0 + rt where the heating rate of the autoclave is r. The autoclave temperature changes linearly with the rate r as well and is described as T„=T0+rt. The width of the fin is assumed to be A and the thickness is 8 . The temperature is assumed to be constant at each cross section through the width of the fin. Here is the system of equations that should be solved to get the temperature distribution in the fin and hence the relation between Tm-Tb and heat inflow (outflow) a2i A / C - ^ - V = C/zA 2 I l = Tb-T0 on nr = 0.0 (2.145) -k hM on £7 = 1.0 . dm In the above equation, I = T-Tai and A and C are the cross sectional area and circumferential length of the fin. Also, m is the non-dimensional length along the fin width. E7 = 0.0 at the base and m = 1.0 at the tip. The solution to the equation above is: 57 i = (r»-r) cosh 1+ ^ ^ t a n h \ h A — — +tanh h A f lKA~ sinh K A~ (2.146) The rate of heat inflow (outflow) from the fin is equal to the heat entering from (to) the spine at the base. Thus Q = AK = {Tb0-T)AKx h C 1+ £ £ t a n h V h A K A~ KA KC —— + tanh V h A ( K A' (2.147) ]h C For the cases where I——A > 2.6 and therefore tanh \KA~ = 1.0, we have (2.148) For the case of nfm fins in unit length we have Qtotal ~ (^ fcO (2.149) If we are investigating unit length of the fin then A -1 x 8 and C = 2 (1 + S). Therefore, n _ QtotalA " (Tb-Tm) K v h j 1 + I A V 8_j -1 (2.150) 58 Notice that h > h i f — > — — which is the case for tool material. Figure 2-37 eq h l + S illustrates Equation (2.150) and shows the efficiency that could be achieved using a reasonable range of variables. 2.10 Summary In this chapter a closed form solution was developed which provides a quick estimate of the temperature distribution through the thickness of a tool-part system during the process. The solutions are exact for the steady state phase. For the transient phase, the best solution was chosen through mathematical manipulation and numerical trial and error studies. These solutions give results with reasonable error range for heat transfer and constant heat generation for single layer materials. The solution is within acceptable error range for the cases of interest in tool-part set but lacks the generality of the single-layer solution. A few examples in the verification and case study chapter (Chapter 4) show the usefulness of approximate solutions developed here. Here we summarize the simple algorithms presented in this chapter and then we present a few recommendations for the future work. 2.10.1 Summary of the Approximate Solution Algorithms 2.10.1-1 One Layer of Material - Symmetric Boundary Conditions The solution in this case is based on using Figure 2-3 for applied hold temperature and Figure 2-4 for applied ramp temperature and/or heat generation. From these diagrams one can read the midplane temperature based on the Fourier number (Fo), Biot number (Bi) and the R number defined below: 59 at F o = i r Bi = — (2.151) K a The midplane temperatures accuracy depends on the resolution of the diagram used. The temperature distribution in the steady state phase of the temperature ramp is presented in Figure 2-5. This diagram can be used as an approximation for the transient case as well. The midplane temperature in the steady state case is given by T -T 1 0 ^ - ^ ^ = 0 . 5 + — . (2.152) R Bi Figure 2-6 illustrates the approximation involved in using the equation above for a given time (Fo number) by showing the percentage of the steady state achieved at that time. 2.10.1-2 One Layer of Material -Asymmetric Boundary Conditions Two approximate strategies are presented in this section. In both these strategies, the location of the maximum (minimum) temperature (adiabatic line), Sp is located with respect to the midplane surface. The through thickness dimension is non-dimensionalized with respect to L, which is the half thickness of the slab. Here is the non-dimensional location of the adiabatic line J 1_ 8P= B \ B t \ • (2.153) 2 + — + — Bir BiB The top and bottom Biot numbers are defined based on the corresponding h values, i.e.: 60 B K h,.-L (2.154) BiT-^ K In the strategy #1, an equivalent Biot number is calculated from the following equation Bi = -„ 1 1 2 + + — BiT BiB j 10 Bij BiB BiT -Big - + 4 4 1 + -BiT.BiB j (2.155) This equivalent Biot number is used to read the midplane temperature from Figure 2-3 and Figure 2-4. In the strategy #2, an equivalent Biot number is determined by Bi. = max(hB,hj.)-Le K (2.156) where L.=L+Sp . (2.157) Figure 2-5 can be used in approximating the temperature distribution for this case as well. Figure 2-6 can again be used for approximating the level of steady state achieved. 2.10.1-3 Composite-Tool Structure The approximate solution for this problem is based on evaluating an equivalent coefficient of convective heat transfer for the interface between the tool and the composite. Here we reiterate the algorithm from Section 2.5.3: 1) Evaluate BiB,BiT,RvR2,Sx,S2 according to: 61 hjL2 ~K7 (2.158) R> = rL2^ z l V « i J '' L2^ V a2 J (**2-r) (2.159) (2.160) 2) Calculate the interface heat transfer coefficient h from: h = hss = BiBBiTS,S2 (tf, (1 +1 / BiB)-R2(\ + \f BiT )) 5,/?, (1 + BiB ) (1 + 2BiT ) + S2R2 ( l + 5 / r ) ( l + 2BiB ) (2.161) 3) If > 0 then /? is the equivalent heat transfer coefficient for the tool, that is: Bi,,.. = —-(2.162) where Bin is the equivalent Biot number for the upper surface of the first layer. The location of the adiabatic line within the tool would then be 1 1 __Bin BiB S a * 1 1 2 + — + — BiTl BiB (2.163) 62 And so we have Bie=max(BiB,BiTi)-(l + \ta\) Re=Rr(\ + \Q2 (2.164) F o e = f i / ( l + \Ca\) These equations are to be used to read the extremum temperature and the temperature distribution in the tool according to Figure 2-3, Figure 2-4 and Figure 2-5. Note that in the case of temperature hold a fictitious heating rate (like r = 1) could be used to locate the adiabatic line. 4) If h < 0 then -h is the equivalent heat transfer coefficient for the composite, or where BiB1 is the equivalent Biot number for the bottom surface of the second layer. The location of the adiabatic line in this layer would be BiT Bit (2.166) ^ 1 1 2 + + — Bij Bil 'B2 And thus we have Bie =max(BiB2,BiT)x(\ + \g-a\) Re=R2x{l+\l\f (2.167) 63 which are to be used to read the extremum temperature and the temperature distribution in the composite according to Figure 2-3, Figure 2-4 and Figure 2-5. Again, in the case of temperature hold a fictitious heating rate should be used to locate the adiabatic line. 2.10.1-4 The Edge Effect In order to evaluate the accuracy of our solutions close to the edges of the slabs we introduced a couple of parameters that illustrate the effect of the second dimension: H M (2.168) Figure 2-26 to Figure 2-33 illustrates the minimum relative distance from the edge i^edge)t0 achieve minimum relative errors equal to x - 0.01,0.02,0.05 in estimating the midplane temperature. These are drawn as a function of Bi for different ratios of through thickness to surface direction conductivities (K/K 2 ). Figure 2-34 shows the maximum value of ged among all values of Biot number and Figure 2-35 demonstrates ged for the case of Dirichlet boundary (Bi = <x>). 2.10.1-5 Equivalent h for a Number of Longitudinal Fins with Rrectangular Profiles Figure 2-37 demonstrates the ratio of equivalent convective heat transfer coefficient to the actual heat transfer coefficient (h /h) on a surface with nfm number of fins each with a thickness of 5. This could be used in estimating the increase in the heat transfer coefficient as a result of using fins. 64 2.10.2 Future Development In order to make the diagrams and closed form solutions even more useful in the process modeling of composite materials a few further steps need to be taken. One is setting up a general strategy for using the simple diagrams developed here and also the Heisler charts for multiple layers. Another important issue to be considered in future research is the change in thermal properties of the solid through the analysis as the charts we have developed here are for constant material properties. A study on the material properties evolution could lead to some modification factors to the results obtained from the diagrams or may set limits on the usefulness of these diagrams. 65 Table 2-1 - Thermal properties for materials of interest. /7(kg/m J) C p ( J / k g K ) kz ( W / m K ) az (ml 1 s) I N V A R 8000 515 11.0 2.67 x lO ' 6 C O M P O S I T E 1580 870 0.69 0.50 x 10"6 A L U M I N U M 2710 896 167 68.9 x 10' 6 STEEL 7860 465 51.9 14.2 x 10 - 6 66 Time Figure 2-1 - Typical autoclave process cycle. Fluid at T^hj. i + L . + 1 - L . - l Fluid at T a , h B Figure 2-2 - Schematic representation of one dimensional transient heat conduction analysis with asymmetric boundaries. 67 1.E-01 1.E+00 1.E+01 1.E+02 Fo 1.E+03 1.E+04 Figure 2-3 — Heisler chart for the temperature in the middle of a slab versus the Fourier number for applied constant temperature at the boundaries. Contours are for various Biot numbers. 1.E+03 1.E+00 1.E+01 1.E+02 1.E+03 Fo 1.E+04 Figure 2-4 - The charts for the temperature in the middle of a slab versus the Fourier number for applied temperature ramp at the boundaries. Contours are for various Biot numbers. 68 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 X=£-8p Figure 2-5 - Temperature distribution through the thickness of the slab in the steady state phase. Contours are for various Biot numbers. 1.E+04 T 1.E-01 A 1.E-02 -\ , , , , , 50 60 70 80 90 100 % of Steady State Achieved Figure 2-6 - The Fourier number at which a percentage of the steady condition state is acievedfor the case of applied ramp temperature at the boundaries. Contours are for various Biot numbers. 69 0.07 T 1 10 100 1000 10000 100000 hT/hB Figure 2-7 — Maximum error of strategy #1 in transient phase for the ramp heating of the autoclave; the error is presented as a function of BiB and hj. I hB. Figure 2-8 - The global maximum errors extracted from curves in Figure 2-7 plotted as a function of BiB. 70 0.09 t 1 10 100 1000 10000 100000 hT/hB Figure 2-9 - Maximum error of strategy #1 in transient phase for the case of temperature hold in the autoclave; the error is presented as a function of BiB and hj I hB. Figure 2-10 - The global maximum of curves in Figure 2-9 as a function of Bit 71 h T/h B Figure 2-11 - A as a function of hj. I hB for 0.001 < BiB < 0.5. . • 10.5 | 1 B i s \ j s ho l i -inn i 1 1 1 , w . 1 1 1 10 100 1000 10000 100000 h T/h B Figure 2-12 - A as a function of hjlhB for 0.5 < BiB < 100. 72 0.06 T hT/hB Figure 2-14 - Maximum error of strategy #2 in transient phase for autoclave temperature ramp; the error is presented in terms of h, I hB for 0.001 < BL < 0.5. 73 Figure 2-15 — Maximum error of strategy #2 in transient phase for autoclave temperature ramp; the error is presented in terms of hj. I hB for 0.5<2?i B<100. 0.045 n Figure 2-16 - The global maximum of curves in Figure 2-15 and Figure 2-17 versus 74 Figure 2-17 - Maximum error of strategy #2 in transient phase for the case of temperature hold in the autoclave; the error is presented in terms of h, I hL for 0.001 < BiB < 0.5. 0.05 -i 1 10 100 1000 10000 100000 hT/hB Figure 2-18 - Maximum error of strategy #2 in transient phase for the case of temperature hold in the autoclave; the error is presented in terms of hr I ht for 0.5 < BiB <100. 75 0.05 i 0.04 0.03 0.02 4 0.01 0.001 0.01 10 100 Figure 2-19 - The global maximum of curves in Figure 2-17 and Figure 2-18 versus BiB • autoclave . Tm,hT L e 1 K2, oc2,H2 / ) layer II: part layer I: tool T Kj, aj y ' z 2 L 2 ± , = . 2 L j ± T h autoclave: 0 0' B Figure 2-20 - Typical two-layer composite studied in this paper. 76 Figure 2-21 — Absolute values of the maximum errors incurred in estimating the temperature at the adiabatic line for the case of applied temperature ramp and for cases when the adiabatic line falls in the second layer (composite part). o • Temperature Lag • Autoclave T + 140 ~ -- 120 *•£ -- 100 I 3 -- 80 H -- 60 -- 40 20 Terr(C) Figure 2-22 - Absolute values of the maximum errors incurred in estimating the temperature at the adiabatic line for the case of applied temperature ramp andfor cases when the adiabatic line falls in the first layer (tool). 11 0.1 0.02 0.04 0.06 0.08 W A T , , old Figure 2-23 - Maximum relative errors incurred in estimating the temperature at the adiabatic line for the case of applied temperature hold and for cases when the adiabatic line falls in the second layer (composite part). Figure 2-24 - Maximum relative errors incurred in estimating the temperature at the adiabatic line for the case of applied temperature hold andfor cases when the adiabatic line falls in the first layer (tool). 78 F l u i d a t F l u i d a t Tx,h . +L,+1 1 - L , - l 1 F l u i d a t T„9h Figure 2-25 - Schematic of the model used to study the edge effect. -/ err=0.0 / t H ) 2 1 / t H ) 5 i I 1 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 Bi Figure 2-26 - 4edge os a function of Biot number for different error tolerances and for K/K2 =0.01. 79 r r r r ferr-QX 7 1 > r 0^02 •—^ 1 1 1 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 Bi Figure 2-27— <^edge as a function of Biot number for different error tolerances andfor K/K2 = 0.05. 25 ~| 1 1 i i i i 20 -15 -o> Ol •a o 10 -5 -0 -1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 Bi Figure 2-28 - ^edge as a function of Biot number for different error tolerances andfor K/K2 =0.1. 80 25 T •D O 10 0 r i r r i — - I— r _ ,— -/err r =0.01 ^ i r r 1 1 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 Bi Figure 2-29 - 4edg(! as a function of Biot number for different error tolerances andfor K/K2 = 0.2. 25 - r 20 15 4 10 -r i — r r i -r r /ferr=0.01 / /om * r ~ ~ r ~ ~ i * — i — 0.05 1 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 Bi Figure 2-30 - 4edge as a function of Biot number for different error tolerances and for K/K2 = 0.5. 81 Figure 2-31 - gedge as a function of Biot number for different error tolerances and for K/K2=\.0. <err=OM /0.1 32 \ 1 —1 1 1 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 Bi Figure 2-32 - g~edge as a function of Biot number for different error tolerances andfor K/K2=\0. 82 Figure 2-33 - 4edge as a function of Biot number for different error tolerances and for K/K2 =100. 25 -t 0 -I 1 \ i 1 0.01 0.1 1 10 100 K / K 2 Figure 2-34 - g~eilge_m!0i as a function of K/K2 for different error tolerances. Markers show the maxima from the previous diagrams and trend lines are drawn through them. 83 Figure 2-35 - ^edge as a function of K/ K2 for different error tolerances for the case of Dirichlet boundaries (Bi = °o ). P A R T T O O L Figure 2-36 - A typical longitudinal fin used to increase the efficient convective heat transfer coefficient. 84 85 3 DEVELOPMENT OF A HIGHER ORDER SHELL E L E M E N T FOR T H E R M O C H E M I C A L ANALYSIS Nomenclature l U I ^ energy norm for volume V «(. , .) typical bilinear form A m , B m , F Z m Parameters in spatial integrations for K Bi Biot number c degree of cure Cp thermal capacity (constant pressure) C capacitance matrix dA differential surface area da differential surface area in parent system dl differential length vector in parent system dL differential length vector dV differential volume dTl differential volume in parent system DKGP thickness of the shell element at each surface integration point erredge estimation of the error at he edge of the shell element / parameter in 2D surface P D E of slab (multilayer) / (C) typical oscillating function / (c, 7/) cure kinetics function F, integral of f(£) over F the force vector for nonlinear procedures F(T) the nonlinear function to be solved for temperature in the next time step g {j; [, rj, typical smooth function 86 GF geometric shape function h coefficient of convective heat transfer hx element size in x-direction h element size in y-direction H exothermic heat generation per unit volume HR resin heat of reaction H force vector H the contribution of convective heat transfer to force vector J Jacobian of transformed coordinate k, k conductivity, conductivity tensor K stiffness matrix K contribution of convective heat rate to stiffness matrix kt, k2 conductivity on both side of interface KIP integration point index L (.) typical linear functional L half thickness of a slab L L lower bound for error Ly upper bound for error L 2 (Q) space of functions over Q equipped with L 2 norm M number of Fourier terms M surface shape functions matrix nT number of trigonometric terms N number of layers N interpolation shape function NG surface geometric shape functions nGP number of Gaussian points in surface nGP_g number of Gaussian points in £ direction nGP-ri number of Gaussian points in 77 direction nIF number of interpolating shape functions 87 nlp number of integration points nT number of trigonometric terms O order of the error q,q heat flux, heat flux vector Q exothermal hear source per unit volume Qs heat flow through the boundary smf (J) surface modification function, defined by Equation (3.64) S surface of the slab - General surface of a three dimensional geometry ST surface with Dirichlet boundary condition S surface with Neumann boundary condition S Jacobian stiffness matrix for nonlinear procedures t time tmnge estimation of the time range T temperature distribution function Trange estimation of the temperature range Tm fluid (autoclave) temperature T G through thickness geometric shape functions u P D E variable with homogeneous boundary condition U internal energy per unit volume v typical distribution function w weight function in Galerkin method W surface P D E variable, after z is separated out W in-plane component of the temperature wG weight function for Gauss quadrature {x,y,z) local coordinate system (X,Y,Z) global coordinate system Y weighted internal energy integral Z ( ^ ) trigonometric shape function, depends on the layers and boundary zf (<^ ) through thickness shape function 88 Greek symbols {4, n,t) parent coordinate system tj through thickness integration points t, through thickness discretization points a diffusivity J3 constant (in multilayer section) 6 interpolation parameter in the 6 -method 8 constant (in multilayer section) <T Stephan-Boltzman constant s emissivity of the surface vR resin volume ratio p density time step in transient analysis A transformation matrix YD modal capacitance matrix c2 modal force vector ¥ time integration of H time integration of H S typical integral n volume in the parent space Q surface area in the parent space fi? surface area in the parent space with Neumann boundary condition Subscripts b boundary B bottom surface c composite i initial value P perturbed value of the parameter R resin T top surface 89 Superscripts approximate value of the parameter + - top and bottom of the interface [/'] layer number ini initial IP the value at the integration point (j) mode number 90 3.1 Introduction The shell finite element solution presented in this chapter is an F E M type of solution that uses the analytical solution of the heat transfer problem to enhance the computational effort, both in terms of meshing and analysis. This shell element is particularly developed for modelling the thermochemical processes in a composite shell during manufacturing. The main assumption in producing this element is that the through thickness direction is the more critical direction, i.e. it has the more complicated distribution function for temperature. For surface direction, multi-linear functions are considered to be good approximation of the distribution. As a result, exact distribution of temperature through the thickness of the shell is investigated and the resulting shape functions are embedded through the thickness to enhance the F E M solution. In addition to using a shell model with hierarchical (Finite Strip) shape functions, multi-scale time stepping and efficient spatial integrations are also used to decrease the solution time. This chapter first describes the constitutive equations that lead to the boundary and initial value partial differential equations. In the next step the Hierarchical Fourier shape functions used in the element are presented. Afterwards, the shape functions are embedded in a general 3D shell element, which allows stacking of elements and thickness variations in an element. Next section describes the challenge in integrating the analytical functions through the thickness and proposes an innovative method to overcome the challenge by performing integrations for each integration point in the first step and using the matrix of solutions throughout the transient analysis. 91 In the next stage, a time integration solution is introduced, which separates the rather quick thermochemical phenomenon from spatial heat transfer. This solution saves a large amount of computational time by avoiding huge matrix inversions in every thermochemical time step. A couple of alternative nonlinear procedures for solving the nonlinear problem are presented next and finally the algorithm that is prepared to add the shell element to the commercial software, A B A Q U S , in the form of a user element (UEL) is presented. 3.2 Constitutive Equations This section presents the general constitutive equation for heat transfer and diffusion in an orthotropic solid. This general formulation is based on the energy balance of a closed system as is stated in heat transfer and thermodynamics text-books, e.g. Cengel and Turner [36]. Then a formulation of exothermic heat generation in thermoset matrix based composites in combination with heat transfer and diffusion is presented ([4]). 3.2.1 Heat Conduction and Diffusion in a Solid Consider a solid of a general shape as is shown in Figure 3-1. Heat transfer in this solid follows the classical Fourier's law of heat transfer, which in one dimensional can be written as ax where q in the equation above is the heat flux, which is the amount of heat transferred per unit area and k is the conductivity. 92 Equation (3.1) can be expressed in three dimensions as q = - k V r (3.2) where k is the conductivity tensor defined by k = kyy kyz kzx kzy k2 (3.3) If conductivity is given in a local coordinate system, we can calculate the global matrix from the following equation k = A k t A r , (3.4) where kL is the conductivity tensor in local coordinates and A is the transformation matrix that is defined by A = X.x X.y X.z Y.x Y.y Y.z Z.x Z.y Z.z (3.5) Here ( X , 7 , Z ) is the global coordinate system and (x,y,z) is the local coordinate system at the point of interest. Hat sign indicates unit vector in the corresponding coordinates. Simple heat equilibrium within the solid of Figure 3-1 is as follows: j^dV = Qs + JQdV (3.6) 93 where U is the internal energy per unit volume and Q is the heat source per unit volume. Also Qs represents heat flow through the boundary and is defined by: Q,=-\qadA = -\v.qdV , (3.7) A V where n is the unit normal to the boundary of interest. The second equality is a result of the Divergence Theorem. Substituting Equation (3.7) into Equation (3.6) we obtain: j [ ^ + V . q - e V = 0 . (3.8) As V is arbitrary, we have dU = - V . q + g . (3.9) dt Then substituting q from Equation (3.2) results in the following equation: dU dt The rate of change of U with temperature is given by dU V ( k V r ) + g . (3.10) = pCp , (3.11) dT where p is the density and Cp is the specific heat capacity under constant pressure. If U is only a function of temperature then we have 94 which, combining with Equation (3.10) leads to the following equation: p C P ^ = V{kVT)+Q (3.13) Or in the case of homogeneous isotropic material 1 3 T - . V T + Q , (3.14) a dt k where a = k/pCP (3.15) is the thermal diffusivity of the material. To solve the above P D E , initial conditions in the following form need to be defined: T(x,y,z,0) = T0(x,y,z) (3.16) Boundary conditions, on the other hand have one of the following formats: T(xb,yb,zb,t) = Tb(t) (3.17) - * ( ? - ) » = ft (0 (3-18) on The above equations are called the boundary conditions of the first kind (Dirichlet) and second kind (Neumann). qB in the above equations is the heat flow out of the surface. The boundary conditions are called homogeneous when the right hand side is equal to zero. Another possible way that heat could be transferred to the surface of a material is through convection which is called boundary condition of the third kind and is formulated as following 95 - A = W r i ) / ) - r . ) ) (3.19) on where h is the convective heat transfer coefficient of the surface and Tx is the temperature of the ambient fluid. This kind of boundary condition is called homogeneous if Tx=0. Other possible boundary conditions are radiation boundary condition: -k(~)b=as[T\rb,t)-Te4] (3.20) and interface boundary condition as defined below: -k-(^-\=-k+(^f-)b (3.21) on on where 8 is the emissivity of the surface and a is the Stefan-Boltzman constant. ky and k2 in Equation (3.21) are the thermal conductivities on both sides of the interface. 3.2.2 Constitutive Equation for Exothermic Heat Generation The constitutive equation of interest for heat transfer in the processing of composites is \p p X -zr=-v-q+vRpRHR— ot at p 2 2 ^ where c is the degree of cure of the resin and represents the advancement of chemical conversion and f{c,T) is a function that describes cure kinetics. (pCP) denotes the average (smeared value) of pCp for the composite while pR is the resin density and v R is resin volume fraction in the composite. The second equation (cure kinetics) however 96 evolves more quickly and therefore needs smaller time increments to track it compared to the first equation which describes the heat transfer. The solution procedure that accounts for this difference in the evolution of T and c is described in the section on time integration (Section 3.6). Different models that can be used for / are described in the literature, e.g. Johnston ([4]). 3.3 Parallelepiped Shell Element Formulation In this section, a set of Fourier shape functions are developed that describe the through thickness distribution of temperature in a shell element. The idea of using non-conventional orthogonal shape functions in Finite Elements has been investigated in two branches. One is the finite strip method where the result of analytical (series) solution to the plate problem in surface direction is used as a shape function in one surface direction to enhance the computation. This method has been widely used for stress analysis of plates by Cheung et al and Zienkiewicz ([37],[38],[39]). Other researchers have used the procedure in modelling heat transfer problems in slabs ([40], [41]). Another closely related concept in the literature is the use of hierarchical shape functions. These functions are usually orthogonal polynomial interpolation functions that are used to interpolate the values between nodal values ([42],[43],[44]). The major usefulness and challenges of using these shape functions are the same as the ones that are encountered in dealing with Fourier shape functions. There is also a vastly developing literature on using this kind of shape functions for multilayer composites. Examples are the works by Babuska ([45]), Surana ([46],[47]) and more recent works by Actis ([48]), Duster ([49]), and Rolfes ([50]). There are also ongoing investigations on error analysis, global-local methods and 97 iterative methods using the hierarchical methods in multilayer composites. A few examples are the works by Fish (e.g. [51]), Reddy ([52]) and Frutos ([53]). In this section and in A P P E N D I X A , Fourier series shape functions are developed for the through thickness analysis of composite shells based on different assumptions on the number of layers and top and bottom boundary types. For practical reasons, only one elements type with embedded shape functions for Dirichlet-Dirichlet boundaries is used for the codes developed in this work and therefore wil l be presented here. Other boundary conditions are presented in A P P E N D I X A . The method that is used here for solution of the heat transfer P D E and determining the shape functions is separation of variables and series expansion in the through thickness direction as is well established in the P D E literature (see for example Zauderer ([54])). Similar techniques and solutions can also be found in the heat transfer books by Holman ([11]) and Ozisik ([22]). 3.3.1 Dirichlet Boundary Conditions Figure 3-2 illustrates the cross section of a typical shell element. Note that (x,y,z) are the local coordinates and the origin is at the centroid of the parallelepiped. The thickness of the parallelepiped element is uniform and equal to 2L . The temperatures on top and bottom surfaces are specified values and the material is assumed to be transversely isotropic with the principal axes parallel to the local element coordinates. The following general 3D formulation holds (conductivities are assumed to be spatially independent): T=T(x,y,z,t) (3.23) 98 pCPT, = kxT + kT + kT + H(x, y, z, t) T(x,y,+L,t) = TT(x,y,t) T(x,y,-L,t) = TB(x,y,t) T(x,y,z,0) = Ti(x,y, z) The P D E above can be rewritten in the following format: T, = aT + aT + aT + (1 / pCP )H(x, y, z, t) (3.24) (3.25) (3.26) (3.27) (3.28) where a x , ay and az are diffusivities in the corresponding directions. In order to change the boundary conditions into homogeneous boundary conditions we apply the following transformation in the variable T s T-r ~T" TQ T-p TB Z T(x,y,z,t) = (^-^ + ^ 2 L ) + u(x,y,z,t) (3.29) U,i = + ayU,yy + UzU,zz ) + 0 1 PCp J V , Z > 0 ^2 + 2L' }__]_ z_. "2 2 V <TBJ -axTBxx )(----) (3.30) u (x ,y ,-L,t) =u [x ,y ,+L,t) = 0 i(x,y,z,0) = Ti(x,y,z)-fT.,.+TB Tr-T0 T ' x f l | T *B 2 L (3.31) (3.32) In order to find the proper shape functions through thickness, the following eigenvalue problem is considered 99 u. — a u + oc u +a u ,i x ,xx y ,yy z ,zz u (jc ,y ,-L,t) =u (x ,y ,+L,t) = 0 (3.33) (3.34) Now consider the following solution based on separation of variables: u(x,y,z,t)=W (x,y,t)Z(z) (3.35) Substituting into Equation(3.33) results in W , , Z = a * W , * * Z + a y W , y y Z + a , W Z " Z(-L) = Z(L) = 0 (3.36) (3.37) Rearranging the terms in the above differential equation leads to w w y w z (3.38) The practical nontrivial solution to the above is when z V ' (3.39) where P is a constant. The solution to the above O D E is Z =A sin(/5— ) + B cos(/5—) . L L (3.40) Applying the boundary conditions in Equation (3.37) we obtain the following set of solutions for Z: cos ( , n z ^ J 2 L sin . n z KJ2Ly 100 j is odd j is even (3.41) And the general solution for u can be written in the following series format: u=^W{J\x,y,t)Z(J\z) 7=1 (3.42) Substituting the above in Equation (3.30) for each j, integrating z from - L to +L and using the orthogonality property of the trigonometric Zy yields W}J) = axWW+ayWV>-az (^-)2W U)+h(j)(x ,y ,t) (3.43) h(j\x,y,t) = -J (1 / pCp )H(x,y,z, t)Z(i) (z)dz +L \zu\z)Zu\z)dz \ \ { T T t - a J ^ - a ^ ^ + ^ + iT^-aJ^-ayTB^-^)jZ{j)(z)d; _ jZ{j\z)Z(j\z)dz 1 1 z , The above P D E is to be solved using the following initial conditions: (3.44) W°\x,y,0) = TT+TB TT TB Z 2 d Z°\z)dz Jt=0, +L \z(i\z)Z(i)(z)dz (3.45) The solution of the above problem can be performed using Finite Element Method (FEM) discretization for each W}, i.e. W^{x,y,t) = ±Nn(x,y)W^{t) (3.46) 101 where N JF is the number of interpolation shape functions in each element and [N n} N IF is the set of interpolation shape functions of the element. In order to apply the Finite Element Method, we need the weak formulation of Equation (3.43), which using the weight function w and integrating by parts wi l l turn out to be = jwhU)(x,y,t)dS , s where w is the weight function and S denotes the in-plane surface area of the element. Using the Galerkin F E M discretization, the shape functions should replace the weight function and applying Equation (3.47) to each element we have 2 jwWJJ)dS + ax \wxW[^dS + ay jwyW(^dS + az s s s (3.47) (3.48) where, (3-49) s (3.50) (3.51) s 102 3.4 General 3D curved shell element formulation with stacking and tapering capabilities The through thickness functions in the previous section ( Z functions) are used in a general shell element in this section to make the element more practical for applications involving complex geometries. In this section, first the geometrical shape functions that are used to approximate the shell surfaces are introduced, then the Z functions are incorporated into the general interpolation functions which are linear in plane and sinusoidal in the through thickness direction. Both these shape functions are then used in a weak formulation of the problem in the next step. It should be noted that the basic shell element here has 20 nodes (8 nodes on each of top and bottom surfaces and one on each one of the surrounding edges). However, only 16 of these nodes are used in describing the geometry through geometrical shape functions as are shown in Figure 3-3. The temperature nodes are the 8 nodes on the vertices and also 4 nodes on the 4 corners which carry the trigonometric terms. These 4 nodes don't have a location and therefore are not shown on Figure 3-3. 3.4.1 Geometrical Shape Functions The geometric shape functions we are using for the shell element introduced here are 8-noded in plane and linear in the through thickness direction so that they can describe surface curvatures (Figure 3-3). The shape functions for the geometry are named and defined so that the following mapping relations hold: 103 "GF x = £ G F , ( £ , , 7 , C ) x , ;=1 "GF ^ = £ G F , ( # , J 7 , 0 y ( ;=1 "OF z = £ G F , . ( ^ / 7 , C ) z , (3.52) Here nGF is the number of geometrical shape functions used in defining the element geometry. (£,77,^) are the coordinates in the parent system and (x,y,z) is the global coordinate system. Now we define the Jacobian J as follows: J = dx/ dx/ dx/ /dt /drj M dy/ dy/ dy/ /dg- /dn /d£ dz/ dz/ dz/ /dt; /dn /dZ (3.53) Substituting Equation (3.52) into Equation (3.53) leads to the following equation for the Jacobian 1=1 y, -V z. 1 1 J {GF, (3.54) J is the ratio between differential vectors in the physical and parent spaces anddet(j) is the ratio between a differential volume in the physical space and the corresponding volume in the parent space. In other words dL = Jdl dV = det(3)dU (3.55) 104 where dL = [dx,dy,dz}T and dV = dx.dy.dz are the differential vector and volume in the physical space and dl = {d<!;,dr],d£}T and dTl = d^-d-q-dC, are the analogous notions in the parent space. In order to get the relation between surfaces we manipulate the above equation as follows: </F = det(j)rfn (dA.dL) = det(j)(da.dl) (3.56) (dA.(Jdl)) = det(j)(da.dl) , dA and da are the differential surface area vectors in the physical and parent spaces. Note that d A. (J dl) = dA, (Jvdlj) = (jydA,) dl} = (j],dA,) dl} = (J r d A).dl det (J) (da.dl) = (det (J) da) .dl . The second equation above holds as det(j) is a scalar. Therefore, further matrix manipulation on Equation (3.56) leads to (j T dA).dl = (det(j)da).dl for arbitrary dl , (3.58) And therefore J T dA = det(j)da dA = det(j)J _ T da (3.59) In the case of our shell element we are interested in surface integrals on top and bottom surfaces as they are the major sources of convective heat transfer. For these surfaces we have 105 da = {0 0 ±dQ}T =±dQ{0 0 l } r , (3.60) which, using Equation (3.59) leads to dA = ± ( d e t ( j ) J Q ) j - T { 0 0 l } r . (3.61) Further matrix manipulation on Equation (3.61) leads to dA = ± (de t ( j ) JQ ) { lNVJ 3 1 INVJ 3 2 INVJ 3 3 } r , (3.62) and taking the absolute value leads to dS = |dA| = ( d e t ( j ) 7 l N V J 3 1 2 +INVJ 3 2 2 +INVJ 3 3 2 J d Q . , (3.63) where INVJ = J" 1 and dQ. is a typical element of surface area on top or bottom surface of the parent element. dS is the corresponding differential surface in the current coordinates. We define the following notion for future use: smf(J) = d e t ( J ) V l N V J 3 1 2 + I N V J 3 2 2 + I N V J 3 3 2 , (3.64) which is the ration between the physical and parent surface areas for top and bottom surfaces. Similarly, for edge surfaces of the element we obtain smf4 (J) = det (J) ^/iNVJn 2 +INVJ 1 2 2 +INVJ 1 3 2 smfv (J) = det (J) A / l N V J 2 1 2 +INVJ 2 2 2 +INVJ which are used as ratios for surfaces normal to dE, and dn . 3.4.2 Interpolation Functions The interpolation functions for the temperature are defined so that: 106 2 23 (3.65) "IF T^N,T, , (3.66) 1=1 where T (scalar) is the spatial temperature distribution function, T = {7]} is the vector of nodal values and nIF is the number of interpolation functions. These interpolation functions are defined within an element as follows: N^M^f^), / = m n = 1 4 + m n = nT+2 , (3.67) 7 + n + (m-\)-nT n = 2,3,...(nr+l) where m = 1,2,3,4 refers to the element edges and nT is the number of trigonometric terms used in defining the through thickness distribution of temperature. Equation (3.67) is our book-keeping equation for degrees of freedom of the element and shows that we first count the bottom surface degrees, then the top surface and at the end we count the trigonometric terms edge by edge. Mm in Equation (3.67) is the surface shape function and for the case of bilinear shape function is defined as follows: M , = l / 4 ( l - < f ) ( l - / 7 ) M 2 = l / 4 ( l + £ ) ( l -77) v M ; (3.68) A / 3 = l / 4 ( l + £ ) ( l + 7) M 4 = l / 4 ( l - < r ) ( l + /7), Here zfn (< )^ is the through thickness shape function and for the case of assumed Dirichlet boundaries, it is defined as follows: 107 2/j(0 = l / 2 ( l - 0 ^ r + 2 ( 0 = l/2(l + 0 (3.69) z/;(C) = 2 ( 0 (O / = 2,3,...,/v+l. Z ( ( ) ( 0 is defined as (Equation (3.41)) 7C cos(i—£) i is odd 2 . (3.70) sin(/'—C) i is even 2 The number of Z ^ (£") used in one step or increment of a transient analysis depends on the accuracy one needs in acquiring temperature distribution. 3.4.3 The Weak Form Applying the Galerkin weighted integral to Equation (3.9) we have Jw dV = -$wV.qdV+ jwQdV , (3.71) v dt v V where V in the equation above is the volume of the system and w is the weight function. Integrating Equation (3.71) by parts results in jw—dV-jVw.qdV = -\wqndS+\wQdV , (3.72) where S is the boundary surface and qn is the normal component of the flux to the boundary. The boundary surface could be broken down into S = ST+Sq , (3.73) 108 where ST and S are the surfaces with Dirichlet and Neumann boundaries correspondingly. Equations (3.71) and (3.72) hold for all permissible w's that are set to zero on Dirichlet boundaries and second integral on the left hand side of Equation (3.72) can be rewritten in the following form f w j r , r\ dw dw dw dV (3.74) The components of q are defined by Equation (3.2). We approximate the temperature using interpolation shape functions as in Equation (3.66) and use the same shape functions as the weight function, as follows: "IF w = Nj j = l,2,...n IF (3.75) Substituting Equation (3.75) into (3.2) leads to q = -£(kVJV,)i; 1=1 (3.76) Substituting (3.75) and (3.76) into (3.72) yields the following equation: dU ^ \Nj—dV + £ j ( k W , ) . V A ^ T^-JN^dS+JNjQdV , (3.77) \v 0* J 1=1 \V ) Sa V which results in "IF . Yj + YjK.Ti=-\N]qndS + H] (3.78) 109 where (3.79) where AU is the change in internal energy from an initial state. On the boundary surface, the general Neumann boundary is defined as follows: qn =-<io+hT Therefore the right hand side of Equation (3.78) can be rewritten as (3.80) "IF -\Njq„dS= \NjqodS-X jhN^jdS "IF (3.81) Finally, Equation (3.78) leads to the following equation in matrix format, which should be solved using a finite difference solver in time Y + ( K + K ) T = H + H . (3.82) Y is evaluated for the case where U is only a function of temperature, as follows: 1 fit dt J 1 i J Pit • J J \ AT dt dT dT dT KdT, dt dV , (3.83) which results in rj=l\lNMpCp)w}tl = flcytl (3.84) The following are the definitions of the above tensors in the physical space and also in the parent space 110 v n Kv = l(kVN,).VNjdV= J ( k J - r V A ^ , ) . ( j - r V i V y ) d e t ( j ) dU v n Kv = jhNtNj dS= \h N,Nj smf(3) dQ (3.85) s, a, Hj = JNjQdV = JNjQdet(j)dU . v n HJ=\q0NjdS= \q0NjSmf(J)dn s, a. Here, Q ? refers to the boundary with Neumann boundary condition in the parent space. These integrals are evaluated several times during the transient analysis as the material properties are generally varying throughout the analysis. In the next section, techniques that are used in evaluating these integrals are presented. 3.5 Efficient Spatial Integration Evaluating the integrals over the domains of interest might become extremely time-consuming i f they are not performed in an efficient way. Specifically in all the noted integrals, we are interested in evaluating the following types of integrals: E=\f(£)g(C)dC . (3.86) 7T f(C) is an oscillating analytical function like / ( £ " ) = sin(/—£0 (especially when while g(C) is a rather smooth function like a lower order polynomial. This type of combination occurs in all integrals in (3.85) and we wi l l further discuss the cases of interest later on in this section. I l l Finding efficient solutions to the above type of integrals is of great interest to researchers in the field of applied mathematics [55]. Some researchers are working on quadrature methods of approximation, e.g. Patterson [56] and Evans [57],[58]. Another class of efficient solutions are referred to as "Asymptotic Expansions of Integrals", such as the method of stationary phase [59]. The asymptotic method would need successive derivatives of the function g(£) which is not generally available. Another common approximation to Fourier coefficients is the Fast Fourier Transform [60], which using the same number of integration points gives less accuracy compared with the more involved quadrature methods. In this work, a particular form of the quadrature method is developed as it provides a rather efficient evaluation of the integrals. This methodology helps avoid large number of integration points in the cases where / (<^ ) is a known analytical function with known closed form solution for the integral. This section also gives details on evaluation of the integrals in (3.85) using the proposed quadrature method. Let the set of integration points be {Ci}"-i w n e r e n j p is the number of integration points through the thickness. The through thickness range in the parent space is [-1,+1] and it is discretized to in a way that the discretization points are placed half-way between two successive integration points. Therefore, the following relation can relate 112 C = ( C + 0 / 2 i = 2,3,...,nIP+\ (3.87) This is schematically presented in Figure 3-5 for the case of three integration points. The above discretization and set of integration points are used to evaluate the integral as follows: s = £ J / ( O s ( 0 ^ = Z J / ( O s ( £ K • (3-88) This simplification is accurate enough if the through thickness discretization is fine for function g so that its value is almost constant in the range • We can further manipulate Equation (3.88) to get S s 2 > ( £ ) jf(t)dt = Zg(£,)F< . (3-89) where ?M \_f{£)dC . (3.90) <r, Here Ft 's are dependent on the number of integration points used in the analysis and can be evaluated a priori. As a matter of fact Ff and other parameters that are not dependent on the element shape (geometry) and its material properties can be evaluated before the start of the transient analysis and saved as data-structures related to the element type. These efficient data structures wi l l be used in the evaluation of integrals during the transient analysis. This type of solution saves in time performing integration in nonlinear 113 problems where the material properties change with time and temperature and hence the stiffness and capacitance matrices should be re-evaluated several times throughout the analysis. It should also be noted that the above solution is consistent in the sense that it can also be used in cases where / (g) is rather smooth in which case the result is close to trapezoidal integration schemes. The practical application of the above integration scheme is shown here for the first and second integrals of (3.85). The rest will follow suit. For the first equation in (3.85) we have C , = fcpCp)NlNJdV= J(/>C,)tf,tfydet(J)dri v n C v = T J 7 (pCp)(M^MH)^{3)(zf^zfH)d4dr,d^ (3.91) *tf.</.0 AO "pp "IP Cjj ~ £ £ WGKGP 8(.4KGP'T1KGP' <?K1P) FK1P ' KGP=\ KIP=\ where F is determined by integrating / as given by Equation (3.90). nGP is the number of Gaussian integration points in the surface integration, which is the preferred method of integration in the surface. wG is the corresponding weight function. Also note that the parameters (w,,^),^^,^) are related to i,j according to the relation in the Equation (3.67). The surface parameters in the Gauss method are calculated from the parameters using the following equations: HGP ~ NGP-i-NGP-N (3 92) wGKGP = wGKGP_rwGKGP_n 114 In the above equations, £,77 denote the variable in the corresponding direction. The second equation in (3.85) is now modified to: K . = l(kVN,).VNjdV = j ( k J " r V i V / ) . ( j - r V i V ; ) d e t ( J ) J n v n k J r \MmJzni Mmfz' C, =+1 n=+l cf=+l * I { J <T=-I <7=-l f=-l f = + 1 ^=+1 <f=+l k J X * 0 X , 0 0 M „ Mm]fz'nj \ r Now we define the following tensors: d e t ( j ) ^ ^ C fa \K J J 0 M 0 i"j,y 0 det(j) d^dnd^ , i = \,2,3,...n,F m = \,2,3,A (3.93) 0 " A m = k . T r 0 0 Mm _ 0 " B m = J r 0 0 m _ F Z " " U J (3.94) (3.95) (3.96) Equation (3.93) can now be written in a more compact form as 115 f=+17=+l f = + 1 K</ = I I j ( A m < F Z " ' ) . ( B m ' F Z " ' ) d e t ( J ) J ^ ; 7 ^ <r=-l,7=-l<?=-l K # f j 7 ( A u 2 F Z " ) ( B ^ F Z j ) d e t ( J ) ^ ^ (3.97) 'l=l '2=1 '3=1 f = - 1 l/=-l f = - 1 > v ' * v ' g(i.i,o no Therefore the element stiffness matrix is computed as follows: K , =iiiTTf(K, ) de t ( J ) (FZ- F Z J ) < W < T ;,=1 / 2 = i / 3 = i ^=_, , ^=_, > v ,< v , g(.tn,o no (3.98) 3 2 2 "cp ";p /,=1 t2=\ / 3 =1 ATG/>=1 K1P=\ The definitions of the notations are the same as ones for Equation (3.91). 3.6 Time integration After developing the F E M solution in space, the resulting ODE's in time should be integrated to obtain the transient solution. The 0 method is used to evaluate the non-analytical integrals in time [61], [62]. The dominant equation is the rather stiff cure kinetics O D E that couples with the heat transfer equation as it evolves much more rapidly compared with the heat transfer equation. On the other hand, the diffusion equation incorporates both spatial and time parameters, which necessitates huge and expensive matrix inversions at each time step. To avoid this challenge, one can think of using different time steps for these ODE's . As a matter of fact, F E M solutions with multiple time steps have long been introduced in the literature. Brackbill and Cohen have reviewed articles on stiff PDE ' s with multiple time scales in different branches of 116 engineering and science [63]. Belytschko et al. have reviewed multi-step and partitioning methods in their review of time integration methods [62]. Smolinski et al. have investigated the subcycling methods in structural dynamics [64],[65]. In his more recent papers, Smolinski has developed multi-time scale solutions and stability conditions for element free methods and first order F E M [66],[67]. Sportisse investigates the error involved in splitting an O D E operator and use of different time-steps for stiff and non-stiff sub-operators [68]. A newer paper by Knoll et al discusses the accuracy of the solutions based on splitting the time integral [69]. More recent papers by Lowrie [70] and Rodriguez-Gomez [71] discuss the accuracy and stability of multi-rate methods. In this section a 0 -method based solution is described for the diffusion equation. In the next step, a multiple-rate solution is introduced for coupling the exothermic heat generation with the diffusion equation, which takes advantage of the slow development of diffusion phenomenon in avoiding cumbersome matrix inversions for each exotherm-time step used in evaluating the degree of cure c and hence the exothermic heat generation H. 3.6.1 Time Integration of Heat Transfer Equation Our linearized equation in space is in the following format ((3.82)): (3.99) We seek integration of this equation in time, i.e. Y , + A , - Y , + J((K + K ) T ) A . - ( ( T l + 4 l - T l ) + ( T l + A / - Y , ) ) = 0 . (3.100) where 117 (3.101) Y , + A , - Y , = JHdt. l+Al H is extremely nonlinear in our case and therefore it is challenging to evaluate JHdt i efficiently and accurately. In the case of thermoset polymer matrices during exothermic heat generation, this integral can be evaluated explicitly in terms of c (Equation (3.22)and Equation (3.85)) as Y , = \Nj(vRpRHRc)tet{S)dTl . (3.102) n H , on the other hand, is calculated from qQ which is part of the user input. Based on equations (3.85) and (3.101) we have <P,= \QoNjdS BIN (3.103) Qo = \ a o d t • The 6 method is used to evaluate the last integral. 0 method is based on the following approximation: ]f(t)dt = {t2-tx){9f(t2) + (\-e)f(tx)) , (3.104) where 6 is an appropriate number in the range [0,1]. Evaluating the integrals in Equation (3.100) by 0 method leads to 118 J-PM - Y,) + [<>(K + K)| T ( + A ( + (1 - 9)(K + K)[ T,) A' (3.105) 3.6.2 Time Integration of Cure Kinetics Equation Equation (3.105) is solved for T ( + A ( , which is the vector of temperatures in the degrees of freedom at the end of the time increment. In order to evaluate (*F,+A, - *P,) we need to calculate (c, + A, - c , ) at each integration point ((3.102)) as (^ + A,-^)=iN ( v ^ A ) ( c , + A , - c , ) d e t ( J ) c/n . (3.106) n As c develops quite rapidly, smaller time increments are needed to track it compared to the time increments used in the diffusion analysis. Let these smaller time increments be denoted as At'p or exotherm-time increment, which is used in contrast with diffusion-time increments, At, in this section. "Time increment" in this document means Diffusion-time increment. If c'0p and are the degree of cure and temperature at the integration point under consideration at t'0p = t, then we evaluate c in the next exotherm - time increments using the explicit integration as follows: cLP=Cl+f(CXP-l)^'P, ™ = U,3 , . .M , (3.107) where the end of final exotherm -time increment must be identical with the end of the diffusion-time increment, i.e. t'p =t + At and m = \,...M are the subinterval notions of the time step At. In order to obtain T1P_X we interpolate between Tt'p and T/pAl to get 119 •IP (3.108) Tt'+Al should be extrapolated from T history for the first iteration. The following extrapolation method is developed by the author and is best suited for the problems of interest in this study: C = C , + [ v ^ J , / ( p C p ) e ] / K - i , C . ) A / 1 P , m = \,2,\.M . (3.109) For the second and higher iterations, the result from previous iterations is used as Equation (3:109) is based on neglecting the diffusion term within one diffusion-time increment in heat transfer Equation ((3.22)), i.e. This extrapolation approximation affects the maximum size of the diffusion-time increment as choosing large time increments may destabilize the solution no matter how small the exotherm-time increments are. In this work this limit is set by trial and error. 3.7 Nonlinear Procedures Equation (3.105) is a nonlinear equation in time as the material properties change with temperature and hence K , K and Y are temperature dependent. The cure kinetics function is also temperature dependent and so are the forcing functions *F and *F. (3.110) (pCP)C^7 = VRPRHR-Z (3.111) 120 Therefore, in order to solve this equation for T, + A , a nonlinear procedure should be adopted. In this section we review two schemes of iterative solution for nonlinear problems; one is Newton-Raphson and the other one is Successive Substitution Method. Both solutions are modified to be used for our problem. Specially, the standard format of successive substitution method is modified in a way that would resemble the Newton-Raphson solution so that it can be used in A B A Q U S as the User Element in this software is devised for Newton-Raphson type of procedure. 3.7.1 Newton-Raphson Iteration The well-known method to solve a set of nonlinear equations is the Newton-Raphson method. This method uses the Taylor series expansion to solve a nonlinear equation F(T, + A , ) = 0 as follows: (3.112) l+M,r+l T<+A(,r) F(rI/+A<,r) •> where T ( + A I , is the r'th approximation for T ( + A , . In our case we have F ( T ( + . ) = - ( Y ( + A ( - Y , ) + (^(K + K ) [ + A ( T , I+Al + ( 1 - 0 ) ( K + K)[T, ) (3.113) { dF Different components of the Jacobian matrix, are calculated as follows: 121 dY, dT. l+Al J i v dU dT NjdV = [NiPCP N.dV = C (3.114) The above is based on the assumption that the density is constant which is compatible with no volume change assumption or assuming the integrals to be evaluated over the initial volume. Also the contribution of the second term on the right hand side of Equation (3.113) to the Jacobian is: dT,, i n ( ( K + R ) L T - ) = ( K , + R . ) L + dk dT (dh_ k dT dV + dS (3.115) The summation £ in taken on all degrees of freedom. The second term is asymmetric dk. dk and is usually neglected as — is small. If — is to be used, it is equal to dT dT dk_dkdc_ dk ~dT~~dcdT 8T (3.116) In order to determine dc ~dT we should solve the following equation in time: d_(dc_) _ (dc\ \ d T j \ d T ) T ' dt (3.117) 122 where/is the cure kinetics function (Equation (3.22)). Equation (3.117) must be solved using the same strategy as is used for c. The contribution of exothermic heat generation to the Jacobian matrix is as follows: S T 1+Al V ^ A , J l(N,vRPRHR c | , + A , ) d V = i N t N j v R p R H R J £ ' \ l t + M / j V V\ U 1 dV I+Al J (3.118) The last term that is contributing to the Jacobian matrix is H , which is differentiated as follows: 5T,, 3(T ( + A , ) y \BIN J&/fys]= . (3.119) Calculation of the above integrals can become very cumbersome except that the calculation of relatively smaller terms is avoided by proper inspection. In case of Equation (3.118) a secant approximation to 8c_ 8T can make life much easier compared to integration of Equation (3.116). In order to get the approximation, the same procedure as described in the section on time integration of exotherm equation can be used to get a perturbation to c| based on a perturbation to T\i+A/. A secant approximation to dc_ dT would then be equal to: ll+Al dc_ dT M ~C\A) (3.120) U+Al 123 ATP\ is the amount of perturbation in temperature at the end of the increment, c l(+A( l(+Ar is the degree of cure calculated based on the perturbed temperature. For the case when the conductivity and thermal capacity do not change quickly with temperature, i.e. with time, we can simplify the Jacobian matrix as follows: S , = i c s + * ( K S + K s ) 4 j At, N,NJvRpRHR dV t+M J (3.121) Therefore, the force vector would be F = - C ( T „ s , - T , ) + ( # ( K + K ) T „ , + ( l - S ) ( K + K)T,) - i ( T « - , , ' ) - i ( ' - - * ' ) -(3.122) As a result, the equation for nonlinear iteration then becomes S ( T r + A r , , + l - T ( + A r , , ) - - F (3.123) 3.7.2 Successive Substitution Method The following procedure, which is based on the successive substitution method needs fewer calculations compared to the previous procedure. First, we rewrite Equation (3.105) as i c ( T ( t , - T , ) + ^ ( K + K ) L T , + , + ( l - ^ ) ( K + K | T , ) J_-At (3.124) where based on the 8 method we have: (3.125) 124 Now regroup terms in Equation (3.124) to get i ^ ( ^ n . J M i e - ( , - ' ) ( K + K ) ' ) T + (3.126, For small enough time increments, the changes in all the parameters are small and so successive substitutive method is convergent. So in iterative format, we have i C + » ( K + K U ) ^ - ( i C _ ( I - ' ) { K + R ) l ) ' l : + (3.127) where r indicates the r'th iteration. It should be noted that for all of the parameters except *P we can use the start of increment temperature for the first iteration. Therefore: T,+ A,,o=T, . (3.128) As A B A Q U S solves a transient problem with Newton-Raphson type of iterations, we should make a small rearrangement of terms in Equation (3.127) to make it look like a Newton iteration: ( ^ » C « L , ) ( U - ' ™ l = ( i c - | l - ( ) ( K ( K l ) „ p i 2 9 ) The extra multiplication in the above equation is actually performed locally within the elements and so it doesn't really add to the computational burden. For the case when conductivity and thermal capacity don't change quickly with temperature, i.e. with time, Equation (3.129) is simplified further to 125 ^ C ^ ( K + K ) ] ( W l - T , + A , „ ) = - F • (3.130) 3.8 Energy Norm In order to make sure that the procedure is convergent, an average norm of the temperature distribution is used as a metric. The best candidate for this measure is the energy norm that is defined by ([72]) IHf =>Mv'v) ' ( 3 - 1 3 1 ) where v is the distribution function and «(. , . ) is the bilinear form in the following type of weak form equation: a(w,v) = Z ( v ) fora l lv (3.132) for a linear functional L (.). Comparing with the weak form of the heat transfer equation (say Equation(3.77)) we note that a ( v ' v ) = Z j ( k W , ) WjdV i=l \ y \viVj.= v1 .K.v . (3.133) Therefore, the norm would be | | v f = V v r X . v . (3.134) In this study, however, the square root is not applied and therefore the Energy norm for a region is defined as elements in V 126 A n F E simulation is deemed to have converged if the change in energy norm as defined above does not exceed a tolerance level after further refinement in time step and/or spatial mesh. In our A B A Q U S subroutine, the energy norm of each element is given as an output and the user can sum up the values in the regions of interest and therefore check the change in energy norm due to changes in the number of terms, spatial mesh and time step and check for convergence. 3.9 Errors and Adaptive Solution Error analysis and its relation with adaptive solutions have been widely discussed in the literature. The classic papers on errors and adaptive solutions are the ones by Zienkiewicz et al ([73], [74], [75]). One can find a survey of the classic works in this area in the work by Oden ([76]). Another paper by Oden ([77]) describes the more recent developments in this area where the benefits of an hp method are investigated. Chapters 10 and 14 of the finite element book by A k i n ([44]) describe the common features of the adaptive methods in F E M and the error estimation methods for elliptic problems. More theoretical and abstract aspects of the problem can be found in the book by Morton ([78]) and also the book and papers by Eriksson et al ([72],[79],[80],[81],[82],[83],[84]). The shell element solution described in this section can accommodate an easy to formulate and implement error estimation and adaptive formulation. As the Fourier terms that are used in the description of the temperature distribution through the thickness of the element are hierarchical and orthogonal, the value of the last term at each corner of the shell element gives a good estimation of the error in that corner. This estimation of 127 the error can be used in determining the number of Fourier terms that should contribute to the global stiffness matrix and damping matrix. Our recommended strategy is to use an a posteriori error estimation to find the necessary number of terms at the end of one time step and use the result for the next time step. If the tolerance is chosen more restrictedly, then the estimation works fine for the next time step. The first step is to have an estimation of the error involved in the analysis. In a Fourier series, the error in estimating a function could be approximated by the last term(s) of the series. Here, we use the same assumption for our mixed formulation and the last two terms are a measure of through thickness error. Therefore we define I h ^ l ^ m ^ d ^ l ' l ^ - i l ) > (3.136) where ||eTeafe(;|| is an error norm for the vertex and |r„ r | and | ^ _ i | are the absolute values of the last two terms on the vertex. This error estimate gives us an estimate of how active the last terms are on the region of interest. The next necessary tools are an upper bound and a lower bound of error which controls adding and dropping of terms. Here we introduce two non-dimensional parameters LL and Lu as non-dimensional lower bound and upper bound for the error to avoid adding and dropping of terms, i.e. Z £ . A ^ Z ^ L < | ^ ^ | | < z , u . A ^ ^ , (3.137) ^ range ^ range 128 where At is the size of the time step and Trange and trange are the estimates of temperature range and time range of the whole analysis. The values of LL and Lv are evaluated by trial and error. If H e r ^ l exceeds the right hand side limit, then more terms are to be added to the edge, in our program and we add two terms. On the other hand, i f |^r e 4 f e e | | is less than the lower bound value on the left hand side of the equation, then two terms are dropped from the edge degrees of freedom. Adding and dropping of terms in our Finite Element code could mean renumbering of degrees of freedom of the system which could be troublesome in bookkeeping and tracking the values of each degree of freedom. Therefore, in our procedure, we set a maximum number of terms that could be reached throughout the analysis and the inactive terms in each increment are dealt with as zero Dirichlet boundaries. The Dirichlet boundaries are removed from the global Jacobian and force vector and hence the size of these matrices changes through time. Again, this could not be implemented in A B A Q U S as we do not have access to global matrices and cannot change the way the program performs assembly. 3.10 Evaluation of Material Properties at Integration Points In order to calculate the stiffness and force matrices we need material properties at the integration points of the element. On the other hand, for the case of the composite materials the inputs to the program are usually the layup material properties and material directions. 129 Part of the input to the program is the layup information in addition to the material properties and the principal material (orthotropic) directions for each layer. The layers are considered to be transversely isotropic as is the case for the materials of interest in this study, mostly unidirectional fiber composites. Three types of input are considered for the input file. In all these types the user should identify one of the surfaces of the element as the reference surface (the surface that includes nodes 1,2,3 and 4 in Figure 3-3). This wi l l establish the local coordinate system at each surface Gaussian point with its first two axes on the surface, all mapped from the corresponding parameters in the parent space (Figure 3-4 (a)). The third axis would be perpendicular to this surface. In the first type of input, the user simply defines the angle between the first local axis and the 'zero vector' of the layup in the reference surface in addition to the layup directions for all layers that is the angles they make with the zero vector. The 'zero vector' is the vector with which material orientations in the layup are defined. Second type of input covers the elements in which the layers are parallel to each other but not necessarily to the reference surface. In this case the user should provide the zero vector of the layup in local coordinates in addition to the layup directions for all layers in terms of angles they make with the zero vector. If the major directions (fiber directions) of layers are not parallel, then the user can input the major direction of the layers in the local coordinate system one after the other. This is the third method of input. 130 The above information shall be repeated for all Gaussian integration points of the element (The numbering of surface integration points is illustrated in Figure 3-4 (a)). If it is mentioned once, the program would consider it to be identical for all surface Gaussian integration points. A tributary length is considered for each one of the through thickness integration points KIP — 1 KIP equal to [DKGP >DKGP ] where DKGP is the thickness of the corresponding nlp n[P Gaussian integration point of the element, nIP is the total number of integration points and KGP and KIP are the indices for the Gaussian point and the integration point under consideration. The material properties of the tributary length are then smeared to get the matching properties of each integration point. 3.11 Program Algorithm In this section the flow of a program that uses the above-mentioned shell element is described. The program consists of a set of modules that are presented separately. The whole algorithm is implemented as a U E L Subroutine in A B A Q U S . The A B A Q U S implementation, however, does not incorporate adaptive procedure described in section 3.8 as the program core cannot be modified to add or drop degrees of freedom throughout the analysis and this is not possible in the current version(s) of A B A Q U S . The modules are introduced in alphabetical order. The main module is described in Section 3.11.14. 131 3.11.1 Adaptivity Module This module requires inspection of the values of trigonometric terms at the shell element edges in every increment to be used in Equation (3.136) for estimation of the error. This error should be passed onto other subroutines that have access to the boundary conditions of the system so that the boundary conditions are added or dropped from the system based on Equation (3.137). This was done in a M A T L A B code that was prepared for two dimensional thermochemical analysis. It could not be implemented in A B A Q U S mainly because the number of degrees of freedom for each boundary condition could not change through one A B A Q U S step, .i.e. it could not adaptively change from one increment (or iteration) to the next. Secondly, in A B A Q U S User Element, one has local access to the degrees of freedom and their values and there is no way of globally marking some degrees of freedom for adding or dropping by changing the number of boundary conditions. The M A T L A B subroutine for adaptivity works as follows: 1) L O O P over all surrounding vertices of the shell elements. 2) READ the active terms of the vertex from the Number of Terms Prediction matrix. Note: The Number of Terms Prediction Matrix is active in the M A T L A B code and contains the number of terms the programs deems necessary for each increment. This matrix is then used in marking the inactive terms as zero boundaries. 132 3) C A L C U L A T E the error norm for the edge ( e r r ^ ) using Equation (3.136). Note that i f only one term is active at the edge, the error norm would then be equal to \T{\. ii n T 4) IF \erredge\ is less than L, • At • (Equation (3.137)), then drop a ^ range couple of terms to the active terms for next increment and update The Number of Terms Prediction Matrix accordingly. T 5) ELSEIF \erredge\ is more than Lv-At-^^- (Equation (3.137)), then add ^ range a couple of terms to the active terms for next increment and update The Number of Terms Prediction Matrix accordingly. 6) END LOOP. The Number of Terms Prediction Matrix is used in marking the Dirichlet Boundary Condition Locations Vector, in which the Dirichlet boundaries are marked as zeros (versus ones for other degrees of freedom). This vector is used in removing the degrees of freedom related to the Dirichlet Boundary from the global stiffness matrix (These concepts could not be implemented in A B A Q U S ) . 3.11.2 Boundary Contribution to the Stiffness and the Force Vector Module (Top and Bottom Surfaces) This module evaluates elements of H and K for the top and bottom surfaces in Equation (3.85) using Gaussian quadrature in the surface; i.e. evaluating 133 "CP Hj = j aoNj smf{i) dQ= £ wGKGP • q0 • AT (%KGP, nKGP ) • smf (J (gKGP, rjKGP, ± l ) ) , o KGP=\ (3.138) and, Kv = lhNlNJsmf(J)dCl (3.139) = Z w G « 3 / > • * • N, {^cr, VKGP ) • Nj {ZKGP > VKGP )' smf ( J (ZKGP » » i 1 ) ) • KGP=\ This subroutine works as follows: 1) L O O P over all the degrees of freedom in the surface. 2) L O O P over all the Gaussian points in the surface (Figure 3-4 (a)). 3) C A L L the Jacobian Calculator Subroutine to compute the inverse of Jacobian of the geometrical shape functions on the top and bottom of the element (refer to Section 3.4.1) 4) C A L C U L A T E smf (J) from Equation (3.64) for top and bottom surfaces at the Gaussian point. 5) C A L L the Surface Shape Function Module. 6) ADD UP the contribution of each Gaussian point to Ktj and Hj (Equation (3.85)). 7) END LOOP. 8) END LOOP. 134 3.11.3 Boundary Contribution to the Stiffness and the Force Vector Module (Edge Surfaces) This module calculates the contribution of convective boundary condition of edge surfaces to stiffness matrix and the force vector. For a surface parallel to the £ axis (like the one in Figure 3-4 (b)), we have Hj = \q0NjSmf(3)da = Z WGKGP-{ • % • M m j {%KGP > ± 0 • j Zfnj d£ -Smfi^J {^KGP-i»±1. ° ) ) KGP-S=\ (3.140) and, Kv = jhNiNjSmfWdQ. KGP-£=\ Z ^ • * - ^ ( ^ , ± l ) - M . j ^ , ± l ) . jzfni-zfndt •smf(j(4KGF^±l,0) U=-> J (3.141) Note that the parameters {mnnl>),{mj,nJ^ are related to i,j according to the relation in the Equation (3.67). The edge surfaces parallel to t] are done in similar fashion. This subroutine works as follows: 1) L O O P over all four edge surfaces. 2) L O O P over all degrees of freedom in the surface. 3) L O O P over all the Gaussian points in the surface. 135 4) C A L L the Jacobian Calculator Subroutine to compute the inverse of Jacobian of the geometrical shape functions at the integration point. 5) C A L C U L A T E smf (J) from Equation (3.64) at the Gaussian point. 6) C A L C U L A T E the exact through thickness integration of shape functions at each Gaussian point. 7) ADD UP the share of each Gaussian point to the Jacobian and force vector (Equation (3.85)). 8) END LOOP. 9) END LOOP. 10) END LOOP. 3.11.4 Cure Module This module consists of two subroutines. The first one calculates the contribution of the exothermic heat generation to the force matrix: 1) LOOP over all the integration points. 2) C A L L Alpha-Change Subroutine (embedded in the same module) to obtain the change in degree of cure and also the perturbation in degree of cure caused by a small change in the end of time segment temperature. 3) END LOOP. 136 The other subroutine calculates the change in degree of cure at a point and the perturbation in this change based on a starting and end temperatures, initial degree of cure and the length of time step. This is based on the theory described in Section 3.6.2 and here is the pseudo-code: 1) READ maximum change of degree of cure and maximum time step for a cure sub-step (At' p, look at Section 3.6.2) and also the variation in the end temperature ( ATP\ in Equation (3.120), Section 3.7.1) to get the perturbation in the degree of cure (cF\ ~c\l+Al m Equation (3.120)) 2) L O O P over temporal sub-steps based on the data read in the above stage. 3) IF it is the first iteration of the increment, then approximate the end-of-sub-step temperature by ^ = C , + [ v ^ A / ( p C , ) e ] / ( C . . C ) A / / ' • (3-142) 4) E L S E interpolate between the subsequent time-steps as follows: T?=Tf+(T£-Tf).{m/M) . (3.143) 5) C A L L the Cure Kinetics Subroutine (embedded in the Database Module) to read the end-of-sub-step degree of cure. 6) END LOOP. The loop described above is repeated for the perturbed end temperature and the result is subtracted from the degree of cure and divided by the variation in temperature to get the perturbation rate (Equation (3.120)) 137 3.11.5 Database Module This module consists of three subroutines. The first one reads the material property based on a material ID, temperature and degree of cure. The second one has the data for cure kinetics and gives out the temporal derivative of degree of cure based on temperature, degree of cure and the cure kinetics model ID. The third subroutine contains the control parameters of the solution and sets the values for number of quadrature points through the thickness (n IP ), number of Gaussian integration points in surface direction ( nGP ), initial degree of cure, cure kinetics model ID, maximum change of degree of cure in a cure sub-step, maximum size of a sub-step, perturbation in temperature to calculate the change in degree of cure (Equation (3.120), Section 3.7.1), the value of 6, method of iteration and the method of sub-stepping. 3.11.6 Element Stiffness and Capacity Matrices Module This module reads in the geometry and material properties of an element and calculates the local stiffness K and capacity C matrices for finite element analysis (Equation (3.85) ). The number of terms and the integral matrix (look at Integral Matrix Module in Section 3.11.9 and its development and equations in Section 3.5) are also input into this module. This module contains four subroutines. Firstly, the Capacity Matrix Subroutine: 1) L O O P over all degrees of freedom in the element. 2) C A L L the Place-Finder Module to see where each degree of freedom sits in the element matrices. 3) C A L L the Surface Shape Function Module. 138 4) C A L L the Jacobian Calculator Subroutine to get the determinant of Jacobian of the geometrical shape functions at the integration points. 5) READ the necessary integral from the Integral Matrix and calculate the capacity matrix. (Equation (3.91)). 6) END LOOP. The stiffness matrix is now calculated with the following subroutine: 1) L O O P over all degrees of freedom in the element. 2) C A L L the Place-Finder Module to see where each degree of freedom sits in the element matrices. 3) C A L L the Surface Shape Function Module. 4) C A L L the Jacobian Calculator Module to get the determinant and inverse of Jacobian of the geometrical shape functions at the integration points. 5) READ the necessary integral from the Integral Matrix. 6) C A L L OSSIL3D Subroutine to elaborate more on the integral matrix and calculate the stiffness matrix. (Equation (3.98)). 7) END LOOP. The third subroutine in this module is the asymmetric part (second term in right hand side of Equation(3.115)). It uses similar procedure as in the symmetric part described above. The forth subroutine in this module is the OSSIL3D subroutine which performs the cumbersome calculations of Equation (3.97). 139 3.11.7 Element Temperature Distribution Module It contains the subroutines that map the temperature values at the degrees of freedom to the values in the integration points using the interpolation shape functions (Equation (3.67) in Section 3.4.2) 3.11.8 Heat Vector Module This module contains a subroutine that calculates H as in Equation (3.85), which is the contribution of heat generation to the force vector. Here is the pseudo-code: 1) L O O P over all degrees of freedom in the element. 2) C A L L the Place-Finder Module to see where each degree of freedom sits in the element matrices. 3) C A L L the Surface Shape Function Module. 4) C A L L the Jacobian Calculator Subroutine to get the determinant of Jacobian of the geometrical shape functions at the integration points. 5) READ the necessary integral from the Integral Matrix and calculate H (Equation (3.106)). 6) END LOOP. 3.11.9 Integral Matrix Module This module includes subroutines that calculate the integral matrices. The input to this module is the number of integration points and number of terms used in the analysis. It evaluates Ft's in the Equation (3.90) for the following cases and saves it an integral matrix: 140 CASE#\:f(C)-CASE#2: f(C) CASE#3:f(C) CASE #4 : / ( £ ) 1 < n, m < nT + 2 (3.144) zf„ is the through thickness shape function as in Equation (3.69). The integrals are saved in a integral matrix to be used throughout the analysis for different g (0 (see Equations (3.86)-(3.90)). 3.11.10 Jacobian Calculator Module The Jacobian Calculator subroutine in this module inputs the nodal coordinates of the spatial point of interest. The Jacobian of the geometry is defined in Equation (3.53) and is calculated using Equation (3.54). 1) C A L L the Surface Geometric Shape Function Subroutine to get the surface shape function and its directional derivatives. 2) C A L L the Thickness Geometric Shape Function Subroutine to get the geometric shape function through thickness and its derivative. 3) L O O P over the parameter kx=1,2,3. element and calculates the Jacobian ( J ) , its determinant ( | j | ) , and inverse ( J 1 ) at the 4) L O O P over the parameter ky=1,2,3. 5) L O O P over the parameter kz=l,2,3. 6) IF kx and ky are not simultaneously equal to 2, T H E N 141 7) C A L L the Geometric Place Finder Subroutine (same module) 8) C A L C U L A T E GF,4, GFiir GFI( (Equation (3.54)) using the in-plane and through thickness geometric shape functions and their derivatives (chain rule). 9) ADD UP the components of the summation from Equation (3.54). 10) END IF. 11) END LOOP. 12) END LOOP. 13) END LOOP. Another subroutine in this module is the Surface Geometric Shape Function that outputs the surface geometric shape functions and their directional derivatives at each node of the surface. The surface geometric shape functions are 142 JVG, = _I(l_^)(1_^)( l 7 + ^ + l) NG2 = +^(l + «7)(l-#)(»7 + #-l) NG3 = +^(l + *)0 + £)07 + £ - l ) NG4 = - I ( l - 7 ) ( l + fl(7-^ + l) NG5 =^0-^)(i+^)0-#) NG6 = I ( l + ? 7 ) ( l + £ ) ( l - £ ) NG7 =^0-^)0+^)0+^) NGS = I ( l - „ ) ( l + £ ) ( l - £ ) . (3.145) The numbering in the above equation is based on the lower surface in Figure 3-3. Another subroutine of this module is the Thickness Geometric Shape Function Subroutine that outputs the through thickness geometric shape functions and their derivatives at each through thickness node. The through thickness geometric shape functions are: \ (3-146) 7u 2 =-( i+<r) . Index 1 in Equation (3.146) refers to the bottom node and index 2 refers to the top node. The last subroutine in this module is the Geometric Place Finder Subroutine. This subroutine correlates the node numbering in the three directions to element node numbering presented in Figure 3-3. 143 3.11.11 Place Finder Module This module has a short subroutine that sets a universal way of arranging the element degrees of freedom in the local matrices and vectors. This is formulated based on Equation (3.67). This subroutine bridges between the way the degrees of freedom are arranged inside the user subroutines and the way they should be arranged to be sent to the A B A Q U S core. Its importance is that this setup can be changed with a quick change in this subroutine. 3.11.12 Smear Module This module includes subroutines that input the laminate properties and output the integration points' properties. This is done based on the argument in Section 3.10. The major subroutine in this module is 'Thickness Smear' Subroutine that smears the laminate properties in the tributary area of every integration point and assigns it to that integration point (look at Figure 3-5). 3.11.13 Superelement Module This module is the main module that deals with constructing the Jacobian matrix and the force vector. The main subroutine in this module is the Superelement Subroutine which retrieves all the input of the problem from the User Element Module and takes the following actions: 1) C A L L Element Temperature Distribution Module to get the temperature at the integration points at the beginning and end of a time step. 2) C A L L the subroutines in the Smear Module to get the material properties at the integration points for the current temperature and degree of cure. 144 3) C A L L the Cure Vector Module to get the change in the degree of cure, rate of change and its perturbation with a variation in the end temperature. 4) C A L L Element Capacity and Stiffness Module to get the capacity matrix (C) , stiffness matrix ( K ) and the asymmetric stiffness matrices (look at Section 3.7; Equation (3.114) and Equation (3.115)). 5) C A L L Boundary Jacobian and Force Module to calculate the contribution of the convective boundary conditions on these matrices. 6) C A L L Jacobian and Force Vector Subroutine (in the current module) to sum up different contributions to these matrices. Jacobian and Force Vector Subroutine is the second subroutine in this module which does the following: 1) R E A D the value of 6 (Section 3.6.1) and the method of nonlinear solution (Section 3.7) from Parameter Control Module. 2) A D D U P the contribution of exotherm and external heating to the force vector (Section 3.6.2). 3) A D D UP the contribution of previous approximation of temperature distribution to the force vector. 4) A D D U P the contribution of C and K matrices to the Jacobian stiffness matrix (Section 3.7). 5) IF the method of nonlinear solution is the Newton-Raphson method, then add the contribution of — ^ (Section 3.7). 145 6) If the non-symmetric contribution of the K matrix to the Jacobian is to be considered, add it as well (Equation (3.116)). 3.11.14 Surface Shape Function Module This Module includes a subroutine that calculates the surface shape function M ( . at any point based on Equation (3.68). The corner id (/) is also input into the subroutine. 3.11.15 User Element Module This module includes two characteristic A B A Q U S subroutines. One is the U E L Subroutine where the program puts in the material properties form the input file and the temperature and state variable vectors from the previous step or iteration and returns the Jacobian and force vector and an updated vector of state variables. The second subroutine initializes the state variables and is called SDVTNI. Here is the flowchart description of the U E L prepared for our shell element. 1) INITIALIZE; set the Jacobian matrix (S in Equation (3.123)) to unit matrix and the force vector (F in Equation (3.123)) to zero. 2) READ boundary temperature, convective heat transfer coefficient and the number of operating degrees of freedom. These are passed in from the input file. 3) IF it is the initializing step, C A L L the Integral Matrix Module and prepare the integral matrix. 4) READ the temperature at the start and end of the time step. 5) READ the layup information that is passed in from the input file. 6) READ the degree of cure that is saved in the state variable vector of the element. 146 7) C A L L the Superemelent Module and get the Jacobian and force vector and updated matrix of degrees of cure. 8) SAVE the update degrees of cure in the state variable vector. 9) C A L C U L A T E the energy of the element (\\T\\E from Section 3.8). 3.12 Summary and Future Work This chapter described a shell model that enhances the heat transfer finite element analysis by using the shape functions of analytical solution in the through thickness direction. A few of the useful shape functions are derived for different boundary conditions although only one of the sets of shape functions (the set derived for Dirichlet type of boundary condition) is practically tested in this thesis. One big challenge in using the exact sinusoidal shape functions is in spatial integration. This is overcome through an innovative efficient spatial integration method which reduces the integration time and decreases the number of integration points. As for the time integral, it was noted that a huge amount of computing time is "wasted" as a result of coupled chemical and thermal analysis because the second one of these analyses has a much larger time scale relative to the first one. The thermal problem on the other hand connects the spatial and temporal dimensions and so the coupled problem would become a big spatial set of equations that should be solved in small time steps. A n innovative solution is proposed in this chapter that neglects the diffusion part of the thermal problem within the small time steps needed in the chemical part of the solution. The diffusion is accounted for within the bigger time steps used for the main iterations. 147 Another feature in the solution proposed here is the straightforward adaptive solution that was easily combined with solution. The numerical solution described in this chapter is very promising in reducing the computational time of thermal analysis. However, there is still room for improvement. Here are two main areas to be considered in future research: • Use of appropriate shape functions based on the boundary conditions and layers. • Separating the through thickness analysis of the shell from the whole mesh using Global-Local Finite Elements concepts. The first concept reduces the number of degrees of freedom and helps with using thicker elements that incorporates different materials. The second concept makes the adaptive solution easier to implement and reduces the size of global stiffness matrix. It might require extra iterations and cause oscillation of the solution which is the main challenge to be addressed in a future research. 148 X Figure 3-1 - Schematic of the notions in heat transfer in a general solid as is presented in Section 3.2.1. The thinner line presents boundary condition of the first kind while the thicker line is the boundary condition of the second or third kind. Figure 3-2-A typical 4-noded heat transfer shell element. 149 x (a) Physical Element Figure 3-3 - General 3D solid-like shell element in the physical coordinate system (top) and its mapped configuration in the parent coordinate system (bottom). (b) Figure 3-4 - Schematic of Gaussian integration points in surface (a) and through the thickness (b) of a shell element. 150 D KGP Tributary area of £ x Figure 3-5 - Schematic of typical layup (right side; full lines), integration points (^, g2, £3 ) and through thickness discretization points (£{, Q2, £3, CJ defined in Section 3.5. This is drawn at integration point KGP with the thickness equal to DKGP. Also illustrated in this drawing is the tributary area of an integration point as is discussed in Section 3.10. 151 4 CASE STUDIES 4.1 Introduction In this chapter a few cases of thermochemical analyses of composite structures during processing in the autoclave are presented to demonstrate the applications of the methodologies that were developed in the previous two chapters. The simple diagrams developed in Chapter 2 and summarized in Section 2.10.1 were used in the cases presented in this chapter to give insight into the capabilities and limitations of the simple diagrams. A l l but one of the F E M modelling with higher order shell elements were performed using A B A Q U S nonlinear finite element analysis software. A B A Q U S linear element (Element ID in A B A Q U S : DCC3D8) for heat transfer analysis and the quadratic element (Element ID in A B A Q U S : DC3D20) for heat transfer analysis were used in the simulations of this chapter. A User Element Subroutine (UEL) was developed based on the concepts presented in Chapter 3 of this thesis for 3D thermochemical analysis of composites. The pseudo-codes for the modules included in this procedure have been presented in Section 3.11. A n A B A Q U S heat generation subroutine (FTEATVAL) was developed to generate heat at the integration points of classic built-in elements of A B A Q U S . This subroutine was based on the Cure Module of Section 3.11.4. A M A T L A B code that was developed by the author was used in demonstrating the adaptivity capabilities of the shell solution presented in Chapter 3. This program performs two dimensional thermochemical analysis of composites and since the assembly 152 and implementation of boundary conditions are performed in the code the adaptive procedure presented in Section 3.9 could be implemented in it. The latter could not be implemented in the current version of the A B A Q U S software as the user does not have access to the core of the program. Eight cases are studied in this chapter. Case I includes different single layer materials studied using the simple diagrams for different boundary conditions and heating rates. In Case II, a composite slab with asymmetric boundary conditions and under a temperature ramp is investigated. This problem is modelled using the simple diagrams and also using the shell element developed here and stacks of linear and quadratic elements. These solutions are compared with an exact solution based on the eigenvalue analysis. In Case III, a composite layer on top of a tool is modelled with asymmetric boundaries and undergoing temperature ramp. The solution to this problem is again carried out using linear and quadratic elements and also the current shell element. It is also solved using the simple diagrams and all solutions are verified against the exact eigenvalue solution. In Case IV, the simple diagrams are applied to the case of a single layer of composite with exotherm using an iterative solution. The error involved in the simple solution are investigated by comparison of the results with the finite element analysis results. In Case V , the adaptive solution presented in this thesis and implemented in a M A T L A B based code is applied to the case of heat transfer and generation in a composite part. Case VI presents a partial implementation of the adaptive solution in A B A Q U S which is based on having some insight into the solution. Case VII investigates the advantages of the subcycling technique developed in this thesis. Specifically, accuracy and time efficiency gained by using the subcycling solution is presented. In Case VIII, the solution is applied 153 to the case of a tool-part assembly similar to the practical procedures in the composites industry in order to demonstrate the practicality of the shell element developed here and also its compatibility with other A B A Q U S elements. 4.2 Case I: Some Case Studies Using Simple Techniques Here, we consider a simple case study where we study the thermal lag in tools of different thicknesses and materials. We consider aluminium, composite, invar, and steel tools, of total thicknesses of 6.25, 12.5, 25, 37.5, and 50 mm, heated in autoclaves with low, intermediate, and high heat transfer coefficient. For each case, we can calculate two temperature lags: first the lag at the start of a hold (assumed to be at 180 °C, having started at 20 °C), and second the steady state lag, which may not have been achieved by the time the hold temperature is reached. The results for aluminium, invar, composite, and steel tooling are shown in Table 4-1, Table 4-2, Table 4-3 and Table 4-4 respectively. In all the tables, the first value is the temperature lag at the start of a hold after a 160 °C temperature increase, and the second value in parentheses is the steady state lag for a sufficiently long ramp. If the values are the same, then there has been enough time during the ramp to the hold temperature of 180 °C to achieve steady state conditions. This is typically the case for the thinner tools and lower ramp rates. Inspection of these results reveals some interesting trends. The convective heat transfer coefficient (and thus the Biot number) is a critical parameter, and at low convective heat transfer coefficient value, the temperature lags are significant even in aluminium tooling of moderate thickness. In fact, on the underside of a tool, air flow plays a very important 154 role in determining the heat transfer coefficient, and this is generally poorly controlled, i f at all. Table 4-3 can be used to evaluate both composite tooling and composite part, since the inputs are the same for both. We observe again that at higher thicknesses, control of the heat transfer coefficient is very important in avoiding thermal lag. 4.3 Case II: A Composite Slab Undergoing a Temperature Ramp with Asymmetric Boundary Conditions A composite slab of thickness 2L -5 (cm) is subjected to air temperature ramp of r = 3 ( ° C / m i n ) for 30 minutes. The initial temperature is T0 = 3 0 ( ° C ) . The convective heat transfer coefficient is hj. -\00(w/m2K) for the top surface and hB = 50(w/m2K) for the bottom surface. In this example we find the maximum temperature lag in the composite slab at the end of the ramp. The material properties for the through thickness direction of the composite slab could be taken from Table 2-1. The following are the parameters needed in the calculation using the closed form solution presented in Chapter 2: az = kz/pCP =0.69/(1580x870) = 0.502 x l0~ 6 (/M2/sec) 1 = 0.05/2 = 0.025 (cm) r = 3.0/60 = 0.05 (°C/sec) BiB =hBL/kz =(50x0.025)/0.69 = 1.8116 BiT = f\L/kz = (100 x 0.025)/0.69 = 3.6232 Fo = azt/L2 = (0.502 x 10"6 x (30 x 60))/0.025 2 =1.44 5 8 R = -rL2/az = -(0.05 x 0.025 2)/0.502x 10"6 = -62.25 (°C) . 155 Exact solution A computer program was developed to solve the eigenvalue problem and compute the temperature distribution. The minimum temperature after 30 minutes was found to be 72.47(°C). The result is based on a 10-term solution with a relative error of T -U- = 4.0 x 10"5. The maximum lag at the same instance is Tlag = -47.53 ( ° C ) . Using simple diagrams Strategy #1 from Section 2.4.4-1 The Fourier number and R number are calculated based on the original half-thickness in this strategy. On the other hand, the equivalent Biot number is calculated from Equation (2.65) to yield Bit = 2.47. Based on this Bie and the parameters calculated above, one can read 9m = 0.76 from Figure 2-4 or its magnified version (Figure 4-1) leading to A T s-62.25x0.76 = -47.3 l ( °C) , and the minimum temperature is 7/ = 7/0 + r x f + A7/ = 30 + 0.05x(30x60)-47.31 = 72.69(°C) . Strategy #2 from Section 2.4.4-2. The location of the minimum temperature at steady state is calculated from Equation (2.59) to be 156 J 1_ Bi Bi = _ 0 0 9 7 6 p „ 1 1 2 + — + — 4 = L (l - ^ ) = 0.025 x (1 + 0.0977) = 0.0274 . Thus resulting in the following parameters: Foe = aztfl2 =(0.502x10^ x(30 x 60))/0.0274 2 =1.2036 Re = -rL2/az = - (0 .05x 0.0274 2)/o.502x 10"* =-74.7769 (°C) Bie =hLjkz = (100x0.0274)/ 0.69 = 3.971. Based on the above values and using Figure 2-4 or its magnified version (Figure 4-1), 6 = 0.64 and therefore we have m AT s -74.78 x 0.64 = -47 .86(°C) T = T0 + r x t + AT = 30 + 0.05 x (30 x 60) - 47.86 = 72.14(°C) . Numerical solutions Current Shell Element A shell member was modelled with the above thickness and material properties. The in-plane size of the element is 10cm x9cm . Table 4-6 shows a convergence study on the element. The converged minimum temperature is 72.43 (°C) which is calculated using 7 terms (total of 9 unknowns). The C P U time for this analysis was 45 (sec). Linear Elements A stack of linear elements (our shell model with analytical terms suppressed) was used to model the thickness of the stack. The same planar size of element as the shell model was used. Table 4-7 shows the convergence study. The minimum through thickness 157 temperature value using 8 elements (total of 9 unknowns) was 72.55 (°C) and the computational time was 49 (sec). Quadratic Elements Similarly, a stack of quadratic elements (user elements) with the same planar size as above was modelled. Table 4-8 shows the convergence study. The minimum through thickness temperature value for 4 elements (total of 9 unknowns) was 72.43 (°C) and the C P U time was 47 (sec). Summary Table 4-9 shows the comparison between different solutions and the error involved for the same number of unknowns in all the solutions (9 unknowns) using the same time step size (At = 3 (Sec)). The F E M solutions take similar amount of C P U time as is expected for the same level of accuracy. The reason that the solution using the current shell element takes slightly less C P U time is that it uses one integration point through the thickness while the stack of linear elements use a total of 16 (two integration points through the thickness of each element) and the stack of quadratic elements uses a total of 12 (three integration points through the thickness of each element). 4.4 Case III: A Two-Layer Material with Asymmetric Boundary Conditions Undergoing a Temperature Ramp In this case a hybrid model was investigated where a composite part with 2L2 = 1 (cm) sits on an invar tool (material properties in Table 2-1) with 2LX -2(cm). The autoclave temperature ramps up from the initial temperature of 20 (°C) to 200 (°C) in 90 minutes 158 which means we have a heating rate of r = 2(°C7min) . The convective heat transfer coefficient is hj = 100(W I m2K) for the top surface and hB - 50(W I m2K) for the bottom surface. Exact Solution A computer program was developed for the eigenvalue analysis of a two layer hybrid system in one dimension. A maximum of 10 eigenvalues were used and the minimum temperature in the composite was found to be 164.61(°C). The relative error involved is T •Ml R = 1.53x10 1 4 (the value of Re on the next page). Approximate solution using the diagrams The parameters involved in the simplified solution are calculated and summarized in Table 4-10. Using these parameters the equivalent convective heat transfer coefficient is calculated from Equation (2.105) to be hss =31.37 (W / m2K) . So the minimum temperature occurs in the first layer in the steady state phase (this is used as an approximation for the transient phase). Therefore, the equivalent Biot number on top of hssL the tool would be Bin = 1 = 0.0285 . Hence, the location of the adiabatic line (where K minimum temperature occurs) is J 1_ C = B \ B \ =0.22 . 2 + — + — Bij BiB Therefore, the equivalent one-layer-symmetric-boundary parameters are calculated to be 159 Bie = max (BiB, BiTl) x ( l + \ga |) = 0.055 ^ = V ( 1 + |C,|) 2=-1.84 (°C) Foe=^/(l + \{a\f=96.7. Using the above parameters we can use Figure 2-4 or its magnified version (Figure 4-2) to read Qm =18.4 and hence the minimum temperature in the tool is estimated to be Tm = R6m + Tx =166.05 ( °C) . The minimum temperature in the composite part occurs at the interface between the tool and the part and we use the distribution diagram (Figure T -T 2-5) to estimate the value of this temperature. The diagram reads — = 0.99 and hence Ts =166.42 (°C) . A look at time to steady state diagram (Figure 2-6) shows that more than 99% of the steady state is achieved and so the steady state solution is a good approximation to the solution. In fact, the steady state solution gives 0 =0.5 + 1/5/ =18.51 . m e This value is only 0.5% higher than the more accurate value of 18.4 considered above. Using this value, we find Tm = 165.94(°C) and Ts =166.28 (°C) . Numerical solutions The problem was first modelled as a shell of composite material stacked on top of a shell of invar (2 cm of invar underneath 1cm of composite). Each shell has a planar dimension of 9cmx\0cm. The edge boundaries are considered to be adiabatic. The convergence was monitored using the energy norm and a maximum tolerance of 0.1% was used as the 160 criteria for convergence. With a constant time step At = 5 sec, this goal was achieved using 5 terms for the part and 3 terms for the tool. The number of terms in each section was constant throughout the analysis. The interface temperature in the converged solution was Ts = 1 6 6 . 4 ( ° C ) . Using linear elements (shell elements with zero number of terms) through the thickness, the same accuracy in energy norm was achieved using 8 elements in the tool and 8 in the part. The interface temperature was calculated to be Ts = 166.4 (°C) at the end of the analysis. Summary A two-layer tool and part hybrid was modelled using the strategy developed in Chapter 2, current shell element and the linear element (shell element with zero terms) in A B A Q U S . A l l the solutions are in close agreement with each other and with the eigenvalue solution of the problem. The solution could also be approximated by the steady state solution as the steady state was 'almost' achieved. The C P U times in all A B A Q U S models had the same order. 4.5 Case IV: Single Layer Behaviour during a Hold -Exotherm A typical autoclave process cycle is constructed from ramp and hold sections (Figure 2-1). However, during a typical composite material curing, the part starts to generate heat as a result of chemical reactions. This phenomenon is called exotherm and it is important to be able to model the exothermic effects using simple diagrams. 161 As a matter of fact, when there is an exotherm, we can evaluate Equation (2.15), and proceed in a similar fashion to the ramp case. Using cure kinetics data for composite (Table 4-5) and the smeared material properties of composite (Table 2-1), we can compare the predictions from the simple approach outlined here with the full 1-D transient solution attained using a transient finite element analysis. We first draw some useful diagrams for isothermal heat generation of this type of composite material. Figure 4-3 shows the rate of change of degree of cure as a function of degree of cure itself during isothermal curing. Figure 4-4, on the other hand, illustrates the heat generated per unit volume as a function of time for isothermal curing. A number of numerical runs were performed for hold temperatures of 170, 175, 180, 185, and 190 °C, for composite parts of total thicknesses of 1, 3, 5, and 7 mm. To obtain the exotherm for these cases the maximum amount of heat generated in the material from Figure 4-4 was used in combination with Figure 2-4. The maximum exotherm from these runs is compared against the predictions from the approximate approach in Figure 4-5. Exotherms range from essentially zero to about 14 °C. In processing of many typical aerospace composite parts, the approximate method, without any modification, does a very good job of predicting small exotherms, but starts to under-predict the higher values of exotherms. The cause of this is that the heat evolution rate is due to the cure kinetics, which is highly temperature sensitive. As the temperature rises due to the exotherm, there is a knock-on effect, and the reaction rate increases even more. This effect can be estimated by taking the first estimate of exotherm in the approximate method and increasing it by a factor corresponding to the temperature increase in the thermally activated term in the cure kinetics equation. This 162 gives us the second prediction line labelled "updated local temperature" in Figure 4-5, which is in very good agreement with the numerical predictions over the range up to 14 ( ° Q . 4.6 Case V: Adaptive Solution using the MATLAB Code In order to examine the adaptive solution presented in Section 3.9, a simple case was studied. The 6.0 (cm) thick model, top and bottom surface temperatures and the autoclave temperature associated with these temperatures are illustrated in Figure 4-6. The smeared material properties for the composite solid are taken from Table 2-1 and the cure kinetics equation and properties are presented in Table 4-5. The surface temperature histories are calculated by a standard finite element analysis of the problem using C O M P R O (a finite element analysis tool for process modelling of composites developed in-house by the U B C Composites Group) and are applied as Dirichlet boundaries to the composite solid. The number of elements for the F E M solution with linear elements was changed manually to arrive at a converged solution for the midplane temperature. This is presented in Figure 4-7 along with the solution from the current shell technique with varying number of terms, based on Equation (3.137). Figure 4-8 illustrates the change in the number of terms and time step size throughout the analysis. A comparison with Figure 4-7 shows that the number of terms reaches a maximum during the exotherm. The size of the time step decreases dramatically during the same period. 163 4.7 Case VI: Implementation of a Solution with Changing Number of Terms in ABAQUS The adaptive solution presented in Section 4.5 could not be implemented in A B A Q U S as this software currently does not allow the user to change the number of degrees of freedom or apply extra boundary conditions during a single step of analysis. However, an experienced user can estimate the necessary number of terms through each stage of the analysis and set up the input file in multiple steps so that the number of active terms is different in each step. As an example to illustrate this technique a case study is considered here in which the "number of terms" is varied through the analysis. This is achieved through setting up seven steps for the analysis and setting zero boundary condition for the terms that are to be suppressed in in each step, thus changing the number of active terms through the analysis. The tool in this example is a 5.0 cm thick invar material (properties in Table 2-1), on which sits a 5.0 cm thick composite with the smeared material properties from Table 2-1. The cure kinetics equation and properties are given in Table 4-5. The boundary conditions are hlop =\00(w / m2K} on the top and hbot =5®(W lm2K} on the bottom surfaces. Figure 4-9 illustrates the autoclave temperature, midplane temperature and degree of cure in a converged solution using shell elements. Figure 4-10 shows different patterns for change in the active number of terms in the part studied here. These patterns have similar shapes and all increase the number of terms in temperature ramps and during the 164 exothermic heat generation. The number of terms in the tool is kept constant and is equal to 3. Figure 4-11, Figure 4-12, and Figure 4-13 show the effect of number of terms in the part on the accuracy of results in terms of energy norm, midplane temperature and final distribution of degree of cure through the thickness of the part. 4.8 Case VII: Subcycling Solution Figure 4-14 shows a simple model that was analyzed to illustrate the efficiency of the subcycling process introduced in Section 3.6.2. This is basically a two-element model of the composite material made with linear elements (shell element with zero terms) with 10 integration points through the thickness of each element. The subcycling solution uses time steps that are small enough to keep the variation in the degree of cure less than 0.001 at each exotherm time increment (refer to Section 3.6.2, Equation (3.107)). The smeared material properties for the composite are given in Table 2-1 and its cure kinetics equation and properties are listed in Table 4-5. The temperature in the autoclave and also the midplane temperature (based on the converged results for the two element model) are illustrated in Figure 4-15. This converged solution is Run No-1 in Table 4-11 (first row), in which both the subcycling and classical time stepping solutions have converged. Columns 2, 3 and 4 of the table show the analysis specification in A B A Q U S which dictate the time stepping for F E M model. The time increment size varies between the minimum and maximum time increment in the table. The program also checks for the maximum change in the temperature in the element to limit the size of the time step. If the maximum temperature change is more than the value specified in the 4 column, the program reduces the time 165 increment and recalculates the temperature at the end of the new time increment. If the temperature change is less than this value, the time increment size increases. Run-time and the accuracy of the final degree of cure are presented both for the subcycling solution and for the classical time-stepping solution. Table 4-11 shows that the amount of time spent for the whole analysis is governed by these limits both for the subcycling method and classical time-stepping procedure. The subcycling scheme, however, yields more accurate solutions. As a matter of fact, i f we can tolerate about 3% error in the analysis, then the subcycling method needs a fraction of one second (rows 7 to 10) for the whole analysis, while the classical solution does not converge. Figure 4-16 to Figure 4-19 illustrate the better accuracy of the solutions gained through subcycling scheme compared with the ones from the classical solution both for the midplane temperature and for the energy norm. 4.9 Case VIII: Application of Higher Order Shell Elements in Modelling a Curved Composite Part on a Tool In this case study a model of a 1 cm thick composite part on a tool is simulated both using the A B A Q U S standard elements and the shell element developed in this work. The tool material is invar (see Table 2-1 for material properties of the tool and the part). The cure kinetics equation and its variables are given in Table 4-5. Figure 4-20 illustrates a 2D schematic of the problem investigated in this section and the dimensions of the tool and the part. Our three dimensional model is an extrusion of this 166 sketch by 1.0 cm. The boundary condition is h - \ 00(w/m2K) on top and right surface of the model in Figure 4-20 and h = 50 [w I m2K) on the bottom and left surfaces of the model. The elements used in the A B A Q U S analysis are quadratic elements for the tool (Element ID in A B A Q U S : DC3D20) and linear elements (Element ID in A B A Q U S : DCC3D8) for the part (Figure 4-21). These elements are tied (constrained) at the interface between the tool and the part and therefore it is not necessary for them to have identical nodes at the interface. In order to simulate the heat generation inside the composite part, the Heat Generation Module of A B A Q U S ( H E A T V A L ) was linked to the program to simulate the curing and the resultant exothermic heat generation. In a separate analysis, the built in quadratic elements in A B A Q U S were used to model the tool again and the user defined shell element developed here was used to model the part (Figure 4-22). The interface nodes are constrained and therefore there is no need for identical nodes on both sides of the interface. One big challenge in practical use of the shell element is in defining material orientation for the user element as it is defined with respect to the bottom surface of the element (the plane made of nodes 1 -4 in Figure 3-3). This is also the reference surface for defining the convective boundary conditions. The Fourier terms are defined in a direction perpendicular to this surface. This is actually the surface that should be put in contact with the tool. With the meshes generated using the current A B A Q U S pre-processor (CAE), however, it is not possible to ensure that this surface is the one in contact with the tool and it is labour intensive to rotate the elements in such a way that they all have this surface in contact with the tool. 167 Figure 4-23 shows the applied autoclave temperature and also the degree of cure at the mid-surface of the composite part for section C (Figure 4-20). These numbers are identical for both analyses. Figure 4-24 shows a snapshot of the temperature distribution in the tool at t=173 min when there is heat generation in the material because of cure as is seen in Figure 4-23 as well. Figure 4-25 and Figure 4-26 illustrate the time history of the temperature at top, bottom and mid points for the cross sections C and A of the composite part (These sections are defined in Figure 4-20). It can be seen that the solution for the current shell element and the A B A Q U S classical elements are identical. Figure 4-27 and Figure 4-28 illustrate the temperature distribution in sections A , B and C as defined in Figure 4-20. These are taken at the end of the first ramp (t=45 min) and during the heat generation in the composite (t=l 73 min). The minimum temperatures in these cases could not be reproduced from the simple diagrams in Chapter 2 as there is a significant edge effect involved here. In fact as the adiabatic line at 45 min occurs in the tool, we can approximate the Biot number of the system by ^ 50x0.031 e K 11 As is seen in Figure 2-31, for £ e d < 1.0, the error is more than 5%. This edge effect in the tool also affects the temperature during exotherm. A total of 15 linear elements were used in A B A Q U S to achieve solution convergence. In the case of the User Element, 5 terms were active throughout the analysis. The analysis time was 710 seconds with A B A Q U S standard element and 684 seconds with the shell 168 element. While these times have the same order, it should be noted that the user element does not need meshing and remeshing through the thickness of the composite shell. 4.10 Summary and Future Work Several cases of thermochemical analysis of composite structures during processing in the autoclave were investigated in this chapter to show the practical efficiencies gained in using the simple diagrams developed in Chapter 2 and also the shell solution developed in Chapter 3 of this thesis. Simple diagrams help with quick estimation of the solution and the shell solution takes away the necessity of meshing and remeshing through the thickness. It is also shown that the shell element can be used along with other classical elements in the same analysis. The future developments that can make our numerical solution more practical is a pre-processor that prepares the meshes of our shell element considering the element orientation of interest. Also, the simulation of more practical and industrial cases should be investigated with the current shell element to ascertain its efficiency and practicality. 169 Table 4-1 - Temperature lag in aluminium tooling for different heat transfer coefficients, autoclave ramp rates, and tool thicknesses - CASE I. f autoclave (°C/min) A7/(°C) 2L = 6.25 mm A 7 (°C) 2L = 12.5 mm A r (°C) 2L = 25 mm A r (°C) 2L = 37.5 mm A r (°C) 21 = 50.0 mm h=20 (W/m 2 K) 0.5 3.2(3.2) 6.3(6.3) 12.7(12.7) 19(19) 25.3(25.3) 1 6.3(6.3) 12.6(12.6) 25.3(25.3) 37.4(38) 48.5(50.6) 2 12.6(12.6) 25.2(25.3) 48.5(50.6) 66.7(75.9) 80.4(101.3) 3 19(19) 37.4(37.9) 66.7(75.9) 86(113.9) 99(151.9) A=100 (W/m 2 K) 0.5 0.6(0.6) 1.3(1.3) 2.5(2.5) 3.8(3.8) 5.1(5.1) 1 1.3(1.3) 2.5(2.5) 5.1(5.1) 7.6(7.6) 10.2(10.2) 2 2.5(2.5) 5.1(5.1) 10.2(10.2) 15.3(15.3) 20.4(20.4) 3 3.8(3.8) 7.6(7.6) 15.2(15.2) 22.9(22.9) 30.4(30.6) h=200 (W/m 2 K) 0.5 0.3(0.3) 0.6(0.6) 1.3(1.3) 1.9(1.9) 2.6(2.6) 1 0.6(0.6) 1.3(1.3) 2.5(2.5) 3.8(3.8) 5.1(5.1) 2 1.3(1.3) 2.5(2.5) 5.1(5.1) 7.7(7.7) 10.3(10.3) 3 1.9(1.9) 3.8(3.8) 7.6(7.6) 11.5(11.5) 15.4(15.4) The first value in the table is temperature lag at the start of a hold after a 160 °C temperature increase and the second value in the parenthesis is the steady state lag. 170 Table 4-2 - Temperature lag in invar tooling for different heat transfer coefficients, autoclave ramp rates, and tool thicknesses - CASE I. f autoclave (°C/min) AT (°C) 21 = 6.25 mm AT (°C) 2L = 12.5 mm AT (°C) 21 = 25 mm AT (°C) 21 = 37.5 mm AT (°C) 21 = 50.0 mm h=20 (W/m 2 K) 0.5 5.4(5.4) 10.8(10.8) 21.7(21.7) 32.5(32.7) 42.8(43.9) 1 10.8(10.8) 21.6(21.6) 42.3(43.4) 59.9(65.5) 73.8(87.8) 2 21.5(21.5) 42.1(43.2) 73.2(86.8) 92.6(130.9) 105.5(175.6) 3 32.1(32.3) 59.3(64.7) 92.3(130.2) 109.8(196.4) 120.5(263.3) A=100 (W/m 2 K) 0.5 1.1(1.1) 2.2(2.2) 4.5(4.5) 7(7) 9.6(9.6) 1 2.2(2.2) 4.4(4.4) 9.1(9.1) 14(14) 19.1(19.1) 2 4.4(4.4) 8.8(8.8) 18.1(18.1) 27.9(27.9) 37.7(38.2) 3 6.5(6.5) 13.2(13.2) 27.1(27.2) 41.1(41.9) 54.1(57.4) h=200 (W/m 2 K) 0.5 0.6(0.6) 1.1(1.1) 2.4(2.4) 3.8(3.8) 5.3(5.3) 1 1.1(1.1) 2.3(2.3) 4.8(4.8) 7.5(7.5) 10.5(10.5) 2 2.2(2.2) 4.5(4.5) 9.6(9.6) 15.1(15.1) 21.1(21.1) 3 3.3(3.3) 6.8(6.8) 14.3(14.3) 22.6(22.6) 31.5(31.6) The first value in the table is temperature lag at the start of a hold after a 160 °C temperature increase and the second value in the parenthesis is the steady state lag. 171 Table 4-3 - Temperature lag in composite tooling for different heat transfer coefficients, autoclave ramp rates, and tool thicknesses - CASE I. f autoclave (°C/min) A r (°c) 21 = 6.25 mm A r (°C) 21 = 12.5 mm A r (°C) 2L = 25 mm A r (°C) 2L = 37.5 mm A r (°C) 21 = 50.0 mm h=20 (W/m 2 K) 0.5 1.9(1.9) 3.9(3.9) 8.5(8.5) 13.7(13.7) 19.5(19.5) 1 3.7(3.7) 7.8(7.8) 16.9(16.9) 27.3(27.3) 38.6(39) 2 7.5(7.5) 15.6(15.6) 33.6(33.8) 52.3(54.7) 69.6(78.1) 3 11.2(11.2) 23.4(23.4) 48.9(50.8) 71.8(82) 90.5(117.1) h=lOO (W/m 2 K) 0.5 0.4(0.4) 1(1) 2.7(2.7) 5.1(5.1) 8.1(8.1) 1 0.9(0.9) 2.1(2.1) 5.5(5.5) 10.1(10.1) 16.1(16.1) 2 1.8(1.8) 4.2(4.2) 10.9(10.9) 20.3(20.3) 32.1(32.2) 3 2.6(2.6) 6.2(6.2) 16.4(16.4) 30.3(30.4) 47.4(48.3) £=200 (W/m 2 K) 0.5 0.3(0.3) 0.7(0.7) 2(2) 4(4) 6.6(6.6) 1 0.5(0.5) 1.4(1.4) 4(4) 8(8) 13.2(13.2) 2 KD 2.7(2.7) 8.1(8.1) 16(16) 26.5(26.5) 3 1.6(1.6) 4.1(4.1) 12.1(12.1) 23.9(24) 39.4(39.7) The first value in the table is temperature lag at the start of a hold after a 160 °C temperature increase and the second value in the parenthesis is the steady state lag. Ill Table 4-4 - Temperature lag in steel tooling for different heat transfer coefficients, autoclave ramp rates, and tool thicknesses - CASE I. f autoclave (°C/min) A r (°c) 2Z = 6.25 mm A r (°C) 2L = 12.5 mm A r (°C) 2L = 25 mm A r (°C) 2L = 37.5 mm A r (°C) 2L = 50.0 mm H=20 (W/m 2 K) 0.5 4.8(4.8) 9.5(9.5) 19.1(19.1) 28.6(28.7) 37.7(38.3) 1 9.5(9.5) 19.1(19.1) 37.6(38.2) 53.8(57.3) 67.1(76.5) 2 19(19) 37.5(38.1) 67(76.3) 86.3(114.6) 99.3(153) 3 28.5(28.6) 53.7(57.2) 86.2(114.5) 104.2(171.9) 115.3(229.5) #=100 (W/m 2 K) 0.5 1(1) 1.9(1.9) 3.9(3.9) 5.8(5.8) 7.8(7.8) 1 1.9(1.9) 3.8(3.8) 7.7(7.7) 11.6(11.6) 15.6(15.6) 2 3.8(3.8) 7.7(7.7) 15.4(15.4) 23.2(23.3) 31(31.2) 3 5.7(5.7) 11.5(11.5) 23.1(23.1) 34.5(34.9) 45.3(46.8) h=200 (W/m 2 K) 0.5 0.5(0.5) 1(1) 1.9(1.9) 3(3) 4(4) 1 K D 1.9(1.9) 3.9(3.9) 5.9(5.9) 8(8) 2 1.9(1.9) 3.9(3.9) 7.8(7.8) 11.8(11.8) 16(16) 3 2.9(2.9) 5.8(5.8) 11.7(11.7) 17.7(17.8) 23.9(23.9) The first value in the table is temperature lag at the start of a hold after a 160 °C temperature increase and the second value in the parenthesis is the steady state lag. 173 Table 4-5 - Composite cure kinetics equation and its variables. Cure Kinetics Equation dc Kcm(\-c)" dt \ + eC{C-(aO+aCTT)} K,=Aie-^,RT Variable Description Values Units c Resin degree of cure. T Resin Temperature °K vf=\-vR Fibre volume fraction 0.573 Total resin heat of reaction (c= 0 to 1) 5.40E5 J/kg A Pre-exponential factor. 1.528E5 /s Activation energy. 6.65E4 J/mol m Equation superscript. 0.8129 _ n Equation superscript. 2.736 _ R Gas Constant 8.3143 J/(mol K) C Diffusion constant 43.09 aco Diffusion constant -1.684 -Diffusion constant 5.475E-3 /°K 174 Table 4-6 - Convergence study for the current shell element - CASE II. Test Number Number of Terms Time Step (sec) Energy Norm (W.C) 1 1 20 453.37 2 3 10 462.12 3 5 5 464.24 4 7 3 464.98 Table 4-7- Convergence study for the linear element - CASE II. Test Number Number of Elements Time Step (sec) Energy Norm (W.C) 1 1 20 372.77 2 3 10 442.58 3 5 5 455.49 4 7 3 460.03 5 20 2 464.59 175 Table 4-8 - Convergence study for the quadratic element - CASE II. Test Number Number of Elements Time Step (sec) Energy Norm (W.C) 1 1 20 461.08 2 2 10 463.55 3 3 5 464.67 4 4 3 465.16 Table 4-9 - Error study for the minimum temperature in the plate from different methods of solution - CASE II. Minimum temperature (C) Error (%) Exact Solution 72.47 0.00 Strategy #1 72.69 0.30 Strategy #2 72.14 0.46 Superelement 72.43 0.06 Linear Elements 72.55 0.11 Quadratic Elements 72.43 0.06 176 Table 4-10 - Parameters used as input for CASE III. Layer 1 Layer II r (C/sec) 0.033 0.033 a z (m 2/sec) 2.67E-06 5.02E-07 L(m) 0.01 0.005 R(C) -1.24E+00 -1.64E+00 BI 0.045455 0.72463768 S (W/m zsec) 1100 138 Table 4-11 - Comparison of the efficiency of the subcycling method with the classical method of identical time-steping parameters - CASE VII. Analysis Specifications Subcyclin g Classical Run No Atmin (sec) Atmax (sec) ATmax (°C) Run Time (sec) Error in final DoC Run Time (sec) Error in final DoC 1 5 5 10 635 0 707 0 2 10 10 10 238 0 248 0.20% 3 20 20 10 84 0 74 0.30% 4 10 100 20 6 0.20% 5.8 1.00% 5 10 200 20 2.4 0.30% 2.7 1.00% 6 10 1000 20 2 0.30% 1.5 2.00% 7 10 1000 100 0.9 2.50% DIVERGENT 8 50 1000 100 0.7 3.00% DIVERGENT 9 100 1000 100 0.6 3.00% DIVERGENT 10 200 1000 100 0.5 2.70% DIVERGENT Atmin: minimum time step used in ABAQUS analysis. AW: maximum time step used in ABAQUS analysis. ATmax: maximum allowable change in temperature. DoC : Degree of Cure 177 Figure 4-1 - Magnified section of Figure 2-4 used to read 9m - Case II. Fo Figure 4-2 - Magnified section of Figure 2-4 used to read Gm - Case III. 178 0.0025 N N 0 0.2 0.4 0.6 0.8 1 Degree of Cure Figure 4-3 - Rate of change of degree of cure as a function of degree of cure for different temperatures during isothermal curing of the composite — CASE IV. 0 600 1200 1800 2400 3000 3600 time (sec) Figure 4-4 - Heat generated per unit volume as a function of time for different temperatures during isothermal curing of the composite - CASE IV. 179 16 0 2 4 ?6J 8 10 12 14 15 E x o t h e r m p r e d i c t i o n f r o m n u m e r i c a l a n a l y s i s ( C ) Figure 4-5 - Error in exotherm prediction for single layer material using simple charts CASE IV. 250 t —200 150 A 3 ra 2.100 E 0) / T / 1tOD Tbot Top Surface temperature Autoclave air temperature profile Bottom'surface temperature 50 100 150 time (min) 200 250 300 Figure 4-6 — Boundary conditions for testing the efficiency of the adaptive solution CASE V. 180 300 g250 £200 3 -*-» SJ150 a §100 50 -\ 0 0 F E M with 10 elements Autoclave air temperature profile\^ Series method with maximum of 11 term^ 50 100 150 time (min) 200 250 300 Figure 4-7 — Comparison of the convergence of the series method with the FEM solution - CASE V. 12 10 (0 E 8 +•» o z 4 2 0 50 100 150 time (min) 200 250 _ I 1 I \ No of Terms J ^ 7 6 5 c" E 4 a v 3 « 2 - i 300 Figure 4-8 - Time step sizes and number of terms used in the series solution - CASE V. 181 300 250 id 200 A i 3 o 150 100 0 3000 6000 9000 12000 15000 18000 21000 time (sec) Figure 4-9 - Autoclave temperature, midplane temperature and degree of cure for Case VI. These are used as a guideline for choosing the patterns specified in Figure 4-10. 0 3000 6000 9000 12000 15000 18000 21000 time (sec) Figure 4-10 - Different patterns of change for number of terms in case VI. These are mainly based on inspecting the solution for autoclave and midplane temperature and change in degree of cure in Figure 4-9. 182 4000 3500 8 3000 5 2500 E2000 o £ 1 5 0 0 E> I 1000 500 0 0 3000 6000 9000 12000 15000 18000 21000 time (sec) Figure 4-11 - Energy norm throughout the analysis for different patterns of number of terms as presented in Figure 4-10 - Case VI. Figure 4-12 - Mid-plane temperature in time for different patterns of number of terms as presented in Figure 4-10 - Case VI. 183 0.94 0.78 -I 1 1 1 1 1 0 1 2 3 4 5 Distance from the interface into the Part (cm) Figure 4-13 - Final distribution of degree of cure through the thickness of the part for different patterns of number of terms as presented in Figure 4-10 - Case VI. h t =100 (W/m 2K) 2.5 (cm) 2.5 (cm) h b o t=50 (W/m 2K) Figure 4-14 - Two element model of a composite layer modeled in Section 4.8 that is used to verify the subcycling solution - CASE VII. 184 Figure 4-15 - Autoclave temperature and mid-surface temperature for the case studied in Section 4.8 - Case VII. Figure 4-16 -A comparison of the errors in the midplane temperature for RUN 4 of Table 4-11 - Case VII. 185 0 2000 4000 6000 8000 10000 12000 14000 time (sec) Figure 4-17 -A comparison of the errors in the energy norm for RUN 4 of Table 4-11 Case VII. 0 -I 1 1 1 1 1 1 1 0 2000 4000 6000 8000 10000 12000 14000 time (sec) Figure 4-18-A comparison of the errors in the midplane temperature for RUN 6 of Table 4-11 - Case VII. 186 1800 1600 % 1400 2 1200 ^ 1000 g 800 P> 600 400 200 0 subcycling i\ classic exact , \ \ \ \ \ \ x 1 1 2000 4000 6000 8000 10000 12000 14000 time (sec) Figure 4-19 -A comparison of the errors in the energy norm for RUN 6 of Table 4-11 Case VII. LOO (cn) : 5.10 (cn) Section B Section A S e c t i o n c Tool -7.20 (cn)-/-R0.46 (cn) / / - R l , 4 6 (cn) n o 3 "O O <+ Figure 4-20 - 2D schematic of the model in Case VIII showing sections along which temperature distributions are sought. 187 Figure 4-21 -ABAQUS thermal elements are used to model both the tool and composite part - Case VIII. Figure 4-23 - Autoclave temperature, degree of cure, and temperature-time history both for cured composite (without exotherm) and curing composite(with exotherm) in the mid-plane of Section C (Figure 4-20) - Case VII. NT 11 - +1 813e+0Z - +1 Bl le+02 - +1 SlOe+OZ - +1 808e+0Z - +1 807e+0Z - +1 80Se+0Z - +1 804etOZ - +1 80Ze+0Z - +1 801e+0Z - +1 ?99e+0Z - +1 798e+02 - +1 796e+0Z - +1 79Se+02 Figure 4-24 — Temperature distribution in the model at exotherm (t=173 min, temperature in degree C) - Case VIII. 189 3 -*-» m 100 a> a. 80 E a H 60 40 20 0 50 100 150 200 time (min) 250 300 Figure 4-25 - Temperature-time history on the top surface, mid surface and bottom surface of the composite from both the current shell element and ABAQUS linear heat transfer elements. These pertain to Section C (as shown in Figure 4-20) - Case VIII. time (min) Figure 4-26 - Temperature-time history on the top surface, mid surface and bottom surface of the composite from both the current shell element and ABAQUS linear heat transfer elements. These pertain to Section A (as shown in Figure 4-20) - Case VIII. 190 Figure 4-27— Through thickness temperature distributions for different sections of the composite ( Figure 4-20) from both the current shell element and ABAQUS linear heat transfer elements at t=45 min - Case VIII Figure 4-28 - Through thickness temperature distributions for different sections of the composite ( Figure 4-20) from both the current shell element and ABAQUS linear heat transfer elements at t=173 min - Case VIII 191 5 SUMMARY, CONCLUSIONS AND FUTURE WORK 5.1 Summary A simple methodology using closed form analytical solutions is developed for thermal analysis of composites during the curing process inside an autoclave. This is based on simple engineering charts developed for a single layer material. The charts developed in this work are applicable to the case of symmetrically applied temperature ramps at the boundaries and constant internal heat generation. Similar classic simple diagrams for the case of temperature hold are adopted from the literature and used in the simple methodology. In the next step, a method is developed for using those diagrams for asymmetric boundary condition and it is shown that this wi l l cause a quantifiable maximum error. This method of solution is also generalized to the case of two-layer tool-part structures and the accuracy of the solution is investigated for the cases of interest. There are some other side issues that need to be addressed so that these solutions become practically useful for process modeling. Some of these side problems that are considered in this work include transition time from steady state of one regime to the next, application of the diagrams to the case of (nonlinear) exothermic heat generation, the edge effect and methods of thermal control for substituting the fin with an equivalent coefficient of heat convection. The second component of this research work is focused on development of a shell element for the thermal simulation of a more general class of composite materials 192 undergoing processing. This shell element uses the shape functions taken from the closed form solutions through the thickness. In this study, however, we only examined the shape functions of one specific set of boundary conditions. The number of these shape functions can vary throughout the analysis leading to adaptive solution features. The main challenge in using these shape functions is the rather cumbersome spatial integration which has been addressed through an efficient spatial integration technique. A n innovative subcycling method has also been used to reduce the computational burden of time-integration in the solution. In the following sections, conclusions of this research, recommendations for future work and a list of the contributions that the author deems worthy of further emphasis are listed in separate sub-sections for both lines of work presented in this research work. 5.2 Conclusion The following are the major conclusions that can be drawn from this research work: Simple Diagrams • Closed form solutions for heat transfer problems and the simple engineering charts generated based on them are practically useful for estimating the temperature distribution in a composite slab undergoing curing in an autoclave. • The simple diagrams for a single layer solid with symmetric boundary conditions can be used for the case of asymmetric boundary conditions and also two layers of material using appropriate strategies developed here. These diagrams could also be used when exothermic heat generation occurs in the material. 193 • Transition time from one steady state solution to the next can also be presented in these simple diagrams. • The above simplified solutions can be used for more general three dimensional cases provided that these solutions are applied to regions that are more than the boundary layer distance away from the edges. Equations needed to estimate the size of the boundary layer region have been established in this study. Higher Order Composite Shell Element • Analysing the thermochemical process during manufacturing using shells provides a seamless transition to structural analysis of composites. • The use of the higher order shell elements with internal degrees of freedom developed here saves time on meshing and remeshing of the spatial domain for processing simulation of composites. • The subcycling technique developed reduces the time integration burden for thermo-chemical analysis involved in process modeling. • The shell elements developed here can render adaptive numerical capability where the number of terms can be activated or de-activated as needed. 5.3 Future work The following are some areas where the author believes further development wi l l be valuable: Simple Diagrams 194 • Setting up a general strategy for using the simple diagrams developed here and also the Heisler charts for multiple layers (more than two layers of material). • Estimation of a general approximation of the error involved in the above mentioned strategy. • Preparing a strategy to deal with tapering solids. • Use of simple solutions in optimization and sensitivity analyses. • Recommendations on size of experimental samples based on the edge effect. Higher Order Composite Shell Element • Use of other boundary conditions than essential (Dirichlet) in developing shape functions for the shell element. • Finding criteria for convergence and number of necessary terms in each simulation in terms of the non-dimensional parameters developed in the chart-form solutions. • Implementing the adaptive solution in A B A Q U S by enabling the software to change the number of degrees of freedom dynamically throughout the analysis. • Developing a special purpose pre-processor for the shell element so that the user has control on defining element direction and implementing convective heat boundaries. 5.4 Contributions The following is a list of noteworthy contributions made in this research work: 195 Simple Diagrams • Development of simple engineering chart for one dimensional thermal analysis of slabs with symmetric convective temperature ramp boundaries and/or constant internal heat generation. • Introduction of new nondimensional parameters for thermal analysis that helps with understanding the heat transfer phenomenon and arriving at simple closed form solutions. • Developing an approximate method for employment of simple diagrams for the case of a slab with asymmetric boundary conditions and estimating the maximum error involved for the general case. • Developing an approximate method for employment of simple diagrams for the case of two-layer slabs and estimating the maximum error involved for cases of interest. • Development of simple methods to analyze the transition time from one regime of heating to the next and also for approximating the edge effects. Higher Order Composite Shell Element • Development of analytical shape functions through the thickness of the shell element. • Development of efficient spatial integration technique. • Development of an innovative and efficient subcycling technique for temporal integration. 196 Implementation of the element as a user defined element in the commercial finite element software, A B A Q U S . 197 REFERENCES [1] Johnston, A., Vaziri, R., and Poursartip, A., "A Plane Strain Model for Process-Induced Deformation of Composite Structures," Journal of Composite Materials, Vol. 35, No. 16, 2001, pp. 1435-1469. [2] Smith, G. D., "Modelling and Experimental Issues in the Processing of Composite Laminates," MAsc. Thesis, The University of British Columbia, 1992. [3] Hubert, P., "Aspects of Flow and Compaction of Laminated Composite Shapes During Cure," PhD Thesis, The University of British Columbia, 1996. [4] Johnston, A., "An Integrated Model of the Development of Process Induced Deformation in Autoclave Processing of Composite Structures," PhD Thesis, The University of British Columbia, 1996. [5] Zobeiry, N., "Viscoelastic Constitutive Models for Evaluation of Residual Stresses in Thermoset Composites during Cure," PhD Thesis, The University of British Columbia, 2006. [6] Arafath, A. A., "Efficient and Simple Numerical Techniques for the Modeling of Process-Induced Stresses and Deformations in Composite Structures," PhD Thesis, The University of British Columbia, 2007. [7] Osooly, A., "Contact Formulation at Tool-Part Interface and Large Strain/Deformation Corotational Stress Based Constitutive Equation in Thermoset Composites during Autoclave Processing," PhD Thesis, The University of British Columbia, 2007. 198 [8] Heisler, M . P., "Temperature Charts for Induction and Constant-Temperature Heating," Transactions of the ASME, Vol. 69, 1947, pp. 227-236. [9] Gurney, H. P. and Lurie, J., "Charts for Estimating Temperature Distributions in Heating or Cooling Solid Shapes," Industrial and Engineering Chemistry, Vol. 15, No. 11, 1923, pp. 1170-1172. [10] Newman, A. B., "Heating and Cooling Rectangular and Cylindrical Solids," Industrial and Engineering Chemistry, Vol. 28, No. 5, 1936, pp. 545-548. [11] Holman, J. P., Heat Transfer, Fifth ed., McGraw-Hill 1981. [12] Suces, J. and Hedge, S., "Transient Conduction in a Slab with Temperature Dependent Thermal Conductivity," Transactions of the ASME, Vol. 100, 1978, pp. 172-174. [13] Fox, E. N., "Temperature Rise in a Heat Evolving Medium," Philosophical Transactions of the Royal Society of London. Vol. 232, 1934, pp. 431-461. [14] Young, F. M . and Savino, C. R., "Time-Temperature Charts for 1-dimensional Conduction with Uniform Internal Heat Generation," Journal of Heat Transfer, Vol. 92, No. 1, 1970, pp. 207-210. [15] Thibault, J., "Analytical Temperature Distribution for Multidimensional Bodies Exposed to a Constant Surface Heat-Flux," Chemical Engineering Communications, Vol. 51, No. 1-6, 1987, pp. 141-151. [16] Zerkle, R. D. and Sunderland, J. E. , "The Transient Temperature Distribution in a Slab Subject to Thermal Radiation," Journal of Heat Transfer, Vol. 87, No. 1, 1965, pp. 117-133. [17] Bowley, W. W. and Koenig, H. A., "Charts for Thermal Transients in Composite Cylinders," Transactions of the ASME, Vol. 93, No. 2, 1971, pp. 248-249. 199 [18] Fasiczka, R. G. and Westermann, E. L . , "Thermal Response of Slabs Subjected to Non-Symmetric Convective Heat Transfer," Journal of Mechanical Design, Transactions of the ASME, Vol. 100, No. 4, 1978, pp. 706-714. [19] Campo, A., "Estimate of the Transient Conduction of Heat in Materials with Linear Thermal-Properties based on the Solution for Constant Properties," Warme undStqffubertragung-Thermo and Fluid Dynamics, Vol. 17, No. 1, 1982, pp. 1-9. [20] Title, C. W., "Boundary Value Problems in Composite Media: Quasi-Orthogonal Functions," Journal of Applied Physics, Vol. 36, No. 4, 1964, pp. 1486-1488. [21] Pagano, N. J., "Generalized Sturm-Liouville Procedure for Composites Domain Anisotropic Transient Conduction Problem," AIAA J., Vol. 12, No. 8, 1974, pp. 1158-1160. [22] Ozisik, M . N., Heat transfer, second ed., JOHN WILEY & SONS, 1993. [23] Haji-Sheikh, A. and Beck, J. V. , "Exact Solution of Heat Conduction in Composite materials and Application to Inverse problems," Journal of Heat Transfer, Vol. 120, No. 3, 1998, pp. 592-599. [24] de Monte, F., "Transient Heat Conduction in One-Dimensional Composite Slab. A Natural Analytic Approach," International Journal of Heat and Mass Transfer, Vol. 43, 2000, pp. 3607-3619. [25] de Monte, F., "An Analytic Approach to the Unsteady Heat Conduction Processes in One-Dimensional Composites Media," International Journal of Heat and Mass Transfer, Vol. 45, 2002, pp. 1333-1343. [26] de Monte, F., "Unsteady heat Conduction in Two Dimensional Two Slab-Shaped Regions. Exact Closed-Form Solution and Results," International Journal of Heat and Mass Transfer, Vol. 46, 2003, pp. 1455-1469. 200 [27] de Monte, F., "Transverse Eigenproblem of Steady-State Heat Conduction for Multi-Dimensional Two-Layered Slabs with Automatic Computation of eigenvalues," International Journal of Heat and Mass Transfer, Vol. 47, 2004, pp. 191-201. [28] Haji-Sheikh, A. and Beck, J. V. , "An Efficient Method of Computing Eigenvalues in Heat Conduction," Numerical Heat Transfer, Part B, Vol. 38, 2000, pp. 133-156. [29] Haji-Sheikh, A. and Beck, J. V., "Temperature Solution in Multi-Dimensional Multi-Layer Bodies," International Journal of Heat and Mass Transfer, Vol. 45, 2002, pp. 1865-1877. [30] Haji-Sheikh, A., Beck, J. V. , and Agonafer, D., "Steady-State Heat Conduction in Multi-Layer Bodies," International Journal of Heat and Mass Transfer, Vol. 46, 2003, pp. 2363-2379. [31] Sun, Y. and Wichman, I. S., "On transient Heat Conduction in a One-Dimensional Composite Slab," International Journal of Heat and Mass Transfer, Vol. 47, 2004, pp. 1555-1559. [32] Lu, X. and Tervola, P., "Transient Heat Conduction in the Composite Slab-Analytical Method," Journal of Physics A: Mathematical and General, Vol. 38, 2005, pp. 81-96. [33] Lu, X., Tervola, P., and Viljanen, M . , "A New Analytical Method to Solve the Heat Equation for a Multi-Dimensional Composite Slab," Journal of Physics A: Mathematical and General, Vol. 38, 2005, pp. 2873-2890. [34] Grober, H. and Erk, S., Fundamentals of Heat Transfer, 2nd ed., McGraw-Hill 1961. [35] Mechanical Engineers' Handbook, 2nd ed., Elsevier 1998. [36] Cengel, Y. A. and Turner, R. H. , Thermal-Fluid Sciences, 2nd ed., McGraw Hill 2005. 201 [37] Cheung, Y . K., "Finite Strip Method Analysis of Elastic Slabs," ASCE Journal of Engineering Mechanics Division, Vol. 94(EM6), 1968, pp. 1365-1378. [38] Cheung, Y. K. and Tham, L. G., Finite Strip Method, CRC Press 1998. [39] Zienkiewicz, O. C. and Taylor, R. L . , Finite Element Method, Volume II, 5 th ed., Elsevier 2000. [40] Manko, Z., "Analysis of Thermal Conduction by the Finite Strip Method," Vol. 2, Pineridge, Swansea, 1985, pp. 1593-1604. [41] Manko, Z., "Thermal Analysis of Engineering Structures by The Finite Strip Method," Canadian Journal of Civil Engineering, Vol. 13, 1986, pp. 761-768. [42] Zienkiewicz, O. C , Gago, J. D. P., and Kelly, D. W., "Hierarchical Concept in Finite-Element Analysis," Computers & Structures, Vol. 16, 1983, pp. 53-65. [43] Zienkiewicz, O. C. and Taylor, R. L . , Finite Element method, Volume I, 5th ed., Elsevier 2000. [44] Akin, J. E . , Finite Elements for Analysis and Design, 1st ed., Academic Press 1994. [45] Babuska, I., Szabo, B. A., and Actis, R. L. , "Hierarchical-Models for laminated Composites," International Journal for Numerical Methods in Engineering, Vol. 33, No. 3, 1992, pp. 503-535. [46] Bose, A. and Surana, S., "Piecewise Hierarchical p-version Curved Shell Finite Element for Heat Conduction in laminated Composites," Computers & Structures, Vol. 49, No. 2, 1992, pp. 283-300. [47] Surana, K. S. and Orith, N. J., "Completely Hierarchical p-Version Axisymmetric Shell element for Nonlinear heat Conduction in Laminated Composites," Computers & Structures, Vol. 46, No. 5, 1992, pp. 777-789. 202 [48] Actis, R. L . , Szabo, B. A., and Schwab, C , "Hierarchical Methods for Laminated Plates and Shells," Computer Methods in Applied Mechanics and Engineering, Vol. 172, 1999, pp. 79-107. [49] Duster, A., Broker, H. , and Rank, E. , "The P-version of Finite Element Method for Three-Dimensional Curved Thin Walled Structures," International Journal for Numerical Methods in Engineering, Vol. 52, 2001, pp. 673-703. [50] Kuhlmann, G. and Rolfes, R., "A Hierarchical 3D Finite Element for Laminated Composites," International Journal for Numerical Methods in Engineering, Vol. 61, 2004, pp. 96-116. [51] Fish, J., Suvorov, A., and Belsky, V., "Hierarchical Composite grid Method for Global-Local Analysis of Laminated Composite Shells," Applied Numerical Mathematics, Vol. 23, 1997, pp. 241-258. [52] Mitchell, J. A. and Reddy, J. N., " A Hierarchical Iterative Procedure for the Analysis of Composite Laminates," Computer Methods in Applied Mechanics and Engineering, Vol. 181, 2000, pp. 237-260. [53] de Frutos, J. and Novo, J., "Element-wise a posteriori Estimates Based on Hierarchical Bases for Non-linear Problems," International Journal for Numerical Methods in Engineering, Vol. 63, 2005, pp. 1146-1173. [54] Zauderer, E. , Partial Differential Equations of Applied Mathematics, second ed., John Wiley & Sons 1998. [55] Davis, P. J. and Rabinowitz, P., Numerical Integration, 1st ed., Blaisdell Pub. Co. 1967. [56] Patterson, T. N. L . , "On High Precision Methods for the Evaluation of Fourier Integrals with Finite and Infinite Limits," Numerical Mathematics, Vol. 27, 1976, pp. 41-52. 203 [57] Evans, G. A., "Two Robust Methods for Irregular Oscillatory Integrals over a Finite Range," Applied Numerical Mathematics, Vol. 14, 1994, pp. 383-395. [58] Evans, G. A. and Webster, J. R., "A Comparison of some Methods for the Evaluation of Highly Oscillatory Integrals," Journal of Computational and Applied Mathematics, Vol. 112, 1999, pp. 55-69. [59] Bleistein, N. and Handelsman, R. A., Asymptotic Expansions of Integrals, 1 st ed., Dover Publications 1986. [60] Briggs, W. L . and Henson, V. E. , The DFT: An Owner's Manual for the Discrete Fourier Transform, lsted. 1995. [61] Heinrich, J. C. and Pepper, D. W., Intermediate Finite Element Method, Fluid Flow and Heat Transfer Applications, Taylor and Francis 1999. [62] Belytschko, T., Engelmann, B. E. , and Liu, W. K., "A Review of Recent Development in Time Integration," State of the Art Surveys on Computational Mechanics A S M E , 1989, pp. 185-199. [63] Brackbill, J. U. and Cohen, B. I., Multiple Time Scales, 1st ed. 1985. [64] Smolinski, P., "Subcycling Integration with Non-Integer Time Steps for Structural Dynamics Problems," Computers & Structures, Vol. 59, No. 2, 1994, pp. 273-281. [65] Wu, Y. S. and Smolinski, P., "A Multi-Time Step Integration Algorithm for Structural Dynamics based on the Modified Trapezoidal Rule," Computer Methods in Applied Mechanics and Engineering, Vol. 187, 2000, pp. 641-660. [66] Smolinski, P. and Pan, J., "Procedures for Multi-time Step Integration of Element-Free Galerkin Methods for Diffusion problems," Computers & Structures, Vol. 77, 2005, pp. 171-183. 204 [67] Smolinski, P. and Wu, Y., "Stability of Explicit Subcycling Time Integration with Linear Interpolation for First-Order Finite lement Semidiscetization," Computer Methods in Applied Mechanics and Engineering, Vol. 151, 2005, pp. 311-324. [68] Sportisse, B., "An Analysis of Operator Splitting techniques in the Stiff Case," Journal of Computational Physics, Vol. 161, 2000, pp. 140-168. [69] Knoll, D. A., Chacon, L . , Margolin, L . G., and Moussaeau, V. A., "On Balanced Approximation for Time Integration of Multiple Time Scale Systems," Journal of Computational Physics, Vol. 185, 2035, pp. 583-611. [70] Lowrie, R. B., "A Comparison of Implicit Time Integration Methods for Nonlinear Relaxation and Diffusion," Journal of Computational Physics, Vol. 196, 2004, pp. 566-590. [71] Rodriguez-Gomez, G., Gonzalez-Casanova, P., and Martinez-Carballido, J., "Computing General Companion Matrices and Stability Regions of Multirate Methods," International Journal for Numerical Methods in Engineering, Vol. 61, 2004, pp. 255-273. [72] Eriksson, K., Estep, D., Hansbo, P., and Johnson, C , Computational Differential Equations, 1st ed., Press Syndicate of the University of Cambridge 1996. [73] Gago, S. R., DE, J. P., Kelly, D. W., and Zienkiewicz, O. C , "A Posteriori Error Analysis and Adaptive Processes in the Finite-Element Method - .2. Adaptive Mesh Refinement," International Journal for Numerical Methods in Engineering, Vol. 19, No. 11, 1983, pp. 1621-1656. [74] Kelly, D. W., Gago, S. R., Zienkiewicz, O. C , and Babuska, I., "A Posteriori Error Analysis and Adaptive Processes in the Finite Element Method: Part I - Error Analysis," International Journal Engineering Science, Vol. 19,1983, pp. 1593-1619. 205 [75] Zienkiewicz, O. C. and Zhu, J. Z., "A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis," International Journal for Numerical Methods in Engineering, Vol. 24, 1987, pp. 337-357. [76] Oden, J. T. and Demkowicz, L . , "A survey of Adaptive Finite Element Methods in Computational Mechanics," State-of-the-art surveys on computational mechanics, edited by A. K. Noor and J. T. Oden ASME, 1989, pp. 441-467. [77] Patra, A. and Oden, J. T., "Computational Techniques for Adaptive hp Finite Element Methods," Finite Elements in Analysis and Design, Vol. 25, 1997, pp. 27-39. [78] Morton, K. W., Numerical Solution of Convective-Diffusion Problems, 1 ed., Chapman & Hall 1996. [79] Eriksson, K. and Johnson, C , "Adaptive Finite-Element Methods for Parabolic Problems .1. A Linear-Model Problem," SIAM Journal on Numerical Analysis, Vol. 28, No. 1, 1991, pp. 43-77. [80] Eriksson, K. and Johnson, C , "Adaptive Finite-Element Methods for Parabolic problems .II. Optimal error-Estimates in L-Infinity -L(2) and L-Infinity-L-Infinity," SIAM Journal on Numerical Analysis, Vol. 32, No. 3, 1995, pp. 706-740. [81] Eriksson, K. and Johnson, C , "Adaptive Finite-Element methods for Parabolic Problems .IV. Nonlinear Problems," SIAM Journal on Numerical Analysis, Vol. 32, No. 6, 1995, pp. 1729-1749. [82] Eriksson, K. and Johnson, C , "Adaptive Finite-Element Methods for parabolic Problems .V. Long Time Integration," SIAM Journal on Numerical Analysis, Vol. 32, No. 6, 1995, pp. 1750-1763. [83] Eriksson, K., Johnson, C , and Larsson, S., "Adaptive finite element methods for parabolic problems - VI: Analytic semigroups," SIAM Journal on Numerical Analysis, Vol. 35, No. 4,1998, pp. 1315-1325. 206 [84] Eriksson, K., Johnson, C , and Logg, A., "Explicit Time-Stepping for Stiff ODEs," SIAM Journal on Scientific Computing, Vol. 25, No. 4, 2003, pp. 1142-1157. [85] Reddy, J. N., Mechanics of laminated Composite Plates, 1st ed., CRC Press 1996. [86] Tamma, K. K. and Yurko, A. A., "Finite-Element Thermal Modeling/Analysis Formulations for Layered Composite Materials," Numerical Heat Transfer, Part B, Vol. 15, 2003, pp. 73-97. [87] Noor, A. K. and Tenek, L. H. , "Steady-State Nonlinear Heat transfer in Multilayered Composite Panels," Journal of Engineering Mechanics, Vol. 118, No. 8, 2003, pp. 1661-1678. [88] Savoia, M . and Reddy, J. N., "Three Dimensional Thermal Analysis of Laminated Composite Plates," International Journal of Solids and Structures, Vol. 32, No. 5, 1995, pp. 593-608. [89] Rolfes, R., Noack, J., and Taeschner, M . , "High Performance 3D-Analysis of Thermo-Mechanically Loaded Composite Structures," Composite Structures, Vol. 46, 1999, pp. 367-379. [90] Noack, J., Rolfes, R., and Tessmer, J., "New Layerwise Theories and Finite Elements for Efficient Thermal Analysis of Hybrid Structures," Computers and Structures, Vol. 81, 20053, pp. 2525-2538. [91] Blanc, M . and Touratier, M . , "A Constrained Discrete layer Model for Heat Conduction in Laminated Composites," Computers & Structures, Vol. 83, 2005, pp. 1705-1718. 207 A P P E N D I X A . E X A C T T H R O U G H T H I C K N E S S E I G E N F U N C T I O N S F O R N O N - D I R I C H L E T B O U N D A R Y CONDITIONS A N D M U L T I L A Y E R PARTS Heat convection on both sides Now consider similar P D E with the boundary condition of the 3 r d kind on both surfaces with the formulation T = T(x,y,z,t) ( A . l ) T,, =ocxTxx +ayTM + azTzz +(l/pCP)H(x,y ,z,t) ( A 2 ) dT h(TT(x,y,t)-T (x, y, +L, t)) = +kz—\ Z=+L h(TB(x,y,t)-T(x,y,-L,t)) = -kz^-\z=_L T(x,y,z,0) = TXx,z) , (AA) where h is the coefficient of convective heat transfer on the surfaces and TT and TB are hd the temperature of the ambient fluid. Defining Biot number as Bi = — , the boundary conditions can be rewritten in the following format: T(x,y,+L,t) +——\z=+L=Tr(x,y,t) Bi dz L dT T(x,y,-L,t)-——\z=_L=TB(x,y,t) . Bi dz (A.5) 208 The following change in the unknown function is necessary to get the homogeneous boundary conditions: T(x, y, z, 0 = fi-^- + . T* ^ 7) + u(x,y,z, t) , 2 2(1 +1 / Bi) L (A.6) which leads to U t = ( a x U xx + ayU yy + «zW zz ) + 0 / J>» Z> 0 — (TT.—GCTT. —GC T,R ) \ T,( TjXx y T,yy )\ ~{^B,1 ~ax^B,xx ~ay^B,yy) 1 —+ -1 2 2(1 + 1/50 L 1 2 2(1 + 1/5/) L (A.7) and the following homogeneous boundary conditions: _ , L du 1 « ( x , ^ , + I , 0 + — — z = + i =0 Bi oz L du 1 u(x,y,-L,t)- ——\Z=_L =0 , Bi oz (A.8) and initial conditions (A.9) to find the shape functions through thickness, we should solve the following eigenvalue problem this time: U , , = a x U , x x + a y U , y y + a : U , * . L du 1 u(x,y,+L,t) +—— z = + i =0 5/ oz r x L du, u(x,y,-L,t)-——\z=_L=Q. Bi dz (A.10) 209 Similar to previous problem, we apply the following solution based on separation of variables: which leads to which means that u(x,y,z,t) = W(x,y,t)Z{z) , ( A . l l ) W,Z = a ^ Z + ayWtyyZ + WZ" W, W Wm Z " ( A - 1 2 ) W x W y W ^ - = - A 2 , (A.13) Z a p is some constant. The solution to the above O D E is: Z = ylsin(/?-) + 5cos(/?4) • ( A - 1 4 ) d d The boundary conditions, on the other hand would translate to B l d z (A. 15) V ' Bi 8z |z" L Applying the boundary conditions we get •0) _ cos sin L , j is odd j is even (A. 16) 210 In order to get { / ^ ] the following equations should be solved in the interval c o t g ( / ? w ) = ^ — j is odd Bi , Aif/K ~Bi . . cot g(/?J>) = —^-j j is even (A. 17) The general solution to u can then be written as u=^WV)(x,y,t)Zv'(z) (A. 18) Inputting the above in equation (A.7) we have W}J) =axW%)+ayWV'-aI 0). v L J W(i) + h{j)(x,y,t) (A.19) J (1 / pCP)H(x, y, z, t) Z(z)dz h(j) (x, y, t) = ^ T l ]z(i\z)Z[i\z)dz -L +L K T T - a TT -a TT ) 1,1 x T,xx y• T,yy J f \ 1 z A + -I 2 2(1 + 1/5/) d Z(z)dz \z[i\z)Z[i\z)dz \ ( T B , , - a X T B ^ - a y T B , y y ) 1 1 Z 2 2(1 + 1/5/) d Z(z)dz \z{i\z)Z{J)(z)dz (A.20) The above problem is to be solved using the following initial conditions: 211 W(x,y,0): +L \ -L TT +TB TT TB z^ v 2 2(\ + \IBi)d)io Z(j)(z)dz +L \zU\z)Z(i\z)dz (A.21) -L Finite Element solution to the P D E in equation (A. 19) is very similar to the one described in the previous section. Dirichlet boundary condition on one side and convection on the other side In this case the P D E and the boundary conditions are as follows: T,t =axT^x +ayT^ +azT^ +(1/pC P)H {x ,y ,z ,t) (A.22) T(x,y,+L,t) +——\z=L=TT(x,y,t) Hi oz (A.23) T(x,y,-L,t) = TB(x,y,t) T(x,y,z,0) = Ti(x,y,z) (A.24) (A.25) TT is the fluid temperature on the top and TT is the essential boundary on the bottom. The following change in variables could be used to change the boundary conditions to homogeneous boundary conditions: T(x,y,z,t) = TT+TB{\ + l/Bi) | TT-TB z 2 + 1/5/ 2 + 1/5/ d + u(x,y,z,t) . (A.26) and the P D E problem changes to 212 U,i ~ \ a x U , x x ~^~ayU,yy ~^~azU,zz ) + pCP)H(x,y,z,t) ~(TT,, ~ axTT,xx ~ ayTT,yy ) ( " + 1 2 + 1/5/ 2 + 1/ Bid 1 z J r - T - T \(l + l / B i V" " ' M 2 + 1/5/ 2 + 1/5/ d (A.27) i <. L du | n u Ix,y ,+L,t) + 1 _., =0 v ; Bi dz u~+l u(x ,y = 0 (A.28) (A.29) i(x,y,z,0) = Ti(x,y,z)-TT+TB{\ + l/Bi) | TT-TB z 2 + 1/ Bi 2 + 1/Bid J,=o (A.30) Separation of Variables, as in previous sections, wi l l lead to Z (A.31) Z(+Z) + - ^ U = 0 Bi dz Z ( - I ) = 0 (A.32) (A.33) The solution to the above O D E is Z =A sin f z ^ l + 5 cos f z^ (A.34) Applying the boundary conditions we have Z 0 ) ( z ) = cos( /? 0 ) )s in( /? 0 ) - ) + sin(/? ( y ) )cos(/? 0 ) —) (A.35) 213 In order to get {/?^] the following equation should be solved in the interval [(2n-l);r,(2n+1)*]: (A.36) Here is the problem that should be solved for each j WV) =axwV} +ayW(jy)-a2 V L J W + hv)(x,y,t) (A.37) hU\x,y,t) = | ( 1 / pCP)H (x, y, z, t) ZiJ) (z)dz -L +L \z(i)(z)Z(i)(z)dz +L 1 - + -2 + 1/5/ 2 + 1/5/ l—L\zW(z)dz MBid) v ' +L \Z(j\z)Z(s)(z)dz -L +i f \(TB>l-axTBxx-ayTByy) -L \ ( I + I/5/3 1 2 + I/5/3 2 + l/5/3 dj Z(j)(z)dz +L jZ°\z)ZU)(z)dz The above problem is to be solved using the following initial conditions: (A.38) W[j\x,y,0) = +L \Tini(x,y,z)--L TT+TB{\ + \IBi) | TT-TB z 2 + 1/ Bi 2 + \/Bid Zu\z)dz +L \z(i\z)Z(i\z)dz (A.39) The rest of the solution is exactly as in previous sections. 214 A multilayer problem with essential boundaries Embedding a multilayer composite inside an element has been of interest to researchers for a while. Reddy has described the way to deal with stress problems in his book: "Mechanics of Laminated Composite Plates" ([85]). Other researchers have studied the multilayer thermal element. Tamma and Yurko ([86]) investigate the problem of multilayer element without implementing the continuity of flux at layer interfaces. Noor and Tenek ([87]) have implemented first order thermal lamination theory to their element and have used the element for nonlinear cases. Savoia and Reddy ([88]) have used a C° continuous shape function for thermal analysis of a laminated composite and have used the result in a thermoelastic analysis. Rolfes et al. ([89], [90]) have also investigated the use of zigzag shape functions in the through layer direction of a composite. Blanc and Touratier, in a recent paper ([91]) have used constrained discrete layer model to describe the through thickness distribution of temperature in a laminated composite. In this section, shape function based on exact through thickness solution of a multilayer composite with essential boundary condition is determined to be used in multilayer shell elements. Multilayered heat problems are also of grave interest in modelling of composites. A 2D cross section of a typical three-layer hybrid parallelepiped is illustrated in Figure A - l . General N-layer problem is mathematically modelled as follows: 215 T, = axT^ + ayTw + azTzz + (l / pCP ) H (x,y, z, t) T(x,y,0,t) = TB(x,y,t) ( N f/l ^ T x,y,^2L['l,t =TT(x,y,t) v 1=1 J T(x,y,z,0) = Ti(x,y,z) (A.40) where the thermal material properties kz, H change at the layer interfaces and TB and TT are the bottom and top essential boundary values. The temperature profile is contentious at the interfaces and its first spatial derivative satisfy the following continuity condition: dT ~dz~ 2 dz @ z = 2L[x\l(d{x] + d[2]p... (A.41) Defining a separate z coordinate for each layer wi l l make the mathematical manipulations easier to grasp so we define z l ' J = z - 2 W=i J (A.42) in which case a layerwise temperature distribution function T ^ [x ,y ,z ^'\t j would describe the temperature distribution through the thickness of the itn layer. The constitutive, boundary and initial value equations ^ ^ would then change to = axT® +ayTl + c$T® + pCP)H^(x,z[i],t) T1(x,y,-L[l],t) = TB(x,y,t) TN (x, y, ,t) = TT (x, y, t) T[i](x,y,z,0) = Ttt(x,y,z[i]). (A.43) The continuity equations would then read 216 V (A.44) for i = l , 2 , . . . , W - l . Going through the same line of arguments as in the previous sections, we try to change the boundary conditions to homogeneous boundaries by the following transformation: 7< 'kx ,2 l^O = ( ^ ( ^ 0 ^ + > 3 O ( ^ ) 2 i ) + « I ' 1 ( * , A / ) • (A.45) Here / r [ , ] and / j ' 1 are linear functions that satisfy the boundary and interface conditions. These functions are defined in terms of the following parameters: Sl>] 2 [ / ] = S - T T n=l S ['] (A.46) which results in J B ~ + • + 1 [N] 5 [ / ] z [ J V ] +1 V " J + (A.47) The P D E would then change to: «!?=(«,«£+«A+«fJ«S)+0 / ^ ) " [ , ] (* • y>zW'') -{TT,,-dxTTxx - « / r , w ) / r [ ' ] - ( ^ , -«?B„)fln . With the following boundary, interface and initial conditions (A.48) 217 u^x,y-L®t)=uW(x,y,+LWt) = 0 ul'](x,+Ll'\t)=ul'+l](x,y,-Ll'+l],t) i =\,2,...,N-\ (A.49) (WlA = ( * J ' + V ' $ ) i=\,2,..,N-\ «I ' l (x ,^zW,0) = r i ( J c , z W ) / x ; x ( A - 5 ° ) - / J ' 1 (* w ) * i (*, * o) - f?] (*") rr (x, y, 0) . We now try to get the Sturm-Liouville solution for the following homogeneous P D E : = a x U , x x + a y U , y y + «z M,zz » = 1, 2 , 3 , N , (A.51) with interface and boundary conditions as defined in (A.49). As it is usually the case with similar problems, we write i/I'l =W (x ,y ,t ) Z 1 ' ] ( z [ , ] ) (A.52) which leads to Z[i]=Ali]cos(r]['](z[i]-S[i]f) . (A.53) Now input(A.52) and (A.53) into (A.51) to get W,=axW^+aF„-ay (IT1'1)2 w, w„ w, — = a + a.. W x W y W which means that (A.54) \2 I I I III a =P2 (A.55) 218 n[i]=plJa? . (A.56) On the other hand, applying the boundary conditions to (A.53) we get: A^cos(^(-d^-S^)) = 0 A^cos(?j^(+d^-S^)) = 0 ^ " c o s ( i 7 [ # 1 - c o s ( ^ ' + , l - = 6 i =\,2,..,N -1 ^ ' ^ ^ V s i n ^ ' 1 ^ ^ ^ ^ i =1,2 -1 or Bi cos(?7 1 J 1 )-C| s'm{r]ldl) = 0.0 BN cos(?iNdN)+CN sm(nNdN) = 0.0 ^ V ' s i n ^ ' V ' ^ ^ (A.58) ^ ' - ^ J W ( A,9 ) ^Ucos-'fflW/i'1). In order for equations in (A.58) to have a non-trivial vector of solution so determinant of the coefficient matrix should be identical to zero, which in combination with equation (A.56) leads to the value for j /?^ '} of different mode shapes which result in j ? / ' ^ j . The eigenvalues of the problem inserted into (A. 59) results in for different mode shapes. Here is the problem that should be solved for each j this time: 219 W}J) =ccxWi?+ayW^-(/3(J))2W + h(j)(x,y,t) (A.60) N +4" h(J)(x,y,t) = J \(\lpCp)H^.Z{i)^ '=' -z!'l The Finite element formulation would be the same as in previous sections. (A.61) . '=' -zl'l '=' ./H 1=1 _jrJ0 The above problem is to be solved using the following initial conditions: W{j)(x,y,0) = ^ ^ -pj . (A.62) 220 1 Layer 3 J T(z) Layer 2 ( Layer 1 AZ V 1 . X 1 p Figure A-l — Typical temperature distribution in a three-layer composite. 221
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Efficient methods for non-linear thermochemical analysis...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Efficient methods for non-linear thermochemical analysis of composite structures undergoing autoclave… Rasekh, Ali 2007
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Efficient methods for non-linear thermochemical analysis of composite structures undergoing autoclave processing |
Creator |
Rasekh, Ali |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | Composite structures are increasingly being used in different industries. Their manufacturing imparts some challenges for the industry: most importantly prediction and control of the process to specification. The usual numerical solutions typically based on the use of the finite element analysis are not very suitable for large parts, especially when there is a need for quick estimation of the results for preliminary design and optimization. Therefore, there is a need both to have enhanced solutions that reduce the modeling effort for computer simulation of large and complex structures and also to simplify the solution and provide easy to use methodologies for quick estimations based on tables and charts. In the present work, a simple methodology is developed to estimate the temperature distribution in a thermoset polymer matrix composite slab placed on a tool and subjected to cycles of temperature ramp and hold leading to the curing of the composite and generation of heat due to the internal chemical reactions. Supplementary diagrams are also generated to set limits on the method. A modified finite element solution for heat transfer is also introduced that reduces the mesh generation and computational effort. This "higher order shell element" uses enhanced shape functions and efficient methods for spatial and temporal integrations in order to reduce the computational run times. The developed methods provide the design engineer with efficient analysis tools for predicting the temperature in a composite part. The simple diagrams and tables can be used for preliminary estimation of the temperature distribution in the part at each stage of the material development. The enhanced finite element methodology developed here can be used to reduce the amount of effort necessary in mesh generation and refinements necessary to achieve accurate solutions for thermochemical modelling of complex composite structures during processing. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-02-03 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0063235 |
URI | http://hdl.handle.net/2429/31050 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2007-267790.pdf [ 7.3MB ]
- Metadata
- JSON: 831-1.0063235.json
- JSON-LD: 831-1.0063235-ld.json
- RDF/XML (Pretty): 831-1.0063235-rdf.xml
- RDF/JSON: 831-1.0063235-rdf.json
- Turtle: 831-1.0063235-turtle.txt
- N-Triples: 831-1.0063235-rdf-ntriples.txt
- Original Record: 831-1.0063235-source.json
- Full Text
- 831-1.0063235-fulltext.txt
- Citation
- 831-1.0063235.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0063235/manifest