Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Efficient methods for non-linear thermochemical analysis of composite structures undergoing autoclave… Rasekh, Ali 2007

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

Efficient Methods for Non-Linear Thermochemical Analysis of Composite Structures Undergoing Autoclave Processing by ALI RASEKH 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 F O R 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 S T U D I E S (Civil Engineering)  The University of British Columbia April 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  T A B L E OF CONTENTS  ABSTRACT  II  T A B L E OF CONTENTS  IV  LIST O F FIGURES  VIII  LIST O F T A B L E S  XIII  ACKNOWLEDGEMENTS  XIV  DEDICATION 1  XV  INTRODUCTION AND B A C K G R O U N D  1  1.1  Autoclave Processing and Process Modeling  2  1.2  Research Objective and Thesis Outline  3  2  SIMPLE D I A G R A M S F O R T H E R M A L ANALYSIS  8  Nomenclature  8  Greek symbols  9  Subscripts  10  2.1  12  Introduction  2.2 One Layer of Material - Symmetric Boundary Conditions 2.2.1 Statement of the Problem 2.2.2 Applied Ramp Temperature 2.2.3 Applied Hold Temperature 2.2.4 Engineering Charts  14 14 14 21 21  2.3  22  Transition  2.4 One Layer of Material - Asymmetric Boundary Conditions 2.4.1 Statement of the Problem 2.4.2 Applied Ramp Temperature 2.4.3 Applied Hold Temperature 2.4.4 Approximate Solution  24 24 25 31 32  2.5 Thermal Analysis of a Composite-Tool Structure Using Simple Diagrams 2.5.1 General Formulation for Ramp Case 2.5.2 Applied Hold Temperature 2.5.3 Approximate Solution  35 .37 44 45  2.6  47  Case studies 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 2.8.1 Case of Hold Boundary 2.8.2 Case of Heat Generation and/or Ramp Boundary 2.8.3 Asymmetric Boundaries and Two Layer Solids 2.9 Methods of Thermal Control 2.9.1 The Longitudinal F i n with a Rectangular Profile  50 51 54 56 56 57  2.10 Summary 2.10.1 Summary of the Approximate Solution Algorithms 2.10.2 Future Development  59 59 65  3 D E V E L O P M E N T O F A H I G H E R ORDER S H E L L E L E M E N T F O R 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  91  Introduction  3.2 Constitutive Equations 3.2.1 Heat Conduction and Diffusion in a Solid 3.2.2 Constitutive Equation for Exothermic Heat Generation  92 92 96  3.3 Parallelepiped Shell Element Formulation 3.3.1 Dirichlet Boundary Conditions  97 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 3.6.1 Time Integration of Heat Transfer Equation 3.6.2 Time Integration of Cure Kinetics Equation  116 117 119  3.7 Nonlinear Procedures 3.7.1 Newton-Raphson Iteration 3.7.2 Successive Substitution Method  120 121 124  3.8  Energy Norm  126  3.9  Errors and Adaptive Solution  127  3.10  Evaluation of Material Properties at Integration Points  3.11 Program Algorithm 3.11.1 Adaptivity Module  129 131 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 4  Summary and Future Work  147  C A S E 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 I V : 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 A B A Q U S 1 6 4  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  169  5  Summary and Future Work S U M M A R Y , CONCLUSIONS AND F U T U R E 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 E I G E N FUNCTIONS F O R NONDIRICHLET B O U N D A R Y 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 Bi and hj I h 70 B  B  Figure 2-8 - The global maximum errors extracted from curves in Figure 2-7 plotted as a function of Bi 70 B  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 Bi and hj I h . B  B  71 Figure 2-10 - The global maximum of curves in Figure 2-9 as a function of Bi  B  Figure 2-11 - A as a function of hj I h for 0.001 < Bi < 0.5 B  Figure 2-12 - A as a function of hj. lh  B  72  B  for 0.5 < Bi < 100  72  B  Figure 2-13 - Maximum values of A for all hj lh  B  71  as a function of Bi  73  B  Figure 2-14 - Maximum error of strategy #2 in transient phase for autoclave temperature ramp; the error is presented in terms of hy I h for 0.001 < Bi < 0.5 73 B  viii  B  Figure 2-15 - Maximum error of strategy #2 in transient phase for autoclave temperature ramp; the error is presented in terms of hj I h for 0.5 < Bi < 100 74 B  B  Figure 2-16 - The global maximum of curves in Figure 2-15 and Figure 2-17 versus Bi . B  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 h for B  0.001 < Bi < 0.5  75  B  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 h for B  0.5<5/  B  < 100  75  Figure 2-19 - The global maximum of curves in Figure 2-17 and Figure 2-18 versus Bi . B  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 Figure 2-26 - %  79  as a function of Biot number for different error tolerances and for  ed  K/K  2  Figure 2-27 - g  =0.01  79  as a function of Biot number for different error tolerances and for  ed  K/K  2  Figure 2-28 - %  =0.05  80  as a function of Biot number for different error tolerances and for  ed  K/K  2  Figure 2-29 - %  =0.1  80  as a function of Biot number for different error tolerances and for  ed  K/K  2  =0.2  81  ix  Figure 2-30 -  as a function of Biot number for different error tolerances and for  %  edge  K/K =0.5  81  2  Figure 2-31 - £,  as a function of Biot number for different error tolerances and for  rf  K/K =\.0  82  2  Figure 2-32 - g  as a function of Biot number for different error tolerances and for  ed  K/K =\0  82  2  Figure 2-33 - g  edge  as a function of Biot number for different error tolerances and for  K/K  =100  2  Figure 2-34 - <^ _ edge  max  83  as a function of K/K  for different error tolerances. Markers show  2  the maxima from the previous diagrams and trend lines are drawn through them 83 Figure 2-35 - g  ed  as a function of K/K  2  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  and through thickness discretization points (  , <^ , £3 > O 2  3  defined in Section  3.5. This is drawn at integration point KGP with the thickness equal to  D . KGP  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 9 - Case II  178  Figure 4-2 - Magnified section of Figure 2-4 used to read 9 - Case III  178  m  m  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 I V 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 I V 179 Figure 4-5 - Error in exotherm prediction for single layer material using simple charts CASE 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 CASE 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 V I . 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 V I 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 V I . 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 V I I 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 p a r t - C a s e 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 o f 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  xii  221  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 D r 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 civil engineering department for their friendship and help. Specially, I would like to thank M r . Ahamed Arafath and M r . 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 tolovely mom/ cwid/ dx&d/  XV  1  INTRODUCTION AND B A C K G R O U N D  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 thermochemical 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 o f 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 o f several years o f 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 o f 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  Bi  Biot number in surface direction: Bi -  2  2  Bi  equivalent Biot number  C  circumferential length of a fin  e  E  hLjK  2  coefficients of the transient series solution  t  E  max  non-dimensional maximum value of the error of maximum lag (lead)  E  G  maximum value of E  Fo  Fourier number:  C  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  h  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  P  ss  max  for  \<h /h <eo T  B  aljL  2  K  thermal conductivity in surface direction  L  half thickness of slab  2  L  e  effective half thickness of the slab  n  number of fins in unit length  r  temperature ramp in the fluid  fm  8  (ti-r)l  R  R number:  S  thermal conductivity per unit length: K/L  T  slab temperature  T  err  if)  m  ja  e difference between the exact and approximate temperature lag in the  part T  lag  (t)  T  err  ^hoia  T  temperature lag in the part T  err  s t e  (t) at the moment when \T (t)/T err  lag  (t)\ maximizes  P change in autoclave temperature  initial temperature of the system  0  T  temperature at the base of the fin  T  initial temperature at the base of the fin  T ,ocia,e  autoclave temperature  T  fluid temperature  T  steady state temperature distribution in the slab  b  b0  au  x  ss  T  non-steady part of the transient solution  u  non-dimensional temperature distribution function through thickness  Tr  u  non-dimensional temperature distribution function in surface direction  z  through thickness coordinate  Z,  eigenfunctions  2  Greek symbols a  diffusivity in through thickness direction  a  diffusivity in surface direction  2  P  parameter in estimating steady state Biot number  p  density of the material  m  non-dimensional length along a fin  r\  non-dimensional eigenvalues for each layer  S.  non-dimensional parameter in eigenfunctions  P  i  9  d  non-dimensional location of the adiabatic line with respect to the midplane  A  width o f a fin  p  X  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  s  £  location of adiabatic line in the first layer  a  £  the coordinate in the second layer in 2 layers section  %  non-dimensional coordinate in surface direction in edge effect section  £,  location of adiabatic line in the second layer  a  <^  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)  9  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 E  (5  non-dimensional parameter  A  error involved in estimating the location of min (max) T by S  ed  m  P  A  i  p  G  Y  maximum value of A for \<hj jh  B  percentage of steady state achieved  Subscripts B  bottom surface of the slab  T  top surface of the slab  10  < °o  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 o f 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 will be introduced in the simple equations developed for the steady state as well as the transient solutions. Secondly, we will 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 w i 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 will 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 2.2.1  One Layer of Material - Symmetric Boundary Conditions 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 8T_  U  2  M  - + H,  = ~{T-T ) A.  \£\<l  (2.1)  o n £ = +l,f>0  m  8T  t>0,  2  (2.2)  hi  = +—(T-T„)  onC =  T = T +rt. 0  oo  14  -\,t>0  (2.3)  Here H is the adiabatic temperature rate, defined as H =  H, H is exothermic heat  generated per unit volume, p is the density and C is the thermal capacity. p  * = Tamodave  T  i n  t h e  a  b  o  v  e  equation is the autoclave temperature, r = t  auloclme  is the  constant autoclave rate of heating (linear temperature ramp), and T is the initial 0  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 pC  P  coefficient, which is identical for both top and bottom surfaces. Finally, C, is the nondimensional 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: T =A{£)t  + B{C)  ss  (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 T  ss  15  can be written as:  57;  = A't + B'  d%SS_  A" A"t + B"  (2.5)  K 8Z ss 1  dt  where A' = dA/dC  and A" =  d A/dC 2  2  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  'a}  (2.7)  B" = A-H  Upon solving Equation (2.7) we obtain: A=CQ+D B = y Cf+y (D-H)£  + E<;  2  6  2  +  F,  (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-T ) o0  A't + B' = +Bi(At + B-T ) x  where Bi is the Biot number and is defined as:  16  on £ = +1 on £ = -l  (2.9)  (2.10)  K  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-T )  on ^ = +1  B' = +Bi(B-T )  on  0  0  (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:  (2.13)  Substituting (2.13) reduces Equations (2.12) to the following: ( 1^ fa +F= +E 1 + Bi, V U2 /  -E  1+ V  1^ Bi  y  +F  —  fa  U  2  \  {R@ + T ) 0  /  (2.14)  \  J  (Re + T ), 0  where R,® are defined as  (H-r)L  2  a 0 = 0.5 + — . Bi  17  (2.15)  The solution to Equations (2.14) is: E=0 a  F=  (2.16)  (R® + T ) 0  Therefore, we have  B = ~RC +R® 2  So that the steady state temperature T  ss  + R\  0  0  becomes  ss  T =T +rt  (2.17)  +T  U +® 2  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 =  T +T =T„+Rss  Tr  1  <Z +® 2  + TTr  (2.19)  Upon substitution into Equations (2.1)-(2.3), we obtain f a \ dT 2  dt 8TTr  %-,\$\<Ut>0  vL j 2  = -Bi-T  (2.20)  on ^ = +1, f >0  Tr  (2.21) = +Bi• T  Tr  on £ = -\ , t>0  18  This problem is to be solved subject to some arbitrary initial condition.  7'  7 >  =£E Z (^)exp(-^ ) /  /  (2.22)  0  1=1  where,  Z, =  [cos(40  i = 1,3,5,7,.  Isin(^)  / = 2,4,6,8,.  (2.23)  and {4} is the set of eigenvalues obtained from the following equations:  Bi +— for A. tan ( 4 ) = r  I  - A f 5/  i = 1,3,5,7,. (2.24) ,- = 2,4,6,8,  or  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 T , that Tr  T  Tr  =T -T„  n  =-R-  a  1 ~2 —C+0  (2.26)  or  -R-  ~C +e] = SE,z,(0 2  (2.27)  J 1=1 Now using the orthogonality of { Z , } , we have  E, = - i ?  — V, +©£•,  2 '  19  (2.28)  where  ^1  J(z,) *T 2  -i  s. =  (2.29) -1  J(z,) <*r 2  Substituting all these into Equation (2.19) leads to:  rp  -I  rjn  00  U, + ©£, Z , ( C ) e x p ( - A , F o ) 2 ' 2  (2.30)  For the case of mid-plane temperature where C, = 0, the equation above is simplified to  - =  0  ~ ^r = -t(-U )  T  jn  L  &  +&£  < (°) e x p ( - A ^ )  Z  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 +  20  Bi  (2.32)  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„=T  (2.33)  X  The solution is well established in heat transfer literature and can easily be calculated  ^ -'o  = 2 > ' ' (C)exp(-A,. Fo) Z  2  (2.34)  ^ 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 will 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) will 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 R , we 2  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:  (2.36)  7>l,=o  J  In the midplane surface we have ( ^ = 0):  (2.37)  Therefore, the solution for E, is modified to: 22  (2.38)  -O, +&£  E,=-(*2-*.)  i  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: 1  T -T^R ®-(R -R ) m  2  2  1  f m  ~ mSS  T  T  V, + &£,  1 '  x  =-{ 2- \) R  R  v 2  ^  V,  +&£.  Z,(0.0)exp(-^ Fo) 2  (2.39) Z, (0.0)exp(-A? Fo) .  1  Comparing with Equation (2.37) we obtain ( Z , (0.0) = 1.0):  1.0-Y/100:  T  m  1  T  — u, + 0 f , e x p ( - ^ Fo) _2 )_ 0  -T  mSS  1  -T  mO  (2.40)  mSS t=0  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: f  Foj -  -1.0  In  \  \  (1.0-Y/100)©  -1.0  0  l n ( l . 0 - Y / 1 0 0 ) + ln  1  A  2  v  2  V, 1  (2.41) +©£•. 1  J  )  We can now consider the extreme cases. For large Biot numbers, from Equations (2.24) f  ^  0  and (2.32) we have \= — , 0 = 0.5, and In ^"2 Therefore 23  u, + 0 f ,  3.1546x10"  Fo =-0.4052847 x l n ( 1 . 0 - Y / 1 0 0 ) + 0.01278521. T  (2.42)  On the other hand, for the cases with small Biot number we have A s yfWi, 0 s — and Bi 4  0  In V 2  1  = 0.0. This gives that 1  Fo = r  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 2.4.1  One Layer of Material - Asymmetric Boundary Conditions 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 h < h unless otherwise mentioned. The solution is obtained for B  r  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: ( a\  dT  dt  ^ d£ ^L  =  V-  dT 2  dt  2  J  1  +  H,\C\<\,t>0  _ M ( _ r ) n 4 - = i,/>o K r  0  (2.44)  +  v  =  hsL( -T )  +  T  on C =  a>  (2.45)  -ht>0  (2.46)  Here h ,hj. are the convective heat transfer coefficients for the top and bottom surfaces. B  This system reaches a steady state condition in which the temperature at each point of the solid changes linearly with time, that is: T =A{t)t  +  ss  (2.47)  B{t).  In order to calculate the spatial functions A and B, T  ss  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) (2.50) B = / C<; +y (D-H)C 3  + E<; + F  2  6  2  where C, D , E , F are integration constants. Similarly, T  ss  is substituted for T in Equations (2.45) to obtain  A't + B' = -Bi {At  + B-T ),  on £ = +1  A't + B' = +Bi (At  + B-T„),  on ^ = -1  T  m  B  (2.51)  which yields A' = -Bi (A-r),  on£" = +l  A' = +Bi (A-r),  on£ = - l .  T  B  Here Bi ,Bi B  T  (2.52)  are Biot numbers corresponding to the bottom and top boundaries and are  defined as  B i  B  Bi =  = ^ K hjL  T  (2.53)  K  Now using Equation (2.49) leads to (2.54)  A =r Another set of equations derived from Equation (2.51) is: B' = -Bi (B-T ),  on £ = +1  B' = +Bi (B-T ),  on £ = -1  T  0  B  0  26  (2.55)  Upon using Equations (2.50) and (2.54), we obtain  1 ^ +F = +E 1 + Bi J V "*T f  -E  (  1  0.5 +  I  T  Bi T  V  \  1+ +F = Bi B J  R 0.5 +  1_ BL  +T  n  J  (2.56)  \  +T  0  J  J  where R = (H -r) — , as defined previously in Equation (2.15). Solving the system of a Equations (2.56), we obtain  R8„  E=  (2.57) T + R 0.5 + — 0,p  F =  where /? , 8  0  J)  are defined as:  1  1  2  + — +  ^  Bi  T  —  Bi  T  B  T  8= P  -  Bi .Bi . i i 2+— + — Bi Bi B  a  \ \ • 2+— + — Bij Bi B  l  (2.58)  B  (2-59)  B  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 + — +  (2.60)  R5 Z-V R?.. P  )  PPJ  V  which can be written as T  -T  (2.61)  0.5 + -  R  PP  The minimum value for the right hand side of Equation (2.61) occurs at C, = 5 . Shifting p  the origin of the coordinate system to this point leads to (£" = 5 + %):  T  -T  SS  <a  1  1  {  _  R  1^ 0.5 + —  (2.62)  which simplifies to  T -T *ss - ° o R  0.5 + ^  + ^- - V PP 2  l  (2.63)  y  2  Now Bi. is defined as  1  1  Bi  - ± PP  g +  - l 2  (  1  + Bi + Bij .Bi Bi ~ i i 2+ +Bi-,. Bir, B  B  which simplifies to  28  1 +—  2  1  1  Bij  Big  r 2+ + — V BL BiB J n  •  (2.64)  Therefore, Equation (2.63) can be rewritten as:  (2.66)  00  R  X in the equation above varies in the range \ -\-5 ,+\-S p  symmetric in % , and 5 changes in the range  p  J . A s Equation (2.66) iis  , 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  [ 0 , l - £ ] to the range above this line. 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:  (2.67)  Substituting the above summation into Equations (2.44)-(2.46) yields 29  87^  8%  (2.68)  dt  571  — = —Bi • T T  Tr  +l,t>0 (2.69)  dT Z- = +Bi -T 7  B  on t =  Tr  on £ =  -\,t>0  with initial condition T (<£", 0) = T . 0  This problem also has a classic solution based on Fourier series, given by  ^=lE,Z,(C)exp(-^Fo) , /=i  (2.70)  Z , = cos (A, ( £ - * , ) ) .  (2.71)  where  Here {/I,} is the set of eigenvalues attained from the following equation:  -A sin {Aj) + Bi cos (A,)  +A cos (A.) + Bi s'm(A )  +Aj sin (A ) - Bi cos (A,)  +A, cos ( A,-) + 5 / sin ( l , )  i  det  r  l  where A e  v  y  2  2  B  i  T  j  = 0  (2.72)  B  and each A is associated with a  from the following  i  equation:  tan(^,) = -  -A sin (A ) + Z?/;. cos (A,) t  t  +A, cos (/I,) + BI sin ( l ) T  (2.73)  (  In order to determine the coefficients {E,} we need to apply the initial condition for TTr 7  This yields,  30  •  Tfr ,  =0  -  T  ss  -  (=0  (2.74)  R • 0.5+—+SX- Ac l  2  Upon using the orthogonality of {Z,.}, we obtain the coefficients E, as  r  E, = -R V  1 ' — v + 0.5 + t  (2,75)  PP  2  where  (2.76)  jz,.z,^ Substituting all these into Equation (2.67) leads to: ryi  rji  1  7?  CO  f  V  i 2 '  r„. i ^ 0.5 + — £ , + < y , Z , ( C ) e x p ( - ^ F o ) v PPJ A  (2.77) The right-hand-side of Equation (2.77) depends on pairs of Biot numbers  (Bi ,Bi ). B  r  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  31  (2.78)  The solution in this case is: T-T.  = £*, ,(0exp(-4 F ) . z  2  O  (2.79)  The parameters and functions involved in this equation are taken from Equations (2.71)(2.77). 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.80)  E,max  32  where T is the minimum (maximum) temperature within the thickness (i.e. at the m  stationary location of the adiabatic line) and T is the approximate value of the same m  state variable using the equivalent Biot number in Figure 2-3. The error is shown for different Bi for a wide range of hj I h . B  B  Figure 2-8 shows the maximums for each Bi or: B  E =max  (E  G  m a x  )  Figure 2-9 and Figure 2-10 show similar error for the case of hold temperature. E  (2.81)  m a x  is  defined slightly different in this case, that is  max (|r - Q^°°  E max  The definition of E  G  -fj) !Z  VI  It*  (2.82)  T*  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 Bi = 0.0. B  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 (S ) of the first trigonometric term; x  ^=005(4(^-4))  (2-83)  In order to establish the relation between these two parameters, we define: A = \S -S\  (2.84)  p  The variation of A with the bottom Biot number and the ratio h I h is investigated in r  B  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, S  p  is first calculated from Equation (2.59). The effective half-thickness,  L , is then estimated as:  L =L(\-S ) e  p  34  (2.86)  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 Bi and hj/hg ratios B  in the case of temperature ramp. Figure 2-16 shows the maximum error for any  hjjh  B  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 h are o f the same order, strategy #1 yields lower errors. Generally B  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 will 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 will 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]. HajiSheikh 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 o f 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 multipledimensions ([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 will 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 o f 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 o f solid with general convective boundaries as shown in Figure 2-20 is as follows:  37  dT  •  2  5T dt  0 < z < 2L  dz  (2.87)  dT 2  a  K dz  +H  2L <z<2(L +L )  2  dz  2  2  x  hB.l _ \ K  =  +  T  T  o  n  =  2  x  2  0  '  K  X  ^—^(T-TJ oz K  (2.88) on z = 2(L  x+  L) 2  2  (2.89)  T„=T +rt 0  The solution must also satisfy the following continuity equation at the interface of the two layers: dT  dz  5z  z=2i,  +  (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:  (2.91)  z-(2L +L ) x  2  and so the P D E and boundary equations transform: 8T _ X  dt ~ dT _  2  x  -\<C<+\  dC  2  (  2  dt ~  dT  M  ^ 2  aj d% ;  d?  (2.92) •+H  n  38  -1<£<+1  8T  X  +Bi {T -T )onC—\ B  i  a  (2.93) --Bi (T -T ) r  84  2  on 4 =  x  where we have  Bi. =  ^  k  (2.94)  1  Bi  T  K  2  Here T , T are the values of the temperature field T in the regions 1 and 2. Therefore the x  2  continuity equations can be rewritten as:  T\  =T\  L 8C '  ( X  x  84 '  ~ L  f=_1  2  2.5.1-1 The steady state (asymptotic) solution We now assume that the steady state solution has the following structure: T =T„ B (t)  = rt  T =T +B {4)  = rt + T  lss  2SS  +  K  l  2  +  T B (<4) 0+  0 +  Substituting the above into Equations (2.92) results in  39  l  B {4). 2  (2.97)  1  B2  -  ^  +*2$  +  fl  where  (2.98)  tf = 2  v  a j  (H -r). 2  2  Using the boundary conditions (Equations (2.93)), we can write (  + e - +Bi x  B  R +e =-Bi 2  2  r  1  ^  V 2 (  1  v  2  (2.99) ^  -~R +e +f 2  2  2  ;  While the continuity Equations (2.95) lead to the following:  ~R +e +f =-^R -e +f l  l  l  S (-R +e ) l  l  i  2  2  2  =  (2.100)  S (R e ) 2  2+  2  where  (2.101) ^2  Equations (2.99) and (2.101) can be rewritten as  40  f  1 ^  f  1 1+ Bij i  2  1 ^ SB  \  e +f =R 2  2  2  2  l  ^2  5/y  ST  J  (2.102)  1  e +f +e -f =-(R -R ) l  \  —+ —  1+Bi,Bj  2  (S )e -(S )e =S R +S R l  l  2  2  l  i  2  2  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  8T  2  {  K l  <r=+i  \L  2  ) 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 oo  -T.  (2.104)  i  Here q can be positive or negative and so is h . If h > 0 then it can be used as an i  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 Bi Bi S S (/?, B  T  x  (1 +1 / Bi )-R (\  2  B  +  2  \lBi )) T  S R (1 + Bi )(1 + 2Bi ) + S R ( l + Bi ) ( l + 2Bi ) ' X  X  B  T  2  2  T  (2.105)  B  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 i  +T  iSS  1  1  -*2  (2.106)  -M7V  T  2SS  27> '  Substituting into Equations (2.92) and (2.93) we obtain aAd TITr  (  2  -1<£<+1  dt  (2.107)  dT2Tr  2Tr  dt  \2  -1<£<+1.  J  L  • = +Bi .T @t B  dt  = -\  m  (2.108)  dT2 _ = -Bi .T @4 dt T  = +l.  2rr  The continuity relations are that T  =T  l  lTr  K  L  dT  x  dt  I  _  ~  |f=1  2  dT  L  2  42  (2.109)  2Tr  K  dt K -  1  The solution has the following form  r™=E%(£)exp(-4 0 2  ;=1  (2.110)  T =f E<Z ({)exp(-A t) 2  Tr2  j  i2  l  Substituting into Equation (2.107) results in:  Z„=0 V  «i  J  2  J  z" +  -1 (2.111)  Z ,=0  -!<£<+!  2  \  a  This system has the following solution:  z  w =cos(j7„(£-<5i )) (  (2.112)  Z . =<r,.cosfa ,(#-£ ,)) 2l  2  2  where  (2.113)  %i=A, L /a 2  2  2  2  Substituting Equation (2.112) into the boundary and continuity Equations (2.108), (2.109) we obtain  Vi,-sin fa,)- 5 / c o s f a „ )  0  cos  ?i,cosfa,) + 5 / s i n f a , )  0  sinfa,,)  0  -7 ,sinfa ,) + i?/ cosfa ,)  f l  det  f l  0  2  2  ^,cosfa,)  r  +  J  2  g/ sinfa,) r  fa,)  -cosfa .) 2l  sinfa ) 2/  .sin  fa,)'  cosfa,) -7/ , sin fa,) 2  -77,, c o s f a , ) (2.114)  43  0  The roots of (2.114) determine the set {A,}, while {£,,}, {S }  are calculated from the  2j  following: >7i,sin(j7„)- BI cos (77,,) B  tan(/7,A) = -  7i,cos(/7„) + 5 / s i n ( % ) B  tan(/7 ,^ ,) = 2  (2.115)  sin (tin ) + T (VI,) /7 cos(?7 ,) + 5/ ,sin(77 ,)  -Vzi  2  B  c  I  2  2/  o  s  2  7  On the other hand, we have: ^ _ cosfaQ.O-^))  (2.116)  cos(^ ,(-1.0-J ,)) " 2  2  The coefficients {£.) are calculated from the following as a result o f the orthogonality of the trigonometric functions: +i  +i  (pC \L \{T -T )cos(rj P  E.=  y  0  xss  {£-S ))d£  u  =1  + a, (pC \  u  L . \(T -T )cos(rj  P  2  -  (pC \ L .jcos x  2SS  (<f-S ))d£  2l  2i  =L 2  P  0  (  Vu  (<T - S )) dt + a  2  u  (C P  ) L .jcos  2  P  2  -i  2  (rj (<f - S )) d^ 2l  2>  -i  (2.117) 2.5.2  Applied Hold Temperature  This is the case where the fluid temperature jumps to a constant temperature T *T . In A  this case the steady state temperature is obviously T = T xss  2SS  0  = T throughout the A  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 will 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 Bi ,Bi ,R R ,S ,S B  T  u  2  l  2  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 is:  ri =  h  h  (2.118)  where Bi  n  is the equivalent Biot number for the upper surface of the first layer. The  location of the adiabatic line is  Bi.ri 2+  BIT l  And so we have  45  Bi,B +  (2.119)  I )  max ( a / , ,  R =R (\ e  \Q  (2.120)  2  r  +  Fo ^/(l e  \C \f.  +  a  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: h  = -h  B2  Bi  =-  B2  hL  (2-121)  2  K  2  where Bi  B2  is the equivalent Biot number for the bottom surface of the second layer. The  location of the adiabatic line in this layer become: J B t  2+  1_ \  B I b  Bi,  +-  \  •  (2-122)  Bi  B2  And so we have: Bi =max(Bi ,Bi )x(] e  B2  R =R x(l e  2  Fo =^/(l L e  + \£ \)  T  a  + \{ \f  (2.123)  a  +  \{ \y. a  2  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 1mm10cm and the autoclave convective heat transfer coefficient varies between 5 WI m K 2  and 200 W  lm K. 2  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 (T ) involved versus the autoclave err  temperature (T  )  auloc!aw  and the temperature lag (T ) for the cases where the (steady state) lag  adiabatic line occurs in the second layer (h < 0). The temperature lag (T ) considered lag  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 (T (t)) is the difference between the err  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 , , ) in these diagrams is recorded at the time when the e  ratio \T {t)JT err  lag  (?)| reaches a maximum, where T  At the same moment, T  lag  lag  and T  autocalve  (t) is the accurate temperature lag.  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 (Ar  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 AT  . These  hold  parameters are recorded in the diagrams for different cases of interest at the moment when \T (t)/T (?)| reaches a maximum. Both diagrams show relative errors of less err  than 10% of AT  lag  hold  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 + 2 5 / , )  -h~L  This equation shows that the composite could be modeled with a Bi  B2  K  in the  2  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 sections. 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 K  2  and a  2  will denote the in-plane  direction conductivity and diffusivity. The other non-dimensional parameters defined here are  i-T  (2.125)  BL = •  These parameters are called the non-dimensional in-plane coordinate and the in-plane Biot number in this section.  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\d T  ( d  2  dt  -+  d?  T  KL j d? 2  = -Bi(T-T )  on £ = +\  = +Bi{T-Tj)  on £ = -1  a  8T_  2  dT_ = +Bi {T-T„)  i  on | = 0  2  d  (2.126)  8T_= 0 on ^ = oo dt; The autoclave temperature T  x  and the initial condition T are assumed to be constant. It 0  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„  = u((Z,t)-u  2  T -T 0  (2.127)  (l,t)  m  N o w we use the following transformation of the variable: T =  T„-(T„-T )-u(£,t)-u (l,t) 0  where u and u satisfy the following PDEs: 2  51  2  (2.128)  du  f  a \ du 2  ~dt'  dt  2  du — = -Biu on C = +i dt du _. — = +Biu on c=—\ dt «(^,0) = 1 „ .  /-  1  (2.129)  and 9  9m,  9/  m.  9£  9m. *2  _  +Bi u 2  2  2  on £ = 0 (2.130)  9w. - - 0 on 4 ~  00  M (^,0) = 1 2  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 Fo)  2 sin (A,.)  2  i=l,3,5,..  X + sin (/l,) cos (/I,)  (2.131)  i  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  u (l,t) = erf  erfc  + exp Bit  2  \  2 J  K  \  K  ZNFo  + Bi Fo 2  \  2 J J  K  2 J  (2.132) In the case of an isotropic material (K = K ), this simplifies further to 2  u (l,t) = erf  tl4F~o  + exp{^Bit)erfc  2  + Bi Fo 2  (2.133)  J)  So, in mathematical terms, the edge effect question would be whether or not u is a valid estimation o f the product uu . The error involved is the difference between the 2  estimation u and the product uu  2  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  Bi,—~,err  4  •• minl^  max (u-uuA<err\  .  (2.134)  J  v  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  err = 0.01,0.02,0.05 and K/K  2  K/K  2  varies as a function of Bi for  = 0.01. In Figure 2-27 to Figure 2-33 the value of  would be 0.05, 0.1, 0.2, 0.5, 1.0, 10, and 100. Figure 2-34 shows the maximum  value of t  ed  among all values of Biot number defined as:  ec/ge-max  = max  53  {tedge}  (2.135)  The values of K/K which are less than one are of special interest in practice as in-plane 2  conductivity of composites is higher. The value of this parameter depends on the layup we are using. Figure 2-35 demonstrates ^ 2.8.2  e d  for the case of Dirichlet boundary (Bi = oo).  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: dT  _(a\d T  8t  \J})dC  dT 2  2  •+  — = -Bi(T-Tj)  d<;  d?  +H  on £ = +l  '  v  ?L = Bi{T-T )  on C = -l  ^ - = +Bi (T-T )  on | = 0  +  a  2  dT_  x  (2.136)  = 0 on £ = 0  where T=T +rt 0  (2.137)  ,  and the initial condition is T and r is the rate of autoclave heating. 0  Now defining T = Rg + Ht + T  0  we can write  54  ,  (2.138)  dg _ d g  1  2  +  dg 2  dC ' f K_\ dl  dFo  2  K  K  2,  ^ - = -Bi(g + Fo) on C = +l ^  = +Bi(g + Fo) on ^ = -1  ^  = +Bi (g + Fo) on £ = 0  (2.139)  2  5<T = 0 on £ = < g ( F o = 0) = 0 .  Now using the Duhamel's Theorem, we form the following auxiliary problem: dL  5 E 2  dFo  1 5 E +' { K\ dig 2  2  \  51  2 J  K  = - 5 / ( S + l ) on £ = +1  — = +5/(2 + 1) on £ = -1  (2.140)  = +fl/ (l + l ) @ | = 0 2  —  5|  = 0 @ # = oo  -  E ( F o = 0) = 0  Comparing Equation (2.140) with Equations (2.129)-(2.130), one can infer that Z = u u -1 . 2  (2.141)  or in the case o f 5 / = 0.0 2  E =«-l 55  (2.142)  And hence the error is identical as both cases are normalized with respect to a maximum absolute value of one. Namely E - E = w - ww  2  (2.143)  The solution of Equation (2.139) could be calculated using the auxiliary Equation (2.140) by  (2.144)  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/K , 2  on the other hand, should be evaluated for the layer o f 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 o f 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 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 T =T + rt where the heating rate of the autoclave is r. b  b0  The autoclave temperature changes linearly with the rate r as well and is described as T„=T +rt. 0  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 T -T m  and heat inflow (outflow)  b  ai 2  A / C - ^ - V = C/zA I 2  l = T -T b  -k  In the above equation, I = T-T  ai  dm  0  on nr = 0.0  (2.145)  hM on £7 = 1.0 .  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  1+  i = (r»-r) cosh  ^^tanh \ h A  sinh  f  — — +tanh h A  K A~  KA~  l  (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  AK Q  =  =  h C  {T -T)AK b0  x  KA  1+  ££tanh V h A  ( KC — — + tanh V h A  h C I——A > 2.6 and therefore tanh \KA~  K A~  (2.147)  K A'  ]  For the cases where  = 1.0, we have  (2.148)  For the case of n  fm  fins in unit length we have  Qtotal ~  (2.149)  (^fcO  If we are investigating unit length of the fin then A -1 x 8 and C = 2 (1 + S). Therefore,  _  n  "  K  QtotalA  (T -T ) b  vh j  m  58  1 V  +  I -1 A  8_j  (2.150)  Notice that h > h i f — > — — which is the case for tool material. Figure 2-37 h l+S eq  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 o f interest in tool-part set but lacks the generality of the singlelayer 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 = — K  (2.151)  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 ^ - ^ ^ R  1 0 =0.5+—. Bi  (2.152)  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), S  p  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 8= P  1_  \ \ • 2+— + — Bi Bi B  B  r  t  (2.153)  B  The top and bottom Biot numbers are defined based on the corresponding h values, i.e.:  60  B  Bi -^  K h,.-L  T  (2.154)  K  In the strategy #1, an equivalent Biot number is calculated from the following equation  „ 1 1 2+ + — Bi Bi j 10 4 4  Bi = -  T  Bij  Bi  B  Bi -Big  (2.155)  B  -+  1 +Bi .Bi T  T  B  j  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(h ,hj.)-L B  e  (2.156)  .  (2.157)  K  where L.=L+S  p  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 Bi ,Bi ,R R ,S ,S B  T  v  2  x  2  according to: 61  (2.158)  hjL  2  ~K7  Lz^l  r  2  V«i J  (2.159)  '' L ^ (** -r) R> = 2  2  V  2  a  J  (2.160)  2) Calculate the interface heat transfer coefficient h from:  h =h  ss  =  Bi Bi S,S B  T  2  (tf, (1 +1 / Bi )-R (\ B  + \f Bi )) T  5,/?, (1 + Bi ) (1 + 2Bi ) + S R ( l + 5 / ) ( l + 2Bi ) B  3) If  2  T  2  2  r  (2.161)  B  > 0 then /? is the equivalent heat transfer coefficient for the tool, that is:  (2.162)  Bi,,.. = —-  where Bi  n  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 Sa  1  __Bi  Bi  n  B  * 1 1 2+— + — Bi Bi Tl  62  B  (2.163)  And so we have  Bi =m x(Bi ,Bi )-(l e  a  R =R (\ e  r  F o  e  B  + \t \)  Ti  +  = f i / ( l  a  (2.164)  \Q  2  \C \)  +  a  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 Bi  B1  is the equivalent Biot number for the bottom surface of the second layer. The  location of the adiabatic line in this layer would be  Bi  Bi  T  t  (2.166)  ^ 1 1 2+ + — 'B2 Bij Bi l  And thus we have Bi =max(Bi ,Bi )x(\ e  B2  T  a  (2.167)  R =R x{l+\l\f e  + \g- \)  2  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 o f parameters that illustrate the effect o f the second dimension:  H (2.168)  M  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 value of g  ed  2  ). Figure 2-34 shows the maximum  among all values of Biot number and Figure 2-35 demonstrates g  ed  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 n  fm  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 o f 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 o f 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 ) 8000 1580 2710 J  INVAR COMPOSITE ALUMINUM STEEL  C (J/kgK)  k (W/mK)  515 870 896 465  11.0 0.69 167 51.9  p  7860  66  z  a (m 1 s) l  z  2.67 0.50 68.9 14.2  x lO' x 10" x 10' x 10  6  6  6  -6  Time Figure 2-1 - Typical autoclave process cycle.  Fluid at T^hj. + L . + 1  i  - 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  1.E+03  1.E+04  Fo 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+02  1.E+01  1.E+03  1.E+04  Fo 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 h /h T  10000  100000  B  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 Bi and hj. I h . B  B  Figure 2-8 - The global maximum errors extracted from curves in Figure 2-7 plotted as a function of Bi . B  70  0.09  t  1  10  100  1000 h /h T  10000  100000  B  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 Bi and hj I h . B  B  Figure 2-10 - The global maximum of curves in Figure 2-9 as a function of Bi  t  71  h /h T  B  Figure 2-11 - A as a function of hj. I h for 0.001 < Bi < 0.5. B  .  B  •  10.5 |  1  Bis\  j  s ho l 1  1  1  i1  10  100  1000 h /h T  Figure 2-12 - A as a function of hjlh  B  -inn .i ,  B  for 0.5 < Bi < 100. B  72  w  1  10000  1  100000  0.06  T  h /h T  B  Figure 2-14 - Maximum error of strategy #2 in transient phase for autoclave temperature ramp; the error is presented in terms of h, I h for B  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 h for B  0.5<2?i <100. B  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 h  L  for 0.001 < Bi < 0.5. B  0.05 -i  1  10  100  1000 h /h T  10000  100000  B  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 h I h r  for 0.5 < Bi  B  <100.  75  t  0.05 i  0.04  0.03  0.02  4  0.01  0.001  10  0.01  100  Figure 2-19 - The global maximum of curves in Figure 2-17 and Figure 2-18 versus Bi • B  autoclave . L  K , oc ,H 2  2  e  T ,h m  1  2  T  )  /  layer II: part layer I: tool  2  ± ,  2 L j ±  T Kj, aj  2 L  z  y '  autoclave:  T 00  '  h 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 andfor cases when the adiabatic line falls in the second layer (composite part).  + 140 ~ -- 120 *•£  o  -- 100  • Temperature Lag • Autoclave T  -- 80  I 3 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.02  0.04  0.06 WAT,,  0.08  0.1  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 T ,h  .  x  +L,+1  1  Fluid at  1  -L,-l  F l u i d a t T„ h 9  Figure 2-25 - Schematic of the model used to study the edge effect.  -  /  1.E-02  err=0.0  1.E-01  /tH)2  1.E+00  /tH)5  1i  I  1.E+01  1.E+02  1  1.E+03  Bi  Figure 2-26 - 4  edge  K/K  2  os a function of Biot number for different error tolerances and for =0.01.  79  r  r  r  r  r  ferr-QX  >  ^002  7  1.E-03  1.E-02  1• — ^  1.E-01  1  1  1.E+01  1.E+02  1  1.E+00  1.E+03  Bi  Figure 2-27— <^  edge  K/K  = 0.05.  2  25 ~|  as a function of Biot number for different error tolerances andfor  1  1  i  i  i  i  1.E-03  1.E-02  1.E-01  1.E+00  1.E+01  1.E+02  20 -  o> •a o  15 -  Ol  10 501.E-04  Bi  Figure 2-28 - ^  as a function of Biot number for different error tolerances andfor  edge  K/K  2  =0.1.  80  25  T  -  —  r  i  r  I—  r  r  r  r  ^ i  1  r  i  _ ,—  •D O  /err=0.01  10  0 1.E-04  1.E-03  1.E-02  1.E-01  1.E+00  1  1.E+01  1.E+02  Bi  Figure 2-29 - 4  as a function of Biot number for different error tolerances andfor  edg(!  K/K  = 0.2.  2  25  -r  r  i —  r  r  r  ~ ~r  i  20 15 4  r  r  ~ ~ i  /ferr=0.01  10 /  /om  * 0.05  1.E-04  * 1.E-03  — i — 1.E-02  1  1.E-01  1.E+00  1.E+01  1.E+02  Bi  Figure 2-30 - 4  edge  K/K  2  as a function of Biot number for different error tolerances and for = 0.5.  81  Figure 2-31 - g  as a function of Biot number for different error tolerances and for  edge  K/K =\.0. 2  <err=OM /0.1 32  \ 1  1.E-05  1.E-04  1.E-03  1.E-02  —11  1.E-01  1  1.E+00  1.E+01  Bi  Figure 2-32 - g~  edge  as a function of Biot number for different error tolerances andfor  K/K =\0. 2  82  Figure 2-33 - 4  edge  K/K  as a function of Biot number for different error tolerances and for =100.  2  25 -t  0 -I 0.01  1  \  i  1  0.1  1  10  100  K/K  Figure 2-34 - g~ _ eilge  m!0i  2  as a function of K/K for different error tolerances. Markers 2  show the maxima from the previous diagrams and trend lines are drawn through them.  83  Figure 2-35 - ^  edge  as a function of K/ K for different error tolerances for the case of 2  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  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  Nomenclature  lUI^  energy norm for volume V  «(.,.)  typical bilinear form  A ,B ,FZ m  m  m  Parameters in spatial integrations for K  Bi  Biot number  c  degree of cure  C  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  D  thickness of the shell element at each surface integration point  p  KGP  ed  err  ge  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(£)  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,  over  typical smooth function 86  GF  geometric shape function  h  coefficient of convective heat transfer  h  element size in x-direction  h  element size in y-direction  H  exothermic heat generation per unit volume  x  H  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  R  k, k  conductivity on both side of interface  KIP  integration point index  L (.)  typical linear functional  L  half thickness of a slab  t  2  L  L  lower bound for error  Ly  upper bound for error  L (Q)  space of functions over Q equipped with L norm  M  number o f Fourier terms  M  surface shape functions matrix  n  number of trigonometric terms  N  number of layers  N  interpolation shape function  NG  surface geometric shape functions  n  number of Gaussian points in surface  2  T  GP  2  n _g  number of Gaussian points in £ direction  GP-ri  number of Gaussian points in 77 direction  n  number of interpolating shape functions  GP  n  IF  87  n  number of integration points  n  number of trigonometric terms  O  order of the error  q,q  heat flux, heat flux vector  Q  exothermal hear source per unit volume  lp  T  Q  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  s  S  surface with Dirichlet boundary condition  S  surface with Neumann boundary condition  S  Jacobian stiffness matrix for nonlinear procedures  t  time  T  t  estimation of the time range  T  temperature distribution function  mnge  T  estimation of the temperature range  T  fluid (autoclave) temperature  TG  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  range  m  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  v  resin volume ratio  p  density  R  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. A s 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, multiscale 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 ( U E L ) 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 = -kVr  (3.2)  where k is the conductivity tensor defined by  k=  yy  k  k  k  zx  (3.3)  yz  k  k  zy  2  If conductivity is given in a local coordinate system, we can calculate the global matrix from the following equation  k = Ak A t  where k  L  r  (3.4)  ,  is the conductivity tensor in local coordinates and A is the transformation  matrix that is defined by X.x  X.y  X.z  A = Y.x  Y.y  Y.z  Z.x  Z.y  Z.z  Here ( X , 7 , Z ) is the global coordinate system and (x,y,z)  (3.5)  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  = Q + JQdV s  93  (3.6)  where U is the internal energy per unit volume and Q is the heat source per unit volume. Also Q represents heat flow through the boundary and is defined by: s  Q,=-\qadA A  = -\v.qdV V  ,  (3.7)  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 dt  = -V.q+g  .  (3.9)  Then substituting q from Equation (3.2) results in the following equation: dU dt  V(kVr)+g  .  (3.10)  The rate of change of U with temperature is given by dU dT where p is the density and C  p  = pC  p  ,  (3.11)  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  a dt  - . V T  +  Q , k  (3.14)  where a = k/pC  (3.15)  P  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) = T (x,y,z) 0  (3.16)  Boundary conditions, on the other hand have one of the following formats: T(x ,y ,z ,t) b  b  b  = T (t) b  - * ( ? - ) » = ft (0 on  (3.17)  (3-18)  The above equations are called the boundary conditions of the first kind (Dirichlet) and second kind (Neumann). q  B  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=Wr /)-r.) i )  on  (3.19)  )  where h is the convective heat transfer coefficient of the surface and T is the x  temperature of the ambient fluid. This kind of boundary condition is called homogeneous if  T =0. x  Other possible boundary conditions are radiation boundary condition:  -k(~) =as[T\r ,t)-T ]  (3.20)  4  b  b  e  and interface boundary condition as defined below:  -k-(^-\=-k (^f-) on  (3.21)  +  b  on  where 8 is the emissivity of the surface and a is the Stefan-Boltzman constant. k and k y  2  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  pX  -zr=- -q+v p H — at v  R  ot  R  R  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. (pC )  average (smeared value) of pC  P  p  for the composite while p  R  denotes the  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 nonconventional 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 will 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)  98  (3.23)  pC T, = k T P  x  + kT  + kT  + H(x, y, z, t)  (3.24)  T(x,y,+L,t)  = T (x,y,t)  (3.25)  T(x,y,-L,t)  = T (x,y,t)  (3.26)  T  B  T(x,y,z,0)  = T (x,y, z)  (3.27)  i  The P D E above can be rewritten in the following format: T, = aT  where a , a x  y  and a  + aT  + aT  + (1 / pC )H(x, y, z, t) P  (3.28)  are diffusivities in the corresponding directions. In order to  z  change the boundary conditions into homogeneous boundary conditions we apply the following transformation in the variable T  T(x,y,z,t)  ,i  U  =  s T-r ~T" T T-p = (^-^ + ^ 2 Q  + y ,yy a  ) + 0  + z ,zz  U  U  U  ^2 x  Bxx  u (x ,y ,-L,t)  i(x,y,z,0)  Z  B  L  ) + u(x,y,z,t)  JV, > 0  Pp  1  (3.29)  C  Z  (3.30)  2L'  }__]_ z_. )(----) "2 2 V  <T -a T BJ  +  T  =  =u [x ,y ,+L,t) = 0  T (x,y,z)-  fT.,.+T  i  T ' fl | x  B  (3.31)  T -T  T *B r  2  0  L  (3.32)  In order to find the proper shape functions through thickness, the following eigenvalue problem is considered  99  ,i  x  (3.33)  + oc u +a u  u. — a u ,xx  u (jc ,y ,-L,t)  y  z ,zz  ,yy  =u (x ,y ,+L,t) = 0  (3.34)  N o w consider the following solution based on separation of variables: u(x,y,z,t)=W  (3.35)  (x,y,t)Z(z)  Substituting into Equation(3.33) results in  W  , ,  Z  =  a  *  W  , * *  + y  Z  a  Z(-L)  W  , y y  +  Z  a  ,  W  Z  "  (3.36)  (3.37)  = Z(L) = 0  Rearranging the terms in the above differential equation leads to  w  w  z  w  y  (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  .nz sin 2Ly J  K  100  j is odd (3.41) j is even  And the general solution for u can be written in the following series format:  u=^W \x,y,t)Z \z) {J  (3.42)  (J  7=1  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)  = a WW+a WV>-a x  y  J (1 / pC )H(x,y,z, h \x,y,t)  =-  +h (x  2  t)Z  U)  (j)  ,y ,t)  (3.43)  (z)dz  (i)  p  (j  (^-) W  z  +L  \z \z)Z \z)dz u  \ \ { T  T  t  u  - a J ^ - a ^ ^ +^  1 1 z, iT^-aJ^-a T ^-^)jZ (z)d;  +  {j)  y  B  _ jZ \z)Z \z)dz {j  (j  (3.44) The above P D E is to be solved using the following initial conditions: TT T  B  W°\x,y,0)  T Z  T  +  B  T  2  =  +L  Z°\z)dz d Jt=0,  (3.45)  \z \z)Z (z)dz (i  (i)  The solution of the above problem can be performed using Finite Element Method (FEM) discretization for each W , i.e. }  W^{x,y,t)  =  ±N (x,y)W^{t) n  101  (3.46)  where N  is the number of interpolation shape functions in each element and [N } N is IF  JF  n  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 will turn out to be 2  jwWJ dS + a \w W ^dS J)  [  x  s  x  s  = jwh (x,y,t)dS U)  + a jw W ^dS  +a  (  y  y  z  (3.47)  s  ,  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  (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 8noded 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 = £GF,(£,,7,C)x, ;=1 "GF  ^ = £GF,(#,J7,0y  (  (3.52)  ;=1 "OF  z = £GF,.(^/7,C)z,  Here n  GF  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. N o w we define the Jacobian J as follows:  dx/ /dt  dx/ /drj  dx/ M  J = dy/  /dg-  dy/ /dn  dy/  dz/ /dt;  dz/ /dn  dz/ /dZ  /d£  (3.53)  Substituting Equation (3.52) into Equation (3.53) leads to the following equation for the Jacobian  (3.54)  y, - { , z. GF  1=1  V  1 1J  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 = J d l dV = det(3)dU  104  (3.55)  where dL = [dx,dy,dz}  and dV = dx.dy.dz are the differential vector and volume in the  T  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, (J dlj) = (jydA,) dl = (j],dA,) dl = (J d A).dl r  v  }  }  det (J) (da.dl) = (det (J) da) . d l .  The second equation above holds as det(j) is a scalar. Therefore, further matrix manipulation on Equation (3.56) leads to ( j d A ) . d l = (det(j)da).dl for arbitrary dl , T  (3.58)  And therefore J dA = det(j)da T  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}  =±dQ{0  T  0  l} ,  (3.60)  r  which, using Equation (3.59) leads to  dA = ± ( d e t ( j ) J Q ) j - { 0  0  T  l} .  (3.61)  r  Further matrix manipulation on Equation (3.61) leads to  dA = ± ( d e t ( j ) J Q ) { l N V J  31  INVJ  INVJ } , r  32  33  (3.62)  and taking the absolute value leads to  dS = |dA| = ( d e t ( j ) 7 l N V J  2 3 1  +INVJ  +INVJ  2 32  JdQ. ,  2 33  (3.63)  where INVJ = J" and dQ. is a typical element of surface area on top or bottom surface 1  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  +INVJ  2 3 1  + INVJ  2 3 2  ,  2 3 3  (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  smf (J) = det (J) ^/iNVJn + I N V J 2  4  f (J) = det (J) / l N V J  sm  v  A  2 2 1  +INVJ  2 12  2 22  +INVJ  2 13  +INVJ 23  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  (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 n  IF  is the number of interpolation functions. These interpolation  functions are defined within an element as follows:  m N^M^f^),  n=1  /= 4+m  n = n +2  ,  T  7 + n + (m-\)-n  n=  T  (3.67)  2,3,...(n +l) r  where m = 1,2,3,4 refers to the element edges and n is the number o f trigonometric T  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. M  m  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 = l / 4 ( l + £)(l-77) 2  v  M  ;  A / = l / 4 ( l + £ ) ( l + 7) 3  M = l / 4 ( l - < r ) ( l + /7), 4  Here zf (<^) is the through thickness shape function and for the case of assumed n  Dirichlet boundaries, it is defined as follows:  107  (3.68)  2/j(0 = l / 2 ( l - 0  ^  r + 2  ( 0 = l/2(l + 0  z/;(C) = 2 ( O  / = 2,3,...,/v+l.  (0  Z  ( ( )  (3.69)  ( 0 is defined as (Equation (3.41))  7C cos(i—£)  i is odd  2  sin(/'—C) 2  .  (3.70)  i is even  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 v  dt  dV = -$wV.qdV+  jwQdV  v  ,  (3.71)  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  = -\wq dS+\wQdV n  ,  (3.72)  where S is the boundary surface and q is the normal component of the flux to the n  boundary. The boundary surface could be broken down into S = S +S T  108  q  ,  (3.73)  where S  T  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  fw  r\  j r ,  dw  dw  dw  (3.74)  dV  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  (3.75) w = Nj  j = l,2,...nIF  Substituting Equation (3.75) into (3.2) leads to  (3.76)  q = -£(kVJV,)i; 1=1  Substituting (3.75) and (3.76) into (3.72) yields the following equation: dU \Nj—dV \v  0*  ^ J  + £  j ( k W , ) . V A ^ T^-JN^dS+JNjQdV  1=1  )  \V  ,  (3.77)  V  S  a  which results in .  "IF  Y  j +  Y K.T =-\N q dS j  i  ]  n  109  +H  ]  (3.78)  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:  q  (3.80)  =-<io+  hT  n  Therefore the right hand side of Equation (3.78) can be rewritten as  "IF  -\N „dS=  "IF  \N dS-X  jq  (3.81)  jhN^jdS  jqo  Finally, Equation (3.78) leads to the following equation in matrix format, which should be solved using a finite difference solver in time (3.82)  Y+(K+K)T = H + H .  Y is evaluated for the case where U is only a function of temperature, as follows:  dT dT 1  fit  dt  J  i  1  J  Pit dt  • J  J  \ AT dT  KdT,  dt  dV ,  (3.83)  which results in  rj=l\l Mp p)w}t N  C  l  = fc t l  y l  (3.84)  The following are the definitions of the above tensors in the physical space and also in the parent space  110  n  v  K = l(kVN,).VNjdV= v  v  K = jhN Nj s, v  t  J ( k J - V A ^ , ) . ( j V i V ) d e t ( j ) dU n dS= \h N,Nj smf(3) dQ a,  Hj = JNjQdV v  H =\q N dS= s, J  0  r  -r  y  = JNjQdet(j)dU n  (3.85)  .  \q N mf(J)dn  j  0  jS  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 timeconsuming if 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  f(C)  .  (3.86)  7T 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 will further discuss the cases of interest later on in this section.  Ill  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  jp  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  /  i = 2,3,...,n +\  2  (3.87)  IP  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/(Os(0^ = Z J/(Os(£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 F 's are dependent on the number of integration points used in the analysis and can t  be evaluated a priori. A s a matter of fact F and other parameters that are not dependent f  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 will be used in the evaluation o f 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 thefirstand second integrals of (3.85). The rest will follow suit. For thefirstequation in (3.85) we have C, = C  J(/>C,)tf,tf det(J)dri  fcpC )N N dV= p  l  J  y  n  v  v = T  J 7 (pC )(M^M )^{3)(zf^zf )d4dr,d^ p  H  *tf.</.0  ££ "pp  Cjj ~  KGP=\  (3.91)  H  AO  "IP  8(.4KGP' 1KGP' ?K1P) <  T  KGP  WG  K1P  '  F  KIP=\  where F is determined by integrating / as given by Equation (3.90). n  GP  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: GP ~ GP-i- GP-  H  N  (3 92)  N  N  wG  KGP  =  wG _ wG _ KGP  114  r  KGP  n  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(kJ" ViV ).(j- ViV )det(J) J n n r  r  /  v  ;  C, =+1 =+l cf=+l n  *  I {J  k J  <T=-I <7=-l f = - l  r  det(j)^^C  \M Jz m  ni  M fz'  M fz'  m  m]  nj  J  f = + 1 ^=+1 <f=+l  kJ  X* X, 0  \ r  0  J 0  0  M  M„  i"j,y  0  0  fa  \K  det(j) d^dnd^  i = \,2,3,...n,  F  m = \,2,3,A (3.93) N o w we define the following tensors: 0 " A =k.T m  (3.94)  0  r  0  _  M m  0 " B  m  =J  0  r  0  FZ"  m_  "UJ  Equation (3.93) can now be written in a more compact form as  115  (3.95)  (3.96)  ,  f=+17=+l f = + 1  </ = I  I  K  j (A <FZ"').(B 'FZ"')det(J)J^;7^ m  m  <r=-l,7=-l<?=-l  K  #  7 ( u  f j  'l=l '2=1 '3=1  f=-1  A  l/=-l  f=-1  F  Z  2  >  ")( ^ F Z j ) d e t ( J ) ^ ^ B  v  '*  g(i.i,o  v  (3.97)  '  no  Therefore the element stiffness matrix is computed as follows:  K,  =iiiTTf(K, ;,=1 / = i / = i ^ _, 2  3  3  2  =  2  "cp  , ^ _, > =  ";p  )det(J)(FZ- F Z J ) < W < T v  g(.tn,o  ,<  v  no  ,  (3.98)  /,=1 t =\ / =1 ATG/>=1 K1P=\ 2  3  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 O D E ' s in time should be integrated to obtain the transient solution. The 0 method is used to evaluate the nonanalytical 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 O D E ' s . A s 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 P D E ' 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 nonstiff 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 exothermtime step used in evaluating the degree of cure c and hence the exothermic heat generation 3.6.1  H.  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  where 117  l + 4 l  -T ) +(T l  l + A /  -Y,)) = 0 .  (3.100)  (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 o f c (Equation (3.22)and Equation (3.85)) as Y , = \Nj(v p H c)tet{S)dTl n R  R  .  R  (3.102)  H , on the other hand, is calculated from q which is part of the user input. Based on Q  equations (3.85) and (3.101) we have <P,=  \QoNjdS BIN  (3.103)  Qo = \ o • a  d t  The 6 method is used to evaluate the last integral. 0 method is based on the following approximation:  ]f(t)dt  = {t -t ){9f(t ) 2  x  2  + (\-e)f(t )) x  where 6 is an appropriate number in the range [0,1]. Evaluating the integrals in Equation (3.100) by 0 method leads to  118  ,  (3.104)  J-PM '  - Y,) + [<>(K + K)|  T  (+A(  + (1 - 9)(K + K)[ T,)  A  3.6.2  (3.105)  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, , - *P,) we need to +A  calculate (c, , - c , ) at each integration point +A  ((3.102)) as  (^ ,-^)=iN(v^A)(c, ,-c,)det(J)c/n + A  + A  .  (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' or exotherm-time increment, which is used in contrast with diffusionp  time increments, At, in this section. "Time increment" in this document means Diffusiontime increment. If c'  and  p 0  are the degree of cure and temperature at the integration  point under consideration at t' = t, then we evaluate c in the next exotherm - time p  0  increments using the explicit integration as follows:  cL =C f(CX - )^' , P  P  l+  ™ = U,3,..M ,  P  l  (3.107)  where the end o f final exotherm -time increment must be identical with the end of the diffusion-time increment, i.e. t' =t + At and m = \,...M p  are the subinterval notions of  the time step At. In order to obtain T _ we interpolate between T ' and T/ 1P  p  X  t  119  p Al  to get  •IP  T '+ should be extrapolated from T t  (3.108)  history for the first iteration. The following  Al  extrapolation method is developed by the author and is best suited for the problems of interest in this study:  C=C,+[v^ ,/(pCp) ]/K-i,C.)A/  1  e  J  P  ,  m = \,2,\.M  .  (3.109)  For the second and higher iterations, the result from previous iterations is used as (3.110)  Equation (3:109) is based on neglecting the diffusion term within one diffusion-time increment in heat transfer Equation ((3.22)), i.e.  (p P) ^7 C  C  =  RPR R-Z  V  H  (3.111)  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.  120  Therefore, in order to solve this equation for T, , a nonlinear procedure should be +A  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 NewtonRaphson 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 , , ) = 0 as follows: +A  (3.112) l+M,r+l  where T  ( + A I  T<+A(,r)  , is the r'th approximation for T  F(T .) = - ( Y ( +  ( + A (  - Y , ) (^(K +  r  , . In our case we have  ( + A  +  F( I/+A<,r) •>  K)[  {  Different components of the Jacobian matrix,  121  +A(  dF  T , I+Al  + (1-0)(K  + K)[T,) (3.113)  are calculated as follows:  dY,  dU  dT.  dT  i v  J  l+Al  = [N C iP  NjdV (3.114)  N.dV = C  P  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,,  K + R  )L -) T  =( , K  dk  + R  .)L  +  dV +  dT dh_  (3.115)  (  in The summation £  k  dS  dT  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  In order to determine  dc ~dT  dk_dkdc_  dk  ~dT~~dcdT  8T  (3.116)  we should solve the following equation in time:  d_(dc_) _  (dc\  dt \ d T j \ d T )  122  T  '  (3.117)  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  l(N,v P H R  V  ^  A  ,  J  '\ t+M/j l  R  c|,  R  + A  dV  ,)dV=iN Njv p H J £ t  V  R  R  R  V\  U  1  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  J&/fys]=  (+A  .  , ) \BIN  (3.119)  y  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|  approximation to  dc_ dT  based on a perturbation to T\  . A secant  i+A/  would then be equal to:  ll+Al  dc_  M  dT  ~ \A) C  U+Al  123  (3.120)  AT \  is the amount of perturbation in temperature at the end of the increment, c  P  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,=ic +*(K s  S  +  K )4j s  dV  N,N v p H J  At,  R  R  R  t+M  (3.121)  J  Therefore, the force vector would be  F = - C ( T „ , - T , ) + ( # ( K + K ) T „ , + ( l - S ) ( K + K)T,) s  (3.122) - i (  T  « -  ,  ,  ' ) - i ( ' - - * ' ) -  As a result, the equation for nonlinear iteration then becomes  S  ( r Ar,, l T  +  3.7.2  +  -  T  (  +  A r , , ) -  -  (3.123)  F  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  J_-  ic(T ,-T,) ^(K At ( t  +  +  K)LT, , (l-^)(K +  +  +  K|T,)  (3.124)  where based on the 8 method we have: (3.125) 124  Now regroup terms in Equation (3.124) to get  i^(^n.JMi 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  iC »( +  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, ,,o=T, +A  .  (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 =(ic-|l- )(K Kl)„ (  (  p  i  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  2  9  )  ^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 =>M ' ) ' v v  (3  -  131)  where v is the distribution function and «(.,.) is the bilinear form in the following type of weak form equation: a(w,v) = Z ( v ) f o r a l l v  (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  ( ' ) =Z v  v  j(  k W  WjdV \v .= v .K.v . 1  ,)  iVj  (3.133)  i=l \ y  Therefore, the norm would be  ||vf = V v X . v r  .  (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. A s 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  Ih^l^m^d^l'l^-il) where ||eT  eafe(;  >  (3.136)  || is an error norm for the vertex and |r„ | and | ^ _ i |  are the absolute values  r  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 L and L as non-dimensional lower bound L  u  and upper bound for the error to avoid adding and dropping of terms, i.e.  Z  £  .  A  ^ Z ^ L < | ^ ^ | | < , . A ^ ^ z  ^ range  u  ^ range  128  ,  (3.137)  where At is the size of the time step and T  range  and t  range  are the estimates of temperature  range and time range of the whole analysis. The values of L and L L  v  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  e4fee  | | 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 will 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  equal to [D  KGP  KIP — 1 >D n  KGP  lp  KIP ] where D n  is the thickness of the corresponding  KGP  [P  Gaussian integration point of the element, n  IP  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 o f 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) R E A D 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  4) IF \err \ edge  n  T 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.  5) E L S E I F \err \ edge  T is more than L -At-^^v  (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 L O O P . 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  j  H  = j o j smf{i) dQ= a  £  N  wG  KGP  • q • AT (% 0  ,n  KGP  KGP  ) • smf (J (g , rj , ± l ) ) KGP  KGP  KGP=\  o  (3.138) and, K =  lhN N smf(J)dCl  v  l  J  (3.139) = Z  w  G  « 3 / >  • * •  N, {^cr, VKGP ) • j {ZKGP > VKGP )' f N  sm  ( (ZKGP » J  »  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 K  tj  (Equation (3.85)). 7) E N D L O O P . 8) END L O O P .  134  and Hj  ,  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 =  =  Z  \q N mf(3)da 0  jS  KGP-{ • % •  WG  M m  j  {%KGP > ± 0 •  j f d£ -Smfi^J {^KGP-i»±1. ° ) ) Z  nj  KGP-S=\  (3.140) and, K = v  jhNiNjSmfWdQ.  Z ^ • * - ^ ( ^ , ± l ) - M . j ^ , ± l ) . KGP-£=\  jzf -zf dt ni  U=->  •smf(j(4 ^±l,0)  n  KGF  J (3.141)  Note that the parameters {m n ),{m ,n ^ are related to i,j according to the relation in >  n  l  j  J  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) E N D L O O P . 9) END L O O P . 10) END L O O P . 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) L O O P 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 L O O P .  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) R E A D maximum change of degree of cure and maximum time step for a cure sub-step (At' , look at Section 3.6.2) and also the variation in the end temperature p  ( AT \ P  in Equation (3.120), Section 3.7.1) to get the perturbation in the degree  of cure (c \  ~\  F  c  m l+Al  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-ofsub-step temperature by  ^=C,+[v^A/(pC,) ]/(C..C)A/ ' /  e  •  (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 L O O P . 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 ), number of Gaussian IP  integration points in surface direction ( n ), initial degree of cure, cure kinetics model GP  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 o f 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) R E A D the necessary integral from the Integral Matrix and calculate the capacity matrix. (Equation (3.91)). 6) END L O O P . 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) R E A D 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 L O O P . 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) R E A D the necessary integral from the Integral Matrix and calculate H (Equation (3.106)). 6) END L O O P . 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 F 's in the Equation (3.90) for the following cases and saves it an integral t  matrix: 140  CASE#\:f(C)-  CASE#2:  f(C) (3.144)  1 < n, m < n + 2 T  CASE#3:f(C)  CASE #4 : / ( £ )  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 element and calculates the Jacobian ( J ) , its determinant ( | j | ) , and inverse ( J ) at the 1  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. 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) I F 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, , GF 4  iir  GF  I(  (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) E N D L O O P . 12) E N D L O O P . 13) END L O O P . 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  l7 +  ^ l) +  NG = +^(l + «7)(l-#)(»7 + #-l) 2  NG = +^(l + *)0 + £)07 + £ - l ) 3  NG = - I ( l - ) ( l 7  4  +  fl( -^ 7  + l) (3.145)  NG =^0-^)(i+^)0-#) 5  NG = I ( l 6  + ? 7  )(l £)(l-£) +  NG =^0-^)0+^)0+^) 7  NG = I ( l - „ ) ( l £ ) ( l - £ ) . +  S  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 =-(i+<r). 2  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 U P 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) R E A D 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) R E A D the temperature at the start and end of the time step. 5) R E A D the layup information that is passed in from the input file. 6) R E A D 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) S A V E 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\\ from Section 3.8). E  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. A s 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 (^, g , £ ) and through thickness discretization points (£ , Q , £ , CJ 2  3  {  2  3  defined in Section 3.5. This is drawn at integration point KGP with the thickness equal to D . KGP  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 : D C C 3 D 8 ) 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 V I 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 o f thickness 2L -5 (cm) is subjected to air temperature ramp o f r = 3 ( ° C / m i n ) for 30 minutes. The initial temperature is T = 3 0 ( ° C ) . The convective 0  heat transfer coefficient is hj. -\00(w/m K)  for the top surface and h =  2  B  50(w/m K) 2  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: a = k /pC z  z  =0.69/(1580x870) = 0.502 x l 0 ~ (/M /sec) 6  P  2  1 = 0.05/2 = 0.025 (cm) r = 3.0/60 = 0.05 (°C/sec) Bi =h L/k  =(50x0.025)/0.69 = 1.8116  Bi = f\L/k  = (100 x 0.025)/0.69 = 3.6232  B  B  z  T  z  Fo = a t/L  2  z  R = -rL /a 2  = (0.502 x 10" x (30 x 60))/0.025 =1.44 5 8 6  2  = -(0.05 x 0.025 )/0.502x 10" = -62.25 ( ° C ) . 2  z  6  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" . The maximum lag at the same instance is T 5  lag  = -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 Bi = 2.47. t  Based on this Bi and the parameters calculated above, one can read 9 = 0.76 from e  m  Figure 2-4 or its magnified version (Figure 4-1) leading to A T s - 6 2 . 2 5 x 0 . 7 6 = -47.3 l ( ° C ) ,  and the minimum temperature is  7/ = 7/ + r x f + A7/ = 30 + 0.05x(30x60)-47.31 = 72.69(°C) . 0  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  „ 1 2+— +  p  =  _  0  0  9  7  6  1 —  4 = L (l - ^ ) = 0.025 x (1 + 0.0977) = 0.0274 .  Thus resulting in the following parameters: Fo = a tfl  =(0.502x10^ x(30 x 60))/0.0274 =1.2036  R = -rL /a  = - ( 0 . 0 5 x 0.0274 )/o.502x 10"* =-74.7769 (°C)  2  e  z  2  e  2  z  Bi =hLjk e  2  z  = (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 = - 4 7 . 8 6 ( ° C ) T = T + r x t + AT = 30 + 0.05 x (30 x 60) - 47.86 = 72.14(°C) . 0  Numerical solutions Current Shell Element A shell member was modelled with the above thickness and material properties. The inplane 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 o f 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 2L = 1 (cm) 2  sits on an invar tool (material properties in Table 2-1) with 2L -2(cm). X  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 m K)  for the top surface and h - 50(W I m K)  2  2  B  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  14  (the value of R on the next page). e  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 h  ss  =31.37 (W / m K) . So the minimum 2  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 hL ss  the tool would be Bi  n  =  1  K  = 0.0285 . Hence, the location of the adiabatic line (where  minimum temperature occurs) is J C=  1_  \ \ =0.22 . 2+— + — Bij Bi B  B  B  Therefore, the equivalent one-layer-symmetric-boundary parameters are calculated to be 159  Bi = max (Bi , Bi ) x ( l + \g |) = 0.055 e  B  Tl  a  ^ = V ( 1 + |C,|) =-1.84 (°C) 2  Fo =^/(l e  \{ \f=96.7.  +  a  Using the above parameters we can use Figure 2-4 or its magnified version (Figure 4-2) to read Q =18.4 and hence the minimum temperature in the tool is estimated to be m  T = R6 + T =166.05 ( ° C ) . The minimum temperature in the composite part occurs at m  m  x  the interface between the tool and the part and we use the distribution diagram (Figure  2-5) to estimate the value of this temperature. The diagram reads  T -T  — = 0.99 and  hence T =166.42 (°C) . s  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 / m  e  =18.51 .  This value is only 0.5% higher than the more accurate value of 18.4 considered above. Using this value, we find T = 165.94(°C) and T =166.28 ( ° C ) . m  s  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 T = 1 6 6 . 4 ( ° C ) . s  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 T = 166.4 (°C) at the end of the s  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. A s 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. A s 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 h  =\00(w / m K} on the top and h 2  lop  =5®(W lm K} 2  bot  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 o f 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. A s 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/m K) 2  on top and right surface  of the model in Figure 4-20 and h = 50 [w I m K) on the bottom and left surfaces of the 2  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 : D C C 3 D 8 ) 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 ^ e  K  50x0.031 11  A s 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 o f 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 preprocessor 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  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  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)  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)  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)  autoclave  (°C/min) h=20 (W/m K) 2  A=100 (W/m K) 2  h=200 (W/m K) 2  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  h=20 (W/m K)  AT (°C) 21 = 50.0 mm  AT (°C) 21 = 25 mm  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)  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)  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)  (°C/min) 0.5  AT (°C) 21 = 6.25 mm  AT (°C) 21 = 37.5 mm  AT (°C) 2L = 12.5 mm  autoclave  2  A=100 (W/m K) 2  h=200 (W/m K) 2  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  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  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)  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)  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)  autoclave  (°C/min) h=20 (W/m K) 2  h=lOO (W/m K) 2  £=200 (W/m K) 2  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  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  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)  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)  0.5  0.5(0.5)  1(1)  1.9(1.9)  3(3)  4(4)  1  KD  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)  autoclave  (°C/min) H=20 (W/m K) 2  #=100 (W/m K) 2  h=200 (W/m K) 2  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. dc dt  Cure Kinetics Equation  Kc (\-c)" \ { -( O+ CT )} m  C  C  a  a  T  + e  K,=A -^  ,RT  ie  Description  Variable c  Resin degree of cure.  T  Resin Temperature  Units  Values  °K  Fibre volume fraction  0.573  Total resin heat of reaction (c= 0 to 1)  5.40E5  J/kg  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  co  Diffusion constant  -1.684  -  Diffusion constant  5.475E-3  /°K  v =\-v f  A  a  R  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 2 3 4  1 3 5 7  20 10 5 3  453.37 462.12 464.24 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 2 3 4 5  1 3 5 7 20  20 10 5 3 2  372.77 442.58 455.49 460.03 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 2 3 4  1 2 3 4  20 10 5 3  461.08 463.55 464.67 465.16  Table 4-9 - Error study for the minimum temperature in the plate from different methods of solution - CASE II.  Exact Solution Strategy #1 Strategy #2 Superelement Linear Elements Quadratic Elements  Minimum temperature (C)  Error (%)  72.47 72.69 72.14  0.00 0.30 0.46 0.06 0.11 0.06  72.43 72.55 72.43  176  Table 4-10 - Parameters used as input for CASE III. Layer 1  Layer II  0.033  0.033  2.67E-06  5.02E-07  R(C) BI  0.01 -1.24E+00 0.045455  0.005 -1.64E+00 0.72463768  S (W/m sec)  1100  138  r (C/sec) a  (m /sec) 2  z  L(m)  z  Table 4-11 - Comparison of the efficiency of the subcycling method with the classical method of identical time-steping parameters - CASE VII.  Analysis Specifications  Run No  Classical  Subcyclin g  Atmin Atmax ATmax Run Time Error in (sec) final DoC (sec) (sec) (°C)  Run Time Error in (sec) 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 9 - Case II. m  Fo  Figure 4-2 - Magnified section of Figure 2-4 used to read G - Case III. m  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  ?6J  4  8  12  10  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  /  T  tOD  1  Top Surface temperature  Tbot  —200  150 A ra 3  Autoclave air temperature profile  2.100 E  Bottom'surface temperature  0)  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  Series method with  g250  maximum of 11 term^  F E M with 10 elements  £200  3 -*-»  Autoclave air temperature profile\^  SJ150 a §100 50 -\ 0 0  50  100  200  150 time (min)  250  300  Figure 4-7 — Comparison of the convergence of the series method with the FEM solution - CASE V.  12 10  (0 E 8 +•»  _  7 6  I  1 I \  5 c" E a v 3 «  4  No of Terms  o z  2 -i  4 2 0  J 50  100  150 time (min)  ^ 200  250  300  Figure 4-8 - Time step sizes and number of terms used in the series solution - CASE V.  181  300 250 id 200 A  3  o i  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 £1500 E>  I  1000 500 0  0  3000  6000  9000 12000 time (sec)  15000  18000  21000  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 f r o m 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  =100 (W/m K) 2  t  2.5 (cm)  2.5 (cm) h =50 (W/m K) 2  bot  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 0  1 2000  1  4000  1 6000  1  1  1  8000  10000  12000  1 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  subcycling  1600 exact  ,  \  2 1200 ^  classic  i\  % 1400  \  1000 \  g 800  \\  P> 600  \ x  400 200 0 2000  4000  1  1  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. 5.10  LOO (cn)  /-R0.46  (cn)  //-Rl,46  (cn)  (cn)  :  Section B n o  Section A  S  e  c  t  i  o  n  c  3 "O O <+  Tool  -7.20  (cn)-  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 - +1 - +1 - +1 - +1 - +1 - +1 - +1 - +1 - +1 - +1 - +1 - +1  813e+0Z Blle+02 SlOe+OZ 808e+0Z 807e+0Z 80Se+0Z 804etOZ 80Ze+0Z 801e+0Z ?99e+0Z 798e+02 796e+0Z 79Se+02  Figure 4-24 — Temperature distribution in the model at exotherm (t=173 min, temperature in degree C) - Case VIII.  189  3 -*-»  100 m a> a. 80 E a H 60 40 20 0 50  100  150  200  250  300  time (min) 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 5.1  SUMMARY, CONCLUSIONS AND F U T U R E WORK  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 will 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 will 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 o f 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 o f the non-dimensional parameters developed in the chartform 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 o f 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 ThermalProperties 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 ofApplied 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 OneDimensional 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 ed., Elsevier 2000.  [40]  Manko, Z., "Analysis of Thermal Conduction by the Finite Strip Method," Vol. 2, Pineridge,  th  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 ThreeDimensional 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. 241258.  [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. 237260.  [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 ofApplied 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. 5569.  [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., D E , 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 A S M E , 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. 13151325.  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., C R C 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  APPENDIX A . E X A C T THROUGH THICKNESS EIGEN FUNCTIONS FOR NON-DIRICHLET  BOUNDARY  C O N D I T I O N S 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 kind on both surfaces rd  with the formulation T = T(x,y,z,t) T,, =oc T x  xx  +a T y  M  + aT z  (A.l)  +(l/pC )H(x,y  zz  ,z,t)  P  (  A  2 )  dT h(T (x,y,t)-T T  (x, y, +L, t)) = +k —\ z  h(T (x,y,t)-T(x,y,-L,t))  =  B  T(x,y,z,0)  = TXx,z)  Z=+L  -k ^-\ =_ z  z  L  ,  (AA)  where h is the coefficient of convective heat transfer on the surfaces and T and T are T  B  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)  +——\ =T (x,y,t) Bi dz L dT T(x,y,-L,t)-——\ _ =T (x,y,t) . Bi dz z=+L  z=  L  B  208  r  (A.5)  The following change in the unknown function is necessary to get the homogeneous boundary conditions:  T(x, y, z, 0 =  +  fi-^2  . * ^ 7) + u(x,y,z, t) , 2(1 +1 / Bi) L T  (A.6)  which leads to  U  t  =(  a  x  U  xx + y a  — (TT.—GCTT.  \  + «z zz ) + 0 /  —GC T,  TjXx  y  —+  -  2  2(1 + 1/50 L  ~ x^B,xx  (A.7)  1  T,yy )\  2 ~{^B,1  Z  1  1  )  R  T,(  J>» > 0  W  yy  U  2(1 + 1/5/) L  ~ y^B,yy)  a  a  and the following homogeneous boundary conditions: _ , L du 1 «(x,^,+I,0 +—— Bi oz L du 1 u(x,y,-L,t)——\ _ Bi oz  z = + i  Z=  =0 (A.8) =0 ,  L  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 +  u(x,y,+L,t)  .  a  y  U  , y y  +  a  :  L du +—— 5/ oz  U  1 z = + i  x L du, u(x,y,-L,t)-——\ _ =Q. Bi dz r  z=  209  L  , *  =0  (A.10)  Similar to previous problem, we apply the following solution based on separation of variables: u(x,y,z,t) = W(x,y,t)Z{z)  ,  (A.ll)  which leads to W,Z = a ^ Z W,  + aW Z  + WZ"  W  Z"  y  W  W  x  tyy  m  W  ( - ) A  1 2  W  y  which means that  ^- =- A Z a  ,  2  (A.13)  p is some constant. The solution to the above O D E is:  Z = ylsin(/?-) + 5 c o s ( / ? 4 ) • d d  ( - ) A  1 4  The boundary conditions, on the other hand would translate to  B l  V  '  d  (A. 15)  z  Bi 8z " |z  L  Applying the boundary conditions we get  j is odd  cos •0) _  (A. 16) sin  L ,  210  j is even  In order to get { / ^ ] the following equations should be solved in the interval  cotg(/? ) =^ — Bi  j is odd  , Aif/K ~Bi cot g(/? >) = —^-j  .. j is even  w  J  (A. 17)  The general solution to u can then be written as  (A. 18)  u=^W (x,y,t)Z '(z) V)  v  Inputting the above in equation (A.7) we have  =a W% +a0). WV'-a  W}  J)  W  (i)  )  x  y  I  v  J (1 / pC )H(x, h  (x, y, t) = ^  {j)  (A.19)  J  y, z, t) Z(z)dz  P  (j)  L  + h (x,y,t)  T l  ]z \z)Z \z)dz (i  [i  -L +L  K  f  T1,1 - a x TT,xx -a y• TT,yy J) T  T  T  -I  \  1  z  2 + 2(1 + 1/5/) d  A  Z(z)dz (A.20)  \z \z)Z \z)dz [i  [i  1 \ (  T  B , , -  a  T X  B ^ -  a  y  T  B , y y )  2  1  Z  2(1 + 1/5/) d  Z(z)dz  \z \z)Z (z)dz {i  {J)  The above problem is to be solved using the following initial conditions:  211  +L  T  2  2(\ +  T  \  v  -L  W(x,y,0):  T +T  B  T  T  z^  B  Z (z)dz (j)  \IBi)d)  io  (A.21)  +L  \z \z)Z \z)dz U  (i  -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, =a T^ t  x  x  +a T^ y  +a T^  T(x,y,+L,t)  +(1/pC )H  z  P  {x ,y ,z ,t)  (A.22)  +——\ =T (x,y,t) Hi oz z=L  T(x,y,-L,t)  T(x,y,z,0)  (A.23)  T  = T (x,y,t)  (A.24)  B  = T (x,y,z)  (A.25)  i  T is the fluid temperature on the top and T is the essential boundary on the bottom. T  T  The following change in variables could be used to change the boundary conditions to homogeneous boundary conditions:  T(x,y,z,t)  =  T +T {\ + l/Bi) T  B  2 + 1/5/  |  T -T T  B  z  2 + 1/5/ d  and the P D E problem changes to  212  + u(x,y,z,t)  .  (A.26)  ,i  ~ \ x , x x ~^~ y ,yy  U  a  ~( T,, T  U  a  ~ x T,xx a  J r V"  ~^~ z ,zz  U  a  ~ y T,yy  T  a  T  )(  T  -  pC )H(x,y,z,t) P  "  l  '  )+ +  2 + 1/5/  \(  T  "  U  +  l  /  B  1  i  M 2 + 1/5/  (A.27)  1  2 + 1/ Bid z  2 + 1/5/ d  i <. L du | u Ix,y ,+L,t) + 1 _., = 0 Bi dz ~  (A.28)  n  v  ;  u  u(x ,y  i(x,y,z,0)  =  T (x,y,z)-  +l  (A.29)  =0  T +T {\ + l/Bi) T  B  i  |  2 + 1/ Bi  T -T T  B  z  2 + 1/Bid  (A.30) J,=o  Separation of Variables, as in previous sections, will lead to  (A.31)  Z  Z(+Z) +  (A.32)  - ^ U = 0 Bi dz  (A.33)  Z(-I) = 0 The solution to the above O D E is  Z =A sin  f  z^l  f  + 5 cos  z^  (A.34)  Applying the boundary conditions we have  Z ( z ) = c o s ( / ? ) s i n ( / ? - ) + sin(/? )cos(/? —) 0 )  0)  0)  213  (y)  0)  (A.35)  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  V)  =a wV}  W  +a W j -a (  x  W + h (x,y,t)  )  y  y  V  | ( 1 / pC )H h \x,y,t)  =  U  J  L  (x, y, z, t) Z  iJ)  P  (A.37)  v)  2  (z)dz  -L +L  \z (z)Z (z)dz (i)  (i)  +L  1  -+ -  —L\zW( )dz  l  z  2 + 1/5/ 2 +MBid) 1/5/  v  '  (A.38)  +L  \Z \z)Z (z)dz (j  (s)  -L  +i \(T -a T -a T ) B>l  x  Bxx  y  f( I + I/5/3  Byy  -L  \  +L  1  Z (z)dz dj (j)  2 + I/5/3 2 + l/5/  3  jZ°\z)Z (z)dz U)  The above problem is to be solved using the following initial conditions: +L  \T (x,y,z)ini  W \x,y,0) [j  =  -L  T T {\ T+  B  + \IBi)  |  2 + 1/ Bi  T -T T  B  z  2 + \/Bid  Z \z)dz u  +L  \z \z)Z \z)dz (i  (i  (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, = a T^ + a T x  y  T(x,y,0,t) (  z  + (l / pC ) H (x,y, z, t)  zz  P  = T (x,y,t) B  N  T x,y,^2L ' ,t [ l  v  + aT  w  f/l ^  T  J  1=1  T(x,y,z,0)  (A.40)  =T (x,y,t)  = T (x,y,z) i  where the thermal material properties k , H change at the layer interfaces and T and z  B  T are the bottom and top essential boundary values. The temperature profile is T  contentious at the interfaces and its first spatial derivative satisfy the following continuity condition: dT  @ z = 2L \l(d [x  ~dz~  dz  2  + d ...  {x]  (A.41)  [2]  p  Defining a separate z coordinate for each layer will make the mathematical manipulations easier to grasp so we define  (A.42)  z ' =z - 2 l  J  W=i  J  in which case a layerwise temperature distribution function T ^[x ,y ,z ^'\t j would describe the temperature distribution through the thickness of the i layer. The tn  constitutive, boundary and initial value equations ^ ^  = a T® +a Tl x  y  T (x,y,-L ,t) [l]  1  T (x, y, N  T (x,y,z,0) [i]  +  c$T®  would then change to  pC )H^(x,z ,t) [i]  +  P  = T (x,y,t) B  ,t) = T (x, y, t) T  =  Ttt(x,y,z ). [i]  The continuity equations would then read 216  (A.43)  (A.44) V 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,2l^O =  Here /  [ , ] r  (^(^0^+>3 (^)2i) O  +« ' (*,A/) I  1  •  (A.45)  and / j ' are linear functions that satisfy the boundary and interface conditions. 1  These functions are defined in terms of the following parameters:  l>]  S  (A.46)  2 =S-TT ['] [/]  n=l  S  which results in  J B  +  ~  [N]  1  • + 5  [/] [JV]  +1  z  (A.47)  + V"  J  The P D E would then change to:  «!?=(«,«£+«A+«f «S) 0 J  -{T ,,-d T T  x  /^  +  -«/r, )/r ' - ( ^ , [  Txx  ]  w  ) "  [ , ]  (*•  zW  -«?B„)fl  With the following boundary, interface and initial conditions 217  y> '')  n  (A.48) .  u^x,y-L®t)=uW(x,y,+LWt)  =0  u ' (x,+L '\t)=u ' (x, ,-L ' ,t)  i =\,2,...,N-\  l ]  l  l +l]  l +l]  y  (WlA  =(*J' V'$)  (A.49)  i=\,2,..,N-\  +  «I'l(x,^zW,0) = r ( c , z W ) i  J  /  x  ;  - / J ' (* ) * i 1  w  x  (  (*, * o) - f? (*") r ]  r  A  -  5  ° )  (x, y, 0) .  We now try to get the Sturm-Liouville solution for the following homogeneous P D E :  = x a  ,xx  U  +  a  y  U  + «z ,zz  » = 1, 2 , 3 , N  M  ,yy  ,  (A.51)  with interface and boundary conditions as defined in (A.49). A s it is usually the case with similar problems, we write  i/I'l =W (x ,y ,t ) Z ' ( z 1  ]  [ , ]  )  (A.52)  which leads to  Z =A cos(r] ' (z -S f) [i]  li]  [ ]  [i]  [i]  .  (A.53)  N o w input(A.52) and (A.53) into (A.51) to get  W,= W^+aF„-ay  (IT ' ) 1 1  ax  w,  —=a  W  x  w„  + a..  W  y  2  (A.54)  w,  W  which means that  III  III  \2  a  218  =P  2  (A.55)  n =plJa?  .  [i]  (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  +  ^"cos(i7  -  [#1  ^ V s i n ^ '  1  c  o  s  (  ^  '  +  ,  l  -  =  6 i =\,2,..,N  ^ ^ ^ ^  -1  ^ ' ^  i =1,2  -1  or B c o s ( ? 7 J ) - C | s'm{r] d ) = 0.0 i  B  N  1  1  cos(?i d )+C  l  N  N  N  l  sm(n d ) N  N  = 0.0  ^ V ' s i n ^ ' V ' ^ ^ (A.58)  ^  '  -  ^  J  W  ^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 different mode shapes. Here is the problem that should be solved for each j this time:  219  for  W}  =cc Wi?+a W^-(/3 ) W  J)  N  h (x,y,t) (J)  =  +  (A.60)  (j)  y  J  .  + h (x,y,t)  (J) 2  x  4"  \(\lpC )H^.Z ^ {i)  p  '=' -z!'l  '=' -zl'l  (A.61)  '=' ./H  1=1 _jrJ0  The above problem is to be solved using the following initial conditions:  W (x,y,0) {j)  =^  ^  -pj  .  The Finite element formulation would be the same as in previous sections.  220  (A.62)  1 J  ( V  Layer 3 Layer 2  T(z)  Layer 1 AZ  1 1  . X p  Figure A-l — Typical temperature distribution in a three-layer composite.  221  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0063235/manifest

Comment

Related Items