z... .I.
.. : L ..

ﬁr 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 . ‘ . . gﬂbvxvpnwtm {0.0va
.égﬁﬁsﬁa4 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 fulﬁllment
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 Afﬁrmative 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 fulﬁllment 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
efﬁciency not available with conventional convection heating. However, due to the
complex interactions between the electromagnetic ﬁelds 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 ﬁeld 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 ﬁeld 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 ﬁnite element method (FEM) is implemented as an

electromagnetic model to ensure the tangential continuity of electric ﬁeld 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 proﬁle across
the polymer sample, as well as the electromagnetic ﬁeld 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 ﬁeld can be speciﬁed a priori to accomplish a desired processing task. The
electric ﬁeld 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 conﬁdence 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 ﬁlled waveguide. ......................................... 45
Figure 4.5 Schematic view of the TEm l-mode bandpass ﬁlter [32]. ................................ 47
Figure 4.6 Meshed geometry of the bandpass ﬁlter .......................................................... 48
Figure 4.7 Calculated magnitude of 821 for the bandpass ﬁlter with FEM ....................... 49
Figure 4.8 Magnitude of S21 for the bandpass ﬁlter 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 modiﬁed
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 ﬁeld 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 coefﬁcient 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 ﬁeld 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 deﬁnition 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 ﬁelds 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 ﬁeld 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, ﬁrst 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 ﬁelds in a cylindrical single mode cavity and microwave/materials
interactions during microwave heating.

In Chapter 4, node-based ﬁnite 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
ﬁeld 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 ﬁeld 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 proﬁle across
the polymer sample, as well as the electromagnetic ﬁeld 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
ﬁeld can be speciﬁed a priori to accomplish a desired processing task. The electric ﬁeld
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
ﬁelds 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 efﬁciency is controlled by the heat transfer coefﬁcient at the
material surface and the thermal conductivity of the material. During thermal heating,
heat ﬂows 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 difﬁculty 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 ﬁeld 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
efﬁcient) 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 ﬁber 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 difﬁculties 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 ﬁeld 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 difﬁcult 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 difﬁcult 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 efﬁcient)
approach to inefﬁcient 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 ﬁrndamental 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 ﬁeld) 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 ﬁeld can be speciﬁed 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 ﬁeld distribution within the applicator during
processing and the resulting temperature proﬁle 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 efﬁcient 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 proﬁle across the polymer sample, as well as the electromagnetic
ﬁeld 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 ﬁnite element method (FEM) with tetrahedral
elements for the applicator and polymer sample since this method offers a high degree of
geometrical and electrical ﬁdelity. In particular, the spatial variation in permittivity can
readily be accommodated and by using “snapshots-in-time”, time-dependant thermal and
composition inﬂuences 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 ﬁeld
oscillation within the applicator.

The ﬁnite 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 efﬁcient manner. Once the ﬁelds 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 ﬁrst-
order ordinary differential equation (ODE) system is obtained and the ﬁnite 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 Teﬂon/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 ﬁelds 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 ﬁeld 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 ﬁxed 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 ﬂat 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 efﬁcient 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.

Speciﬁcally, this microwave applicator is designed with the goal of:

10

1. Reducing the processing cycle time.

2. Maximizing the processing efﬁciency (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 ﬁeld intensity to
achieve high ﬁeld strength in the regions of the applicator requiring high ﬁelds (and
hence regions where heating is desirable) and low ﬁeld 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 ﬁeld strength and distribution patterns are essential factors that
inﬂuence microwave heating efﬁciency 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 ﬁelds,
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 useﬁrl 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 ﬁeld intensity, H is the magnetic ﬁeld intensity, D is the electric
displacement density or electric ﬂux density, B is the magnetic ﬂux density, J is the

electric current density, and p is the charge density. D is deﬁned 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 ﬁeld intensity, and M is the
volume magnetization, the measure of the density of magnetic dipoles in the material. In
a simple isotropic medium, these ﬁeld 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 ﬁelds.

Maxwell’s equations are ﬁrst-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 ﬁelds 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 ﬁeld
components are transverse, and the magnetic ﬁeld components are parallel to the
direction of wave propagation, which is the axial direction. In TM modes, the electric
ﬁeld components are parallel, and the magnetic ﬁeld 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
Longthofcavltﬂcm)

[+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 ﬁelds inside the cavity can be derived from
Maxwell's equations and boundary conditions.

For TE modes, the electromagnetic ﬁeld components inside an empty cavity are

[17]:

 

 

 

_ 1 ‘32qu E ___1_3"1ﬂ
” 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 ﬁeld 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<ﬂz) (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 ﬁeld
inside the cavity. The most widely used numerical techniques include the method of
moments (MOM), the ﬁnite element method (FEM), and the ﬁnite difference time

domain (FDTD).

3.2 Microwave/Materials Interactions

3.2.1 Mechanisms of Microwave/Materials Interactions

Non-magnetic materials are classiﬁed into conductors, semiconductors and
dielectrics according to their electric conductivity. Conductors contain free charges,
which are conducted inside the material with alternating electric ﬁelds so that a
conductive current is induced. Electromagnetic energy is dissipated into the materials
while the conduction current is in phase with the electric ﬁeld inside the materials.
Dissipated energy is proportional to conductivity and the square of the electric ﬁeld
strength. Conduction requires long-range transport of charges.

In dielectric materials, electric dipoles are created when an external electric ﬁeld

is applied and they will rotate until they are aligned in the direction of the ﬁeld.

19

Therefore, the normal random orientation of the dipoles becomes ordered. These ordered

polar segments tend to relax and oscillate with the ﬁeld. 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 ﬁeld. 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 ﬁeld. When no electric ﬁeld is
applied, the center of positive charges (nucleus) coincides with the center of
negative charges (electron cloud). When an external electric ﬁeld 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 ﬁeld is applied, the positive charges move in the

direction of the ﬁeld 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 ﬁeld. At the
presence of external electric ﬁeld, the dipoles will rotate until they are aligned
in the direction of the ﬁeld. The dipolar rotation of molecules is accompanied
by intermolecular friction, which is responsible for heat generation. Orientation
polarization is ﬁrndamentally different from electronic and atomic polarization.
The latter is due to the fact that the external ﬁeld induces dipole moments and
exerts displacing force on the electrons and atoms, while the orientation
polarization is because of the torque action of the ﬁeld 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 ﬁeld.

21

3.2.2 Dielectric Properties

Most polymers and composites are non-magnetic materials. The electromagnetic
energy loss is only dependent on the electric ﬁeld. Incident electromagnetic ﬁelds 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 deﬁned 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 ﬁeld, materials will interact with the
oscillating electromagnetic ﬁeld at the molecular level. Different materials will have
different responses to the microwave irradiation. Microwave heating of conductive
materials, such as carbon ﬁbers and acid solutions, is mainly due to the interaction of the
motion of ions or electrons with the electric ﬁeld. However, conductors with high
conductivity will reﬂect the incident microwaves at their surface and accordingly can not
be effectively heated.

The ﬁelds 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. Deﬁned as the distance into the sample, at

which the electric ﬁeld strength is reduced to He, the skin depth is given by [17]:

ds = 1 (3.21)

1 1/2
[Ewﬂo/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 ﬁeld is reﬂected and no electric ﬁeld 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 inﬁnite. 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 ﬁbers,
and Kevlar ﬁbers, is mainly due to the interaction of the motion of dipoles with the
alternating electric ﬁeld. 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 ﬁbers. When non-conducting ﬁbers, such as glass,
are used, microwaves will selectively heat the polymer matrix. When conducting ﬁbers
like graphite are used, energy is preferably absorbed by the conductive ﬁbers and heat is
locally conducted from the ﬁber to the matrix. In this case, loss factor is mainly due to
the ﬁber 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 ﬁeld strength and power
absorption in the materials. The change in electric ﬁeld 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 ﬁeld 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 ﬁnite element method (FEM) is widely used in computational
electromagnetics since it is very efﬁcient 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 brieﬂy 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 ﬁnite 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 ﬁnite 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 coefﬁcients. 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 coefﬁcients 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 coefﬁcients 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 ﬁrst-order elements)

ue(x,y,z)=ae +bex+cey+dez (4.7)

29

3

Figure 4.1 Linear tetrahedral element.

The coefﬁcients 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 coefﬁcients 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 ﬁelds. 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 difﬁcult to represent conducting and dielectric edges and

comers due to ﬁeld singularities associated with these structures [25].

32

Edge-based ﬁnite elements, whose degrees of freedom are associated with the
edges and the faces of the ﬁnite element mesh, have been shown to be free of the above
shortcomings. Edge basis ﬁrnctions 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 ﬁnite element basis functions that span a curl conforming space with degrees

of freedom associated with the edges, faces, and elements of a ﬁnite 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 deﬁned according to Table 4.1, then edge bases for a triangular element are deﬁned 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 ﬁeld 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 ﬁeld along the ith edge. It can be easily shown that the

edge-based functions deﬁned 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 ﬁrnctions for these
elements follow the same pattern as that for triangular vector basis functions. Consider
the tetrahedron shown in Fig. 4.1 and deﬁne 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 ﬁeld within the element can then be expanded as

6
e _ e e
E _ ZEiWi (4.15)

i =1
where the coefﬁcients E 1.9 represent the average value of the ﬁeld along the ith edge of

th
the e element.

Table 4.2 Edge deﬁnition 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 ﬁeld turns around the axis 3-4 and is normal to planes containing nodes 3 and
4. The ﬁeld thus has only tangential continuity across element faces. Edge elements can
also be described as Whitney elements of degree one and can be broadly classiﬁed as
belonging to the H 0(curl) space.

When the edge-based basis functions are employed for a three—dimensional ﬁnite
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 deﬁned 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 sufﬁciently 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 speciﬁc
conﬁguration 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 inﬁnitely long, the

region of interest extends to inﬁnity. To apply the ﬁnite element method to such an

37

unbounded region, it is ﬁrst necessary to truncate the inﬁnite region to a ﬁnite region.

This can be accomplished by placing ﬁctitious planes on each of the two sides of the

obstacles. These two ﬁctitious 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 deﬁne 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 sufﬁciently 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 ﬁeld can be expressed as a superposition of an incident TE10
wave and a reﬂected 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 reﬂection

coefﬁcient, 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(ﬁxE)=Uinc (4.23)
with

7 = 1k210a Umc = -21'kz10Ei"c-

Similarly, the ﬁeld 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 coefﬁcient. From this, the following condition is

obtained

39

x(VxE)+yﬁx(ﬁxE)=0 (4.25)
valid for the ﬁeld 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 ﬁeld within the region of interest and the boundary condition
at the waveguide walls, can be solved uniquely using the ﬁnite 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)-(ﬁxE)]dS

(4.27)

For the numerical discretization of (4.27), ﬁrst, the volume V is subdivided into

small volume elements and within each element the ﬁeld 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 ﬁeld 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; = ﬁxWiS , 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 ﬁelds

coincident with the waveguide wall, for the edge ﬁelds inside Vas well as on S 1 and $2.

Once (4.33) is solved for the ﬁelds at S 1 and Sz, the reﬂection and transmission

coefﬁcients can be obtained from (4.21) and (4.24), respectively. Speciﬁcally, from

(4.21), the reﬂection coefﬁcient 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 coefﬁcient can be obtained as

E(x, y, 2)
E0810 (x. y)e"jkz‘oz

 

T:

 

2= 22 (4-35)

Ideally, the reﬂection and transmission coefﬁcients are not functions of (x,y).
However, in a numerical solution, when S 1 and S2 are not placed sufﬁciently far away
from the obstacle one ﬁnds 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 reﬂection and transmission coefﬁcients 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 inﬂuence of higher-

42

order modes excited by the obstacle on the calculation of the reﬂection and transmission

coefﬁcients 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 ﬁxed 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 ﬁxed 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 ﬁlled 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 ﬁrst test involves reversing the extraction process to
generate scattering parameters for the special case of a completely ﬁlled 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 ﬁnite element method for electromagnetic model, the completely ﬁlled waveguide
case is studied with the FEM code and is compared with mode-matching results from
[31].

The conﬁguration 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 ﬁlls 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 ﬁeld as well as the reﬂection and transmission coefﬁcients
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' _._,ﬁ_‘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 ﬁlled 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 ﬁlters at millimeter wavelengths. Direct
coupled ﬁlters 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 ﬁlters.

In [32], a mode-matching technique was presented to analyze accurately this kind
of ﬁlters. 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 ﬁlter 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 ﬁlter [32].

Ii

 

 

 

The studied conﬁguration 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 ﬁlter.

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 ﬁlter 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 ﬁlter 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
signiﬁcant after the curing process starts.

In this chapter, a self-consistent 3D marching-in-time multiphysics model is
presented, which includes electromagnetic ﬁeld 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 ﬁeld 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 proﬁle across the polymer sample, as well as the electromagnetic
ﬁeld 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 ﬁeld strength and the

perrnittvity as below:
1 ,, It
P = Ewsogr (a,T)E - E (5.2)
where E is the electric ﬁeld 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 ﬁeld is

52

determined by solving the vector wave equation using ﬁnite 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 efﬁcient steady-state analysis method whose results are used in solving the
energy balance equation for heat transfer problem. Speciﬁcally, in the heat transfer model,
after the node-based FEM for spatial variations is performed, a ﬁrst-order ordinary
differential. equation (ODE) system is obtained and the ﬁnite 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 ﬁeld 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 proﬁle 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 ﬁts 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,
speciﬁcally, 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
4Iﬁ+ + 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-Tﬁhfﬁ'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-ﬁ—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 ﬁrst 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 deﬁned 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 vitriﬁcation, 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 ﬁt 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. Speciﬁcally, 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 ﬁt 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—ﬁ' 1+1 rjﬁ' 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 modiﬁed such that they follow the same trend

as the dielectric constant and are shown in Fig. 5.8. Speciﬁcally, 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 difﬁcult 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 modiﬁed 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 modiﬁed 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" )—Modiﬁed 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 modiﬁed

Davidson-Cole model for a wider temperature range.

5.5 Electromagnetic Model

The edge-based ﬁnite element method with tetrahedral elements for the applicator
and polymer sample is used since this method offers a high degree of geometrical and
electrical ﬁdelity. In particular, the spatial variation in permittivity can readily be
accommodated and by using “snapshots-in-time”, time-dependant thermal and
composition inﬂuences 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 ﬁeld

propagation within the applicator.

66

The ﬁnite 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 efﬁcient manner. Once the ﬁelds are determined
throughout the cavity for a given excitation, local power absorption rate can be
determined for the heat transfer model.

Speciﬁcally, the conﬁguration 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 ﬁgure. 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 ﬁeld distribution for TMmz
mode is shown in Fig. 5.10. The epoxy sample is centered at the middle of the cavity,
where the electric ﬁeld 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 Teﬂon holder with the wall

thickness of 0.2 cm and the base thickness of 0.3 cm. The Teﬂon has the complex

permittivity a: = 2.1—j0.0004.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

‘ 7mm
PointC , F "X'—
t 7mm
PointB > —}%
6mm
PointA F %F 7
mm
Teﬂon _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 ﬁeld distribution for TMon mode inside the

cavity.

Regarding the mesh generation for computational purpose by ﬁnite 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
ﬂexibility on geometry.

The electric ﬁeld, inside the cavity and at each point within the epoxy sample
being processed, can be computed using edge-based ﬁnite element method. The set of
ﬁnite 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 ﬁmction is denoted by W,- and the vector expansion function for
the total electric ﬁeld 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 ﬁeld.
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, ﬁxE = 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 ﬁeld 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. Speciﬁcally, 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 ﬁrst-order ordinary differential equation (ODE) system
is obtained and the ﬁnite difference scheme for temporal variations will be used for time-
marching process.

Recall the energy balance equation for the epoxy sample and Teﬂon 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 Teﬂon holder are listed in Table 5.5. For convenience, there is no

difference in symbols between the epoxy resins and the Teﬂon 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
Teﬂon
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 Teﬂon 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—gGAHﬂNidV

+ LpCpNiNj

In the above equation, the natural air free convection boundary condition is

already applied between the Teﬂon/air or epoxy/air interfaces, as
—(IE-VT)-ﬁ=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 coefﬁcient 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 inﬂuence of its variation on numerical results.

For each tetrahedral element, the elemental matrix is given by
e
[Cel{%}+[1<e]{re}={ge} (5.12)

where

e _ . . 1
Ci]. -- LpCpN,NJdP

71

 

 

 

6N'5N. 6N~ . aN- .
K€=L kx—L ’+ky—i§N—'+kz J 6N’ V+j hNiNj-dS
'1 0x 0y 6y 62 62 S

Q? = [lg—meoaxafmymdm JShTwNidS+ Lp%%(—AHr)NidV

(5.13)
Performing the integrations and assembling, the linear system representing the

ﬁrst-order ordinary differential equations (ODE) is given as
dT
[01(3} + [K1{T} = {Q} (5.14)

The ﬁnite difference scheme known as the O-method [59] is used to solve the
system and to obtain the updated temperature for the sample and Teﬂon holder. The
dielectric property of the epoxy sample is updated with the new temperature and extent-
of-cure, then the electric ﬁeld will be solved again with edge-based FEM to get the new
power absorption rate and solve for the heat transfer equation. This process is repeated
and the temperature proﬁle can be obtained across the sample for any given time.

Figure 5.11 shows the ﬂow chart for this coupled electro-, thermo- and kinetic
model. Essentially, the material property (complex permittivity for the DGEBA/DDS) is
updated with the current local temperature and extent-of-cure according to the modiﬁed
Davidson-Cole model with a sufﬁciently ﬁne timestep as the thermoset curing is
proceeding. At the same time, the time derivative of extent-of-cure (da/dt), which will be
used in the energy balance equation, also is updated according to the kinetics model and
by the interpolation between curves with different heating temperature. The extent-of-
cure distribution also is updated by performing the integration of the da/dt over the

corresponding timestep. Therefore, in this multiphysics model, the node-based FEM for

72

heat transfer model together with the ﬁnite difference scheme, and the edge-based FEM
for the electromagnetic model, will be called alternatively during the time-marching

process.

Take the thermal model as main function with initial
conditions and boundary conditions, starting mesh, at time t

 

 

J
l
\
Call EM model as a function, calculating electric ﬁeld
distribution as well as power, with up dated 6," (T,a)
_.a

 

I

Increase time by At, solve energy balance equation with up dated
power absorption rate and daldt, get up dated temperature and 1“

extent—of—cure distribution at t' = t + At _/

 

From modified Davidson-Cole model, get up dated S,* (T,o.) at
t'=t+At, and get up date time derivative of extent—of-cure
(da/dt) from kinetics curves _2

‘ Yes
E Is time t'>tem1 E

No
Yes

No
‘ ﬂ: 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-ﬁeld 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 ﬂuoroptic temperature sensing device
(Luxtron 750), which does not absorb electromagnetic energy and minimally interferes
with the electromagnetic ﬁeld, is used to on-line monitor the material temperature. The
ﬂuoroptic probe is protected by a 3 mm O.D. pyrex capillary tube. A cylindrical Teﬂon
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 puriﬁcation. 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 ﬁve minutes). Finally, the resins were degassed at 0.02 bar at
100°C for ﬁve minutes.

Before the epoxy sample held by the Teﬂon 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 Teﬂon
holder. The dimensions of the epoxy resins and the Teﬂon 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 Teﬂon holder with a ﬂuorOptic 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 Teﬂon holder
such that the sample is located at the place where the electric ﬁeld is the strongest.

The single frequency microwave curing and diagnostic system was used to heat
the liquid samples in the Teﬂon 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 speciﬁed 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. Thereaﬁer,

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 ﬁxed at 5
Watts and the temperature was recorded as it increased. The temperature/time proﬁle
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 ﬁg. 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 coefﬁcient, which was used in
the convective boundary condition. It can be seen that the calculated temperature proﬁle
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 ﬁnal 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 coefﬁcient 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 signiﬁcant

. . . . . . 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 ﬁrnction 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 modiﬁed 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 speciﬁcally 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 ﬁn

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 :___ﬂ__ . l l i ' '

o 200 400 600 800 1000 1200 14001500
time(seconds)

 

 

 

Figure 5.18 Parametric study on the heat transfer coefﬁcient 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 ﬁeld can be speciﬁed a priori to accomplish a desired processing task. The
electric ﬁeld 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
efﬁcient alternative energy source compared to inefﬁcient 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 speciﬁc 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 conﬁgured to accomplish a variety of processing
tasks, is therefore very important.

Commonly used applicators for materials processing can be classiﬁed 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 ﬁeld 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. Difﬁculties for multimode cavity analysis of electric ﬁeld 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-ﬁeld 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 ﬁeld 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 difﬁcult to predict the electric ﬁeld 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 efﬁciently
provide high ﬁeld strength at mode-speciﬁc locations within the cavity since single
frequency systems can be tuned for maximum throughput. Although single-mode
applicators have a high efﬁciency relative to multimode and traveling wave applicators,
sometimes it is very difﬁcult 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 ﬁxed 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 efﬁcient 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 ﬁeld intensity can be tailored to achieve high ﬁeld strength in the
regions of the applicator requiring high ﬁelds (and hence regions where heating is
desirable) and low ﬁeld intensities in regions where heating is not desired. Multi-port,
multi-feed applicators will be studied where the ﬁeld spatial distribution can be speciﬁed
a priori to accomplish a desired processing task. Single frequency operation is preferred
since the over—all system efﬁciency is generally higher for such “tuned” systems.

Speciﬁcally, 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 ﬁeld strength distribution obtained at speciﬁc 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 sufﬁcient 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 signiﬁcant

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 ﬁeld 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 ﬁeld pattern
(TMozo) at the center as well as some distributed spots (TEzn or TE311) around it. TM
mode has the electric ﬁeld 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 ﬁeld polarization to the bottom plate of the applicator, which is best
suited for heating ﬂat panels or disks.

For each combination of TM and TE mode with ﬁxed 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 ﬁg. 6.2 for the combination of TMozo and TE” case. To investigate the electric ﬁeld
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 veriﬁed 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 ﬁeld distribution for each or , where (c) is the same as that in ﬁg. 6.1(a) but
with different scale. It was determined that a = 0.55 yielded the most uniform electric
ﬁeld for curing purpose in the center plane of the applicator. A uniform ﬁeld is desirable

since that will allow the curing of the largest diameter part.

95

2.5

 

 

   

(b) TMozo + T133”

Figure 6.1 Electric ﬁeld 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 inﬂuences 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 ﬁelds inside the cavity can be derived from Maxwell's

equations and the boundary conditions [1 7]. For TMozo mode, the ﬁelds 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 ﬁelds are given by

 

 

 

___1_.‘9_W_2n__ ﬂéﬂ __
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 ﬁeld generated by azimuthal slot for TMozo mode, E p and E ¢ are
zero, and E2 is very difﬁcult to be coupled out through the axial slot; similarly, for the

ﬁeld generated by axial slot for TEZH mode, E2 is zero and it is very difﬁcult for E p and

E ¢ to be coupled out through the azimuthal slot. So it is seen that with these two speciﬁc

modes, a very low port-to-port coupling can theoretically be obtained, which is essential
for achieving improved efﬁciency.

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 ﬁne-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 signiﬁcantly due to the loaded epoxy sample. In practice, in

order to adjust the reflection coefﬁcients 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 2IIIhIIIrIII2IIILIIIklllFlIlFllrﬁﬁ

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 .\
.ﬁ 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..ﬁ
.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 .
ﬂ 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 .ﬁ! 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 ﬁeld 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 ﬁeld 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 proﬁle across
the polymer sample, as well as the electromagnetic ﬁeld 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
ﬁeld can be speciﬁed a priori to accomplish a desired processing task. It has been shown
that the electric ﬁeld 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 difﬁcult 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 modiﬁcation 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 , ﬁrst 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 ﬁllers 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/ﬂ’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 ﬁlters,” 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