z... .I. .. : L .. fir 25.5. . ' $13 an]? {I 3“,.” arms. .3, 3: £313..) F .4 Q- 2 {Evita a ‘2 3}.“ x3 . .3... A33. 5~u an . t . ‘ . . gflbvxvpnwtm {0.0va .égfifisfia4 a gig“? . ‘ 1 .v: gettin 3% .i: I A ‘ UBRARY 1004 Michigsm State University This is to certify that the dissertation entitled COUPLED ELECTROMAGNETIC THERMAL AND KINETIC MODELING FOR MICROWAVE PROCESSING OF POLYMERS WITH TEMPERATURE- AND CURE- DEPENDENT PERMITI'IVITY USING 3D FEM presented by RENSHENG SUN has been accepted towards fulfillment of the requirements for the Ph.D. degree in Electrical EngineerinL 2,1 /M Major Professor’s Signature 08/22/05 Date MSU is an Affirmative Action/Equal Opportunity Institution PLACE IN REI‘URN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 2/05 c:/CIRC/Dlt00\n.lndd-p.15 COUPLED ELECTROMAGNETIC THERMAL AND KINETIC MODELING FOR MICROWAVE PROCESSING OF POLYMERS WITH TEMPERATURE- AND CURE-DEPENDENT PERMITTIVITY USING 3D FEM By Rensheng Sun A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Electrical and Computer Engineering 2005 ABSTRACT COUPLED ELECTROMAGNETIC THERMAL AND KINETIC MODELING FOR MICROWAVE PROCESSING OF POLYMERS WITH TEMPERATURE- AND CURE-DEPENDENT PERMITTIVITY USING 3D FEM By Rensheng Sun Polymer processing is an important application of microwave heating in industry. Microwave assisted curing of thermoset polymers has the advantage of heating the polymer precursor materials volumetrically and hence can lead to superior cure with efficiency not available with conventional convection heating. However, due to the complex interactions between the electromagnetic fields and the material, achieving the promise of microwave assisted curing is challenging. This is due to the fact that electrical properties (e.g. the complex permittivity) of the material change non-linearly with temperature and composition during the curing process. Hence, the field distribution within the cavity applicator changes as a function of the extent-of-cure and local temperature of the materials being processed. It is vital in modeling the curing process that the microwave power deposition, heat transfer, and polymer curing kinetics be coupled together. In this work, we develop a self-consistent 3D marching-in-time multiphysics model, which includes electromagnetic field distribution, microwave power absorption, heat transfer, and polymer curing kinetics. Temperature- and cure-dependent permittivity and curing kinetics for DGEBA/DDS based on experimental data are explicitly included in the model. An edge-based finite element method (FEM) is implemented as an electromagnetic model to ensure the tangential continuity of electric field and divergence-free condition for source free region, while node-based FEM is used in thermal model to solve for the temperature distribution. The numerical results can be used to determine the time-dependent temperature distribution and curing profile across the polymer sample, as well as the electromagnetic field distribution within the cavity applicator. With the help of this numerical model, robust control strategies can be developed for the polymer curing process. The numerical results are compared with the measured data and a good agreement is achieved. In this research, the electromagnetic modeling of a novel adaptable multi-feed multimode cylindrical cavity applicator is performed, where the spatial distribution of the electric field can be specified a priori to accomplish a desired processing task. The electric field intensity inside the cavity can be tailored by varying the power delivered to each port, and mode-switching can be realized without mechanically adjusting the cavity dimensions. An orthogonal feeding mechanism is developed to reduce the cross coupling between the ports. Numerical simulations are performed for the cavity applicator to verify the theoretical analysis. Copyright by RENSHENG SUN 2005 To my parents. ACKNOWLEDGEMENTS Many peeple have assisted me on this journey towards the completion of the degree. First of all, I am especially grateful to my advisor Dr. Leo Kempel for his constant encouragement, support, and guidance throughout the years. He has been an outstanding mentor and friend. Appreciation is also given to the other committee members, Dr. Martin Hawley, Dr. Edward Rothwell, and Dr. Andre Benard, for their invaluable guidance during the research. Special thanks are given to Dr. Shanker Balasubramaniam for his insightful suggestions and valuable discussions. I wish to thank Dr. Liming Zong, Dr. Shuangjie Zhou, Mr. Pedro Barba, and Dr. Gregory Kobidze for their countless help on research topics. The assistance and friendship of my other fellow graduate students from different departments can not be overstated. I wish to express my sincere thanks to Jun Gao, Chuan Lu, Jun Yuan, Jorge Villa-Giron, Sandra Soto-Caban, Dilip Manda], Nikki Sgn'ccia, and Susan Farhat. I would like to extend my deepest gratitude to my parents, Shouzhi and J ingrong, and my brother, Renhui, for their encouragement and confidence in me during so many years. Finally, but certainly not last, I am especially grateful to my wife Wei for all of her love, support, encouragement, and patience with me throughout my graduate studies. This work was supported by the National Science Foundation under contract number DMI-0200346. vi TABLE OF CONTENTS LIST OF FIGURES .................................................................................. ix LIST OF TABLES ................................................................................... xii CHAPTER 1 INTRODUCTION ...................................................................................... 1 CHAPTER 2 RESEARCH OBJECTIVES ....................................................................... 4 2.1 MICROWAVE HEATING VS. CONVENTIONAL HEATING ...................... 4 2.2 CETK MULTIPHYSICS MODELING .............................................................. 7 2.2.1 Electromagnetic Modeling .......................................................................... 8 2.2.2 Thermal Modeling ...................................................................................... 9 2.2.3 Curing Kinetics Modeling ........................................................................... 9 2.3 ADAPTABLE APPLICATOR DESIGN ......................................................... 10 CHAPTER 3 BACKGROUND ON MICROWAVE HEATING ................................... 12 3.1 ELECTROMAGNETIC FIELDS IN A MICROWAVE ENCLOSURE ......... 12 3.1.1 Maxwell’s Equations ................................................................................ 12 3.1.2 Resonant Modes in a Cylindrical Single Mode Cavity ............................. 15 3.1.3 Electromagnetic Fields in a Cylindrical Single Mode Cavity .................. 18 3.2 MICROWAVE/MATERIALS INTERACTIONS ........................................... 19 3.2.1 Mechanisms of Microwave/Materials Interactions ................................... 19 3.2.2 Dielectric Properties .................................................................................. 22 CHAPTER 4 FINITE ELEMENT METHOD ................................................................ 26 4.1 BASIS FUNCTIONS F OR NODE-BASED FEM ........................................... 26 4.1.1 Triangular Elements .................................................................................. 27 4. l .2 Tetrahedral Elements ................................................................................ 29 4.2 BASIS FUNCTIONS FOR EDGE-BASED FEM ............................................ 32 vii 4.2.1 Triangular Elements .................................................................................. 33 4.2.2 Tetrahedral Elements ................................................................................ 34 4.3 CODE VALIDATION WITH WAVEGUIDE APPLICATIONS ................... 37 4.3.1 Loaded Rectangular Waveguide ............................................................... 43 4.3.2 Waveguide Bandpass Filters ..................................................................... 46 CHAPTER 5 CETK MULTIPHYSICS MODELING .................................................... 51 5.1 OVERVIEW ..................................................................................................... 51 5.2 ENERGY BALANCE EQUATION ................................................................. 52 5.3 CURING KINETICS ........................................................................................ 54 5.4 DIELECTRIC PROPERTIES ........................................................................... 57 5.5 ELECTROMAGNETIC MODEL .................................................................... 66 5.6 HEAT TRANSFER MODEL ........................................................................... 70 5.7 EXPERIMENTAL SETUP ............................................................................... 74 5.8 NUMERICAL RESULTS VS. MEASURED DATA ...................................... 77 CHAPTER 6 ADAPTABLE MULTIMODE APPLICATOR MODELING ................. 90 6.1 INTRODUCTION ............................................................................................ 90 6.2 MODE-SWITCHIN G ....................................................................................... 93 6.3 PORT-TO-PORT COUPLING ANALYSIS .................................................... 99 CHAPTER 7 CONCLUSIONS ..................................................................................... 103 CHAPTER 8 RECOMMENDATIONS FOR FUTURE WORK ................................. 104 APPENDIX EXPRESSIONS FOR INTEGRAL TERM WITH DOT PRODUCT ....... 106 REFERENCES ............................................................................................................... 108 viii LIST OF FIGURES Figure 3.1 TM mode diagram for an empty cavity with a radius of 10.75 cm. ................ 16 Figure 3.2 TB mode diagram for an empty cavity with a radius of 10.75 cm. ................. 17 Figure 4.1 Linear tetrahedral element. .............................................................................. 30 Figure 4.2 Rectangular waveguide loaded with an obstacle [22]. .................................... 38 Figure 4.3 Loaded rectangular waveguide. ....................................................................... 44 Figure 4.4 S-parameters for the completely filled waveguide. ......................................... 45 Figure 4.5 Schematic view of the TEm l-mode bandpass filter [32]. ................................ 47 Figure 4.6 Meshed geometry of the bandpass filter .......................................................... 48 Figure 4.7 Calculated magnitude of 821 for the bandpass filter with FEM ....................... 49 Figure 4.8 Magnitude of S21 for the bandpass filter from [32]. ........................................ 50 Figure 5.1 Extent-of-cure vs. time and temperature for DGEBA/DDS (markers represent experimental data and lines represent calculated data from the model) ........ 56 Figure 5.2 Time derivative of extent-of-cure vs. extent-of-cure. ..................................... 56 Figure 5.3 Measured dielectric properties for the DGEBA/DDS. .................................... 59 Figure 5.4 DGEBA resin curing mechanism. ................................................................... 60 Figure 5.5 Comparison between the experimental and calculated permittivity of DGEBA/DDS. ............................................................................................... 62 Figure 5.6 Calculated dielectric constant for DGEBA/DDS based on Davidson-Cole model for a wider temperature range ............................................................. 64 Figure 5.7 Calculated dielectric loss factor for DGEBA/DDS based on Davidson-Cole model for a wider temperature range ............................................................. 65 Figure 5.8 Calculated dielectric loss factor for DGEBA/DDS based on modified Davidson-Cole model for a wider temperature range. .................................. 66 Figure 5.9 Coax-fed single-mode cavity with the epoxy sample loaded. ......................... 67 Figure 5.10 Pattern of the total electric field distribution for TMmz mode inside the cavity. ............................................................................................................ 68 Figure 5.11 Flow chart for the coupled electro-, thermo- and kinetic model. .................. 73 Figure 5.12 Schematic illustration of the microwave processing and diagnostic system. 75 Figure 5.13 Microwave processing and diagnostic system. ............................................. 75 Figure 5.14 Calculated and measured temperature vs. time for the epoxy sample. ......... 80 Figure 5.15 Calculated extent-of-cure vs. time for the epoxy sample. ............................. 81 Figure 5.16 Calculated dielectric constant vs. time for the epoxy sample. ...................... 82 Figure 5.17 Calculated dielectric loss factor vs. time for the epoxy sample. ................... 83 Figure 5.18 Parametric study on the heat transfer coefficient h for the epoxy sample ..... 84 Figure 5.19 Simulation with and without reaction heat for the epoxy sample. ................ 85 Figure 5.20 Simulation on the variation effects of the complex permittivity for the epoxy sample. ........................................................................................................... 86 Figure 5.21 Temperature distribution at cross section of point B at t = 100 seconds ....... 87 Figure 5.22 Temperature distribution at cross section of point B at t = 300 seconds ....... 87 Figure 5.23 Temperature distribution at cross section of point B at t = 400 seconds ....... 88 Figure 5.24 Temperature distribution at cross section of point B at t = 500 seconds ....... 88 Figure 5.25 Temperature distribution at cross section of point B at t = 700 seconds ....... 89 Figure 5.26 Temperature distribution at cross section of point B at t = 1000 seconds ..... 89 Figure 6.1 Electric field distribution for different combination of TM and TE modes. 96 Figure 6.2 Parametric study on the weighting factor or for the combination of the two modes ............................................................................................................. 98 Figure 6.3 The two-port cavity fed with waveguides via irises. ..................................... 101 Figure 6.4 TB mode (Sn and S2,) for the two-port cavity. ............................................. 102 Figure 6.5 TM mode (S22 and 812) for the two-port cavity. ............................................ 102 xi LIST OF TABLES Table 2.1 Comparison between microwave and thermal heating methods ........................ 5 Table 4.1 Edge numbering for a triangular element ......................................................... 34 Table 4.2 Edge definition for a tetrahedral element ......................................................... 35 Table 5.1 Epoxy curing reaction constants and parameters .............................................. 55 Table 5.2 Properties of the epoxy resin and the curing agent ........................................... 5 8 Table 5 .3 DSC results of the curing DGEBA/DDS system .............................................. 59 Table 5.4 Values of the parameters for the curing DGEBA/DDS system ........................ 63 Table 5.5 Parameters of the material properties in the heat transfer model ..................... 70 xii CHAPTER 1 INTRODUCTION Polymer processing is an important application of microwave heating in industry. Microwave assisted curing of thermoset has the advantage of heating the polymer precursor materials volumetrically and hence can lead to superior cure with efficiency not available with conventional convection heating. However, due to the complex interactions between the electromagnetic fields and the material, achieving the promise of microwave assisted curing is challenging. This is due to the fact that electrical prOperties (e.g. the complex permittivity) of the material change non-linearly with temperature and composition during the curing process. Hence, the field distribution within the cavity applicator changes as a function of the extent-of-cure and local temperature of the materials being processed. It is vital that in modeling the curing process microwave power deposition, heat transfer, and polymer curing kinetics be coupled together. In this work, first in Chapter 2, the research objectives of this interdisciplinary project are discussed in a little more detail, with the emphasis on the multiphysics modeling and the adaptable multimode applicator design. In Chapter 3, a review of the background on microwave processing of materials is performed, including the electromagnetic fields in a cylindrical single mode cavity and microwave/materials interactions during microwave heating. In Chapter 4, node-based finite element method (FEM) and edge-based FEM are presented since they will be used for heat transfer model and electromagnetic model separately in the multiphysics model in Chapter 5. Some numerical results generated from EM model are compared with existing data in literature for the purpose of validation of the code. In Chapter 5, which is the main body of the dissertation, we develop a self- consistent 3D marching-in-time multiphysics model, which includes electromagnetic field distribution, microwave power absorption, heat transfer, and polymer curing kinetics. Temperature- and cure-dependent permittivity and curing kinetics for diglycidyl ether of bisphenol A (DGEBA) / diaminodiphenyl sulfone (DDS) based on experimental data are explicitly included in the model. Edge-based FEM is implemented for electromagnetic model to ensure the tangential continuity of electric field and divergence-free condition for source free region, while node-based FEM is used in thermal model to solve for the temperature distribution. The numerical results can be used to determine the time-dependent temperature distribution and curing profile across the polymer sample, as well as the electromagnetic field distribution within the cavity applicator. With the help of this numerical model, robust control strategies will be developed for the polymer curing process. The numerical results are compared with the measured data and a good agreement is achieved. In Chapter 6, we present the electromagnetic modeling of a novel adaptable multi- feed multimode cylindrical cavity applicator, where the spatial distribution of the electric field can be specified a priori to accomplish a desired processing task. The electric field intensity inside the cavity can be tailored by varying the power delivered to each port, and mode-switching can be realized without mechanically adjusting the cavity dimensions. An orthogonal feeding mechanism is developed to reduce the cross coupling between the ports. Numerical simulations are performed for the cavity applicator to verify the theoretical analysis. Finally, in Chapter 7 and Chapter 8, we summarize the conclusions and outline the future work, respectively. CHAPTER 2 RESEARCH OBJECTIVES 2.1 Microwave Heating vs. Conventional Heating Microwaves are one form of electromagnetic radiation. It is a wave motion associated with electric and magnetic forces. Microwave refers to electromagnetic waves in a frequency range from 300MHz to 300GHz or a free-space wavelength range from 1m to 1mm. Heating is one of the major non-communication applications of microwaves. The fundamental electromagnetic constitutive property of nonmagnetic materials for microwave heating and diagnosis is complex dielectric constant (8* = 8' - js"). The real part of the complex dielectric constant is the dielectric constant, which is related to the microwave energy stored in the materials. The imaginary part is dielectric loss factor, which is related to microwave energy dissipated as heat in materials. The dielectric loss factor of materials is generally due to contributions from the motion of dipoles and charges, conductivity, etc. Polymers have polar groups to interact with electromagnetic fields and exhibit dielectric relaxation. These polar groups can absorb microwave energy directly and the localized heating on the reactive polar sites can initiate or promote polymerizations that require heat. Microwave processing of materials has been studied as an attractive alternative to conventional thermal processing. Thermal heating is a surface-driven, non-selective process. The heating efficiency is controlled by the heat transfer coefficient at the material surface and the thermal conductivity of the material. During thermal heating, heat flows from the surface to the interior of the material. This tends to cause remarkable temperature gradients in thick materials. Residual thermal stresses resulting from large temperature gradients can impair the physical and mechanical properties of the materials. In addition, the production cycle is long because of the difficulty in heating poor thermal conductors such as polymers. Microwave heating offers a number of advantages over thermal heating in a wide range of heating applications. A comparison between microwave and thermal heating is summarized in Table 2.1. Microwave heating is selective, instantaneous, and volumetric with heat loss at the boundaries while thermal heating is nonselective and depends on temperature gradient. Microwave heating can be easily controlled by rapid changes in the applied electric field whereas thermal heating is characterized with long lag times and difficult composite cure control. The heat source of microwave heating can be readily removed to prevent thermal excursion. Microwave processing has potential for rapid processing of thick-section and complex-shaped composites. Table 2.1 Comparison between microwave and thermal heating methods Thermal heating Microwave Heating Heat conduction/convection Energy coupling/transport Surface heating Molecular level coupling Slow Fast Surface Volumetric Non-selective Selective Less property dependent Material property dependent Surface temperature control Intelligent control Established technology Emerging technology In recent years, microwave has been investigated as an attractive alternative (and efficient) energy source for material processing, including polymers and composites. Results have been reported on enhanced polymerization rates [1-2], increased glass transition temperatures of cured epoxy [1], improved interfacial bonding between graphite fiber and polymer matrix [3], and increased mechanical properties of the composites [4]. Practical application of microwave heating in industry is however still limited due to difficulties in achieving optimal temperature control for microwave heating. Although microwaves can penetrate materials, thus resulting in volumetric heating, it is not easy to achieve desired field distribution within the material. This is especially true for thermoset polymers. The processing of these polymers involves an autocatalytic reaction (cure) with a large heat of reaction. Curing is difficult to control and uneven cure has detrimental effect on material properties such as mechanical strength. A uniform or controlled temperature history is also desired for reducing residual stress in polymeric composites. This is difficult to achieve without active cooling. It is however possible to initiate cure at precise locations or almost uniformly. In addition to temperature uniformity (or a temperature distribution that minimizes residual stress buildup), curing temperature level is desired to be precisely determined in order to obtain optimal product qualities. For example, a precise temperature cycle should be applied to carbon/epoxy prepregs in order to obtain satisfactory manufactured components in which the matrix crosslinking is almost complete [5]. Hence, the industrial use of microwave as an alternative (energy efficient) approach to inefficient pressurized ovens has been impeded by the lack of proper applicator design, process modeling, and control/monitoring methods. The objective of this interdisciplinary research therefore includes three firndamental elements: 1. Development of a coupled e1ectromagnetic-thermal-kinetic (CETK) multiphysics model that will be used to understand the curing process and to develop robust control strategies. This combined effects model is unique in that the permittivity (and hence electromagnetic field) variations attributed to temperature and composition changes is explicitly included. 2. Design of an adaptable multi-feed multimode microwave applicator, where the spatial distribution of the electric field can be specified a priori to accomplish a desired task. 3. Development of a process control methodology that utilizes non-invasive process monitoring. In this dissertation, the multiphysics process modeling and the adaptable multimode applicator design are emphasized. The process control strategy will be developed later based on the numerical results from the multiphysics model and the new applicator system. Next, the numerical modeling and the applicator design will be discussed in more detail. 2.2 CETK multiphysics modeling To gain a fundamental understanding of polymer/microwave interactions, to reduce risk, and to help develop control/diagnosis scheme, a multiphysics process simulation package is required. This simulation will include all the relevant physical and chemical properties required to assess the field distribution within the applicator during processing and the resulting temperature profile of the polymer sample. The simulation is supported by three parts: electromagnetic, thermal and curing kinetics modeling. The critical aspect of our modeling approach relies on the fact that electromagnetic interactions occur much faster, 0(ns), than thermal effects, 0(ns), or the kinetic/composition effects, 0(5). Hence, the electromagnetic analysis can be performed with an efficient steady-state (e.g. time-harmonic) analysis method whose results are used in solving the energy balance equation for heat transfer problem. The numerical results from this process model can be used to determine the time-dependent temperature distribution and curing profile across the polymer sample, as well as the electromagnetic field distribution within the cavity applicator. With the help of this numerical model, robust control strategies can be developed for the polymer curing process. 2.2.1 Electromagnetic Modeling We are using the edge-based finite element method (FEM) with tetrahedral elements for the applicator and polymer sample since this method offers a high degree of geometrical and electrical fidelity. In particular, the spatial variation in permittivity can readily be accommodated and by using “snapshots-in-time”, time-dependant thermal and composition influences on the local permittivity can be included in the model. Recall that such quasi-stationary analysis approaches are valid since the evolution of thermal and composition effects are relatively slow compared to the rate of electromagnetic field oscillation within the applicator. The finite element method is well-suited for resonant cavity simulation since it can incorporate the presence of a sample as well as determine the response of the system to single frequency excitation in an efficient manner. Once the fields are determined throughout the cavity for a given excitation, local power absorption rate can be determined for the heat transfer model. 2.2.2 Thermal Modeling Node-based FEM with tetrahedral elements is used for heat transfer model. The temperature distribution will be obtained by solving the energy balance equation, which includes microwave power absorption, heat transfer and heat generation from chemical reactions. After the node-based FEM representing spatial variations is performed, a first- order ordinary differential equation (ODE) system is obtained and the finite difference scheme for temporal variations will be used for time-marching process. The primary heat transfer mechanism within the sample is conduction, while the free convection boundary condition will be applied between the Teflon/air or epoxy/air interfaces. 2.2.3 Curing Kinetics Modeling As the polymer heats, chemical reactions occur that result in composition changes. Such composition changes lead to different dielectric properties with consequential impact on the electromagnetic fields and power absorption (input power) within the applicator and material. In addition, reaction heat needs to be included in the energy balance equation. For polymer curing, existing reaction kinetics models include an nth-order reaction model [6-7] and an autocatalyzed reaction model [8-10]. The curing reactions for an epoxy are characteristic of autocatalyzed reactions. The reactions between amines and epoxide are autocatalyzed by the hydroxide group formed in the reactions. The generalized form of the autocatalyzed reaction kinetics for epoxy curing was derived from the reaction kinetic mechanism [11]. The extent of reaction can be measured experimentally as a function of time and temperature. The reaction constants and parameters can be determined with regression of experimental data. 2.3 Adaptable Applicator Design Microwave heating is nearly instantaneous, volumetric, and uneven in nature. Heating is most effective when resonant modes are excited inside the microwave applicator, characterized by low power reflectance. Different resonant modes have different electric field distributions, which are typically not uniform. In a multi-mode oven, the applicator is over-moded, time averaged, and heating uniformity can be obtained by frequency sweeping [12]. Single-mode cavities can provide uniform heating using the mode-switching method, in which several modes with complementary heating patterns are alternatively excited [13-14]. With a fixed frequency microwave power source, mode switching can be achieved by mechanically adjusting the volume of the cavity. This mechanical process slows the response of the system to temperature changes. With a variable frequency power source, modes can be changed by changing frequency. As a result of the instantaneous variable frequency mode switching, not only the speed of the process but also the controllability of the process is much improved. Variable frequency mode switched heating of flat graphite/epoxy panels has been presented in a recent paper [15] and dissertation [16]. Results showed the success of controlled uniform heating. However, variable frequency sources are inherently less efficient than the single frequency versions; the equipment is also very expensive. In this research effort, with our novel adaptable multi-feed cavity applicator, mode switching can be obtained by varying the power delivered to multiple ports, hence eliminating the need for mechanical applicator adjustments. Specifically, this microwave applicator is designed with the goal of: 10 1. Reducing the processing cycle time. 2. Maximizing the processing efficiency (e.g. maximize power delivered to the material at desired locations and at the desired moment in time). 3. Improving processed material uniformity. The applicator developed in this effort will be unique since: 1. It is a multi-port, multi-mode applicator that operates at a single frequency using a coherent high-power source. 2. Since the input power is divided between several ports, the applicator can be constructed using components rated for low power application. The use of such components will reduce applicator construction cost. In this adaptable multimode cavity applicator, we can tailor the field intensity to achieve high field strength in the regions of the applicator requiring high fields (and hence regions where heating is desirable) and low field intensities in regions where heating is not desired. And as we mentioned earlier, this can be done by just varying the power delivered to each port, and the mode-switching can be realized without mechanically adjusting the cavity dimensions. An orthogonal feeding mechanism is developed to reduce the cross coupling between the ports. Numerical simulations are performed for the cavity applicator to verify the theoretical analysis. 11 CHAPTER 3 BACKGROUND ON MICROWAVE HEATING 3.1 Electromagnetic Fields in a Microwave Enclosure Electromagnetic field strength and distribution patterns are essential factors that influence microwave heating efficiency and uniformity. They are determined by microwave operating conditions, applicator dimensions, and material properties. To understand microwave heating characteristics, fundamentals of microwave processing are reviewed. 3.1.1 Maxwell’s Equations The basic laws governing electromagnetic wave phenomena are Maxwell's equations, which describe the relations and variations of the electric and magnetic fields, charges, and currents associated with electromagnetic waves. Maxwell's equations can be written in either differential or integral form. The differential form is widely used to solve electromagnetic boundary-value problems and is usefirl for FEM formulations. These coupled partial differential equations (PDE) are given by V x E = 3%]; (Faraday's law) (3.1) V x H = J + gt)- (Ampere's law) (3.2) V - D = p (Gauss law) (3.3) V -B = 0 (Gauss law — magnetic) (3.4) 12 where E is the electric field intensity, H is the magnetic field intensity, D is the electric displacement density or electric flux density, B is the magnetic flux density, J is the electric current density, and p is the charge density. D is defined as: D = 50E + P (3.5) where so is the dielectric constant of free space, P is the volume polarization, the measure of the density of electric dipoles. B can be expressed as: B = p0(H + M) (3.6) where in) is the permeability of free space, H is the magnetic field intensity, and M is the volume magnetization, the measure of the density of magnetic dipoles in the material. In a simple isotropic medium, these field quantities are related as follows: D = 5E (3.7) B = pH (3.8) where a is the dielectric constant, and p. is the permeability. In addition to the Maxwell's Equations, the Equation of Continuity holds due to the conservation of electric charge: a V-J+— =0 3.9 atp ( ) In Maxwell's equations, only two of the PDEs are independent. Usually Equations 3.1 and 3.2 are used with Equation 3.9 to solve for electromagnetic fields. Maxwell’s equations are first-order differential equations with E and B coupled. They can be converted into decoupled second-order wave equations through mathematical manipulations: S 86 {B}: (a) ,u at (3-10) — ,uV x J s where o is the conductivity, and J’ is the source current term. In a source free region, Equations 3.10 becomes: 2 2 (3 a E V — a—— 5— = 0 3.11 l .U 6t 1“ 6:2]IB} ( ) Equations 3.1-3.4 are the time-domain representation of Maxwell's equations. If the source functions, J(r, t) and p(r, t), oscillate with a constant angular frequency a) in the system, all the fields will oscillate at the same frequency. The Maxwell’s equations can be written in time-harmonic form: rV x E(r) = —ja)B(r) < V x H(r) = J(r) + ij(r) V - 1)(r) = p(r) [V - B(r) = 0 (3.12) For the time-harmonic case, Equation 3.10 becomes the vector Helmholtz equation and Equation 3.11 becomes the vector Helmholtz equation in source-free region: [v2 + wzyg‘]{§} = 0 (3.13) 5‘ = 50—15;) (3.14) where 8* is the complex dielectric constant. Note that isotropic, non-magnetic materials are assumed throughout this thesis. l4 3.1.2 Resonant Modes in a Cylindrical Single Mode Cavity In a cylindrical single mode cavity, there are two types of resonant modes, transverse electric (TE) and transverse magnetic (TM). In TE modes, the electric field components are transverse, and the magnetic field components are parallel to the direction of wave propagation, which is the axial direction. In TM modes, the electric field components are parallel, and the magnetic field components are transverse to the direction of wave propagation. Three subscripts, n, p, and q, are used to describe the corresponding mode in an empty cavity, i.e. TEnpq and TMan. The subscript n denotes the number of the periodicity in the circumferential direction, n=0,1,2...; p denotes the number of zeroes of the radial wave function, p=1,2,3...; q denotes the number of half wavelengths along the length of the equivalent circular waveguide, q=0,1,2... for TM modes and q=1,2,3. .. for TE modes. Theoretically the relationship between the frequency and cavity diameter and length for a given resonant mode can be calculated in closed form. The equations for TE and TM modal resonant frequency for an empty cylindrical single mode cavity [17] are: 2 (f)npq= 2m :/—‘/x'np HUM) (3.15) 2 (f)npq= 2m ——:/=\/x,,p 2 +[q—h—m) (3.16) where f is frequency, a is cavity radius, 11 is cavity height, x and x' are the tabulated "P "P zeroes of the Bessel’s function and the derivative of the Bessel’s function, respectively. Fig. 3.1 and Fig. 3.2 show the calculated resonant frequency as a function of cavity length in an empty cavity with a radius of 10.75 cm for TM and TE modes, respectively. 15 TM modes In the a = 10.75 cm empty cavity 4.5 3.5 Frequency (GHz) I" u- 4 6 8 10 12 14 16 18 Longthofcavltflcm) [+TM010 +TM011 ~e—TM012 ——TM013 -a—TM014 +TM020 l—t—TM021—o—TM022 ——TM023 +TM110 TM111—a—TM112 >+TM113 +TM114 ., TM120 +TM121 —a—TM210 —--TM211 --o—TM212 I >TM213 +TM310 —9—TM311 +TM312 Figure 3.1 TM mode diagram for an empty cavity with a radius of 10.75 cm. I 1 TE modes inthe a= 10.75 cmenptycavity Frequency (GI-1:) 4 6 8 10 12 14 16 18 Cavity length (cm) :9—TE011 +TE01’2' TEO13 T715014 +1E021 -e—TE111 +1511} +15113 —-—TE114 +TE121 —-+—TE122 —x—TE123 ,, - TE211 7+.TE212 +1153“ +TE411 +TE412 —o—TE413 —‘-—TE511 -—TE512 Figure 3.2 TB mode diagram for an empty cavity with a radius of 10.75 cm. l I e» TE213 TE214 +TE221 —-—TE222 --e—TE311 I TE312 +TE31:I 3.1.3 Electromagnetic Fields in a Cylindrical Single Mode Cavity In a homogeneous, source-free cylindrical single mode cavity with perfectly conducting walls, the electromagnetic fields inside the cavity can be derived from Maxwell's equations and boundary conditions. For TE modes, the electromagnetic field components inside an empty cavity are [17]: _ 1 ‘32qu E ___1_3"1fl ” jaw apaz " p as 2 11¢: 1 la V’npq E¢=aw_”-p?_ (3.17) jaw/9 (M62 6p 2 l a V’npq 2 H = +k E =0 2 jaw 622 Wnpq) z where (p, (t), z) are the cylindrical coordinates, k2 = (02,11 6 , and 1;!an is the wave potential for TEnpq modes: J (x), Wm" '6’” ) (318) = — s1n —z . Wnpq n a '0 cosn¢ h where a is the radius, and h is the height of the cavity. In TM modes, the field components inside an empty cavity are [17]: 2 = 1 a V’npq H :lail’npq p jars apaz p p 5;» I 1621/an aWnpq E¢= . —- H¢=— (3.19) ngp W52 5,0 18 2 E = l (6 V’npq “ jw‘9 622 + kzu/npq) Hz = o where l/lnpq is the wave potential for Tanq modes: x . mg =1” (—"’—’-p>{sm"¢}cos<flz) (3.20) a cos n h When the cavity is loaded with materials, equations (3.17)-(3.20) are no longer applicable and numerical techniques can be used to solve for the electromagnetic field inside the cavity. The most widely used numerical techniques include the method of moments (MOM), the finite element method (FEM), and the finite difference time domain (FDTD). 3.2 Microwave/Materials Interactions 3.2.1 Mechanisms of Microwave/Materials Interactions Non-magnetic materials are classified into conductors, semiconductors and dielectrics according to their electric conductivity. Conductors contain free charges, which are conducted inside the material with alternating electric fields so that a conductive current is induced. Electromagnetic energy is dissipated into the materials while the conduction current is in phase with the electric field inside the materials. Dissipated energy is proportional to conductivity and the square of the electric field strength. Conduction requires long-range transport of charges. In dielectric materials, electric dipoles are created when an external electric field is applied and they will rotate until they are aligned in the direction of the field. 19 Therefore, the normal random orientation of the dipoles becomes ordered. These ordered polar segments tend to relax and oscillate with the field. The energy used to hold the dipoles in place is dissipated as heat into the material while the relaxation motion of dipoles is out of phase with the oscillation of the electric field. Both the conduction and the electric dipole movement cause losses and are responsible for heat generation during microwave processing. The contribution of each loss mechanism largely depends on the types of materials and operating frequencies. Generally, conduction loss is dominant at low frequencies while polarization loss is important at high frequencies. Most dielectric materials can generate heat via both loss mechanisms. There are mainly four different kinds of dielectric polarization: 1. Electron or optical polarization occurs at extremely high frequencies, close to ultraviolet range of electromagnetic spectrum [18, 19]. It refers to the displacement of the electron cloud center of an atom relative to the center of the nucleus, caused by an external electric field. When no electric field is applied, the center of positive charges (nucleus) coincides with the center of negative charges (electron cloud). When an external electric field is present, the electrons are pushed away from their original orbits and electric dipoles are created. Atomic polarization is also referred to as ionic polarization. It occurs in the infrared region of the electromagnetic spectrum. This type of polarization is usually observed in molecules consisting of two different kinds of atoms. When an external electric field is applied, the positive charges move in the direction of the field and the negative ones move in the opposite direction. This 20 mainly causes the bending and twisting motion of molecules. Atomic polarization can occur in both non-ionic and ionic materials. The magnitudes of atomic polarization in non-ionic materials are much less than that in ionic or partially ionic materials. Orientation or dipole alignment polarization occurs in the microwave range of the electromagnetic spectrum. It is the dominant polarization mechanism in microwave processing of dielectrics. Orientation polarization is usually observed when dipolar or polar molecules are placed in an electric field. At the presence of external electric field, the dipoles will rotate until they are aligned in the direction of the field. The dipolar rotation of molecules is accompanied by intermolecular friction, which is responsible for heat generation. Orientation polarization is firndamentally different from electronic and atomic polarization. The latter is due to the fact that the external field induces dipole moments and exerts displacing force on the electrons and atoms, while the orientation polarization is because of the torque action of the field on the pre-existing permanent dipole moments of the molecules. Interfacial or space charge polarization occurs at low frequencies, e.g. radio frequency (RF). It is a fundamental polarization mode in semiconductors. This type of polarization is caused by the migration of charges inside and at the interface of dielectrics under a large scale field. 21 3.2.2 Dielectric Properties Most polymers and composites are non-magnetic materials. The electromagnetic energy loss is only dependent on the electric field. Incident electromagnetic fields can interact with conductive and nonconductive materials. The fundamental electromagnetic property of non-magnetic materials for microwave heating and diagnosis is the complex dielectric constant: £=E'—j£ (3.17) The real part of the complex dielectric constant is dielectric constant. The larger the polarizability of a molecule is, the larger its dielectric constant. The imaginary part is dielectric loss factor, which is related to energy dissipated as heat in the materials. Usually, the relative values with respect to the dielectric constant of free space are used: 5:80(5r_jgeff) (3-18) where so is the dielectric constant of free space, a", is the relative dielectric constant, and 82/]: is the effective relative dielectric loss factor. Typically, the loss factor of materials consists of both polarization and conduction loss. The polarization loss is further contributed by all four polarization mechanisms mentioned earlier. The effective relative loss factor is expressed as [19]: egg...) = g;(w)+g;(w)+g; +5; +—"—- (3.19) 800) where the subscripts d, e, a, and 5 refer to dipolar, electronic, atomic and space charge polarization, respectively. The loss factor is a function of material structures, compositions, angular frequency, temperature, pressure, etc. 22 The ratio of the effective loss factor to the dielectric constant is defined as the loss tangent, which is also commonly used to describe dielectric losses: 6" tan 5,17 = if- (3.20) gr When introduced into a microwave field, materials will interact with the oscillating electromagnetic field at the molecular level. Different materials will have different responses to the microwave irradiation. Microwave heating of conductive materials, such as carbon fibers and acid solutions, is mainly due to the interaction of the motion of ions or electrons with the electric field. However, conductors with high conductivity will reflect the incident microwaves at their surface and accordingly can not be effectively heated. The fields attenuate towards the interior of the material due to skin effect, also. For conducting materials, the conducting electrons are limited in the skin area to some extent, which is called the skin depth, d3. Defined as the distance into the sample, at which the electric field strength is reduced to He, the skin depth is given by [17]: ds = 1 (3.21) 1 1/2 [Ewflo/lraj where a) is the angular frequency of the electromagnetic waves in rad/sec, 110 (= 41t><10'7 H/m) is the permeability of the free space, ,u, is the relative permeability, and 0" is the conductivity of the conductor in mhos/m. For example, graphite has 0' = 7x104 mhos/m and d5. = 38.4 ,u m at 2.45 GHz. The skin depth decreases as frequency increases. For a perfect conductor, the electric field is reflected and no electric field is induced inside a 23 perfect conductor at any frequency. Therefore, no electromagnetic energy will be dissipated even though the conductivity of the perfect conductor is infinite. Since the vast majority of the microwave power delivered to the cavity in this work is dissipated in the polymer part rather than the cavity brass, perfect electric conductor will be assumed. Microwave heating of nonconductive materials, such as polymers, glass fibers, and Kevlar fibers, is mainly due to the interaction of the motion of dipoles with the alternating electric field. Microwave processing of thermosets is self-adjusting. As the crosslinking occurs, the mobility of dipoles decreases because of the “trapping” or reaction and the dielectric loss factor decreases. Energy absorbed by crosslinking molecules decreases and microwaves are concentrated in unreacted molecules. During microwave processing of thermoplastics, the dielectric loss factor usually increases with temperature and thermal runaway is likely to occur. Thermal runaway can be prevented by decreasing or even turning off power at a temperature close to thermal excursion. Microwave heating selectivity of polymer composites depends on the magnitude of dielectric loss factor of polymers and fibers. When non-conducting fibers, such as glass, are used, microwaves will selectively heat the polymer matrix. When conducting fibers like graphite are used, energy is preferably absorbed by the conductive fibers and heat is locally conducted from the fiber to the matrix. In this case, loss factor is mainly due to the fiber conductivity and can not be used to diagnose the curing process of the low loss matrix materials. Dielectric measurement of epoxy curing systems has shown that generally both the dielectric constant and dielectric loss factor increase with temperature and decrease with extent of reaction [20]. This dependence on temperature and extent of reaction is 24 non-linear. During microwave processing, the dielectric properties of materials change as a result of heating and reaction. This affects the electric field strength and power absorption in the materials. The change in electric field and power absorption in turn affects the temperature and extent of reaction inside the materials. Thus, the modeling of microwave heating is a coupled non-linear problem, which involves Maxwell’s equations for solving the electric field strength, a heat transfer equation for obtaining the temperature distribution inside the material, and a reaction kinetic equation for calculating extent of reaction. 25 CHAPTER 4 FINITE ELEMENT METHOD The finite element method (FEM) is widely used in computational electromagnetics since it is very efficient on modeling full three-dimensional volumes that have complex geometrical details and material inhomogeneities. Since the edge- based FEM for electromagnetics and node-based FEM for heat transfer, respectively, are used herein in Chapter 5, it is necessary to briefly review the basis functions in 2D and 3D cases in this chapter. More information on the general procedure and detailed formulation on FEM can be easily found in some classic textbooks [21, 22]. In this chapter, some numerical results generated from the EM model are compared with existing data in literature for the purpose of validation of the code. 4.1 Basis Functions for N ode-based FEM The finite element method is used for modeling a wide class of problems by breaking up the computational domain into elements of simple shapes. Suitable interpolation polynomials (commonly referred to as shape or basis functions) are used to approximate the unknown function within each element. Once the shape functions are chosen, it is possible to program a computer to solve complicated geometries by specifying the basis functions, geometry, and forcing functions. Node-based shape functions have been used extensively in civil and mechanical engineering applications, for example in heat transfer problems. In node-based finite elements, the form of the sought function in the element is controlled by the function values at its nodes. The approximate function can then be expressed as a linear combination of basis functions weighted by the nodal coefficients. If 26 the function values uie at the nodes are taken as nodal variables, then the approximating function for a two-dimensional element e with p nodes has the form P ue(x,y) = Zquie (x,y) (4.1) i=1 Since the expression (4.1) must be valid for any nodal variable uie , the basis function N i6 (x, y) must be unity at node i and zero for all remaining nodes within the element. 4.1.1 Triangular Elements First-order triangular elements are popular in 2D case because they can model arbitrary geometries very easily. The unknown function it within each element is approximated as ue(x,y) 2 ae + bex + cey (4.2) where ae, be, and ce are constant coefficients to be determined, and e is the element number. For a linear triangular element, there are three nodes located at the vertices of the triangle. Assume that the nodes are numbered counterclockwise (to satisfy a right hand rule) by numbers 1, 2, and 3, with the corresponding values of u denoted by ule , u§ , u§ , respectively. Enforcing (4.2) at the three nodes, the following expressions can be obtained uf 2 ae +bexle +cey1e e_e ee ee u2—a +b x2+c y2 27 e_e ee ee u3—a +b x3+c y3 e where x J and y; (j = 1, 2, 3) denote the coordinate values of the jth node in the em element. Solving for the constant coefficients ae , be , and ce in terms of ue. and J substituting them back into (4.2) yields 3 ue(x, y) = 2N? (x, y)uj (4.3) j=1 where N13. (x, y) are the interpolation or basis functions given by 1 e —.._ e. c: e. -= Nj(x,y)—2Ae(aj+bjx+cjy) 1 1,2,3 (4.4) inwhich e_ e e e e, e_ e e, e_ e e “1‘x2y3“y2x3’ b1 ‘yz‘yr €1‘x3‘x2 e_ee_ee, e_ e_ e, e_ e_e “2—x3y1 y3x1’ b2"y3 Y1, C241 x3 e_ee_ee, e_ e_ e, e_ e_e “3"‘1y2 y1x2’ b3‘y1 yz’ C3“"2 x1 and 1 x13 yle 1 1 Ae=51 x; yze =—2-(blec2e—b§cf)=area ofthe eth element. e e 1 x3 y3 It can be easily shown that these basis functions have the property 1 i=j. e e e _ ..: Ni(xj,yj)—5,J {0 we)“, (4.5) 28 As a result, ue in (4.3) reduces to its nodal value uie at node i. Another important feature of N; (x, y) is that it vanishes when the observation point (x,y) is on the element side . .th . . opposrte the 1 node. Therefore, the value of ue at an element s1de rs not related to the value of u at the opposite node, but rather, it is determined by the values at the two endpoints of its associated side. This important feature guarantees continuity of the solution across the element sides. A very useful integration formula in terms of the area coordinates — N1, N 2, N3 — over a triangular domain is given by [23] IIAeiNfIINiinINilndW- M"! Me (4.6) —(l+m+n+2)! 4.1.2 Tetrahedral Elements The three-dimensional analogue of a two-dimensional triangle is a tetrahedron (four-faced element). Once again, special coordinates, called volume coordinates or simplex coordinates, can be used to simplify the derivation of shape functions. Let us consider the tetrahedral element illustrated in Fig. 4.1. Within the element, the unknown function u can be approximated as (for first-order elements) ue(x,y,z)=ae +bex+cey+dez (4.7) 29 3 Figure 4.1 Linear tetrahedral element. The coefficients are , be, ce ,and de can be determined by enforcing (4.7) at the four nodes of the element assuming a given value at the vertices. Thus, denoting the .th . . . value of u at the 1 node as ue- , the followmg expressrons can be obtained uf = ae + bexf + ceyf + .1er e_e ee ee ee uz—a +b x2+c y2+d 22 e_e ee ee ee u3—a +b x3+c y3+d Z3 e_e ee ee ee u4—a +b x4+c y4+d 24 from which we obtain 2 Q § >1 wmwmwmwm R t H V: hm hm hm hm O\ V s: O\ V m __ \t «Hm—Web‘s: ‘< H mmNmNmNm ‘< N N N N 30 1 1 1 1 1 ul u2 “3 ”4 1 e _ — e e e e e e e e b — 6V9 yle y; y; y: — 6Ve(b1ul +b.2“2 +b3“3 +194144) e e e e 21 22 Z3 Z4 1 l 1 1 e 6 e e 1 x1 x2 x3 x4 1 “e = 6,. ur a; u; u: 2 We «5.5.1.3.; +c§u§ +czuz> e e e e Z1 Z2 Z3 Z4 1 1 1 1 e e e e e — ____ e e e e e e e e d - 6V3 yle y; y; y: — 6Ve(d1u1 +d2u2 +d3u3 +“14%) e e e e "1 ”2 ”3 “4 where 1 1 l 1 e e e e e 1 x1 x2 x3 x4 V =__ e e e e = volume of the element. 6 yl Y2 y3 y4 e e e e 21 22 Z3 Z4 The coefficients a; , b? , ce. 1 j , and d; can be determined from expansion of the determinants. Substituting the expressions for ae , be , Ca ,and de back into (4.7), (4.7) becomes 4 e e e u (xaysz) = Zle(xayaz)uj (4'8) 1: where the interpolation functions N j (x, y, z) are given by 31 e _ e e e e Nj(x,y,z)— (aj+bjx+cjy+djz). (4.9) 61/9 Similar to the two-dimensional case, it can be shown that the interpolation functions have the property 1 '= '. ' 1 (4.10) Ni(xjayjazj) 5:} {0 iij. and furthermore that N 37 (x, y, z) vanishes when the observation point is on the surface of . .th . . . . the tetrahedron oppos1te the j node. As a result, facral mter-element contmurty of the interpolated function is guaranteed. Similarly to triangular elements, volume coordinates greatly simplify integration over tetrahedral elements. A useful formula for integrating over the volume of a tetrahedron is [23] IILeINfIkiNiiiNseiniNirdV= mm”! 6V8 (4.11) (k+I+m+n+3)! 4.2 Basis Functions for Edge-based FEM In electromagnetics, serious problems are encountered when node-based elements are employed to represent vector electric or magnetic fields. First, spurious modes are observed when modeling cavity problems using node-based elements [24], which is generally attributed to lack of enforcement of the divergence condition. Second, it is inconvenient to impose boundary conditions at material interfaces as well as at conducting surfaces. Third, it is difficult to represent conducting and dielectric edges and comers due to field singularities associated with these structures [25]. 32 Edge-based finite elements, whose degrees of freedom are associated with the edges and the faces of the finite element mesh, have been shown to be free of the above shortcomings. Edge basis firnctions were described by Whitney [26] over 35 years ago and have been revived by Bossavit and Verite [27], Nedelec [28], and Hano [29] in the recent past. It was Nedelec’s landmark paper [28] that laid down the guidelines for constructing finite element basis functions that span a curl conforming space with degrees of freedom associated with the edges, faces, and elements of a finite element mesh. 4.2.1 Triangular Elements Since the edges of an arbitrary triangular element are not parallel to the x- or y- axis, it is not easy to guess the form of the edge-based basis function by inspection. Therefore, the edge-based basis for a triangular element will be expressed in terms of its area coordinates, L‘f , L; , L; . These are the Whitney elements. If the local edge numbers are defined according to Table 4.1, then edge bases for a triangular element are defined as e _ e __ e e e e e . ._ Wi —Wij.—liJ.(LiVLj—LjVLi), 1,]-1,2,3 (4.12) where Wie denotes the basis function for the ith edge of the eth element and 1,)- is the length of the edge formed by nodes i and j. The vector field inside the triangular element can, therefore, be expanded as 3 e e e E : ZEI Wi (4.13) i=1 where E 1.6 denotes the tangential field along the ith edge. It can be easily shown that the edge-based functions defined in (4.12) have the following properties within the element 33 e_ V-Wij—O e _ e e e VxWi]. —21ijVLi xVLJ- If e] is the unit vector pointing from node 1 to node 2, then é] oVLl =—1/11 and él -VL2 = 1/11 . Since L1 is a linear function that varies from unity at node 1 and zero at node2 and L 2 is unity at node 2 and zero at node 1, it is clear that éyw§=q+g=1 along the entire length of edge 1. This implies that W1":2 has a constant tangential component along edge 1. Moreover, since LT vanishes along edge 2, VLf is normal to edge 2, L; vanishes along edge 3, and VB; is normal to edge 3, W182 has no tangential component along these edges. Similar observations apply to W333 and W381 . Thus, tangential continuity is preserved across inter—element boundaries but normal continuity is not maintained. These are the conditions necessary for vector EM problems. Table 4.] Edge numbering for a triangular element Edge No. Node i 1 Node 1'; 1 1 2 2 2 3 3 3 1 4.2.2 Tetrahedral Elements Tetrahedra are, by far, the most popular element shapes to be employed for three- dimensional applications. This is because the tetrahedral element is the simplest 34 tessellation shape capable of modeling arbitrary three-dimensional geometries and is also well suited for automatic mesh generation. The derivation of shape firnctions for these elements follow the same pattern as that for triangular vector basis functions. Consider the tetrahedron shown in Fig. 4.1 and define the edge numbers according to Table 4.2, it is clear that e_ e_e e e e e .._ w, ~W1j —lij(LiVLj —LjVLi), 1,)- 1,2, 3,4 (4.14) where LT, Li, L; , L: are the basis functions for the node-based tetrahedral elements. The vector field within the element can then be expanded as 6 e _ e e E _ ZEiWi (4.15) i =1 where the coefficients E 1.9 represent the average value of the field along the ith edge of th the e element. Table 4.2 Edge definition for a tetrahedral element Edge No. Node i 1 Node i 2 1 l 2 2 1 3 3 1 4 4 2 3 5 4 2 6 3 4 An elegant explanation of the physical character of the edge-based interpolation function is given by Bossavit [30]. Let us consider edge number 1 connecting nodes 1 and 2 in Fig. 4.1. Since VL': is orthogonal to facet {134} and VB? is orthogonal to facet 35 {234}, the field turns around the axis 3-4 and is normal to planes containing nodes 3 and 4. The field thus has only tangential continuity across element faces. Edge elements can also be described as Whitney elements of degree one and can be broadly classified as belonging to the H 0(curl) space. When the edge-based basis functions are employed for a three—dimensional finite element discretization of a vector wave equation, after applying Galerkin testing and the vector identity, the resulting elemental matrices contain the following two integrals [22]: e E; =[[L8(wai)-(waj.)dV (4.16) e _ e e Ft]. _m/ewi -WjdV. (4.17) These two integrals can be evaluated analytically for tetrahedral element. If we adopt the expression given in (4.14), the curl of Wig is given by V x wf = 21f VL‘? x VL‘? I1 12 21.8 =——'—[(c?d.e —d.ec.e japan)? 4rd? )y+[b?c? —c.eb? )2] (6V"’)2 '1 '2 '1 ’2 '1 '2 '1 '2 '1 I2 11 12 (4.18) where bf , c? and die are defined in Section 4.1.2. Substituting (4.18) into (4.16), (4.16) I 9 becomes 8 e e git-=4” ’1'V [(csd? _d€c€)(ce.de we? ) J (6Ve)4 ’1 ’2 ‘1 ’2 11 12 11 12 4.1.86? _6?d.e)(debe 416.18. ]+(6.ec.e —c?b?)(b€c€ -ce.be. )] '1 ’2 ll '2 1] 12 11 12 ll '2 ll '2 11 12 1l 12 (4.19) 36 For the integral in (4.17), using the expression in (4.14), the dot product is given by e e e.e_i1L+LeeHeeH wi wj —(6Ve)2[Lei1 lefizsz tel szfizle i2+Lj1f1112 Lizszf‘lll] (4.20) where 11'j =bi 8b; + cie c: + d ede Substituting this into (4. 17) and utilizing the general integration formula given in (4.11), each entry in the elemental matrix can be obtained analytically, as shown in Appendix for completeness. 4.3 Code Validation with Waveguide Applications In chapter 5, the multiphysics model for polymer processing will be studied, where the heat transfer model is implemented with node-based elements while the electromagnetic model is implemented with edge-based elements. Before hand, the EM code is validated for a couple waveguide applications in this section. These problems have solution data in literature and are sufficiently general as to validate the code for this thesis. Regarding the mode-matching boundary condition on the two ports of the waveguide, the formulations in [22] are used and the details are cited below for completeness of the section. The geometry of the problem is illustrated in Fig. 4.2, whose specific configuration consists of a dielectric and/or metallic obstacle placed in a rectangular waveguide. As usual, analysis of this problem amounts to solving Maxwell’s equations in the region of interest, which for this problem is the region bounded by the waveguide walls. Since the waveguide under consideration is assumed to be infinitely long, the region of interest extends to infinity. To apply the finite element method to such an 37 unbounded region, it is first necessary to truncate the infinite region to a finite region. This can be accomplished by placing fictitious planes on each of the two sides of the obstacles. These two fictitious planes are denoted by 51 and S2, respectively, in Fig. 4.2. I[ ,, LIL—+1 Figure 4.2 Rectangular waveguide loaded with an obstacle [22]. To define the boundary-value problem uniquely for the region bounded by S1 , $2 and the waveguide walls, it is necessary to prescribe a boundary condition for each of 51 and $2. For this, let us assume that the waveguide is operating at a frequency at which only the dominant mode, TEIO, can propagate without attenuation. This is a practical assumption. Further, let us assume that S1 and 52 are placed sufficiently far from the obstacle so that higher-order modes excited by the obstacle die out before they reach S1 38 and 52- Thus, at S 1 the total field can be expressed as a superposition of an incident TE10 wave and a reflected TElO wave. That is, E(x, y, z) = Einc (x, y, z) + Eref (x, y, z) = E0810(x, y)e—jk2102 + RE0e10(x, y)ejk2102 (4.21) where E 0 denotes the magnitude of the incident wave, R denotes the reflection coefficient, and em and kzlo are given by 2 A . 75$ 2 ”- e10(w)= ysm? kle = k0 {2) with a being the width of the waveguide. From (4.21), the following relationship is obtained a x (V x E): —2 x (V x E): —jkzloEi"C + jkleEref = J'kzloE - ijzroEmc (422) which can also be written as hx(VxE)+7hx(fixE)=Uinc (4.23) with 7 = 1k210a Umc = -21'kz10Ei"c- Similarly, the field at 52 can be expressed as E(x, y, z) =E’m'” (x, y, z) = TE e (x, y)e_jk2102 (4.24) 0 10 where T denotes the transmission coefficient. From this, the following condition is obtained 39 x(VxE)+yfix(fixE)=0 (4.25) valid for the field at $2. The relations (4.23) and (4.25) can be regarded as the boundary conditions applicable at S1 and S2, respectively. These, together with the well-known governing differential equation for the field within the region of interest and the boundary condition at the waveguide walls, can be solved uniquely using the finite element method. In accordance with the generalized variational principle, the equivalent variational problem is given by {67703) = O (4.26) 1'1 x E : 0 on waveguide wall where F(E) Z—III/[i' (V)xE )(VxE)—kggrE-E]dV + Hg][§(axn).(axn)+n.uinc]25+ [L2[22’—(th)-(fixE)]dS (4.27) For the numerical discretization of (4.27), first, the volume V is subdivided into small volume elements and within each element the field is expressed as (4.15). For a discretization of the surface integrals, the planar surfaces S1 and S2 are subdivided into small surface elements. Note that the surface and volume discretizations must be compatible; that is, an element edge resulting from the volume discretization must also be an element edge in the surface discretization. Actually, this is always the case because 40 once the volume is discretized, its surface is also discretized. Within each surface element, we may approximate the field as A S "S S S S S S T S an = ZSiEi = = (4.28) i=1 where the superscript 5 denotes the surface element number, n3 represents the number of edges comprising the element, and the S IS denote the vector basis functions. Obviously, S; = fixWiS , where Wis are the vector basis functions that have a unit tangential component at edge i. Substituting (4.15) and (4.28) into (4.27), we obtain 1 M T 1 MS T MS] T bile} #:1547214} Issllssl- 215:} lbs} («29> 2e=l e=l e=l where M is the total number of volume elements, MS is the number of surface elements on 51 and 52, and Ms] is the number of surface elements on S 1 only. Also, [1.4: lite [i.{wwel {Wei 3.3.4.4121]... mm [3‘]: [[5, ;ss .SSdS (4.31) {65 j: [[5, 5" .(U’"0 x mas. (4.32) Carrying out the summation, replacing the local edge numbers with the corresponding global edge numbers, and applying the Ritz procedure, we obtain the system of equations given by [K]{E} = {b} (4.33) 41 where [K] is assembled from [K8] and [3"], and {b} is assembled from {b5}. This system can be solved, after imposing the Dirichlet boundary condition to zero out the edge fields coincident with the waveguide wall, for the edge fields inside Vas well as on S 1 and $2. Once (4.33) is solved for the fields at S 1 and Sz, the reflection and transmission coefficients can be obtained from (4.21) and (4.24), respectively. Specifically, from (4.21), the reflection coefficient can be obtained as _ E(xayx) - Eo‘?10(x,y)e_jkz102 15081006 10842102 R 2-2, (4.34) and from (4.24), the transmission coefficient can be obtained as E(x, y, 2) E0810 (x. y)e"jkz‘oz T: 2= 22 (4-35) Ideally, the reflection and transmission coefficients are not functions of (x,y). However, in a numerical solution, when S 1 and S2 are not placed sufficiently far away from the obstacle one finds that they vary slightly with (x,y). This variation becomes higher when S 1 and 52 are brought closer to the obstacle. This observation provides a good check on whether S1 and 52 are placed far enough from the obstacle. Better accuracy for the reflection and transmission coefficients can be obtained by utilizing the orthogonality properties of waveguide modes. Since the dominant mode is orthogonal to all other modes, we may take a scalar product of (4.21) and (4.24) with em and integrate over the cross section of the waveguide. Doing this, the influence of higher- 42 order modes excited by the obstacle on the calculation of the reflection and transmission coefficients can be effectively removed. The new expressions for R and T then become -'k z R=—— Ex, 1,2 -e x, dS—e 121021 4.36 0on [[1 < i 1) 10( y) ( > 2617621022 =———— Ex, ,z -e x, as. 4.37 0on [[52 ( y 2) 10( y) ( ) Another good check for the solution accuracy is based on the conservation of power: If the waveguide and obstacle are lossless, the equation |R|2+|7]2 = 1 must hold exactly. 4.3.1 Loaded Rectangular Waveguide Rectangular waveguides are commonly used to facilitate material characterization since they have good operational bandwidth for single-mode propagation, reasonably low attenuation, and good mode stability for the fundamental mode of propagation [17]. In addition, for various bands of operation, appropriate standard waveguides are ubiquitous in a well-equipped electromagnetic test facility. Consider a rectangular waveguide with a TE“) mode wave incident on a sample at z = 0 as shown in Figure 4.3. The waveguide has fixed dimensions for width —a/ 2 S x S a/ 2 and height y = b, where the width 0 is in general greater than twice the height b. The length L of the sample waveguide region (S) may vary. The sample, which has one fixed dimension, height y = b, and two varying dimensions, width —d/ 2 S x S d/ 2 and length 2 = L , will be centered in the waveguide at x=0. 43 x=a/2 y=b l ——> ‘— i =-a/2 j y= x=-d/ 2 x=d/ 2 F0 z=L x=-a/2 x=a/2 Top view Cross-sectional view Figure 4.3 Loaded rectangular waveguide. Andrew Bogle’s thesis [31] demonstrated, through the use of a modal-analysis technique, how using a partially filled rectangular waveguide cross-section allows for better transmission responses to extract the complex constitutive parameters. In this work, before testing any samples experimentally, a few simple validation tests were done to verify the extraction process. The first test involves reversing the extraction process to generate scattering parameters for the special case of a completely filled waveguide cross-section. This special case allows the generated scattering parameters from mode- matching method to be compared with results for well known transmission line theory and experimentally measured data. Here, in order to test the numerical code developed with finite element method for electromagnetic model, the completely filled waveguide case is studied with the FEM code and is compared with mode-matching results from [31]. The configuration studied consists of the standard X-band rectangular waveguide (a = 0.9 inch and b = 0.4 inch), and a sample of lossless, non-magnetic acrylic with effective permittivity of 2.5 that fills the cross-section of the waveguide with d = a. The total length of the waveguide is 3cm while the length of the sample is L = 0.75 cm. The FEM with edge-based tetrahedral elements described above in this section is used to calculate the electric field as well as the reflection and transmission coefficients according to (4.34) and (4.35). The results are compared with data from [31] using the mode-matching method for the frequency range of 8 — 12 GHz and shown in Fig. 4.4. It can be seen that a good agreement is achieved. S - Parameters 1' _._,fi_‘rii, ’—* "—“‘ r r 1 f 0.9 -------------------------- i | ...... i ' I I I 0.8 -—----—-- -— ————————————— 4'. : ._ ..... -811(FEM): 0.7 -------------------- 14 0.6 -------- ,—-.~ .__J—.—S11(MM) l? (as; 1 ' - . ' ‘ 5 ;_ 0.5 " . -~ ------~----.' (D!- 1 z o . - . I 0.4 1 : - s._ ............ 0.3 ------- g ------------------------------ .. 0'2 ——————— L —————— 1 ——————————— 1 ______ 1 J: —————— I- '\ : : : : ' " 0.1»----—-i— ~~~~~~~ 'r—————— ------ I ------ i —————— 3 0 l_ _ ..__-_.1.____. -_ __. l__ _ _ _._ _._¥_ ___IL_ _- __,,I _._. _ , ._I_S _Sl _. _ . -.I 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 Frequency (Hz) 10 Figure 4.4 S-parameters for the completely filled waveguide. 45 4.3.2 Waveguide Bandpass Filters The high unloaded Q-factor of TEon-mode in cylindrical cavity is very attractive for the realization of low-loss narrow-band filters at millimeter wavelengths. Direct coupled filters employing this kind of resonators can be designed to have a very selective passband with steep skirts, together with a very low midband insertion loss. One further advantage is the higher midband pulse power capacity, when compared to other waveguide filters. In [32], a mode-matching technique was presented to analyze accurately this kind of filters. That procedure accurately took into account the effects of the thick coupling apertures, the irises angular offset 20, the spurious responses and the higher mode interaction between adjacent resonators, overcoming the limitations of available approximate models [33], [34]. The main building block of the TEon-mode bandpass filter is the cylindrical resonator shown in Fig. 4.5, coupled both with the external rectangular waveguides and with the adjacent cavities through two rectangular apertures (coupling irises) which are placed on the cavity sidewall with an angular offset equal to 20. In this part, the FEM code developed herein is checked with the results from mode-matching method [32] for a simple case where 20 = 180°. 46 Resonator '| External Waveguide q. . Coupling Irises Figure 4.5 Schematic view of the TEon-mode bandpass filter [32]. Ii The studied configuration has the cavity with diameter D = 34 mm and height h = 20.5 mm, coupled by two identical irises with dimensions a,- = 8 mm, bi = 4 mm and ti = 3 mm to an external R120 waveguide ( aw] = aw; = 19.05 mm, and bw1 = bw2 = 9.525 mm) forming an angle 20 = 180°. For the FEM code, the geometry and mesh information are extracted from FEMLAB [35], as shown in Fig. 4.6. Calculated values of S21 using the FEM code is shown in Fig. 4.7, while the measured and computed S21 with mode- matching from [32] are shown in Fig. 4.8 for a wide frequency band. It can be seen that a good agreement with each other is obtained. 47 Figure 4.6 Meshed geometry of the bandpass filter. 48 S - Parameters 1.5 —-—--———+—————--4--——-—-— __--—-____._.— lIIIIIl 1.3 1.4 1.2 Frequency (Hz) I. 1.1 0.9 x 1010 Figure 4.7 Calculated magnitude of 821 for the bandpass filter with FEM. 49 ° 1 1 ""rr 1 ' u 1 TEm E 1'15er ‘ TEorr ' 1 . T3311 I I T . '. ’ -3] ~— ~_-o I Computed ’3. -40 I. ' ' U r r" I °°°°° I m T I ' I ' I I ; ' I Measured _w b .. . 2 . , . _ II ’ I i f § ‘ m i E -_l l _ I w“I I - _l__ 8 9 10 ‘II 12 13 14 15 Frequency [GHz] Figure 4.8 Magnitude of S2] for the bandpass filter from [32]. 50 CHAPTER 5 CETK MULTIPHYSICS MODELING 5.1 Overview Coupled EM-thermal problems have been studied for various applications [36-46], for example, microwave drying of capillary porous materials [36], microwave epoxy curing [37], microwave processing of foods [38] [39], microwave sintering of ceramics [40], microwave heating of laminar material [41], and RF heating for non-ablative cutaneous therapy [42]. The multiphysics numerical models in these papers were implemented using either FDTD or FEM to simulate the time-marching process, considering both electromagnetic and heat transfer effects. Among these efforts, articles [39] and [40] include the consideration of the sample dielectric property variation as function of temperature during heating process, and update the material property from time to time when necessary. Article [37] reported a two-dimensional EM-thermal modeling of epoxy curing process for two operating frequencies, but did not include the reaction heat due to the thermoset curing in the energy balance equation, which is very significant after the curing process starts. In this chapter, a self-consistent 3D marching-in-time multiphysics model is presented, which includes electromagnetic field distribution, microwave power absorption, heat transfer, and polymer curing kinetics. Temperature-dependent permittivity and curing kinetics for diglycidyl ether of bisphenol A (DGEBA) / diaminodiphenyl sulfone (DDS) based on experimental data are explicitly included in the model. Edge-based FEM is implemented for electromagnetic model to ensure the tangential continuity of electric field and divergence-free condition for source free region, while node-based FEM is used in thermal model to solve for the temperature distribution. 51 The numerical results can be used to determine the time-dependent temperature distribution and curing profile across the polymer sample, as well as the electromagnetic field distribution within the cavity applicator. With the help of this numerical model, robust control strategies can be developed for the polymer curing process. The numerical results are compared with the measured data and a good agreement is achieved. 5.2 Energy Balance Equation During the epoxy thermoset curing process, the energy balance for the material being processed within the applicator is governed by the following equation 6T - a pCp—a—t-=P+V-[k(T)-VT]+p—a%(—AH,) (5.1) where p is density (kg/m3), C17 is heat capacity (J/kg/K), T is temperature (K), t is time (second), It is the anisotropic thermal conductivity tensor (W/K/m), a is exten-of-cure, and — AH , is reaction heat (J/kg). On the right hand side in (5.1), P is microwave power absorption rate (W/m3) inside the epoxy sample, and is expressed in terms of the electric field strength and the perrnittvity as below: 1 ,, It P = Ewsogr (a,T)E - E (5.2) where E is the electric field strength inside the epoxy material (V/m), wis the radial frequency (rad/sec), £0 =8.854><10_12 farad/m is the free-space permittivity. 5" is the effective relative loss factor and is the function of local temperature and extent-of-cure for epoxy sample. More details will be discussed in section 5.4. The electric field is 52 determined by solving the vector wave equation using finite element method, and will be presented in details in section 5.5. The term V-[k(T)-VT] represents the conduction heat transfer, which is the primary heat transfer mechanism, inside the epoxy sample. The free convection boundary condition is applied on all surfaces, and will be discussed in detail in section 5.6. Since the epoxy thermoset curing process is exothermic, the chemical reaction heat is included in the energy balance equation (5.1), which is the term paa—a(— AH r). Each thermoset system has a particular extent-of-cure vs. time characteristic that is also dependent on temperature. More details on the curing kinetics model will be presented in section 5.3. Hence, the energy balance equation (5.1) readily demonstrates the interrelationship between electromagnetic, composition, thermal, and kinetic modeling. The critical aspect of the modeling approach used herein relies on the fact that electromagnetic interactions occur much faster, 0(ns), than thermal effects, 0(ns), or the kinetic/composition effects, 0(s). Hence, the electromagnetic analysis can be performed with an efficient steady-state analysis method whose results are used in solving the energy balance equation for heat transfer problem. Specifically, in the heat transfer model, after the node-based FEM for spatial variations is performed, a first-order ordinary differential. equation (ODE) system is obtained and the finite difference scheme for temporal variations will be used for the time-marching process. The dielectric property of the epoxy sample is updated with the new temperature and extent-of-cure, then the electric field will be solved again to get the new power absorption rate and solve for the heat transfer equation. This process is repeated and the temperature profile can be 53 obtained across the sample for any given time. More details on the procedure of this process model will be given in section 5 .6. 5.3 Curing Kinetics For polymer curing, existing reaction kinetics models include an n'h-order reaction model [6-7] and an autocatalyzed reaction model [8-10]. The curing reactions for an epoxy are characteristic of autocatalyzed reactions. The reactions between amines and epoxide are autocatalyzed by the hydroxide group formed in the reactions. The generalized form of the autocatalyzed reaction kinetics for epoxy curing was derived from the reaction kinetic mechanism [1 l]. The extent of reaction can be measured experimentally as a function of time and temperature. The reaction constants and parameters can be determined with regression of experimental data. Reaction kinetics for a simple epoxy resin, diglycidyl ether of bisphenol A (DGEBA) / diaminodiphenyl sulfone (DDS), at 145°C, 165°C, and 185°C, has been studied at Michigan State University [47]. The reaction heat, — AHr , was 413.9i5.1 J/g as determined with Differential Scanning Calorimetry (DSC). The curing reactions for an epoxy are characteristic of autocatalyzed reactions. The reactions between amines and epoxide are autocatalyzed by the hydroxide group formed in the reactions. The following phenomenological model was used to describe the autocatalyzed reaction for a stoichiometric reactant mixture: id?- 2 (k1 (T) + k2(T)am)(1— a)" (5.3) 54 where k1 is the non-catalytic polymerization reaction rate constant, k2 is the autocatalytic polymerization reaction rate constant, m is the autocatalytic polymerization reaction order, and n is the non-catalytic polymerization reaction order. The epoxy was cured with microwave power in a cylindrical single mode cavity. The extent of reaction was measured experimentally as a function of time and temperature. The reaction constants and parameters were determined with regression of the experimental data and presented in Table 5.1 [47]. Using these parameters, the extent of reaction can be regenerated with the kinetics model. Fig. 5.1 shows the experimental data (markers) and the calculated values (lines) [47]. The kinetic model fits the experimental data very well. Fig. 5.2 shows the calculated time derivative of extent-of- cure (do/dt) as a function of extent-of-cure for different heating temperature, which will be used in the numerical model on the right hand side (RHS) considering reaction heat. In the code, the extent-of-cure distribution across the sample can be obtained by the integration of the time derivative of extent-of-cure for the timestep locally for each tetrahedral element. The time derivative of extent-of-cure will be updated during curing process for each element with the local heating temperature and extent-of-cure. Table 5.1 Epoxy curing reaction constants and parameters Microwave Microwave Microwave Curing at Curing at Curing at 145°C 165°C 185°C Non-catalytic reaction order n 1.57 Autocatalytic reaction order m 0.94 Non-catalytic reaction rate constant k1 0.016 0.040 0.089 Autocatalytic reaction rate constant k2 0.084 0.16 0.24 55 Extent of Cure T=l85 30 Time (min) 40 60 Figure 5.1 Extent-of-cure vs. time and temperature for DGEBA/DDS (markers represent experimental data and lines represent calculated data from the model). x 10'3 time derivative of extent-of-cure (da/dt) 1.8, I" r I I 1 "' l | | I 777 r T i i 7%, ..... ~-~. . . . 145 °Cj ‘ 7 .’-7 7 ‘.~77I 7 7 j 777777 1 "61 x’ 1 f 4,‘ L : ----- 165 °cI I I ". i 1 I _ ...... o 1_4 1 . I .44777‘1". i 44 SI 185 CII I I I ‘0 I l : : ' ~- 1 ‘ 1 2 1 7777777 I 7 ?°‘77 7 7 7 71 7 ' I I ‘7 I p y I I ‘7 I '1: l : 1 I l x 1 : 8 11L——-——>-————|- ————— o————+ ———————— ‘0‘ ————— —1 ——————————————— —~ 8 1 1 1 ; l '\. I ‘19, 1 I, —-1‘—-_-1~~-~ 1 \.\ l g 0.8—~;"',--A--‘---4-%”3‘3; ---------- \ul -------------- - \ j I 1 1 ~\ 1 \. 1 8 ‘1’ : 1 \, ‘1 extent-of—cure (a) i1,7 .7 Figure 5.2 Time derivative of extent-of-cure vs. extent-of-cure. 5.4 Dielectric Properties The work on assessing the dielectric properties of the curing systems of DGEBA epoxy resin with the curing agent DDS as a function of temperature and extent-of-cure at 2.45 GHz has been primarily conducted by Dr. Liming long at MSU [20, 48]. Generally, the real and the imaginary parts of the complex permittivity increase as the temperature increases and decrease as the extent-of-cure increases. The two reasons accounting for the evolution of the dielectric quantities are a decrease in the number of the dipolar groups in the reactants and an increase in the viscosity. The Davidson-Cole model can be used to describe the experimental data. The details will be discussed below in this section. The complex permittivity of materials is expressed as follows: 8* = 6' -—je" = 86(4— jet) (5.4) where s' is the dielectric constant, s" is the dielectric loss factor, so (=8.85 x 10'12 F/m) is the free-space permittivity, s'r is the relative dielectric constant, and s"r is the relative dielectric loss factor. Although s'r is referred to as the dielectric constant in many papers, it is a function of temperature, pressure, humidity, and other conditions. For convenience, specifically, in this section hereafter, s' refers to s'.r while 8" refers to s",. The epoxy resin used in this study was DGEBA and the curing agent is DDS. The chemical structures and properties of DGEBA and DDS are shown in Table 5 .2. Stoichiometric DGEBA and DDS were mixed before each experiment. A single- frequency microwave processing and diagnostic system was used to heat the epoxy materials and to assess the dielectric properties. The heating mode was TM 012 and the frequency was 2.45 GHz. The curing temperature was 145°C. After heating, the single frequency microwave curing system was switched to become a low-power swept 57 frequency diagnostic system. Measurements of temperature and dielectric properties using the swept frequency method were made during free convective cooling of the samples. The cooled samples were analyzed with a Differential Scanning Calorimeter (DSC) to determine the residual heat of reaction per gram and, thus, the extents of cure. More details on the experimental system, sample preparation, and measurement procedures can be found in [48], or section 5.7 in this chapter. Table 5.2 Properties of the epoxy resin and the curing agent Name Chemical Structure Epoxy/amine Density equivalent at 25°C weight (g/ml) CH, DGEBA 0H0 0 171-175 1.16 O |>—‘/ CH, HZN NH2 76 62 1.33 The measured dielectric constant and dielectric loss factor of reacting DDS 203:0 O DGEBA/DDS epoxy resins at different temperatures for different extent-of-cure states are shown in Figure 5.3. The extent-of-cure was calculated from DSC data, which are shown in Table 5.3. Due to the non-uniformity of microwave heating, a distribution of extent-of-cure exists in the samples. From Figure 5.3, it is found that both dielectric constant and dielectric loss factor increase as the temperature increases and decrease as the curing reaction proceeds. 58 o 0% 6 3 ‘ I o 9% D ‘ ° I A 14°/ :1 . . D l o o 0 - ‘ .u '3‘ o + + I o 220/0 “5 1 0.20‘ o + j o 3 003.00 7 + 7 . ' . +37/o ‘ “O I X ' o 4Ifi+ + I+..XX I-42/0 * + 9: X. x x A A A l X64% ‘1‘A. A A A A A I o 3 -: 'T-Tfihffi'T"' "r”“r 1' 1.11... I r _T '1 _g _T 7 A 79/0 20 4O 60 80 100 120 Terrperature (°C) 0.8 A? *7 _SSSS_ __ 5 o 0% _ 7 _ 0.6 '1 Q . D ‘ i D 9% : .’ 6 ° . o 414% 300.4; 0 no A 0 0° +;022% I ’ D ‘ o + ; 1 . g, D A ‘0 + + 7 I + 37% 02 " .33 Q, o +."'e ' . I42% . 7 o + E 4%+oj'I+- .X x xx X XA1X64o/o 741 .Ax XA XA A A A A A : 0.0 “I‘T‘T'TTTf‘r—T—‘W—r‘r—j-fi—j~rw—r~~—T THESE, A 79% 20 40 6O 80 100 120 Temperature (0 C) Figure 5.3 Measured dielectric properties for the DGEBA/DDS. Table 5.3 DSC results of the curing DGEBA/DDS system Reaction Heat (J/ g) 432 387 395 324 249 168 203 72 448 436 415 350 253 299 153 109 419 403 360 362 269 262 58 63 439 370 330 370 294 222 57 64 429 365 363 290 308 297 319 157 Average 433 392 372 339 275 250 158 93 Extent of Cure 0% 9% 14% 22% 37% 42% 64% 79% Standard Error 1% 3% 3% 3% 3% 6% 11% 4% 59 The stoichiometric curing reaction of epoxy resins with amines is illustrated in Figure 5.4 [49, 50]. In the first step, an epoxy group reacts with a primary amine to form a hydroxyl group and a secondary amine. In the second step, the formed secondary amine reacts with another epoxy group to produce a hydroxyl group and a tertiary amine. The progress of the reaction is defined in terms of extent-of—cure. The viscosity of the reacting systems rise as the polymerization goes on and two distinguishable transitions, i.e. the gelation and the vitrification, were observed. Step 1 R—'NH R1 R—NH2 + I>—R, ——> 0 OH Step 2 OH R—NH R1 f—(R + l>—R2 ——> R—N R2 O 1 OH H OH Figure 5.4 DGEBA resin curing mechanism. The unreacted epoxy groups, -O-, and amine groups, -NH2-, within the reaction system are the main driving forces for the apparent combined 7 relaxation observed in this study. The disappearance of the epoxy and amine groups during the reaction is one reason accounting for the changes of the dielectric properties. An increase in the viscosity during the reaction is the other reason. The increasing viscosity of the reacting systems 60 should hinder the mobility of the dipolar groups associated with the relaxation and cause the relaxation time to increase. The Davidson-Cole model [51] can be used to describe the dielectric behavior of DGEBA epoxy resins: 8*-8oo : (£0 _goo) (1+jwr)" , + (so — soo )cos(n 19) n (1 + (4455 (5.5) ,, 7 (so — s00 )sin(n 6) (new)? 8 6 = arc tan(a)r) where t is the relaxation time in second; so is the static frequency dielectric constant; 800 is the high frequency dielectric constant, and n is the shape parameter with a range of 0 S nS l. The dominant y relaxation time should fit the Arrhenius rate law [52]: r = A x 1%] (5.6) where Ea is the activation energy in J/mol, R (=8.314 J/mol-K) is the molar gas constant, T is the temperature in K, and A is the relaxation time in the high temperature limit in second. The values of the parameters, n, sac, and E0, and a relation between (so-so.) and A, can be calculated using equations (5.5) and (5.6) [48], while the values of (so-soc) was taken from the literature [53] and then A was determined. Specifically, all the values of 61 the parameters for the curing DGEBA/DDS system with different extent—of-cure are listed in Table 5.4. Figure 5.5 shows the comparison of the experimental and calculated data of the reacting DGEBA/DDS system, where the curves represent the calculated data from the Davidson-Cole expression and the markers represent experimental data. It can be shown that the data from Davidson-Cole model fit the measured data very well within the temperature range between 20 °C and 130 °C (note that the thermoset curing for DGEBA/DDS starts once the temperature goes beyond about 125 °C). 7 , I o 0% ‘ t" a 91V 6 1 7,07)" 7° ° “ .’ 13/4” ,0” A 14% — /, D //" '16 ’o, 5 4 .3;ng 7/+ o 22% 4” 4.5‘9136/ ’ r! .Jt/x I42% 2 /+//I'/::/” M « W x 64% - 7r- 3 :F:A:1 r—fi' 1+1 rjfi' t—r—r 1 A 79% 20 40 60 80 100 120 Temperature (° C) 0.8 o 0% u 9% 0.6 — « / A 14% . 04 1 ° 0 22% ° ' j + 37% - 420/ 0.2 — ”V ' ° ' W x 64% : Amt-1494’“ “379% 0_0 ll] 1T1 T‘FTTTT‘I’I 11 20 40 60 80 100 120 Temperature (0 C) Figure 5.5 Comparison between the experimental and calculated permittivity of DGEBA/DDS. 62 Table 5.4 Values of the parameters for the curing DGEBA/DDS system Extent of Cure n so a.o (so-s00) Ea (kJ/Mol) A (s) 0% 0.17 9.10 3.29 5.81 110 4.2E-25 9% 0.15 8.67 3.33 5.34 110 1.7E-24 14% 0.14 8.47 3.30 5.16 112 2.1E-24 22% 0.13 8.12 3.43 4.69 107 5.2E-23 37% 0.11 7.45 3.11 4.34 119 8.6E-24 42% 0.10 7.19 2.89 4.30 133 7.1E-25 64% 0.08 6.24 2.93 3.31 129 5.4E-23 79% 0.07 5.57 2.73 2.83 121 7.2E-20 In the numerical code for the process model, the complex permittivity of the epoxy sample will be updated with the local temperature and extent-of-cure as the microwave heating is proceeding. In order to get the values of the complex permittivity that are beyond the temperature range of the measured data and to do the interpolation between different extents of cure, the dielectric constant and loss factor are calculated for the temperature range of 0 °C and 300 °C based on Davidson-Cole model, and shown in Fig. 5 .6 and 5.7, respectively. However, it is found that the dielectric loss factor calculated from Davidson-Cole model is not accurate for relatively high temperature, say, beyond 130°C for DGEBA/DDS in this case. Considering Fig. 5.7, since interpolation is necessary for the values of the dielectric loss factor between different extents of cure in the numerical model, at the high temperature the values of the loss factor will change very abruptly. This is not the real situation as observed from experimental value obtained in [54], where the loss factor changes smoothly. In order to overcome this problem, the values of the dielectric loss factor calculated from Davidson-Cole model are modified such that they follow the same trend as the dielectric constant and are shown in Fig. 5.8. Specifically, after the dielectric loss 63 factor reaches its maximum value for each extent-of-cure, it remains constant at higher temperature. In practice, it is very difficult and sometimes even impossible to obtain accurate measured data at high temperature for the low extent-of-cure, since the curing is actively going on at high temperature and both the dielectric loss factor and extent-of- cure keep changing during the measurement. Therefore, even this modified model may not be true in practice; however, it gives a very good way to do the interpolation and to provide reasonable approximation values for the loss factor, and as will be seen in section 5.8, good results are obtained based on this modified Davidson-Cole model. It is hoped more experimental efforts and theoretical research can be done in the near future to provide more accurate model on the complex permittivity for epoxy resins. Dielectric Constant (s' )-—Davidson-Cole Model 10 1 T ‘1 I I I ---- cure=0°/o I — cure=9% 7 __ 7é 9| -'-"Cure=14% ’,7_ ‘3. l + cure=22% ‘ I, d...“- = 8 I O cure=37% I I” W : o a i" - cure=42% ,’ .1 7y :‘3 II . cure=64% ’, I.’ +++ . ......eeeeeeeeeeeeeee E 7II A cure=79% I ’ +° ~1“I'“mmnm' '4 m I I .4» .0 +‘ 8 l ’ 9’ * e. l.- U 1 [I 0’ +++ .°. I I go it) 6'— ’ .’ If... 0.. 34an ..O .... 8 ” . +++ 7.. 7‘1‘ ...00 E 1’ if”. x", “v". D 5 r 1’ 0’43“ .-°.. 4"“ "0'... T ’ 4’... e'. .+'. '0'... ..e 7 K" ....O. 4 h ‘ xx’:::e eeeeee 7 o... ..‘ta... 9.9.!” .. - - - ' . . . - - ' ' g 1 1 _L l 0 50 100 150 200 250 300 Temperature ( °C) Figure 5.6 Calculated dielectric constant for DGEBA/DDS based on Davidson-Cole model for a wider temperature range. Dielectric Loss Factor (3" )——Davidson-Cole Model 1.4 , i . met—— —-—»—A-__.~~_. ---- cure=0% ! — cure=9% I l 1.2 -----cure=14% a, + cure=22% o cure=37% cure=42% I \ - cure=64% I 0.3 A cure=79% , a Dielectric Loss Factor (8" ) Temperature ( °C) Figure 5.7 Calculated dielectric loss factor for DGEBA/DDS based on Davidson-Cole model for a wider temperature range. 65 Dielectric Loss Factor (3" )—Modified Davidson—Cole Model 1.4 v l l l - - -- cure=0% ‘ — cure=9% 1.2 ----- cure=14% + cure=22% : o cure=37% x cure=42% ‘ - cure=64% 0.8 A cure=79% Dielectric Loss Factor (8" ) , , , ,W, , . 7 , l,,,,, _—._ 0 50 100 150 200 250 300 Temperature ( °C) Figure 5.8 Calculated dielectric loss factor for DGEBA/DDS based on modified Davidson-Cole model for a wider temperature range. 5.5 Electromagnetic Model The edge-based finite element method with tetrahedral elements for the applicator and polymer sample is used since this method offers a high degree of geometrical and electrical fidelity. In particular, the spatial variation in permittivity can readily be accommodated and by using “snapshots-in-time”, time-dependant thermal and composition influences on the local permittivity can be included in the model. Recall that such quasi-stationary analysis approaches are valid since the evolution of thermal and composition effects are relatively slow compared to the rate of electromagnetic field propagation within the applicator. 66 The finite element method is well-suited for resonant cavity simulation since it can incorporate the presence of a sample as well as determine the response of the system to single frequency excitation in an efficient manner. Once the fields are determined throughout the cavity for a given excitation, local power absorption rate can be determined for the heat transfer model. Specifically, the configuration of the applicator loaded with an epoxy sample in the numerical model is shown in Fig. 5.9, with an enlarged view at the epoxy sample on the right side in the figure. The applicator is a coax-fed single-mode cylindrical cavity with the diameter of 15.24 cm and the height of 16.0 cm, such that TMmz mode is excited inside the cavity at 2.45 GHz. The pattern of the total electric field distribution for TMmz mode is shown in Fig. 5.10. The epoxy sample is centered at the middle of the cavity, where the electric field strength is the strongest. The diameter of the epoxy sample is 1.0 cm and the height is 2.7 cm. The sample is held by a Teflon holder with the wall thickness of 0.2 cm and the base thickness of 0.3 cm. The Teflon has the complex permittivity a: = 2.1—j0.0004. ‘ 7mm PointC , F "X'— t 7mm PointB > —}% 6mm PointA F %F 7 mm Teflon _L, ' I 4L— holder Epoxy sample Figure 5.9 Coax-fed single-mode cavity with the epoxy sample loaded. 67 Figure 5.10 Pattern of the total electric field distribution for TMon mode inside the cavity. Regarding the mesh generation for computational purpose by finite element method, SKY/Mesh2 mesh generator [55] was used to get the 2D mesh with triangular elements. The triangular elements are extruded to form prism elements, and the latter elements are subdivided to get tetrahedral elements [56], which provide increased flexibility on geometry. The electric field, inside the cavity and at each point within the epoxy sample being processed, can be computed using edge-based finite element method. The set of finite element equations, derived from the vector wave equations and using Galerkin testing, is given by 21% VxWi -Vij - [—————- 435w,- .wj] dV = -jk0Z0 j'Vwi .J‘de (5.7) 68 where the vector testing fimction is denoted by W,- and the vector expansion function for the total electric field is denoted as W j- Tetrahedral elements with edge-based basis functions, which have been discussed in chapter 4, are used to solve for the electric field. The term on the right hand side (RHS) accounts for the excitation from coaxial cable at the top, in the way demonstrated by [21]. The Dirichlet boundary condition, fixE = 0 where I“: is the surface normal, is applied on the inner conducting surfaces of the cylindrical cavity. Regarding the data storage scheme for the sparse matrices generated from assembling, compressed sparse row (CSR) format [57] was used to save the memory and CPU time. For the solving process, a biconjugate gradient (BiCG) iterative solver [58] was used to solve the linear equations. Once the values of the electric field are obtained, the power absorption rate within the epoxy sample can be obtained with the expression 1 , * P = §w£05,(a,T)E - E (5.8) Since the complex permittivity of the DGEBA/DDS is the function of the local temperature and extent-of-cure, both the dielectric constant and the loss factor keep changing during the thermoset curing process. Therefore, in order to simulate the process accurately, the material property must be updated. Specifically, the complex permittivity of the epoxy sample, in the system equation (5.7) and the power absorption rate expression (5.8), must be updated with an appropriate timestep in the time-marching process. This procedure will be discussed in more details in section 5.6. 69 5.6 Heat Transfer Model In the heat transfer model, node-based FEM with tetrahedral elements is used since the temperature is a scalar variable. The temperature distribution can be obtained by solving the energy balance equation (5.1), which includes microwave power deposition, heat transfer and heat generation due to thermoset curing. After the node-based FEM for spatial variations is performed, a first-order ordinary differential equation (ODE) system is obtained and the finite difference scheme for temporal variations will be used for time- marching process. Recall the energy balance equation for the epoxy sample and Teflon holder mpg:P+V~[k(T)-VT]+p%(—AH,) (5.9) in the numerical model, the parameters of the material properties for DGEBA/DDS epoxy resins and the Teflon holder are listed in Table 5.5. For convenience, there is no difference in symbols between the epoxy resins and the Teflon holder. Table 5.5 Parameters of the material properties in the heat transfer model Parameter Symbol Value Unit DGEBA/DDS Density p 1 160.0 kg/m3 Heat capacity Cp 200.0 J/kg/K Thermal conductivity k 0.24 W/m/K Teflon Density p 2190.0 kg/m3 Heat capacity Cp 995.0 J/kg/K Thermal conductivity k 0.24 W/m/K 70 Although the same mesh from the electromagnetic model is used, only the region consisting of the epoxy resins and the Teflon holder needs to be considered for the heat transfer problem. After applying Galerkin testing (assuming diagonal tensor for thermal conductivity), (5.9) becomes 6N- . aN. . aN- . [[19, 1%{T}+k ’ aN'{T}+kz 621%{T}}V 6x 6x y 73767 6{T} ——a—t—dV + IShNiNJ-{TMS (5.10) 1 ,, 6 = Lawgogr(a,T)E-E*NidV+ [Sham-(15+ Lp—gGAHflNidV + LpCpNiNj In the above equation, the natural air free convection boundary condition is already applied between the Teflon/air or epoxy/air interfaces, as —(IE-VT)-fi=h(T—T,o) (5.11) where I; is the anisotropic thermal conductivity tensor, ii is the surface normal, and T 00 (= 20 °C) is the surrounding air temperature. The convective heat transfer coefficient h is a very important parameter and needs to be considered very carefully. Typical value of h for natural air free convection ranges between 5—25 W/mZ/K. A parametric study on h will be done in section 5 .8 to see the influence of its variation on numerical results. For each tetrahedral element, the elemental matrix is given by e [Cel{%}+[1tem1 E No Yes No ‘ fl: Is ABS(e,* (t +At) - a,‘ (t)) > 6? ‘F 1 Figure 5.11 Flow chart for the coupled electro-, thermo- and kinetic model. 73 5.7 Experimental Setup The microwave diagnostic and processing system at MSU was developed by Dr. Jinder Jow. The details of the system and the single mode perturbation method can be found in his dissertation [54]. A microwave switch between the heating and diagnostic systems was developed by Mr. Gregory Charvat [60]. Here is a brief description of the system. The microwave diagnostic and processing system consists of a microwave external circuit (an energy source, transmission lines, and the coupling probe), diagnostic elements (i.e. an X-Y oscilloscope, power meters, the E-field diagnostic probe, and a temperature sensor device), and the loaded cavity (i.e. the cavity and the loaded material). The schematic illustration of the system is shown in Fig. 5.12 while the system with the equipments is shown in Fig. 5.13. A single-mode circular cylindrical cavity is used to process epoxy resins at a frequency of 2.45 GHz. A fluoroptic temperature sensing device (Luxtron 750), which does not absorb electromagnetic energy and minimally interferes with the electromagnetic field, is used to on-line monitor the material temperature. The fluoroptic probe is protected by a 3 mm O.D. pyrex capillary tube. A cylindrical Teflon sample holder is also required to contain the liquid samples. The cylindrical geometry of both the cavity and loaded materials is selected in order to facilitate diagnosis, modeling, and theoretical analysis. 74 Fluoroptic Temperature S ‘ t Power Meter ensmg SYS em Microwave , Loaded Pt = P1 — Pr . . Material Circulator Coax1al Coupling Cable Probe /_ T Cavity Directional Coupler Pb E Field Diagnostic Reference X Y 20 dB Attenuator Probe Signal Crystal Detector X-Y Osc11105cope Pbo Power Meter Power Meter Figure 5.12 Schematic illustration of the microwave processing and diagnostic system Figure 5.13 Microwave processing and diagnostic system. 75 As mentioned earlier, in this thesis, the epoxy resin is diglycidyl ether of bisphenol A (DGEBA, DER 332 by Dow Chemical) and the curing agent is 3,3- diaminodiphenyl sulfone (DDS by TCI America). All of the materials were used as received without further purification. In preparing neat DGEBA/DDS epoxy resins, stoichiometric DGEBA and DDS (2.79:1 by weight) were mixed in a glass beaker. The mixtures were well stirred by hand in a 130°C oil bath until the DDS was completely dissolved (in approximately five minutes). Finally, the resins were degassed at 0.02 bar at 100°C for five minutes. Before the epoxy sample held by the Teflon holder was loaded in the single-mode cavity, the height of the cavity was adjusted such that the mom mode was excited inside the cavity at 2.45 GHz. The degassed liquid epoxy resins were poured into the Teflon holder. The dimensions of the epoxy resins and the Teflon holder have been described in section 5.5 and in Fig. 5.9, which are the same as here since the numerical model follows the experimental system. Then the Teflon holder with a fluorOptic probe and epoxy resin sample was centrally loaded in the middle of the cavity. Careful and small adjustments need to be done on both the height of the cavity and the position of the Teflon holder such that the sample is located at the place where the electric field is the strongest. The single frequency microwave curing and diagnostic system was used to heat the liquid samples in the Teflon holder. Usually for the purpose of dielectric property measurement, the procedure is as follows. The fresh DGEBA/DDS samples were heated to react at 145°C for specified reaction time periods, e.g. 1, 5, and 20 minutes, with the exception of those for the 0% cured epoxy resin, which were heated to 100°C. Thereafier, 76 the single frequency microwave curing system was switched to a low-power swept frequency diagnostic system by changing the switch position. Measurements of temperature and dielectric properties using the swept frequency method were made during free convective cooling of the samples. The cooled samples were analyzed with a Differential Scanning Calorimeter (DSC) to determine the residual heat of reaction per gram and thus the extents of cure. However, in this study, in order to test the numerical model with the experimental data, only the curing mode (heating mode) was used. The input power was fixed at 5 Watts and the temperature was recorded as it increased. The temperature/time profile generated from the numerical model is compared with the experimental data in the next section. 5.8 Numerical Results vs. Measured Data Fig. 5.14 shows the comparison between the calculated and measured temperature as a function of time during thermoset curing process. For the positions of point A, B, and C, please refer to fig. 5.9. The measured temperature was taken at the center of the epoxy sample. In the numerical model, At = 5 seconds was used as the timestep and h = 15 W/mZ/K was used as the free convective heat transfer coefficient, which was used in the convective boundary condition. It can be seen that the calculated temperature profile is in good agreement with the measured data. Figs. 5.15, 5.16, and 5.17 show the calculated extent-of—cure, dielectric constant, and dielectric loss factor as a function of time for the three points during microwave heating. From Fig. 5.16 and Fig. 5 .17, it can be seen that both the dielectric constant and 77 the dielectric loss factor increase as the temperature increases at the beginning, before the curing starts at 125°C. Once the temperature reaches 125°C, the curing process will begin with chemical reactions. As the extent-of-cure increases, both the dielectric constant and the dielectric loss factor decrease. As shown in the expression (5.2), the power absorption rate is related with the dielectric loss factor. As the dielectric loss factor becomes smaller and smaller, the microwave power energy absorbed by the epoxy sample becomes less and less. From Fig. 5.2, it can be seen that the chemical reaction heat reaches maximum at the extent-of-cure around 20% — 30%, and becomes smaller and smaller after that. Fig. 5.14 shows that the temperature reaches the maximum at 250°C and then starts to cool down. However, the extent-of-cure still increases since the temperature is still above 125°C. The final extent-of-cure of the epoxy resins stops at about 80% and the curing process reaches the end. Fig. 5.18 shows the parametric study on the heat transfer coefficient h with three different values of 5, 15, and 25 (W/m2/1(), and it is clear that the value of h plays a very important role on the accuracy of the numerical model. Since the epoxy thermoset curing process is exothermic, which means that heat energy is released during the chemical reactions, the term p6_a(_ M) in the energy balance equation (5.9) must be used to account for this. In order to see how significant . . . . . . 6a . the chemical reaction heat rs, the srmulatron was run wrthout the term pE—(—AH ) 1n r (5.9), and the result is shown in Fig. 5.19. Apparently, the chemical reaction heat can not be neglected for an accurate model of the process. 78 In this numerical model, the complex permittivity of the DGEBA/DDS was updated as a frmction of both local temperature and extent-of-cure during heating process. Another two cases were also performed to see how important it is to include these variations. In case 1, the complex permittivity is kept constant (5: =4— j0.2 ), which is the initial value at the room temperature before the heating process begins. In case 2, the complex permittivity as the firnction of the local temperature, but not the function of the extent-of-cure, is used. This means the curve with 0% extent-of-cure in the modified Davidson-Cole model is used. All three cases are shown in Fig. 5.20 with the measured data. It is clear that updating the complex permittivity as a function of both temperature and extent-of-cure gives most accurate results. The temperature distributions at the horizontal cross-section layers containing points A, B, and C respectively are specifically recorded in the numerical model for each timestep during the calculation. As the examples to illustrate the temperature variations in both spatial and temporal domains, the snap-shots of the temperature distribution at the cross section of point B are listed in Figs. 21-26 for the heating time of 100, 300, 400, 500, 700, and 1000 seconds. 79 250 T 200 ~~~~~~ ,; —————— I—. -‘ --------- . y150e—"—~---E---I _______________________ : q ‘6 i x’ : :5. i 1’ I E i . a I : , 5100+ ..... ’ _______ ., . . "-1 Q) ' It i i I . I , _ u i a r ‘ A Ii. i i i I I ’ s,» )- I 3 i i i l : °~.‘ 50 W, i ----- simulated data at point A ” ‘ “ ‘3 """" simulated data at point e """""" ' ' i - ------ simulated data at point C I i —-1— measured data at center 0 L l i i L .3 r 0 200 400 600 800 1000 1200 14001500 time(seconds) Figure 5.14 Calculated and measured temperature vs. time for the epoxy sample. 80 14001500 1200 ' _0 ...................... O r . A. _0 _1 ” ..ABC _ttt .m.m.m.} ................. m ammm. e ~.m “a.m “ . aaa“ ttt .aaa. rddd ............................... m www _ 4 mmrla ” UUU _ fmmm U 0 +588 _ -0 W. a _ 2 . _ n . n _ _ . . . H H H ” “H41~.l.i L. _i _. lll . —. . h . r 0 8. 7. 6. 5. A 3. 2 1. 0 0 0 0 0 0 0 0 0 0.9 fin 2:32:26 time(seconds) Figure 5.15 Calculated extent-of-cure vs. time for the epoxy sample. 81 simulated data at point B 3 -----~ simulated data at point A ‘1.-||lr' lllllll lllrillll lllilil Ecumcoo 050290 'lll’lll iiiiiii 1200 14001500 1 000 800 time(seconds) Figure 5.16 Calculated dielectric constant vs. time for the epoxy sample. 82 0 . w r 1 ABC - m wmm M mmm mmm .............. :m 1383 2 ddd 1 mmm U U U m ..m .m ..m ..................... 7.50 S SS 1 . u _ . . m . . 0 .. _- rim 0 rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr 0 6 _ r i--- ............. hm 4 0 ....... o, 110 _ 2 _ _ _ . _ _ . _ _ « ill II. _ i i..l 0 2 1 8. 6 4. 2 0 1. 0 o 0 0 .28.— mmo. 05856 time(seconds) Figure 5.17 Calculated dielectric loss factor vs. time for the epoxy sample. 83 400 ; i r r I l r r r t I I --..-.-i_-- .-_L_.-___¢.-_ L .1 _______ r.__. ..... '._a., 350. . ,4 . ~. . r r . l r I ,_ r temperature( ° C) _-_..t E , simulated with h=15 at point B r . 50 - - - - i - - H; ------- simulated with h=25 at point B — - -j ------- 3- - —~ ' i —1-— measured data at center ] ; ' 0 :___fl__ . l l i ' ' o 200 400 600 800 1000 1200 14001500 time(seconds) Figure 5.18 Parametric study on the heat transfer coefficient h for the epoxy sample. 84 temperature( ° C) 501. _.__J ....... L ...... J _______________ i ............ .. .., , ‘ ----- simulated w/o reaction heat at point B i ' I simulated with reaction heat at point B i , —-1— measured data at center J] ‘ i 0 r ’1" gif’ ' ' ”T:_—, -- e_J 0 200 400 600 800 1000 1200 14001500 time(seconds) Figure 5.19 Simulation with and without reaction heat for the epoxy sample. 85 300 :- 1 r' '1 _,.— —.~_- ‘1 r 1‘ : 11" i ~~.’l'-‘ . ‘ l i ,"l l 1 "mg,” I I .’ I l I l “mg“ 250i , ’f“;"' ‘4 '5 - “are: i l r i r I l t I j I t I 1 I I l r 1’ : 200' ------------- ,'i ~ 4 ------- . _._-. t 77777 5 .’ I I . l a, ,-’ ' : 1 : : ; 93 . : 1 :i ........... . ......... _,.. 7L6 150 i I l I 8 I : : : E I . . . 2 I l I 100 -------------- --------------- ; ----------- - -— l I C ' ‘ ‘,.H- ‘-'-'_- l r *\ ’,r"‘- i ----- simulated with constant a = 4 — 10.2 50 " ,7," ' ‘ """" T ‘ ' ' ' simulated with iarying e i"— ; ' """ simulated with varying g at 0% cure : —1— measured data at center 0 .-,, , i r . e . r ..f r:::::::‘_.4 0 200 400 600 800 1000 1200 14001500 time(seconds) Figure 5.20 Simulation on the variation effects of the complex permittivity for the epoxy sample. temperature distribution ( °C) _8 . . . . . . . . -8 —6 41 -2 0 2 4 6 8 Figure 5.21 Temperature distribution at cross section of point B at t = 100 seconds. temperature distribution ( °C) -8 -6 41 -2 0 2 4 6 8 Figure 5.22 Temperature distribution at cross section of point B at t = 300 seconds. temperature distribution ( °C) _8 . . . . . —8 -6 4 -2 O 2 4 6 8 Figure 5.23 Temperature distribution at cross section of point B at t = 400 seconds. temperature distribution ( °C) ‘8 L 4 A A A A r a e .4 -2 o 2 4 s a Figure 5.24 Temperature distribution at cross section of point B at t = 500 seconds. 88 temperature distribution ( °C) -8 -6 -4 -2 0 2 4 6 8 Figure 5 .25 Temperature distribution at cross section of point B at t = 700 seconds. temperature distribution ( °C) -8 —6 Figure 5.26 Temperature distribution at cross section of point B at t = 1000 seconds. 89 CHAPTER 6 ADAPTABLE MULTIMODE APPLICATOR MODELING In this chapter, the electromagnetic modeling of a novel adaptable multi-feed multimode cylindrical cavity applicator is presented, where the spatial distribution of the electric field can be specified a priori to accomplish a desired processing task. The electric field intensity inside the cavity can be tailored by varying the power delivered to each port, and mode-switching can be realized without mechanically adjusting the cavity dimensions. An orthogonal feeding mechanism is developed to reduce the cross coupling between the ports. Numerical simulations are performed for the cavity applicator to verify the theoretical analysis. 6.1 Introduction As mentioned earlier, microwave has been investigated as an attractive and efficient alternative energy source compared to inefficient pressurized ovens for material processing, including polymers and composites. However, industrial use of microwave processing has been impeded by the lack of proper applicator design, modeling, and control/monitoring methods. Most existing applicators are for specific applications only and the applicator design has to be performed over and over again for new processes, mostly by trial and error without the assistance of a model. The development of adaptable applicators, which can be configured to accomplish a variety of processing tasks, is therefore very important. Commonly used applicators for materials processing can be classified into three basic types: waveguide, multimode and single-mode applicators. A waveguide applicator is a hollow conducting pipe with either a rectangular or a circular cross-section. The 90 wave inside a waveguide applicator is fundamentally different from that inside multimode and single-mode applicator. The former is a traveling wave and the latter is a standing wave. Energy from the microwave generator travels through the waveguide and is partially absorbed by the material being processed. The remainder of the energy is directed to a terminating load. Traveling wave applicators are primarily used for continuous processing of high-loss materials; low-loss materials require an excessively long waveguide or a long processing period to absorb the necessary energy. The most popular applicator type is the overrnoded or multimode cavity where the electric field distribution is given by the sum of all the modes excited at a particular frequency. The frequent use of multimode cavity applicators is a result of their low cost, simplicity of construction, and adaptability to many different heating loads. This kind of applicator is very versatile in that it can accept a wide range of material loads of different dielectric losses, size and shape [61]. However, that may limit product quality, particularly with regard to the uniformity of temperature distribution in processed materials. Difficulties for multimode cavity analysis of electric field distributions result essentially from coexistence of many resonant modes. By rotating the sample and/or using metal stirrers, it is possible to improve the E-field uniformity and, thus, the heating uniformity inside the multimode applicator [62]. A mode-stirrer is a fan within a multimode cavity designed to change the resonance of multiple modes within a 4 MHz band near 2.45 GHz. As the mode-stirrer rotates, the resonant conditions of the cavity change, focusing the energy into different field patterns. The rotation of the fan cycles the cavity through the different resonant modes with a pattern that accepts whatever random heating occurs across the sample. This may improve the heating uniformity to a certain 91 degree, but makes it more difficult to predict the electric field distribution and is not suitable for precision material processing. The single-mode resonant applicator is designed to support only one resonant mode, therefore resulting in highly localized heating. These applicators can efficiently provide high field strength at mode-specific locations within the cavity since single frequency systems can be tuned for maximum throughput. Although single-mode applicators have a high efficiency relative to multimode and traveling wave applicators, sometimes it is very difficult for them to provide desired uniform heating across a large sample. To obtain uniform heating under these conditions, a technique called mode- switching was developed, in which several modes with complementary heating patterns are alternatively excited [63]. With a fixed frequency microwave power source, mode-switching can be achieved by mechanically adjusting the length of the cavity. This mechanical process slows the response of the system to temperature changes. With a variable frequency power source, modes can be changed by changing frequency. As a result of the instantaneous variable frequency mode switching, not only the speed of the process but also the controllability of the process is much improved [64]. However, variable frequency sources are inherently less efficient than the single frequency versions; the equipment is also very expensive. In this research effort, mode-switching can be obtained by varying the power delivered to multiple ports, hence eliminating the need for mechanical applicator adjustments. 92 In this case, the field intensity can be tailored to achieve high field strength in the regions of the applicator requiring high fields (and hence regions where heating is desirable) and low field intensities in regions where heating is not desired. Multi-port, multi-feed applicators will be studied where the field spatial distribution can be specified a priori to accomplish a desired processing task. Single frequency operation is preferred since the over—all system efficiency is generally higher for such “tuned” systems. Specifically, in our modeling, a cylindrical cavity is excited via two ports to get a TM mode and a TB mode simultaneously. With these separately controllable TM and TE modes, not only is the desired field strength distribution obtained at specific locations inside the cavity, but also the ability to heat along a preferred direction can be achieved; this is particularly useful when anisotropic materials are processed. Taking advantages of TM and TE modes, it is easy for us to get low port-to-port coupling, especially when a lossy load is present. This will be shown later in both theoretical analysis and numerical results. 6.2 Mode-Switching Typically, inside a cylindrical cavity, TM modes are excited with a coaxial probe while TE modes are excited with sidewall-mounted loops. Although these are useful methods for feeding, they may not lead to sufficient decoupling of the two ports. If the feeds are ideally orthogonal, the mutual coupling between the ports is zero (i.e., S2, = 0 = S”). This is a desirable condition since a non-zero transmission term indicates that power from one port will exit the cavity through the other port doing no significant heating of the polymer. Our cavity design has two feed ports. The feed mechanism 93 chosen, to minimize mode cross-coupling, are an axial slot (to excite 3 TE mode) and an azimuthal slot (to excite a TM mode). These slots are electrically very narrow to avoid mode cross-coupling. By observing the mode diagrams in Fig. 3.1 and Fig. 3.2, which show the resonant frequency as a function of cavity length for TM and TE modes in an empty cavity with radius a = 10.75 cm, it is seen that when the cavity height is adjusted, we can have different combinations of TM and TE modes resonant at 2.45 GHz, and thus get different electric field distribution inside the cavity. When the radius a = 10.75 cm and the height h = 7.34 cm, we can get the combination of mom and TE“, which is shown in Fig. 6.1(a); while a = 10.75 cm and h = 9.50 cm, the combination of mm and TE311 is shown in Fig. 6.1(b). For both cases, we have an axisymmetric electric field pattern (TMozo) at the center as well as some distributed spots (TEzn or TE311) around it. TM mode has the electric field orientation parallel to the central axis of the cylinder, which is best suited for heating rods and cylinders located along the axis of the cylinder; TE mode has the electric field polarization to the bottom plate of the applicator, which is best suited for heating flat panels or disks. For each combination of TM and TE mode with fixed cavity height, by varying the power delivered to each port, mode-switching can be obtained at a single frequency without mechanically adjustment of the cavity dimensions. This transition can be shown in fig. 6.2 for the combination of TMozo and TE” case. To investigate the electric field distribution for simultaneous mode excitation, a parameter a is introduced. A pure TE mode excitation has a = 0 while a pure TM mode excitation is given by or =1. The cavity model was verified for these two cases using the appropriate slot excitation methods. The 94 numerical simulations using COMSOL’s FEMLAB code [35] were performed for different or , which means different relative power delivered to each port. Fig. 6.2 illustrates the field distribution for each or , where (c) is the same as that in fig. 6.1(a) but with different scale. It was determined that a = 0.55 yielded the most uniform electric field for curing purpose in the center plane of the applicator. A uniform field is desirable since that will allow the curing of the largest diameter part. 95 2.5 (b) TMozo + T133” Figure 6.1 Electric field distribution for different combination of TM and TE modes. (b) a. = 0.25 97 (d)a=l Figure 6.2 Parametric study on the weighting factor or for the combination of the two modes. 98 6.3 Port-to-Port Coupling Analysis As mentioned earlier, the mutual coupling between the two ports is a very important factor that influences the performance of the applicator. To get better performance, lowest possible coupling is needed in order to avoid energy leakage through the ports. In a homogeneous, source-free cylindrical cavity with perfectly conducting walls, the electromagnetic fields inside the cavity can be derived from Maxwell's equations and the boundary conditions [1 7]. For TMozo mode, the fields are given by Ep=E¢=O 2 k 5.52 5.52 Ez= —Jo( p]=— WNO( p] (6.3) jam 0 Hp=Hz=0 H : _awozo _5.52J1(5.52p] ¢ 6p a a For TE211 mode, the fields are given by ___1_.‘9_W_2n__ fléfl __ Ep — p 6,, J2( p)sin(2¢)sin(’: ) __ 61/1211 .._ . 2 3.054 _ 2a 3.054 19,4 ‘—ap cos(2¢)sm( h)[J1("—P) 054p J2(—ap)] E, = 0 (6.4) _ 1 621/1211_ Ir/h _ 3054 _ Za 3.054 p— jaw apaz my COS(2¢)COS()[J1(_P) 3054p J2(—'—p)] 1 _l_62 _27r/h 3. 054 ”’2“: ——sin(2¢)cos est—)J2(—‘—p) jam 1) 64402 jaw H¢= 99 2 + kzy/211)= —1- ((02,118 — [E] )JZ (m p) COS(2¢) Sin(£) Jaw h h h 2 l 6 Hz: . ( Will 10):“ 62 Considering the field generated by azimuthal slot for TMozo mode, E p and E ¢ are zero, and E2 is very difficult to be coupled out through the axial slot; similarly, for the field generated by axial slot for TEZH mode, E2 is zero and it is very difficult for E p and E ¢ to be coupled out through the azimuthal slot. So it is seen that with these two specific modes, a very low port-to-port coupling can theoretically be obtained, which is essential for achieving improved efficiency. In practice, the microwave power is delivered into the cylindrical cavity by rectangular waveguides (for this work, WR-284) via iris-coupling (shown in Fig. 6.3 as the meshed geometry in F EMLAB). The dominant mode inside the waveguide is TE") mode. The dimensions of the irises need to be fine-tuned such that the S“, 821, 822, and 812 are minimized. Under these conditions, the maximum power is delivered to the load for processing purpose. Fig. 6.4 and Fig. 6.5 show the simulated S-parameters for this empty two-port cavity, where the TE-port is denoted as port 1 and TM-port is shown as port 2. The resonant frequency shifted down a little bit from 2.45 GHz due to the existence of the irises. After the lossy epoxy samples are loaded inside the cavity, the 821 and 812 are expected to be even smaller. Since the complex permittivity of epoxy resins varies as a function of temperature and extent—of-cure during microwave heating, the Sn and 822 are expected to vary significantly due to the loaded epoxy sample. In practice, in order to adjust the reflection coefficients at the two ports, a brass plunger can be 100 introduced lying very close to each coupling iris [65], such that the S]. and 822 can be retuned in real-time. Currently, the iris-coupled two-port cavity is being built and the feed network that includes the coax-to—waveguide transition is also being implemented. Experiments will be carried out to verify the model and the control strategies will be developed for microwave thermoset curing of epoxy resins. .- V Figure 6.3 The two-port cavity fed with waveguides via irises. 101 S-parameters 5 .2... 2 2 2 . 2 2 I1 . -li.-l 4““ 2 2 . _ 2 _ _ _ 2 _ 1 1.2 2 2 2 2 2 _ _ .2 .2 1 2 2 2 _ 2 2 _ 2 2 ...» S S 2 2 2 2 2 _ 2 2 h ._ 2 . 2 _ _ .. 2 2 2 M2 2+2 x2; _ 2 2 2 2 . 2 2 ..2 V X .1\_ .2 *7. m 2 _ _ . _ 2 _ _ .\J 2, . _ 2 2 2 2 2 2 _ 13 2. III .. 2. IIIIIII 2IIIhIIIrIII2IIILIIIklllFlIlFllrfifi 2 2 2 2 2 2 . 2 .v . .2 . 2 2 2 2 2 2 2 £2 4 2 2 _ 2 _ 2 2 _ .u 2 2 2 2 2 _ 2 _ 2 2 . 22 . 2 2 2 2 . 2 _ 2 .\ .fi 2 2 . 2 2 2 2 2 2 \. . 2 2 2 2 2 2 _ _ 2 .2 22. 2 2 2 2 _ 2 . 2 _ .31 22 . 2 2 2 2 2 2 2 2 H2 2‘III2II112|II+IIITIII2IIIJII14I|I2IIiI2III..fi .2. 2 _ 2 2 2 2 2 2 2 ... Z 2 2 2 2 _ 2 _ _ 2 2 _ .2 2 2 2 _ 2 2 2 2 _ w _ 2 2 2 2 2 _ 2 2 2 2 2 2 2 2 2 2 2 _ 2.x 2 2 2 2 2 2 _ 2 2 2 2 2 2 2 2 2 2 2 2 2. 22.. 2 2 _ 2 2 2 . 2 2 2X9 ...-I,_|--12.-y:2| r-yu_|;ut.r--2-t|r»--21 -. .2. 2 2 2 2 2 2 2 2 AM A _ 2 2 2 2 2 _ _ 2 m2 / 2 2 2 _ 2 2 2 2 .... It..- 2 2 2 2 _ _ 2 \n\t\\ . ...»... .P 2. 2 _ 2 2 22‘...\. _ 2 2 2 .2. 2.-.. Y\ 2 22L.“ _ .... _ _ _ _ 0.2 .Y 2 2 2 2 2 2 _ . 2.... 2...._ 2 2 _ 2 _ 2 2 .N7 IRI|2101JIII41I11I|12 IIIIIII 4IIIIIII2III.M ‘ _ 2 . u 2 2 2 2 _ . 2 2 2. 2 2 2 2 2 _ . 2 2 _ . . 2 2 2 2 2 2 _ 2 2 f . 2 \ 2 2 2 2 2 2 2 2 2 .272 x 2 2 2 _ 2 2 2 2 2 2% x. 2 2 2 2 2 2 2 2 2 22A 2 2 2 2 2 2 2 2 . fl 2 2 2 2 _ _ 2 m5 2 . 2 2 . F 2 .-.. 2 F111. M 1‘4 . 1 9 8. 7 6 5. 4 3. 2 1. OZ 0 0 0 O 0 0 0 0 0 _.Nw ._._.m Frequency (Hz) "" 43.} .IuSV Figure 6.4 TB mode (Sn and $21) for the two-port cavity. S-paramee rs 5 . fit .fi! 2 .2-T 2 .. 2 .../M .v. 2 2 2 _ 2 o 2 2* 2 _ 2 2 a2 82 81.2 2 2 2 2 _ .2 .22 2 _ . _ . _ ._ .2 _ _ . _ 2 2 ,\. .2 . 2.1. ‘ _ 2 2 _ 2 . ... , 2 2 2 2 ~. 2 2 2 2 2 2 m ¥ ...2...-..L|i.—.I IFIIL ILIlo .2 2 2 2 2 2 2 2H. . . 2 _ 2 2 2 2 ...22 M. 2 2 2 2 2 2 . 2 2 2 2 2 2 _. w 2 2 2 2 . 2 Mr 2 _ 2 2 2 2 g, . 2 2 2 _ 2 2 2 2 2 2 2 _ 2 ...1 2 2 2 2 2 2 ., A n.2il-.2| -«uiu11412..-u2-| % _ 2 2 2 _ 2 . 2 2 _ _ 2 2 2 U. 2 2 2 2 2 2 .. 2 2 2 2 2 2 . W _ _ 2 _ 2 2 .. 2 2 2 _ _ 2 . _ 2 2 2 2 2 \ . _ 2 2 _ . 2 . . .4 9 ...---2- Lwtlruuu..lxu2iun2n 2 2 .... I 2 2 2 2 _ 2 2 _ 2 ......2 M. X. 2 2 2 _ 2 2 2 _ 2 2..., 2 222:: 2 2 2 2 2 2 2 2 _>.\ 2 .4! ..2.. _ 2 2 2 2 .. 7.1. 2 2 _ 2 .....2. .... .\2 2 ...... 2. _ 2 2 2 + .2 _ 2 2 , z.. . _ . 2X ..2. . 2 _ 2 2 2 _. ......2 .-2 . 2 v... 2 _ 2 2 2 2 2 _ 2... h 2\ 2 2 2 2 2 _ 2 2 . ..7 IT..I.II..JII»1|II1|IJIIIJIII4|x11:111....42“ x 2 2 2 2 2 2 2 2 2 .. . 2. 2 2 2 2 2 2 2 2 . 2.2.2 ... 2 2 2 2 2 2 _ _ 2 :2. ‘ 2 2 2 2 2 2 2 2 Am . _ _ . _ 2 _ . 2 ..A. 2 2 2 2 2 _ _ . 2 X 2 2 _ 2 2 2 2 _ 2 _ +2. _ 2 2 2 _ 2 2 2 _ J5 ...-llrl. ILllIl-rl..||c1LlltlllLll.lul2 L *M 1 9. 8. 7. 6. 5. 4 3. 2 1. 02. 0 0 0 0 0 0 0 0 0 9 x10 Frequency (Hz) Figure 6.5 TM mode (822 and 812) for the two-port cavity. 102 CHAPTER 7 CONCLUSIONS In this study, a self-consistent 3D marching-in-time multiphysics model has been developed, which includes electromagnetic field distribution, microwave power absorption, heat transfer, and polymer curing kinetics. The temperature- and cure- dependent permittivity and curing kinetics for DGEBA/DDS based on experimental data were explicitly included in the model. Edge-based FEM was implemented for electromagnetic model to ensure the tangential continuity of electric field and divergence-free condition for source free region, while node-based FEM was used in thermal model to solve for the temperature distribution. The numerical results can be used to determine the time-dependent temperature distribution and curing profile across the polymer sample, as well as the electromagnetic field distribution within the cavity applicator. With the help of this numerical model, robust control strategies can be developed for the polymer curing process. The numerical results have been compared with the measured data and a good agreement was achieved. The electromagnetic modeling of a novel adaptable multi-feed multimode cylindrical cavity applicator was presented, where the spatial distribution of the electric field can be specified a priori to accomplish a desired processing task. It has been shown that the electric field intensity inside the cavity can be tailored by just varying the power delivered to each port, and the mode-switching can be realized without mechanically adjusting the cavity dimensions. An orthogonal feeding mechanism was developed to reduce the cross coupling between the ports. Numerical simulations were performed for the cavity applicator to verify the theoretical analysis. 103 CHAPTER 8 RECOMMENDATIONS FOR FUTURE WORK The recommendations for the future work include the following: 1. The dielectric property measurements need to be done at higher temperature, up to 250°C for DGEBA/DDS epoxy resins with different extents of cure at 2.45 GHz, in order to get accurate values of the complex permittivity. This may be difficult since the extent of cure keeps changing at high temperature. A better model than Davidson-Cole model should be investigated that has the ability to account for the high temperature and different extent of cure. 2. The automatic microwave processing and diagnostic system should be developed such that the dielectric constant and dielectric loss factor can be measured and recorded real time as the heating process is proceeding. 3. In this study, the input power delivered to the cavity was kept constant and recorded the temperature as it increased during the thermoset curing process. The main purpose of this is to test the numerical model. In practice, the heat temperature is usually kept constant during the curing process while the level of the input power is being adjusted. So in order to make the code useful for practical applications, especially on developing feedback control strategies, the modification should be done in the model accordingly. 4. After the two-port multimode cavity is built, experiments are expected to be carried out to verify the model and the control strategies should be developed for microwave thermoset curing of epoxy resins. A non- 104 invasive online monitoring technique should be developed as well, to provide feedback to the control system. 105 APPENDIX EXPRESSIONS FOR INTEGRAL TERM WITH DOT PRODUCT In section 4.2.2, to evaluate the integral F if = ”Le Wie -Wje.dV , first the dot product is obtained as [.9 5? w,e.we.= '1 LeLe L‘?Le Le. ' J (6Ve)2 Li111f’21211 12 £2le 12 11f’1j2+Li2L121’111 where j,-j =b.eb;+ cieci+dfdi Substituting this into the integral and utilizing the general integration formula given in (4.11), each entry in the elemental matrix can be obtained analytically [22]. For completeness, the expressions for these entries are listed below 8 (If? F11= e(f22-f12+f11) 360V ’f’§ F182= 8 (NB -f21-f13 +f11) 720V e e 1 Fle3= 38(2f24 -f21-f14 +f11) 720V [e e 1 4 F13= 720V V9 (f23 -f22 -2f13 +f12) e e 15 8.. 1[715‘ Ve(f22- f24- f12+2f14) le 1e F12 = 10: (f24 -f23 -f14 +f13) 720V F e -—(I§)2 (f f f ) = — + 22 360Ve 33 13 11 106 88 l 2 3 23: 8(2f34-f13-f14+f11) 720V 8 [3’14 6 (f33 -f23 -f13 +2f12) 720V 8 1515 e (f23 -f34 -f12 + f14) 720V 8 8 I 2 6e (f13 -f33 -2f14 +f34) 720V (13812 360Ve 88 = 3 4e (f34 -f24 -f13 +f12) 720V (f44-f14+f11) 8 8 I = 3 Se (f24 -f44 —2/12 +f14) 720V 88 = 3 6e (f44 -f34 -f14 +2f13) 720V (1:)?- 360Ve (f33 -f23 + f22) = 4 Se (f23 -2f34 -f22 +f24) 720V 8 8 l 4 6 e (f34 - f33 - 2f24 + f23) 720V (1:12 55: (f22—f24+f44) 360Ve 88 Fe —£1—6—(f —2f — + -720Ve 24 23 f44 f34) (1g )2 66: (f44—f34+f33)° 360Ve 107 10. 11. 12. REFERENCES J. Wei, M. C. Hawley, and J. D. Delong, “Comparison of Microwave and Thermal Cure of Epoxy Resins,” Polymer Engineering and Science, vol. 33, no. 17, pp. 1132- 1140, 1993. D. A. Lewis, J. C. Hedrick, T. C. Ward, and J. E. McGrath, “The Accelerated Curing of Epoxy Resins Using Microwave Radiation,” Polymer Preprints, vol. 28, no. 2, pp. 330-331, 1987. R. Agrawal and LT. Drzal, “Effects of Microwave Processing on Fiber-Matrix Adhesion in Composites,” Journal of Adhesion, vol. 29, pp. 63-79, 1989. J. Wei, J. Jow, J. D. DeLong, and M. C. Hawley, “Microwave Processing of Crossply Continuous Graphite Fiber/Epoxy Composites,” SAWE Journal, vol. 27, no. 1, pp. 33-39, 1991. M. LeGrand and V. Bellenger, “The Cure Optimisation of Carbon/Epoxy Prepregs,” Composite Science and Technology, vol. 58, no. 5, pp. 639-644, 1998. M. A. Acitelly, R. B. Prime, and E. Sacher, “Kinetics of Epoxy Cure: (1) The System Bisphenol-A Diglycidyl Ether/m-phenylene Diamine,” Polymer, vol. 12, no. 5, pp. 335-343, 1971. G. L. Hagnauer, P. J. Pearce, B. R. Laliberte, and M. E. Roylance, in T hermorheolog/ of Thermosetting Polymers, pp. 25-47, C. A. May, ed., Washington, DC, 1983. A. Dutta and M. B. Ryan, “Effect of fillers on Kinetics of Epoxy Cure,” Journal of Applied Polymer Science, vol. 24, no. 3, pp. 635-649, 1979. J. Mijovic, “Cure Kinetics of Neat Versus Reinforced Epoxies,” Journal of Applied Polymer Science, vol. 31, no. 5, pp. 1177-1187, 1986. A. Moroni, J. Mijovic, E. M. Pearce, and C. C. Foun, “Cure Kinetics of Epoxy Resins and Aromatic Diamines,” Journal of Applied Polymer Science, vol. 32, no. 2, pp. 3761-3773, 1986. J. Wei, M. C. Hawley and M. T. Demeuse, “Kinetics Modeling and Time- Temperature-Transforrnation Diagram of Microwave and Thermal Cure of Epoxy Resins,” Polymer Engineering and Science, vol. 35, no. 6, pp. 461-470, 1995. Z. Fathi, D. A. Tucker, W. A. Lewis, and J. B. Wei, “Industrial Applications of Variable Frequency Microwave Energy in Materials Processing,” Microwave Processing of Materials V, Materials Research Society Symposium Proceedings, vol. 430, pp. 21-28, 1996. 108 13. 14. 15. l6. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. V. Adegbite, “Automation and Control of the Microwave Processing of Composite Materials,” Ph.D. Dissertation, Michigan State University, 1995. Y. Qiu, V. Adegbite, and MC. Hawley, “Intelligent Composite Processing Using Variable Frequency Method: Preliminary Heating Results," 27'” International SAA/fl’E Technical Conference Proceedings, pp.173-183, 1995. Y. Qiu and M. C. Hawley, “Automated Variable Frequency Microwave Processing of Graphite/Epoxy Composite in a Single Mode Cavity,” Proceedings of the 13’ Annual Technical Conference, American Society for Composites, 1998. Y. Qiu, “Variable Frequency Microwave Processing and Microwave Process Control for Polymer Composites,” Ph.D. Dissertation, Michigan State University, 2000. R. F. Harrington, T ime—Harmonic Electromagnetic Fields, McGraw-Hill, New York, 1987. C. C. Ku, and R. Liepins, Electrical Properties of Polymers: Chemical Principles, Hanser, Munich, 1987. C. Saltiel and A. K. Datta, “Heat and Mass Transfer in Microwave Processing,” Advances in Heat Transfer, vol. 33, pp. 1-94, 1999. L. Zong, S. Zhou, R. Sun, L. C. Kempel, and M. C. Hawley, “Dielectric Analysis of a Crosslinking Epoxy Resin at a High Microwave Frequency,” Journal of Polymer Science Part B: Polymer Physics, vol. 42, no. 15, pp. 2871-2877, 2004. J. L. Volakis, A. Chatterjee, and L. C. Kempel, Finite Element Method for Electromagnetics: Antennas, Microwave Circuits, and Scattering Application, New York: IEEE Press, 1998. Jianming Jin, The Finite Element Method in Electromagnetics, second edition, Wiley-Interscience, New York, 2002. P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers (2nd edition). Cambridge: Cambridge University Press, 1990. Z. J. Csendes and P. Silvester, “Numerical Solution of Dielectric Loaded Waveguides: I —— Finite Element Analysis,” IEEE Trans. Microwave Theory Tech., vol. 18, no. 12, pp. 1124-1131,1970. X. Yuan, D. R. Lynch, and K. Paulsen, “Importance of Normal Field Continuity in Inhomogeneous Scattering Calculations,” IEEE Trans. Microwave Theory T ech., vol. 39, no. 4, pp. 638-642, 1991. H. Whitney, Geometric Integration Theory. Princeton University Press, NJ, 1957. 109 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. A. Bossavit and J. C. Verite, “A Mixed FEM-BIEM Method to Solve 3D Eddy Current Problems,” IEEE Trans. Magnetics, vol. 18, no. 2, pp. 431-435, 1982. J. C. Nedelec, “Mixed Finite Elements in R3,” Numer. Math, vol. 35, pp. 315-341, 1980. M. Hano, “Finite Element Analysis of Dielectric-Loaded Waveguides,” IEEE Trans. Microwave Theory T ech., vol. 32, no. 10, pp. 1275-1279, 1984. A. Bossavit, “Whitney Forms: A Class of Finite Elements for Three-Dimensional Computations in Electromagnetism,” IEE Proceedings, vol. 135, pt. A, no. 8, pp. 493-500, 1988. A. E. Bogle, “Electromagnetic Material Characterization Using a Partially Filled Rectangular Waveguide,” Master Thesis, Michigan State University, 2004. A. Melloni, M. Politi, and G. G. Gentili, “Mode-Matching Analysis of TEon-Mode Waveguide Bandpass Filters,” IEEE Trans. Microwave Theory T ech., vol. 43, no. 9, pp. 2109-2116, 1995. G. L. Matthaei, L. Young, and E. M. T. Jones, Microwave Filters Impedance- Matching Networks and Coupling Structures. New York: McGraw-Hill, 1964. F. Shnurer, “Design of aperture-coupled filters,” IRE Trans. Microwave Theory T ech., pp. 238-243, Oct. 1957. COMSOL, Inc., FEMAB User ’s Manual, version 3.0a, 2004. P. Ratanadecho, K. Aoki, and M. Akahori, “Experimental Validation of a Combined Electromagnetic and Thermal Model for a Microwave Drying of Capillary Porous Materials inside a Rectangular Wave Guide,” Journal of Microwave Power & Electromagnetic Energy, vol. 37, no. 1, pp. 15-40, 2002. L. Pichon and 0. Meyer, “Coupled Thermal-Electromagnetic Simulation of a Microwave Curing Cell,” IEEE Transactions on Magnetics, vol. 38, no. 2, pp. 977- 980, 2002. B. Wappling-Raaholt, et al., “A Combined Electromagnetic and Heat Transfer Model for Heating of Foods in Microwave Combination Ovens,” Journal of Microwave Power & Electromagnetic Energy, vol. 37, no. 2, pp. 97-111, 2002. H. Zhang and A. Datta, “Coupled Electromagnetic and Thermal Modeling of Microwave Oven Heating of Foods,” Journal of Microwave Power & Electromagnetic Energy, vol. 35, no. 2, pp. 71-85, 2000. Z. Huang, M. Iskander, J. Tucker, and H. Kimrey, “FDTD Modeling of Realistic Microwave Sintering Experiments,” Materials Research Society Symposium Proceedings, vol. 347, pp. 331-345, 1994. 110 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. V. Soriano, C. Devece, and E. De los Reyes, “A Finite Element and Finite Difference Formulation for Microwave Heating Laminar Material,” Journal of Microwave Power & Electromagnetic Energy, vol. 33, no. 2, pp. 67-76, 1998. L. Pham and K. Pope, “3D Finite Element Model of RF Heating: Novel Non- Ablative Cutaneous Therapy,” Proceedings of SPIE, vol. 4949, pp. 22-31, 2003. J. Haala and W. Wiesbeck, “Modeling Microwave and Hybrid Heating Processes Including Heat Radiation Effects,” IEEE Trans. Microwave Theory T ech., vol. 50, no. 5, pp. 1346-1354, 2002. T. V. C. T. Chan, J. Tang and F. Younce, “3-Dimensiona1 Numerical Modeling of an Industrial Radio Frequency Heating System Using Finite Elements,” Journal of Microwave Power & Electromagnetic Energy, vol. 39, no. 2, pp. 87-105, 2004. Q. Tang, N. Tummala, S. K. S. Gupta, and L. Schwiebert, “Communication Scheduling to Minimize Thermal Effects of Implanted Biosensor Networks in Homogeneous Tissue,” IEEE Trans. Biomedical Engineering, vol. 52, no. 7, pp. 1285-1294, 2005. K. G. Ayappa, H. T. Davis, G. Crapiste, E. A. Davis, and J. Gordon, “Microwave Heating: an Evaluation of Power Formulations,” Chemical Engineering Science, vol. 46, no. 4, pp. 1005-1016, 1991. 8. Zhou, “Interactions in Microwave Adhesive Bonding of Polymers and Composites,” Ph.D. Dissertation, Michigan State University, 2002. L. Zong, “Microwave Processing of Epoxy Resins and Synthesis of Carbon Nanotubes by Microwave Plasma Chemical Vapor Deposition,” Ph.D. Dissertation, Michigan State University, 2005. H. Lee and K. Neville, Epoxy Resins: Their Applications and Technology, McGraw- Hill, New York, 1957. J. 0. Simpson and S. A. Bidstrup, “Rheological and Dielectric Changes During Isothermal Epoxy-Amine Cure,” Journal of Polymer Science Part B: Polymer Physics, vol. 33, no. 1, pp. 55-62, 1995. D. W. Davidson and R. H. Cole, “Dielectric Relaxation in Glycerol, Propylene Glycol, and n-Propanol,” Journal of Chemical Physics, vol. 19, pp. 1484-1490, 1951. W. Kauzmann, “Dielectric Relaxation as a Chemical Rate Process,” Reviews of Modern Physics, vol. 14, no. 1, pp. 12-44, 1942. G. Gallone, S. Capaccioli, G. Levita, P.A. Rolla, S. Corezzi, “Dielectric Analysis of the linear Polymerization of an Epoxy Resin,” Polymer International, vol. 50, pp. 545-551, 2001. 111 54 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. J. Jow, “Microwave Processing and Dielectric Diagnosis of Polymers and Composites Using a Single-Mode Resonant Cavity Technique,” Ph.D. Dissertation, Michigan State University, 1988. Skyblue Systems, SK Y/Mesh2 user ’s manual, version 2.33, 1999. J. Meese, L. Kempel, and S. Schneider, “Subdividing Distorted Prisms into Tetrahedra,” 19th Annual Review of Progress in Applied Computational Electromagnetics, Monterey, CA, March 2003. D. R. Kincaid and T. C. Oppe, “ITPACK on Supercomputers,” Numerical Methods: Lecture Notes in Mathematics, vol. 1005, pp. 151-161, 1982. C. Lanczos, “Solution of Systems of Linear Equations by Minimized Iterations,” J. Res. Natl. Bur. Stand, vol. 49, pp. 33-53, 1952. J. Heinrich and D. Pepper, Intermediate Finite Element Method: Fluid Flow and Heat Transfer. PA: Taylor & Francis, 1999. L. Zong, G. L. Charvat, L. C. Kempel, and M. C. Hawley. “Automated method for characterizing temperature dependent dielectric materials,” 26th AMTA Meeting & Symp., Stone Mountain, GA, Oct., 2004. J. E. Gerling, “Microwave Oven Power: a Technical Review,” Journal of Microwave Power and Electromagnetic Energy, vol. 22, no. 4, pp. 199-207, 1987. A. C. Metaxas and R. J. Meredith, Industrial Microwave Heating, pp. 137-138, IEE Power Engineering Series, Peter Peregrinus Ltd., London, 1983. L. A. Fellows and M. C. Hawley, “Microwave Heating of Composites with Complex Geometries in Cylindrical Single-Mode Resonant Microwave Cavities,” Proceedings of the 29’” Microwave Power Symposium, pp. 91-94, Chicago, IL, July 1994. Y. Qiu and M. C. Hawley, “Computer Controlled Variable Frequency Microwave Processing of Complex-Shaped Graphite/Epoxy Composite,” Proceedings of the 14th American Society for Composites Annual Technical Conference, pp. 52-61, Dayton, OH, September 1999. L. Chen, C. K. Ong, and B. T. G. Tan, “A Resonant Cavity for High-Accuracy Measurement of Microwave Dielectric Properties,” Measurement Science and Technology, vol. 7, no. 9, pp. 1255-1259, September 1996. 112 l111111111111111111111111111l